Problem in calculating the derivative of a cost function

1 view (last 30 days)
Felippe Trigueiro Angelo
Felippe Trigueiro Angelo on 2 Sep 2019
Answered: darova on 3 Sep 2019
Hi, I'm having a problem with the calculation of the following cost function:
, where the vector T of displacement is calculated according the parametric function B-Splines (see below) and F and R are 3D images :
.
I have to calculate the derivative of the SSD function according with the parameters phi. The expression for this derivative is given by:
, where the derivatives inside of the function are given by:
Sorry for the amount of mathematical expressions.
My problem is: The expression of the derivative 1 should give me a vector with the number of parameters of phi, right? But my implementation is not returning that number of parameters (see below). So, what am I doing wrong?
[Gx, Gy, Gz] = imgradientxyz(ft, 'central');
Dt_Dtheta = g_transformation3_gradient(ft, delta);
Dssd_Dphi = cat(4, Gx.*2.*(ft - fr).*Dt_Dtheta, Gy.*2.*(ft - fr).*Dt_Dtheta, ...
Gz.*2.*(ft - fr).*Dt_Dtheta);
function grad = g_transformation3_gradient(f, delta)
%f --> image to be transformed
%c --> control points matrix
%delta --> spacing in pixels, between the control points
[X, Y, Z] = size(f);
x0 = 0:(X - 1);
y0 = 0:(Y - 1);
z0 = 0:(Z - 1);
u = x0./delta - floor(x0./delta);
v = y0./delta - floor(y0./delta);
w = z0./delta - floor(z0./delta);
Bu = zeros(4, numel(u));
Bv = zeros(4, numel(v));
Bw = zeros(4, numel(w));
%Computing the vector of B-splines
Bu(1, :) = ((1 - u).^3)/6;
Bu(2, :) = (3*u.^3 - 6*u.^2 + 4)/6;
Bu(3, :) = (-3*u.^3 + 3*u.^2 + 3*u + 1)/6;
Bu(4, :) = (u.^3)/6;
Bv(1, :) = ((1 - v).^3)/6;
Bv(2, :) = (3*v.^3 - 6*v.^2 + 4)/6;
Bv(3, :) = (-3*v.^3 + 3*v.^2 + 3*v + 1)/6;
Bv(4, :) = (v.^3)/6;
Bw(1, :) = ((1 - w).^3)/6;
Bw(2, :) = (3*w.^3 - 6*w.^2 + 4)/6;
Bw(3, :) = (-3*w.^3 + 3*w.^2 + 3*w + 1)/6;
Bw(4, :) = (w.^3)/6;
T0 = Bu(1, :).*(Bv(1, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])) + ...
Bv(2, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])) +...
Bv(3, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])) + ...
Bv(4, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])));
%matr = reshape(Bu(1, :), [1 1 Z]);
%T0 = Bu(1, :).*(Tw0);
%clear Tw0 ;
T1 = Bu(2, :).*(Bv(1, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])) + ...
Bv(2, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])) + ...
Bv(3, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])) + ...
Bv(4, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])));
%matr = reshape(Bu(2, :), [1 1 Z]);
%T1 = Bu(2, :).*(Tw1);
%clear Tw1;
T2 = Bu(3, :).*(Bv(1, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])) + ...
Bv(2, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])) + ...
Bv(3, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])) + ...
Bv(4, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])));
%matr = reshape(Bu(3, :), [1 1 Z]);
%T2 = Bu(3, :).*(Tw2);
%clear Tw2;
T3 = Bu(4, :).*(Bv(1, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])) + ...
Bv(2, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])) + ...
Bv(3, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])) + ...
Bv(4, :)'.*(reshape(Bw(1, :), [1 1 Z]) + ...
reshape(Bw(2, :), [1 1 Z]) + ...
reshape(Bw(3, :), [1 1 Z]) + ...
reshape(Bw(4, :), [1 1 Z])));
%matr = reshape(Bu(4, :), [1 1 Z]);
%T3 = Bu(4, :).*(Tw3);
%clear Tw3;
grad = T0 + T1 + T2 + T3;

Sign in to answer this question.

Accepted Answer

darova
darova on 3 Sep 2019
Try this
[U,V,W] = meshgrid(u,v,w);
%Computing the vector of B-splines
Bu(:,:,:,1) = ((1 - U).^3)/6;
Bu(:,:,:,2) = (3*U.^3 - 6*U.^2 + 4)/6;
Bu(:,:,:,3) = (-3*U.^3 + 3*U.^2 + 3*U + 1)/6;
Bu(:,:,:,4) = (U.^3)/6;
Bv(:,:,:,1) = ((1 - V).^3)/6;
Bv(:,:,:,2) = (3*V.^3 - 6*V.^2 + 4)/6;
Bv(:,:,:,3) = (-3*V.^3 + 3*V.^2 + 3*V + 1)/6;
Bv(:,:,:,4) = (V.^3)/6;
Bw(:,:,:,1) = ((1 - W).^3)/6;
Bw(:,:,:,2) = (3*W.^3 - 6*W.^2 + 4)/6;
Bw(:,:,:,3) = (-3*W.^3 + 3*W.^2 + 3*W + 1)/6;
Bw(:,:,:,4) = (W.^3)/6;
Tgrad = zeros(size(U));
for i = 1:4
for j = 1:4
for k = 1:4
Tgrad = Tgrad + Bu(:,:,:,i).*Bv(:,:,:,j).*Bw(:,:,:,k);
end
end
end

More Answers (0)

Categories

Find more on Geometric Transformation and Image Registration in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!