% AbelInversion_Resample_v1
%
%    Syntax: AbelInversion_Resample(FigureSwitch, PhaseShiftFilename, ParameterFilename)
%
%    Siurce data is resampled to be 10 times data points.
%
%    FigureSwitch = 1: figure on
%    FigureSwitch = 0: figure off
%
%    parameter file:
%       LeftCenter_RowIndex
%       RightCenter_RowIndex
%       ZeroRange_RowIndex
%       pixel (length calibration, unit: m)
%
%    output files:
%       (1) xxx_PhaseShift_Symmetrized.mat
%              Matlab binary file (.mat)
%              structure.PhaseShift_Symmetrized: 2-D matrix, symmetrized phase shift data (unit: rad)
%              structure.Position              : 1-D row vector, calibrated position (unit: micrometer)
%       (2) xxx_PhaseShift_Symmetrized.bmp
%              symmetrized phase shift image, 8-bit BMP file.
%       (3) xxx_dPdx.mat
%              Matlab binary file (.mat)
%              structure.dPdx    : 2-D matrix, dPhase/dx (unit: rad/m)
%              structure.Position: 1-D row vector, calibrated position (unit: micrometer)
%       (4) xxx_dPdx.bmp
%              dPhase/dx image, 8-bit BMP file.
%
%    returned value:
%       dPhase/dx, 2-D matrix (unit: rad/m)
%
%    Example:
%       a = AbelInversion_Resample(1,'IF-1-01_PhaseShift.mat', 'IF-1-01_Abel-parameters.txt');
%
%    Author: Hsu-hsin Chu (2003/5/30)
%

function result = AbelInversion_Resample(FigureSwitch, PhaseShiftFilename, ParameterFilename)

% importation of the parameter file
fid = fopen(ParameterFilename,'r');
while feof(fid) == 0
   line = fgetl(fid);
   eval(line);
end
fclose(fid);
Center_RowIndex = round((LeftCenter_RowIndex + RightCenter_RowIndex)/2);

% Loading of the phase-shift file
Source       = load(PhaseShiftFilename);
PhaseShift_1 = Source.PhaseShift;
[RowNumber_1, ColumnNumber_1] = size(PhaseShift_1);

% Rotation of the phase shift image
angle        = (RightCenter_RowIndex - LeftCenter_RowIndex)/ColumnNumber_1 * 180 / pi;
PhaseShift_2 = imrotate(PhaseShift_1, angle, 'bicubic');  % rotated phase shift data
[RowNumber_2, ColumnNumber_2] = size(PhaseShift_2);

% Average of the upper-half part and the lower-half part (symmetrize the phase shift image)
RowNumber_2_half  = min(Center_RowIndex, RowNumber_2 - Center_RowIndex);
PhaseShift_2_up   = PhaseShift_2(Center_RowIndex - RowNumber_2_half + 1:Center_RowIndex,:);
PhaseShift_2_down = PhaseShift_2(Center_RowIndex + 1:Center_RowIndex + RowNumber_2_half,:);
PhaseShift_3      = (flipud(PhaseShift_2_up) + PhaseShift_2_down) / 2;
%PhaseShift_3      = flipud(PhaseShift_2_up);
%PhaseShift_3      = PhaseShift_2_down;
[RowNumber_3, ColumnNumber_3] = size(PhaseShift_3);

% Noise reducing (cutting the phase image in the supposed-to-be-zero region)
PhaseShift_4 = PhaseShift_3(1:ZeroRange_RowIndex,:);
[RowNumber_4, ColumnNumber_4] = size(PhaseShift_4);

% Data resample
n_resample  = 10;
y_original  = [0:RowNumber_4-1]';
y_resample  = linspace(0, RowNumber_4-1, RowNumber_4 * n_resample)';
PhaseShift_5 = zeros(RowNumber_4 * n_resample, ColumnNumber_4);
for i = 1:ColumnNumber_4
    PhaseShift_5(:,i) = spline(y_original,PhaseShift_4(:,i),y_resample);
end
pixel = pixel/n_resample;

% Solving dPhase/dx (dPdx) by Abel inversion
dPdx_half = Abel(PhaseShift_5) / pixel; % dPdx unit: rad/m, pixel: length calibration

% Data exportation
PhaseShift_Symmetrized = [flipud(PhaseShift_5);PhaseShift_5];   % (unit: rad)
dPdx                   = [flipud(dPdx_half);dPdx_half];         % (unit: rad/m)
Position               = [1:ColumnNumber_4] * pixel * 10^6;     % (unit: micrometer)

result = dPdx;

temp = strtok(PhaseShiftFilename,'_');
OutputFilename_1 = strcat(temp,'_PhaseShift_Resampled_Symmetrized.mat');
OutputFilename_2 = strcat(temp,'_PhaseShift_Resampled_Symmetrized.bmp');
OutputFilename_3 = strcat(temp,'_Resampled_dPdx.mat');
OutputFilename_4 = strcat(temp,'_Resampled_dPdx.bmp');

save(OutputFilename_1,'PhaseShift_Symmetrized','Position');

max_PS = max(max(PhaseShift_Symmetrized));
min_PS = min(min(PhaseShift_Symmetrized));
PhaseShift_image = (PhaseShift_Symmetrized - min_PS)./(max_PS - min_PS)*255;
imwrite(PhaseShift_image,jet(256),OutputFilename_2,'bmp');

save(OutputFilename_3,'dPdx','Position');

max_dPdx = max(max(dPdx));
min_dPdx = min(min(dPdx));
dPdx_image = (dPdx - min_dPdx)./(max_dPdx - min_dPdx)*255;
imwrite(dPdx_image,jet(256),OutputFilename_4,'bmp');



if FigureSwitch
    x = [1:RowNumber_2]';
    y1 = PhaseShift_2(:,round(ColumnNumber_2/4));
    y2 = PhaseShift_2(:,round(ColumnNumber_2/4*2));
    y3 = PhaseShift_2(:,round(ColumnNumber_2/4*3));
    y_max = max([y1;y2;y3])+(max([y1;y2;y3])-min([y1;y2;y3]))/5;
    y_min = min([y1;y2;y3])-(max([y1;y2;y3])-min([y1;y2;y3]))/5;
    x_c = [Center_RowIndex, Center_RowIndex];
    y_c = [y_min,y_max];
    
    subplot(3,2,1), imagesc(PhaseShift_1), title('original phase shift');
    PhaseShift_2_image = PhaseShift_2;
    PhaseShift_2_image(Center_RowIndex:Center_RowIndex + 1,:) = min(min(PhaseShift_2));
    subplot(3,2,3), imagesc(PhaseShift_2_image),  title('rotated phase shift');
    subplot(3,2,5), plot(x,y1,x,y2,x,y3,x_c,y_c), title('vertical line profile'),
    axis([-inf, inf, y_min, y_max]), xlabel('(pixel)'), ylabel('(rad)');
    
    PhaseShift_3_image = PhaseShift_3;
    PhaseShift_3_image(ZeroRange_RowIndex,:) = min(min(PhaseShift_3));
    subplot(3,2,2), imagesc(PhaseShift_3_image),  title('symmetrized phase shift');
    subplot(3,2,4), plot(PhaseShift_3(ZeroRange_RowIndex,:)), xlabel('(pixel)'), ylabel('(rad)'),
    title(['horizontal profile, row ' num2str(ZeroRange_RowIndex)]),;
    subplot(3,2,6), imagesc(PhaseShift_4),        title('noise cutted phase shift');
    
    figure;
    subplot(3,1,1), imagesc(PhaseShift_Symmetrized), title('symmetrized and resampled phase shift');
    subplot(3,1,2), imagesc(dPdx),                   title('dPhase/dx');
    subplot(3,1,3), contour(flipud(dPdx),15),        title('dPhase/dx');
end



%-------------------------------------------------------------------------------------------------------
%  Abel inversion
%-------------------------------------------------------------------------------------------------------

function result = Abel(PhaseShift)

[RowNumber, ColumnNumber] = size(PhaseShift);
r  = [1:RowNumber]';
rr = cumsum(ones(RowNumber));
yy = rr';
%     1         1          1 1 1...1          1 2 3...N
%     2         2          2 2 2...2          1 2 3...N
% r = 3,    y = 3,    rr = 3 3 3...3,    yy = 1 2 3...N
%     :         :              :                  :
%     N         N          N N N...N          1 2 3...N

% Abel inversion kernal:
%    A(i,j) = 1 / sqrt( y(j)^2 - r(i)^2 ) for j > i, A(i,j) = 0 for j <= i

A  = yy.^2 - rr.^2;   % set A(i,j) = y(j)^2 - r(i)^2
A  = A - (A == 0);    % set zero element to become -1, make 1./A possible
A  = 1 ./ A;          % set A(i,j) = 1 / ( y(j)^2 - r(i)^2 )
A  = A .* (A > 0);    % set A(i,j) = 0 for j <= i, keep other A(i,j)
A  = sqrt(A);         % set A(i,j) = 1 / sqrt( y(j)^2 - r(i)^2 )

% Abel inversion delta phase shift:
%    DeltaPS(j,k) = PhaseShift(j+1,k)-PhaseShift(j,k)

DeltaPS = zeros(RowNumber, ColumnNumber);
DeltaPS(1:RowNumber-1,:) = PhaseShift(2:RowNumber,:) - PhaseShift(1:RowNumber-1,:);

result = -1/pi*A*DeltaPS;









