TKD_DED

Adding routines for transmission kikuchi diffraction / direct electron diffraction analysis paper.

Work by Tianbi Zhang and Ben Britton (UBC)

gen/loading/h5_pixet.m

@@ -0,0 +1,113 @@
+function [rawdata,numfiles] = h5_pixet(h5filenames,path1,multframe_option,multframe_num)
+%H5_PIXET read patterns from a h5 pixet output
+% (c) Tianbi Zhang and Ben Britton, 2023
+
+[~,numfiles] = size(h5filenames);
+
+rawdata = zeros(254^2,numfiles);
+
+for i=1:numfiles
+    filename = fullfile(path1, h5filenames{i});
+    if exist(filename,'file') == false
+        error(['The input file ' filename ' cannot be found']);
+    end
+
+    if multframe_option == true
+        rawdata(:,i) = rawhdfmult(filename, multframe_num);
+    else
+        if i~=numfiles
+            rawdata(:,i) = rawhdfint(filename) ./ multframe_num; % due to data formet issue, use a special function here
+        else
+            rawdata(:,i) = rawhdfint_special(filename) ./ multframe_num;
+        end
+    end
+
+end
+
+end
+
+function event_data_vector = rawhdfint(h5filename)
+
+% Read width and height values from the HDF5 file. Both should be 256
+% because of the DED geometry
+width = h5read(h5filename,'/Frame_0/Width');
+height = h5read(h5filename,'/Frame_0/Height');
+
+% Extract pixel-by-pixel intensity from the "Event" and "iTOT" data.
+% Change "Event" to "iToT", "ToT", "ToA" etc. based on data structure.
+% Please check the file under HDF5 viewer or use h5info() first.
+event_data = h5read(h5filename,'/Frame_0//SubFrames/Event/Data');
+event_data_matrix = double(reshape(event_data,width,height));
+
+event_data_matrix = event_data_matrix';
+event_data_matrix = imcrop(event_data_matrix, [2 2 253 253]);
+
+%TZ QUERY
+%what are these pixels?
+pix_cleanxval=[11,26,240,44];
+pix_cleanyval=[60,208,22,47];
+event_data_matrix=pix_clean(event_data_matrix,pix_cleanxval,pix_cleanyval);
+
+event_data_vector = double(reshape(event_data_matrix, [],1));
+
+end
+
+function event_data_vector = rawhdfint_special(h5filename)
+% Read width and height values from the HDF5 file. Both should be 256
+% because of the DED geometry
+width = h5read(h5filename,'/Frame_0/Width');
+height = h5read(h5filename,'/Frame_0/Height');
+
+% Extract pixel-by-pixel intensity from the "Event" and "iTOT" data.
+% Change "Event" to "iToT", "ToT", "ToA" etc. based on data structure.
+% Please check the file under HDF5 viewer or use h5info() first.
+event_data = h5read(h5filename,'/Frame_0/Data');
+event_data_matrix = double(reshape(event_data,width,height));
+
+event_data_matrix = event_data_matrix';
+event_data_matrix = imcrop(event_data_matrix, [2 2 253 253]);
+
+pix_cleanxval=[11,26,240,44];
+pix_cleanyval=[60,208,22,47];
+event_data_matrix=pix_clean(event_data_matrix,pix_cleanxval,pix_cleanyval);
+
+event_data_vector = double(reshape(event_data_matrix, [],1));
+end
+
+function h5matrix=pix_clean(h5matrix,xval,yval)
+%clean specific pixels fro the array to 0
+for n=1:numel(xval)
+    h5matrix(xval(n),yval(n)) = 0;
+end
+
+end
+
+
+function event_data_vector = rawhdfmult(h5filename, framenum)
+width = h5read(h5filename,'/Frame_0/Width');
+height = h5read(h5filename,'/Frame_0/Height');
+
+event_data = double(zeros(65536,1));
+
+for i=0:(framenum-1)
+    frame_event_key = strcat('/Frame_',num2str(i),'/SubFrames/Event/Data');
+    event_data_this_frame = double(h5read(h5filename, frame_event_key));
+
+    event_data = event_data + event_data_this_frame;
+end
+
+event_data_matrix = double(reshape(event_data,width,height));
+
+event_data_matrix = event_data_matrix';
+event_data_matrix = imcrop(event_data_matrix, [2 2 253 253]);
+
+%TZ QUERY
+pix_cleanxval=[11,26,240,44];
+pix_cleanyval=[60,208,22,47];
+event_data_matrix=pix_clean(event_data_matrix,pix_cleanxval,pix_cleanyval);
+
+
+event_data_matrix = event_data_matrix / framenum;
+
+event_data_vector = double(reshape(event_data_matrix, [],1));
+end

modules/ded_tkd/TKD_cor_saturated.m

@@ -0,0 +1,19 @@
+function [fillingdomains, satdomains] = TKD_cor_saturated(rawdata,numfiles,sat_index)
+%TKD_COR_SATURATED Summary of this function goes here
+%   Detailed explanation goes here
+%% Identify saturated domains 
+satdomains = cell(1,numfiles);
+for i=1:numfiles
+    satdomains{i} = find(rawdata(:,i) >= sat_index);
+end
+
+fillingdomains = cell(1,numfiles);
+for i=1:numfiles
+    if i==1
+        fillingdomains{i} = find(rawdata(:,i) < sat_index);
+    else
+        fillingdomains{i} = setdiff(satdomains{i-1}, satdomains{i});
+    end
+end
+end
+

modules/ded_tkd/TKD_flatfield.m

@@ -0,0 +1,17 @@
+function [fused_flat,fusedpattern,normfactor, ref_curve] = TKD_flatfield(rawdata,normalizeddata,angle_matrix,fillingdomains,numfiles)
+%TKD_FLATFIELD fit to the radial function, and then flat field with it
+
+
+[normfactor, ref_curve] = TKD_referencecurvefit(rawdata,normalizeddata,angle_matrix,numfiles);
+
+% Stitch the pattern
+[fusedpattern] = TKD_patternfuse(rawdata,normalizeddata,fillingdomains);
+
+%% Flat field
+npix=sqrt(size(rawdata,1));
+
+normfactor = reshape(normfactor, [npix npix]);
+fused_flat = fusedpattern ./ normfactor;
+
+end
+

modules/ded_tkd/TKD_normalize.m

@@ -0,0 +1,29 @@
+function [normalizeddata] = TKD_normalize(fillingdomains,rawdata,numfiles)
+%TKD_NORMALIZE Summary of this function goes here
+%   Detailed explanation goes here
+
+normalizeddata = 0*rawdata;
+
+% Calculate normalization factors
+eventratio = zeros(length(fillingdomains{1}),numfiles-1);
+for j=2:numfiles
+    reference_data = rawdata(:,1);
+    current_data = rawdata(:,j);
+    eventratio(:,j-1) = reference_data(fillingdomains{1})./ current_data(fillingdomains{1});
+end
+
+eventratio(isinf(eventratio)) = NaN; % account for dead pixels
+
+clear reference_data current_data;
+
+scale_factor = mean(eventratio, 1 ,"omitnan");
+
+for i=1:numfiles
+    if i==1
+        normalizeddata(:,i) = rawdata(:,i);
+    else
+        normalizeddata(:,i)= rawdata(:,i) .* scale_factor(i-1);
+    end
+end
+end
+

modules/ded_tkd/TKD_patternfuse.m

@@ -0,0 +1,20 @@
+function [fusedpattern] = TKD_patternfuse(rawdata,normalizeddata,fillingdomains)
+%TKD_PATTERNFUSE Summary of this function goes here
+%   Detailed explanation goes here
+numfiles=size(rawdata,2);
+
+fusedpattern = zeros(254^2,1);
+
+for i=1:numfiles
+    current_data = normalizeddata(:,i);
+    fusedpattern(fillingdomains{i}) = current_data(fillingdomains{i});
+end
+
+fusedpattern = reshape(fusedpattern, 254,254);
+stitched_mask = medfilt2(fusedpattern, [3 3]);
+fusedpattern(rawdata(:,1)==0) = stitched_mask(rawdata(:,1)==0);
+
+
+
+end
+

modules/ded_tkd/TKD_referencecurvefit.m

@@ -0,0 +1,11 @@
+function [normfactor, ref_curve] = TKD_referencecurvefit(rawdata,normalizeddata,angle_matrix,numfiles)
+%fit to the angular data to get the scattering profile
+
+%% Fit reference curve to a smoothing function
+notdeadpixel = find(rawdata(:,1) ~=0);
+
+ref_curve = normalizeddata(:,numfiles);
+
+normfactor = smooth(angle_matrix(:), ref_curve,200,'sgolay',2);
+
+end

modules/ded_tkd/TKD_scattergrid.m

@@ -0,0 +1,9 @@
+function [angle_matrix,xgrid,ygrid] = TKD_scattergrid(PCx,PCy,DD,npix)
+%TKD Scattering Grid
+%calculates the scattered grid pixel positions
+
+[xgrid,ygrid] = meshgrid(1:npix,1:npix);
+angle_matrix = atan(sqrt((xgrid - PCx).^2 + (ygrid - PCy).^2) * 55 / 1000 ./ DD) / pi * 180;
+
+end
+

modules/ded_tkd/fitedge.m

@@ -0,0 +1,6 @@
+function [fitlinex, fitliney] = fitedge(x1,x2,y1,y2)
+fitslope = (y2 - y1) / (x2 - x1);
+
+fitlinex = linspace(x1,x2,(abs(x2 - x1) + 1));
+fitliney = y1 + fitslope .* (fitlinex - x1);
+end

modules/ded_tkd/normalizeto1.m

@@ -0,0 +1,4 @@
+function normalized_matrix = normalizeto1(matrix_in)
+normalized_matrix = matrix_in - min(matrix_in(:)); % make the lowest 0
+normalized_matrix = normalized_matrix ./ max(normalized_matrix(:)); % Normalize
+end

modules/ded_tkd/normalizeto16bit.m

@@ -0,0 +1,7 @@
+function normalized_matrix = normalizeto16bit(matrix_in)
+% normalized_matrix = matrix_in;
+normalized_matrix = matrix_in - min(matrix_in(:)); % make the lowest 0
+normalized_matrix = normalized_matrix ./ max(normalized_matrix(:)); % Normalize
+normalized_matrix = uint16(normalized_matrix * (2^16 - 1)); % Convert to 16 bit
+end
+

modules/rtm/bin/EBSP_gen.m

@@ -1,4 +1,4 @@
-function [ EBSP_out ] = EBSP_gen( EBSP_av,rotmat,screen_int,isHex )
+function [ EBSP_out,r2] = EBSP_gen( EBSP_av,rotmat,screen_int,isHex )
 %EBSP_GEN Generate a pattern for a specific orientation matrix
 
 % This code is copyright Alex Foden and Ben Britton 09/04/2019

start_AstroEBSD.m

@@ -88,6 +88,9 @@
 
 addpath([Astro_FP,link, 'outputs']);
 addpath([Astro_FP,link, 'plot']);
+
+addpath([Astro_FP,link, 'modules',link,'ded_tkd']);
+
 addpath([Astro_FP]);
 
 disp('AstroEBSD file paths loaded');