Hi

First, use the MATLAB tools for Nifti to load the file:

http://uk.mathworks.com/matlabcentral/fileexchange/8797-tools-for-nifti-and-analyze-image

To load the image type:

nii_file = load_untouch_nii(file_path)

This yields a nifti struct nii_file. Now you have to look into this struct for the quaternion information to recover the axis rotation matrix. This information is strored in nii_file.hdr.hist

There are 6 fields in nii_file.hdr.hist you're interested in

quatern_b : b value of quaternion

quatern_c : c value of quaternion

quatern_d : d value of quaternion

qoffset_x : x component of translation vector

qoffset_y : y component of translation vector

qoffset_z : z component of translation vector

Now you have to recover the rotation matrix from the quaternion values a, b, c, and d where a = sqrt(1.0-(b^2+c^2+d^2))

R =  [ a*a+b*b-c*c-d*d        2*b*c-2*a*d          2*b*d+2*a*c    ;

        2*b*c+2*a*d          a*a+c*c-b*b-d*d        2*c*d-2*a*b    ;

        2*b*d-2*a*c              2*c*d+2*a*b       a*a+d*d-c*c-b*b ]

And the translation vector is

 t =  [qoffset_x; qoffset_y; qoffset_z]

Now for each image location [i; j; k], the world coordinate [x; y; z] is

[x; y; z] = R * [i; j; k] + t

I realised this needs a small correction, since the calculation before works only for images with 1x1x1mm voxel size

You need to first multiply each component of the [i; j; k] vector by the corresponding voxel sizes pixdim(2), pixdim(3) and pixdim(4) (Note you have to consider the 2nd to 4th element of the vector, as the first one is just a flag), like this

[i2; j2; k2] = [i * pixdim(2); j * pixdim(3); k * pixdim(4)] (Note you can do this in vector form using pixdim(2:4)). You can find the pixdim vector in nii_file.hdr.dime

finally the world coordinates are computed as before with this new vector

[x; y; z] = R * [i2; j2; k2] + t