Posted by Graham French on June 04, 1998 at 09:34:01:
The equation I'm trying to solve is
H[i,j] = 1/4pi * M[i,j,k,l]g[k,l]
where i,j,k,l are indices running from 1..N where N is the number of
pixels in the image.
H is the measured field at a point i,j in the detector film. g is the
magnetisation at point [k,l] in the sample being measured. M is a
geometric function.
I need to determine g. One way is to invert this matrix equation i.e
g = M^-1 x H but this requires too much memory for a large image.
Therefore I've decided to solve iterativley.
For a given sample geometry M has only to be calculated once and stored
for subsequent use in the iterations.
As you can see, the amount of memory required to do this will be
proportional to N^4. Therefore at present I am calculating M at each
point in every iteration which slows the whole process quite
considerably.
According to the papers M is a Toeplitz block Toeplitz matrix with only
N^2 DIFFERENT elements. Hence it is possible for the memory requirements
to be reduced.
the problem is that I don't know how to map the N^2 different elements
to the full N^4 array for the calculation.
Here is M in full
-arctan((sqr(a)*(k-i+0.5)*(l-j+0.5))/((d+t)*sqrt(sqr(a*(k-i+0.5))+sqr(a*
(l-j+0.5))+sqr(d+t))))
+
arctan((sqr(a)*(k-i+0.5)*(l-j+0.5))/(d*sqrt(sqr(a*(k-i+0.5))+sqr(a*(l-j+
0.5))+sqr(d))))
+
arctan((sqr(a)*(k-i+0.5)*(l-j-0.5))/((d+t)*sqrt(sqr(a*(k-i+0.5))+sqr(a*(
l-j-0.5))+sqr(d+t))))
-
arctan((sqr(a)*(k-i+0.5)*(l-j-0.5))/(d*sqrt(sqr(a*(k-i+0.5))+sqr(a*(l-j-
0.5))+sqr(d))))
+
arctan((sqr(a)*(k-i-0.5)*(l-j+0.5))/((d+t)*sqrt(sqr(a*(k-i-0.5))+sqr(a*(
l-j+0.5))+sqr(d+t))))
-
arctan((sqr(a)*(k-i-0.5)*(l-j+0.5))/(d*sqrt(sqr(a*(k-i-0.5))+sqr(a*(l-j+
0.5))+sqr(d))))
-
arctan((sqr(a)*(k-i-0.5)*(l-j-0.5))/((d+t)*sqrt(sqr(a*(k-i-0.5))+sqr(a*(
l-j-0.5))+sqr(d+t))))
+
arctan((sqr(a)*(k-i-0.5)*(l-j-0.5))/(d*sqrt(sqr(a*(k-i-0.5))+sqr(a*(l-j-
0.5))+sqr(d))))
where i,j,k,l are the indices as before and a,d,t are constants.
For N=2, M has the following results:
k=1 k=2
i=1 i=2 i=1 i=2
j=1 7.63 -9.67 | -9.67 7.63
l=1 |
j=2 -9.67 -3.01 | -3.01 -9.67
--------------------------------|--------------------
j=1 -9.67 -3.01 | -3.01 -9.67
l=2 |
j=2 7.63 -9.67 | -9.67 7.63
As you can see there are only N^2 different elements, those in the top
left sub matrix.. all other sub matrices are reflections.
I was wondering if there is a way to only store a N^2 array such as the
top left one and then map those values for different [k,l].
If it is any help, apparantly M is a symmetric semipositive definite
matrix, and the matrix equation above can be solved iteratively using
the conjugate gradient method.
I would be extremely grateful for any help you can give on this one.
Many thanks
Graham French