diff --git a/examples/line_based_frc.py b/examples/line_based_frc.py new file mode 100644 index 0000000..980a7f4 --- /dev/null +++ b/examples/line_based_frc.py @@ -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() diff --git a/src/cdtools/tools/analysis/analysis.py b/src/cdtools/tools/analysis/analysis.py index 9844edb..38b1306 100644 --- a/src/cdtools/tools/analysis/analysis.py +++ b/src/cdtools/tools/analysis/analysis.py @@ -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 diff --git a/tests/tools/test_analysis.py b/tests/tools/test_analysis.py index f0a1cf1..265a5ea 100644 --- a/tests/tools/test_analysis.py +++ b/tests/tools/test_analysis.py @@ -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