Rigid ein 2D-Bild zu einem 3D-Volumen mit guter Anfangsrate für affine Transformation

Question

Ich habe ein 3D-Volumen und ein 2D-Bild und eine ungefähre Zuordnung (affine Transformation ohne Skwewing, bekannte Skalierung, Rotation und Übersetzung etwa bekannt und muss passen) zwischen den beiden. Weil bei diesem Mapping ein Fehler aufgetreten ist und ich das 2D-Bild weiter auf das 3D-Volumen registrieren möchte. Ich habe noch keinen Code für Registrierungszwecke geschrieben, aber weil ich keine Programme oder Codes finden kann, um das zu lösen, möchte ich versuchen, dies zu tun. Ich glaube, der Standard für die Registrierung ist zu optimieren mutual information. Ich denke, das wäre auch hier passend, weil die Intensitäten der beiden Bilder nicht gleich sind. Also ich denke, ich sollte eine Funktion für die Transformation, eine Funktion für die gegenseitige Information und eine Funktion für die Optimierung machen.

Ich habe Matlab-Code auf einem Mathworks thread vor zwei Jahren gefunden, basierend auf einem article. Die OP berichtet, dass sie es geschafft hat, dass der Code funktioniert, aber ich verstehe nicht, wie sie das genau gemacht hat. Auch im IP-Paket für Matlab gibt es eine implementation, aber ich habe dieses Paket nicht und es scheint kein Äquivalent für octave zu sein. SPM ist ein Programm, das Matlab verwendet und hat Registrierung implemented, aber nicht mit 2D-Registrierung zu bewältigen. Auf dem Dateiaustausch gibt es eine Brute-Force method, die zwei 2D-Bilder mit gegenseitigen Informationen registriert.

Was sie tut, besteht darin, eine multiplanare Rekonstruktionsfunktion und eine Ähnlichkeits-/Fehlerfunktion in einen Minimierungsalgorithmus zu überführen. Aber die Details verstehe ich nicht ganz. Vielleicht wäre es besser, neu zu beginnen:

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

Ok ich hatte gehofft, dass Ansatz jemand anderes, das ist wahrscheinlich besser gewesen wäre, als ich, wie ich noch nie irgendeine Optimierung oder Registrierung getan haben. Also wartete ich darauf, dass Knedlsepps Kopfgeld fast zu Ende ging. Aber ich habe einen Code, der jetzt funktioniert. Es wird ein lokales Optimum finden, so dass die anfängliche Schätzung gut sein muss. Ich habe selbst einige Funktionen geschrieben, einige Funktionen aus dem Dateiaustausch übernommen und einige andere Funktionen aus dem Dateiaustausch ausgiebig bearbeitet. Jetzt, da ich den ganzen Code zusammenarbeite, um hier als Beispiel zu arbeiten, sind die Rotationen ausgeschaltet, werden versuchen, das zu korrigieren. Ich bin mir nicht sicher, wo der Unterschied im Code zwischen dem Beispiel und meinem ursprünglichen Code ist, muss einen Tippfehler beim Ersetzen einiger Variablen und Datenladeschemas gemacht haben.

Was ich mache ist, ich nehme die anfängliche affine Transformationsmatrix, zerlege sie in eine orthogonale Matrix und eine obere Dreiecksmatrix. Ich nehme dann an, dass die orthogonale Matrix meine Rotationsmatrix ist, also berechne ich daraus die Euler-Winkel. Ich nehme die Übersetzung direkt aus der affinen Matrix, und wie in der Aufgabe angegeben, nehme ich an, dass ich die Skalierungsmatrix kenne und dass es keine Scherung gibt. So habe ich alle Freiheitsgrade für die affine Transformation, die meine Optimierungsfunktion verändert und daraus eine neue affine Matrix konstruiert, auf das Volumen anwendet und die gegenseitige Information berechnet. Die Matlab-Optimierungsfunktion Mustersuche minimiert dann 1-MI/MI_max.

Was ich bei der Verwendung meiner realen Daten, die multimodale Bilder des Gehirns sind, bemerkte, ist, dass es viel besser auf Gehirn extrahiert Bilder funktioniert, so mit dem Schädel und Gewebe außerhalb des Schädels entfernt.

%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 
    [q r] = 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