Hard-register a 2D image into a 3D volume with a good initial premise for affine transformation

Question

I have a 3D volume and a 2D image, as well as an approximate display (affine transformation without twisting, known scaling, rotation and translation, approximately known and need to be adjusted) between them. Because there is an error in this display, and I would like to additionally register the 2D image in the 3D volume. I didn’t write code for registration before, but since I can’t find any programs or code to solve this problem, I would like to try and do it. I believe that the registration standard is the optimization of mutual information . I think this would also be suitable here, because the intensities are not equal between the two images. Therefore, I think that I should make a function for conversion, a function for mutual information, and a function for optimization.

I found Matlab code in mathworks thread from two years ago based on an article . The OP reports that she managed to get the code to work, but I don’t understand how she did it. There is also an implementation in the IP package for Matlab, but I do not have this package, and it does not seem to be equivalent for octave . SPM is a program that uses Matlab and has registration implemented , but does not cope with the registration of 2d and 3d. In file sharing, there is a brute force method that registers two 2D images using mutual information.

What she does is pass the multifaceted recovery function and the similarity / error function to the minimization algorithm. But I don’t quite understand the details. It might be better to start a new life:

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 

Answer 1

Well, I was hoping someone else would do, it would probably be better than mine, as I had never done any optimization or registration before. So I waited until the Knedlsepps bouncer almost finished. But I have some code that works now. He will find a local optimum, so the initial guess should be good. I wrote some functions myself, took some functions from file sharing as it is, and I widely edited some other functions from file sharing. Now, when I give all the code together to work as an example, the rotation is disabled, it will try to fix it. I'm not sure where the difference in code between the example and my source code was to make a typo when replacing some variables and the data loading scheme.

What I do, I take the original matrix of affine transformations, decompose it into an orthogonal matrix and an upper triangular matrix. Then I assume that the orthogonal matrix is ​​my rotation matrix, so I calculate Euler angles. I directly take the translation from the affine matrix, and, as indicated in the problem, I assume that I know the scaling matrix and there is no shift. So, I have all degrees of freedom for the affine transformation, which my optimization function changes and builds a new affine matrix, applies it to that and calculates mutual information. The matlab patternsearch optimization function then minimizes 1-MI / MI_max.

What I noticed when using it on my real data, which are multimodal images of the brain, is that it works much better on images extracted from the brain, therefore it is removed from the skull and tissue outside the skull.

 %data load mri; volume = double(squeeze(D)); %transformation parameters phi = 3; theta = 1; psi = 5; %some small angles tx = 1; ty = 1; tz = 3; % 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); %image center and size dims = size(volume); p0 = [round(dims(1)/2);round(dims(2)/2);round(dims(3)/2)]; %slice coordinate ranges range_x = 1:dims(1)/dx; range_y = 1:dims(2)/dy; range_z = 1; %rotation R = dofaffine([0;0;0], r, [1,1,1]); T1 = ( eye(3)-R(1:3,1:3) ) * p0(1:3); %rotate about p0 %scaling S = eye(4); S(1,1) = dx; S(2,2) = dy; S(3,3) = dz; %translation T = [[eye(3), T1 + t]; [0 0 0 1]]; %affine A = T*R*S; % first guess for A r00 = [1,1,1]*pi/180; R00 = dofaffine([0;0;0], r00, [1 1 1]); t00 = T1 + t + ( eye(3) - R00(1:3,1:3) ) * p0; A0 = dofaffine( t00, r00, [dx, dy, dz] ); [ t0, r0, s0 ] = dofaffine( A0 ); x0 = [ t0.', r0, s0 ]; %the transformation slice = affine3d(volume, A, range_x, range_y, range_z, 'cubic'); guess = affine3d(volume, A0, range_x, range_y, range_z, 'cubic'); %initialisation Dt = [1; 1; 1]; Dr = [2 2 2].*pi/180; Ds = [0 0 0]; Dx = [Dt', Dr, Ds]; %limits LB = x0-Dx; UB = x0+Dx; %other inputs ref_levels = length(unique(slice)); Qref = imquantize(slice, ref_levels); MI_max = MI_GG(Qref, Qref); %patternsearch options options = psoptimset('InitialMeshSize',0.03,'MaxIter',20,'TolCon',1e-5,'TolMesh',5e-5,'TolX',1e-6,'PlotFcns',{@psplotbestf,@psplotbestx}); %optimise [x2, MI_norm_neg, exitflag_len] = patternsearch(@(x) AffRegOptFunc(x, slice, volume, MI_max, x0), x0,[],[],[],[],LB(:),UB(:),options); %check p0 = [round(size(volume)/2).']; R0 = dofaffine([0;0;0], x2(4:6)-x0(4:6), [1 1 1]); t1 = ( eye(3) - R0(1:3,1:3) ) * p0; A2 = dofaffine( x2(1:3).'+t1, x2(4:6), x2(7:9) ) ; fitted = affine3d(volume, A2, range_x, range_y, range_z, 'cubic'); overlay1 = imfuse(slice, guess); overlay2 = imfuse(slice, fitted); figure(101); ax(1) = subplot(1,2,1); imshow(overlay1, []); title('pre-reg') ax(2) = subplot(1,2,2); imshow(overlay2, []); title('post-reg'); linkaxes(ax); function normed_score = AffRegOptFunc( x, ref_im, reg_im, MI_max, x0 ) t = x(1:3).'; r = x(4:6); s = x(7:9); rangx = 1:size(ref_im,1); rangy = 1:size(ref_im,2); rangz = 1:size(ref_im,3); ref_levels = length(unique(ref_im)); reg_levels = length(unique(reg_im)); t0 = x0(1:3).'; r0 = x0(4:6); s0 = x0(7:9); p0 = [round(size(reg_im)/2).']; R = dofaffine([0;0;0], r-r0, [1 1 1]); t1 = ( eye(3) - R(1:3,1:3) ) * p0; t = t + t1; Ap = dofaffine( t, r, s ); reg_im_t = affine3d(reg_im, A, rangx, rangy, rangz, 'cubic'); Qref = imquantize(ref_im, ref_levels); Qreg = imquantize(reg_im_t, reg_levels); MI = MI_GG(Qref, Qreg); normed_score = 1-MI/MI_max; end function [ varargout ] = dofaffine( varargin ) % [ t, r, s ] = dofaffine( A ) % [ A ] = dofaffine( t, r, s ) if nargin == 1 %affine to degrees of freedom (no shear) A = varargin{1}; [T, R, S] = decompaffine(A); r = GetEulerAngles(R(1:3,1:3)); s = [S(1,1), S(2,2), S(3,3)]; t = T(1:3,4); varargout{1} = t; varargout{2} = r; varargout{3} = s; elseif nargin == 3 %degrees of freedom to affine (no shear) t = varargin{1}; r = varargin{2}; s = varargin{3}; R = GetEulerAngles(r); R(4,4) = 1; S(1,1) = s(1); S(2,2) = s(2); S(3,3) = s(3); S(4,4) = 1; T = eye(4); T(1,4) = t(1); T(2,4) = t(2); T(3,4) = t(3); A = T*R*S; varargout{1} = A; else error('incorrect number of input arguments'); end end function [ T, R, S ] = decompaffine( A ) %I assume A = T * R * S T = eye(4); R = eye(4); S = eye(4); %decompose in orthogonal matrix q and upper triangular matrix r %I assume q is a rotation matrix and r is a scale and shear matrix %matlab 2014 can force real solution [qr] = qr(A(1:3,1:3)); R(1:3,1:3) = q; S(1:3,1:3) = r; % A*S^-1*R^-1 = T*R*S*S^-1*R^-1 = T*R*I*R^-1 = T*R*R^-1 = T*I = T T = A*inv(S)*inv(R); t = T(1:3,4); T = [eye(4) + [[0 0 0;0 0 0;0 0 0;0 0 0],[t;0]]]; end function [varargout]= GetEulerAngles(R) assert(length(R)==3) dims = size(R); if min(dims)==1 rx = R(1); ry = R(2); rz = R(3); R = [[ cos(ry)*cos(rz), -cos(ry)*sin(rz), sin(ry)];... [ cos(rx)*sin(rz) + cos(rz)*sin(rx)*sin(ry), cos(rx)*cos(rz) - sin(rx)*sin(ry)*sin(rz), -cos(ry)*sin(rx)];... [ sin(rx)*sin(rz) - cos(rx)*cos(rz)*sin(ry), cos(rz)*sin(rx) + cos(rx)*sin(ry)*sin(rz), cos(rx)*cos(ry)]]; varargout{1} = R; else ry=asin(R(1,3)); rz=acos(R(1,1)/cos(ry)); rx=acos(R(3,3)/cos(ry)); if nargout > 1 && nargout < 4 varargout{1} = rx; varargout{2} = ry; varargout{3} = rz; elseif nargout == 1 varargout{1} = [rx ry rz]; else error('wrong number of output arguments'); end end end