How to convert 2d images to a 3d image (MRI)?

Hello everybody,

I have MRI scan file with 9 slices...how can I make a 3d image with this 9 slices?

Here the data I have (Aspect,ImageOrientation,ImagePosition):

I have used the following code and I got an error:

load data.mat
im = data(1).Img;

for i = 1 : size(im,2)
    for j = 1 : size(im,1)
        di = data(1).Aspect(1);
        dj = data(1).Aspect(2);
        Xxyz = data(1).ImageOrientation(1:3);
        Yxyz = data(1).ImageOrientation(4:6);
        Sxyz = data(1).ImagePosition;
        X(i,j) = Xxyz(1)*di*(i-1) + Yxyz(1)*dj*(j-1) + Sxyz(1);
        Y(i,j) = Xxyz(2)*di*(i-1) + Yxyz(2)*dj*(j-1) + Sxyz(2);
        Z(i,j) = Xxyz(3)*di*(i-1) + Yxyz(3)*dj*(j-1) + Sxyz(3);
    end
end
figure(1);
imReal.X = X;
imReal.Y = Y;
imReal.Z = Z;
figure(1);
surf(imReal.X', imReal.Y', imReal.Z', double(im), 'EdgeColor', 'none'); hold on;

error:

Error using surf (line 69)
Data dimensions must agree.

Error in TestTwoDtoThreeD (line 24)
surf(imReal.X', imReal.Y', imReal.Z', double(im), 'EdgeColor', 'none'); hold on;

Thanks in advance and I greatly appreciate your help, the data is attached!

Ivan

@Walter Roberson and @Image Analyst

[Oliver Woodford]

What kind of 3D plot do you want? E.g. an isocontour, or a voxel array?

[Ivan Shorokhov]

Hello Oliver,

I would like to get any of those, but prferably a voxel array.

Thank you for your response.

[Image Analyst]

What do you have? It doesn't look like parallel slices, but a random collection of slices (cutaway views) taken at a bunch of oddball locations and orientations. If that's what you have, I don't think you can reconstruct all the "missing" data into a nice 3D rectangular volumetric image, at least not easily and without some guessing. For a simple example, let's say you have 4 slices like this:

and you want to fill in the missing voxels so that you have voxel data at every point in the rectangular block. I just would not trust data there even if you could get it.

[Ivan Shorokhov]

@Image Analyst

So what I have is:

%       Orientation:        The most prominent image volume orientation
%       Aspect:             The pixel spacing in all 3 directions
%       ImageOrientation:   The direction cosines (compare DICOM standard)
%       ImagePosition:      The coordinates of the upper left corner of
%                           each image slice (compare DICOM standard)

The data is correct, according to clinicians, but in the current data I have replaced original Cardiac MRI's with single brain snapshot image, as I'm not allowed to disclose original images.

[Ivan Shorokhov]

Here is the code for solving my problem:

If anyone found this question useful please vote for it.

clc; clear all; close all;
loaded = load('data.mat');
for ni = 1:size(loaded.data,2);
    topl = 1;
    topr = 2;
    botl = 3;
    botr = 4;

      X = 1;
      Y = 2;
      Z = 3;
      fighandle = [];
      transparency = 1;

      %save the current path

          %%Read DICOM header and image data
  %     dinfo= dicominfo( dicom_filename ); %
  %     imdata = dicomread(  dicom_filename ); % loaded.data(1).Img
  %    
  imdata =  loaded.data(ni).Img;
      %%Calculate slice corner positions from the DICOM header info
      % Get the top left corner position in XYZ coordinates

      %pos    = dinfo.ImagePositionPatient; % loaded.data(ni).ImagePosition

      pos    =  loaded.data(ni).ImagePosition;

%     nc = double(dinfo.Columns);
%     nr = double(dinfo.Rows);
    nc = size(imdata,1);
    nr = size(imdata,2);

      % Get the row and column direction vercors in XYZ coordinates

      row_dircos(X) = loaded.data(ni).ImageOrientation(1);
      row_dircos(Y) = loaded.data(ni).ImageOrientation(2);
      row_dircos(Z) = loaded.data(ni).ImageOrientation(3);
      col_dircos(X) = loaded.data(ni).ImageOrientation(4);
      col_dircos(Y) = loaded.data(ni).ImageOrientation(5);
      col_dircos(Z) = loaded.data(ni).ImageOrientation(6);

      % % Check normality and orthogonality of the row and col vectors
      % Crownorm = dot(row_dircos, row_dircos);
      % Ccolnorm = dot(col_dircos, col_dircos);
      % Cdotprod = dot(row_dircos, col_dircos);
      %
      % if abs(Cdotprod) > 1e-5
      %     warnstr = sprintf('Possible dicominfo error: the dotproduct of the row and col vectors is %f should be 0',Cdotprod );
      %     disp(warnstr)
      % end

      % Calculate image dimensions
      % row_length = dinfo.PixelSpacing(1) * nr;% loaded.data(ni).Aspect(1)
      % col_length = dinfo.PixelSpacing(2) * nc;% loaded.data(ni).Aspect(2)

      row_length = loaded.data(ni).Aspect(1)  * nr;% loaded.data(ni).Aspect(1)
      col_length = loaded.data(ni).Aspect(2) * nc;% loaded.data(ni).Aspect(2)

      %%Set up the corner positions matrix in XYZ coordinates
      % Top Left Hand Corner
      corners( topl, X) = pos(X);
      corners( topl, Y) = pos(Y);
      corners( topl, Z) = pos(Z);

      % Top Right Hand Corner
      corners( topr, X) = pos(X) + row_dircos(X) * row_length;
      corners( topr, Y) = pos(Y) + row_dircos(Y) * row_length;
      corners( topr, Z) = pos(Z) + row_dircos(Z) * row_length;

      % Bottom Left Hand Corner
      corners( botl, X) = pos(X) + col_dircos(X) * col_length;
      corners( botl, Y) = pos(Y) + col_dircos(Y) * col_length;
      corners( botl, Z) = pos(Z) + col_dircos(Z) * col_length;

      % Bottom Right Hand Corner
      corners( botr, X) = pos(X) + row_dircos(X) * row_length + col_dircos(X) * col_length;
      corners( botr, Y) = pos(Y) + row_dircos(Y) * row_length + col_dircos(Y) * col_length;
      corners( botr, Z) = pos(Z) + row_dircos(Z) * row_length + col_dircos(Z) * col_length;

      %%Prepare the figure
      % Select active figure, set hold on to alow multiple slices in one figure

          colormap( gray );

      figure(1);
      hold on;

      %Tidy up the figure
      % set aspect ratio
      daspect( [1,1,1]);
      set( gca, 'color', 'none')

      %%Display slice
      %  normalize image data
      imdata = double( imdata );
      imdata = imdata / max( imdata(:));
      % scale the image
      I = imdata*255;
      % create an alternative matrix for corner points
      A( 1,1 , 1:3 ) = corners( topl, : );
      A( 1,2 , 1:3 ) = corners( topr, : );
      A( 2,1 , 1:3 ) = corners( botl, : );
      A( 2,2 , 1:3 ) = corners( botr, : );
      % extract the coordinates for the surfaces
      x = A( :,:,X );
      y = A( :,:,Y );
      z = A( :,:,Z );

      % plot surface
      surfh = surface('XData',x,'YData',y,'ZData',z,...
          'CData', I,...
          'FaceColor','texturemap',...
          'EdgeColor','none',...
          'LineStyle','none',...
          'Marker','none',...
          'MarkerFaceColor','none',...
          'MarkerEdgeColor','none',...
          'CDataMapping','direct');
      %set transparency level
      set( surfh, 'FaceAlpha', transparency );
      % label axes and optimize figure
      xlabel('RL');
      ylabel('AP');
      zlabel('FH');
       axis tight

end

      %%Optional: rotate to show all
      % edit this to create a movie
      do_movie = 0;
      show_rotate = 1;

      if do_movie
          % open avifile
          aviobj = avifile('slice3Drotate.avi');
          show_rotate = 1;
      end

      if show_rotate
          % tilt and rotate
          el = 22;
          for azplus = 0:10:360
              az = mod( 45+azplus, 360);
              view( [az, el] );
              if do_movie
                  % add farem to the movie
                  frame = getframe(fh);
                  aviobj = addframe(aviobj,frame);
              else
                  % needed to see the intemediate steps
                  drawnow;
              end
          end
      end

      if do_movie
          % close avifile
          close( aviobj );
      end

Original function taken form: http://uk.mathworks.com/matlabcentral/fileexchange/27869-mri-slice-viewer

Answer 1

[kian ghotbi]

hi, what is version of your matlab?