flexible algorithms for image registration

% ======================================================================== % function By = mfBy(yc,para,flag); % (c) Jan Modersitzki 2009/04/08, see FAIR.2 and FAIRcopyright.m. % % MatrixFree B*u for elastic, diffusion, curvature % computes B*zy or B'*yc (depends on flag, works for dim=2,3) % where B is implicitly given by the parameters para. %========================================================================= function By = mfBy(yc,para,flag); % default is B*yc if ~exist('flag','var'), flag = 'By'; end; % extract parameters omega = para.omega; m = para.m; dim = length(omega)/2; h = (omega(2:2:end)-omega(1:2:end))./m; if ~isfield(para,'mu'), para.mu = 0; end; if ~isfield(para,'lambda'), para.lambda = 0; end; a = sqrt(para.mu); b = sqrt(para.mu+para.lambda); nc = prod(m); ns = prod(ones(length(m),1)*m+eye(length(m)),2); flag = [para.regularizer,'-',flag,'-',int2str(dim)]; switch flag % By for elastic-grad % yc -> [u1,u2,u3] % % dim == 3, dim == 2, % | d11 | % | d21 | % | d31 | | u1 | | d11 | % | d12 | | | | d21 | | u1 | % y= | d22 | * | u2 | y = | d12 | * | | % | d32 | | | | d22 | | u2 | % | d13 | | u3 | | d11 d22 | % | d23 | % | d33 | % | d11 d22 d33 | case 'mfElastic-By-2' y1 = reshape(yc(1:ns(1)),m(1)+1,m(2)); y2 = reshape(yc(ns(1)+(1:ns(2))),m(1),m(2)+1); p1 = nc; % 1:p1 d11 y1 p2 = p1+(m(1)+1)*(m(2)-1); % p1+1:p2 d21 y1 p3 = p2+(m(1)-1)*(m(2)+1); % p2+1:p3 d12 y2 p4 = p3+nc; % p3+1:p4 d22 y2 p5 = p4+nc; % p4+1:p5 div y By = zeros(p5,1); d1 = @(Y) reshape(Y(2:end,:)-Y(1:end-1,:),[],1)/h(1); d2 = @(Y) reshape(Y(:,2:end)-Y(:,1:end-1),[],1)/h(2); By(1:p4) = [d1(y1);d2(y1);d1(y2);d2(y2)]; % div y = d11 y1 + d22 y2 By(p4+1:p5) = By(1:p1) + By(p3+1:p4); By(1:p4) = a*By(1:p4); By(p4+1:p5) = b*By(p4+1:p5); case 'mfElastic-BTy-2' % the input yc is the result of B*y, therefore it can be splitted into 5 % parts: d11u1, d21u1,d12u2,d22u2,Divy p1 = nc; % 1:p1 d11 y1 p2 = p1+(m(1)+1)*(m(2)-1); % p1+1:p2 d21 y1 p3 = p2+(m(1)-1)*(m(2)+1); % p2+1:p3 d12 y2 p4 = p3+nc; % p3+1:p4 d22 y2 p5 = p4+nc; % p4+1:p5 div y By = zeros(sum(ns),1); d1T = @(Y) reshape([-Y(1,:);Y(1:end-1,:)-Y(2:end,:);Y(end,:)],[],1)/h(1); d2T = @(Y) reshape([-Y(:,1),Y(:,1:end-1)-Y(:,2:end),Y(:,end)],[],1)/h(2); By(1:ns(1)) = a*d1T(reshape(yc( 1:p1),m)) ... + a*d2T(reshape(yc(p1+1:p2),m(1)+1,m(2)-1)) ... + b*d1T(reshape(yc(p4+1:p5),m)); By(1+ns(1):end) = a*d1T(reshape(yc(p2+1:p3),m(1)-1,m(2)+1)) ... + a*d2T(reshape(yc(p3+1:p4),m)) ... + b*d2T(reshape(yc(p4+1:p5),m)); case 'mfElastic-By-3' y1 = reshape(yc(1:ns(1)),m(1)+1,m(2),m(3)); y2 = reshape(yc(ns(1)+(1:ns(2))),m(1),m(2)+1,m(3)); y3 = reshape(yc(ns(1)+ns(2)+(1:ns(3))),m(1),m(2),m(3)+1); p1 = nc; % 1:p1 d11 y1 p2 = p1+(m(1)+1)*(m(2)-1)*m(3); % p1+1:p2 d21 y1 p3 = p2+(m(1)+1)*m(2)*(m(3)-1); % p2+1:p3 d31 y1 p4 = p3+(m(1)-1)*(m(2)+1)*m(3); % p3+1:p4 d12 y2 p5 = p4+nc; % p4+1:p5 d22 y2 p6 = p5+m(1)*(m(2)+1)*(m(3)-1); % p5+1:p6 d32 y2 p7 = p6+(m(1)-1)*m(2)*(m(3)+1); % p6+1:p7 d13 y3 p8 = p7+m(1)*(m(2)-1)*(m(3)+1); % p7+1:p8 d23 y3 p9 = p8+nc; % p8+1:p9 d33 y3 p10 = p9+nc; % p9+1:p10 div y By = zeros(p10,1); d1 = @(Y) reshape(Y(2:end,:,:)-Y(1:end-1,:,:),[],1)/h(1); d2 = @(Y) reshape(Y(:,2:end,:)-Y(:,1:end-1,:),[],1)/h(2); d3 = @(Y) reshape(Y(:,:,2:end)-Y(:,:,1:end-1),[],1)/h(3); By(1:p9) = [... d1(y1);d2(y1);d3(y1);... d1(y2);d2(y2);d3(y2);... d1(y3);d2(y3);d3(y3) ]; % div y = d11 y1 + d22 y2 + d33 y3 By(p9+1:p10) = By(1:p1) + By(p4+1:p5)+ By(p8+1:p9); By( 1:p9 ) = a* By(1:p9); By(p9+1:p10) = b* By(p9+1:p10); case 'mfElastic-BTy-3' % the input yc is the result of B*y, therefore it can be splitted into 10 % parts: d11y1, d21y1,d31y1,d12y2,d22y2,d32y2,d13y3,d23y3,u33y3,Divy p1 = nc; % 1:p1 d11 y1 p2 = p1+(m(1)+1)*(m(2)-1)*m(3); % p1+1:p2 d21 y1 p3 = p2+(m(1)+1)*m(2)*(m(3)-1); % p2+1:p3 d31 y1 p4 = p3+(m(1)-1)*(m(2)+1)*m(3); % p3+1:p4 d12 y2 p5 = p4+nc; % p4+1:p5 d22 y2 p6 = p5+m(1)*(m(2)+1)*(m(3)-1); % p5+1:p6 d32 y2 p7 = p6+(m(1)-1)*m(2)*(m(3)+1); % p6+1:p7 d13 y3 p8 = p7+m(1)*(m(2)-1)*(m(3)+1); % p7+1:p8 d23 y3 p9 = p8+nc; % p8+1:p9 d33 y3 p10 = p9+nc; % p9+1:p10 div y By = zeros(sum(ns),1); d1T = @(Y) reshape(d1t(Y),[],1)/h(1); d2T = @(Y) reshape(d2t(Y),[],1)/h(2); d3T = @(Y) reshape(d3t(Y),[],1)/h(3); By(1:ns(1)) = a*d1T(reshape(yc( 1:p1),m)) ... + a*d2T(reshape(yc(p1+1:p2),m(1)+1,m(2)-1,m(3))) ... + a*d3T(reshape(yc(p2+1:p3),m(1)+1,m(2),m(3)-1)) ... + b*d1T(reshape(yc(p9+1:p10),m)); By(ns(1)+(1:ns(2))) = ... a*d1T(reshape(yc(p3+1:p4),m(1)-1,m(2)+1,m(3))) ... + a*d2T(reshape(yc(p4+1:p5),m)) ... + a*d3T(reshape(yc(p5+1:p6),m(1),m(2)+1,m(3)-1)) ... + b*d2T(reshape(yc(p9+1:p10),m)); By(ns(1)+ns(2)+1:end) = ... a*d1T(reshape(yc(p6+1:p7),m(1)-1,m(2),m(3)+1)) ... + a*d2T(reshape(yc(p7+1:p8),m(1),m(2)-1,m(3)+1)) ... + a*d3T(reshape(yc(p8+1:p9),m)) ... + b*d3T(reshape(yc(p9+1:p10),m)); case {'mfCurvature-By-2','mfCurvature-BTy-2'}, % % By = | \Delta y^1 0 | = | d_1^2 y^1 + d_2^2 y^1 | % | 0 \Delta y^2 | = | d_1^2 y^2 + d_2^2 y^2 | yc = reshape(yc,[m,2]); d2 = @(i) D2(i,omega,m); yc(:,:,1) = d2(1)*yc(:,:,1) + yc(:,:,1)*d2(2); yc(:,:,2) = d2(1)*yc(:,:,2) + yc(:,:,2)*d2(2); By = yc(:); case {'mfCurvature-By-3','mfCurvature-BTy-3'}, % % | \Delta y^1 0 0 | = | d_1^2 y^1 + d_2^2 y^1 + d_3^2 y^1 | % By = | 0 \Delta y^2 0 | = | d_1^2 y^2 + d_2^2 y^2 + d_3^2 y^2 | % | 0 0 \Delta y^3 | = | d_1^2 y^3 + d_2^2 y^3 + d_3^2 y^3 | % efficient implementation of Bcurvature*Y n = prod(m); % number of voxels yc = reshape(yc,[m,3]); % note that now yc(:,:,:,ell) is the ell-th % component of Y=(Y^1,Y^2,Y^3) % of size m(1)-by-m(2)-by-m(3) By = zeros(numel(yc),1); % allocate memory for the output d2 = @(i) D2(i,omega,m); % this is a shortcut to the discrete second derivative % the following line is a shortcut for % - permuting the 3d-array using the permutation J % - reshape it to a 2D-array of size q-by-prod(m)/q, where q=m(J(1)) % - multiply by A (which is q-by-q) % - undo the reshape, i.e. make the result to m(J(1))-by-m(J(2))-by-m(J(3)) % - undo the permutation operate = @(A,z,J) ipermute(reshape(A*reshape(permute(z,J),m(J(1)),[]),m(J)),J); % run over all compunents y^ell of Y=(Y^1,Y^2,Y^3) for ell=1:3, % compute % (I_3\otimes I_2\otimes d2(1) + I_3\otimes d2(2)\otimes I_1 ... % + d2(3)\otimes I_2\otimes I_1) y^ell for k=1:3, z = operate(d2(k),yc(: