Matlab: help with biot savart law

[arronslacey]

Hi, if anyone good with MATLAB and knows biot savart law then i hope you can help. I have the following program:

 clear all
 Img = imread('littlecircle.png');
 Img = Img(:,:,1);
 Img = double(Img);
 w = size(Img,1);               % width size
 h = size(Img,2);               % height size
 [Ix,Iy] = gradient(Img);       %gradient of image
 
 i=1;     %iteration for magnetic field loop
 b=0;     %initialize b to zero
 
 % Magnetic Field
 for pxRow = 1:h % fixed pixel row
 for pxCol = 1:w % fixed pixel column
 
 for r = 1:h % row of distant pixel
 for c = 1:w % column of distant pixel
 
 R(c,r) = sqrt((r-pxRow)^2 + (c-pxCol)^2);                               % pythagoras theorem to get distance to each pixel
 O(c,r) = atan(Iy(c,r)./Ix(c,r));                         % direction of current
 If(c,r) = sqrt((Ix(c,r)).^2 + (Iy(c,r)).^2);                            % magnitude of current 
 Rxs(c,r) = R(c,r)./norm(R(c,r));                                        % unit vector from x to s                     
 b(c,r) = If(c,r).*O(c,r).*(Rxs(c,r)./(R(c,r)).^2);                      % b field = If(s)O(s) x Rxs/R.^2  BIOT SAVART LAWRxs/R.^2  BIOT SAVART LAW
 end
 end
 B(i) = {b}; % filling a cell array with results. read below
 i = i+1;
 end
 end

all lines in red are currently not calculating what I want them to. Can anyone help? Thanks!

 

[arronslacey]

by the way, here is the reference paper http://www.cs.swan.ac.uk/~csjason/papers/xxmm-pami2008.pdf

 

[gsal]

Debugging hints/help:

They are not calculating what you want them to?
What ARE they calculating? and
what do YOU expect them to be?

Are you sure your starting numbers are correct (Ix,Iy) ?

What does the function 'norm' do for you?
You are taking the norm of scalar R(c,r)...what is that? Its length? Its magnitude? Is it always a positive number?
By the way, you just calculated R(c,r) above as the sqrt() of something...so, R(c,r) IS a scalar and IS positive; also, it is just a scalar...not a matrix.
So...what do you expect Rxs(c,r) = R(c,r)./norm(R(c,r)) to be? You call it a unit vector...but I don't think it is a vector at all...it is just a number.

It seems to me that Rxs(c,r) will always be 1.0 ?

 

[arronslacey]

I am struggling with getting to grips with vectors in MATLAB i'll admit. but here's what i need to do:

1. read in an image
2. compute image gradient vectors [Ix,Iy], and then magnitude sqrt(Ix.^2+Iy.^2) = If(c,r)
3. compute the edge orientation of the image by rotating the gradient vectors 90 degrees clockwise or anti clockwise (O(c,r))
4. take each pixel in the image, and find its distance to all other pixels (R(c,r))
5. compute the unit vectors of the pixel distances Rxs(c,r)
6. Compute b field with Biot-Savart law B = If.*O cross (Rxs/R.^2)

and yes my understanding of how MATLAB handles vectors isn't helping the issue.

 

[gsal]

I hear you...but you did not answer ANY of my questions. :-(

 

[arronslacey]

sorry i'll have a go -

only R(c,r) calculates what I want. it calculates the distance one pixel to all others and stores into an array.

O(c,r) is supposed to be the current direction, ie taking the gradient vectors Ix(c,r) and Iy(c,r) and rotating them by 90 degrees. When I quiver plot these rotated vectors it looks as expected, however I am told that the actual result of these are just scalars, ie numbers - but they need to be vectors.

Rxs(c,r) is supposed to be the unit vector of R(c,r). Matlab computes unit vectors with the norm function:

Rxs(c,r) = R(c,r)/norm(R(c,r)) give me the unit vectors, and they are indeed all supposed = 1

Where I could be going wrong, as you say is that none of the above are vectors in which I have a problem. I am new to MATLAB and when you do your sums on the blackboard you take the nature of vectors for granted and just punch them in the calculator. but I'm having difficulty representing them here in matlab.

The clue "you are taking the norm of a scalar R(c,r)" as you say...well how do I turn R(c,r) into a vector? If i knew this I am sure the rest of the program would fall into place. I have been relying on matlabs quiver function (vector plot) to show me the direction of these numbers, and they do indeed have direction, hence why I'm assuming I'm returning vectors...argh!

 

[gsal]

Yeap...you are struggling.

When on the blackboard you take the nature of vectors granted?
What do vectors look like?

Recall that a vector has a magnitude AND a direction (angle).

So, step back a bit and take a good look at your own code as if you had not written it...just because you had something in mind and had certain intent, that does not mean that it went through...take a look at your calculations and you will see that you are just doing a bunch of scalar operations...nothing there is a vector...

That is not necessarily bad, i.e., keeping track of your own magnitudes and your own angles separately is o.k...if you know that you are doing it!.

Other than that, you say that O(c,r) is a current direction (angle)? I guess it is, after all, it is the arctan() of something; so it is an angle, but just an angle, not a vector...are these the angles of your vectors? Where are the corresponding magnitudes? If()?

And Rxs() is not giving you unit vectors, it is just giving you units!



Alternatively...if you want to take vectors for granted in matlab, than you first need to find out how to handle vectors in matlab...really, vectors, you know, those things that do have magnitude and angle...can MATLAB handle this quantity in a single entity, variable? I think it might but I have not used MATLAB in about 20 years...so, take a few well invested minutes and do a google search or read the user's manual...I guess I could do it, too; but it is best if you learn this by yourself.