- 10

- 0

Code:

```
function J=ImageDerivatives2D(I,sigma,type)
% Gaussian based image derivatives
%
% J=ImageDerivatives2D(I,sigma,type)
%
% inputs,
% I : The 2D image
% sigma : Gaussian Sigma
% type : 'x', 'y', 'xx', 'xy', 'yy'
%
% outputs,
% J : The image derivative
%
% Make derivatives kernels
[x,y]=ndgrid(floor(-3*sigma):ceil(3*sigma),floor(-3*sigma):ceil(3*sigma));
switch(type)
case 'x'
DGauss=-(x./(2*pi*sigma^4)).*exp(-(x.^2+y.^2)/(2*sigma^2));
case 'y'
DGauss=-(y./(2*pi*sigma^4)).*exp(-(x.^2+y.^2)/(2*sigma^2));
case 'xx'
DGauss = 1/(2*pi*sigma^4) * (x.^2/sigma^2 - 1) .* exp(-(x.^2 + y.^2)/(2*sigma^2));
case {'xy','yx'}
DGauss = 1/(2*pi*sigma^6) * (x .* y) .* exp(-(x.^2 + y.^2)/(2*sigma^2));
case 'yy'
DGauss = 1/(2*pi*sigma^4) * (y.^2/sigma^2 - 1) .* exp(-(x.^2 + y.^2)/(2*sigma^2));
end
J = conv2(I,DGauss,'same');
```

Code:

```
% Squared magnitude of force field
Fx= Fext(:,:,1);
Fy= Fext(:,:,2);
% Calculate magnitude
sMag = Fx.^2+ Fy.^2;
% Set new vector-field to initial field
u=Fx; v=Fy;
for i=1:Iterations,
% First order image derivatives
Uxx=ImageDerivatives2D(u,Sigma,'x');
Uyy=ImageDerivatives2D(u,Sigma,'y');
Vxx=ImageDerivatives2D(v,Sigma,'x');
Vyy=ImageDerivatives2D(v,Sigma,'y');
% Compute curl and update vector field
u = u + cross(Uxx,Uyy,3);
v = v + cross(Vxx,Vyy,3);
end
Fext(:,:,1) = u;
Fext(:,:,2) = v;
```