From 55ffd316aada595d56f799ed2406a4fc898ae77e Mon Sep 17 00:00:00 2001 From: benjaminbritton Date: Sun, 25 Jun 2023 00:11:42 -0700 Subject: [PATCH] TKD_DED Adding routines for transmission kikuchi diffraction / direct electron diffraction analysis paper. Work by Tianbi Zhang and Ben Britton (UBC) --- decks/AstroEBSD_PatternMatch_2021.m | 494 ++++++++++++++++++ decks/TKD_HDR_pattern_astro.m | 374 +++++++++++++ decks/TKD_HDR_stereo_reproject.m | 668 ++++++++++++++++++++++++ gen/loading/h5_pixet.m | 113 ++++ modules/ded_tkd/TKD_cor_saturated.m | 19 + modules/ded_tkd/TKD_flatfield.m | 17 + modules/ded_tkd/TKD_normalize.m | 29 + modules/ded_tkd/TKD_patternfuse.m | 20 + modules/ded_tkd/TKD_referencecurvefit.m | 11 + modules/ded_tkd/TKD_scattergrid.m | 9 + modules/ded_tkd/fitedge.m | 6 + modules/ded_tkd/normalizeto1.m | 4 + modules/ded_tkd/normalizeto16bit.m | 7 + modules/rtm/bin/EBSP_gen.m | 2 +- start_AstroEBSD.m | 3 + 15 files changed, 1775 insertions(+), 1 deletion(-) create mode 100644 decks/AstroEBSD_PatternMatch_2021.m create mode 100644 decks/TKD_HDR_pattern_astro.m create mode 100644 decks/TKD_HDR_stereo_reproject.m create mode 100644 gen/loading/h5_pixet.m create mode 100644 modules/ded_tkd/TKD_cor_saturated.m create mode 100644 modules/ded_tkd/TKD_flatfield.m create mode 100644 modules/ded_tkd/TKD_normalize.m create mode 100644 modules/ded_tkd/TKD_patternfuse.m create mode 100644 modules/ded_tkd/TKD_referencecurvefit.m create mode 100644 modules/ded_tkd/TKD_scattergrid.m create mode 100644 modules/ded_tkd/fitedge.m create mode 100644 modules/ded_tkd/normalizeto1.m create mode 100644 modules/ded_tkd/normalizeto16bit.m diff --git a/decks/AstroEBSD_PatternMatch_2021.m b/decks/AstroEBSD_PatternMatch_2021.m new file mode 100644 index 0000000..1fd8990 --- /dev/null +++ b/decks/AstroEBSD_PatternMatch_2021.m @@ -0,0 +1,494 @@ +clear +home +close all + +% This script was created to kick off discussions at the 2019 EBSD meeting +% Deck to show how pattern matching can work with MTEX + +%% Toolbox locations for AstroEBSD and MTEX + +%these are Ben's settings +InputUser.Astro_loc='C:\Users\tbritton\Documents\GitHub\AstroEBSD_v2\'; %Change this to your AstroEBSD location +% InputUser.Astro_loc='C:\Users\bbrit\Documents\GitHub\AstroEBSD_v2'; %Change this to your AstroEBSD location + +% location_mtex='C:\Users\bbrit\Documents\mtex-5.4.0'; %Change this to where you keep your MTEX folder +location_mtex='E:\Communal_MatlabPlugins\mtex-5.4.0'; + +run(fullfile(InputUser.Astro_loc,'start_AstroEBSD.m')); +run(fullfile(location_mtex,'startup_mtex.m')); + + +%% Load the experimental data we want to work with +InputUser.Phase_Folder = fullfile(InputUser.Astro_loc,'phases\phasefiles'); %location of the pha files +% InputUser.HDF5_folder='C:\Users\bbrit\Documents\Issues\BigFile\'; %Change this to the file location in whch you have saved the example data +InputUser.HDF5_folder='E:\Ben\'; + +InputUser.HDF5_file= 'Demo_Ben.h5'; +InputUser.Phase_Input = {'Ferrite'}; + + +%% Set up the RTM + +%setttings for RTM +RTM.screensize = 128; %size of the library patterns and the resize of the raw EBSPs +RTM.Sampling_Freq=8; %Set the SO(3) sampling freq in degrees +RTM.iterations = 4;%Set the number of iterations to do in the refinement step +RTM.LPTsize = 128; %LPT size used in pixels + +%From AstroEBSD +%background correction +Settings_CorX.gfilt=1; %use a high pass filter (do you mean high pass?) +Settings_CorX.gfilt_s=5; %low pass filter sigma +Settings_CorX.radius=0; %use a radius mask +Settings_CorX.radius_frac=0.85; %fraction of the pattern width to use as the mask +Settings_CorX.hotpixel=0; %hot pixel correction +Settings_CorX.hot_thresh=1000; %hot pixel threshold +Settings_CorX.resize=1; %resize correction +Settings_CorX.size=RTM.screensize; %image height +Settings_CorX.RealBG=0; %use a real BG +Settings_CorX.EBSP_bgnum=30; %number of real pattern to use for BG +Settings_CorX.SquareCrop = 1; %make square the EBSP +Settings_CorX.SplitBG=1; %deal with a split screen + +%% 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 + +[ SettingsXCF, correction, SettingsXCF2 ] = FFT_Filter_settings( RTM.screensize, RTM.LPTsize ); + +%Define all rotation matrices needed in the code +RTM.Rz=@(theta)[cos(theta) sin(theta) 0;-sin(theta) cos(theta) 0;0 0 1]; %z rotation +RTM.Rx=@(theta)[1 0 0;0 cos(theta) sin(theta);0 -sin(theta) cos(theta)]; %x rotation +RTM.Ry=@(theta)[cos(theta) 0 sin(theta);0 1 0; -sin(theta) 0 cos(theta)]; %y rotation + + +%% Read data +[ MapData,MicroscopeData,PhaseData,EBSPData ] = bReadHDF5( InputUser ); + +%load the data from h5 +[MapInfo.MapData,MicroscopeData,PhaseData,MapInfo.EBSPData ]=bReadHDF5( InputUser ); +[MapInfo.Data_InputMap] = EBSD_Map(MapInfo.MapData,MicroscopeData); + +%adjust the pattern centre for the square cropping - this corrects PCx +[ MapInfo.Data_InputMap ] = PC_square( MapInfo.EBSPData, MapInfo.Data_InputMap,Settings_CorX ); + +%get stuff ready for XCF +Detector_tilt = RTM.Rx(MicroscopeData.TotalTilt); + +%% create a reference pattern from the Dynamics bin file +[ ~,~,~,~,~, RTM_info ] = Phase_Builder_RTM( {InputUser.Phase_Input{1}},RTM.Phase_Folder); +[screen_int] = Cube_Generate(RTM_info.bin_file,RTM_info.isHex); + +%% Plot the sphere (because we can) +%create a spherical pattern +[rx,ry,rz]=sphere(1000); +%if you rotate rx, ry and rz you can rotate the pattern +[i_data] = Cube_Sample(rx(:),ry(:),rz(:),screen_int,RTM_info.isHex); + +figure; +surf(rx,ry,rz,reshape(i_data,size(rx,1),size(rx,2)),'EdgeColor','none'); +colormap('gray'); +axis off; axis equal; +axis vis3d; + +%% Generate a simulated pattern +% create a geometry +PC_simulation=[0.5,0.5,0.4]; %try varying this +[EBSD_simulation ] = EBSP_Gnom( RTM,PC_simulation); %you can change PC_in if you want + +%create an orientation +eangs=[25 20 0]*pi/180; %phi1, PHI, phi2 +gmatrix=RTM.Rz(eangs(3))*RTM.Rx(eangs(2))*RTM.Rz(eangs(1)); +[ Pat_sim_eang ] = EBSP_gen( EBSD_simulation,gmatrix*Detector_tilt,screen_int); %generate the EBSP for this iteration + +figure; +subplot(1,1,1); +I1=pPattern(Pat_sim_eang,EBSD_simulation); +title('Simulated pattern'); + +%% Now generate a pattern from the Bruker orientation data + +P_test=10; %the 10th pattern captured +[P_r,P_c]=fPointFind(P_test,MapInfo.Data_InputMap.PMap); + +%read the experimental pattern +[ Pat_exp ] = bReadEBSP(EBSPData,P_test); +[ Pat_exp_cor ] = EBSP_BGCor( Pat_exp,Settings_CorX); + +%Read the PC from Bruker data +PC_in=[MapInfo.Data_InputMap.PCX(P_r,P_c),MapInfo.Data_InputMap.PCY(P_r,P_c),MapInfo.Data_InputMap.DD(P_r,P_c)]; +[EBSD_geom ] = EBSP_Gnom( RTM,PC_in); %you can change PC_in if you want + +%read the Euler angles +eangs=[MapInfo.Data_InputMap.phi1(P_r,P_c),MapInfo.Data_InputMap.PHI(P_r,P_c),MapInfo.Data_InputMap.phi2(P_r,P_c)]*pi/180; +gmatrix=RTM.Rz(eangs(3))*RTM.Rx(eangs(2))*RTM.Rz(eangs(1)); + +%simulate +[ Pat_sim_B] = EBSP_gen( EBSD_geom,gmatrix*Detector_tilt,screen_int); + +figure; +subplot(1,2,1); +I1=pPattern(Pat_exp_cor,EBSD_geom); +title('Experimental pattern'); + +subplot(1,2,2); +I1=pPattern(Pat_sim_B,EBSD_geom); +title('Simulated pattern'); + +%% Refine the Bruker solution + +%prepare the experimental pattern for refinement +[Pat_Ref_r,XCF_data_fill] = refine_prep(Pat_exp_cor,SettingsXCF,RTM); + +%refine +[G_Refined,regout_R] = refine5(Pat_Ref_r,EBSD_geom,EBSD_geom.PC,gmatrix*Detector_tilt,SettingsXCF,screen_int,RTM_info.isHex,RTM); +%[orientation,XCF values]=refine5(Pattern,mean-geometry,updatedPC,orientation,SettingsXCF,screen_int,isHex,RTM) +%XCF values = [x shift, y shift, ~ , normalised chi] + +%update simulation +[ Pat_sim_B_ref ] = EBSP_gen( EBSD_geom,G_Refined,screen_int,screen_int.isHex ); %generate the EBSP for this iteration + +figure; +subplot(2,2,1); +I1=pPattern(Pat_exp_cor,EBSD_geom); +title('Experimental pattern'); + +subplot(2,2,2); +I1=pPattern(Pat_sim_B-Pat_sim_B_ref,EBSD_geom); +title('Difference'); + +subplot(2,2,3); +I1=pPattern(Pat_sim_B,EBSD_geom); +title('Refined Simulation pattern'); + +subplot(2,2,4); +I1=pPattern(Pat_sim_B_ref,EBSD_geom); +title('Bruker Simulation pattern'); + +%% Now do SO(3) searching +%populate fundamental zone of SO(3) +cs_phase=loadCIF(RTM_info.cif_file); +[ library_G ] = SO3_rotmat_gen( cs_phase,RTM.Sampling_Freq); + +%generate the library +[ template_library ] = Library_Gen(EBSD_geom,screen_int,RTM_info.isHex,library_G,2,SettingsXCF); + +%compare with the library +[G_SO3,Library_PH]=fLibraryTest(template_library,library_G,Pat_Ref_r.FFT,SettingsXCF,XCF_data_fill,1); + +%refine +[G_Refined_SO3,regout_R] = refine5(Pat_Ref_r,EBSD_geom,EBSD_geom.PC,G_SO3,SettingsXCF,screen_int,RTM_info.isHex,RTM); + +%update simulation +[ Pat_sim_SO3 ] = EBSP_gen( EBSD_geom,G_SO3,screen_int,screen_int.isHex ); %generate the EBSP for this iteration +%update simulation +[ Pat_sim_SO3r ] = EBSP_gen( EBSD_geom,G_Refined_SO3,screen_int,screen_int.isHex ); %generate the EBSP for this iteration + +figure; +subplot(2,2,1); +I1=pPattern(Pat_exp_cor,EBSD_geom); +title('Experimental pattern'); + +subplot(2,2,2); +I1=pPattern(Pat_sim_SO3-Pat_sim_SO3r,EBSD_geom); +title('Difference'); + +subplot(2,2,3); +I1=pPattern(Pat_sim_SO3,EBSD_geom); +title('SO(3) Template Found'); + +subplot(2,2,4); +I1=pPattern(Pat_sim_SO3r,EBSD_geom); +title('SO(3) Refined'); + +%% Prepare data from a h5 to run as map + +pat_list=1:MicroscopeData.NPoints; +%read the the patterns and prepare them for XCF + +%prepare the data fill for the faster FFT +[~,XCF_data_fill] = refine_prep(zeros(RTM.screensize),SettingsXCF,RTM); + +%create empty arrays in the structure +Pat_Ref_r = struct('FFT', zeros(numel(XCF_data_fill), numel(XCF_data_fill)), 'logp', zeros(RTM.LPTsize,RTM.LPTsize),'pat_in',zeros(RTM.screensize)); + +%build the structure with the input data +parfor p=1:numel(pat_list) + [ Pat_exp ] = bReadEBSP(EBSPData,p); + [ Pat_exp_cor ] = EBSP_BGCor( Pat_exp,Settings_CorX); + [ Pat_Ref_r(p) ] = refine_prep(Pat_exp_cor,SettingsXCF,RTM); %note this stores every pattern, the FFT and the log-polar +end + +%% Extract useful mapdata +map_pcx=MapInfo.MapData.PCX; +map_pcy=MapInfo.MapData.PCY; +map_pcz=MapInfo.MapData.DD; + +map_phi1=map_pcx; +map_PHI=map_pcx; +map_phi2=map_pcx; + +g_bruker=zeros(3,3,numel(map_pcx)); +g_bruker_det=zeros(3,3,numel(map_pcx)); + +for p=1:numel(map_pcx) + %find the point in the map + [P_r,P_c]=fPointFind(p,MapInfo.Data_InputMap.PMap); + + %extract the PC list + map_pcx(p)=MapInfo.Data_InputMap.PCX(P_r,P_c); + map_pcy(p)=MapInfo.Data_InputMap.PCY(P_r,P_c); + map_pcz(p)=MapInfo.Data_InputMap.DD(P_r,P_c); + + %extract the orientations + map_phi1(p)=MapInfo.Data_InputMap.phi1(P_r,P_c); + map_PHI(p)=MapInfo.Data_InputMap.PHI(P_r,P_c); + map_phi2(p)=MapInfo.Data_InputMap.phi2(P_r,P_c); + + %calculate the orientation matrix + eangs=[map_phi1(p),map_PHI(p),map_phi2(p)]*pi/180; %bruker stores in degrees + g_bruker(:,:,p)=RTM.Rz(eangs(3))*RTM.Rx(eangs(2))*RTM.Rz(eangs(1)); + g_bruker_det(:,:,p)=g_bruker(:,:,p)*Detector_tilt; %convert to the detector +end + +%% Check that this is reasonalbe + +p=100; + +%generate the geometry of the pattern +[EBSD_geom ] = EBSP_Gnom( RTM,[map_pcx(p) map_pcy(p) map_pcz(p)]); %you can change PC_in if you want + +%update simulation for this pattern +[ Pat_sim_bruker ] = EBSP_gen( EBSD_geom,g_bruker_det(:,:,p),screen_int,screen_int.isHex ); %generate the EBSP for this iteration + +%perform a refinement +tic +[G_refined_bruker_test,~] = refine5(Pat_Ref_r(p),EBSD_geom,EBSD_geom.PC,g_bruker_det(:,:,p),SettingsXCF,screen_int,RTM_info.isHex,RTM); +toc + +%update simulation for the refined pattern +[ Pat_sim_brukerref ] = EBSP_gen( EBSD_geom,G_refined_bruker_test,screen_int,screen_int.isHex ); %generate the EBSP for this iteration + +figure; +subplot(2,2,1); +I1=pPattern(Pat_Ref_r(p).pat_in,EBSD_geom); +title('Experimental pattern'); + +subplot(2,2,2); +I1=pPattern(Pat_sim_bruker,EBSD_geom); +title('Simulation - Bruker measurement'); + +subplot(2,2,3); +I1=pPattern(Pat_sim_brukerref,EBSD_geom); +title('Simulation - Bruker measurement'); + +subplot(2,2,4); +I1=pPattern(Pat_sim_brukerref-Pat_sim_bruker,EBSD_geom); +title('Difference'); + +%% Bruker refinement of patterns + +%remove the original patterns - good for memory, but plotting will need +%them to be read again.. +if isfield(Pat_Ref_r,'pat_in') + Pat_Ref_r=rmfield(Pat_Ref_r,'pat_in'); +end + +%preallocate variables +G_refined=zeros(3,3,numel(pat_list)); +regout_Grefined=zeros(4,numel(pat_list)); +disp(['Matching started']); +tic +parfor p=1:numel(pat_list) + + %update geometry + [EBSD_geom ] = EBSP_Gnom( RTM,[map_pcx(p) map_pcy(p) map_pcz(p)]); %you can change PC_in if you want + + %perform a refinement + [G_refined(:,:,p),regout_Grefined(:,p)] = refine5(Pat_Ref_r(p),EBSD_geom,EBSD_geom.PC,g_bruker_det(:,:,p),SettingsXCF,screen_int,RTM_info.isHex,RTM); + + %{ + if p/100 == round(p/100) + disp(['Pattern ' int2str(p) ' run, of ' int2str(numel(pat_list))]); + end + %} +end +toc +disp(['All Patterns Matched']); + +%% create the maps from this data + +%rotate the detector data into the sample frame +G_refined_sample=zeros(3,3,numel(pat_list)); + +for p=1:numel(pat_list) + G_refined_sample(:,:,p)=inv(G_refined(:,:,p)*inv(Detector_tilt)); +end + +%create the MTEX variables +setMTEXpref('xAxisDirection','west'); +setMTEXpref('zAxisDirection','outOfPlane'); + +prop.x = double(MapInfo.MapData.XSample); %have to transpose - not sure why... +prop.y = double(MapInfo.MapData.YSample); +prop.quality=regout_Grefined(4,:); %PH from refinement +prop.quality=prop.quality(:); + +ori_Brefined = ... + rotation('Matrix',G_refined_sample(:,:,:)); +ebsd_Bruker = EBSD(ori_Brefined, ones(size(ori_Brefined)),{'notIndexed',cs_phase},prop); + +figure; +colorKey = ipfHSVKey(cs_phase); + +colorKey.inversePoleFigureDirection = xvector; +plot(ebsd_Bruker,colorKey.orientation2color(ebsd_Bruker.orientations)) +nextAxis + +colorKey.inversePoleFigureDirection = yvector; +plot(ebsd_Bruker,colorKey.orientation2color(ebsd_Bruker.orientations)) +nextAxis + +colorKey.inversePoleFigureDirection = zvector; +plot(ebsd_Bruker,colorKey.orientation2color(ebsd_Bruker.orientations)) +nextAxis + +%{ +%% +%get the mean PC for the map +Mean_PC=[mean(MapInfo.Data_InputMap.PCX(:)),mean(MapInfo.Data_InputMap.PCY(:)),mean(MapInfo.Data_InputMap.DD(:))]; +%generate the library for this PC +[ template_library ] = Library_Gen(EBSD_geom,screen_int,RTM_info.isHex,library_G,2,SettingsXCF); + + + + + +%% Map based matching +G_SO3_map=zeros(3,3,numel(pat_list)); +G_Refined_SO3_map=zeros(3,3,numel(pat_list)); +G_Refined_Bruker_map=zeros(3,3,numel(pat_list)); +Library_PH=zeros(numel(pat_list),1); +regout_R=zeros(4,numel(pat_list)); +regout_Bruker=zeros(4,numel(pat_list)); + + + + +%% Refine the Bruker map + +for p=1:numel(pat_list) + + G_Bruker_map=gmatrix; + PC_in=[map_pcx(p),map_pcy(p),map_pcz(p)]; + [G_Refined_Bruker_map(:,:,p),regout_Bruker(:,p)] = refine5(Pat_Ref_r(p),EBSD_geom,PC_in,gmatrix*Detector_tilt,SettingsXCF,screen_int,RTM_info.isHex,RTM); +end + + + +%% SO3 Search +parfor p=1:numel(pat_list) + %compare with the library + [G_SO3_map(:,:,p),Library_PH(p)]=fLibraryTest(template_library,library_G,Pat_Ref_r(p).FFT,SettingsXCF,XCF_data_fill,1); + + PC_in=[map_pcx(p),map_pcy(p),map_pcz(p)]; + + %refine + [G_Refined_SO3_map(:,:,p),regout_R(:,p)] = refine5(Pat_Ref_r(p),EBSD_geom,PC_in,G_SO3_map(:,:,p),SettingsXCF,screen_int,RTM_info.isHex,RTM); +end + +%% Image one pattern and check solutions + +p=4000; %pattern number + +%extract the PC for this pattern +PC_in=[map_pcx(p),map_pcy(p),map_pcz(p)]; +[EBSD_geom ] = EBSP_Gnom( RTM,PC_in); %you can change PC_in if you want + +%extract the pattern +Pat_EXPp=Pat_Ref_r(p).pat_in; + +%create SO3 simulation +[ Pat_simp_SO3 ] = EBSP_gen( EBSD_geom,G_SO3_map(:,:,p),screen_int,screen_int.isHex ); %generate the EBSP for this iteration +%create refined simulation +[ Pat_simp_SO3r ] = EBSP_gen( EBSD_geom,G_Refined_SO3_map(:,:,p),screen_int,screen_int.isHex ); %generate the EBSP for this iteration + +%plot the figure +f1=figure; +subplot(2,2,3); +I1=pPattern(Pat_EXPp,EBSD_geom); +title('Experimental pattern'); + +subplot(2,2,2); +I1=pPattern(Pat_simp_SO3-Pat_simp_SO3r,EBSD_geom); +title('Difference'); + +subplot(2,2,1); +I1=pPattern(Pat_simp_SO3,EBSD_geom); +title('SO(3) Template Found'); + +subplot(2,2,4); +I1=pPattern(Pat_simp_SO3r,EBSD_geom); +title('SO(3) Refined'); +f1.Position=[100 100 800 400]; + +%% Convert to MTEX orientation data +G_Samp_ref=zeros(3,3,numel(pat_list)); +G_Samp_SO3=zeros(3,3,numel(pat_list)); + +% prop.quality_so3=Library_PH; %PH from SO3 +% +%convert to the sample frame +for p=1:numel(pat_list) + G_Samp_ref(:,:,p)=inv(G_Refined_SO3_map(:,:,p)*inv(Detector_tilt)); + G_Samp_SO3(:,:,p)=inv(G_SO3_map(:,:,p)*inv(Detector_tilt)); +end + +setMTEXpref('xAxisDirection','west'); +setMTEXpref('zAxisDirection','outOfPlane'); + +% build the coordinate maps +prop.x = double(MapInfo.MapData.XSample); %have to transpose - not sure why... +prop.y = double(MapInfo.MapData.YSample); +prop.quality=regout_R(4,:); %PH from refinement +prop.quality_so3=Library_PH; %PH from SO3 + +ori_ref = ... + rotation('Matrix',G_Samp_ref(:,:,:)); +ori_so3 = ... + rotation('Matrix',G_Samp_SO3(:,:,:)); + +ebsd_ref = EBSD(ori_ref, ones(size(ori_ref)),{'notIndexed',cs_phase},'Options',prop); +ebsd_so3 = EBSD(ori_so3, ones(size(ori_ref)),{'notIndexed',cs_phase},'Options',prop); + +%% Plot some data +figure; +colorKey = ipfHSVKey(cs_phase); +colorKey.inversePoleFigureDirection = xvector; + +plot(ebsd_so3,colorKey.orientation2color(ebsd_so3.orientations)) +nextAxis +plot(ebsd_ref,colorKey.orientation2color(ebsd_ref.orientations)) + +figure; +colorKey = ipfHSVKey(cs_phase); +colorKey.inversePoleFigureDirection = yvector; + +plot(ebsd_so3,colorKey.orientation2color(ebsd_so3.orientations)) +nextAxis +plot(ebsd_ref,colorKey.orientation2color(ebsd_ref.orientations)) + +figure; +colorKey = ipfHSVKey(cs_phase); +colorKey.inversePoleFigureDirection = zvector; + +plot(ebsd_so3,colorKey.orientation2color(ebsd_so3.orientations)) +nextAxis +plot(ebsd_ref,colorKey.orientation2color(ebsd_ref.orientations)) + +figure; +plot(ebsd_so3,ebsd_so3.prop.quality_so3); +nextAxis; +plot(ebsd_ref,ebsd_so3.prop.quality); + +%} \ No newline at end of file diff --git a/decks/TKD_HDR_pattern_astro.m b/decks/TKD_HDR_pattern_astro.m new file mode 100644 index 0000000..073a041 --- /dev/null +++ b/decks/TKD_HDR_pattern_astro.m @@ -0,0 +1,374 @@ +%% Preamble +% TKD_HDR_Pattern.m - a script in AstroEBSD Package to generate a fused +% high dynamic range on-axis transmission kikuchi diffraction pattern (TKP) +% +% Originally scripted by Tianbi Zhang in November 2022 with major rework in +% April 2023 (structure and variables) and June 2023 (functions) +% Annotated in April-June 2023 +% For how this script works, please refer to our paper +% +% Pre-requisite: MATLAB Image Processing & Statistics and Machine Learning +% Toolbox, Signal Processing Toolbox +% +% Additional requisites (does not affect this code, but helps you check +% file structure of the h5) +% (a) h5 file viewer to check h5 file structure +% (b) A few new functions, especially yaxis rotation, is optimized in +% MATLAB R2023a. Not using this version though should not affect its +% functionality. +% +% Inputs: on-axis transmission Kikuchi patterns captured by a TimePIX3 +% direct detector in .h5 format. The patterns can be saved as individual +% frames within a single .h5 file, or as a single integrated frame. There +% are two different functions to import the data to the workspace. +% +% Outputs: fused high dynamic range pattern (raw and flatfielded), +% corresponding 2D FFT spectra (modulus) and assorted plots. +% +% The standard version of the script generates similar plots as Figure 3-5 +% of the original paper (you will need a simulated pattern for Fig.5.) +% +% You can customize this script to produce other plots - please work on a +% copy of this code for your own pattern, and acknowledge the original +% work if possible. +% Fig. 6-8 are generated by a separate reprojection code which can be found +% in the AstroEBSD repository as well. + +%% +home; +close all; +clear; + +%% Load astro +% location_astro='C:\Users\billy\Documents\GitHub\AstroEBSD_v2\'; %Change this to your AstroEBSD location +location_astro='C:\Users\benja\OneDrive\Documents\GitHub\AstroEBSD_v2'; +run(fullfile(location_astro,'start_AstroEBSD.m')); + +%% Settings +export_option = false; % select true if you want image outputs +multframe_option = false; % true if the data are stored in individual frames within a h5 file +multframe_num = 50; % the number of frames in a multi-frame h5 file, +% or the number of frames integrated if your raw data is collected as sum of frames. + +% InputUser.Astro_loc='C:\Users\billy\Documents\MATLAB\AstroEBSD_v2-live'; +% run(fullfile(InputUser.Astro_loc,'start_AstroEBSD.m')); % run AstroEBSD + +%% Initial PC info +DD = 10.535; %% CL in mm - you need to measure this in the SEM. +PCx = 136; %%PCx, PCy in pixel units - you can read from the image +PCy = 146; % 149 + +%Saturation of the camera, in counts +sat_index = 1022; + +% I added an npix definition +npix = 254; % final dimension of the pattern (a square) in pixel units - we remove the edge two pixels, as they are larger than the rest + +Settings_Cor.gfilt=1; %gaussian (high frequency filter) - 1 = on, 0 = off +%you can switch off the gaussian high pass filter if you need/want +Settings_Cor.gfilt_s=15; %low pass filter sigma, in pixels, e.g. remove things that vary bigger than 15 pixels +Settings_Cor.size = [npix, npix]; +%% Import Data + +% Change file path based on your settings +path1= 'C:\Users\benja\OneDrive\Documents\MATLAB\Al_HDR\Demo'; +% path1= 'C:\Users\billy\OneDrive - UBC\PhD\TKD\Al_HDR\Demo'; + +% Please fill in the file names in decreasing order of exposure time. The +% file corresponding to the baseline exposure time t0 (the longest in the set) should be the +% first one. +h5filenames = {'spot1_50f_01s.h5','spot1_50f_001s.h5', 'spot1_50f_0005s.h5','spot1_50f_0001s_Event.h5'}; +sampleID = ["0.1s Exposure", "0.01s Exposure", "0.005s Exposure", "0.001s Exposure"]; + +[rawdata,numfiles] = h5_pixet(h5filenames,path1,multframe_option,multframe_num); + +% Signal Saturation Correction & Normalization +[fillingdomains, satdomains] = TKD_cor_saturated(rawdata,numfiles,sat_index); +% added sat_index input in function dfefinition +% added satdomains in function output + +[normalizeddata] = TKD_normalize(fillingdomains,rawdata,numfiles); + +% note: for the "demo" data, due to a mistake in collection, the event and +% iToT subframes are separated into 2 files. Thi resuslts in the extra +% if-else on lines 17-23 of h5_pixet. For other data, you can just use line +% 19 and discard the rest of the sub if-loop. + +%% Plot raw signal versus scattering angle +PC_TKD=[PCx/npix,1-PCy/npix,1000*DD/(55*npix)]; % TZ - I think this should be right - can you have a think whether this should be 256 of 254? + +% Detector.screensize = npix; +[ TKD_Geometry ] = EBSP_Gnom( Settings_Cor, PC_TKD ); + +angle_matrix=atand(sqrt(TKD_Geometry.ypts_screen.^2+TKD_Geometry.xpts_screen.^2)); +% angle_matrix(find(angle_matrix == min(angle_matrix(:)))) = 0; + +% added 2 additional outputs so the fitted curve will pe properly plotted. +[normfactor, ref_curve] = TKD_referencecurvefit(rawdata,normalizeddata,angle_matrix,numfiles); + +% Stitch the pattern +[fusedpattern] = TKD_patternfuse(rawdata,normalizeddata,fillingdomains); + +%% Flat field +% normfactor = feval(fitresult2, angle_matrix); +normfactor = reshape(normfactor, [npix npix]); +fused_flat = fusedpattern ./ normfactor; + +% If your fused pattern still shows uneven background, consider adding a +% Gaussian filter, for example: +fused_gaussian = imgaussfilt(fused_flat, 15); +fused_flat_gaus = fused_flat ./ fused_gaussian; +% Flat field using the fitted radial scattering profile +[fused_flat,fusedpattern,normfactor, ref_curve] = TKD_flatfield(rawdata,normalizeddata,angle_matrix,fillingdomains,numfiles); + +% [fused_flat,Settings_Cor ] = EBSP_BGCor( fused_flat,Settings_Cor); +[~,Settings_Cor ] = EBSP_BGCor( fused_flat,Settings_Cor); +% This step is completely optional + +%% Visualization and FFT + +fused_fft = fft2(fusedpattern); +fused_flat_fft = fft2(fused_flat); +fused_gaus_fft = fft2(fused_flat_gaus); + +figure; +subplot(2,2,1); +imagesc(fusedpattern); axis xy; axis image; colormap('gray');axis off; +title("Raw fused"); +subplot(2,2,2); +imagesc(fused_flat); axis xy; axis image; colormap('gray');axis off; +title("Flattened fused"); +subplot(2,2,3); +imagesc(log10(abs(fftshift(fused_fft)))); axis xy; axis image; colormap('gray');axis off; +subplot(2,2,4); +imagesc(log10(abs(fftshift(fused_flat_fft)))); axis xy; axis image; colormap('gray');axis off; + +%% Assorted Plots + +figure; +subplot(1,3,1); %Figure 1: raw data as a function of scattering angle +hold on; +for i=1:numfiles +scatter(angle_matrix(:), rawdata(:,i), 1, 'o', 'filled'); +end +grid on; +xlabel("Scattering Angle (degrees)"); +ylabel(["Event"; "Count"],'Rotation',0); +legend(sampleID); +subplot(1,3,2); %Figure 2: normalized data and the conformity to the reference curve +hold on; +for i=1:numfiles +scatter(angle_matrix(:), normalizeddata(:,i), 1, 'o', 'filled'); +end +grid on; +xlabel("Scatter Angle (degrees)"); +ylabel(["Weighted"; "Event"; "Count"],'Rotation',0); +% legend(sampleID); +subplot(1,3,3); %Figure 3: Plot fit with the reference curve +plot(angle_matrix(:), ref_curve,'.','Color',[0 0 1],'MarkerSize',10); +hold on; +plot(angle_matrix(:), normfactor(:),'-', 'LineWidth',2, 'Color',[1 0 0]); +legend('Reference Curve', 'Smoothed Reference Curve', 'Location', 'NorthEast'); +% Label axes +xlabel( 'Scatter Angle (degree)', 'Interpreter', 'none' ); +ylabel( ["Weighted"; 'Event'; 'Count'], 'Rotation',0); +grid on; + +%% Export images: 16-bit TIFF + +fusedpattern_16 = normalizeto16bit(fusedpattern); +fused_flat_16 = normalizeto16bit(fused_flat); +fused_gaus_16 = normalizeto16bit(fused_flat_gaus); + +% note: below are log10 fft spectrum +fused_fft_16 = normalizeto16bit(log10(abs(fftshift(fused_fft)))); +fused_flat_fft_16 = normalizeto16bit(log10(abs(fftshift(fused_flat_fft)))); +fused_gaus_fft_16 = normalizeto16bit(log10(abs(fftshift(fused_gaus_fft)))); + + +if export_option == true + imwrite(uint16(flipud(fusedpattern_16)), fullfile(path1,'fused_raw.tif')); + imwrite(uint16(flipud(fused_flat_16)), fullfile(path1,'fused_flat.tif')); + imwrite(uint16(flipud(fused_gaus_16)), fullfile(path1,'fused_flat_gaus.tif')); + + imwrite(uint16(flipud(fused_fft_16)), fullfile(path1,'fused_raw_fft.tif')); + imwrite(uint16(flipud(fused_flat_fft_16)), fullfile(path1,'fused_flat_fft.tif')); + imwrite(uint16(flipud(fused_gaus_fft_16)), fullfile(path1,'fused_flat_gaus_fft.tif')); +end + +%% Side scripts for other plots +%%% SideScript A: patterns vs exposure time (not normalized) +numpats = numfiles+1; +sampleID{end+1} = 'Raw fused'; +%this will be exposure time longest to shortest, then fused + +%% Figure: patterns at different exposure times, short to long, and gaussian +figure; +hold on; +for k=1:numfiles +pattern = reshape(rawdata(:,k),[254 254]); +subplot(1,numpats,k); +imagesc(pattern); axis xy; axis image; axis off; colormap('gray'); colorbar; +title(sampleID{k}); +end +subplot(1,numpats, numpats); +imagesc(fusedpattern); axis xy; axis image; axis off; colormap('gray'); colorbar; +title(sampleID{end}); + +figure; +hold on; +for k=1:numfiles +pattern = reshape(rawdata(:,k),[254 254]); +pattern_gaus = pattern ./ normfactor; +pattern_gaus(satdomains{k})=0; + +subplot(1,numpats,k); +imagesc(pattern_gaus); axis xy; axis image; axis off; colormap('gray'); +end +subplot(1,numpats,numpats); +imagesc(fused_flat); axis xy; axis image; axis off; colormap('gray'); + +%% intensity histogram (not in the paper) +figure; +histogram(fusedpattern_16, 200,'EdgeColor','none','Normalization','count'); +hold on; +histogram(fused_gaus_16, 200,'EdgeColor','none','Normalization','count'); +xlabel("Normalized grayscale"); +ylabel("Pixel count"); +legend("Raw fused","Background corrected"); +grid on; + +%% angle contours + corrected intensity plot +f3 = figure; +h2 = axes; +p2 = pPattern(fused_flat,TKD_Geometry); +h3 = axes; +[p3,h4] = contour(TKD_Geometry.xpts_screen,TKD_Geometry.ypts_screen,angle_matrix,[10 20 30 40 50],'ShowText','on', 'linewidth', 2); +set(h3,'color','none','visible','off','fontsize',22); +axis equal; +clabel(p3,h4,'color','yellow', 'fontsize',20); +set(h3,'ydir', 'normal'); +linkaxes([h2 h3]); + + +figure; +plot(angle_matrix(:),fused_flat(:),'.','MarkerSize',10,'Color',"#77AC30"); +xlabel("Scattering Angle (degrees)"); +ylabel(["Corrected"; "Intensity"; "(a.u.)"], 'Rotation',0); +grid on; + +%% filling zone contours (for easier visualization) +fillingcontour = zeros(254^2,1); + +for i=1:numfiles + fillingcontour(fillingdomains{i}) = i; +end + +fillingcontour = reshape(fillingcontour, 254,254); + +figure; +contour(fillingcontour); +axis xy; axis image; + +%% Side Script 6: line profile across a band - need to plot corrected images +% and relative intensity; almost perfectly perpendicular profiling +% This section will not work if you don't have a simulated pattern as the reference. +% The simulated pattern for the demo pattern is included in the raw +% dataset. +newsampleID = sampleID; +newsampleID(end) = "Fused pattern"; +extendedsampleID = newsampleID; +extendedsampleID(end+1) = "Simulated Pattern"; + +% read the simulated pattern +opts = delimitedTextImportOptions("NumVariables", 254, "Encoding", "UTF-8"); +% Specify range and delimiter +opts.DataLines = [1, Inf]; +opts.Delimiter = "\t"; +opts.VariableNames = ["VarName1", "VarName2", "VarName3", "VarName4", "VarName5", "VarName6", "VarName7", "VarName8", "VarName9", "VarName10", "VarName11", "VarName12", "VarName13", "VarName14", "VarName15", "VarName16", "VarName17", "VarName18", "VarName19", "VarName20", "VarName21", "VarName22", "VarName23", "VarName24", "VarName25", "VarName26", "VarName27", "VarName28", "VarName29", "VarName30", "VarName31", "VarName32", "VarName33", "VarName34", "VarName35", "VarName36", "VarName37", "VarName38", "VarName39", "VarName40", "VarName41", "VarName42", "VarName43", "VarName44", "VarName45", "VarName46", "VarName47", "VarName48", "VarName49", "VarName50", "VarName51", "VarName52", "VarName53", "VarName54", "VarName55", "VarName56", "VarName57", "VarName58", "VarName59", "VarName60", "VarName61", "VarName62", "VarName63", "VarName64", "VarName65", "VarName66", "VarName67", "VarName68", "VarName69", "VarName70", "VarName71", "VarName72", "VarName73", "VarName74", "VarName75", "VarName76", "VarName77", "VarName78", "VarName79", "VarName80", "VarName81", "VarName82", "VarName83", "VarName84", "VarName85", "VarName86", "VarName87", "VarName88", "VarName89", "VarName90", "VarName91", "VarName92", "VarName93", "VarName94", "VarName95", "VarName96", "VarName97", "VarName98", "VarName99", "VarName100", "VarName101", "VarName102", "VarName103", "VarName104", "VarName105", "VarName106", "VarName107", "VarName108", "VarName109", "VarName110", "VarName111", "VarName112", "VarName113", "VarName114", "VarName115", "VarName116", "VarName117", "VarName118", "VarName119", "VarName120", "VarName121", "VarName122", "VarName123", "VarName124", "VarName125", "VarName126", "VarName127", "VarName128", "VarName129", "VarName130", "VarName131", "VarName132", "VarName133", "VarName134", "VarName135", "VarName136", "VarName137", "VarName138", "VarName139", "VarName140", "VarName141", "VarName142", "VarName143", "VarName144", "VarName145", "VarName146", "VarName147", "VarName148", "VarName149", "VarName150", "VarName151", "VarName152", "VarName153", "VarName154", "VarName155", "VarName156", "VarName157", "VarName158", "VarName159", "VarName160", "VarName161", "VarName162", "VarName163", "VarName164", "VarName165", "VarName166", "VarName167", "VarName168", "VarName169", "VarName170", "VarName171", "VarName172", "VarName173", "VarName174", "VarName175", "VarName176", "VarName177", "VarName178", "VarName179", "VarName180", "VarName181", "VarName182", "VarName183", "VarName184", "VarName185", "VarName186", "VarName187", "VarName188", "VarName189", "VarName190", "VarName191", "VarName192", "VarName193", "VarName194", "VarName195", "VarName196", "VarName197", "VarName198", "VarName199", "VarName200", "VarName201", "VarName202", "VarName203", "VarName204", "VarName205", "VarName206", "VarName207", "VarName208", "VarName209", "VarName210", "VarName211", "VarName212", "VarName213", "VarName214", "VarName215", "VarName216", "VarName217", "VarName218", "VarName219", "VarName220", "VarName221", "VarName222", "VarName223", "VarName224", "VarName225", "VarName226", "VarName227", "VarName228", "VarName229", "VarName230", "VarName231", "VarName232", "VarName233", "VarName234", "VarName235", "VarName236", "VarName237", "VarName238", "VarName239", "VarName240", "VarName241", "VarName242", "VarName243", "VarName244", "VarName245", "VarName246", "VarName247", "VarName248", "VarName249", "VarName250", "VarName251", "VarName252", "VarName253", "VarName254"]; +opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double", "double"]; + +% Specify file level properties +opts.ExtraColumnsRule = "ignore"; +opts.EmptyLineRule = "read"; +% Import the data +simpattern_table = fullfile(path1, 'sim_pattern.csv'); + +simpattern = readtable(simpattern_table, opts); +simpattern = table2array(simpattern); + +simpattern = simpattern + 0.5; % offset to move the band intensities up +clear opts; + +figure; +imagesc(simpattern); axis xy; axis image; colormap('gray'); axis off; + +[line_1_x, line_1_y] = fitedge(28, 57, 34, 31); % these are read from the pattern manually +[line_2_x, line_2_y] = fitedge(200,222,167,178); + +figure; +subplot(2,2,[1 3]); +hold on; +imagesc(fused_flat); axis xy; axis image; colormap('gray'); axis off; +hold on; +plot(line_1_x, line_1_y,'Color','cyan','LineWidth',2); +plot(line_2_x, line_2_y,'Color','magenta','LineWidth',2); +hold off; + + +for k=1:numfiles +thispattern = reshape(normalizeddata(:,k),[254 254]); +bgcorpattern = thispattern ./ normfactor; +subplot(2,2,2); +hold on; +lineprof1 = improfile(bgcorpattern, line_1_x, line_1_y, 'nearest'); +plot(lineprof1,'LineWidth',2); +subplot(2,2,4); +hold on; +lineprof2 = improfile(bgcorpattern, line_2_x, line_2_y, 'nearest'); +plot(lineprof2,'LineWidth',2); +end +subplot(2,2,2); +lineprof = improfile(fused_flat, line_1_x,line_1_y, 'nearest'); +plot(lineprof,'LineWidth',2); +lineprof = improfile(simpattern, line_1_x,line_1_y, 'nearest'); +plot(lineprof,'LineWidth',2); +legend(extendedsampleID); +grid on; +xlabel("Pixel Position"); +ylabel("Corrected Intensity (a.u.)"); +title("Line 1 (Cyan)"); +subplot(2,2,4); +lineprof = improfile(fused_flat, line_2_x, line_2_y, 'nearest'); +plot(lineprof,'LineWidth',2); +lineprof = improfile(simpattern, line_2_x-2,line_2_y, 'nearest'); +plot(lineprof,'LineWidth',2); +legend(extendedsampleID); +grid on; +xlabel("Pixel Position"); +ylabel("Corrected Intensity (a.u.)"); +title("Line 2 (Magenta)"); + +%% Helper functions + + +% write a dead pixel correction program + +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 + +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 + +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 diff --git a/decks/TKD_HDR_stereo_reproject.m b/decks/TKD_HDR_stereo_reproject.m new file mode 100644 index 0000000..aaf7e3d --- /dev/null +++ b/decks/TKD_HDR_stereo_reproject.m @@ -0,0 +1,668 @@ +%% stereo_reproject.m +% Created by T. Ben Britton in June 2023 +% Edited and annotated by Tianbi Zhang + +% This script takes a Kikuchi pattern, index it, reproject it onto the +% diffraction sphere and create a stereogram. + +% This script accompanies the following article: +% "Multi-exposure diffraction pattern fusion applied to enable wider-angle +% transmission Kikuchi diffraction with direct electron detectors" +% Tianbi Zhang, T. Ben Britton + +% Requirements: +% (1) MATLAB toolboxes: image processsing, statistics and machine learning, +% parallel computing +% (2) AstroEBSD package - this script is available as a part of the latest +% AstroEBSSD distribution. +% (3) MTEX toolbox (https://mtex-toolbox.github.io/) + +%% Start of the script +close all; +clear; +home; + +%% start astro & mtex - please edit the paths +location_astro='C:\Users\benja\OneDrive\Documents\GitHub\AstroEBSD\'; +location_mtex='C:\Users\benja\OneDrive\Documents\MATLAB\mtex-5.4.0'; +path1= 'C:\Users\benja\OneDrive\Documents\MATLAB\Al_HDR\Demo'; + + +run(fullfile(location_astro,'start_AstroEBSD.m')); +run(fullfile(location_mtex,'startup.m')); + + +%% load the pattern + Figure 1 + + +pat1=double(flipud(imread(fullfile(path1,'stitched_flat_gaus.tif')))); % this is the demo TKP +figure; +imagesc(pat1); axis image; axis xy; colormap('gray'); + +%% Enter pattern information and set up indexing + +% Euler angles (Bunge) and PC - they must be very accurate for the +% reprojection. +eangs_pat1=[140.4739,23.8782,246.8339];%demo +pc_pat1=[0.5454366,0.4350097,0.7578597];% demo + +% eangs_pat1=[277.5706,40.1078,68.0877];%ap1 +% pc_pat1=[0.5308719,0.4267852,0.6428683];% ap1 + +% eangs_pat1=[351.0974,35.4311,350.7168];%ap2 +% pc_pat1=[0.5402837,0.4316793,0.6424058];% ap2 + +MicroscopeData.TotalTilt=0; +%Set up the radon transform peak finder +Settings_Rad.theta_range=[-10 180 1]; %theta min, theta max, theta step - in degrees + +%peak huntfor the Hough map +Settings_Rad.max_peaks=12; %max number of peaks to return +Settings_Rad.num_peak=20; %number of peaks to search for - peaks will be rejected +Settings_Rad.theta_search_pix=6; %search size in theta steps +Settings_Rad.rho_search_per=0.2; %radon search in fractions +Settings_Rad.min_peak_width=0.002; %min rseperation of the peak width, in pixels + + +% Normalise intensities +[EBSP_One.PatternIn,Settings_Cor ] = EBSP_BGCor( pat1,[]); + +%Define all rotation matrices needed in the code +RTM.Rz=@(theta)[cos(theta) sin(theta) 0;-sin(theta) cos(theta) 0;0 0 1]; %z rotation +RTM.Rx=@(theta)[1 0 0;0 cos(theta) sin(theta);0 -sin(theta) cos(theta)]; %x rotation +RTM.Ry=@(theta)[cos(theta) 0 sin(theta);0 1 0; -sin(theta) 0 cos(theta)]; %y rotation + +%% index in AstroEBSD (Hough) + +EBSP_Fe_One.PatternIn=pat1; +EBSP_Fe_One.PC=pc_pat1; + +InputUser.Phase_Folder = fullfile(location_astro,'phases\phasefiles'); +InputUser.Phase_Input = {'Al'}; +[ Crystal_UCell,Crystal_Family,Crystal_LUT,Settings_LUT,Phase_Num ] = Phase_Builder( InputUser.Phase_Input,InputUser.Phase_Folder ); + +%Radon transform +[ EBSP_Fe_One.Peak_Centre,EBSP_Fe_One.Single.Peak_Set_All,EBSP_Fe_One.Peak_Set_All,... + EBSP_Fe_One.R_EBSP,EBSP_Fe_One.R_Edge,EBSP_Fe_One.R_rho,EBSP_Fe_One.R_theta ] ... + = EBSP_RadHunt( EBSP_Fe_One.PatternIn,Settings_Rad); + +% Convert the bands to normal space +[ EBSP_Fe_One.nhat_gnom] = EBSP_NormConv( EBSP_Fe_One.Peak_Centre,size(EBSP_Fe_One.PatternIn),EBSP_Fe_One.PC); + +R_x=@(theta)[1 0 0;0 cos(theta) sin(theta);0 -sin(theta) cos(theta)]; +rot_det=R_x(MicroscopeData.TotalTilt); + +% Index this pattern +[EBSP_Fe_One.rotdata{1},EBSP_Fe_One.banddata{1}]=EBSP_Index(EBSP_Fe_One.nhat_gnom,Crystal_LUT{1},Settings_LUT{1}.thresh_trig,Crystal_UCell{1},rot_det); + +%generate the geometry +[ EBSP_Fe_One.PatternGeometry ] = EBSP_Gnom( Settings_Cor,EBSP_Fe_One.PC ); + +EBSP_OneFigure=Plot_SinglePattern(EBSP_Fe_One,Crystal_UCell,Crystal_LUT,1); + +%% Simulate the pattern + +eangs=eangs_pat1*pi/180; +PC_simulation=pc_pat1; +pattern_info=struct; +pattern_info.size=size(pat1); +Detector_tilt = RTM.Rx(MicroscopeData.TotalTilt); + +% InputUser.Phase_Folder +[ ~,~,~,~,~, RTM_info ] = Phase_Builder_RTM( {InputUser.Phase_Input{1}},[location_astro 'phases']); +[screen_int] = Cube_Generate(RTM_info.bin_file,RTM_info.isHex); + +[EBSD_simulation ] = EBSP_Gnom( pattern_info,PC_simulation); %you can change PC_in if you want +gmatrix=RTM.Rz(eangs(3))*RTM.Rx(eangs(2))*RTM.Rz(eangs(1)); +[ Pat_sim_eang ] = EBSP_gen( EBSD_simulation,gmatrix*Detector_tilt,screen_int); %generate the EBSP for this iteration + +[Pat_sim_astroind]= EBSP_gen( EBSD_simulation,EBSP_Fe_One.rotdata{1}.detector,screen_int); %generate the EBSP for this iteration + +%% refine + Figure 2 + +RTM.screensize = 254; %size of the library patterns and the resize of the raw EBSPs +RTM.Sampling_Freq=8; %Set the SO(3) sampling freq in degrees +RTM.iterations = 10;%Set the number of iterations to do in the refinement step +RTM.LPTsize = 254; %LPT size used in pixels +[ SettingsXCF, correction, SettingsXCF2 ] = FFT_Filter_settings( RTM.screensize, RTM.LPTsize ); + +[ Pat_Ref_r ] = refine_prep(EBSP_Fe_One.PatternIn,SettingsXCF,RTM); +%perform a refinement +[G_refined,regout_Grefined] = refine5(Pat_Ref_r,EBSD_simulation,EBSD_simulation.PC,gmatrix,SettingsXCF,screen_int,RTM_info.isHex,RTM); +[ Pat_sim_ref ] = EBSP_gen( EBSD_simulation,G_refined,screen_int); %generate the EBSP for this iteration + +figure; +subplot(1,2,1); +pPattern(Pat_sim_ref,EBSD_simulation); +subplot(1,2,2); +pPattern(EBSP_Fe_One.PatternIn,EBSD_simulation); + +%% Figure 3 - fitted/indexed patterns +figure; + +subplot(1,3,1); +I1=pPattern(pat1,EBSD_simulation); +title('Experimental pattern'); + +subplot(1,3,2); +I1=pPattern(Pat_sim_eang,EBSD_simulation); +title('Simulated pattern - Dynamics'); + +subplot(1,3,3); +I1=pPattern(Pat_sim_astroind,EBSD_simulation); +title('Simulated pattern - AstroH'); + +%% Now try to plot a sterogram of the dynamically simulated pattern + Figure 4 +num_stereo=2001; +stereo_range=1; +stereo_xl=linspace(-stereo_range,stereo_range,num_stereo); +stereo_yl=linspace(-stereo_range,stereo_range,num_stereo); + +[stereo_x,stereo_y]=meshgrid(stereo_xl,stereo_yl); + +stereo_r=stereo_x.^2+stereo_y.^2; + +stereo_nu=2./(stereo_r+1); + +stereo_vx=stereo_nu.*stereo_x; +stereo_vy=stereo_nu.*stereo_y; +stereo_vz=-1+stereo_nu; + +stereo_vl=sqrt(stereo_vx.^2+stereo_vy.^2+stereo_vz.^2); +stereo_vx=stereo_vx./stereo_vl; +stereo_vy=stereo_vy./stereo_vl; +stereo_vz=stereo_vz./stereo_vl; + + +%sample from the sphere +[stereo_pat] = Cube_Sample(stereo_vx(:),stereo_vy(:),stereo_vz(:),screen_int,0); +% figure; imagesc(stereo_xl,stereo_yl,reshape(stereo_pat,num_stereo,num_stereo)); +% axis image; axis xy; colormap('gray'); + +%% Now take the experimental pattern and generate a look up, so we can do the spherical sampling... +% Figure 5 +% rotmat=EBSP_Fe_One.rotdata{1}.detector; +rotmat=gmatrix*Detector_tilt; +[Pat_sim_astroind,r_exp]= EBSP_gen( EBSD_simulation,rotmat,screen_int); %generate the EBSP for this iteration + +r_lambda=1./(r_exp(:,3)+1); + +%convert from r_exp to the stereogram +P_rx=r_exp(:,1).*r_lambda; +P_ry=r_exp(:,2).*r_lambda; + +pat1_plot=pat1; +pat1_plot=pat1_plot/std(pat1(:)); +pat1_plot=pat1_plot-mean(pat1_plot(:)); + +% create a intensity normalized stereogram +stereo_pat_plot=stereo_pat; +stereo_pat_plot=stereo_pat/std(stereo_pat_plot(:)); +stereo_pat_plot=stereo_pat_plot-mean(stereo_pat_plot(:)); +% +% figure; +% hist(pat1_plot(:),100,'r'); +% hold on; +% hist(stereo_pat_plot(:),100,'b'); + +figure; imagesc(stereo_xl,stereo_yl,reshape(stereo_pat_plot,num_stereo,num_stereo)); +axis image; axis xy; colormap('gray'); axis off; +hold on; +% here the pattern is plotted as a scatter plot whose marker color scale +% corresponds to the grayscale value. This is still OK at low mag. +% scatter(P_rx,P_ry,1,pat1_plot(:)*0.8); + +%% use the interpolant to generate the pattern +% create the interpolant +stero_interp_exp=scatteredInterpolant(P_rx(:),P_ry(:),pat1_plot(:),'natural','none'); + +%read the intepolant +stereo_remap=stero_interp_exp(stereo_x(:),stereo_y(:)); + +stereo_remap_2d=reshape(stereo_remap,num_stereo,num_stereo); +stereo_mask_2d=isnan(stereo_remap_2d); + +figure; imagesc(stereo_xl,stereo_yl,reshape(stereo_pat_plot,num_stereo,num_stereo)); +axis image; axis xy; colormap('gray'); +hold on; +i1=imagesc(stereo_xl,stereo_yl,stereo_remap_2d); +axis image; axis xy; colormap('gray'); +i1.AlphaData=1-stereo_mask_2d; + +% plot the outline of the detector (experimental pattern) +line_top=[EBSD_simulation.xpts_screen(:,1),EBSD_simulation.ypts_screen(:,1)]; +line_bottom=flipud([EBSD_simulation.xpts_screen(:,end),EBSD_simulation.ypts_screen(:,end)]); +line_right=[EBSD_simulation.xpts_screen(end,:);EBSD_simulation.ypts_screen(end,:)]'; +line_left=flipud([EBSD_simulation.xpts_screen(1,:);EBSD_simulation.ypts_screen(1,:)]'); + +line_allx=[line_top(:,1);line_right(:,1);line_bottom(:,1);line_left(:,1)]; +line_ally=[line_top(:,2);line_right(:,2);line_bottom(:,2);line_left(:,2)]; +line_allz=line_ally*0+1; + +%rotate these points around the sphere +line_r = [line_allx(:), line_ally(:), line_allz(:)].*1./sqrt((line_allx(:).^2+line_ally(:).^2+1)); +line_r2 = line_r*rotmat'; + +%convert these points to the stereogram +line_r_lambda=1./(line_r2(:,3)+1); + +%convert from r_exp to the stereogram +line_P_rx=line_r2(:,1).*line_r_lambda; +line_P_ry=line_r2(:,2).*line_r_lambda; + +%plot the blue box +plot(line_P_rx,line_P_ry,'b','LineWidth',2); + +% +%do the fundamental zone +%write three axes, and construct the lines that go between one and the next +ax1=[0 0 1]; +ax2=[1 1 2]; +ax3=[1 0 1]; + +%normalize length +ax1=ax1./norm(ax1); +ax2=ax2./norm(ax2); +ax3=ax3./norm(ax3); +nl=100; + +%inscribe along the great circles and equispace +ang_12=acos(dot(ax1,ax2)); +ax_12=cross(ax1,ax2); +ax_12=ax_12/norm(ax_12); +ax_2p=cross(ax_12,ax1); +ax_2p=ax_2p/norm(ax_2p); +sample_spacing12=linspace(0,ang_12,nl); +line_12=ax1'*cos(sample_spacing12)+ax_2p'*sin(sample_spacing12); + +ang_23=acos(dot(ax2,ax3)); +ax_23=cross(ax2,ax3); +ax_23=ax_23/norm(ax_23); +ax_3p=cross(ax_23,ax2); +ax_3p=ax_3p/norm(ax_3p); +sample_spacing23=linspace(0,ang_23,nl); +line_23=ax2'*cos(sample_spacing23)+ax_3p'*sin(sample_spacing23); + +ang_31=acos(dot(ax3,ax1)); +ax_31=cross(ax3,ax1); +ax_31=ax_31/norm(ax_31); +ax_1p=cross(ax_31,ax3); +ax_1p=ax_1p/norm(ax_1p); +sample_spacing31=linspace(0,ang_31,nl); +line_31=ax3'*cos(sample_spacing31)+ax_1p'*sin(sample_spacing31); + +%projected on stereogram +line_set=[line_12';line_23';line_31']; +line_set_lambda=1./(line_set(:,3)+1); +line_set_rx=line_set(:,1).*line_set_lambda; +line_set_ry=line_set(:,2).*line_set_lambda; + +plot(line_set_rx(1:end),line_set_ry(1:end),'y', 'LineWidth',2); + + +% %plot the individual segments +% plot(line_set_rx(1:end),line_set_ry(1:end),'r'); +% scatter(line_set_rx(201),line_set_ry(201),'m'); +% plot(line_set_rx(nl+(1:nl)),line_set_ry(nl+(1:nl)),'g'); +% plot(line_set_rx(2*nl+(1:nl)),line_set_ry(2*nl+(1:nl)),'b'); + +%% lets inflate this up to the full pattern +nline=20; +line_P_rxn=zeros(4*nline,24); +line_P_ryn=zeros(4*nline,24); +stereo_remap_maskp=zeros(num_stereo,num_stereo,24); +stereo_remap_2dnp=stereo_remap_maskp; + +cs_al=loadCIF('Al-Aluminum.cif'); + +EBSP_av=EBSD_simulation; +[Pat_sim_astroind,r_exp]= EBSP_gen( EBSD_simulation,rotmat,screen_int); %generate the EBSP for this iteration + +r_n=zeros(size(r_exp,1),size(r_exp,2),24); +n_sym=24; +for n=1:n_sym + sym_matrix=cs_al.rot(n).matrix; + r_n(:,:,n) = r_exp*sym_matrix; +end + +stereo_x_ok=stereo_x(:); +stereo_y_ok=stereo_y(:); + +% calculate and plot each symmetrically equivalent pattern in the +% stereogram +parfor p=1:24 + + x_lookup=EBSP_av.xpts_screen; + y_lookup=EBSP_av.ypts_screen; + + x_lookup=x_lookup(:,2:end-1); + y_lookup=y_lookup(:,2:end-1); + pat1_plot_n=pat1_plot(:,2:end-1); + + x_lookup=x_lookup(2:end-1,1); + y_lookup=y_lookup(2:end-1,1); + pat1_plot_n=pat1_plot_n(2:end-1,:); + + r = [x_lookup(:), y_lookup(:), y_lookup(:)*0+1].*1./sqrt((x_lookup(:).^2+y_lookup(:).^2+1)); + r_exp2 = r*rotmat'; + + r_3=r_n(:,3,p); + r_3(r_3<0)=-r_3(r_3<0); + r_lambda=1./(r_3+1); + + %convert from r_exp to the stereogram + P_rxn=r_n(:,1,p).*r_lambda; + P_ryn=r_n(:,2,p).*r_lambda; + + % interpolate + stero_interp_exp_n=scatteredInterpolant(P_rxn(:),P_ryn(:),pat1_plot(:),'natural','none'); + + %plot the outline of the detector + line_top=[EBSD_simulation.xpts_screen(:,1),EBSD_simulation.ypts_screen(:,1)]; + line_bottom=flipud([EBSD_simulation.xpts_screen(:,end),EBSD_simulation.ypts_screen(:,end)]); + line_right=[EBSD_simulation.xpts_screen(end,:);EBSD_simulation.ypts_screen(end,:)]'; + line_left=flipud([EBSD_simulation.xpts_screen(1,:);EBSD_simulation.ypts_screen(1,:)]'); + + s2=linspace(1,254,nline); + s2(2:end-1)=round(s2(2:end-1)); + + line_allx=[line_top(s2,1);line_right(s2,1);line_bottom(s2,1);line_left(s2,1)]; + line_ally=[line_top(s2,2);line_right(s2,2);line_bottom(s2,2);line_left(s2,2)]; + line_allz=line_ally*0+1; + + %rotate these points around the sphere + line_r = [line_allx(:), line_ally(:), line_allz(:)].*1./sqrt((line_allx(:).^2+line_ally(:).^2+1)); + line_r2 = line_r*rotmat'*cs_al.rot(p).matrix; + line_r2_3=line_r2(:,3); + line_r2_3(line_r2_3<0)=-line_r2_3(line_r2_3<0); + %convert these points to the stereogram + line_r_lambda=1./(line_r2_3+1); + + %convert from r_exp to the stereogram + line_P_rx=line_r2(:,1).*line_r_lambda; + line_P_ry=line_r2(:,2).*line_r_lambda; + + %create the blank frame + stereo_remap_2dn=zeros(num_stereo,num_stereo); + stereo_remap_mask=ones(num_stereo,num_stereo); +num_el=numel(stereo_remap_mask); + % %find the points in the detector frame + % [in,on] = inpolygon(stereo_x_ok,stereo_y_ok,line_P_rx,line_P_ry); + in=true(num_el,1); + stereo_remap_mask_0=stereo_remap_mask*0; + %read the intepolant + stereo_remap_n=stero_interp_exp_n(stereo_x_ok(in),stereo_y_ok(in)); + stereo_remap_2dn(in)=stereo_remap_n; + stereo_remap_mask(in)=stereo_remap_mask_0(in); + + % stereo_remap_2dn=reshape(stereo_remap_n,num_stereo,num_stereo); + % stereo_mask_2dn=isnan(stereo_remap_2dn); + + % Figure 6,7 + %plot the base dynamical pattern + figure; imagesc(stereo_xl,stereo_yl,reshape(stereo_pat_plot,num_stereo,num_stereo)); + axis image; axis xy; colormap('gray'); + hold on; + + %plot the remapped experiment + i1=imagesc(stereo_xl,stereo_yl,stereo_remap_2dn); + axis image; axis xy; colormap('gray'); + i1.AlphaData=1-stereo_remap_mask; + + %plot the blue box + plot(line_P_rx,line_P_ry,'b'); + xlim([-1 1]); + ylim([-1 1]); + + title(int2str(p)); + + line_P_rxn(:,p)=line_P_rx; + line_P_ryn(:,p)=line_P_ry; + stereo_remap_maskp(:,:,p)=stereo_remap_mask; + stereo_remap_2dnp(:,:,p)=stereo_remap_2dn; +end + +%% Individual plots of the 24 equivalent patterns reprojected +% (this is commented out by default to reduce window count) + +% for p=1:24 +% %plot the base dynamical pattern +% figure; imagesc(stereo_xl,stereo_yl,reshape(stereo_pat_plot,num_stereo,num_stereo)); +% axis image; axis xy; colormap('gray'); +% hold on; +% +% %plot the remapped experiment +% i1=imagesc(stereo_xl,stereo_yl,stereo_remap_2dnp(:,:,p)); +% axis image; axis xy; colormap('gray'); +% i1.AlphaData=1-stereo_remap_maskp(:,:,p); +% +% %plot the blue box +% plot(line_P_rxn(:,p),line_P_ryn(:,p),'b'); +% xlim([-1 1]); +% ylim([-1 1]); +% title(int2str(p)); +% end + +%% Plot all outlines of the reprojected experimental patterns +% on a simulated pattern + +figure; imagesc(stereo_xl,stereo_yl,reshape(stereo_pat_plot,num_stereo,num_stereo)); +axis image; axis xy; colormap('gray'); +hold on; + +for p=1:24 + %plot the blue box + plot(line_P_rxn(:,p),line_P_ryn(:,p),'b', 'LineWidth',1.5); + xlim([-1 1]); + ylim([-1 1]); + title(int2str(p)); +end + +%% Plot the sterogram and try to work out the points for this projection + +%read the experimental pattern +pat1_plot2=pat1_plot; +pat_c=pat1_plot2(:); + +pat_maps=zeros(2001,2001,24); +pat_mapn=pat_maps; + +%pick the symemtry +for p=1:24 + %extract the points of this pattern + %rotate these points of the detector around the sphere + box_x=EBSD_simulation.xpts_screen(:); + box_y=EBSD_simulation.ypts_screen(:); + box_z=box_y*0+1; + + % box_r=[box_x,box_y,box_z]; + box_r = [box_x(:), box_y(:), box_z(:)].*1./sqrt((box_x(:).^2+box_y(:).^2+1)); + + box_r2=box_r*rotmat'*cs_al.rot(p).matrix; + box_r2x=box_r2(:,1); + box_r2y=box_r2(:,2); + box_r2z=box_r2(:,3); + + %now find those which are in the positive and negative hemisphere ('p' + %and 'n' points) based upon their z coords + box_r2z_pz=find(box_r2z>0); + box_r2z_nz=find(box_r2z<0); + + %subselect the XYZ positives + box_r2xp=box_r2x(box_r2z_pz); + box_r2yp=box_r2y(box_r2z_pz); + box_r2zp=box_r2z(box_r2z_pz); + + %subselect the XYZ negatives + box_r2xn=-box_r2x(box_r2z_nz); + box_r2yn=-box_r2y(box_r2z_nz); + box_r2zn=-box_r2z(box_r2z_nz); + + %extract the pattern data (we are going to use scatter for a lazy plot) + pat_cp=pat_c(box_r2z_pz); + pat_cn=pat_c(box_r2z_nz); + + %sterographic projection + box_r2_lambdap=1./(box_r2zp+1); + box_r2_xp=box_r2xp.*box_r2_lambdap; + box_r2_yp=box_r2yp.*box_r2_lambdap; + + %find the boundary of the point cloud + boundary_r2_p=boundary(box_r2_xp,box_r2_yp); + boundary_r2_px=box_r2_xp(boundary_r2_p); + boundary_r2_py=box_r2_yp(boundary_r2_p); + + %sterographic projection + box_r2_lambdan=1./(box_r2zn+1); + box_r2_xn=box_r2xn.*box_r2_lambdan; + box_r2_yn=box_r2yn.*box_r2_lambdan; + + %find the boundary of the point cloud + boundary_r2_n=boundary(box_r2_xn,box_r2_yn); + boundary_r2_nx=box_r2_xn(boundary_r2_n); + boundary_r2_ny=box_r2_yn(boundary_r2_n); + + %find the pounts within each boundary + sx=stereo_x; + sy=stereo_y; + + [in_p]=inpolygon(sx,sy,boundary_r2_px(:),boundary_r2_py(:)); + [in_n]=inpolygon(sx,sy,boundary_r2_nx(:),boundary_r2_ny(:)); + + x_lookup=EBSD_simulation.xpts_screen; + y_lookup=EBSD_simulation.ypts_screen; + + r = [x_lookup(:), y_lookup(:), y_lookup(:)*0+1].*1./sqrt((x_lookup(:).^2+y_lookup(:).^2+1)); + r_exp2 = r*rotmat'*cs_al.rot(p).matrix; + %create the interpolant + r_3=r_exp2(:,3); + r_lambda=1./(r_3+1); + + %convert from r_exp to the stereogram + P_rxn=r_exp2(:,1).*r_lambda; + P_ryn=r_exp2(:,2).*r_lambda; + + % generate interpolate + stero_interp_exp_s=scatteredInterpolant(P_rxn(:),P_ryn(:),pat1_plot(:),'natural','none'); + + stero_interp_exp_sp=scatteredInterpolant(box_r2_xp(:),box_r2_yp(:),pat_cp(:),'natural','none'); + stero_interp_exp_sn=scatteredInterpolant(box_r2_xn(:),box_r2_yn(:),pat_cn(:),'natural','none'); + + P_stero_xp=sx(in_p); + P_stero_yp=sy(in_p); + Pat_interp_p=stero_interp_exp_sp(P_stero_xp,P_stero_yp); + + P_stero_xn=sx(in_n); + P_stero_yn=sy(in_n); + Pat_interp_n=stero_interp_exp_sn(P_stero_xn,P_stero_yn); + + pat_map=0*sx; + pat_num=0*sx; + pat_num(in_p)=pat_num(in_p)+1; + pat_num(in_n)=pat_num(in_n)+1; + + pat_map(in_p)=Pat_interp_p; + pat_map(in_n)=Pat_interp_n; + + pat_maps(:,:,p)=pat_map; + pat_mapn(:,:,p)=pat_num; +end + +%% Normalize the overlaid stereogram and plot +pat_maps_num=sum(pat_mapn,3); +pat_maps_sum=sum(pat_maps,3); +pat_maps_norm=pat_maps_sum./pat_maps_num; + +%% final stereogram +figure; +subplot(1,3,1); +imagesc(stereo_xl,stereo_yl,pat_maps_sum); % total stereogram +axis image; axis xy; colormap('gray'); axis off; +subplot(1,3,2); +imagesc(stereo_xl,stereo_yl,pat_maps_norm); % stereogram normalized by weight +axis image; axis xy; colormap('gray'); axis off; +% sum/weight hemisphere +subplot(1,3,3); imagesc(stereo_xl,stereo_yl,pat_maps_num); % overlapped patterns +axis image; axis xy; colormap('gray'); axis off; + +%% Difference between simulated and reprojected patterns +% normalize the reprojected dynamics pattern +pat_maps_dyn=reshape(stereo_pat_plot,num_stereo,num_stereo); +pat_maps_dynn=(pat_maps_dyn-mean(pat_maps_dyn(:)))./std(pat_maps_dyn(:)); +pat_maps_normn=(pat_maps_norm-mean(pat_maps_norm(:)))./std(pat_maps_norm(:)); + +figure; imagesc(stereo_xl,stereo_yl,pat_maps_norm-pat_maps_dyn); +axis image; axis xy; colormap('gray'); +hold on; +clim([-6 6]); + +%% Write some sterograms to tifs + + +imwrite(uint16(flipud(normalizeto16bit(pat_maps_norm))), fullfile(path1,'results','stereo_stack_exp.tif')); +imwrite(uint16(flipud(normalizeto16bit(pat_maps_dyn))), fullfile(path1,'results','stereo_stack_dyn.tif')); + + +%% +line_r = [line_allx(:), line_ally(:), line_allz(:)].*1./sqrt((line_allx(:).^2+line_ally(:).^2+1)); +line_r2 = line_r*rotmat'*cs_al.rot(p).matrix; +line_r2_3=line_r2(:,3); + +line_p=line_r2_3*0; +line_y1=line_r2_3*0+1; +line_p(line_r2_3<0)=line_y1(line_r2_3<0); +line_p=logical(line_p); +line_n=logical(1-line_p); + +line_r2_3(line_r2_3<0)=-line_r2_3(line_r2_3<0); + + +%convert these points to the stereogram +line_r_lambda=1./(line_r2_3+1); + +%convert from r_exp to the stereogram +line_P_rx=line_r2(:,1).*line_r_lambda; +line_P_ry=line_r2(:,2).*line_r_lambda; + +%plot the base dynamical pattern +figure; imagesc(stereo_xl,stereo_yl,reshape(stereo_pat_plot,num_stereo,num_stereo)); +axis image; axis xy; colormap('gray'); +hold on; + +%plot the remapped experiment +i1=imagesc(stereo_xl,stereo_yl,stereo_remap_2dnp(:,:,p)); +axis image; axis xy; colormap('gray'); +i1.AlphaData=1-stereo_remap_maskp(:,:,p); + +%plot the blue box +plot(line_P_rx,line_P_ry,'b'); + +%plot the positive +line_P_rx_p=line_P_rx(line_p); +line_P_ry_p=line_P_ry(line_p); +plot(line_P_rx_p,line_P_ry_p,'y'); + +%plot the negative +line_P_rx_n=line_P_rx(line_n); +line_P_ry_n=line_P_ry(line_n); +plot(line_P_rx_n,line_P_ry_n,'m'); + +[in_l]=inpolygon(line_P_rx_n,line_P_ry_n,line_P_rx_p,line_P_ry_p); +line_P_rx_n1=line_P_rx_n(in_l); +line_P_ry_n1=line_P_ry_n(in_l); + +line_P_rx_n1_r=sqrt(line_P_rx_n1.^2+line_P_ry_n1.^2); + +line_P_rx_n1=line_P_rx_n1./line_P_rx_n1_r; +line_P_ry_n1=line_P_ry_n1./line_P_rx_n1_r; +plot(line_P_rx_n1,line_P_ry_n1,'w','linewidth',3); + +xlim([-1 1]); +ylim([-1 1]); +title(int2str(p)); + +% plot a circle in green +thetac=linspace(0,2*pi,360); +x_c=1*sin(thetac); +y_c=1*cos(thetac); +plot(x_c,y_c,'g') \ No newline at end of file diff --git a/gen/loading/h5_pixet.m b/gen/loading/h5_pixet.m new file mode 100644 index 0000000..b7a05ed --- /dev/null +++ b/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 \ No newline at end of file diff --git a/modules/ded_tkd/TKD_cor_saturated.m b/modules/ded_tkd/TKD_cor_saturated.m new file mode 100644 index 0000000..7d0eb7a --- /dev/null +++ b/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 + diff --git a/modules/ded_tkd/TKD_flatfield.m b/modules/ded_tkd/TKD_flatfield.m new file mode 100644 index 0000000..e1dd632 --- /dev/null +++ b/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 + diff --git a/modules/ded_tkd/TKD_normalize.m b/modules/ded_tkd/TKD_normalize.m new file mode 100644 index 0000000..8e88c87 --- /dev/null +++ b/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 + diff --git a/modules/ded_tkd/TKD_patternfuse.m b/modules/ded_tkd/TKD_patternfuse.m new file mode 100644 index 0000000..3ec9c4a --- /dev/null +++ b/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 + diff --git a/modules/ded_tkd/TKD_referencecurvefit.m b/modules/ded_tkd/TKD_referencecurvefit.m new file mode 100644 index 0000000..f2a63cf --- /dev/null +++ b/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 \ No newline at end of file diff --git a/modules/ded_tkd/TKD_scattergrid.m b/modules/ded_tkd/TKD_scattergrid.m new file mode 100644 index 0000000..3b4f133 --- /dev/null +++ b/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 + diff --git a/modules/ded_tkd/fitedge.m b/modules/ded_tkd/fitedge.m new file mode 100644 index 0000000..d1b2af9 --- /dev/null +++ b/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 diff --git a/modules/ded_tkd/normalizeto1.m b/modules/ded_tkd/normalizeto1.m new file mode 100644 index 0000000..a1b4480 --- /dev/null +++ b/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 diff --git a/modules/ded_tkd/normalizeto16bit.m b/modules/ded_tkd/normalizeto16bit.m new file mode 100644 index 0000000..60b887b --- /dev/null +++ b/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 + diff --git a/modules/rtm/bin/EBSP_gen.m b/modules/rtm/bin/EBSP_gen.m index 623351f..1cf2d04 100644 --- a/modules/rtm/bin/EBSP_gen.m +++ b/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 diff --git a/start_AstroEBSD.m b/start_AstroEBSD.m index 7e8ba67..f4fd09f 100644 --- a/start_AstroEBSD.m +++ b/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');