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)*(ki+0.5)*(lj+0.5))/((d+t)*sqrt(sqr(a*(ki+0.5))+sqr(a*
(lj+0.5))+sqr(d+t))))
+
arctan((sqr(a)*(ki+0.5)*(lj+0.5))/(d*sqrt(sqr(a*(ki+0.5))+sqr(a*(lj+
0.5))+sqr(d))))
+
arctan((sqr(a)*(ki+0.5)*(lj0.5))/((d+t)*sqrt(sqr(a*(ki+0.5))+sqr(a*(
lj0.5))+sqr(d+t))))

arctan((sqr(a)*(ki+0.5)*(lj0.5))/(d*sqrt(sqr(a*(ki+0.5))+sqr(a*(lj
0.5))+sqr(d))))
+
arctan((sqr(a)*(ki0.5)*(lj+0.5))/((d+t)*sqrt(sqr(a*(ki0.5))+sqr(a*(
lj+0.5))+sqr(d+t))))

arctan((sqr(a)*(ki0.5)*(lj+0.5))/(d*sqrt(sqr(a*(ki0.5))+sqr(a*(lj+
0.5))+sqr(d))))

arctan((sqr(a)*(ki0.5)*(lj0.5))/((d+t)*sqrt(sqr(a*(ki0.5))+sqr(a*(
lj0.5))+sqr(d+t))))
+
arctan((sqr(a)*(ki0.5)*(lj0.5))/(d*sqrt(sqr(a*(ki0.5))+sqr(a*(lj
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