Skip to content

Updating cluster module and ACOM for cluster based strain mapping #732

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
smribet merged 6 commits into from
Nov 3, 2025
Merged

Updating cluster module and ACOM for cluster based strain mapping #732

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
fc7fb17
Enhance similarity calculation with optional Gaussian smoothing and i…
serinlee1065 Sep 2, 2025
e3e918a
Replacing np.integer with np.int32
Sep 4, 2025
2ca5901
updating cluster module and ACOM to use cluster dataset for strain ma…
Sep 23, 2025
c21eca4
doc strings update and delete extra code
smribet Sep 29, 2025
da82837
Adding lines to check the r_space_mask condition
serinlee1065 Oct 28, 2025
e59820c
a few small changes
smribet Oct 31, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
129 changes: 66 additions & 63 deletions py4DSTEM/process/diffraction/crystal_ACOM.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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")
Expand All @@ -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, :]
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -2330,71 +2332,72 @@ 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(
Expand Down
Loading