Curl of Gaussian image derivatives

Hi guys, I am trying to create a magnetic field from an image contour. I have a function which takes the image derivative via a Guassian:

    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');

and now i need to take the curl of those image derivatives to get the pseudo magnetic field. having a bit of trouble with that though:

   % 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;

anyone have any tips as the second script does not work. i can to the divergence and laplacian easily as they are just first and second derivatives, but the curl is proving a bit of a problem. Any help would be great thanks!

[Andrew Newell]

How exactly does it not work?

Answer 1

[arron lacey]

Hi Andrew - I get the following error:

??? Error using ==> cross at 43 A and B must be of length 3 in the dimension in which the cross product is taken.

Answer 2

[David Young]

Note that if you take the cross-product of 2 vectors in the plane, only the Z-component of the result is non-zero.

The following two suggestions are both untested.

Use cross:

U = cat(3, Uxx, Uyy, zeros(size(Uxx)));
V = cat(3, Vxx, Vyy, zeros(size(Vxx)));
crs = cross(U, V, 3);
crsZ = crs(:, :, 3);

Directly from the formula:

crsZ = Uxx .* Vyy - Uyy .* Vxx;

(but do check I've got the sign right).