Image processing using matrix multiplication and interpolation of image data
Thomas.Dane
I am having trouble performing a correction to some 2D X-ray diffraction data I have obtained. The data is stored as a text file containing a matrix of intensity values (z-values). The camera used has 2048 x 2048 pixels (so quite a large data set!). In my experimental setup the camera is oriented at 10 degrees to the direct x-ray beam and as such there is a slight distortion to the collected data.
The correction for this data was presented by P. Boesecke (2007) (for full details see this article http://scripts.iucr.org/cgi-bin/paper?aj6013 - important operations are highlighted below) and involves applying rotation matrices to shift the pixel positions relative to the angles rotated about various axes.
<br />
R1 = 1 0 0<br />
0 cos(p1) -sin(p1)<br />
0 sin(p1) cos(p1)<br />
<br />
R2 = cos(p2) 0 sin(p2)<br />
0 1 0<br />
-sin(p2) 0 cos(p2)<br />
<br />
R3 = cos(p3) -sin(p3) 0<br />
sin(p3) cos(p3) 0<br />
0 0 1<br />
<br />
R1 = 1 0 0<br />
0 cos(p1) -sin(p1)<br />
0 sin(p1) cos(p1)<br />
<br />
R2 = cos(p2) 0 sin(p2)<br />
0 1 0<br />
-sin(p2) 0 cos(p2)<br />
<br />
R3 = cos(p3) -sin(p3) 0<br />
sin(p3) cos(p3) 0<br />
0 0 1<br />
<br />
where p = rotation angle, in my case (p1) = 10, (p2) = 0 and (p3) = 0
The pixel coordinate system is defined as a matrix:
<br />
x = f(n1)<br />
f(n2)<br />
-L<br />
x = f(n1)<br />
f(n2)<br />
-L<br />
where f(n1) and f(n2) are the physical distances on the detector in the x and y directions for each pixel (units of micrometers) from a specific pixel on the detector defined to be the centre (0) and -L is the distance from the sample to the detector (~200 milimeters).
The rotation matrix is then applied to the pixel coordinate as:
x' = R3.R2.R1.x
My problems are as follows;
Firstly, as I said the data is recorded as a 2048 x 2048 point pixel matrix and as such there is no x and y position information at the moment. I would like to know how to assign x and y values to each of the data points as these will then become the f(n1) and f(n2) values which will be shifted by the rotation matrix. In doing this I would like to define an arbitrary pixel as point 0 (for example if top left hand side of the image is (0,0) then pixel (1024,1800) should become (0,0) and all other pixels scaled accordingly) and set the units of x and y as physical distances based on the fact that a pixel has the dimensions 80 x 80 um. I don't know if this would involve concatenating every column of z data and coupling with a new two columns for x and y values or if there is a simpler solution.
The second, problem is then the matrix mathematics. I will need to compose the x matrix using the f(n1) (x data) f(n2); (y data) and -L (200 mm); then I will need to create the R rotation matrix and multiply the two. I am not certain what the result of such an operation, I assume it will be another matrix? In the paper it says "The scattering angle 2theta and the azimuth of the scattered beam can be found by expressing x' in polar coordinates." I would rather extract the data in terms of physical distance along the detector relative to the data zero as was input to the pixel coordinate matrix x than in terms of scattering angles.
Finally, the resultant output should have 3 data values per pixel; x and y coordinates which have been shifted by the rotatation matrix and intensity values. I assume the position of the x and y coordinates will no longer lie on a regular matrix grid and as such the image will need to be interpolated so that an Igor image can be generated from the data set. I have read through the Image Interpolate options and it seems that this stage would not be too difficult but any suggestions would be welcome.
Apologies for such a long request but I am really stuck! Thank you.
The first problem that you identified is handled (usually) through wave scaling. You can define your origin as well as the X and Y dimensions. To learn more about this execute the command:
DisplayHelpTopic "SetScale"
I do not have access to the reference you cited so I will proceed based on the expressions you posted which seem to represent the product of three rotations. In general, one would create the rotation matrices (3x3 double-precision waves) and multiply them together using MatrixOP, e.g.,
MatrixOP/O rotMatrix=R3 x R2 x R1
In your case with p2=p3=0 you really only have R1 to worry about (the others are identity matrices) so you don't even need to compute this product.
The form of your rotation matrix is not what I am used to so I imagine we are working with different definitions/notations. It is possible that what you need is something as simple as ImageRotate. Alternatively you could transform your data using ImageRegistration and providing the offset and rotation in the /USER=pwave format.
If you have a problem implementing any of these approaches contact support@wavemetrics.com with full details regarding the definitions of the various axes and rotations.
A.G.
WaveMetrics, Inc.
May 11, 2011 at 01:00 pm - Permalink