<jupyter_start><jupyter_text>2021 PARADIM Summer School  4D STEM Data Analysis T5: Zhaslan, <NAME>  Table of Contents [1. Introduction](intro)[2. Import packages and data](load)[3. Browse 4D STEM data](browse)&emsp;[a) Interactive/Dynamic Viewer](3a)&emsp;&emsp;[checkpoint 1](cp1)&emsp;[b) Static Viewer](3b)&emsp;[c) ADF image](3c) [4. Exit Wave Power Cepstrum (EWPC)](EWPC)&emsp;[a) Perform cepstral transform](ct)&emsp;&emsp;[checkpoint 2](cp2)&emsp;[b) Select EWPC peaks](peak)&emsp;[c) Calculate spot maps](spotmap) [5. Spot maps](bspot) - Map of length and angle of the EWPC spots with respect to the center.&emsp;&emsp;[checkpoint 3](cp3)[6. Strain maps](bstrain) - Map of the strain matrix components.&emsp;&emsp;[checkpoint 4](cp4)[7. (Optional) Segmentation](segment) - Segment the region of interest.  0. MAKE A COPY OF THIS NOTEBOOKBefore working through this notebook, please make a copy and save it with the same name plus your initials. Continue this tutorial using your copy.  1. Introduction Typical imaging methods in STEM sum up all the electrons falling on the detector area to form the image. Depending on range of scattered angles you collect on the detector, you can get different types of contrast on your image - phase contrast, diffraction contrast, Z contrast, etc. With recent developments in pixelated detectors, it has become possible to efficiently collect the whole diffraction pattern at each scan position. The resulting datasets have 4 dimensions - 2 dimensions in real space along the scan directions, and a 2 dimensional diffraction pattern for each real space pixel, hence the technique is popularly called 4D-STEM.However, 4D-STEM datasets are huge (a single 128 x 128 scan on EMPAD gives 1 GB data) and hence it is important to be familiar with the different tools to analyse the data. This notebook is aimed at introducing you to the different kinds of imaging modes possible with 4D-STEM data, and what information these different modes provide about the sample.As an example of another application of 4D-STEM, we also demonstrate strain mapping on PTO/DSO sample, using the Exit Wave Power Cepstral (EWPC) transform method developed at Cornell - [Padgett et al., Ultramicroscopy 214, 112994 (2019)](https://doi.org/10.1016/j.ultramic.2020.112994). In summary, this notebook introduces: &emsp;(1) different imaging modes and&emsp;(2) strain mappingNote that this notebook is color-coded. When your input is needed, the markdown cell will indicate the part with this orange color. Checkpoints that will guide you through the notebook will be shown with this turquoise color. [go to top](toc) 2. Import packages and functions<jupyter_code>#Custom functions used in this script, which are stored in fun.py
from fun import * #* means import all functions within the script

#*%matplotlib notebook* mode makes the cells to be interactive.
#If cells don't seem to be interactive, restart the kernel and run this cell first.
%matplotlib notebook<jupyter_output><empty_output><jupyter_text>**Optional**: To learn what kind of input a function needs or see what the function does, type the name of the function and a question mark, then run the cell like following:```pythonfunctionName?```Or you can open the `fun.py` and check the details.  Load 4D-STEM Data (.raw format)Locate the 4D-STEM data you want to look at in a file browser. Copy and paste the directory of the file  to assign the variable ***`inputDir`*** and the file name  as ***`fname`***. (The path and filename for the sample 4D data is already filled out).<jupyter_code>inputDir='../Data/' #directory where the 4D STEM data file is.
fname= 'PTODSO_small_128x130x64x64_float_x64_y64.raw' #name of the 4D-STEM data
saveDir='./'#./ means current directory; where this Jupyter notebook is.spotmaps file will be saved here.

data4d=load_raw_to_dp(inputDir+fname,130,128)<jupyter_output><empty_output><jupyter_text>The above code directly reads in the 4D STEM data from the raw file output acquired with EMPAD 1. The file name needs to end with the scan size in real space - for this file, 64 steps in x and y directions are indicated by '...x64_y64.raw'. For EMPAD 1, the data for each scan position is saved as 130 x 128 pixels with some metadata. If you need to work with 4D STEM data acquired on other detectors, you will need to change the arguments passed to the function `load_raw_to_dp` in `fun.py`. Check the dimensions<jupyter_code>(rx,ry,sx,sy)=data4d.shape #order: N_x1,N_x2,N_k1,N_k2
(rx,ry,sx,sy)<jupyter_output><empty_output><jupyter_text>[go to top](toc) 3. Browse the 4D-STEM Data  3a) Interactive/Dynamic version Run the following cell. On the real space image (image on the left), click and drag over the region of interest (from where you want to get the strain map). Then, the average diffraction pattern (DP) from the region inside the shaded red box will be shown on the right. Diffraction space image will ideally show clear diffraction pattern. Draw a box over a diffraction spot. Then, the real space image will show the dark field image corresponding to the selected diffraction spot. Check which region lights up in the real space image. Here, you can see if the corresponding part is tilted or has a different orientation from others, etc. Also, by selecting a diffraction spot, you can distinguish your sample from vacuum/amorphous region easily.**Goals**: Browse the 4D-STEM Data to (1) check the DP throughout the sample; (2) identify the coordinates of the central beam (coordinates are shown at the bottom right of the figure window); (3) identify diffraction spots that give more information about the sample.<jupyter_code>plt.close('all')
guis=browser(data4d)<jupyter_output><empty_output><jupyter_text>CHECKPOINT 1**1)** Draw a box on the real space image - what do you see on the diffraction space image? Move the box around, what changes and why? **2)** Now do the same on the diifraction pattern - try to interpret the different images you get as you move around in diffraction space. **3)** What can we learn from selecting particular diffraction disks to form the image? **4)** Also, notice how there is a contrast reversal between the BF and DF images. Why?  [go to top](toc) 3b) Static version **Optional**: You can also get snapshots of the 4D data by specifying both real and diffraciton space coordinates.To do so, assign `stptReal, stptDiff` with desired coordinates. `statDisp` function will be useful for comparing images.(Format: `statDisp(dataset,dataset, coordinates in diffraction space, coordinates in real space)`)<jupyter_code>stptDiff = int(rx/2); stptReal = int(sx/2);#starting points
statDisp(data4d,[60,60],[20,50]);<jupyter_output><empty_output><jupyter_text>[go to top](toc) c) Virtual ADF image **Optional**: Specify the radii (`r0` & `r1`) for the virtual ADF mask. The virtual ADF image generated during this process will be used for selecting the region of interest for the EWPC strain mapping.**Optional Task** : Change the radii values to form a BF and ABF image.<jupyter_code>r0=0 #inner radius
r1=1 #outer radius

### display diffraction pattern at real space position 
i=30; j=30 #coordinates in real space to get the diffraction image shown on the right
disp_haadf(data4d,[r0,r1],i,j)<jupyter_output><empty_output><jupyter_text>[go to top](toc) 4. EWPC 4a) Perform cepstral transform Theory behind: [Padgett et al., Ultramicroscopy 214, 112994 (2019)](https://doi.org/10.1016/j.ultramic.2020.112994).<jupyter_code>ewpc=convert_dp_to_ewpc(data4d)<jupyter_output><empty_output><jupyter_text>[go to top](toc) 4b) Select EWPC peaks for strain map Just like you did with the 4D browser,1) Draw a rectangle in the leftmost window to indicate the area of interest2) Draw another box around a diffraction peak in the reciprocal space (cepstral space in this case). Make sure the peak is close to the center of the rectangle.3) The rightmost window will show a zoomed in image of the box drawn in the middle window. Select just the peak (the brightest spot within the diffraction disk) area, by drawing another box.4) Click on `Save Results` button to add the region of interest and window location for the strain analysis.5) Repeat steps 2-4 at least once (you need at least two g vectors to calculate strain)<jupyter_code>guis,wins,rois=browser_with_peak_selection(ewpc)<jupyter_output><empty_output><jupyter_text>[go to top](toc) CHECKPOINT 2**1)** Remember that EWPC picks out the periodicity in the diffraction patterns. Based on your understanding of the method, try to think about the optimal experimental conditions for collecting the 4D-STEM data - think about convergence angle and camera length.**2)** What do the different peaks in cepstral space correspond to?**Optional**: Since you already know how to do simulations from Tutorial 3 - you can also try to simulate diffraction patterns for this sample, then calculate the EWPC pattern, and try to correlate the different projected lattice spacings to the EWPC peaks.**3)** Similar to the 4D browser, try to select different areas on the real space image and see how the EWPC peaks change. And try to understand the images you get when you form images from different cepstral peaks.**4)** Compare the cepstral DF images to the DF images you obtained from the diffraction pattern. What differences do you see? **Hint**: Think about diffraction contrast, tilt and thickness artefacts.  Confirm selection on the dark field cepstral STEM imageCheck the region you selected.<jupyter_code>roi=rois[0]
wins=np.array(wins)
show_roi(ewpc,roi,wins)<jupyter_output><empty_output><jupyter_text>Make sure the selected peaks are shown below within the red boxes.<jupyter_code>data4d_roi,ewpc_img=show_wins(data4d,wins,roi)<jupyter_output><empty_output><jupyter_text>[go to top](toc) 4c) Calculate spot maps Calculate the peak position of the above selected EWPC spots with sub-pixel precision using the Nelder Mead algorithm. After the peak position is found, the length and angle of the peak from the center is found.**Optional**: lower the ***`tol`***, which is tolerance, to get better precision at a cost of increased processing time. ***`tol`*** of 1e-3 to 1e-4 is a good place to start.Available peak finding ***`method`*** from fastest to slowest : ***`L-BFGS-B`***, ***`Powell`***,***`Nelder-Mead`*** (default),***`TNC`***<jupyter_code>spotMaps_upd = get_spotMaps(data4d_roi,wins,tol=1e-4,method='L-BFGS-B')<jupyter_output><empty_output><jupyter_text>Save ***`spotMaps`***<jupyter_code>spotMaps_upd['wins']=np.array(wins)
spotMaps_upd['roi']=rois[0]
pickle.dump(spotMaps_upd,open(saveDir+'/spotMaps.dat','wb'))<jupyter_output><empty_output><jupyter_text>[go to top](toc) 5. Spot Maps Plot of the vector length and angle between the peak positions of the selected EWPC spots and the center of the DP. The vector length is in units of pixels and vector angle in units of radians.<jupyter_code>spotMaps_upd=pickle.load(open(saveDir+'/spotMaps.dat','rb'))
wins=np.array(spotMaps_upd['wins'])
roi=spotMaps_upd['roi']<jupyter_output><empty_output><jupyter_text>***`plot_ids`***: Select indices of EWPC peaks from which you would like to get the maps. If you selected only two, enter `[1,2]`**Optional**: For the following figures, adjust the ***`figureSize`*** if needed.<jupyter_code>peak_ids=[1,2]

plotSpotMaps(wins[np.array(peak_ids)-1,:],np.log(ewpc_img),trim_spotMaps(spotMaps_upd,peak_ids),figureSize=(8,4),plot_ids=peak_ids)<jupyter_output><empty_output><jupyter_text>CHECKPOINT 3PTO has a tetragonal unit cell with its long axis in the direction of ferroelectric polarization. The lattice mismatch with the DSO substrate causes the PTO film to form domains to accomodate the epitaxial strain.What are the stripes you see here? What does this tell you about the domain structure in PTO - think about where the substrate is and try to characterize the polarization direction in the different domains based on the distance of the EWPC spots from the center. [go to top](toc) 6. EWPC strain map<jupyter_code>new_spotMaps=trim_spotMaps(spotMaps_upd,peak_ids)
spotRef = makeRelativeSpotReference( new_spotMaps, [0,64,0,64]) #uses the average position of the EWPC spots in the ROI as reference for strain mapping
strainComponents = calculateStrainMap(new_spotMaps, spotRef,latticeCoords=1)<jupyter_output><empty_output><jupyter_text>The affine transform connecting the position of the EWPC peak position at each pixel to the reference position is calculated, and decoupled into in-plane rotation and lattice parameter change (strain tensor).If latticeCoords = 1, the strain map is calculated taking the direction of the two EWPC spots as the basis vectors. If set to 0, the usual x and y coordinates in Matlab are used as the basis vectors.The strain tensor components describes the lattice expansion/contraction in Cartesian coordinates: x- ($\epsilon_{11}$), y- ($\epsilon_{22}$) and diagonal xy-direction($\epsilon_{12}$) stretches. $\theta$ is the angle of in-plane lattice rotation.<jupyter_code>plotStrainTensor(strainComponents,figureSize=(8,8))<jupyter_output><empty_output><jupyter_text>Save ***`strainComponents`***<jupyter_code>pickle.dump(strainComponents,open(saveDir+'/strainComponents.dat','wb'))<jupyter_output><empty_output><jupyter_text>[go to top](toc) CHECKPOINT 4  Now try the strain mapping with a different set of EWPC spots. In the earlier case, we choose EWPC spots with a good SNR to form the strain maps.**1)** How does this compare with the earlier results? Explore using different EWPC spots and see how the results compare. This should help develop intuition about which EWPC spots to choose to get optimal results.**2)** What happens if you choose 2 EWPC spots that are the inverse of each other, i.e spots centered around (x,y) and (-x,-y) pixels? Try to explain why, thinking about how the EWPC pattern is mathematically defined. In general, the EWPC spots must be robust over the entire ROI, and should have a good SNR. Noisy spots $\implies$ Noisy strain maps.Also make sure, the window around the EWPC spot is such that it is: i) big enough to captures any shifts due to changes in lattice spacing ii) small enough that there is no interference from other nearby EWPC spots. [go to top](toc) 7. Optional: Segmentation using PCA decomposition and K-means clustering   This part is for advanced users and is optional. You can segment your data into your sample region and background using PCA decomposition and K-means clustering.<jupyter_code>n_components=20
radii=[5,64]

center_beam_mask=create_haadf_mask((ewpc.shape[2],ewpc.shape[3]),radii)
pca,scores=pca_decomposition(ewpc,n_components,center_beam_mask)
## save pca model and scores
pickle.dump(pca, open(saveDir+'pca.pkl','wb'))
np.save(saveDir+'/scores.npy',scores)
plot_false_color_img(pca,scores,center_beam_mask,cmap='gray')
clusters_range=[2,5]
cut_off=6

kmeans_dict=perform_kmeans(scores,cut_off,clusters_range)
plot_kmeans_dict(kmeans_dict)
pickle.dump(kmeans_dict,open(saveDir+'/kmeans.dat','wb'))
# kmeans_dict=pickle.load(open(saveDir+'/kmeans.dat','rb'))

### choose the number of clusters(3 looks good), and select the index of the region of interest

mask=kmeans_dict[3]['label']==0<jupyter_output><empty_output><jupyter_text>Use mask from pca and kmeans for spotmaps<jupyter_code>spotMaps=get_spotMaps(data4d,wins,valid=mask,tol=1e-4)<jupyter_output><empty_output><jupyter_text>you can go back and plot these spotmaps and calculate strain  compare your kmeans segmentation to the dark field threshold based segmentation (choose the threshold manually)<jupyter_code>thresh=2000
df_cep=calculate_DF(ewpc,wins)
segment_manually(df_cep,thresh=thresh)
valid=df_cep>thresh
spotMaps= get_spotMaps(data4d,wins,valid=valid,tol=1e-4,method='L-BFGS-B') 
plotSpotMaps(np.array(wins),np.log(ewpc_img),spotMaps,figureSize=(6,4))
new_spotMaps=trim_spotMaps(spotMaps,peak_ids)
spotRef = makeRelativeSpotReference( new_spotMaps, [0,64,0,64]) #uses the average position of the EWPC spots in the ROI as reference for strain mapping
strainComponents = calculateStrainMap(new_spotMaps, spotRef,latticeCoords=1)
plotStrainTensor(strainComponents,figureSize=(5,5))<jupyter_output><empty_output><jupyter_text>ETCYou may change the functions within ***`fun.py`***. If you do so, you can import the updated functions by running the following cell without restarting the kernel.<jupyter_code>import sys, importlib
importlib.reload(sys.modules['fun'])
from fun import *<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08537472714043172, 'fraction_numerical': 0.017341741450400194, 'mean_word_length': 6.0755040755040755, 'pattern_counts': {'":': 0, '<': 86, '<?xml version=': 0, '>': 87, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '354198', 'license': 'MIT', 'n_tokens_mistral': 5110, 'n_tokens_neox': 4717, 'n_words': 2043, 'path': 'Tutorial 5 - 4D STEM/T5-4DSTEM.ipynb', 'repo_name': 'paradimdata/Cornell_EM_SummerSchool_2021', 'size': '16497'}