cSAXS/cdtools
0
0
Fork 0
mirror of https://github.com/cdtools-developers/cdtools.git synced 2026-10-12 19:10:26 +02:00
Code Issues Packages Projects Releases Wiki Activity

Add line-based FRC function and corresponding tests, and example file

Browse Source 
- Implemented `line_based_frc` function for computing line-based Fourier Ring Correlation with multiple threshold methods.
- Added unit tests for `line_based_frc` to validate functionality with both PyTorch tensors and NumPy arrays.
- Created example script to demonstrate usage of `line_based_frc` with generated test images.
This commit is contained in:
gnzng
2025-07-10 14:24:24 -07:00
parent 603d47417b
commit 6de3491928
3 changed files with 356 additions and 2 deletions
Download Patch File Download Diff File Expand all files Collapse all files
examples/line_based_frc.py
+117
View File
@@ -0,0 +1,117 @@
import numpy as np
import matplotlib.pyplot as plt
from cdtools.tools.analysis import line_based_frc
def create_stripe_test_pattern(shape, stripe_spacing, noise_level=0.1):
"""Create test pattern with anisotropic stripes"""
h, w = shape
# Create vertical stripes (varying along x-axis)
x = np.arange(w)
pattern = np.sin(2 * np.pi * x / stripe_spacing)
# Broadcast to full image
image = np.tile(pattern, (h, 1))
# Add some modulation along y-axis for realism
y = np.arange(h)
y_mod = 1 + 0.3 * np.sin(2 * np.pi * y / (h // 3))
image = image * y_mod[:, np.newaxis]
# make the lower part 0
image[7 * h // 8:] = 0
# Add noise
image += noise_level * np.random.randn(h, w)
return image
# Create test images
shape = (200, 300)
stripe_spacing = 8 # pixels
unit = "nm"
pixel_size = 19.0 # nm, for example
# Create two slightly different versions to simulate repeated measurements
np.random.seed(42)
image1 = create_stripe_test_pattern(
shape,
stripe_spacing,
noise_level=0.05,
)
np.random.seed(43) # Different seed for second image
image2 = create_stripe_test_pattern(
shape,
stripe_spacing,
noise_level=0.05,
)
# Display test images
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5))
ax1.imshow(image1, cmap="gray")
ax1.set_title("Test Image 1")
ax1.axis("off")
ax2.imshow(image2, cmap="gray")
ax2.set_title("Test Image 2")
ax2.axis("off")
plt.tight_layout()
frc_thresholds = {"1/2": 0.5, "1/4": 0.25, "1/7": 1 / 7}
# Compute binned FRC
frequencies, frc_curve, resolution_dict = line_based_frc(
image1, image2, axis=1, n_bins=100, thresholds=frc_thresholds,
pixel_size=pixel_size, unit=unit
)
###########################
# Plot FRC curve with threshold crossings
###########################
# Color map for different thresholds
colors = ["red", "blue", "green", "orange", "purple", "brown", "pink", "gray"]
fig, ax1 = plt.subplots(1, 1, figsize=(6, 6))
# FRC curve with all thresholds
ax1.plot(frequencies, frc_curve, "k-", linewidth=3, label="FRC", zorder=10)
for i, (method, threshold_val) in enumerate(frc_thresholds.items()):
color = colors[i % len(colors)]
ax1.axhline(
y=threshold_val,
color=color,
linestyle="--",
alpha=0.8,
label=f"{method} ({resolution_dict[method]:.3f} {unit})",
)
# Mark threshold frequency using the resolution value
resolution = resolution_dict[method]
if resolution is not None and resolution > 0:
freq = pixel_size / resolution
# Find the closest frequency in the array
ax1.scatter(
# threshold_val * np.ones_like(freq),
freq,
threshold_val,
color=color,
s=80,
edgecolor="black",
zorder=20,
)
ax1.set_xlabel(f"Spatial Frequency (1/{unit})")
ax1.set_ylabel("FRC")
ax1.set_title("Line-based FRC Curve")
ax1.legend()
plt.tight_layout()
plt.show()
src/cdtools/tools/analysis/analysis.py
+145 -1
View File
@@ -34,6 +34,7 @@ __all__ = [
'standardize_reconstruction_set',
'standardize_reconstruction_pair',
'calc_spectral_info',
'line_based_frc',
]
@@ -1525,4 +1526,147 @@ def calc_spectral_info(dataset, nbins=50):
pattern_snr = sum_spectrum_sq / sum_spectrum
return sum_pattern, frc_bins[:-1], mean_spectrum
def line_based_frc(image1, image2, axis=1, n_bins=None, thresholds={}, pixel_size=1.0, unit="pixels"):
"""
Compute line-based FRC with multiple threshold methods
Parameters
----------
image1, image2: np.ndarray or torch.Tensor
Input images (should be same size). Must be 2D arrays.
axis: int
axis along which to extract lines (1 for horizontal lines, 0 for vertical)
n_bins: int
number of frequency bins (if None, no binning is applied)
thresholds: dict
Dictionary of threshold methods and values
pixel_size: float
Physical size of one pixel
unit: str
Unit for the pixel size and resolution output
Returns
-------
frequencies: np.ndarray
Spatial frequency array
frc_curve: np.ndarray
FRC values as function of frequency
resolution_dict: dict
Resolution at different thresholds
"""
# transform all torch tensors to np.arrays
if isinstance(image1, t.Tensor):
image1 = image1.numpy()
if isinstance(image2, t.Tensor):
image2 = image2.numpy()
assert image1.shape == image2.shape, "Images must have same shape"
if image1.ndim != 2:
raise ValueError("Images must be 2D")
if image2.ndim != 2:
raise ValueError("Images must be 2D")
if axis == 1:
# Extract horizontal lines (perpendicular to vertical stripes)
lines1 = image1
lines2 = image2
n_freq = image1.shape[1] // 2
else:
# Extract vertical lines (perpendicular to horizontal stripes)
lines1 = image1.T
lines2 = image2.T
n_freq = image1.shape[0] // 2
# Compute 1D FFT for each line
fft1 = np.fft.fft(lines1, axis=1)
fft2 = np.fft.fft(lines2, axis=1)
# Only keep positive frequencies
fft1 = fft1[:, :n_freq]
fft2 = fft2[:, :n_freq]
# Get the raw frequency array
raw_frequencies = np.fft.fftfreq(lines1.shape[1])[:n_freq]
if n_bins is None:
# No binning - use original approach
numerator = np.sum(fft1 * np.conj(fft2), axis=0)
denominator = np.sqrt(np.sum(np.abs(fft1)**2, axis=0) * np.sum(np.abs(fft2)**2, axis=0))
frc_curve = np.abs(numerator) / denominator
frequencies = raw_frequencies
else:
# Binning approach - similar to 2D FRC
max_freq = raw_frequencies[-1]
freq_bins = np.linspace(0, max_freq, n_bins + 1)
# Compute binned FRC
frc_curve = []
frequencies = []
for i in range(len(freq_bins) - 1):
# Find frequency indices in this bin
mask = (raw_frequencies >= freq_bins[i]) & (raw_frequencies < freq_bins[i+1])
n_freq_components = np.sum(mask)
if n_freq_components > 0:
# Extract FFT values in this frequency bin
fft1_bin = fft1[:, mask] # Shape: (n_lines, n_freq_components)
fft2_bin = fft2[:, mask]
# Compute FRC for this bin (sum over both lines and frequencies)
numerator = np.sum(fft1_bin * np.conj(fft2_bin))
denominator = np.sqrt(np.sum(np.abs(fft1_bin)**2) * np.sum(np.abs(fft2_bin)**2))
frc_value = np.abs(numerator) / denominator
frc_curve.append(frc_value)
frequencies.append((freq_bins[i] + freq_bins[i+1]) / 2)
# Convert lists to numpy arrays
frc_curve = np.array(frc_curve)
frequencies = np.array(frequencies)
# Find resolution at each threshold using interpolation
resolution_dict = {}
for method, threshold_val in thresholds.items():
# Find where FRC crosses the threshold
above_threshold = frc_curve >= threshold_val
# Find all crossings (from above to below threshold)
crossings = np.where(np.diff(above_threshold.astype(int)) == -1)[0]
if crossings.size > 0:
# Use the first crossing (lowest frequency where FRC drops below threshold)
idx = crossings[0]
freq1 = frequencies[idx]
freq2 = frequencies[idx + 1]
frc1 = frc_curve[idx]
frc2 = frc_curve[idx + 1]
# Linear interpolation: find frequency where FRC == threshold_val
resolution_frequency = freq1 + (threshold_val - frc1) * (freq2 - freq1) / (frc2 - frc1)
else:
# No crossing found, check if all above or all below
if np.all(above_threshold):
# FRC never drops below threshold, use max frequency
resolution_frequency = frequencies[-1]
else:
# FRC never reaches threshold
resolution_frequency = 0
if resolution_frequency > 0:
resolution_pixels = 1 / resolution_frequency
resolution_physical = resolution_pixels * pixel_size
else:
resolution_physical = np.inf
resolution_dict[method] = resolution_physical
return frequencies, frc_curve, resolution_dict
tests/tools/test_analysis.py
+94 -1
View File
@@ -2,9 +2,9 @@ import numpy as np
from scipy import linalg as la
from scipy.sparse import linalg as spla
import torch as t
from itertools import combinations
from cdtools.tools import analysis, initializers
from cdtools.tools.analysis import line_based_frc
def test_product_svd():
@@ -519,3 +519,96 @@ def test_calc_generalized_rms_error():
assert (analysis.calc_generalized_rms_error(fields_1, fields_2, dims=1).shape == t.Size([3]))
def test_line_based_frc():
# First test some basic input with two fields of the same size
field_1 = t.rand(30, 17, dtype=t.complex128)
field_2 = t.rand(30, 17, dtype=t.complex128)
# Compute FRC
freq, frc, res_dict = line_based_frc(field_1, field_2,
axis=1, n_bins=50,
thresholds={"1/7": 1/7},
)
# check if length of frc is longer than 0
assert len(frc) > 0
# check if freq is same length as frc
assert len(freq) == len(frc)
# check if res_dict is a dictionary
assert isinstance(res_dict, dict)
# Now check if it works with numpy arrays
field_1_np = field_1.numpy()
field_2_np = field_2.numpy()
freq_np, frc_np, res_dict_np = line_based_frc(field_1_np, field_2_np)
# Do the same checks as above
assert len(frc_np) > 0
assert len(freq_np) == len(frc_np)
assert isinstance(res_dict_np, dict)
# Create test images from the provided example
def create_stripe_test_pattern(shape, stripe_spacing, noise_level=0.1):
"""Create test pattern with anisotropic stripes"""
h, w = shape
# Create vertical stripes (varying along x-axis)
x = np.arange(w)
pattern = np.sin(2 * np.pi * x / stripe_spacing)
# Broadcast to full image
image = np.tile(pattern, (h, 1))
# Add some modulation along y-axis for realism
y = np.arange(h)
y_mod = 1 + 0.3 * np.sin(2 * np.pi * y / (h // 3))
image = image * y_mod[:, np.newaxis]
# make the lower part 0
image[7 * h // 8:] = 0
# Add noise
image += noise_level * np.random.randn(h, w)
return image
shape = (200, 300)
stripe_spacing = 8 # pixels
unit = "nm"
pixel_size = 19.0 # nm, for example
# Create two slightly different versions to simulate repeated measurements
np.random.seed(42)
image1 = create_stripe_test_pattern(
shape,
stripe_spacing,
noise_level=0.05,
)
np.random.seed(43) # Different seed for second image
image2 = create_stripe_test_pattern(
shape,
stripe_spacing,
noise_level=0.05,
)
frc_thresholds = {"1/2": 0.5, "1/4": 0.25, "1/7": 1 / 7}
# Compute binned FRC
frequencies, frc_curve, res_dict = line_based_frc(
image1, image2, axis=1, n_bins=100, thresholds=frc_thresholds,
pixel_size=pixel_size, unit=unit
)
# run through the above checks
assert len(frc) > 0
assert len(freq) == len(frc)
assert isinstance(res_dict, dict)
# check if the values are within expected ranges
assert res_dict["1/2"] > 70.0 and res_dict["1/2"] < 75.0
assert res_dict["1/4"] > 52.0 and res_dict["1/4"] < 57.0
assert res_dict["1/7"] > 45.0 and res_dict["1/7"] < 48.0