title: Rigidly register a 2D image to a 3D volume with good initial guess for affine transformation
body: <p>I have a 3D volume and a 2D image and an approximate mapping (affine transformation with no skwewing, known scaling, rotation and translation approximately known and need fitting) between the two. Because there is an error in this mapping and I would like to further register the 2D image to the 3D volume. I have not written code for registration purposes before, but because I can't find any programs or code to solve this I would like to try and do this. I believe the standard for registration is to optimize <a href="http://en.wikipedia.org/wiki/Mutual_information" rel="nofollow noreferrer">mutual information</a>. I think this would also be suitable here, because the intensities are not equal between the two images. So I think I should make a function for the transformation, a function for the mutual information and a function for optimization.</p>
  
  <p>I did find some Matlab code on a mathworks <a href="http://www.mathworks.com/matlabcentral/newsreader/view_thread/312687" rel="nofollow noreferrer">thread</a> from two years ago, based on an <a href="http://www.ncbi.nlm.nih.gov/pubmed/17278510" rel="nofollow noreferrer">article</a>. The OP reports that she managed to get the code to work, but I'm not getting how she did that exactly. Also in the IP package for matlab there is an <a href="http://nl.mathworks.com/help/images/examples/registering-multimodal-3-d-medical-images.html?prodcode=IP" rel="nofollow noreferrer">implementation</a>, but I dont have that package and there does not seem to be an equivalent for <a href="http://wiki.octave.org/Image_package" rel="nofollow noreferrer">octave</a>. <a href="http://www.fil.ion.ucl.ac.uk/spm/" rel="nofollow noreferrer">SPM</a> is a program that uses matlab and has registration <a href="https://www.youtube.com/watch?v=s-91ulNhtHc&amp;x-yt-ts=1421914688&amp;x-yt-cl=84503534" rel="nofollow noreferrer">implemented</a>, but does not cope with 2d to 3d registration. On the file exchange there is a brute force <a href="http://www.mathworks.com/matlabcentral/fileexchange/4145-automatic-image-registration-using--normalized--mutual-information-for-users-of-ip-toolbox" rel="nofollow noreferrer">method</a> that registers two 2D images using mutual information.</p>
  
  <p>What she does is pass a multi planar reconstruction function and an similarity/error function into a minimization algorithm. But the details I don't quite understand. Maybe it would be better to start fresh:</p>
  
  <pre><code>load mri; volume = squeeze(D);
  
  phi = 3; theta = 2; psi = 5; %some small angles
  tx = 1; ty = 1; tz = 1; % some small translation
  dx = 0.25, dy = 0.25, dz = 2; %different scales
  t = [tx; ty; tz];
  r = [phi, theta, psi]; r = r*(pi/180);
  dims = size(volume);
  p0 = [round(dims(1)/2);round(dims(2)/2);round(dims(3)/2)]; %image center
  S = eye(4); S(1,1) = dx; S(2,2) = dy; S(3,3) = dz;
  
  Rx=[1 0 0 0;
      0 cos(r(1)) sin(r(1)) 0;
      0 -sin(r(1)) cos(r(1)) 0;
      0 0 0 1];
  Ry=[cos(r(2)) 0 -sin(r(2)) 0;
      0 1 0 0;
      sin(r(2)) 0 cos(r(2)) 0;
      0 0 0 1];
  Rz=[cos(r(3)) sin(r(3)) 0 0;
      -sin(r(3)) cos(r(3)) 0 0;
      0 0 1 0;
      0 0 0 1];
  R = S*Rz*Ry*Rx;
  %make affine matrix to rotate about center of image
  T1 = ( eye(3)-R(1:3,1:3) ) * p0(1:3);
  T = T1 + t; %add translation
  A = R;
  A(1:3,4) = T;
  Rold2new = A;
  Rnew2old = inv(Rold2new);
  
  %the transformation
  [xx yy zz] = meshgrid(1:dims(1),1:dims(2),1:1);
  coordinates_axes_new = [xx(:)';yy(:)';zz(:)'; ones(size(zz(:)))'];
  coordinates_axes_old = Rnew2old*coordinates_axes_new;
  Xcoordinates = reshape(coordinates_axes_old(1,:), dims(1), dims(2), dims(3));
  Ycoordinates = reshape(coordinates_axes_old(2,:), dims(1), dims(2), dims(3));
  Zcoordinates = reshape(coordinates_axes_old(3,:), dims(1), dims(2), dims(3));
  
  %interpolation/reslicing
  method = 'cubic'; 
  slice= interp3(volume, Xcoordinates, Ycoordinates, Zcoordinates, method);
  %so now I have my slice for which I would like to find the correct position
  
  % first guess for A
  A0 = eye(4); A0(1:3,4) = T1; A0(1,1) = dx; A0(2,2) = dy; A0(3,3) = dz; 
  % this is pretty close to A
  % now how would I fit the slice to the volume by changing A0 and examining some similarity measure?
  % probably maximize mutual information?
  % http://www.mathworks.com/matlabcentral/fileexchange/14888-mutual-information-computation/content//mi/mutualinfo.m
  </code></pre>