Geometry comparison scripts

2025-02-18 update

This update includes:

1. Assorted scripts used for the manuscript "Comparison of Kikuchi Diffraction Geometries in Scanning Electron Microscope" (https://arxiv.org/abs/2411.13018). They can be found `modules/ded_geometry_comparison` with a more detailed `readme`. [^1]

2. Nomenclature change for dynamic templates:
  - The folder `phases/masterpattern` is renamed `dynamic_templates`. Example scripts are changed accordingly.
  - Dynamic patterns simulated from a few different software packages [^2] are now supported. A few examples (using Si) were provided. You may refer to the example .pha files to build phase files on your own, or check out the short guide under `tutorial_docs`.


[^1]: Reconstructed reference spheres with experimental data are uploaded separately to Zenodo. See: 10.5281/zenodo.14111773
.
[^2]: Espirit DynamicS, EMsoft, and the package by Winkelmann et al. (hereby called "bwkd", outlined in doi:10.1111/jmi.13051).

Author: TianbiZhang

.gitignore

@@ -1,9 +1,7 @@
 /utils/gmail.m
 /testing/*
 /outputs/*
-/phases/cifs/*
-/phases/masterpatterns/*
-/phases/phasefiles/*
+
 
 # Logs and databases #
 ######################
@@ -29,3 +27,15 @@ Thumbs.db
 #Matlab Files
 *.asv
 *.fig
+
+# Phase folders
+/phases/cifs
+/phases/cifs/*
+/phases/dynamic_templates
+/phases/dynamic_templates/*
+/phases/phasefiles
+/phases/phasefiles/*
+
+# Phase folders - deprecated
+/phases/masterpatterns
+/phases/masterpatterns/*

decks/AstroEBSD_PatternMatch_2021.m

@@ -52,7 +52,7 @@
 
 %% Low level setting up stuff - you shouldn't need to change this
 RTM.Phase_Folder = fullfile(InputUser.Astro_loc,'phases'); %location of the AstroEBSD phases super-folder
-RTM.Bin_loc = fullfile(RTM.Phase_Folder,'masterpatterns'); %location of the binary files used for RTM
+RTM.Bin_loc = fullfile(RTM.Phase_Folder,'dynamic_templates'); %location of the binary files used for RTM
 
 [ SettingsXCF, correction, SettingsXCF2 ] = FFT_Filter_settings( RTM.screensize, RTM.LPTsize );
 

decks/PCA/PCA_deck.m

@@ -47,7 +47,7 @@
 
 %These need to be the names of '.pha' files in the '\phases\phasefiles'
 %folder. There need to be corresponding '.cif' and '.BIN' files in 'cifs'
-%and 'masterpatterns' folders. - Only relevant to RTM analyses.
+%and 'dynamic_templates' folders. - Only relevant to RTM analyses.
 InputUser.Phases={'Ni','ZrC','M6C'};
 
 % Run PCA on EBSD, EDX, or both

decks/RTM/RTM_multiphase_deck.m

@@ -62,7 +62,7 @@
 
 %% now run the code
 RTM.Phase_Folder = fullfile(InputUser.Astro_loc,'phases'); %location of the AstroEBSD phases super-folder
-RTM.Bin_loc = fullfile(RTM.Phase_Folder,'masterpatterns'); %location of the binary files used for RTM
+RTM.Bin_loc = fullfile(RTM.Phase_Folder,'dynamic_templates'); %location of the binary files used for RTM
 RTM.parsearch=2;
 
 [MapInfo.MapData,MicroscopeData,PhaseData,MapInfo.EBSPData ]=bReadHDF5( InputUser );

gen/loading/EBSP_BGCor.m

@@ -11,22 +11,22 @@
 % %background correction
 % Settings_Cor.gfilt=1; %use a low pass filter
 % Settings_Cor.gfilt_s=5; %low pass filter sigma
-%
+% 
 % %radius mask
 % Settings_Cor.radius=1; %use a radius mask
 % Settings_Cor.radius_frac=0.9; %fraction of the pattern width to use as the mask
-%
+% 
 % %hold pixel
 % Settings_Cor.hotpixel=1; %hot pixel correction
 % Settings_Cor.hot_thresh=1000; %hot pixel threshold
-%
+% 
 % %resize
 % Settings_Cor.resize=1; %resize correction
 % Settings_Cor.size=150; %image width
-%
+% 
 % Settings_Cor.RealBG=0; %use a real BG
 % Settings_Cor.EBSP_bgnum=30; %number of real pattern to use for BG
-%
+% 
 % Settings_Cor.SquareCrop=1; %crop to a square
 %
 %Settings_Cor
@@ -45,13 +45,12 @@
 EBSP2=EBSP;
 
 %check fields exist & create if needed - this is ordered in the order of
-%operations to aid with debugging & adding new correction routines
+%operations to aid with debugging & adding new correction routines 
 %as needed
 
-if ~isfield(Settings_Cor,'tkd_onaxis')
-    Settings_Cor.tkd_onaxis = [0 0 0 0 0]; %[1/0, xc,yc,ro,ri]
+if ~isfield(Settings_Cor,'NaNClean')
+   Settings_Cor.NaNClean=0;
 end
-
 if ~isfield(Settings_Cor,'MeanCentre')
     Settings_Cor.MeanCentre=0;
 end
@@ -117,63 +116,88 @@
 end
 
 if ~isfield(Settings_Cor,'LineError')
-    Settings_Cor.LineError=0;
+   Settings_Cor.LineError=0; 
 end
 
 if ~isfield(Settings_Cor,'MeanCentre')
-    Settings_Cor.LineError=0;
+   Settings_Cor.LineError=0; 
 end
 
-%% Start the corrections
-if Settings_Cor.RealBG == 1
-    if size(EBSP2) == size(Settings_Cor.EBSP_bg) %this if statement helps if the debug isn't working properly
-        EBSP2=EBSP2./Settings_Cor.EBSP_bg;
-    end
+if ~isfield(Settings_Cor,'rkd')
+    Settings_Cor.rkd=0;
 end
 
-if Settings_Cor.tkd_onaxis(1) == 1 %blank the central beam and crop the edge
-    [xgrid,ygrid]=meshgrid(1:size(EBSP2,2),1:size(EBSP2,1));
-    
-    tkd_set=Settings_Cor.tkd_onaxis(2:5);%[3,3,30,80];
-    
-    xc=floor(mean(xgrid(:)))+tkd_set(1);
-    yc=floor(mean(ygrid(:)))+tkd_set(2);
-    ri=tkd_set(3);
-    ro=tkd_set(4);
-    pgrid=sqrt((xgrid-xc).^2+(ygrid-yc).^2);
-    grid_ok=false(size(xgrid));
-    grid_ok(pgrid > ri & pgrid < ro) = true;
-    pat_filtered = zeros(size(EBSP2));
-    pat_filtered(grid_ok)=EBSP2(grid_ok);
-    %     EBSP2=pat_filtered;
-    
-    x_new=xc+(-ro:ro);
-    y_new=yc+(-ro:ro);
-    EBSP2=pat_filtered(y_new,x_new);
-    mean_tkd=mean(EBSP2(EBSP2~=0));
-    std_tkd=std(EBSP2(EBSP2~=0));
-    
-    EBSP2(EBSP2~=0)=(EBSP2(EBSP2~=0)-mean_tkd)/std_tkd;
-    
+%% Start the corrections
+
+if Settings_Cor.rkd == 1 %RKD
+
+    pat2=EBSP2(2:end-1,2:end-1);
+
+    pat2_q1=pat2(1:253,1:218);
+    pat2_q2=pat2(1:218,225:end);
+    pat2_q3=pat2(260:end,1:253);
+    pat2_q4=pat2(225:end,260:end);
+
+    pat_std=std([pat2_q1(:);pat2_q2(:);pat2_q3(:);pat2_q4(:)]);
+    pat_mean=mean([pat2_q1(:);pat2_q2(:);pat2_q3(:);pat2_q4(:)]);
+    EBSP2=randn(size(pat2))*pat_std+pat_mean;
+    EBSP2(1:253,1:218)=pat2_q1;
+    EBSP2(1:218,225:end)=pat2_q2;
+    EBSP2(260:end,1:253)=pat2_q3;
+    EBSP2(225:end,260:end)=pat2_q4;
+
 end
 
 if Settings_Cor.SquareCrop == 1 %crop the image to a square
     EBSP_size=min(size(EBSP2));
     crop.centre = [size(EBSP2, 2)/2 size(EBSP2, 1)/2];
-    s=1:EBSP_size; s=s-mean(s);
-    sy=s+crop.centre(2)+0.5;
+    s=1:EBSP_size; s=s-nanmean(s);
+
+    sy=s+crop.centre(2)+0.5; 
     sx=s+crop.centre(1)+0.5;
+    
+    %helps deal with issue with when the patterns are not even in size to start
+    sx=floor(sx);
+    sy=floor(sy);
     EBSP2=EBSP2(sy,sx);
 end
 
+if Settings_Cor.NaNClean == 1 %NaN correction
+    [y,x]=find(isnan(EBSP2)==1);
+    
+    [threexthreey,threexthreex]=meshgrid(-1:1,-1:1);
+    
+    if numel(y) > 0
+        for s=1:numel(y)
+            y_v=y(s)+threexthreey;
+            x_v=x(s)+threexthreex;
+            
+            y_v(y_v<0)=[];
+            x_v(x_v<0)=[];
+            
+            y_v(y_v<0)=[];
+            x_v(x_v<0)=[];
+            
+            y_v(y_v>size(EBSP2,1))=[];
+            x_v(x_v>size(EBSP2,2))=[];
+                        
+            sv=sub2ind(size(EBSP2),y_v(:),x_v(:));
+            EBSP2(y,x)=median(EBSP2(sv),'omitnan');
+        end
+    end
+    
+end
+
 if Settings_Cor.SatCor == 1 %saturation correction
     I_noise= imnoise(zeros(size(EBSP2)),'gaussian');
-    I_noise=mean(EBSP2(:))+std(EBSP2(:))*I_noise;
+    I_noise=nanmean(EBSP2(:))+nanstd(EBSP2(:))*I_noise;
     mvals=find(EBSP2 == max(EBSP2(:)));
     EBSP2(mvals)=I_noise(mvals);
 end
 
-
+if Settings_Cor.RealBG == 1
+    EBSP2=EBSP2./Settings_Cor.EBSP_bg;
+end
 
 if Settings_Cor.SplitBG == 1
     EBSPw=size(EBSP2,2);
@@ -186,7 +210,7 @@
     EBSP2b_g = imgaussfilt(EBSP2b,gf);
 
     EBSP2=[EBSP2a./EBSP2a_g EBSP2b./EBSP2b_g];
-    
+   
 end
 
 %cor the pattern for hot pixels
@@ -197,8 +221,8 @@
 
 if Settings_Cor.LineError==1
 
-    meanval=mean(EBSP2(:));
-    patstd=std(EBSP2(:));
+    meanval=nanmean(EBSP2(:));
+    patstd=nanstd(EBSP2(:));
 
     loc1=floor(0.9667*size(EBSP,1));
     loc2=ceil(0.98*size(EBSP,1));
@@ -209,7 +233,7 @@
     end
 end
 
-%fix the mean and std
+%fix the mean and nanstd
 EBSP2=fix_mean(EBSP2);
 
 %resize the image
@@ -259,9 +283,13 @@
     r_thresh=Settings_Cor.radius_frac*4/3*cs(1)/2;
 
     [xgrid,ygrid]=meshgrid(1:cs(2),1:cs(1));
-    r_grid=sqrt((xgrid-size(EBSP2,2)/2).^2+(ygrid-size(EBSP2,1)/2).^2);
-    EBSP2(r_grid>=r_thresh) = 0;
-    EBSP2(r_grid<r_thresh)=  EBSP2(r_grid<r_thresh)-mean(EBSP2(r_grid<r_thresh));
+    xgrid=int32(xgrid);
+    ygrid=int32(ygrid);
+    root_values=(xgrid-size(EBSP2,2)/2).^2+(ygrid-size(EBSP2,1)/2).^2;
+    r_grid=sqrt(double(root_values));
+    %%
+    EBSP2(r_grid>=double(r_thresh)) = 0;
+    EBSP2(r_grid<double(r_thresh))=  EBSP2(r_grid<double(r_thresh))-nanmean(EBSP2(r_grid<double(r_thresh)));
 else
     EBSP2=fix_mean(EBSP2);
 end
@@ -275,12 +303,12 @@
         hpix=find(abs(EBSP_med-EBSP2) > 0.8);
         EBSP2(hpix)=EBSP_med(hpix);
     end
-    %     EBSP2=fix_mean(EBSP2);
+%     EBSP2=fix_mean(EBSP2);
 end
 
 if Settings_Cor.MeanCentre==1
     %taken from fix fixmean function, without the make positive part
-    EBSP2=(EBSP2-mean(EBSP2(:)))/std(EBSP2(:));
+    EBSP2=(EBSP2-nanmean(EBSP2(:)))/nanstd(EBSP2(:));
 end
 
 Settings_Cor_out=Settings_Cor;
@@ -304,7 +332,7 @@
 
 function EBSP2=fix_mean(EBSP)
 %zero mean & fix stdev
-EBSP2=(EBSP-mean(EBSP(:)))/std(EBSP(:));
+EBSP2=(EBSP-nanmean(EBSP(:)))/nanstd(EBSP(:));
 %make positive
 EBSP2=EBSP2-min(EBSP2(:))+1;
 end

gen/loading/H5OINA_Convert.m

@@ -0,0 +1,86 @@
+function [ebsdtemp] = H5OINA_Convert(ebsd_data,header_data,file_dset,phase_names,phase_colors,ebsd_nonpm)
+%H5OINA_CONVERT Convert the output from the HFOINA_Read function into a
+%MTEX EBSD container
+%
+% Ben Britton - Feb 2024
+%
+% Adapted from load_h5oina
+
+
+slice_name=['s' file_dset];
+
+CS{1}='notIndexed';
+num_phases=numel(header_data.(slice_name).phase);
+
+for n=1:num_phases
+    phaseN=1+n;
+    Space_Group=double(header_data.(['s' file_dset]).phase{n}.Space_Group);
+    Lattice_Dimensions=double(header_data.(['s' file_dset]).phase{n}.Lattice_Dimensions);
+    Lattice_Angles=double(header_data.(['s' file_dset]).phase{n}.Lattice_Angles);
+
+    if nargin > 5 && ~isempty(phase_names)
+        Mineral=phase_names{n};
+        CS{phaseN} = crystalSymmetry('SpaceId',Space_Group, ...
+            Lattice_Dimensions,...
+            Lattice_Angles,...
+            'Mineral',Mineral,'Color',phase_colors{n});
+    else
+        Mineral=char(header_data.(['s' file_dset]).phase{1}.Phase_Name);
+
+        CS{phaseN} = crystalSymmetry('SpaceId',Space_Group, ...
+            Lattice_Dimensions,...
+            Lattice_Angles,...
+            'Mineral',Mineral);
+    end
+end
+
+rc = rotation.byEuler(double(header_data.(slice_name).Specimen_Orientation_Euler(:)')*degree); % what definition? Bunge?
+
+% set up EBSD data
+rot = rc*rotation.byEuler(ebsd_data.(slice_name).Euler');
+phase = ebsd_data.(slice_name).Phase;
+opt=struct;
+
+% read some fields
+EBSD_fieldnames=fieldnames(ebsd_data.(slice_name));
+
+num_fields=size(EBSD_fieldnames,1);
+for n = 1: num_fields
+    s=ebsd_data.(slice_name).(EBSD_fieldnames{n});
+    if size(s,2) == 1 && size(s,1) == numel(phase)
+        try
+            opt.(EBSD_fieldnames{n})=double(s);
+        catch
+            try
+                opt.(EBSD_fieldnames{n})=s;
+            catch
+            end
+        end
+
+    end
+end
+
+%not put the Euler angles in the options too
+if isfield(ebsd_data.(slice_name),'Euler')
+    opt.euler1=double(ebsd_data.(slice_name).Euler(1,:)');
+    opt.euler2=double(ebsd_data.(slice_name).Euler(2,:)');
+    opt.euler3=double(ebsd_data.(slice_name).Euler(3,:)');
+end
+
+try
+    opt.x=opt.Beam_Position_X;
+    opt.y=opt.Beam_Position_Y;
+catch
+    try
+        opt.x=ebsd_nonpm.prop.x;
+        opt.y=ebsd_nonpm.prop.y;
+    catch
+        error('Beam Positions Not Loaded - was this a pattern match file?')
+    end
+end
+
+ebsdtemp = EBSD(rot,phase,CS,opt,'unitCell',calcUnitCell([opt.x,opt.y]));
+ebsdtemp.opt.Header = header_data.(slice_name);
+
+end
+