3D interpolation Dicom data

I have a series of dicom images which i would like to extend by interpolation, the old series is a [256 256 166] matrix with a voxel size of [0,9375 0,9375 1,200]. I would like to interpolate this in all 3 dimensions to create smaller voxels. When i try to apply the function interp3() i immediately get an error: Sample values must be a single or double array. But when i rewrite my data to single or double array i will lose all the information in my image.

Could anybody tell me how to apply this method or which method i should use to interpolate without losing the intensity values of my original data?

Answer 1

[Alex Taylor]

If you have the Image Processing Toolbox, I would make this problem easier on yourself by using imwarp and possibily imref3d. You can represent what you are trying to do as a geometric transformation of your input image with a scale transformation. It sounds like your intention is to choose a scale factor > 1 to increase the grid size of your image.

Using imref3d will help you define the world location of your input image from the metadata in your dicom image.

S = load('mri');
I = squeeze(S.D);
RI = imref3d(size(I),0.9375 0.9375 1.200);

scaleFactor = 2;
tformResize = affine3d([scaleFactor 0 0 0; 0 scaleFactor 0 0; 0 0 scaleFactor 0; 0 0 0 1]); 
[Inew,Rnew] = imwarp(I,RI,tformResize);

Conceptually, it should be possible to do what you are trying to do with interp3. My guess is that interp3 is returning a single/double image outside the range [0 1] because you are probably just using simple cast single(I) or double(I) instead of rescaling and casting (im2double). Many MATLAB/IPT functions that require a dynamic range assume that single/double data is normalized to the range [0 1]. I bet you have not actually lost information if you are interpolating correctly. You either need to adjust the dynamic range used in visualizing your floating point data:

imshow(I(:,:,10),[]);

Move your output image back to its original datatype, for example uint8, or rescale your output image to the range [0 1]:

I = I./max(I(:));

[Hugo]

Thank a lot, i will explore this possibilities. Furthermore could you explain why you use squeeze (S.D) instead of squeeze (S), what does this .D do?

[Hugo]

I tried to apply the imref method, suggested above, but when i apply the last line (imwarp) i get the next error messages:

Error using - Out of memory. Type HELP MEMORY for your options.

Error in images.spatialref.internal.SpatialDimensionManager/intrinsicToWorld (line 220) worldCoordinate = self.StartCoordinateInWorld + (intrinsicCoordinate-0.5)*self.DeltaNumerator/self.DeltaDenominator;

Error in imref3d/intrinsicToWorld (line 214) zw = self.Dimension.Z.intrinsicToWorld(zIntrinsic);

Error in imwarp>remapAndResampleGeneric3d (line 214) [dstXWorld, dstYWorld, dstZWorld] = outputRef.intrinsicToWorld(dstXIntrinsic,dstYIntrinsic,dstZIntrinsic);

Error in imwarp>remapPointsAndResample (line 195) outputImage = remapAndResampleGeneric3d(inputImage,R_A,tform,outputRef,method,fillValues);

Error in imwarp (line 175) outputImage = remapPointsAndResample(inputImage,R_A,tform,outputRef,method,fillValues);

Can anybody explain me what is going worng, and which adjustments i should make to make it work?

[Alex Taylor]

Hi Hugo,

Looks like you are running out of memory. There are two possibilities:

1) The syntax you are using to call imwarp is somehow incorrect and causing imwarp to attempt to generate a huge volumetric image.

2) The syntax is correct, and you simply do not have enough memory to perform the 3D resize on your machine.

To be able to help you, Can you provide the output of:

disp(size(I))
disp(RI)
disp(tform.T)

Where RI is the spatial referencing object for your input image, I is your input image, and tform is the geometric transformation that you are passing to imwarp.

[Hugo]

Hi Alex,

I resolved the problem with my memory already, i was still running a 32-bit version of matlab, on the 64-bit version it works fine. I didn't realize the other errors were caused by the same memory error.

Thanks a lot for your help, the syntax works fine now!

[Hugo]

Finally got the interp3 working also (used im2double to convert uint16 to double), it saves a lot of calculation time compared to imwarp, so i would like to use this one. But is there a fast way to rescale the interpolated (class double) volume to it's original intensity values (uint16)? when i just convert it using uint16 my whole volume gets the same intensity

I tried to rescale using MRIi = MRIi./max(MRIi(:)); like you suggested, but it doesn't work yet. I get the next errors:

Subscript indices must either be real positive integers or logicals.

Error using ./ Integers can only be combined with integers of the same class, or scalar doubles.

Does anyone have some ideas for me?

[Daniel Ribeiro de Sousa]

Hi,

Could you explain me why the first and the last figures are empty? I mean they are only black without any interpolation data. I should have 27*2=54 pictures but there are only 52 images.

Thank you for explanation.