MATLAB: 3D interpolation Dicom data

dicominterp3?interpolation

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?

Best Answer

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(:));
Related Question