Merge pull request #742 from serinlee1065/diffraction

From diffraction (fork) to dev for clustering module

Author: smribet

py4DSTEM/__init__.py

@@ -25,7 +25,6 @@
 from py4DSTEM import io
 from py4DSTEM.io import import_file, read, save
 
-
 ### basic data classes
 
 # data
@@ -40,7 +39,6 @@
 # datacube
 from py4DSTEM.datacube import DataCube, VirtualImage, VirtualDiffraction
 
-
 ### visualization
 
 from py4DSTEM import visualize
@@ -64,15 +62,12 @@
 # diffraction
 from py4DSTEM.process.diffraction import Crystal, Orientation
 
-
 # ptycho
 from py4DSTEM.process import phase
 
-
 # polar
 from py4DSTEM.process.polar import PolarDatacube
 
-
 # strain
 from py4DSTEM.process.strain.strain import StrainMap
 
@@ -89,7 +84,6 @@
 from py4DSTEM import preprocess
 from py4DSTEM import process
 
-
 ### utilities
 
 # config

py4DSTEM/braggvectors/diskdetection_parallel_new.py

@@ -21,7 +21,6 @@
 from emdfile import PointListArray, PointList
 from py4DSTEM.braggvectors.diskdetection import _find_Bragg_disks_single_DP_FK
 
-
 #### SERIALISERS ####
 # Define Serialiser
 # these are functions which allow the hdf5 objects to be passed. May not be required anymore

py4DSTEM/data/propagating_calibration.py

@@ -3,7 +3,6 @@
 
 import warnings
 
-
 # This is the abstract pattern:
 
 

py4DSTEM/datacube/virtualimage.py

@@ -13,7 +13,6 @@
 from py4DSTEM.preprocess import get_shifted_ar
 from py4DSTEM.visualize import show
 
-
 # Virtual image container class
 
 

py4DSTEM/io/__init__.py

@@ -3,6 +3,5 @@
 from py4DSTEM.io.read import read
 from py4DSTEM.io.save import save
 
-
 # google downloader
 from py4DSTEM.io.google_drive_downloader import gdrive_download, get_sample_file_ids

py4DSTEM/io/google_drive_downloader.py

@@ -2,7 +2,6 @@
 import os
 import warnings
 
-
 ### File IDs
 
 # single files

py4DSTEM/io/legacy/legacy13/v13_emd_classes/io.py

@@ -6,7 +6,6 @@
 from numbers import Number
 from emdfile import tqdmnd
 
-
 # Define the EMD group types
 
 EMD_group_types = {

py4DSTEM/io/legacy/legacy13/v13_py4dstem_classes/io.py

@@ -18,7 +18,6 @@
     _read_metadata,
 )
 
-
 # Calibration
 
 

py4DSTEM/io/legacy/legacy13/v13_to_14.py

@@ -3,7 +3,6 @@
 import numpy as np
 from emdfile import tqdmnd
 
-
 # v13 imports
 
 from py4DSTEM.io.legacy.legacy13.v13_emd_classes import (
@@ -25,7 +24,6 @@
     BraggVectors as BraggVectors13,
 )
 
-
 # v14 imports
 
 from emdfile import Root, Metadata, Array, PointList, PointListArray

py4DSTEM/process/diffraction/crystal_ACOM.py

@@ -335,7 +335,8 @@ def orientation_plan(
     )
     self.orientation_zone_axis_steps = (
         np.round(step / self.orientation_refine_ratio) * self.orientation_refine_ratio
-    ).astype(np.integer)
+    ).astype(np.int32)
+    # ).astype(np.integer)
 
     if self.orientation_fiber and self.orientation_fiber_angles[0] == 0:
         self.orientation_num_zones = int(1)
@@ -370,7 +371,8 @@ def orientation_plan(
             (self.orientation_zone_axis_steps + 1)
             * (self.orientation_zone_axis_steps + 2)
             / 2
-        ).astype(np.integer)
+        ).astype(np.int32)
+        # ).astype(np.integer)
         self.orientation_vecs = np.zeros((self.orientation_num_zones, 3))
         self.orientation_vecs[0, :] = self.orientation_zone_axis_range[0, :]
         self.orientation_inds = np.zeros((self.orientation_num_zones, 3), dtype="int")
@@ -379,7 +381,8 @@ def orientation_plan(
         # or circular arc SLERP for fiber texture
         for a0 in np.arange(1, self.orientation_zone_axis_steps + 1):
             inds = np.arange(a0 * (a0 + 1) / 2, a0 * (a0 + 1) / 2 + a0 + 1).astype(
-                np.integer
+                np.int32
+                # np.integer
             )
 
             p0 = pv[a0, :]
@@ -617,7 +620,8 @@ def orientation_plan(
 
     # Solve for number of angular steps along in-plane rotation direction
     self.orientation_in_plane_steps = np.round(360 / angle_step_in_plane).astype(
-        np.integer
+        np.int32
+        # np.integer
     )
 
     # Calculate -z angles (Euler angle 3)
@@ -2208,8 +2212,6 @@ def calculate_strain(
     deformation tensor which transforms the simulated diffraction pattern
     into the experimental pattern, for all probe positons.
 
-    TODO: add robust fitting?
-
     Parameters
     ----------
     bragg_peaks_array (PointListArray):
@@ -2330,71 +2332,75 @@ def calculate_strain(
                     inds_match[a0] = ind_min
                     keep[a0] = True
 
-            # Get all paired peaks
-            qxy = np.vstack((p.data["qx"][keep], p.data["qy"][keep])).T
-            qxy_ref = np.vstack(
-                (p_ref.data["qx"][inds_match[keep]], p_ref.data["qy"][inds_match[keep]])
-            ).T
+            if np.sum(keep) >= min_num_peaks:
+                # Get all paired peaks
+                qxy = np.vstack((p.data["qx"][keep], p.data["qy"][keep])).T
+                qxy_ref = np.vstack(
+                    (
+                        p_ref.data["qx"][inds_match[keep]],
+                        p_ref.data["qy"][inds_match[keep]],
+                    )
+                ).T
 
-            # Fit transformation matrix
-            # Note - not sure about transpose here
-            # (though it might not matter if rotation isn't included)
-            if intensity_weighting:
-                weights = np.sqrt(p.data["intensity"][keep, None]) * 0 + 1
-                m = lstsq(
-                    qxy_ref * weights,
-                    qxy * weights,
-                    rcond=None,
-                )[0].T
-            else:
-                m = lstsq(
-                    qxy_ref,
-                    qxy,
-                    rcond=None,
-                )[0].T
-
-            # Robust fitting
-            if robust:
-                for a0 in range(5):
-                    # calculate new weights
-                    qxy_fit = qxy_ref @ m
-                    diff2 = np.sum((qxy_fit - qxy) ** 2, axis=1)
-
-                    weights = np.exp(
-                        diff2 / ((-2 * robust_thresh**2) * np.median(diff2))
-                    )[:, None]
-                    if intensity_weighting:
-                        weights *= np.sqrt(p.data["intensity"][keep, None])
-
-                    # calculate new fits
+                # Fit transformation matrix
+                # Note - not sure about transpose here
+                # (though it might not matter if rotation isn't included)
+                if intensity_weighting:
+                    weights = np.sqrt(p.data["intensity"][keep, None]) * 0 + 1
                     m = lstsq(
                         qxy_ref * weights,
                         qxy * weights,
                         rcond=None,
                     )[0].T
+                else:
+                    m = lstsq(
+                        qxy_ref,
+                        qxy,
+                        rcond=None,
+                    )[0].T
 
-            # Set values into the infinitesimal strain matrix
-            strain_map.get_slice("e_xx").data[rx, ry] = 1 - m[0, 0]
-            strain_map.get_slice("e_yy").data[rx, ry] = 1 - m[1, 1]
-            strain_map.get_slice("e_xy").data[rx, ry] = -(m[0, 1] + m[1, 0]) / 2.0
-            strain_map.get_slice("theta").data[rx, ry] = (m[0, 1] - m[1, 0]) / 2.0
-
-            # Add finite rotation from ACOM orientation map.
-            # I am not sure about the relative signs here.
-            # Also, maybe I need to add in the mirror operator?
-            if orientation_map.mirror[rx, ry, 0]:
-                strain_map.get_slice("theta").data[rx, ry] += (
-                    orientation_map.angles[rx, ry, 0, 0]
-                    + orientation_map.angles[rx, ry, 0, 2]
-                )
-            else:
-                strain_map.get_slice("theta").data[rx, ry] -= (
-                    orientation_map.angles[rx, ry, 0, 0]
-                    + orientation_map.angles[rx, ry, 0, 2]
-                )
+                # Robust fitting
+                if robust:
+                    for a0 in range(5):
+                        # calculate new weights
+                        qxy_fit = qxy_ref @ m
+                        diff2 = np.sum((qxy_fit - qxy) ** 2, axis=1)
+
+                        weights = np.exp(
+                            diff2 / ((-2 * robust_thresh**2) * np.median(diff2))
+                        )[:, None]
+                        if intensity_weighting:
+                            weights *= np.sqrt(p.data["intensity"][keep, None])
+
+                        # calculate new fits
+                        m = lstsq(
+                            qxy_ref * weights,
+                            qxy * weights,
+                            rcond=None,
+                        )[0].T
+
+                # Set values into the infinitesimal strain matrix
+                strain_map.get_slice("e_xx").data[rx, ry] = 1 - m[0, 0]
+                strain_map.get_slice("e_yy").data[rx, ry] = 1 - m[1, 1]
+                strain_map.get_slice("e_xy").data[rx, ry] = -(m[0, 1] + m[1, 0]) / 2.0
+                strain_map.get_slice("theta").data[rx, ry] = (m[0, 1] - m[1, 0]) / 2.0
+
+                # Add finite rotation from ACOM orientation map.
+                # I am not sure about the relative signs here.
+                # Also, maybe I need to add in the mirror operator?
+                if orientation_map.mirror[rx, ry, 0]:
+                    strain_map.get_slice("theta").data[rx, ry] += (
+                        orientation_map.angles[rx, ry, 0, 0]
+                        + orientation_map.angles[rx, ry, 0, 2]
+                    )
+                else:
+                    strain_map.get_slice("theta").data[rx, ry] -= (
+                        orientation_map.angles[rx, ry, 0, 0]
+                        + orientation_map.angles[rx, ry, 0, 2]
+                    )
 
-        else:
-            strain_map.get_slice("mask").data[rx, ry] = 0.0
+            else:
+                strain_map.get_slice("mask").data[rx, ry] = 0.0
 
     if rotation_range is not None:
         strain_map.get_slice("theta").data[:] = np.mod(