Commit aeda86a4 authored by krempasky's avatar krempasky
Browse files

tar files removed, folders created

parent 664994cd
CC = g++
DOUBLE_PRECISION = defined
ifdef DOUBLE_PRECISION
PREC = -DUSE_DOUBLE_PRECISION
FFTW = -L /home/optics/Software/fftw/lib -lfftw3
else
PREC =
FFTW = -lfftw3f
endif
CFLAGS = -O3 -std=c++0x $(PREC) -DHAVE_TIFF -DHAVE_HDF5 -I/usr/include/hdf5/serial -I/home/optics/Software/fftw/include
LDFLAGS = -lm -ltiff -lpng -lhdf5 $(FFTW)
src = $(wildcard *.cpp)
obj = $(src:.cpp=.o)
all: $(obj)
$(CC) -o main $(obj) $(LDFLAGS)
%.o: %.cpp
$(CC) $(CFLAGS) -c -o $@ $<
.PHONY: clean
clean:
rm $(obj)
figure(102)
subplot(2,4,1);
imagesc(imread('roi.tif'));colormap(flipud(bone)); title('moire image')%freezeColors;
colormap default
ax_ft=subplot(2,4,2); plot(log(fftshift(imread('fft.tif'))));set(ax_ft, 'YTick','');title('1D FFT')
%subplot_tight(2,4,3,.06); imagesc(amp);
cifimg_re = imread('ifftre.tif');
cifimg_im = imread('ifftim.tif');
subplot(2,4,3); imagesc(abs(cifimg_re + 1j*cifimg_im));title('1st order FFT filter')
cwrap = imread('dphi_w.tif');
subplot(2,4,4); imagesc(cwrap);title('wrapped fringe phase');
%subplot(2,4,5); imagesc(unwrap);
cwrap = imread('alpha2D.tif');
subplot(2,4,5); imagesc(cwrap);title('fringe deflection \alpha');
cpolys = imread('polys.tif').';
%ax1= subplot(2,4,6); plot(cpolys(1,:));title('\delta\alpha/\delta x')%xlim(ax1,[0 172])
cwpa = imread('wpa.tif');
ax2= subplot(2,4,6); plot(cwpa);title('1D \alpha'), ylabel('radians')%xlim(ax2,[0 0.0005])
%subplot(2,4,7); imagesc(out.alpha2D);
%ax2= subplot(2,4,8); plot(out.hp*1E9);title('height profile');ylabel('nanometers')%xlim(ax2,[0 0.0005])
cdphi = imread('hp.tif');
subplot(2,4,7); plot(cdphi*1E9);title('hp');
cphiy = imread('phi_y.tif');
ax2 = subplot(2,4,8); plot(cphiy*1E9);title('spherical wavefront');ylabel('nanometers')
This diff is collapsed.
#ifndef XGI_H
#define XGI_H
#ifdef USE_DOUBLE_PRECISION
#define FLOATING double
#define C_ONE 1.0
#define C_EPSILON 1e-11
#define C_PI 3.14159265358979323846
#define FMAX fmax
#define FMIN fmin
#define ROUND round
#define ABS fabs
#define COS cos
#define ATAN2 atan2
#define SQRT sqrt
#define FFTPLAN fftw_plan
#define FFTPLAN1D fftw_plan_dft_1d
#define FFTPLANMANY fftw_plan_many_dft
#define FFTPLAND fftw_destroy_plan
#define FFTC fftw_complex
#define FFTMALLOC fftw_malloc
#define FFTFREE fftw_free
#define FFTEXEC fftw_execute
#define FFTEXECA fftw_execute_dft
#define FFTTINIT fftw_init_threads
#define FFTTPLAN fftw_plan_with_nthreads
#define FFTTCLEANUP fftw_cleanup_threads
#define FLOATING_SF "%lf"
#else
#define FLOATING float
#define C_ONE 1.0f
#define C_EPSILON 1e-5
#define C_PI 3.1415926f
#define FMAX fmaxf
#define FMIN fminf
#define ROUND roundf
#define ABS fabsf
#define COS cosf
#define ATAN2 atan2f
#define SQRT sqrtf
#define FFTPLAN fftwf_plan
#define FFTPLAN1D fftwf_plan_dft_1d
#define FFTPLANMANY fftwf_plan_many_dft
#define FFTPLAND fftwf_destroy_plan
#define FFTC fftwf_complex
#define FFTMALLOC fftwf_malloc
#define FFTFREE fftwf_free
#define FFTEXEC fftwf_execute
#define FFTEXECA fftwf_execute_dft
#define FFTTINIT fftwf_init_threads
#define FFTTPLAN fftwf_plan_with_nthreads
#define FFTTCLEANUP fftwf_cleanup_threads
#define FLOATING_SF "%f"
#endif
#define SOF sizeof(FLOATING)
#define C_DISCONT (C_ONE*C_PI)
#include <math.h>
struct Rect
{
int x_lo;
int x_hi;
int y_lo;
int y_hi;
Rect(int _x_lo,int _x_hi,int _y_lo, int _y_hi):x_lo(_x_lo),x_hi(_x_hi),y_lo(_y_lo),y_hi(_y_hi){};
Rect():x_lo(0),x_hi(0),y_lo(0),y_hi(0){};
};
struct Orientation
{
char tpose;
char fliplr;
char flipud;
Orientation(char _tpose,char _fliplr,char _flipud):tpose(_tpose),fliplr(_fliplr),flipud(_flipud){};
};
struct Matrix
{
int width;
int height;
FLOATING *data;
Matrix();
Matrix(int w,int h);
~Matrix();
void free();
int resize(int w,int h);
int copy(Matrix *tgt);
int zero(int w,int h);
int add(Matrix *other);
int sub(Matrix *other);
int mul(Matrix *other);
int mul(FLOATING factor);
int crop(Matrix *tgt,Rect roi, Orientation orientation = Orientation(0,0,0));
int unwrap(Matrix *tgt);
int save_tiff(const char *fn);
int save_png(const char *fn,FLOATING data_low = INFINITY, FLOATING data_high = - INFINITY);
int save_txt(const char *fn);
int load_txt(const char *fn);
void print();
int load_h5_ushort(const char *fn,const char *ds, int slice=0);
static int add(Matrix *a, Matrix *b, Matrix *c);
static int sub(Matrix *a, Matrix *b, Matrix *c);
static int mul(Matrix *a, Matrix *b, Matrix *c);
static int mul(Matrix *a, FLOATING b, Matrix *c);
};
struct Xgi_ip
{
//roi pixel positions: <lo,hi)
Rect roi;
Orientation cam_orientation;
FLOATING beta_fix;
FLOATING beta_rot;
//period of grating 1 [m]
FLOATING period_g1;
//period of grating 2 [m]
FLOATING period_g2;
FLOATING distance;
int talbot_order;
FLOATING energy;
FLOATING pixel_size;
//zero padding, power of 2
int n_zeropadding;
//search region in fourier space
int excl_para;
//filter region in fourier space
int filter_para;
//crop after fft
int ncrop;
//motor position not used
//FLOATING motor_pos;
//threads for fftw
int fftw_thread;
//verbose option
int verbose;
Xgi_ip()
:roi(505,827,190,525)
,cam_orientation(0,0,1)
,beta_fix(0.0271)
,beta_rot(-0.0185)
,period_g1(3.75e-6)
,period_g2(2.000e-6)
,distance(144e-3)
,talbot_order(11)
,energy(9000)
,pixel_size(2.857e-6)
,n_zeropadding(16)
,excl_para(13)
,filter_para(3)
,ncrop(50)
//,motor_pos(18000)
,fftw_thread(1)
,verbose(0)
{};
};
struct Xgi_out
{
Matrix roi;
Matrix fft;
Matrix amp;
Matrix fft_re;
Matrix fft_im;
Matrix ifimg_re;
Matrix ifimg_im;
//wrapped phase
Matrix dphi_wrap;
//unwrapped phase
Matrix dphi_unwrap;
//phase with removed average horizontal slope
Matrix dphi;
Matrix polys;
Matrix alpha2D;
Matrix alpha;
Matrix laxis;
Matrix wpa;
Matrix phi_y;
Matrix yaxis;
Matrix hp;
FLOATING roc;
FLOATING droc;
Xgi_out(){};
~Xgi_out(){};
};
class Xgi
{
public:
Xgi();
~Xgi();
int compute(Matrix *put_in, Xgi_ip *ip, Xgi_out *put_out);
};
#endif //XGI_H
clear all;
dir_path = struct;
dir_path.expdata = '.\';
%parameter initializations (input comes from moire XGI calibration)
beta_fix=0.0271; %fixed absolute angle on phase grating
beta_rot=-0.0185; %rotated absolte angle on absorption grating
%reading moire interferograms and dark images (we have two)
img = read_img(dir_path,87,'sfb_');
dark = read_img(dir_path,61:62,'sfb_');
img=img-mean(dark,3);
roi = struct;
roi.x1= 506;roi.x2= 827;roi.y1= 191;roi.y2= 525; %this roi is found to well center the moire interf.
img=img(roi.y1:roi.y2,roi.x1:roi.x2);
ip = struct;
ip.beta_fix=beta_fix;
ip.beta_rot=beta_rot;
% period of G1 [m] phase grating
p1=3.75e-6;
ip.p1 = p1;
% period of G2 [m] absorption grating
p2 = 2.000e-6;
ip.p2 = p2;
% intergrating (G1-G2) distance [m]
ip.d = 144e-3;
%ip.d = 152e-3;
% % fractional Talbot order
ip.order = 11;
% X-ray energy [keV]
ip.kev = 9;
% X-ray wavelength [m]
ip.lambda = 12.3984191/ip.kev*1E-10; %*power(10,-10);
% magnification
ip.M0 =p2/p1*2;
% camera pixel size [m]
ip.pixel_size = 2.857e-6;
% some convenience parameters for the FFT filter
% zeropadding parameter
ip.n_zeropadding = 16; % n_zeropadding = 4; % must be power of 2
% parameter for search region in fourier space
ip.excl_para = 13;
% parameter for filter region in fourier space
ip.filter_para = 3;
% cropping after FFT
ip.ncrop = 50;
M=img;
[out] = xgi_func2(M,ip)
%show data output
ft=fftshift(out.fft);
amp = abs(out.amp);
ifimg = abs(out.ifimg);
wrap=out.dphi_wrap;
unwrap=out.dphi_unwrap;
polys = out.polys;
figure
subplot(1,6,1);
imagesc(M);colormap(flipud(bone)); title('moire image')
colormap default
ax_ft=subplot(1,6,2); plot(log(ft));set(ax_ft, 'YTick','');title('1D FFT')
subplot(1,6,3); imagesc(ifimg);title('1st order FFT filter')
subplot(1,6,4); imagesc(wrap);title('wrapped fringe phase')
ax2= subplot(1,6,5); plot(out.hp*1E9);title('height profile');ylabel('nanometers')
ax2 = subplot(1,6,6); plot(out.phi_y*1E9);title('spherical wavefront');ylabel('nanometers')
%%
function [ img_stack ] = read_img(dir_path,img_nb,file_base)
file_ext = 'h5';
cam_orientation = [ 1 0 1 ];
groupname = '/entry/instrument/detector/data' ;
img_stack = zeros(1024,1024,length(img_nb),'double');
for i_imgnb = 1:length(img_nb)
i_img = img_nb(i_imgnb);
filename = [ dir_path.expdata file_base sprintf('%04.0f.', i_img) file_ext ]
% put the array into the preallocated stack
if exist(filename,'file') == 2
img = double(h5read(filename,groupname));
if cam_orientation(1)
img = permute(img, [ 2 1 3 ]);
end
if cam_orientation(2)
img = flipdim(img, 2);
end
if cam_orientation(3)
img = flipdim(img, 1);
end
% img_stack(:,:,i_img-img_nb(1)+1) = img;
img_stack(:,:,i_imgnb) = img;
else
[ dir_path.expdata file_base sprintf('%04.0f.', i_img) file_ext ]
error('The file referred to in the above line does not exist')
end
end
end
function [out] = xgi_func2(img, ip)
ROC=0;
laxis=[];
wf=[]; hp=[];
aux=[];
out=struct;
[cols, rows] = size(img);
beta = -(ip.beta_rot - ip.beta_fix)/2;
dim_ana = 2;dim_ana_t=1;
% separate Hann window for each row
%hanning_window = repmat(hanning(size(img,2)),1,size(img,1))';
%s = myhanning(cols);
nn=1:1:cols;
HANN = .54 - 0.46*cos(2*pi.*nn/(cols+1));
HANN=HANN';
hanning_window = repmat(HANN,1,rows);%here removes the ' operator for Matlab2018
% zero padded, by a factor of 16
% Fourier transform along dimension 2 (rows) with hanning window and zero-padding
fimg = fft((img).*hanning_window,ip.n_zeropadding*size(img,2),dim_ana);
% projection of the FFT image on the orthogonal axis system
% mean of Fourier transform along the direction of interest (zero-padded by
% factor of 16)
mean_proj = mean(abs(fimg),dim_ana_t);
out.fft = mean_proj;
% get the position of the maximum value -> DC component position
[val_M, pos_ind_M] = max(mean_proj);
% get the position where the intensity is half the maximum value in order to have an estimate for the FWHM (Gaussian approximation of the DC component)
[~,pos_ind_fwhm] = min(abs(mean_proj-val_M/2));
% calculate from the FWHM the width of the Gaussian
G_width = 2*abs(pos_ind_M - pos_ind_fwhm) / 2.35;
%ensure max G_width
fprintf(2, 'pos_ind_M=%f pos_ind_fwhm=%f G_width=%f ', pos_ind_M, pos_ind_fwhm, G_width);
if G_width>10
G_width=10;
end
% find the correct Fourier peak
[val_m,pos_ind_p] = max(mean_proj(pos_ind_M+round(ip.excl_para*G_width):length(mean_proj)/2));
pos_ind_p = pos_ind_p+pos_ind_M+round(ip.excl_para*G_width)-1;
%shift
freq = pos_ind_p-pos_ind_M;
N2 = size(fimg,2);
period = N2*ip.pixel_size/freq * 10^6;
fprintf(2, 'period=%f \n',period)
% finge visibility
visib = val_m/val_M;
fprintf(2, 'freq=%f fringe visib=%.3f \n',freq, visib);
% mask which will be the filter
line_mask = zeros(size(fimg,1), size(fimg,2));
filter_width = round(ip.filter_para*G_width);
% 0-order and first order component processing
% filtering out the frequency of interest
line_mask(:, 1:filter_width+1) = 1;
line_mask(:, end-filter_width:end) = 1;
% mask out zero order
fimg_amp = fimg;
fimg_amp(~line_mask) = 0;
% shifts fourier transform (complex) so that frequency of interest is
% at zero
fimg = circshift(fimg, [0 -round(freq)]);
% sets everything outside mask to zero
fimg(~line_mask) = 0;
%%%%%%%%%%%%% changed here %%%%%%%%%%%%
fimg = circshift(fimg, [0 round(freq)]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 0-rder -> intensity on camera
amp = abs(ifft(fimg_amp,[],dim_ana));
% crop back to original image size
amp = amp(1:size(img,1), 1:size(img,2));
amp = amp ./ (hanning_window);
out.amp = amp;
%amp = sum(sum(amp(ncrop:end-ncrop+1,ncrop:end-ncrop+1),2));
%fprintf(2, 'intensity on image: %.0f \n',amp);
% inverse Fourier transform
ifimg = ifft(fimg,[],dim_ana);
% cropping in the inverse Fourier transformed image (back to original size)
ifimg = ifimg(1:size(img,1), 1:size(img,2));
out.ifimg = ifimg;
% differential fringe phase; returns the phase angle in radians
dphi_wrapped = angle(ifimg); %this is atan -> values netween -pi,pi -> need to unwrap
out.dphi_wrap = dphi_wrapped;
dphi_unwrap = unwrap(dphi_wrapped,pi,2);
dphi_unwrap = unwrap(dphi_unwrap,pi,1);
%dphi_unwrap = double(Miguel_2D_unwrapper(single(dphi_wrapped)));
out.dphi_unwrap = dphi_unwrap;
% coordinates in units of pixels
[ Y, X ] = ndgrid(size(img,1):-1:1, 1:size(img,2));
% fit a first order polynom through the unwrapped phase
% which is due to the fringe fundamental frequency,
%doing this is equvialent to mTakeda1984's shifting in fourier space
polys = zeros(2,size(img,1));
for ii = 1:size(img,1)
poly_tmp = polyfit(X(ii,:),dphi_unwrap(ii,:),1);
polys(1,ii) = poly_tmp(1);
polys(2,ii) = poly_tmp(2);
end
dphi = dphi_unwrap;
% alpha wavefront propagation angle; Eq. 2.74 Rutishauser thesis(calibration params)
alpha = (ip.M0*cos(ip.beta_fix)/cos(ip.beta_rot) - 1)*Y*ip.pixel_size/ip.d + ip.M0*ip.p2*cos(ip.beta_fix)/(2*pi*ip.d*(cos(ip.beta_rot))^2) * dphi;
% first order polynomial fit of the propagation angle of the wavefront through the image
polys = zeros(1,size(alpha,1)-2*ip.ncrop+2);
%%%%%%%% changed here %%%%%%%%%%%
for ii = ip.ncrop:size(alpha,1)-ip.ncrop+1
poly_tmp = polyfit(X(ii,:),alpha(ii,:),1);
polys(1,ii-ip.ncrop+1) = poly_tmp(1);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
out.polys = polys;
out.alpha2D = alpha;
%%
% run-time radius of curvature ROC determination
% connected through a geometric consideration to the propagation angle of the wavefront; Rutishauser thesis Eq. 2.26
ROC = ip.pixel_size/abs(mean(polys));
dROC = ip.pixel_size/(abs(mean(polys))^2)*std(polys);
%
fprintf(2, 'ROC %.3f +/- %.3f m \n', ROC, dROC);
out.ROC=ROC;
out.dROC=dROC;
%% find center of mass the image
X_hist=sum(img,2)'; X=1:size(img,1);
fprintf(2, 'X_hist %f X %f', size(X_hist),size(X))
i_x=round(sum(X.*X_hist)/sum(X_hist));
laxis = (1:size(img,1))*ip.pixel_size;
out.laxis=laxis;
%slope_data = alpha(:, i_x)';
%%%%%%%%%%%%%%% changed here %%%%%%%%%%%%%%%
slope_data = alpha(i_x,:);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
slope_data = slope_data-mean(slope_data);
out.wpa=slope_data;
phi_y1 = cumsum( slope_data'*ip.pixel_size );
%yaxis = (1:size(img,1))*ip.pixel_size;
%%%%%%%%%%%%%%% changed here %%%%%%%%%%%%%%
yaxis = (1:size(img,2))*ip.pixel_size;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% second order polynomial fit (spherical component)
py_phi = polyfit(yaxis', phi_y1, 2);
phi_y_aspherical = phi_y1 - polyval(py_phi, yaxis');
%% height profile
dphi_quant = ip.p2 / (2*pi * ip.d);
% scale for the image [m]
xaxis = (1:size(img,2))*ip.pixel_size;
% scale for the image [m]
yaxis = (1:size(img,1))*ip.pixel_size;
slope_data_y=slope_data';
% creating a rectangular grid with same dimensions as the data for the propagation angles
[ Yoffset_y, Xoffset_y ] = ndgrid(1:size(slope_data_y,1), 1:size(slope_data_y,2));
nbin=1;
[ Yoffset_y_all, Xoffset_y_all ] = ndgrid(1:size(slope_data_y,1), 1:size(slope_data_y,2));
phi_y_quant = ip.p2 / (ip.lambda * ip.d) * ip.pixel_size;
phi_y_surface_profile_quant = phi_y_quant * ip.lambda/(2*pi); %lambda is in meters
phi_y = cumsum( slope_data_y / dphi_quant, 1);
out.phi_y=phi_y * phi_y_surface_profile_quant;
%out.phi_y=cumsum( slope_data_y, 1);
% second order polynomial fit (spherical component)
py_phi = polyfit(xaxis', phi_y, 2);
% aspherical
%phi_y_aspherical = phi_y - repmat_size(polyval(py_phi, yaxis'), size(phi_y));
phi_y_aspherical = phi_y - polyval(py_phi, xaxis');
% metric dimension on the image
eff_pix_size = ip.pixel_size;
out.yaxis=yaxis;
% height profile hp
hp=phi_y_aspherical*phi_y_surface_profile_quant;
hp = hp - mean(hp);
out.hp = hp
INSTRUCTIONS
============
The XGI.tar contains C++ code for wavefront analysis (XGI) of an image saved in sfb_0087.h5 file. The code is subtacting an average of two dark images ( sfb_0061.h5 and sfb_0062.h5 ) before proceeding with the analysis. The code processing the XGI is generating images in order to follow and verify the processing pipeline [1]. The final output is a text file containig the wavefront profile in nanometers.
C-implementation
================
This folder contains C++ code for wavefront analysis (XGI) of an image saved in sfb_0087.h5 file. The code is subtacting an average of two dark images ( sfb_0061.h5 and sfb_0062.h5 ) before proceeding with the analysis. The code processing the XGI is generating images in order to follow and verify the processing pipeline, as described in Ref.[1]. The final output is a text file containig the wavefront profile in nanometers.
The XGI.7v is a Matlab equivalent of the XGI.tar. It should be noted that the C++ code was not generated by Matlab. Therefore both XGI.tar and XGI.7v are two independent software implementations producing the same results
Matlab
======
This folder contains Matlab scripts equivalent to the C++ implementation. It should be noted that the C++ code was not generated by Matlab.
References:
==========
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment