interferogramphaseshift.m

干涉仪相位分析程序

% InterferogramPhaseShift_v7
% 
%    Syntax: InterferogramPhaseShift(FigureSwitch, SourceFilename, ReferenceFilename, ParameterFilename)
%
%       FigureSwitch:
%          FigureSwitch = 0: no figures
%          FigureSwitch = 1: generate figures
%
%       source & reference file:
%          BMP, TIFF or FITS... image file, vertical fringes.
%
%       parameter file (ASCII-type file):
%          element 1: Row_min
%          element 2: Row_max
%          element 3: Column_min
%          element 4: Column_max
%          element 5: Mask_range_up
%          element 6: Mask_range_down
%
%       returned value:
%          phase shift (2-D matrix)
%
%       output files:
%          (1) Matlab binary file, 2-D matrix (.mat)
%          (2) 8-bit BMP file
%
%       example:
%          a = InterferogramPhaseShift(1, 'IF-1-04.fit', 'IF-1-01.fit', 'IF-parameters.txt');
%
%    Author: Hsu-hsin Chu (2003/5/1)

function result = InterferogramPhaseShift_v7(FigureSwitch, SourceFilename, ReferenceFilename, ParameterFilename)

% Loading the parameter file
ProgramParameters = load(ParameterFilename);
Rmin = ProgramParameters(1);
Rmax = ProgramParameters(2);
Cmin = ProgramParameters(3);
Cmax = ProgramParameters(4);
filter_range_up   = ProgramParameters(5);
filter_range_down = ProgramParameters(6);

% importation of the interferogram
SourceData    = OpenImageFile(SourceFilename);
SourceData    = double(SourceData(Rmin:Rmax,Cmin:Cmax));
ReferenceData = OpenImageFile(ReferenceFilename);
ReferenceData = double(ReferenceData(Rmin:Rmax,Cmin:Cmax));
[RowNumber, ColumnNumber] = size(SourceData);

% rescaling, Fourier transformation, and filtering
filter_range      = [filter_range_up, filter_range_down];
SourceSpectrum    = Rescaling_FT_Filtering(SourceData, filter_range);
ReferenceSpectrum = Rescaling_FT_Filtering(ReferenceData, filter_range);

Source_FT      = SourceSpectrum.original;
Source_FT_F    = SourceSpectrum.filted;
Reference_FT   = ReferenceSpectrum.original;
Reference_FT_F = ReferenceSpectrum.filted;


% inverse Fourier transform
Source_back         = ifft(Source_FT_F')';
Source_back_real    = real(Source_back);
Reference_back      = ifft(Reference_FT_F')';
Reference_back_real = real(Reference_back);

if FigureSwitch
    figure;
    colormap(gray(256));        
    RowIndex = double(int8(RowNumber/2));
    temp = strcat(' spectrum (row ', num2str(RowIndex), ')');
    subplot(4,2,1), imagesc(SourceData),                   title('source');
    subplot(4,2,2), imagesc(ReferenceData),                title('reference');
    subplot(4,2,3), plot(abs(Source_FT(RowIndex,:))),      title(strcat('source',temp));
    subplot(4,2,4), plot(abs(Reference_FT(RowIndex,:))),   title(strcat('reference',temp));
    subplot(4,2,5), plot(abs(Source_FT_F(RowIndex,:))),    title(strcat('filted',temp));
    subplot(4,2,6), plot(abs(Reference_FT_F(RowIndex,:))), title(strcat('filted',temp));
    subplot(4,2,7), imagesc(Source_back_real),             title('filtered source');
    subplot(4,2,8), imagesc(Reference_back_real),          title('filtered reference');
end

% phase shift calculation
PhaseShift = -angle(Source_back ./ Reference_back);
PhaseShift = PhaseShift(:,1:ColumnNumber);
PhaseShift = PhaseCompensation(PhaseShift);
PhaseShift = BackgroundSubstration(PhaseShift);

PhaseShift(:,ColumnNumber-14:ColumnNumber) = 0;

if FigureSwitch
    figure;
    colormap(gray(256));
    colormap(jet);
    imagesc(PhaseShift), title(strtok(SourceFilename,'.'));
end

% data exportation
temp = strtok(SourceFilename,'.');
OutputFilename  = strcat(temp,'_PhaseShift.mat');
save(OutputFilename,'PhaseShift');
OutputImageName = strcat(temp,'_PhaseShift.bmp');
MaxValue = max(max(PhaseShift));
MinValue = min(min(PhaseShift));
OutputImage = (PhaseShift - MinValue)/(MaxValue - MinValue)*255;
colormap(gray(256));
colormap(jet);
imwrite(OutputImage,colormap,OutputImageName,'bmp');


result = PhaseShift;



%--------------------------------------------------------------------------------------------------
% OpenImageFile
%    Open an image file. ('FITS','TIF','BMP','JPEG','PNG',and 'HDF'... formats)
%    Return a 2D matrix.
%
%    author: Hsu-hsin Chu (2003/1/1)
%--------------------------------------------------------------------------------------------------

function result = OpenImageFile(Filename)

[SourceFilename_1 SourceFilename_2] = strtok(Filename,'.');
SourceFilename_2 = SourceFilename_2(2:length(SourceFilename_2));

if strcmpi(SourceFilename_2, 'FITS')
   source = fitsread(Filename);
elseif strcmpi(SourceFilename_2, 'FIT')
   source = fitsread(Filename);
elseif strcmpi(SourceFilename_2, 'FTS')
   source = fitsread(Filename);
else
   source = imread(Filename, SourceFilename_2);
end

dimension = size(source,3);
if dimension == 3
   source = RGB2gray(source);
   image = double(source);
else
   image = double(source);
end
result = image;


%--------------------------------------------------------------------------------------------------
% Rescaling_FT_Filtering
%
%    Rescaling, Fourier transformation and Filtering of the input interferogram
%--------------------------------------------------------------------------------------------------

function result = Rescaling_FT_Filtering(input, filter_range)

[RowNumber, ColumnNumber] = size(input);

% rescaling
temp = double(input);
MaxValue = max(max(temp));
MinValue = min(min(temp));
temp = (temp - MinValue)/(MaxValue - MinValue);

% Fourier transformation and filtering
temp_FT = fft(temp')';

center = (filter_range(2) + filter_range(1))/2;   % super Gaussian center
radius = (filter_range(2) - filter_range(1))/2;   % super Gaussian radius
x = ones(RowNumber, ColumnNumber);
x = cumsum(x')';
mask = exp(-(x-center).^4/radius^4);            % (super Gaussian mask)
%mask = zeros(RowNumber, ColumnNumber);
%mask(:,filter_range(1):filter_range(2)) = 1;

temp_FT_F = temp_FT .* mask;

spectrum.original = temp_FT;
spectrum.filted   = temp_FT_F;

result = spectrum;


%--------------------------------------------------------------------------------------------------
% PhaseCompensation
%
%    Compensate the phase jump from Pi to -Pi and from -Pi to Pi.
%--------------------------------------------------------------------------------------------------

function result = PhaseCompensation(PhaseShift)

[RowNumber, ColumnNumber] = size(PhaseShift);

diff_col = zeros(RowNumber, ColumnNumber);
diff_col(:,2:ColumnNumber) = PhaseShift(:,1:ColumnNumber-1)-PhaseShift(:,2:ColumnNumber);
flag_1 = (diff_col > 3) - (diff_col < -3);
flag_2 = cumsum(flag_1,2);

PhaseShift_2 = PhaseShift + 2 * pi * flag_2;

diff_row = zeros(RowNumber, 1);
diff_row(2:RowNumber,1) = PhaseShift_2(1:RowNumber-1,1) - PhaseShift_2(2:RowNumber,1);
%diff_row(2:RowNumber,ColumnNumber) = PhaseShift_2(1:RowNumber-1,ColumnNumber) - PhaseShift_2(2:RowNumber,ColumnNumber);
flag_3 = (diff_row > 3) - (diff_row < -3);
flag_4 = cumsum(flag_3);
flag_4(:,2:ColumnNumber) = 0;
flag_4 = cumsum(flag_4,2);

result = PhaseShift_2 + 2 * pi * flag_4;


%--------------------------------------------------------------------------------------------------
% BackgroundSubstration
%
%    Substract of the linear background phase
%--------------------------------------------------------------------------------------------------

function result = BackgroundSubstration(PhaseShift)

[RowNumber, ColumnNumber] = size(PhaseShift);

x = [1:ColumnNumber];
Background = mean(PhaseShift([11:20,RowNumber-19:RowNumber-10],:));
LinearFit_coefficient     = polyfit(x, Background,1);
LinearFit_BackgroundPhase = polyval(LinearFit_coefficient, x);
BackgroundPhase = zeros(RowNumber, ColumnNumber);
for i = 1:ColumnNumber
    BackgroundPhase(:,i) = LinearFit_BackgroundPhase(i);
end

result = PhaseShift - BackgroundPhase;