% main file to run SLIC using the configuration specified in SLICconfig.mat
function BaSiC_main_AW(configurefile)
clearvars -except configurefile;
close all; clc;
msgbox('Starting BaSiC, just sit back and Relax. Wait for output in Console');
load(configurefile);
for wavelength_num = 1:length(FluorWL_group)
    FluorWL = char(FluorWL_group(wavelength_num));
    temp = textscan(FluorWL,'%s %s','delimiter', '.');
    wavelength = temp{1};
   
    %% Step I: Generation of flatfield_ap (a base flatfield for all positions)
    % output: flatfield_ap in outputDirectory; if exist, please set runStepI to
    % false
    if ~runStepI
    else
        if runAllPositionsI
            files = dir(inputDirectory);
            dirFlags = [files.isdir];
            subFolders = files(dirFlags);  % extract only those that are directories
            %subFolders = subFolders(arrayfun(@(x)~strcmp(x.name(1),'.'),subFolders));
            subFolders = subFolders(arrayfun(@(x)~isempty(strfind(x.name,sprintf('%s_p',movieID))),subFolders));
    %         IndexC = strfind(subFolders.name, sprintf('%s_p',movieID)); % position subfolder 
    %         Index = find(not(cellfun('isempty', IndexC)));
            positionVectorI = zeros(1,length(subFolders));
            for i = 1:length(subFolders)
                positionVectorI(i) = sscanf(subFolders(i).name, [sprintf('%s_p',movieID) '%d']);
            end
        end
        % read through all positionVectorI
        totalImageNum = 0; 
        incomplete=1;%some experiments do not acquire wavelengts for all positions AW
        for i = 1:length(positionVectorI)
            positionDirectory = sprintf('%s%s_p%.4d/',inputDirectory,movieID,positionVectorI(i));
            imagefiles = dir([positionDirectory '*' FluorWL]);
            totalImageNum = totalImageNum+length(imagefiles);
            empty=isempty(imagefiles); %some experiments do not acquire wavelengts for a positions AW
            if incomplete && ~empty %AW
                image = imread([positionDirectory '/' imagefiles(1).name]);
                sizeX = size(image,1);
                sizeY = size(image,2);
                incomplete=0; %AW
            end
        end
        
        sampleNumPerPosition = ceil(min(totalImageNum,500)./length(positionVectorI)); 
        if ~incomplete %AW
            fi = [];
            for i = 1:length(positionVectorI)
                positionDirectory = sprintf('%s/%s_p%.4d/',inputDirectory,movieID,positionVectorI(i));
                imagefiles = dir([positionDirectory '*' FluorWL]);
                sample = randperm(length(imagefiles),min(length(imagefiles),sampleNumPerPosition));
                fi_temp = zeros(sizeX, sizeY, length(sample));
                for j = 1:length(sample)
                    fi_temp(:,:,j) = imread([positionDirectory '/' imagefiles(sample(j)).name]);
                end
                fi = cat(3,fi,fi_temp);
                clear fi_temp;
            end
            % run SLIC and save estimated shading profile
            [flatfield_ap,~]=BaSiC(double(fi),'darkfield','true','lambda',lambdaFlatfield_allposition,'lambda_darkfield',lambdaFlatfield_allposition*0.5);
             save([outputDirectory 'flatfield_allposition_' char(wavelength) '.mat'],'flatfield_ap')
             figure(2); imagesc(flatfield_ap);colorbar;title(['flatfield-allposition-' char(wavelength)]);
             saveas(gcf,[outputDirectory 'flatfield_allposition_' char(wavelength)],'tif');
        end  %AW

    end


    %% step II: Generation of position-wise flatfield and darkfield 
    % output: flatfield and darkfield in position subdirectory under
    % outputDirectory as well as a overview csv file
    % runStepII = 'true'; % bool, if true then run SLIC to generate new position-wise flatfield & darkfield and please specify the following parameters
    % runAllPositionsII = 'true'; %bool, if true run across all positions in the movie folder. If false, specity positionVector
    % positionVectorII = [1 2 3 10 15 20]; % or in the format of 1:36

    if ~runStepII
    else
        if useFlatfield_allposition
            flatfield_ap_file = [outputDirectory 'flatfield_allposition_' char(wavelength) '.mat'];
            if exist(flatfield_ap_file,'file');
                flatfield_ap = importdata(flatfield_ap_file);
            else
                msgbox(sprintf('%s%s%s','There is no flatfield_allposition file, please run Step I and include positions where wavelength ', char(wavelength), ' was included'));
                msgbox('The BaSiC Run will terminate now. Either include positions with the given wavenlength or leave this wavelength out');
                error('%s%s%s','There is no flatfield_allposition file, please run Step I and include positions where wavelength ', char(wavelength), ' was included');
            end
        else
            flatfield_ap = [];
        end
        if runAllPositionsII
            files = dir(inputDirectory);
            dirFlags = [files.isdir];
            subFolders = files(dirFlags);  % extract only those that are directories
            %subFolders = subFolders(arrayfun(@(x)~strcmp(x.name(1),'.'),subFolders));
            subFolders = subFolders(arrayfun(@(x)~isempty(strfind(x.name,sprintf('%s_p',movieID))),subFolders));
    %         IndexC = strfind(subFolders.name, sprintf('%s_p',movieID)); % position subfolder 
    %         Index = find(not(cellfun('isempty', IndexC)));
            positionVectorII = zeros(1,length(subFolders));
            for i = 1:length(subFolders)
                positionVectorII(i) = sscanf(subFolders(i).name, [sprintf('%s_p',movieID) '%d']);
            end
        end
        % prepare log file
        %LastName = {'Smith';'Johnson';'Williams';'Jones';'Brown'};
        goodestimation_measure1 = false(length(positionVectorII),1);
        goodestimation_measure2 = false(length(positionVectorII),1);
        % read through all positionVectorI
        incomplete=1;%some experiments do not acquire wavelengts for all positions AW
        for i = 1:length(positionVectorII)
            positionDirectory = sprintf('%s/%s_p%.4d/',inputDirectory,movieID,positionVectorII(i));
            outputPositionDirectory = sprintf('%s/%s_p%.4d/',outputDirectory,movieID,positionVectorII(i));
            mkdir(outputPositionDirectory);
            imagefiles = dir([positionDirectory '*' FluorWL]);
            if isempty(temporalRangeII)
                sample = randperm(length(imagefiles),min(length(imagefiles),300));
            else
                timepoint = zeros(length(imagefiles),1);
                for k = 1:length(imagefiles)
                    time_str = regexp(imagefiles(k).name,'t(?<t>\d+)','match');
                    timepoint(k) = sscanf(time_str{1},'t%d');
                end
                valid_idx = find(timepoint>=temporalRangeII(1)&timepoint<=temporalRangeII(2));
                sample = valid_idx(randperm(length(valid_idx),min(length(valid_idx),300)));
            end
            
            empty=isempty(imagefiles); %some experiments do not acquire wavelengts for a positions AW
            if incomplete && ~empty %AW
                image = imread([positionDirectory imagefiles(1).name]);
                sizeX = size(image,1);
                sizeY = size(image,2);
                incomplete=0; %AW
            end
            fi = zeros(sizeX, sizeY, length(sample));
            for j = 1:length(sample)
                fi(:,:,j) = imread([positionDirectory imagefiles(sample(j)).name]);
            end
            if ~isempty(fi) %AW
                if and(isempty(lambdaFlatfield),isempty(lambdaDarkfield)) 
                    [flatfield,darkfield]=BaSiC(double(fi),'darkfield','true','initial_flat',flatfield_ap);
                     %adjust parameter if the obtained flatfield is not smooth enough
                    flatfield_dct = dct2(flatfield);
                    %darkfield_cmp = importdata('darkfield_test.mat');
                    darkfield_norm = (darkfield-mean(darkfield(:)))./mean(darkfield(:));
                    darkfield_dct = dct2(darkfield_norm);
                    nonsmoothindex = max(sum(abs(flatfield_dct(:))./max(abs(flatfield_dct(:)))>1e-4)./numel(flatfield_dct),sum(abs(darkfield_dct(:))./max(abs(darkfield_dct(:)))>1e-4)./numel(darkfield_dct));
                    if nonsmoothindex > 0.002
                        lambda_dark = (nonsmoothindex-0.002)/0.01+0.05;
                        lambda_dark = min(lambda_dark,2.0);
                        lambda_flat = lambda_dark./0.4;
                        [flatfield,darkfield]=BaSiC(double(fi),'est_dark','true','initial_flat',flatfield_ap,'lambda',lambda_flat,'lambda_darkfield',lambda_dark);
                    end
                else
                    lambda_dark = lambdaDarkfield(i);
                    lambda_flat = lambdaFlatfield(i);
                    [flatfield,darkfield]=BaSiC(double(fi),'est_dark','true','initial_flat',flatfield_ap,'lambda',lambda_flat,'Lambda_darkfield',lambda_dark);
                end
                if isempty(temporalRangeII)
                    flatfield_file = [outputPositionDirectory 'flatfield_' char(wavelength)];
                    darkfield_file = [outputPositionDirectory 'darkfield_' char(wavelength)];
    %                 save([outputPositionDirectory 'flatfield_' char(wavelength) '.mat'],'flatfield');
    %                 save([outputPositionDirectory 'darkfield_' char(wavelength) '.mat'],'darkfield'); 
                else
                    flatfield_file = [outputPositionDirectory 'flatfield_' char(wavelength) '_t' num2str(temporalRangeII(1)) 'to' num2str(temporalRangeII(2)) ];
                    darkfield_file = [outputPositionDirectory 'darkfield_' char(wavelength) '_t' num2str(temporalRangeII(1)) 'to' num2str(temporalRangeII(2)) ];
    %                 save([outputPositionDirectory 'flatfield_' char(wavelength) '_t' num2str(temporalRangeII(1)) 'to' num2str(temporalRangeII(2)) '.mat'],'flatfield');
    %                 save([outputPositionDirectory 'darkfield_' char(wavelength) '_t' num2str(temporalRangeII(1)) 'to' num2str(temporalRangeII(2)) '.mat'],'darkfield'); 
                end
                saveToTif(flatfield, flatfield_file);
                saveToTif(darkfield, darkfield_file);


                if useFlatfield_allposition
                    error1 = sum(abs(flatfield(:)-flatfield_ap(:)))./sum(abs(flatfield(:))); % error between flatfield and the common flatfield_ap
                    if error1(1)<0.005
                        goodestimation_measure1(i) = true;
                    end
                end
                 % estimate flatfield and darkfield with segmentation considered
                 if existSegmentation
                    fi_seg = false(sizeX, sizeY, length(sample));
                    segDirectory = sprintf('%s/%s_p%.4d/',segmentationDirectory,movieID,positionVectorII(i));
                    segfiles = dir([segDirectory '*mask.png']);
                    segfilesname = {segfiles.name}';
                    matchsegfilesname = cell(1,length(sample));                  
                    temp_name = segfiles(sample(1)).name;
                    expression = 't(?<t>\d+)';
                    for  j = 1:length(sample)
                        temp_str = regexp(imagefiles(sample(j)).name,'t(?<t>\d+)','match'); 
                        matchsegfilesname{j} = regexprep(temp_name,expression,temp_str);
                    end
                    matchindex = find(ismember(matchsegfilesname,segfilesname));            
                    if isempty(segfiles)
                        error('Cannot find segmentation files, please check the path to segmentation files is correct');
                    else
                        for j = 1:length(sample(matchindex))
                            fi_seg(:,:,matchindex(j)) = im2bw(imread([segDirectory '/' segfiles(sample(matchindex(j))).name]),1e-6);
                        end
                        if and(exist('lambda_dark','var'),exist('lambda_flat','var'))
                            [flatfield_seg, darkfield_seg] = BaSiC(fi,'est_dark','true','initial_flat',flatfield_ap,'segmentation',fi_seg,'lambda',lambda_flat,'Lambda_darkfield',lambda_dark);
                        else
                            [flatfield_seg, darkfield_seg] = BaSiC(fi,'est_dark','true','initial_flat',flatfield_ap,'segmentation',fi_seg);
                        end
                        if useSegmentation  % overwrite flatfield and darkfield without using segmenation
                            if isempty(temporalRangeII)
                                flatfield_file = [outputPositionDirectory 'flatfield_' char(wavelength)];                            
                                darkfield_file = [outputPositionDirectory 'darkfield_' char(wavelength)];
    %                             save([outputPositionDirectory 'flatfield_' char(wavelength) '.mat'],'flatfield_seg');
    %                             save([outputPositionDirectory 'darkfield_' char(wavelength) '.mat'],'darkfield_seg'); 
                            else
                                flatfield_file = [outputPositionDirectory 'flatfield_' char(wavelength) '_t' num2str(temporalRangeII(1)) 'to' num2str(temporalRangeII(2))];
                                darkfield_file = [outputPositionDirectory 'darkfield_' char(wavelength) '_t' num2str(temporalRangeII(1)) 'to' num2str(temporalRangeII(2))]; 
                            end  
                            saveToTif(flatfield_seg, flatfield_file);
                            saveToTif(darkfield_seg, darkfield_file);
                        end
                    end
                    error2 = sum(abs(flatfield(:)-flatfield_seg(:)))./sum(abs(flatfield(:))); % error between flatfield and the flatfield_seg      
                    if error2(1)<0.005
                        goodestimation_measure2(i) = true;
                    end        
                 end
          end %AW
        end
        % check the reliabiblity of estimated flatfield and darkfield and
            % write to a table file
         
        qualitymeasurefile =  [outputDirectory 'qualitymeasure_' char(wavelength) '.txt'];   
        if exist(qualitymeasurefile,'file')
            temp_txt = regexp( fileread(qualitymeasurefile), '\n', 'split');
            temp_txt = temp_txt(~cellfun('isempty',temp_txt)); 
            for i = 1:length(positionVectorII)
                temp_IndexC = strfind(temp_txt, sprintf('%s_p%.4d',movieID,positionVectorII(i)));
                temp_Index = find(not(cellfun('isempty', temp_IndexC)));
                if not(isempty(temp_Index))
                    temp_txt{temp_Index} = sprintf('%s_p%.4d,%d,%d',movieID,positionVectorII(i),goodestimation_measure1(i),goodestimation_measure2(i));
                else
                    temp_txt{length(temp_txt)+1} = sprintf('%s_p%.4d,%d,%d',movieID,positionVectorII(i),goodestimation_measure1(i),goodestimation_measure2(i));
                end
            end
            fid = fopen(qualitymeasurefile,'w');
            fprintf(fid,'%s\n',temp_txt{:});
            fclose(fid);
        else
            fileID = fopen(qualitymeasurefile,'w');
            fprintf(fileID,'%s,%s,%s\n', 'Position','Meansure1','Measure2');
            for i = 1:length(positionVectorII)
                fprintf(fileID,'%s_p%.4d,%d,%d\n',movieID,positionVectorII(i),goodestimation_measure1(i),goodestimation_measure2(i));
            end
            fclose(fileID);
        end
        
    end

    %% step III: Generation of baseline fluorescence to correct background bleaching
    % output: basefluor in position subdirectory under
    % outputDirectory
    % runStepIII = 'true'; %bool, if true then run BaSiC to generate new basefluor and please specify the following parameters
    % runAllPositionsIII = 'true'; %bool, if true run across all positions in the movie folder. If false, specity positionVector
    % positionVectorIII = [1 2 3 10 15 20]; % or in the format of 1:36
    if ~runStepIII
    else
        if runAllPositionsIII
            files = dir(inputDirectory);
            dirFlags = [files.isdir];
            subFolders = files(dirFlags);  % extract only those that are directories
            %subFolders = subFolders(arrayfun(@(x)~strcmp(x.name(1),'.'),subFolders));
            subFolders = subFolders(arrayfun(@(x)~isempty(strfind(x.name,sprintf('%s_p',movieID))),subFolders));
    %         IndexC = strfind(subFolders.name, sprintf('%s_p',movieID)); % position subfolder 
    %         Index = find(not(cellfun('isempty', IndexC)));
            positionVectorIII = zeros(1,length(subFolders));
            for i = 1:length(subFolders)
                positionVectorIII(i) = sscanf(subFolders(i).name, [sprintf('%s_p',movieID) '%d']);
            end
        end
        % read through all positionVectorIII
        for i = 1:length(positionVectorIII)
            positionDirectory = sprintf('%s/%s_p%.4d/',inputDirectory,movieID,positionVectorIII(i));
            outputPositionDirectory = sprintf('%s/%s_p%.4d/',outputDirectory,movieID,positionVectorIII(i));
            imagefiles = dir([positionDirectory '*' FluorWL]);
            blocksamplenum = 1000;    % divide the images into several blocks %AW180418 TYHIS IS shit, make it according to timepoints
            if isempty(temporalRangeIII) 
                blocknum = ceil(length(imagefiles)./blocksamplenum);
                valid_idx = 1:length(imagefiles);
            else          
                timepoint = zeros(length(imagefiles),1);
                for k = 1:length(imagefiles)
                    temp_str = regexp(imagefiles(k).name,'t(?<t>\d+)','match');
                    timepoint(k) = sscanf(temp_str{1},'t%d');
                end
                %valid_idx = find(timepoint>=temporalRangeIII(1)&timepoint<=max(temporalRangeIII(2));
                valid_idx = find(timepoint>=min(temporalRangeIII)&timepoint<=max(temporalRangeIII)); %AW180417
                blocknum = ceil(length(valid_idx)./blocksamplenum);
            end
            % read in the flatfield and darkfield
            if isempty(temporalRangeIII)
                flatfield_file = [outputPositionDirectory 'flatfield_' char(wavelength) '.tif'];
                darkfield_file = [outputPositionDirectory 'darkfield_' char(wavelength) '.tif'];
%                 flatfield_file = [outputPositionDirectory 'flatfield_' char(wavelength) '.mat'];
%                 darkfield_file = [outputPositionDirectory 'darkfield_' char(wavelength) '.mat'];
            else
                flatfield_file = [outputPositionDirectory 'flatfield_' char(wavelength) '_t' num2str(temporalRangeIII(1)) 'to' num2str(temporalRangeIII(2)) '.tif'];
                darkfield_file = [outputPositionDirectory 'darkfield_' char(wavelength) '_t' num2str(temporalRangeIII(1)) 'to' num2str(temporalRangeIII(2)) '.tif'];
%                 flatfield_file = [outputPositionDirectory 'flatfield_' char(wavelength) '_t' num2str(temporalRangeIII(1)) 'to' num2str(temporalRangeIII(2)) '.mat'];
%                 darkfield_file = [outputPositionDirectory 'darkfield_' char(wavelength) '_t' num2str(temporalRangeIII(1)) 'to' num2str(temporalRangeIII(2)) '.mat'];
            end 
            if ~and(exist(flatfield_file,'file'),exist(darkfield_file,'file')) && isempty(imagefiles) %AW
               fprintf('No image acquired at this wavelength at this position');
            elseif ~and(exist(flatfield_file,'file'),exist(darkfield_file,'file')) && ~isempty(imagefiles) %Aw
            error('There is no suitable flatfield or darkfield file, please run Step II');
            else
                flatfield = importdata(flatfield_file);
                darkfield = importdata(darkfield_file);
                sizeX = size(flatfield,1);
                sizeY = size(flatfield,2); 
                fi_base = zeros(1,length(valid_idx));
                flat_amp = zeros(1,length(valid_idx));
                time_str = cell(1,length(valid_idx));               
                for block = 1:blocknum
                    if isempty(temporalRangeIII)
                        sample = (block-1)*blocksamplenum+1:min(block*blocksamplenum,length(imagefiles));
                    else
                        sample = valid_idx((block-1)*blocksamplenum+1:min(block*blocksamplenum,length(valid_idx)));
                    end
                    fi = zeros(sizeX, sizeY, length(sample));
                    for j = 1:length(sample)
                        fi(:,:,j) = imread([positionDirectory imagefiles(sample(j)).name]);
                        time_str{j} = regexp(imagefiles(sample(j)).name,'t(?<t>\d+)','match');
                    end
                     if existSegmentation
                        fi_seg = false(sizeX, sizeY, length(sample));
                        segDirectory = sprintf('%s/%s_p%.4d/',segmentationDirectory,movieID,positionVectorIII(i));
                        segfiles = dir([segDirectory '*mask.png']);
                        segfilesname = {segfiles.name}';
                        matchsegfilesname = cell(1,length(sample));                  
                        temp_name = segfiles(sample(1)).name;
                        expression = 't(?<t>\d+)'; 
                        for  j = 1:length(sample)
                            matchsegfilesname{j} = regexprep(temp_name,expression,time_str{j});
                        end
                        matchindex = find(ismember(matchsegfilesname,segfilesname));
                        if isempty(segfiles)
                            error('Cannot find segmentation files, please check the path to segmentation files is correct');
                        else
                            for j = 1:length(sample(matchindex))                                                   
                                fi_seg(:,:,matchindex(j)) = im2bw(imread([segDirectory matchsegfilesname{matchindex(j)}]),1e-6);
                            end 
                        end
                        [fi_base_temp,flat_amp_temp] = BasefluorEst_AW(fi,flatfield,darkfield,fi_seg);
                     else
                        [fi_base_temp,flat_amp_temp] = BasefluorEst_AW(fi,flatfield,darkfield);
                     end
                     fi_base(sample) = fi_base_temp;
                     flat_amp(sample) = flat_amp_temp;
                end
                if isempty(temporalRangeIII)
                    fileID = fopen([outputPositionDirectory 'basefluor_' char(wavelength) '.txt'],'w');
                else
                    fileID = fopen([outputPositionDirectory 'basefluor_' char(wavelength) '_t' num2str(temporalRangeIII(1)) 'to' num2str(temporalRangeIII(2)) '.txt'],'w');
                end
                fprintf(fileID,'meanbasefluor:%.8f\n', mean(fi_base));
                for p = 1:length(time_str)
                    fprintf(fileID,'%s:%.8f\n',char(time_str{p}),fi_base(p));
                end
                fclose(fileID);
                if isempty(temporalRangeIII)
                    fileID = fopen([outputPositionDirectory 'ratioflat_' char(wavelength) '.txt'],'w');
                else
                    fileID = fopen([outputPositionDirectory 'ratioflat_' char(wavelength) '_t' num2str(temporalRangeIII(1)) 'to' num2str(temporalRangeIII(2)) '.txt'],'w');
                end
                fprintf(fileID,'meanratioflat:%.8f\n', mean(flat_amp));
                for p = 1:length(time_str)
                    fprintf(fileID,'%s:%.8f\n',char(time_str{p}),flat_amp(p));
                end
                fclose(fileID);
                 figure(2); subplot(221);imagesc(flatfield);colorbar;
                 %saveas(gcf,[save_dir 'flatfield_' char(wavelength)],'png');
                 if ~isempty(darkfield)
                    subplot(222); imagesc(darkfield); colorbar;
                    %saveas(gcf,[save_dir 'darkfield_' char(wavelength)],'png');
                 end
                 subplot(223); plot(fi_base);
                 subplot(224); plot(flat_amp);
                 if isempty(temporalRangeIII)
                    figure(2); saveas(gcf,[outputPositionDirectory 'shading_' char(wavelength)],'tif');
                 else
                     figure(2); saveas(gcf,[outputPositionDirectory 'shading_' char(wavelength) '_t' num2str(temporalRangeIII(1)) 'to' num2str(temporalRangeIII(2))],'tif');
                 end
                     
                 %% run correction
                 if ~runCorrection
                 else
                     saveCorrectionDirectory = [outputPositionDirectory 'Correction/'];
                     mkdir(saveCorrectionDirectory);
                     if ~adjustFlatfield
                         flat_amp = ones(length(fi_base),1);
                     end
                     if ~floatOutput
                         for j = 1:length(valid_idx)
                            temp = strsplit(imagefiles(valid_idx(j)).name,'.');
                            correctedimagefile = [saveCorrectionDirectory temp{1}];
                            correction_modeLS(fi(:,:,j),correctedimagefile,flatfield,darkfield,fi_base(valid_idx(j)),flat_amp(valid_idx(j)));
                         end
                     else
                         for j = 1:length(valid_idx)
                            temp = strsplit(imagefiles(valid_idx(j)).name,'.');
                            correctedimagefile = [saveCorrectionDirectory temp{1}];
                            correction_mode2(fi(:,:,j),correctedimagefile,flatfield,darkfield,fi_base(valid_idx(j)),flat_amp(valid_idx(j)));
                         end                         
                     end
                 end
            end
        end
    end  
end

function correction_mode1(imag,saveimagename,flatfield,darkfield,fi_base,flat_amp)
imag_corr = (imag-darkfield)./(1.0+(flatfield-1.0)*flat_amp)-fi_base;
imag_corr(imag_corr<0) = 0; imag_corr(imag_corr>255) = 255;
imwrite(uint8(imag_corr),[saveimagename '.png']);

function correction_modeLS(imag,saveimagename,flatfield,darkfield,fi_base,flat_amp)
imag_corr = (imag-darkfield)./(1.0+(flatfield-1.0)*flat_amp)-fi_base;
% imag_corr = (imag-darkfield)./flatfield -fi_base
imag_corr = imag_corr.*10000; imag_corr = imag_corr+10000; imag_corr(imag_corr<0) = 0; imag_corr(imag_corr>(2^16-1)) = 2^16-1;
imwrite(uint16(imag_corr),[saveimagename '.png']);

function correction_mode2(imag,saveimagename,flatfield,darkfield,fi_base,flat_amp)
imag_corr = (imag-darkfield)./(1.0+(flatfield-1.0)*flat_amp)-fi_base;
saveToTif(imag_corr, saveimagename);

function saveToTif(imag,saveimagename)
t = Tiff([saveimagename '.tif'], 'w'); 
tagstruct.ImageLength = size(imag, 1); 
tagstruct.ImageWidth = size(imag, 2); 
tagstruct.Compression = Tiff.Compression.None; 
tagstruct.SampleFormat = Tiff.SampleFormat.IEEEFP; 
tagstruct.Photometric = Tiff.Photometric.MinIsBlack; 
tagstruct.BitsPerSample = 32; 
tagstruct.SamplesPerPixel = 1; 
tagstruct.PlanarConfiguration = Tiff.PlanarConfiguration.Chunky; 
t.setTag(tagstruct); 
t.write(single(imag)); 
t.close();



