Contents
Clear Variables. Turn off Obnoxious warnings...
if S.GUI_flag == 1
else
clear all;
close all;
S.GUI_flag = 0;
end
warning('off','MATLAB:nearlySingularMatrix');
warning('off','MATLAB:polyfit:PolyNotUnique');
warning('off','Images:initSize:adjustingMag');
warning('off','curvefit:fit:equationBadlyConditioned');
warning('off','curvefit:prepareSurfaceData:removingNaNAndInf');
warning('off','curvefit:prepareSurfaceData:swapXAndYForTable');
warning('off','MATLAB:MKDIR:DirectoryExists');
Undefined function or variable 'S'.
Error in ==> Nano6_Flatten_surf_poly_GUI at 6
if S.GUI_flag == 1 % program was called from the GUI
Set relevant input paramters. Scalebar sizes, fitting limits and so forth
ScaleBar = 1000;
if S.GUI_flag == 1
FontSize = 11;
FontName = 'Arial';
else
FontSize = 16;
FontName = 'Times';
end
print_dpi = 150;
amplitude_fract = 0.95;
crop_fraction = 0.0;
channels_nb = 0;
if S.GUI_flag
existingChannels = zeros(14, 1);
for k = 1:7
channel_name = get(S.channelProcessText(k), 'string');
channel_name = channel_name{1};
switch channel_name
case 'Height'
existingChannels(1) = 1;
case 'Deflection'
existingChannels(2) = 1;
case 'Phase'
existingChannels(3) = 1;
case 'Amplitude'
existingChannels(4) = 1;
case 'Amplitude Error'
existingChannels(5) = 1;
case 'Potential'
existingChannels(6) = 1;
case 'Peak Force Error'
existingChannels(7) = 1;
case 'DMTModulus'
existingChannels(8) = 1;
case 'LogDMTModulus'
existingChannels(9) = 1;
case 'Adhesion'
existingChannels(10) = 1;
case 'Deformation'
existingChannels(11) = 1;
case 'Dissipation'
existingChannels(12) = 1;
case 'Input1'
existingChannels(13) = 1;
case 'Input2'
existingChannels(14) = 1;
end
end
channels_indices = find(existingChannels);
for k = 1:14
Med_corr(k) = get(S.Med_corr(k), 'value');
Med_corr_2(k) = get(S.Med_corr_2(k), 'value');
vert_med_corr(k) = get(S.vert_med_corr(k), 'value');
poly_order_horiz(k) = get(S.poly_order_horiz(k), 'value')-1;
poly_order_vert(k) = get(S.poly_order_vert(k), 'value')-1;
z_height_reference(k) = str2num(get(S.z_height_reference(k), 'string'));
min_spacing_stdev_user(k) = str2num(get(S.min_spacing_stdev_user(k), 'string'));
max_spacing_stdev_user(k) = str2num(get(S.max_spacing_stdev_user(k), 'string'));
min_z(k) = str2num(get(S.min_z(k), 'string'));
max_z(k) = str2num(get(S.max_z(k), 'string'));
peaks_fit(k) = get(S.peaks_fit(k), 'value')-1;
peak_reference_height(k) = get(S.peak_reference_height(k), 'value')-1;
peak_to_flatten(k) = get(S.peak_to_flatten(k), 'value')-1;
Calc_flatten_lim(k) = get(S.Calc_flatten_lim(k), 'value');
overwrite_xlim(k) = get(S.overwrite_xlim(k), 'value');
sort_peaks(k) = get(S.sort_peaks(k), 'value');
end
Med_corr = Med_corr(existingChannels~=0);
Med_corr_2 = Med_corr_2(existingChannels~=0);
vert_med_corr = vert_med_corr(existingChannels~=0);
poly_order_horiz = poly_order_horiz(existingChannels~=0);
poly_order_vert = poly_order_vert(existingChannels~=0);
peaks_fit = peaks_fit(existingChannels~=0);
peak_reference_height = peak_reference_height(existingChannels~=0);
peak_to_flatten = peak_to_flatten(existingChannels~=0);
Calc_flatten_lim = Calc_flatten_lim(existingChannels~=0);
overwrite_xlim = overwrite_xlim(existingChannels~=0);
sort_peaks = sort_peaks(existingChannels~=0);
z_height_reference = z_height_reference(existingChannels~=0);
min_spacing_stdev_user = min_spacing_stdev_user(existingChannels~=0);
max_spacing_stdev_user = max_spacing_stdev_user(existingChannels~=0);
min_z = min_z(existingChannels~=0);
max_z = max_z(existingChannels~=0);
for k = 1:7
channels_flatten(k) = get(S.channelProcessSelect(k), 'value');
if channels_flatten(k) == 1
channels_nb = channels_nb + 1;
end
end
Med_corr = Med_corr(channels_flatten~=0);
Med_corr_2 = Med_corr_2(channels_flatten~=0);
vert_med_corr = vert_med_corr(channels_flatten~=0);
poly_order_horiz = poly_order_horiz(channels_flatten~=0);
poly_order_vert = poly_order_vert(channels_flatten~=0);
peaks_fit = peaks_fit(channels_flatten~=0);
peak_reference_height = peak_reference_height(channels_flatten~=0);
peak_to_flatten = peak_to_flatten(channels_flatten~=0);
Calc_flatten_lim = Calc_flatten_lim(channels_flatten~=0);
overwrite_xlim = overwrite_xlim(channels_flatten~=0);
sort_peaks = sort_peaks(channels_flatten~=0);
z_height_reference = z_height_reference(channels_flatten~=0);
min_spacing_stdev_user = min_spacing_stdev_user(channels_flatten~=0);
max_spacing_stdev_user = max_spacing_stdev_user(channels_flatten~=0);
min_z = min_z(channels_flatten~=0);
max_z = max_z(channels_flatten~=0);
channels_indices = channels_indices(channels_flatten~=0);
Offset_data = zeros(channels_nb, 1);
fit_hist = zeros(channels_nb, 1);
for k = 1:channels_nb
if z_height_reference(k) == 0
Offset_data(k) = 0;
else
Offset_data(k) = 1;
end
if peaks_fit(k) == 0
fit_hist(k) = 0;
else
fit_hist(k) = 1;
end
end
histogram_bins = str2num(get(S.histogram_bins, 'string'));
smoothing_level = str2num(get(S.smoothing_level, 'string'));
thresh_range = str2num(get(S.thresh_range, 'string'));
min_peak_height = str2num(get(S.min_peak_height, 'string'));
min_peak_distance = str2num(get(S.min_peak_distance, 'string'));
peak_threshold = str2num(get(S.peak_threshold, 'string'));
else
Med_corr = [1,0,0,1,1,1,1];
Med_corr_2 = [1,0,0,1,1,1,1];
vert_med_corr = [0,0,0,0,0,0,0];
initial_poly_order_horiz = [2, 0, 0, 0, 0, 0, 0];
initial_poly_order_vert = [2, 0, 0, 0, 0, 0, 0];
poly_order_horiz = [3, 0, 0, 0, 0, 1, 1];
poly_order_vert = [3, 0, 0, 0, 0, 1, 1];
Offset_data = [1, 0, 0, 0, 0, 0, 0];
histogram_bins = 256;
smoothing_level = 6;
z_height_reference = [ 5.0e-0, 0.0e-0, 0.0e-0, 0.0e-0, 0.0e-0, 0.0e-0, 0.0e-0];
Calc_flatten_lim = [1, 1, 1, 1, 1, 1, 1];
thresh_range = 3;
min_spacing_stdev_user = [ 5.0e-2, 2.0e+0, 2.0e-1, 5.0e-2, 1.0e+1, 5.0e-3, 1.0e-3];
max_spacing_stdev_user = [ 5.0e-1, 2.0e+1, 2.0e+0, 1.0e-0, 4.0e+1, 2.0e-2, 5.0e-2];
overwrite_xlim = [1, 1, 1, 1, 1, 1, 1];
min_z = [0.0e+0, -3.0e+1, -5.0e+0, -1.0e+1, -1.0e+2, -1.0e-2, -1.0e-1];
max_z = [ 10.0e+0, 3.0e+1, 5.0e+0, 1.0e+1, 1.0e+2, 3.0e-2, 1.0e-1];
channels_flatten = [1 0 0 0 0 0 0];
channels_indices = find(channels_flatten);
channels_nb = length(find(channels_flatten));
fit_hist = [1, 1, 1, 1, 1, 1, 1];
peaks_fit = [2,1,1,2,1,2,1];
peak_reference_height = [2, 1, 1, 1, 1, 1, 1];
peak_to_flatten = [2, 1, 1, 1, 1, 1, 1];
min_peak_height = 1e2;
min_peak_distance = 5;
peak_threshold = 5;
sort_peaks = [1, 1, 1, 1, 1, 1, 1];
end
disp_percent = 1;
percent_search = [1 10 90 99];
sequential_flattening = 1;
plot_intermediate_steps = 0;
plot_intermediate_image = 0;
peak_pos_tolerance = 5.0e-1;
Select Files to read
if S.GUI_flag
Filename = evalin('base', 'Filename');
ImageDir = get(S.selectFolderText, 'string');
ImageDir = [ImageDir,filesep];
if S.processAll == 1
Filename_length = length(Filename);
else
Filename = Filename{1};
Filename = {Filename, Filename};
Filename_length = 1;
end
else
[Filename, ImageDir] = uigetfile( {'*.*','Nanoscope Files (*.*)'}, ...
'Pick a file',...
'C:\Users\LBNI\Desktop\Test ata\',...
'MultiSelect','on');
Filename_length = length(Filename);
end
Prep for file output
mkdir([ImageDir,'Matlab',filesep,'Figure Jpegs',filesep]);
mkdir([ImageDir,'Matlab',filesep,'Figure Fig',filesep]);
Loop the import and processing routiene for the number of selected files
load('Sky','mycmap_sky')
searchstring(1).label='@2:Z scale:';
searchstring(2).label='Samps';
searchstring(3).label='Lines:';
searchstring(4).label='Data offset';
searchstring(5).label='Image Data:';
searchstring(6).label='Scan Size:';
searchstring(7).label='@Sens. Zsens:';
searchstring(8).label='@Z magnify:';
searchstring(9).label='Scan Rate:';
searchstring(10).label='@Sens. DeflSens: V';
searchstring(11).label='@2:CantFrequency: V';
searchstring(12).label='@2:CantDrive: V';
searchstring(13).label='@2:TMSetAmplitude: V';
Init_poly_coef = cell(length(Filename),length(channels_flatten));
for index = 1:Filename_length
file_name = [ImageDir,Filename{index}];
Read in DI Image and image propreties This code is mostly from Jaco De Groot at University College London
fid = fopen(file_name,'r');
[message,errnum] = ferror(fid);
if(errnum)
fprintf(1,'I/O Error %d \t %s',[errnum,message]);
end
header_end=0; eof = 0; counter = 1; byte_location = 0;
nstrings=size(searchstring,2);
parcounter = ones(nstrings,1);
paramters = struct('trace',cell(1,nstrings),'values',cell(1,nstrings),'channel',cell(1,nstrings));
for ij=1:nstrings
parcounter(ij)=1;
parameters(ij).trace=0;
end;
while( and( ~eof, ~header_end ) )
byte_location = ftell(fid);
line = fgets(fid);
for ij=1:nstrings
if strfind(line,searchstring(ij).label)
if (extract_num(line))
b=strfind(line, 'LSB');
if (b>0)
parameters(ij).values(parcounter(ij))=extract_num(line(b(1):end));
else parameters(ij).values(parcounter(ij))=extract_num(line);
end;
else
b= strfind(line,'"');
if (b>0)
parameters(ij).channel(parcounter(ij)).name=line(b(1)+1:b(2)-1);
else
if (strfind(line,'Trace')>0)
parameters(ij).trace=1;
end;
end;
end;
parcounter(ij)=parcounter(ij)+1;
end
end;
if( (-1)==line )
eof = 1;
end
if ~isempty( strfind( line, '\*File list end' ) )
header_end = 1;
end
counter=counter+1;
end
fclose(fid);
param = parameters;
scaling = param(1).values;
spl = param(2).values;
linno = param(3).values;
image_pos = param(4).values;
ScanSize = param(6).values(1);
Z_Sensitivity = param(7).values;
Z_Magnification = param(8).values;
Scan_Rate = param(9).values(1);
Defl_Sens = param(10).values;
Drive_Freq = param(11).values;
Drive_Amp = param(12).values;
Amp_Set_Point = param(13).values;
[tok remain] = strtok(Filename{index},'.');
file_number = str2double(remain(2:4));
Image_Data_File = [ImageDir,'Matlab',filesep,'Imaging_Parameters.txt'];
if index == 1
fid = fopen(Image_Data_File, 'a');
fprintf(fid, [ 'File Number','\t',...
'Samples per Line','\t',...
'Lines per Image','\t',...
'Scan Rate (Hz)','\t',...
'Deflection Sensitivity (nm/V)','\t',...
'Drive Frequency (kHz)','\t',...
'Drive Amplitude (mV)','\t',...
'Amplitude Setpoint (mV)','\n'],...
'char');
fprintf(fid, [ num2str(file_number,'%0.0f'),'\t',...
num2str(spl(1),'%0.0f'),'\t',...
num2str(linno,'%0.0f'),'\t',...
num2str(Scan_Rate,'%0.3f'),'\t',...
num2str(Defl_Sens,'%0.3f'),'\t',...
num2str(Drive_Freq,'%0.2f'),'\t',...
num2str(Drive_Amp,'%0.2f'),'\t',...
num2str(Amp_Set_Point,'%0.2f'),'\n'],...
'char');
fclose(fid);
else
fid = fopen(Image_Data_File, 'a');
fprintf(fid, [ num2str(file_number,'%0.0f'),'\t',...
num2str(spl(1),'%0.0f'),'\t',...
num2str(linno,'%0.0f'),'\t',...
num2str(Scan_Rate,'%0.3f'),'\t',...
num2str(Defl_Sens,'%0.3f'),'\t',...
num2str(Drive_Freq,'%0.2f'),'\t',...
num2str(Drive_Amp,'%0.2f'),'\t',...
num2str(Amp_Set_Point,'%0.2f'),'\n'],...
'char');
fclose(fid);
end
L = length(image_pos);
channel_info = struct('Trace',cell(L,1),'Name',cell(L,1),'Finalscaling',cell(L,1),'Unit',cell(L,1));
finalscaling = zeros(L,1);
for im=1:L
channel_info(im).Trace=char((param(6).trace)*'Trace '+(1-param(6).trace)*'Retrace');
channel_info(im).Name=param(5).channel(im).name;
channel_info(im).Finalscaling(im)=param(8).values(im)*param(1).values(im);
mkdir([ImageDir,'Matlab',filesep,'Figure Jpegs',filesep,channel_info(im).Name,filesep]);
mkdir([ImageDir,'Matlab',filesep,'Figure Jpegs',filesep,channel_info(im).Name,filesep,'hist',filesep]);
mkdir([ImageDir,'Matlab',filesep,'Figure Jpegs',filesep,channel_info(im).Name,filesep,'grey',filesep]);
switch channel_info(im).Name
case 'Height'
channel_info(im).Unit='nm' ;
channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);
finalscaling(im)=(scaling(im)*Z_Sensitivity)/(2^16);
case 'Deflection'
channel_info(im).Unit='V';
finalscaling(im)=(scaling(im))/(2^16);
case 'Phase'
channel_info(im).Unit='Degree';
finalscaling(im)=(scaling(im))/(2^16);
case 'Amplitude'
channel_info(im).Unit='mV';
channel_info(im).Finalscaling(im)=param(8).values(im)*param(1).values(im);
finalscaling(im)=(scaling(im))/(2^16);
case 'Amplitude Error'
channel_info(im).Unit='mV';
channel_info(im).Finalscaling(im)=param(8).values(im)*param(1).values(im);
finalscaling(im)=(scaling(im))/(2^16);
case 'Potential'
channel_info(im).Unit='V';
channel_info(im).Finalscaling(im)=param(8).values(im)*param(1).values(im);
finalscaling(im)=(scaling(im))/(2^16);
case 'Peak Force Error'
channel_info(im).Unit='nN' ;
channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);
finalscaling(im)=(scaling(im)*Z_Sensitivity)/(2^16);
case 'DMTModulus'
channel_info(im).Unit='MPa' ;
channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);
finalscaling(im)=(scaling(im)*Z_Sensitivity)/(2^16);
case 'LogDMTModulus'
channel_info(im).Unit='log(MPa)' ;
channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);
finalscaling(im)=(scaling(im)*Z_Sensitivity)/(2^16);
case 'Adhesion'
channel_info(im).Unit='nN ' ;
channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);
finalscaling(im)=(scaling(im)*Z_Sensitivity)/(2^16);
case 'Deformation'
channel_info(im).Unit='nm' ;
channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);
finalscaling(im)=(scaling(im)*Z_Sensitivity)/(2^16);
case 'Dissipation'
channel_info(im).Unit='eV' ;
channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);
finalscaling(im)=(scaling(im)*Z_Sensitivity)/(2^16);
case 'Input1'
channel_info(im).Unit='nm' ;
channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);
finalscaling(im)=(scaling(im)*Z_Sensitivity*Defl_Sens*1.50)/(2^16);
case 'Input2'
channel_info(im).Unit='nm' ;
channel_info(im).Finalscaling(im)=param(7).values*param(8).values(im)*param(1).values(im);
finalscaling(im)=(scaling(im)*Z_Sensitivity*22*Defl_Sens)/(2^16);
otherwise
channel_info(im).Unit='V';
channel_info(im).Finalscaling(im)=param(8).values(im)*param(1).values(im);
finalscaling(im)=(scaling(im))/(2^16);
end;
end;
if S.GUI_flag
if S.processAll == 0
disp({['Image ',num2str(index),' of 1'];['Z Sensitivity ',num2str(Z_Sensitivity)];['Hard Scaling ',num2str(scaling)];['Final Scaling ',num2str(finalscaling(1))]});
else
disp({['Image ',num2str(index),' of ',num2str(length(Filename))];['Z Sensitivity ',num2str(Z_Sensitivity)];['Hard Scaling ',num2str(scaling)];['Final Scaling ',num2str(finalscaling(1))]});
end
else
disp({['Image ',num2str(index),' of ',num2str(length(Filename))];['Z Sensitivity ',num2str(Z_Sensitivity)];['Hard Scaling ',num2str(scaling)];['Final Scaling ',num2str(finalscaling(1))]});
end
fid = fopen(file_name,'r');
images = zeros(linno, spl(1), L);
for i = 1:L
fseek(fid,image_pos(i),-1);
A = fread(fid, [spl(i) linno],'int16');
images(:,:,i) = rot90(finalscaling(i)*A);
end;
Crop Images. Helps deal with some of the turn around ripple
images = images(ceil(crop_fraction/2*size(images,1))+1:floor((1-crop_fraction/2)*size(images,1))-1,ceil(crop_fraction/2*size(images,2))+1:floor((1-crop_fraction/2)*size(images,2))-1,:);
images = images(:,:, channels_flatten == 1);
channel_info = channel_info(channels_flatten == 1);
This section actuall does all of the image processing. Different blocks will be explained below Code Loops through each of the channels and applies all of the corrections
for j = 1:channels_nb
disp({'Channel:'; channel_info(j).Name})
try
if S.GUI_flag
if S.processAll == 0 && j == 1
raw_images = zeros(size(images,1), size(images,2), 14);
med_corr_1_images = zeros(size(images,1), size(images,2), 14);
poly_flat_1_images = zeros(size(images,1), size(images,2), 14);
med_corr_2_images = zeros(size(images,1), size(images,2), 14);
poly_flat_2_images = zeros(size(images,1), size(images,2), 14);
final_images = zeros(size(images,1), size(images,2), 14);
end
end
if S.GUI_flag
if S.processAll == 0
raw_images(:,:,channels_indices(j)) = images(:,:,j);
end
end
max_peak_ampl = 0;
ScanZ = max_z(j)-min_z(j);
if plot_intermediate_steps == 1
if index == 1 && j == channels_flatten(1)
summary_fig = figure('Name','Summary', 'OuterPosition', [1 1 1800 800]);
else
close(summary_fig)
summary_fig = figure('Name','Summary', 'OuterPosition', [1 1 1800 800]);
end
subplot(2, 7, 1);
plot_nice_image;
title('Initial image')
subplot(2, 7, 8);
hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins)
title ('Initial histogram')
end
[height_ini position_ini] = hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins);
hist_std_ini = std(position_ini);
Does an initial Median correction to fix basic line skips. Don't use on grids or other similar structures
saved_image_no_med_corr = images(:,:,j);
if Med_corr(j) == 1
disp('--- First median correction ----')
[images(:,:,j), hist_std, improv_param] = medianCorrection(images(:,:,j), hist_std_ini, 4, histogram_bins);
if S.GUI_flag
if S.processAll == 0
med_corr_1_images(:,:,channels_indices(j)) = images(:,:,j);
end
end
if plot_intermediate_steps == 1
figure(summary_fig)
subplot(2, 7, 2);
plot_nice_image;
title('After 1st median correction')
subplot(2, 7, 9);
hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins)
title('Histogram')
if (improv_param == 1)
xlabel('This step improved the image');
else
xlabel('This step did NOT improve the image');
end
end
else
[height_ini position_ini] = hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins);
hist_std = std(position_ini);
end
if vert_med_corr(j) == 1
[images(:,:,j), hist_std, improv_param] = medianCorrection(images(:,:,j)', hist_std, 4, histogram_bins);
images(:,:,j) = images(:,:,j)';
end
final_flattening_image = images(:,:,j);
Does an initial tilt correction to fix some of the initial sample tilt. Needed before doing any histogrammed flattening % Determines the number of coeffecients for the user selected polynomial order
saved_image = images(:,:,j);
initial_poly_order_horizontal = 1;
initial_poly_order_vertical = 1;
coef_length = nchoosek(max([initial_poly_order_horizontal,initial_poly_order_vertical])+2,2);
[xInput, yInput, zOutput] = prepareSurfaceData(1:size(images(:,:,j),1),1:size(images(:,:,j),2) , images(:,:,j));
poly_type = fittype( ['poly',num2str(initial_poly_order_vertical),num2str(initial_poly_order_horizontal)] );
poly_opts = fitoptions( poly_type );
if index == 1
poly_opts.Lower = -Inf*ones(coef_length,1);
poly_opts.Upper = Inf*ones(coef_length,1);
else
poly_opt.Lower = coeffvalues(Init_poly_coef{index-1,j})-0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
poly_opt.Upper = coeffvalues(Init_poly_coef{index-1,j})+0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
end
[poly_result] = fit( [xInput, yInput], zOutput, poly_type, poly_opts );
images(:,:,j) = images(:,:,j) - reshape(poly_result([xInput,yInput]),size(images(:,:,j),1),size(images(:,:,j),2));
Init_poly_coef{index,j} = poly_result;
saved_image_1st_order_fit_med_cor = images(:,:,j);
images(:,:,j) = saved_image_no_med_corr;
initial_poly_order_horizontal = 1;
initial_poly_order_vertical = 1;
coef_length = nchoosek(max([initial_poly_order_horizontal,initial_poly_order_vertical])+2,2);
[xInput, yInput, zOutput] = prepareSurfaceData(1:size(images(:,:,j),1),1:size(images(:,:,j),2) , images(:,:,j));
poly_type = fittype( ['poly',num2str(initial_poly_order_vertical),num2str(initial_poly_order_horizontal)] );
poly_opts = fitoptions( poly_type );
if index == 1
poly_opts.Lower = -Inf*ones(coef_length,1);
poly_opts.Upper = Inf*ones(coef_length,1);
else
poly_opt.Lower = coeffvalues(Init_poly_coef{index-1,j})-0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
poly_opt.Upper = coeffvalues(Init_poly_coef{index-1,j})+0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
end
[poly_result] = fit( [xInput, yInput], zOutput, poly_type, poly_opts );
images(:,:,j) = images(:,:,j) - reshape(poly_result([xInput,yInput]),size(images(:,:,j),1),size(images(:,:,j),2));
Init_poly_coef{index,j} = poly_result;
saved_image_1st_order_fit = images(:,:,j);
images(:,:,j) = saved_image;
initial_poly_order_horizontal = 2;
initial_poly_order_vertical = 2;
coef_length = nchoosek(max([initial_poly_order_horizontal,initial_poly_order_vertical])+2,2);
[xInput, yInput, zOutput] = prepareSurfaceData(1:size(images(:,:,j),1),1:size(images(:,:,j),2) , images(:,:,j));
poly_type = fittype( ['poly',num2str(initial_poly_order_vertical),num2str(initial_poly_order_horizontal)] );
poly_opts = fitoptions( poly_type );
if index == 1
poly_opts.Lower = -Inf*ones(coef_length,1);
poly_opts.Upper = Inf*ones(coef_length,1);
else
poly_opt.Lower = coeffvalues(Init_poly_coef{index-1,j})-0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
poly_opt.Upper = coeffvalues(Init_poly_coef{index-1,j})+0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
end
[poly_result] = fit( [xInput, yInput], zOutput, poly_type, poly_opts );
images(:,:,j) = images(:,:,j) - reshape(poly_result([xInput,yInput]),size(images(:,:,j),1),size(images(:,:,j),2));
Init_poly_coef{index,j} = poly_result;
saved_image_2nd_order_fit_med_cor = images(:,:,j);
images(:,:,j) = saved_image_no_med_corr;
initial_poly_order_horizontal = 2;
initial_poly_order_vertical = 2;
coef_length = nchoosek(max([initial_poly_order_horizontal,initial_poly_order_vertical])+2,2);
[xInput, yInput, zOutput] = prepareSurfaceData(1:size(images(:,:,j),1),1:size(images(:,:,j),2) , images(:,:,j));
poly_type = fittype( ['poly',num2str(initial_poly_order_vertical),num2str(initial_poly_order_horizontal)] );
poly_opts = fitoptions( poly_type );
if index == 1
poly_opts.Lower = -Inf*ones(coef_length,1);
poly_opts.Upper = Inf*ones(coef_length,1);
else
poly_opt.Lower = coeffvalues(Init_poly_coef{index-1,j})-0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
poly_opt.Upper = coeffvalues(Init_poly_coef{index-1,j})+0.5*abs(coeffvalues(Init_poly_coef{index-1,j}));
end
[poly_result] = fit( [xInput, yInput], zOutput, poly_type, poly_opts );
images(:,:,j) = images(:,:,j) - reshape(poly_result([xInput,yInput]),size(images(:,:,j),1),size(images(:,:,j),2));
Init_poly_coef{index,j} = poly_result;
saved_image_2nd_order_fit = images(:,:,j);
tilt_corr_names = ['saved_image_1st_order_fit '; ...
'saved_image_2nd_order_fit '; ...
'saved_image_1st_order_fit_med_cor'; ...
'saved_image_2nd_order_fit_med_cor'];
for ii = 1:4
eval(['images(:,:,j) = ', num2str(tilt_corr_names(ii,:)), ';']);
[hcdf xcdf] = ecdf(reshape(images(:,:,j),size(images,1)*size(images,2),1));
k = 1;
while hcdf(k) <= 0.01
k = k+1;
end
x_1percent = xcdf(k);
index_1percent = floor((x_1percent-min(min(images(:,:,j))))/(max(max(images(:,:,j)))-min(min(images(:,:,j)))) * histogram_bins);
index_1percent = max([1,index_1percent]);
k = length(hcdf);
while hcdf(k) >= 0.99
k = k-1;
end
x_99percent = xcdf(k);
index_99percent = floor((x_99percent-min(min(images(:,:,j))))/(max(max(images(:,:,j)))-min(min(images(:,:,j)))) * histogram_bins);
index_99percent = min([index_99percent,length(xcdf)]);
[height position] = hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins);
height = smooth(height,smoothing_level,'loess');
[peak_height peak_index]= findpeaks(height,'MinpeakDistance', min_peak_distance, 'Sortstr','descend');
peak_index = peak_index(1:min(length(peak_index),peaks_fit(j)));
peak_height = height(peak_index);
tilt_hist_std(ii) = std(position(index_1percent:index_99percent));
end
if plot_intermediate_steps == 1
figure(summary_fig)
end
disp('--- Tilt correction ----')
if (min(tilt_hist_std) > hist_std)
images(:,:,j) = saved_image;
disp('Neither the 1st order nor the 2nd order inital tilt correction did improve the image')
if plot_intermediate_steps == 1
subplot(2, 7, 10);
xlabel('This step did NOT improve the image')
end
else
[minimum_std best_tilt_corr] = min(tilt_hist_std);
images_size = size(saved_image_1st_order_fit, 1);
all_tilt_corr_im = [saved_image_1st_order_fit; saved_image_2nd_order_fit; saved_image_1st_order_fit_med_cor; saved_image_2nd_order_fit_med_cor];
images(:,:,j) = all_tilt_corr_im(images_size*(best_tilt_corr-1)+1:images_size*best_tilt_corr, :);
titles = ['tilt 1st order without med cor'; 'tilt 2nd order without med cor'; 'tilt 1st order with med cor'; 'tilt 2nd order with med cor'];
disp(['best tilt correction: ', titles(best_tilt_corr, :)])
if plot_intermediate_steps == 1
subplot(2, 7, 10);
xlabel(['best corr: ', titles(best_tilt_corr, :)])
end
end
if plot_intermediate_steps == 1
figure(summary_fig)
subplot(2, 7, 3);
plot_nice_image;
title('After inital tilt correction')
subplot(2, 7, 10);
hold on;
end
images(:,:,j) = images(:,:,j)-median(median(images(:,:,j)));
min_hist = min(min(images(:,:,j))) + 0.05*min(min(images(:,:,j)));
max_hist = max(max(images(:,:,j))) - 0.05*max(max(images(:,:,j)));
if S.GUI_flag
if S.processAll == 0
ini_tilt_images(:,:,channels_indices(j)) = images(:,:,j);
end
end
Generates initial fitting bound estimates
image_inside_range = images(:,:,j);
image_inside_range((images(:,:,j) < min_hist)) = NaN;
image_inside_range((images(:,:,j) > max_hist)) = NaN;
[height_est position_est] = hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);
height_est = smooth(height_est,smoothing_level,'loess');
[peak_height_est peak_index_est] = findpeaks(height_est,'MinpeakDistance', min_peak_distance, 'MinpeakHeight', min_peak_height, 'Threshold', peak_threshold, 'Sortstr','descend');
if sort_peaks(j) == 1
peak_index_est = sort(peak_index_est(1:min([peaks_fit(j),3, length(peak_height_est)])),'ascend');
else
peak_index_est = peak_index_est(1:min([peaks_fit(j),3, length(peak_height_est)]));
end
peak_height_est = height_est(peak_index_est);
clear std_est
if length(peak_index_est) == 1
std_est(1) = std(reshape(images(:,:,j),size(images,1)*size(images,2),1));
elseif length(peak_index_est) == 2
std_est(1) = std( images(images(:,:,j) <= position_est(round((peak_index_est(1)+peak_index_est(2))/2))));
std_est(2) = std( images(images(:,:,j) >= position_est(round((peak_index_est(1)+peak_index_est(2))/2))));
elseif length(peak_index_est) == 3
std_est(1) = std( images(images(:,:,j) <= position_est(round((peak_index_est(1)+peak_index_est(2))/2))));
std_est(2) = std( images(position_est(round((peak_index_est(1)+peak_index_est(2))/2)) <= images(:,:,j) <= position_est(round((peak_index_est(2)+peak_index_est(3))/2))));
std_est(3) = std( images(images(:,:,j) >= position_est(round((peak_index_est(2)+peak_index_est(3))/2))));
else
end
if Calc_flatten_lim(j) == 1
max_spacing_stdev(j) = 1.2*max(std_est);
min_spacing_stdev(j) = 0.01*min(std_est);
else
min_spacing_stdev = min_spacing_stdev_user;
max_spacing_stdev = max_spacing_stdev_user;
end
if plot_intermediate_steps == 1
figure(summary_fig)
subplot(2, 7, 10);
title('Histogram')
image_inside_range = images(:,:,j);
image_inside_range((images(:,:,j) < min_hist)) = NaN;
image_inside_range((images(:,:,j) > max_hist)) = NaN;
hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);
end
if plot_intermediate_image == 1
save_intermediate_image;
end
[peak_height_est peak_index_est] = findpeaks(height_est,'MinpeakDistance', min_peak_distance, 'MinpeakHeight', min_peak_height, 'Threshold', peak_threshold, 'Sortstr','descend');
peak_index_est = peak_index_est(1:min([peaks_fit(j),3, length(peak_height_est)]));
ref_peak_pos = position_est(peak_index_est);
if plot_intermediate_steps == 1
figure(summary_fig)
subplot(2, 7, 10)
plot(position_est(peak_index_est), height_est(peak_index_est), 'k^', 'MarkerEdgeColor','k', 'MarkerFaceColor', 'g')
hold off;
end
This section histograms the image, identifies up to 3 distinct levels in the historam, fits them to Gaussians, and uses this to determin an approriate threshold range to use for a masked polynomial approach. This leads to truly flat planes that are not biased by steps and terraces in the image data
saved_image = images(:,:,j);
image_inside_range = images(:,:,j);
image_inside_range((images(:,:,j) < min_hist)) = NaN;
image_inside_range((images(:,:,j) > max_hist)) = NaN;
[height_1 position_1] = hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);
height_1 = smooth(height_1,smoothing_level,'loess');
[peak_height_1 peak_index_1] = findpeaks(height_1,'MinpeakDistance', min_peak_distance, 'MinpeakHeight', min_peak_height, 'Threshold', peak_threshold, 'Sortstr','descend');
peak_index_1 = find_closest(ref_peak_pos, peak_index_1(1:min(7, length(peak_index_1))), position_1);
if sort_peaks(j) == 1
peak_index_1 = sort(peak_index_1(1:min(peaks_fit(j), length(peak_index_1))),'ascend');
else
peak_index_1 = peak_index_1(1:min(peaks_fit(j), length(peak_index_1)));
end
peak_height_1 = height_1(peak_index_1);
func_form = [];
t_hist_lower_a = zeros(1,min(3,length(peak_height_1)));
t_hist_lower_s = zeros(1,min(3,length(peak_height_1)));
t_hist_lower_x = zeros(1,min(3,length(peak_height_1)));
t_hist_upper_a = zeros(1,min(3,length(peak_height_1)));
t_hist_upper_s = zeros(1,min(3,length(peak_height_1)));
t_hist_upper_x = zeros(1,min(3,length(peak_height_1)));
t_hist_start_a = zeros(1,min(3,length(peak_height_1)));
t_hist_start_s = zeros(1,min(3,length(peak_height_1)));
t_hist_start_x = zeros(1,min(3,length(peak_height_1)));
for i = 1:min(3,length(peak_height_1))
t_hist_lower_a(i) = (1-amplitude_fract)*peak_height_1(i);
t_hist_lower_s(i) = min_spacing_stdev(j);
t_hist_lower_x(i) = position_1(peak_index_1(i))-3.5/histogram_bins;
t_hist_upper_a(i) = (1+amplitude_fract)*peak_height_1(i);
t_hist_upper_s(i) = max_spacing_stdev(j);
t_hist_upper_x(i) = position_1(peak_index_1(i))+3.5/histogram_bins;
t_hist_start_a(i) = 1.00*peak_height_1(i);
t_hist_start_s(i) = (min_spacing_stdev(j)+max_spacing_stdev(j))/2;
t_hist_start_x(i) = position_1(peak_index_1(i));
func_form = [func_form,'a',num2str(i),'*exp(-(x-x',num2str(i),')^2/(2*(s',num2str(i),')^2))'];
if i < min(3,length(peak_height_1))
func_form = [func_form,'+'];
end
end
t_hist_lower = [t_hist_lower_a,t_hist_lower_s,t_hist_lower_x];
t_hist_upper = [t_hist_upper_a,t_hist_upper_s,t_hist_upper_x];
t_hist_start = [t_hist_start_a,t_hist_start_s,t_hist_start_x];
t_hist = fitoptions('Method','NonlinearLeastSquares',...
'Lower', t_hist_lower,...
'Upper', t_hist_upper,...
'Startpoint', t_hist_start,...
'MaxIter', 1e2,...
'MaxFunEvals',3e2);
f_Number_hist = fittype(func_form, 'options' ,t_hist);
[peak_1_gaussian,peak_1_gaussian_gof] = fit( position_1',...
height_1,...
f_Number_hist);
peaks_nb = min(3,length(peak_height_1));
for i = 1:min(3,length(peak_height_1))
eval(['peak_ampl = peak_1_gaussian.a', num2str(i), ';']);
if peak_ampl < min_peak_height
peaks_nb = peaks_nb - 1;
end
end
peak_1_gaussian_ampl = [];
peak_1_gaussian_pos = [];
peak_1_gaussian_std = [];
for i = 1:peaks_nb
eval(['peak_1_gaussian_ampl(', num2str(i), ')= peak_1_gaussian.a', num2str(i), ';']);
eval(['peak_1_gaussian_pos(', num2str(i), ')= peak_1_gaussian.x', num2str(i), ';']);
eval(['peak_1_gaussian_std(', num2str(i), ')= peak_1_gaussian.s', num2str(i), ';']);
end
if peaks_nb == 1
peak_gaussian_pos = peak_1_gaussian_pos;
peak_gaussian_std = peak_1_gaussian_std;
peak_gaussian_ampl = peak_1_gaussian_ampl;
else
[peak_gaussian_pos, ind] = sort(peak_1_gaussian_pos, 'ascend');
for i = 1:peaks_nb
peak_gaussian_std(i) = peak_1_gaussian_std(ind(i));
peak_gaussian_ampl(i) = peak_1_gaussian_ampl(ind(i));
end
end
for i = 1:peaks_nb
eval(['thresh_low_', num2str(i), ' = peak_gaussian_pos(', num2str(i), ') - thresh_range*peak_gaussian_std(', num2str(i), ');']);
eval(['thresh_high_', num2str(i), ' = peak_gaussian_pos(', num2str(i), ') + thresh_range*peak_gaussian_std(', num2str(i), ');']);
end
if peaks_nb == 1
else
for i = 1:peaks_nb-1
fct = @(x)(peak_gaussian_ampl(i)*exp(-(x-peak_gaussian_pos(i)).^2./(2*(peak_gaussian_std(i))^2))- peak_gaussian_ampl(i+1)*exp(-(x-peak_gaussian_pos(i+1)).^2./(2*(peak_gaussian_std(i+1))^2)));
eval(['thresh_start_x = (thresh_low_', num2str(i), '+ thresh_high_', num2str(i+1), ')/2;']);
[intersect_x, fval, exit_flag] = fzero(fct, thresh_start_x);
if intersect_x > peak_gaussian_pos(i)+ 5*peak_gaussian_std(i) || intersect_x < peak_gaussian_pos(i) - 5*peak_gaussian_std(i)
exit_flag = 0;
end
if exit_flag == 1
if i == 1
if thresh_high_1 > intersect_x
thresh_high_1 = intersect_x;
end
if thresh_low_2 < intersect_x
thresh_low_2 = intersect_x;
end
end
if i == 2
if thresh_high_2 > intersect_x
thresh_high_2 = intersect_x;
end
if thresh_low_3 < intersect_x
thresh_low_3 = intersect_x;
end
end
end
end
end
[peak_gaussian_ampl, ind] = sort(peak_gaussian_ampl, 'ascend');
for i = 1:peaks_nb
eval(['thresh_low(', num2str(i), ') = thresh_low_', num2str(ind(i)), ';']);
eval(['thresh_high(', num2str(i), ') = thresh_high_', num2str(ind(i)), ';']);
end
if peaks_nb >= peak_to_flatten(j)
eval(['thresh_low = thresh_low_', num2str(peak_to_flatten(j)), ';']);
eval(['thresh_high = thresh_high_', num2str(peak_to_flatten(j)), ';']);
else
thresh_low = thresh_low_1;
thresh_high = thresh_high_1;
end
image_test = images(:,:,j);
image_test((images(:,:,j) < thresh_low)) = NaN;
image_test((images(:,:,j) > thresh_high)) = NaN;
First_Thresh_File = [ImageDir,'Matlab',filesep,'first_threshold.txt'];
if index == 1 && j == 1
fid = fopen(First_Thresh_File, 'a');
fprintf(fid, [ 'Image','\t',...
'Channel','\t',...
'Thresh Low','\t',...
'Thresh High','\n'],...
'char');
fprintf(fid, [ num2str(file_number,'%0.0f'),'\t',...
channel_info(j).Name,'\t',...
num2str(thresh_low,'%0.0f'),'\t',...
num2str(thresh_high,'%0.3f'),'\n'],...
'char');
fclose(fid);
else
fid = fopen(First_Thresh_File, 'a');
fprintf(fid, [ num2str(file_number,'%0.0f'),'\t',...
channel_info(j).Name,'\t',...
num2str(thresh_low,'%0.0f'),'\t',...
num2str(thresh_high,'%0.3f'),'\n'],...
'char');
fclose(fid);
end
coef_length = nchoosek(max([poly_order_horiz(j),poly_order_vert(j)])+2,2);
[xInput_poly, yInput_poly, zOutput_poly] = prepareSurfaceData(1:size(images(:,:,j),1),1:size(images(:,:,j),2) , image_test);
poly_type = fittype( ['poly',num2str(poly_order_vert(j)),num2str(poly_order_horiz(j))] );
poly_opts = fitoptions( poly_type );
if index == 1
poly_opts.Lower = -Inf*ones(coef_length,1);
poly_opts.Upper = Inf*ones(coef_length,1);
else
poly_opt.Lower = coeffvalues(Second_poly_coef{index-1,j})-0.5*abs(coeffvalues(Second_poly_coef{index-1,j}));
poly_opt.Upper = coeffvalues(Second_poly_coef{index-1,j})+0.5*abs(coeffvalues(Second_poly_coef{index-1,j}));
end
poly_result = fit( [xInput_poly, yInput_poly], zOutput_poly, poly_type, poly_opts );
images(:,:,j) = images(:,:,j) - reshape(poly_result([xInput,yInput]),size(images(:,:,j),1),size(images(:,:,j),2));
Second_poly_coef{index,j} = poly_result;
images(:,:,j) = images(:,:,j)-median(median(images(:,:,j)));
if S.GUI_flag
if S.processAll == 0
poly_flat_1_images(:,:,channels_indices(j)) = images(:,:,j);
end
end
if plot_intermediate_steps == 1
figure(summary_fig)
subplot(2, 7, 4);
plot_nice_image;
title('After 1st polynomial flattening')
subplot(2, 7, 11);
image_inside_range = images(:,:,j);
image_inside_range((images(:,:,j) < min_hist)) = NaN;
image_inside_range((images(:,:,j) > max_hist)) = NaN;
hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);
plot_peaks;
title('Histogram')
xlabel('This step improved the image');
end
disp('--- First Poly Flattening ----')
[improv_param, peak_2_gaussian_std, peak_2_gaussian_ampl] = check_improvement(images(:,:,j), ...
peak_1_gaussian_std, peak_1_gaussian_ampl,...
min_spacing_stdev(j), max_spacing_stdev(j), ...
amplitude_fract, min_peak_distance, min_peak_height, peak_threshold, ...
position_est, peak_height_est, peak_index_est, ...
peak_pos_tolerance, histogram_bins, smoothing_level,...
min_hist, max_hist ,ref_peak_pos, peaks_fit(j));
if max(peak_gaussian_ampl) > max_peak_ampl
max_peak_ampl = max(peak_gaussian_ampl);
end
if improv_param < 0
images(:,:,j) = saved_image;
peak_2_gaussian_std = peak_1_gaussian_std;
peak_2_gaussian_ampl = peak_1_gaussian_ampl;
disp('The first polymomial flattening did NOT improve the image')
if plot_intermediate_steps == 1
xlabel('This step did NOT improve the image');
end
else
end
if improv_param == 0
if plot_intermediate_steps == 1
xlabel('No change')
end
end
[height position] = hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins);
image_inside_range = images(:,:,j);
image_inside_range((images(:,:,j) < min_hist)) = NaN;
image_inside_range((images(:,:,j) > max_hist)) = NaN;
[height position] = hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);
height = smooth(height,smoothing_level,'loess');
[peak_height peak_index] = findpeaks(height,'MinpeakDistance', min_peak_distance, 'MinpeakHeight', min_peak_height, 'Threshold', peak_threshold, 'Sortstr','descend');
peak_index = find_closest(ref_peak_pos, peak_index(1:min(7, length(peak_index))), position);
ref_peak_pos = position(peak_index);
After the first round of thresholded plane correction this does another median correction, respecting the above thresholds to correct scaring
if Med_corr_2(j) == 1
saved_image(:,:) = images(:,:,j);
disp('--- Second median correction ----')
[images(:,:,j), peak_gaussian_std, peak_gaussian_ampl, improv_param] = medianCorrection2(images(:,:,j), 2.5, ...
peak_2_gaussian_std, peak_2_gaussian_ampl, ...
min_spacing_stdev(j), max_spacing_stdev(j), ...
amplitude_fract, min_peak_distance, min_peak_height, peak_threshold, ...
position_est, peak_height_est, peak_index_est, ...
peak_pos_tolerance, histogram_bins, smoothing_level, ...
min_hist, max_hist, ref_peak_pos, peaks_fit(j));
images(:,:,j) = images(:,:,j)-median(median(images(:,:,j)));
if S.GUI_flag
if S.processAll == 0
med_corr_2_images(:,:,channels_indices(j)) = images(:,:,j);
end
end
if plot_intermediate_steps == 1
figure(summary_fig)
subplot(2, 7, 5);
plot_nice_image;
title('After 2nd median correction')
subplot(2, 7, 12);
image_inside_range = images(:,:,j);
image_inside_range((images(:,:,j) < min_hist)) = NaN;
image_inside_range((images(:,:,j) > max_hist)) = NaN;
hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);
plot_peaks;
title('Histogram')
end
if max(peak_gaussian_ampl) > max_peak_ampl
max_peak_ampl = max(peak_gaussian_ampl);
end
if improv_param < 0
images(:,:,j) = saved_image(:,:);
peak_gaussian_std = peak_2_gaussian_std;
peak_gaussian_ampl = peak_2_gaussian_ampl;
if plot_intermediate_steps == 1
xlabel('This step did NOT improve the image');
end
elseif improv_param == 0
if plot_intermediate_steps == 1
xlabel('No change')
end
else
if plot_intermediate_steps == 1
xlabel('This step improved the image');
end
end
end
[height position] = hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),histogram_bins);
This does a second round of thresholding and polynomial flattening to improve on the results from before.
saved_image = images(:,:,j);
image_inside_range = images(:,:,j);
image_inside_range((images(:,:,j) < min_hist)) = NaN;
image_inside_range((images(:,:,j) > max_hist)) = NaN;
[height_final position_final] = hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);
height_final = smooth(height_final,smoothing_level,'loess');
[peak_height_final peak_index_final] = findpeaks(height_final,'MinpeakDistance', min_peak_distance, 'MinpeakHeight', min_peak_height, 'Threshold', peak_threshold, 'Sortstr','descend');
peak_index_final = find_closest(ref_peak_pos, peak_index_final(1:min(peaks_fit(j)+1, length(peak_index_final))), position_final);
ref_peak_pos = position_final(peak_index_final);
if sort_peaks(j) == 1
peak_index_final = sort(peak_index_final(1:min(peaks_fit(j), length(peak_index_final))),'ascend');
else
peak_index_final = peak_index_final(1:min(peaks_fit(j), length(peak_index_final)));
end
peak_height_final = height_final(peak_index_final);
func_form = [];
t_hist_lower_a = zeros(1,min(3,length(peak_height_final)));
t_hist_lower_s = zeros(1,min(3,length(peak_height_final)));
t_hist_lower_x = zeros(1,min(3,length(peak_height_final)));
t_hist_upper_a = zeros(1,min(3,length(peak_height_final)));
t_hist_upper_s = zeros(1,min(3,length(peak_height_final)));
t_hist_upper_x = zeros(1,min(3,length(peak_height_final)));
t_hist_start_a = zeros(1,min(3,length(peak_height_final)));
t_hist_start_s = zeros(1,min(3,length(peak_height_final)));
t_hist_start_x = zeros(1,min(3,length(peak_height_final)));
for i = 1:min(3,length(peak_height_final))
t_hist_lower_a(i) = (1-amplitude_fract)*peak_height_final(i);
t_hist_lower_s(i) = min_spacing_stdev(j);
t_hist_lower_x(i) = position_final(peak_index_final(i))-3.5/histogram_bins;
t_hist_upper_a(i) = (1+amplitude_fract)*peak_height_final(i);
t_hist_upper_s(i) = max_spacing_stdev(j);
t_hist_upper_x(i) = position_final(peak_index_final(i))+3.5/histogram_bins;
t_hist_start_a(i) = 1.00*peak_height_final(i);
t_hist_start_s(i) = (min_spacing_stdev(j)+max_spacing_stdev(j))/2;
t_hist_start_x(i) = position_final(peak_index_final(i));
func_form = [func_form,'a',num2str(i),'*exp(-(x-x',num2str(i),')^2/(2*(s',num2str(i),')^2))'];
if i < min(3,length(peak_height_final))
func_form = [func_form,'+'];
end
end
clear peak_gaussian_std peak_gaussian_ampl peak_gaussian_pos;
t_hist_lower = [t_hist_lower_a,t_hist_lower_s,t_hist_lower_x];
t_hist_upper = [t_hist_upper_a,t_hist_upper_s,t_hist_upper_x];
t_hist_start = [t_hist_start_a,t_hist_start_s,t_hist_start_x];
t_hist = fitoptions('Method','NonlinearLeastSquares',...
'Lower', t_hist_lower,...
'Upper', t_hist_upper,...
'Startpoint', t_hist_start,...
'MaxIter', 1e2,...
'MaxFunEvals',3e2);
f_Number_hist = fittype(func_form, 'options' ,t_hist);
[peak_final_gaussian,peak_final_gaussian_gof] = fit( position_final',...
height_final,...
f_Number_hist);
peaks_nb = min(3,length(peak_height_final));
for i = 1:min(3,length(peak_height_final))
eval(['peak_ampl = peak_final_gaussian.a', num2str(i), ';']);
if peak_ampl < min_peak_height
peaks_nb = peaks_nb - 1;
end
end
peak_final_gaussian_ampl = [];
peak_final_gaussian_pos = [];
peak_final_gaussian_std = [];
for i = 1:peaks_nb
eval(['peak_final_gaussian_ampl(', num2str(i), ')= peak_final_gaussian.a', num2str(i), ';']);
eval(['peak_final_gaussian_pos(', num2str(i), ')= peak_final_gaussian.x', num2str(i), ';']);
eval(['peak_final_gaussian_std(', num2str(i), ')= peak_final_gaussian.s', num2str(i), ';']);
end
if peaks_nb == 1
peak_gaussian_pos = peak_final_gaussian_pos;
peak_gaussian_std = peak_final_gaussian_std;
peak_gaussian_ampl = peak_final_gaussian_ampl;
else
[peak_gaussian_pos, ind] = sort(peak_final_gaussian_pos, 'ascend');
for i = 1:peaks_nb
peak_gaussian_std(i) = peak_final_gaussian_std(ind(i));
peak_gaussian_ampl(i) = peak_final_gaussian_ampl(ind(i));
end
end
for i = 1:peaks_nb
eval(['thresh_low_', num2str(i), ' = peak_gaussian_pos(', num2str(i), ') - thresh_range*peak_gaussian_std(', num2str(i), ');']);
eval(['thresh_high_', num2str(i), ' = peak_gaussian_pos(', num2str(i), ') + thresh_range*peak_gaussian_std(', num2str(i), ');']);
end
if peaks_nb == 1
else
for i = 1:peaks_nb-1
fct = @(x)(peak_gaussian_ampl(i)*exp(-(x-peak_gaussian_pos(i)).^2./(2*(peak_gaussian_std(i))^2))- peak_gaussian_ampl(i+1)*exp(-(x-peak_gaussian_pos(i+1)).^2./(2*(peak_gaussian_std(i+1))^2)));
eval(['thresh_start_x = (thresh_low_', num2str(i), '+ thresh_high_', num2str(i+1), ')/2;']);
[intersect_x, fval, exit_flag] = fzero(fct, thresh_start_x);
if intersect_x > peak_gaussian_pos(i)+ 5*peak_gaussian_std(i) || intersect_x < peak_gaussian_pos(i) - 5*peak_gaussian_std(i)
exit_flag = 0;
end
if exit_flag == 1
if i == 1
if thresh_high_1 > intersect_x
thresh_high_1 = intersect_x;
end
if thresh_low_2 < intersect_x
thresh_low_2 = intersect_x;
end
end
if i == 2
if thresh_high_2 > intersect_x
thresh_high_2 = intersect_x;
end
if thresh_low_3 < intersect_x
thresh_low_3 = intersect_x;
end
end
end
end
end
[peak_gaussian_ampl, ind] = sort(peak_gaussian_ampl, 'ascend');
for i = 1:peaks_nb
eval(['thresh_low(', num2str(i), ') = thresh_low_', num2str(ind(i)), ';']);
eval(['thresh_high(', num2str(i), ') = thresh_high_', num2str(ind(i)), ';']);
end
if peaks_nb >= peak_to_flatten(j)
eval(['thresh_low = thresh_low_', num2str(peak_to_flatten(j)), ';']);
eval(['thresh_high = thresh_high_', num2str(peak_to_flatten(j)), ';']);
else
thresh_low = thresh_low_1;
thresh_high = thresh_high_1;
end
image_test = images(:,:,j);
image_test((images(:,:,j) < thresh_low)) = NaN;
image_test((images(:,:,j) > thresh_high)) = NaN;
Second_Thresh_File = [ImageDir,'Matlab',filesep,'second_threshold.txt'];
if index == 1 && j == 1
fid = fopen(Second_Thresh_File, 'a');
fprintf(fid, [ 'Image','\t',...
'Channel','\t',...
'Thresh Low','\t',...
'Thresh High','\n'],...
'char');
fprintf(fid, [ num2str(file_number,'%0.0f'),'\t',...
channel_info(j).Name,'\t',...
num2str(thresh_low,'%0.0f'),'\t',...
num2str(thresh_high,'%0.3f'),'\n'],...
'char');
fclose(fid);
else
fid = fopen(Second_Thresh_File, 'a');
fprintf(fid, [ num2str(file_number,'%0.0f'),'\t',...
channel_info(j).Name,'\t',...
num2str(thresh_low,'%0.0f'),'\t',...
num2str(thresh_high,'%0.3f'),'\n'],...
'char');
fclose(fid);
end
coef_length = nchoosek(max([poly_order_horiz(j),poly_order_vert(j)])+2,2);
[xInput_poly, yInput_poly, zOutput_poly] = prepareSurfaceData(1:size(images(:,:,j),1),1:size(images(:,:,j),2) , image_test);
poly_type = fittype( ['poly',num2str(poly_order_vert(j)),num2str(poly_order_horiz(j))] );
poly_opts = fitoptions( poly_type );
if index == 1
poly_opts.Lower = -Inf*ones(coef_length,1);
poly_opts.Upper = Inf*ones(coef_length,1);
else
poly_opt.Lower = coeffvalues(Third_poly_coef{index-1,j})-0.5*abs(coeffvalues(Third_poly_coef{index-1,j}));
poly_opt.Upper = coeffvalues(Third_poly_coef{index-1,j})+0.5*abs(coeffvalues(Third_poly_coef{index-1,j}));
end
poly_result = fit( [xInput_poly, yInput_poly], zOutput_poly, poly_type, poly_opts );
images(:,:,j) = images(:,:,j) - reshape(poly_result([xInput,yInput]),size(images(:,:,j),1),size(images(:,:,j),2));
Third_poly_coef{index,j} = poly_result;
images(:,:,j) = images(:,:,j) - median(median(images(:,:,j)));
if S.GUI_flag
if S.processAll == 0
poly_flat_2_images(:,:,channels_indices(j)) = images(:,:,j);
end
end
if plot_intermediate_steps == 1
figure(summary_fig)
subplot(2, 7, 6);
plot_nice_image;
title('After 2nd polynomial flattening')
subplot(2, 7, 13);
image_inside_range = images(:,:,j);
image_inside_range((images(:,:,j) < min_hist)) = NaN;
image_inside_range((images(:,:,j) > max_hist)) = NaN;
hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);
plot_peaks;
title('Histogram')
xlabel('This step improved the image');
end
disp('--- Second Poly Flattening ----')
[improv_param, peak_gaussian_std, peak_gaussian_ampl] = check_improvement(images(:,:,j), ...
peak_final_gaussian_std, peak_final_gaussian_ampl,...
min_spacing_stdev(j), max_spacing_stdev(j), ...
amplitude_fract, min_peak_distance, min_peak_height, peak_threshold, ...
position_est, peak_height_est, peak_index_est, ...
peak_pos_tolerance, histogram_bins, smoothing_level, ...
min_hist, max_hist, ref_peak_pos, peaks_fit(j));
if max(peak_gaussian_ampl) > max_peak_ampl
max_peak_ampl = max(peak_gaussian_ampl);
end
if improv_param < 0
images(:,:,j) = saved_image;
peak_gaussian_std = peak_final_gaussian_std;
peak_gaussian_ampl = peak_final_gaussian_ampl;
disp('The second polymomial flattening did NOT improve the image')
if plot_intermediate_steps == 1
figure(summary_fig)
xlabel('This step did NOT improve the image')
end
else
end
if improv_param == 0 && plot_intermediate_steps == 1
figure(summary_fig)
xlabel('No change')
end
Computes the final mask for direct one-step correction of the initial image (image after first median correction)
if sequential_flattening == 1
else
image_test = final_flattening_image(:,:);
image_test((images(:,:,j) < thresh_low)) = NaN;
image_test((images(:,:,j) > thresh_high)) = NaN;
coef_length = nchoosek(max([poly_order_horiz(j),poly_order_vert(j)])+2,2);
[xInput_poly, yInput_poly, zOutput_poly] = prepareSurfaceData(1:size(images(:,:,j),1),1:size(images(:,:,j),2) , image_test);
poly_type = fittype( ['poly',num2str(poly_order_vert(j)),num2str(poly_order_horiz(j))] );
poly_opts = fitoptions( poly_type );
if index == 1
poly_opts.Lower = -Inf*ones(coef_length,1);
poly_opts.Upper = Inf*ones(coef_length,1);
else
poly_opt.Lower = coeffvalues(Third_poly_coef{index-1,j})-0.5*abs(coeffvalues(Third_poly_coef{index-1,j}));
poly_opt.Upper = coeffvalues(Third_poly_coef{index-1,j})+0.5*abs(coeffvalues(Third_poly_coef{index-1,j}));
end
poly_result = fit( [xInput_poly, yInput_poly], zOutput_poly, poly_type, poly_opts );
images(:,:,j) = final_flattening_image(:,:) - reshape(poly_result([xInput,yInput]),size(images(:,:,j),1),size(images(:,:,j),2));
images(:,:,j) = images(:,:,j) - median(median(images(:,:,j)));
end
Applies an offset to the whole image to put the reference peak at a given height reference_height
if Offset_data(j) == 1
image_inside_range = images(:,:,j);
image_inside_range((images(:,:,j) < min_hist)) = NaN;
image_inside_range((images(:,:,j) > max_hist)) = NaN;
[height_offset_calc position_offset_calc] = hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);
[peak_height_hist peak_index_hist] = findpeaks(smooth(height_offset_calc,smoothing_level,'loess'),'MinpeakDistance', min_peak_distance, 'MinpeakHeight', min_peak_height, 'Threshold', peak_threshold, 'Sortstr','descend');
peak_index_hist = find_closest(ref_peak_pos, peak_index_hist(1:min(peaks_fit(j)+1, length(peak_index_hist))), position_offset_calc);
peak_index_hist = sort(peak_index_hist(1:min([3, peaks_fit(j), length(peak_index_hist)])),'ascend');
peak_height_hist = height_offset_calc(peak_index_hist);
func_form = [];
t_hist_lower_a = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_lower_s = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_lower_x = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_upper_a = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_upper_s = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_upper_x = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_start_a = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_start_s = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_start_x = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
for i = 1:min(peaks_fit(j),length(peak_height_hist))
t_hist_lower_a(i) = (1-amplitude_fract)*peak_height_hist(i);
t_hist_lower_s(i) = min_spacing_stdev(j);
t_hist_lower_x(i) = position_offset_calc(max([1,peak_index_hist(i)-4]));
t_hist_upper_a(i) = (1+amplitude_fract)*peak_height_hist(i);
t_hist_upper_s(i) = max_spacing_stdev(j);
t_hist_upper_x(i) = position_offset_calc(min([peak_index_hist(i)+4,length(position_offset_calc)]));
t_hist_start_a(i) = 1.00*peak_height_hist(i);
t_hist_start_s(i) = (min_spacing_stdev(j)+max_spacing_stdev(j))/2;
t_hist_start_x(i) = position_offset_calc(peak_index_hist(i));
func_form = [func_form,'a',num2str(i),'*exp(-(x-x',num2str(i),')^2/(2*(s',num2str(i),')^2))'];
if i < min(peaks_fit(j),length(peak_height_hist))
func_form = [func_form,'+'];
end
end
t_hist_lower = [t_hist_lower_a,t_hist_lower_s,t_hist_lower_x];
t_hist_upper = [t_hist_upper_a,t_hist_upper_s,t_hist_upper_x];
t_hist_start = [t_hist_start_a,t_hist_start_s,t_hist_start_x];
t_hist = fitoptions( 'Method','NonlinearLeastSquares',...
'Lower', t_hist_lower,...
'Upper', t_hist_upper,...
'Startpoint', t_hist_start,...
'MaxIter', 1e2,...
'MaxFunEvals',3e2);
f_Number_hist = fittype(func_form, 'options' ,t_hist);
[peak_hist_gaussian,peak_hist_gaussian_gof] = fit( position_offset_calc',...
smooth(height_offset_calc,smoothing_level,'loess'),...
f_Number_hist);
peak_hist_gaussian_ampl = [];
peak_hist_gaussian_pos = [];
peak_hist_gaussian_std = [];
for i = 1:length(peak_height_hist)
eval(['peak_hist_gaussian_ampl(', num2str(i), ')= peak_hist_gaussian.a', num2str(i), ';']);
eval(['peak_hist_gaussian_pos(', num2str(i), ')= peak_hist_gaussian.x', num2str(i), ';']);
eval(['peak_hist_gaussian_std(', num2str(i), ')= peak_hist_gaussian.s', num2str(i), ';']);
end
[peak_hist_pos, ind] = sort(peak_hist_gaussian_pos, 'ascend');
for i = 1:length(peak_hist_pos)
peak_hist_ampl(i) = peak_hist_gaussian_ampl(ind(i));
end
if length(peak_hist_pos) >= peak_reference_height(j)
eval(['peak_hist_gaussian_x_ref = peak_hist_pos(', num2str(peak_reference_height(j)), ');']);
elseif peak_reference_height(j) > 1 && length(peak_height_hist) >= peak_reference_height(j)-1
eval(['peak_hist_gaussian_x_ref = peak_hist_pos(', num2str(peak_reference_height(j)-1), ');']);
else
peak_hist_gaussian_x_ref = peak_hist_pos(1);
end
images(:,:,j) = images(:,:,j) + z_height_reference(j) - peak_hist_gaussian_x_ref;
else
end
image_inside_range = images(:,:,j);
image_inside_range((images(:,:,j) < min_hist)) = NaN;
image_inside_range((images(:,:,j) > max_hist)) = NaN;
[height_final position_final] = hist(reshape(image_inside_range(:,:),size(image_inside_range,1)*size(image_inside_range,2),1),histogram_bins);
[peak_height_final peak_index_final] = findpeaks(smooth(height_final,smoothing_level,'loess'),'MinpeakDistance', min_peak_distance, 'MinpeakHeight', min_peak_height, 'Threshold', peak_threshold, 'Sortstr','descend');
peak_index_final = sort(peak_index_final(1:min([3, peaks_fit(j), length(peak_index_final)])),'ascend');
peak_height_final = height_final(peak_index_final);
peak_position_final = position_final(peak_index_final);
while length(peak_height_final) < max(peaks_fit)
peak_height_final = [peak_height_final NaN];
peak_position_final = [peak_position_final NaN];
end
if plot_intermediate_steps == 1
figure(summary_fig)
subplot(2, 7, 10)
axis([min_hist max_hist 0 max_peak_ampl*1.5])
subplot(2, 7, 11)
axis([min_hist max_hist 0 max_peak_ampl*1.5])
subplot(2, 7, 12)
axis([min_hist max_hist 0 max_peak_ampl*1.5])
subplot(2, 7, 13)
axis([min_hist max_hist 0 max_peak_ampl*1.5])
end
all_images(index, :, :, j) = images(:,:,j);
Computes a final histogram and CDF of the corrected image and fits Gaussians (like before) Displays the final histogram and fits and saves plots
if S.GUI_flag == 0
if index == 1 && j == channels_flatten(1)
figure_hand_hist = figure('Name','Hist');
else
figure(figure_hand_hist);
end
else
axes(S.ax2)
end
[tok remain] = strtok(Filename{index},'.');
file_number = str2double(remain(2:4));
if index == 1 && ~overwrite_xlim(j) == 1
[h_lim x_lim] = ecdf(reshape(images(:,:,j),size(images,1)*size(images,2),1));
min_z(j) = floor(x_lim(find(h_lim >= 0.01,1)));
max_z(j) = ceil(x_lim(find(h_lim >= 0.99,1)));
span_z = max_z(j) - min_z(j);
min_z(j) = ceil(min_z(j)-0.05*span_z);
max_z(j) = floor(max_z(j)+0.05*span_z);
else
end
[y1_data x_data] = hist(reshape(images(:,:,j),size(images,1)*size(images,2),1),min_z(j):(max_z(j)-min_z(j))/(histogram_bins-1):max_z(j));
y2_data = zeros(size(y1_data,1),1);
x_data = x_data';
y1_data = y1_data';
y1_data = 100*y1_data/sum(y1_data);
for i = 1:size(y1_data,1)
y2_data(i,1) = sum(y1_data(1:i));
end
[AX,H1,H2] = plotyy(x_data,y1_data,x_data,y2_data,'plot');
set(AX(1),'YColor','k');
set(AX(2),'YColor','r');
axis auto;
set(AX(1),'YLim',[-0.002*ceil(max(y1_data)) 1.002*ceil(max(y1_data))],'YTick',0:ceil(max(y1_data))/10:ceil(max(y1_data)),'YTickLabel',0:ceil(max(y1_data))/10:ceil(max(y1_data)))
set(AX(2),'YLim',[-0.2 100.2],'YTick',0:10:100,'YTickLabel',0:10:100)
axis_lim = [min_z(j),max_z(j)];
set(AX(1),'XLim',axis_lim(1:2),'XTick',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2),'XTickLabel',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2));
set(AX(2),'XLim',axis_lim(1:2),'XTick',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2),'XTickLabel',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2));
set(AX(1),'yticklabel',sprintf('%0.1f |',get(AX(1),'ytick')'));
set(AX(2),'yticklabel',sprintf('%0.0f |',get(AX(2),'ytick')'));
set(get(AX(1),'Ylabel'),'String','% Area','FontSize', FontSize, 'FontName', FontName,'Color','k');
set(get(AX(2),'Ylabel'),'String','% Below Given Value','FontSize', FontSize, 'FontName', FontName,'Color','r');
set(H1,'Color','k','LineWidth',2);
set(H2,'Color','r','LineWidth',2);
xlabel([channel_info(j).Name,' (',channel_info(j).Unit,')'],'FontSize', FontSize, 'FontName', FontName );
set(AX(1),'FontSize', FontSize, 'FontName', FontName);
set(AX(2),'FontSize', FontSize, 'FontName', FontName);
if disp_percent == 1
percent_val = zeros(length(percent_search),1);
text_disp = cell(length(percent_search),1);
for i = 1:length(percent_search)
index_height = find(y2_data >= percent_search(i),1);
percent_val(i) = x_data(index_height);
text_disp{i,1} = [num2str(percent_search(i)),'% below ',num2str(percent_val(i),'%02.2f'),' ',channel_info(j).Unit];
end
else
end
x_range = get(AX(1),'XLim');
y_range = get(AX(1),'YLim');
text(x_range(2) -.02*(x_range(2)-x_range(1)),y_range(2)-.02*(y_range(2)-y_range(1)),text_disp,'FontSize', 0.75*FontSize, 'FontName', FontName,'VerticalAlignment','top','HorizontalAlignment','right' );
if fit_hist(j) == 1
[peak_height_hist peak_index_hist] = findpeaks(smooth(y1_data,smoothing_level,'loess'),'MinpeakDistance', 8, 'Sortstr','descend');
func_form = [];
t_hist_lower_a = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_lower_s = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_lower_x = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_upper_a = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_upper_s = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_upper_x = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_start_a = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_start_s = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_start_x = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
for i = 1:min(peaks_fit(j),length(peak_height_hist))
t_hist_lower_a(i) = (1-amplitude_fract)*peak_height_hist(i);
t_hist_lower_s(i) = min_spacing_stdev(j);
t_hist_lower_x(i) = x_data(max([1,peak_index_hist(i)-3]));
t_hist_upper_a(i) = (1+amplitude_fract)*peak_height_hist(i);
t_hist_upper_s(i) = max_spacing_stdev(j);
t_hist_upper_x(i) = x_data(min([peak_index_hist(i)+3,length(x_data)]));
t_hist_start_a(i) = 1.00*peak_height_hist(i);
t_hist_start_s(i) = (min_spacing_stdev(j)+max_spacing_stdev(j))/2;
t_hist_start_x(i) = x_data(peak_index_hist(i));
func_form = [func_form,'a',num2str(i),'*exp(-(x-x',num2str(i),')^2/(2*(s',num2str(i),')^2))'];
if i < min(peaks_fit(j),length(peak_height_hist))
func_form = [func_form,'+'];
end
end
t_hist_lower = [t_hist_lower_a,t_hist_lower_s,t_hist_lower_x];
t_hist_upper = [t_hist_upper_a,t_hist_upper_s,t_hist_upper_x];
t_hist_start = [t_hist_start_a,t_hist_start_s,t_hist_start_x];
t_hist = fitoptions( 'Method','NonlinearLeastSquares',...
'Lower', t_hist_lower,...
'Upper', t_hist_upper,...
'Startpoint', t_hist_start,...
'MaxIter', 1e2,...
'MaxFunEvals',3e2);
f_Number_hist = fittype(func_form, 'options' ,t_hist);
[peak_hist_gaussian,peak_hist_gaussian_gof] = fit( x_data,...
smooth(y1_data,smoothing_level,'loess'),...
f_Number_hist);
hold on
plot(x_data,peak_hist_gaussian(x_data),'.b','LineWidth',2);
for i = 1:min(peaks_fit(j),length(peak_height_hist))
text(1.05*eval(['peak_hist_gaussian.x',num2str(i)]),...
0.85*eval(['peak_hist_gaussian.a',num2str(i)]),...
{ [num2str(sum(eval(['peak_hist_gaussian.a',num2str(i),'*exp(-(x_data-peak_hist_gaussian.x',num2str(i),').^2/(2*(peak_hist_gaussian.s',num2str(i),')^2))'])),'%02.2f'),'% @'];...
[num2str(eval(['peak_hist_gaussian.x',num2str(i)]),'%02.2f'),' +- ',num2str(eval(['peak_hist_gaussian.s',num2str(i)]),'%02.2f'),' ',channel_info(j).Unit]},...
'FontSize', 0.75*FontSize, 'FontName', FontName,'VerticalAlignment','bottom','HorizontalAlignment','left','color','b' );
end
else
end
if S.GUI_flag == 0
print('-loose', '-djpeg' , [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/hist/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_hist_cdf.jpg']);
saveas(gcf, [ImageDir,'Matlab',filesep,'Figure Figs/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_hist_cdf.fig']);
hold off;
else
if S.processAll == 1
end
end
hold off;
Plots the final corrected image with user set color and scale bars and saves the output
ScanZ = max_z(j)-min_z(j);
if S.GUI_flag == 0
if index == 1 && j == channels_flatten(1)
figure_hand_image = figure('Name','Image');
else
figure(figure_hand_image);
end
else
axes(S.ax1)
end
imshow(imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),[min_z(j),max_z(j)])
if S.GUI_flag == 0
print('-dtiff' , [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/grey/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_grey.tiff']);
else
end
colormap(mycmap_sky);
hcb= colorbar('YTick',min_z(j):ScanZ/10:max_z(j),...
'YTickLabel',{ [num2str(min_z(j)+0*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+1*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+2*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+3*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+4*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+5*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+6*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+7*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+8*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+9*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+10*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit]},...
'FontSize', FontSize, 'FontName', FontName);
set(hcb,'YTickMode','manual','FontSize', FontSize, 'FontName', FontName);
DistFactor = 0.05;
Length = ScaleBar;
Pos = [0.5*max(size(images(:,:,j))),0.95*max(size(images(:,:,j)))];
Scale = ScanSize/(size(images,2)*(1+crop_fraction));
UnitsName = ' nm';
NextPlot = get(gca,'NextPlot');
hold on;
XLim = get(gca,'XLim');
YLim = get(gca,'YLim');
Xdiff = abs(diff(XLim));
Ydiff = abs(diff(YLim));
LineX = Pos(1) + 0.5.*Length./Scale.*[-1;+1];
LineY = Pos(2).*[+1;+1];
DistXdir = 0;
DistYdir = -1;
plot(LineX,LineY,'w-','LineWidth',2);
DistX = DistXdir.*DistFactor.*Xdiff;
DistY = DistYdir.*DistFactor.*Ydiff;
Htext = text(Pos(1)+DistX,Pos(2)+DistY,sprintf('%5.1f %s',Length,UnitsName));
set(Htext,'HorizontalAlignment','center','Color', 'w', 'FontSize', FontSize, 'FontName', FontName);
set(gca,'XLim',XLim);
set(gca,'YLim',YLim);
set(gca,'NextPlot',NextPlot);
set(gca, 'XTickMode', 'manual');
if S.GUI_flag == 0
print('-loose', '-djpeg' , [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'.jpg']);
hold off;
else
end
hold off;
if S.GUI_flag
if S.processAll == 0
final_images(:,:,channels_indices(j)) = images(:,:,j);
end
end
pause(0.1)
catch exception
msgbox([{'An error occured while processing image:'}, {file_name}, {'Channel: ', channel_info(j).Name}], 'Error', 'error');
if S.GUI_flag
raw_images(:,:,channels_indices(j)) = zeros(size(images, 1), size(images, 2));
med_corr_1_images(:,:,channels_indices(j)) = zeros(size(images, 1), size(images, 2));
poly_flat_1_images(:,:,channels_indices(j)) = zeros(size(images, 1), size(images, 2));
med_corr_2_images(:,:,channels_indices(j)) = zeros(size(images, 1), size(images, 2));
poly_flat_2_images(:,:,channels_indices(j)) = zeros(size(images, 1), size(images, 2));
final_images(:,:,channels_indices(j)) = zeros(size(images, 1), size(images, 2));
all_images(index, :, :, j) = zeros(size(images, 1), size(images, 2));
if S.processAll == 0
if exist('min_spacing_stdev', 'var') == 0
min_spacing_stdev = min_spacing_stdev_user;
max_spacing_stdev = max_spacing_stdev_user;
end
if exist('min_hist', 'var') == 0
min_hist = 0;
max_hist = 0;
end
if exist('std_est', 'var') == 0
std_est = 0;
end
fit_hist(j) = 0;
pause(2)
end
end
end
end
if S.GUI_flag
if S.processAll == 0
save('Nano6_variables', 'channel_info', 'Filename', ...
'final_images', 'raw_images', 'med_corr_1_images', 'poly_flat_1_images', 'med_corr_2_images', 'poly_flat_2_images', ...
'min_spacing_stdev', 'max_spacing_stdev', 'mycmap_sky', 'min_z', 'max_z', 'histogram_bins', 'ScanSize', ...
'smoothing_level', 'min_hist', 'max_hist', 'std_est');
break;
end
end
end
In case of "Process all", prints all the images and all the histograms in jpg and matlab figures
if S.GUI_flag
if S.processAll == 1
FontSize = 16;
FontName = 'Times';
for index = 1:length(Filename)
for j = 1:length(find(channels_flatten == 1))
if max(max(all_images(index, :, :, j))) == 0 && min(min(all_images(index, :, :, j))) == 0
break;
end
--- Prints all histograms --- %%%
if index == 1 && j == 1
figure_hand_hist = figure('Name','Saving all histograms...');
else
figure(figure_hand_hist);
end
[tok remain] = strtok(Filename{index},'.');
file_number = str2double(remain(2:4));
if index == 1 && ~overwrite_xlim(j) == 1
[h_lim x_lim] = ecdf(reshape(all_images(index, :, :, j),size(all_images,2)*size(all_images,3),1));
min_z(j) = floor(x_lim(find(h_lim >= 0.01,1)));
max_z(j) = ceil(x_lim(find(h_lim >= 0.99,1)));
span_z = max_z(j) - min_z(j);
min_z(j) = ceil(min_z(j)-0.05*span_z);
max_z(j) = floor(max_z(j)+0.05*span_z);
else
end
[y1_data x_data] = hist(reshape(all_images(index, :, :, j),size(all_images,2)*size(all_images,3),1),min_z(j):(max_z(j)-min_z(j))/(histogram_bins-1):max_z(j));
y2_data = zeros(size(y1_data,1),1);
x_data = x_data';
y1_data = y1_data';
y1_data = 100*y1_data/sum(y1_data);
for i = 1:size(y1_data,1)
y2_data(i,1) = sum(y1_data(1:i));
end
[AX,H1,H2] = plotyy(x_data,y1_data,x_data,y2_data,'plot');
set(AX(1),'YColor','k');
set(AX(2),'YColor','r');
axis auto;
set(AX(1),'YLim',[-0.002*ceil(max(y1_data)) 1.002*ceil(max(y1_data))],'YTick',0:ceil(max(y1_data))/10:ceil(max(y1_data)),'YTickLabel',0:ceil(max(y1_data))/10:ceil(max(y1_data)))
set(AX(2),'YLim',[-0.2 100.2],'YTick',0:10:100,'YTickLabel',0:10:100)
axis_lim = [min_z(j),max_z(j)];
set(AX(1),'XLim',axis_lim(1:2),'XTick',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2),'XTickLabel',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2));
set(AX(2),'XLim',axis_lim(1:2),'XTick',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2),'XTickLabel',axis_lim(1):(axis_lim(2)-axis_lim(1))/10:axis_lim(2));
set(AX(1),'yticklabel',sprintf('%0.1f |',get(AX(1),'ytick')'));
set(AX(2),'yticklabel',sprintf('%0.0f |',get(AX(2),'ytick')'));
set(get(AX(1),'Ylabel'),'String','% Area','FontSize', FontSize, 'FontName', FontName,'Color','k');
set(get(AX(2),'Ylabel'),'String','% Below Given Value','FontSize', FontSize, 'FontName', FontName,'Color','r');
set(H1,'Color','k','LineWidth',2);
set(H2,'Color','r','LineWidth',2);
xlabel([channel_info(j).Name,' (',channel_info(j).Unit,')'],'FontSize', FontSize, 'FontName', FontName );
title([channel_info(j).Name,'(',num2str(j),')',' ',num2str(file_number)],'FontSize', FontSize*1.5, 'FontName', FontName );
set(AX(1),'FontSize', FontSize, 'FontName', FontName);
set(AX(2),'FontSize', FontSize, 'FontName', FontName);
if disp_percent == 1
percent_val = zeros(length(percent_search),1);
text_disp = cell(length(percent_search),1);
for i = 1:length(percent_search)
index_height = find(y2_data >= percent_search(i),1);
percent_val(i) = x_data(index_height);
text_disp{i,1} = [num2str(percent_search(i)),'% below ',num2str(percent_val(i),'%02.2f'),' ',channel_info(j).Unit];
end
else
end
x_range = get(AX(1),'XLim');
y_range = get(AX(1),'YLim');
text(x_range(2) -.02*(x_range(2)-x_range(1)),y_range(2)-.02*(y_range(2)-y_range(1)),text_disp,'FontSize', 0.75*FontSize, 'FontName', FontName,'VerticalAlignment','top','HorizontalAlignment','right' );
if fit_hist(j) == 1
[peak_height_hist peak_index_hist] = findpeaks(smooth(y1_data,smoothing_level,'loess'),'MinpeakDistance', 8, 'Sortstr','descend');
func_form = [];
t_hist_lower_a = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_lower_s = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_lower_x = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_upper_a = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_upper_s = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_upper_x = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_start_a = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_start_s = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
t_hist_start_x = zeros(1,min(peaks_fit(j),length(peak_height_hist)));
for i = 1:min(peaks_fit(j),length(peak_height_hist))
t_hist_lower_a(i) = (1-amplitude_fract)*peak_height_hist(i);
t_hist_lower_s(i) = min_spacing_stdev(j);
t_hist_lower_x(i) = x_data(max([1,peak_index_hist(i)-3]));
t_hist_upper_a(i) = (1+amplitude_fract)*peak_height_hist(i);
t_hist_upper_s(i) = max_spacing_stdev(j);
t_hist_upper_x(i) = x_data(min([peak_index_hist(i)+3,length(x_data)]));
t_hist_start_a(i) = 1.00*peak_height_hist(i);
t_hist_start_s(i) = (min_spacing_stdev(j)+max_spacing_stdev(j))/2;
t_hist_start_x(i) = x_data(peak_index_hist(i));
func_form = [func_form,'a',num2str(i),'*exp(-(x-x',num2str(i),')^2/(2*(s',num2str(i),')^2))'];
if i < min(peaks_fit(j),length(peak_height_hist))
func_form = [func_form,'+'];
end
end
t_hist_lower = [t_hist_lower_a,t_hist_lower_s,t_hist_lower_x];
t_hist_upper = [t_hist_upper_a,t_hist_upper_s,t_hist_upper_x];
t_hist_start = [t_hist_start_a,t_hist_start_s,t_hist_start_x];
t_hist = fitoptions( 'Method','NonlinearLeastSquares',...
'Lower', t_hist_lower,...
'Upper', t_hist_upper,...
'Startpoint', t_hist_start,...
'MaxIter', 1e2,...
'MaxFunEvals',3e2);
f_Number_hist = fittype(func_form, 'options' ,t_hist);
[peak_hist_gaussian,peak_hist_gaussian_gof] = fit( x_data,...
smooth(y1_data,smoothing_level,'loess'),...
f_Number_hist);
hold on
plot(x_data,peak_hist_gaussian(x_data),'.b','LineWidth',2);
for i = 1:min(peaks_fit(j),length(peak_height_hist))
text(1.05*eval(['peak_hist_gaussian.x',num2str(i)]),...
0.85*eval(['peak_hist_gaussian.a',num2str(i)]),...
{ [num2str(sum(eval(['peak_hist_gaussian.a',num2str(i),'*exp(-(x_data-peak_hist_gaussian.x',num2str(i),').^2/(2*(peak_hist_gaussian.s',num2str(i),')^2))'])),'%02.2f'),'% @'];...
[num2str(eval(['peak_hist_gaussian.x',num2str(i)]),'%02.2f'),' +- ',num2str(eval(['peak_hist_gaussian.s',num2str(i)]),'%02.2f'),' ',channel_info(j).Unit]},...
'FontSize', 0.75*FontSize, 'FontName', FontName,'VerticalAlignment','bottom','HorizontalAlignment','left','color','b' );
end
end
print(gcf, '-loose', '-djpeg' , [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/hist/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_hist_cdf.jpg']);
saveas(gcf, [ImageDir,'Matlab',filesep,'Figure Fig/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_hist_cdf.fig']);
hold off;
Does Kumar's Grain Analysis if turned on
grain_analysis = get(S.disp_grain_analysis, 'value');
image_test = [];
image_particles = [];
if grain_analysis == 1
image_test = all_images(index, :, :, j);
image_particles(1:size(image_test,2),1:size(image_test,3)) = image_test(1,1:size(image_test,2),1:size(image_test,3));
image_particles = wiener2(image_particles,[5,5]);
image_test_min = (min(min(image_particles)));
image_test_span = (max(max(image_particles))-min(min(image_particles)));
part_thresh = (peak_hist_gaussian.a1+peak_hist_gaussian.s1 - image_test_min)/image_test_span;
image_particles = (image_particles - image_test_min)./image_test_span;
image_particles = im2bw(image_particles, part_thresh);
B = bwconncomp(image_particles);
disp(B);
stats = regionprops(B, 'Area');
idx = find([stats.Area] > 25);
BW2 = ismember(labelmatrix(B), idx);
graindata = regionprops(BW2,'basic');
BW2 = bwconncomp(BW2);
disp(BW2);
grain_areas = [graindata.Area]*(ScanSize/size(image_test,2))^2;
grain_centroids = [graindata.Centroid];
[y1_area_data x_area_data] = hist(grain_areas,[min(grain_areas):(max(grain_areas)-min(grain_areas))/(length(grain_areas)*5):max(grain_areas)]);
y2_area_data = zeros(size(y1_area_data,1),1);
x_area_data = x_area_data';
y1_area_data = y1_area_data';
y1_area_data = 100*y1_area_data/sum(y1_area_data);
for i = 1:size(y1_area_data,1)
y2_area_data(i,1) = sum(y1_area_data(1:i));
end
[AX,H1,H2] = plotyy(x_area_data,y1_area_data,x_area_data,y2_area_data,'plot');
set(AX(1),'YColor','k');
set(AX(2),'YColor','r');
axis auto;
set(AX(1),'YLim',[-0.002*ceil(max(y1_area_data)) 1.002*ceil(max(y1_area_data))],'YTick',0:ceil(max(y1_area_data))/10:ceil(max(y1_area_data)),'YTickLabel',0:ceil(max(y1_area_data))/10:ceil(max(y1_area_data)))
set(AX(2),'YLim',[-0.2 100.2],'YTick',0:10:100,'YTickLabel',0:10:100)
set(AX(1),'yticklabel',sprintf('%0.1f |',get(AX(1),'ytick')'));
set(AX(2),'yticklabel',sprintf('%0.0f |',get(AX(2),'ytick')'));
set(get(AX(1),'Ylabel'),'String','% of total Particles','FontSize', FontSize, 'FontName', FontName,'Color','k');
set(get(AX(2),'Ylabel'),'String','% Below Given Value','FontSize', FontSize, 'FontName', FontName,'Color','r');
set(H1,'Color','k','LineWidth',2);
set(H2,'Color','r','LineWidth',2);
xlabel('Particle area [nm^2]','FontSize', FontSize, 'FontName', FontName );
title( 'Particle area [nm^2]','FontSize', FontSize*1.5, 'FontName', FontName );
set(AX(1),'FontSize', FontSize, 'FontName', FontName);
set(AX(2),'FontSize', FontSize, 'FontName', FontName);
x_range = get(AX(1),'XLim');
y_range = get(AX(1),'YLim');
print(gcf, '-loose', '-djpeg' , [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/hist/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_part_area_hist_cdf.jpg']);
part_vol = zeros(length(grain_areas),1);
part_height = zeros(length(grain_areas),1);
image_test = all_images(index, :, :, j);
image_particles(1:size(image_test,2),1:size(image_test,3)) = image_test(1,1:size(image_test,2),1:size(image_test,3));
image_particles = wiener2(image_particles,[5,5]);
for i = 1:length(grain_areas)
grain = false(size(BW2));
grain(BW2.PixelIdxList{i}) = true;
disp(size(grain));
disp(size(image_particles));
end
disp(part_height);
else
end
end
end
close(figure_hand_hist)
--- Prints all images --- %%%
for index = 1:length(Filename)
for j = 1:length(find(channels_flatten == 1))
[tok remain] = strtok(Filename{index},'.');
file_number = str2double(remain(2:4));
ScanZ = max_z(j)-min_z(j);
if index == 1 && j == 1
figure_hand_image = figure('Name','Saving all images...');
else
figure(figure_hand_image);
end
if max(max(all_images(index, :, :, j))) == 0 && min(min(all_images(index, :, :, j))) == 0
break;
end
images(:,:, j) = all_images(index, :, :, j);
imshow(imresize(images(:,:,j),[max(size(images(:,:,j))),max(size(images(:,:,j)))]),[min_z(j),max_z(j)])
print(gcf, '-dtiff' , [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/grey/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'_grey.tiff']);
color_string = get(S.ColorBarSelect, 'string');
color_value = get(S.ColorBarSelect, 'value');
switch color_string{color_value}
case 'Autumn'
colormap(autumn);
case 'Bone'
colormap(bone);
case 'ColorCube'
colormap(colorcube);
case 'Cool'
colormap(cool);
case 'Copper'
colormap(copper);
case 'Flag'
colormap(flag);
case 'Gray'
colormap(gray);
case 'Hot'
colormap(hot);
case 'HSV'
colormap(hsv);
case 'Jet'
colormap(jet);
case 'Lines'
colormap(lines);
case 'Pink'
colormap(pink);
case 'Prism'
colormap(prism);
case 'Spring'
colormap(spring);
case 'Summer'
colormap(summer);
case 'Winter'
colormap(winter);
case 'Sky'
load('Sky','mycmap_sky')
colormap(mycmap_sky);
end
hcb= colorbar('YTick',min_z(j):ScanZ/10:max_z(j),...
'YTickLabel',{ [num2str(min_z(j)+0*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+1*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+2*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+3*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+4*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+5*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+6*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+7*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+8*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+9*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit],...
[num2str(min_z(j)+10*ScanZ/10,'%0.2f'),' ', channel_info(j).Unit]},...
'FontSize', FontSize, 'FontName', FontName);
set(hcb,'YTickMode','manual','FontSize', FontSize, 'FontName', FontName);
disp_scale_bar = get(S.disp_scale_bar, 'value');
if disp_scale_bar == 1
DistFactor = 0.05;
Length = ScaleBar;
Pos = [0.5*max(size(all_images(index, :, :, j))),0.95*max(size(all_images(index, :, :, j)))];
Scale = ScanSize/(size(all_images,3)*(1+crop_fraction));
UnitsName = ' nm';
NextPlot = get(gca,'NextPlot');
hold on;
XLim = get(gca,'XLim');
YLim = get(gca,'YLim');
Xdiff = abs(diff(XLim));
Ydiff = abs(diff(YLim));
LineX = Pos(1) + 0.5.*Length./Scale.*[-1;+1];
LineY = Pos(2).*[+1;+1];
DistXdir = 0;
DistYdir = -1;
plot(LineX,LineY,'w-','LineWidth',2);
DistX = DistXdir.*DistFactor.*Xdiff;
DistY = DistYdir.*DistFactor.*Ydiff;
Htext = text(Pos(1)+DistX,Pos(2)+DistY,sprintf('%5.1f %s',Length,UnitsName));
set(Htext,'HorizontalAlignment','center','Color', 'w', 'FontSize', FontSize, 'FontName', FontName);
set(gca,'XLim',XLim);
set(gca,'YLim',YLim);
set(gca,'NextPlot',NextPlot);
set(gca, 'XTickMode', 'manual');
else
end
print(gcf, '-loose', '-djpeg' , [ImageDir,'Matlab',filesep,'Figure Jpegs/',channel_info(j).Name,'/',channel_info(j).Name,'(',num2str(j),')','_',num2str(file_number),'.jpg']);
hold off;
end
end
close(figure_hand_image)
end
end