XTI differential phase image calculation using sinus fit

This program calculates the differential phase image from X-ray Talbot or grating interferometry. It uses sinusoidal fitting.

Here, I post simple algorithm to calculate differential phase image from M=5 foreground and M=5 background images. The attached foreground and background images are raw 16 bit unsiged integer images. The attached images were measured by X-ray Talbot inteferometry (XTI) also called grating inteferometry. The attached images corerspond to polymer polypropylene sample. In the X-ray measurement, there were phase and amplitude gratings. The amplitude grating was used for scanning. The amplitiude grating in fringe scanning technique was moved five times per grating pitch. Therefore, M parameter is equal to 5. The measured images have size 1536 x 512 pixels (H x V). In other words, 1536 columns and 512 rows. The foreground and background images are encoded as unsigned 16 bit images with little endian byte order.

Here, I use another approach to calculate differential phase image. The phase image for background is subtratced from phase image for foreground, and differential phase image is calculated. The phase image for foreground is caculated from M=5 foreground images. The number M=5 means that ampllitude grating was displaced equally 5 times with step d/5 where d is the grating pitch. The phase image for background is calculated from M=5 background images. For all M=5 foreground images, we have in single pixel 5 intensity values. These 5 intensity values in single pixel lie on sinusoidal curve. You can fit sinus function for eaxh pixel and calculate phase shift of sinusoidal curve or phase image for foreground. You can repeat the same procedure for calculation of phase background image.

The phase values or angles for foreground are stored in phase_image_fg. The phase values or angles for background are stored in phase_image_bg. The subtraction of phase_image_fg and phase_image_bg generates differential phase image (dph_image). The final differential phase image is stored as double or 64bit real number image with little endian byte order as well. The raw images attached to this contribution were measured by laboratory of Prof. Atsushi Momose from Tohoku University in Sendai, Japan. You can try to calculate absorbtion or visibility (dark field) image from the image data.

Here, the posted algorithm uses for sinusoidal fitting Levenberh-Marquardt algorithm modified for sinus function. There is levmar_sinus function in attached m-sript on the 201 line. The phase data are generated equidistantially in the lines 95 and 102 and phase values are stored in the phase vector. The sinus curve in single pixel is stored in y_data vector. The initial parameters for foreground sinusoidal fiting as well as background sinusoidal fitting is stored in x0=[x01 x02 x03] vector. The x01 is initial amplitude, x02 is intial phase and x03 is initial offset. The estimation of x01, x02, x03 from y_data vector is happening between 128 and 136 lines for foreground. The estimation of x01, x02, x03 from y_data vector is happening between 162 and 170 lines for background. The x02 parameter for foreground and background is initialized to pi/2 value. The best results were achieved here for x02=pi/2, but you can experiement with x02 value by yourself. The other parameters used in Levenberh-Marquardt algorithm are:

% maximum number of iterations

k_max=1000;

% epsilon 1

epsilon_1=1.0e-8;

% epsilon 2

epsilon_2=1.0e-8;

% tau

tau=1.0e-3;

These values were taken from document posted by Technical University of Denmark written by Madsen, Nielsen and Tingleff on the web-page <http://www2.imm.dtu.dk/pubdb/edoc/imm3215.pdf>. The maximnum number of iterations in fitting procedure is set to k_max=1000 iterations.

The fitting reults are stored in fit_result variable in the 137 and 171 lines for foreground and background, respectively.

The phase for foreground is stored in the following line

phase_image_fg(index_0,index_1)=fit_result(1,2);

and for background in the

phase_image_bg(index_0,index_1)=fit_result(1,2);

Finally, differential pahse iamge is calcualted in the line 179:

dph_image=phase_image_fg(:,:)-phase_image_bg(:,:);

I wish you good experiments with raw images for foreground and background attached to this email. You can try also to calculate absorption image or visibility (dark field image) by yourself.

The algorithm for calculation of differential phase image is using calculation on CPU. It takes almost 300 seconds to calculate differentail pahse image on the CPU. I recommend to you to do calculations on GPU using MATLAB to speed up calculations. You can do experiments with your computer as well.