diff --git a/packages/scratch-core/src/conversion/preprocess_impression/preprocess_impression.py b/packages/scratch-core/src/conversion/preprocess_impression/preprocess_impression.py index 7272661f..247e4428 100644 --- a/packages/scratch-core/src/conversion/preprocess_impression/preprocess_impression.py +++ b/packages/scratch-core/src/conversion/preprocess_impression/preprocess_impression.py @@ -22,7 +22,7 @@ ) from conversion.preprocess_impression.tilt import apply_tilt_correction from conversion.preprocess_impression.utils import update_mark_data, Point2D -from conversion.resample import get_scaling_factors, resample_image_array +from conversion.resample import get_scaling_factors, resample_array_2d def preprocess_impression_mark( @@ -77,7 +77,7 @@ def preprocess_impression_mark( factors = get_scaling_factors( scales=original_scales, target_scale=params.pixel_size ) - fitted_surface = resample_image_array(fitted_surface, factors=factors) + fitted_surface = resample_array_2d(fitted_surface, factors=factors) # Stage 7: High-pass filter if params.highpass_cutoff is not None: diff --git a/packages/scratch-core/src/conversion/profile_correlator/__init__.py b/packages/scratch-core/src/conversion/profile_correlator/__init__.py new file mode 100644 index 00000000..bf49ad63 --- /dev/null +++ b/packages/scratch-core/src/conversion/profile_correlator/__init__.py @@ -0,0 +1,66 @@ +""" +Profile correlator module for striated mark comparison. + +This module provides functions for comparing 1D profiles of striated marks +(scratch marks, toolmarks, etc.) using global brute-force search and +correlation analysis. + +The main entry point is :func:`correlate_profiles`, which handles the complete +comparison workflow including: + +- Sampling distance equalization +- Global brute-force search over all shift positions and scale factors +- Selection of alignment with maximum cross-correlation +- Computation of comparison metrics + +Submodules +---------- +- data_types: Core data structures (Profile, AlignmentParameters, StriationComparisonResults) +- correlator: Main entry point function (global brute-force search) +- transforms: Resampling operations for pixel scale equalization +- statistics: Statistical metrics (correlation, roughness, overlap ratio) +""" + +# Core data types +from conversion.profile_correlator.data_types import ( + AlignmentParameters, + StriationComparisonResults, + Profile, +) + +# Main entry point +from conversion.profile_correlator.correlator import correlate_profiles + +# Transform functions +from conversion.profile_correlator.transforms import equalize_pixel_scale + +# Statistics functions +from conversion.profile_correlator.statistics import ( + compute_cross_correlation, + compute_overlap_ratio, + compute_roughness_sa, + compute_roughness_sq, + compute_normalized_square_based_roughness_differences, +) + +# Re-export cutoff_to_gaussian_sigma from conversion.filter for convenience +from conversion.filter.gaussian import cutoff_to_gaussian_sigma + +__all__ = [ + # Main entry point + "correlate_profiles", + # Data types + "Profile", + "AlignmentParameters", + "StriationComparisonResults", + # Transforms + "equalize_pixel_scale", + # Statistics + "compute_cross_correlation", + "compute_overlap_ratio", + "compute_roughness_sa", + "compute_roughness_sq", + "compute_normalized_square_based_roughness_differences", + # Filtering (re-exported from conversion.filter) + "cutoff_to_gaussian_sigma", +] diff --git a/packages/scratch-core/src/conversion/profile_correlator/correlator.py b/packages/scratch-core/src/conversion/profile_correlator/correlator.py new file mode 100644 index 00000000..e341a572 --- /dev/null +++ b/packages/scratch-core/src/conversion/profile_correlator/correlator.py @@ -0,0 +1,274 @@ +""" +Main entry point for profile correlation. + +This module provides the primary interface for comparing striated mark profiles. +It uses a brute-force approach: try all possible shifts and scales between +profiles and select the one with maximum cross-correlation, subject to a +minimum overlap constraint. + +All length and height measurements are in meters (SI units). + +Algorithm Overview +------------------ +The algorithm uses a global brute-force search strategy: + +1. **Equalizes pixel scales** between profiles (downsample to common resolution) +2. **Tries multiple scale factors** (e.g., 0.95, 0.97, ..., 1.03, 1.05) +3. **For each scale, tries all shifts** with at least min_overlap_distance overlap +4. **Computes correlation** at each shift position +5. **Selects maximum** correlation as the optimal alignment + +This approach is guaranteed to find the global maximum correlation, which may +be at a position far from zero shift for repetitive patterns. +""" + +import numpy as np + +from container_models.base import FloatArray1D +from conversion.profile_correlator.data_types import ( + AlignmentParameters, + AlignmentResult, + StriationComparisonResults, + Profile, + RoughnessMetrics, +) +from conversion.profile_correlator.transforms import equalize_pixel_scale +from conversion.profile_correlator.statistics import ( + compute_cross_correlation, + compute_overlap_ratio, + compute_roughness_sa, + compute_roughness_sq, + compute_normalized_square_based_roughness_differences, +) +from conversion.resample import resample_array_1d + + +def correlate_profiles( + profile_ref: Profile, + profile_comp: Profile, + params: AlignmentParameters = AlignmentParameters(), +) -> StriationComparisonResults | None: + """ + Compare two striated mark profiles and compute similarity metrics. + + This function performs profile alignment using a brute-force approach: + it tries multiple scale factors and all possible shifts between the profiles + (subject to a minimum overlap constraint) and selects the combination with + maximum cross-correlation. + + The workflow is: + + 1. **Equalize sampling distances**: If the two profiles have different pixel + sizes, the higher-resolution one is downsampled to match the lower resolution. + + 2. **Brute-force search**: Try all shift/scale combinations and compute + cross-correlation at each position. For each scale factor in the range + defined by params.max_scaling, determine valid shifts that maintain at + least min_overlap_distance of overlap, and select the shift/scale with + maximum correlation. + + 3. **Compute metrics**: Calculate comprehensive comparison metrics including + correlation coefficient, roughness parameters, and signature differences. + + All measurements are in meters (SI units). + + :param profile_ref: Reference profile to compare against. + :param profile_comp: Compared profile to align to the reference. + :param params: Alignment parameters. Key parameters: + - max_scaling: Maximum scaling deviation (e.g., 0.05 for ±5%) + - min_overlap_distance: Minimum overlap distance in meters (default 200 μm) + :returns: StriationComparisonResults containing all computed metrics, or None if no + valid alignment could be found. + """ + # Step 1: Equalize pixel scales + profile_ref_eq, profile_comp_eq = equalize_pixel_scale(profile_ref, profile_comp) + pixel_size = profile_ref_eq.pixel_size + + heights_ref = profile_ref_eq.heights + heights_comp = profile_comp_eq.heights + len_ref = len(heights_ref) + len_comp = len(heights_comp) + + # Minimum overlap in samples + min_overlap_samples = int(np.ceil(params.min_overlap_distance / pixel_size)) + + # Early exit if either profile is too short + if len_ref < min_overlap_samples or len_comp < min_overlap_samples: + return None + + # Generate scale factors to try (7 steps gives ~1.7% intervals for ±5%) + scale_factors = np.linspace( + 1.0 - params.max_scaling, 1.0 + params.max_scaling, params.n_scale_steps + ) + + # make scaling symmetric to what you choose as ref or comp + scale_factors = np.unique(np.concatenate((scale_factors, 1 / scale_factors))) + + # Step 2: Find best alignment + alignment = _find_best_alignment( + heights_ref, heights_comp, scale_factors, min_overlap_samples + ) + if alignment is None: + return None + + # Step 3: Compute and return metrics + return _compute_metrics(alignment, pixel_size, len_ref, len_comp) + + +def _calculate_idx_parameters( + shift: int, len_small: int, len_large: int +) -> tuple[int, int, int]: + """ + Find starting idx for both striations, and compute overlap length + """ + + if shift >= 0: + idx_large_start = shift + idx_small_start = 0 + overlap_length = min(len_large - shift, len_small) + else: + idx_large_start = 0 + idx_small_start = -shift + overlap_length = min(len_large, len_small + shift) + + return idx_small_start, idx_large_start, overlap_length + + +def _find_best_alignment( + heights_ref: FloatArray1D, + heights_comp: FloatArray1D, + scale_factors: FloatArray1D, + min_overlap_samples: int, +) -> AlignmentResult | None: + """ + Find the best alignment between two profiles using brute-force search. + + Tries all combinations of scale factors and shifts, returning the one + with maximum correlation. + + :param heights_ref: Reference profile heights. + :param heights_comp: Comparison profile heights. + :param scale_factors: Array of scale factors to try. + :param min_overlap_samples: Minimum required overlap in samples. + :returns: Best alignment result, or None if no valid alignment found. + """ + len_ref = len(heights_ref) + + best_correlation = -np.inf + best_shift = None + best_scale = None + + for scale in scale_factors: + heights_comp_scaled = resample_array_1d(heights_comp, scale) + len_comp = len(heights_comp_scaled) + # Calculate shift range (ensure minimum overlap) + min_shift = -(len_comp - min_overlap_samples) + max_shift = len_ref - min_overlap_samples + + for shift in range(min_shift, max_shift + 1): + idx_comp_start, idx_ref_start, overlap_length = _calculate_idx_parameters( + shift, len_comp, len_ref + ) # Calculate overlap region for this shift + + if overlap_length < min_overlap_samples: + continue + + partial_ref = heights_ref[idx_ref_start : idx_ref_start + overlap_length] + partial_comp = heights_comp_scaled[ + idx_comp_start : idx_comp_start + overlap_length + ] + + correlation = compute_cross_correlation(partial_ref, partial_comp) + + if correlation and correlation > best_correlation: + best_correlation = correlation + best_shift = shift + best_scale = scale + + if best_correlation is None: + return None + + # Redo computations for best_cale and best_shift (instead of copying partial_ref and partial_comp above multiple times. This saves time.) + heights_comp_scaled = resample_array_1d(heights_comp, best_scale) + len_comp = len(heights_comp_scaled) + idx_comp_start, idx_ref_start, overlap_length = _calculate_idx_parameters( + best_shift, len_comp, len_ref + ) + + best_ref_overlap = heights_ref[idx_ref_start : idx_ref_start + overlap_length] + best_comp_overlap = heights_comp_scaled[ + idx_comp_start : idx_comp_start + overlap_length + ] + + return AlignmentResult( + correlation=best_correlation, + shift=best_shift, + scale=best_scale, + ref_overlap=best_ref_overlap, + comp_overlap=best_comp_overlap, + ) + + +def _compute_metrics( + alignment: AlignmentResult, + pixel_size: float, + len_ref: int, + len_comp: int, +) -> StriationComparisonResults: + """ + Compute comparison metrics from an alignment result. + + :param alignment: The best alignment found. + :param pixel_size: Pixel size in meters. + :param len_ref: Length of reference profile in samples. + :param len_comp: Length of comparison profile in samples. + :returns: Full comparison results. + """ + ref_overlap = alignment.ref_overlap + comp_overlap = alignment.comp_overlap + + # Convert to meters + position_shift = alignment.shift * pixel_size + overlap_length = len(ref_overlap) * pixel_size + ref_length = len_ref * pixel_size + comp_length = len_comp * pixel_size + + overlap_ratio = compute_overlap_ratio(overlap_length, ref_length, comp_length) + + # Roughness metrics + sa_ref = compute_roughness_sa(ref_overlap) + mean_square_ref = compute_roughness_sq(ref_overlap) + sa_comp = compute_roughness_sa(comp_overlap) + mean_square_comp = compute_roughness_sq(comp_overlap) + + # Difference profile roughness + diff_profile = comp_overlap - ref_overlap + sa_diff = compute_roughness_sa(diff_profile) + mean_square_of_difference = compute_roughness_sq(diff_profile) + + # Signature differences + roughness = RoughnessMetrics( + mean_square_ref=mean_square_ref, + mean_square_comp=mean_square_comp, + mean_square_of_difference=mean_square_of_difference, + ) + signature_diff = compute_normalized_square_based_roughness_differences(roughness) + + return StriationComparisonResults( + pixel_size=pixel_size, + position_shift=position_shift, + scale_factor=alignment.scale, + similarity_value=alignment.correlation, + overlap_length=overlap_length, + overlap_ratio=overlap_ratio, + correlation_coefficient=alignment.correlation, + sa_ref=sa_ref, + mean_square_ref=mean_square_ref, + sa_comp=sa_comp, + mean_square_comp=mean_square_comp, + sa_diff=sa_diff, + mean_square_of_difference=mean_square_of_difference, + ds_roughness_normalized_to_reference=signature_diff.roughness_normalized_to_reference, + ds_roughness_normalized_to_compared=signature_diff.roughness_normalized_to_compared, + ds_roughness_normalized_to_reference_and_compared=signature_diff.roughness_normalized_to_reference_and_compared, + ) diff --git a/packages/scratch-core/src/conversion/profile_correlator/data_types.py b/packages/scratch-core/src/conversion/profile_correlator/data_types.py new file mode 100644 index 00000000..ded2d7e7 --- /dev/null +++ b/packages/scratch-core/src/conversion/profile_correlator/data_types.py @@ -0,0 +1,177 @@ +""" +Data types for profile correlation. + +This module defines the core data structures used for profile correlation analysis +of striated marks. It follows the patterns established in preprocess_impression. + +The main types are: +- Profile: Container for 1D profile data with metadata +- AlignmentParameters: Configuration for the brute-force alignment algorithm +- StriationComparisonResults: Full comparison metrics for striated mark analysis + +All length and height measurements are in meters (SI units). +""" + +from dataclasses import dataclass + +from container_models.base import FloatArray1D + + +@dataclass(frozen=True) +class Profile: + """ + Container for a 1D profile with metadata. + + Profiles represent 1D height measurements along a scratch/striation mark. + All measurements are in meters (SI units). + + :param heights: Height values as a 1D array of shape (N,). + :param pixel_size: Physical distance between samples in meters. + """ + + heights: FloatArray1D + pixel_size: float + + @property + def length(self) -> int: + """ + Get the number of samples in the profile. + + :returns: Number of samples in heights. + """ + return len(self.heights) + + +@dataclass(frozen=True) +class AlignmentParameters: + """ + Configuration parameters for profile alignment. + + This dataclass contains parameters for the global brute-force alignment + algorithm. The algorithm tries all possible shift positions and scale + factors, selecting the combination with maximum cross-correlation. + + All length parameters are in meters (SI units). + + :param max_scaling: Maximum allowed scaling deviation as a fraction. + E.g., 0.05 means scaling can vary from 0.95 to 1.05. + :param n_scale_steps: Number of steps along each scale factor. The scaling procedure is made symmetric resulting in (2n - 1) different scaling elements. + :param min_overlap_distance: Minimum required overlap between profiles + in meters. Alignments with less overlap are rejected. + """ + + max_scaling: float = 0.05 + n_scale_steps: int = 7 + min_overlap_distance: float = 350e-6 # 350 μm + + +@dataclass(frozen=True) +class RoughnessMetrics: + """ + Container for square-based roughness metrics of a profile pair. + + Uses ISO 25178 naming conventions: + - Sq: Quadratic mean roughness ('S' = surface, 'q' = quadratic mean) + + :param mean_square_ref: mean square 'roughness' of the reference profile. + :param mean_square_comp: mean square 'roughness' of the comparison profile. + :param mean_square_of_difference: mean square 'roughness' of the difference profile. + """ + + mean_square_ref: float + mean_square_comp: float + mean_square_of_difference: float + + +@dataclass(frozen=True) +class NormalizedSquareBasedRoughnessDifferences: + """ + Container for normalized square-based roughness difference metrics. + + These metrics quantify the difference between profiles normalized by + their roughness, providing dimensionless measures of dissimilarity. + + :param roughness_normalized_to_reference: square-based roughness difference normalized to reference, + computed as (mean_square_of_difference / mean_square_ref)^2. + :param roughness_normalized_to_compared: square-based roughness difference normalized to comparison, + computed as (mean_square_of_difference / mean_square_comp)^2. + :param roughness_normalized_to_reference_and_compared: square-based roughness difference using geometric mean + normalization, computed as mean_square_of_difference^2 / (mean_square_ref * mean_square_comp). + """ + + roughness_normalized_to_reference: float + roughness_normalized_to_compared: float + roughness_normalized_to_reference_and_compared: float + + +@dataclass(frozen=True) +class AlignmentResult: + """Result from the alignment search.""" + + correlation: float + shift: int + scale: float + ref_overlap: FloatArray1D + comp_overlap: FloatArray1D + + +@dataclass(frozen=True) +class StriationComparisonResults: + """ + Full comparison metrics for striated mark analysis. + + This immutable dataclass contains all metrics computed during profile comparison, + including registration parameters, roughness measurements, and + signature differences. + + All length and height measurements are in meters (SI units). + + Roughness parameters use ISO 25178 naming conventions: + - Sa: Arithmetic mean roughness ('S' = surface, 'a' = arithmetical mean) + - Sq: Root-mean-square roughness ('S' = surface, 'q' = quadratic mean) + + :param pixel_size: Pixel separation of profile in meters. + :param position_shift: Registration shift of compared profile relative to + reference in meters. + :param scale_factor: Registration scale factor applied to compared profile + (1.0 means no scaling). + :param similarity_value: Optimized value of the similarity metric used + during registration. + :param overlap_length: Length of the overlapping region after registration + in meters. + :param overlap_ratio: Ratio of overlap length to the length of the shorter + profile. + :param correlation_coefficient: Pearson cross-correlation coefficient + between the aligned profiles. + :param sa_ref: Arithmetic mean roughness (Sa) of the reference profile in meters. + :param mean_square_ref: Root-mean-square roughness (Sq) of the reference profile in meters. + :param sa_comp: Arithmetic mean roughness (Sa) of the compared profile in meters. + :param mean_square_comp: Root-mean-square roughness (Sq) of the compared profile in meters. + :param sa_diff: Arithmetic mean roughness (Sa) of the difference profile + (comparison minus reference) in meters. + :param mean_square_of_difference: Root-mean-square roughness (Sq) of the difference profile + (comparison minus reference) in meters. + :param ds_roughness_normalized_to_reference: Signature difference normalized by reference Sq, + computed as (mean_square_of_difference / mean_square_ref)^2. + :param ds_roughness_normalized_to_compared: Signature difference normalized by compared Sq, + computed as (mean_square_of_difference / mean_square_comp)^2. + :param ds_roughness_normalized_to_reference_and_compared: Combined signature difference using geometric mean + normalization, computed as mean_square_of_difference^2 / (mean_square_ref * mean_square_comp). + """ + + pixel_size: float + position_shift: float + scale_factor: float + similarity_value: float + overlap_length: float + overlap_ratio: float + correlation_coefficient: float + sa_ref: float + mean_square_ref: float + sa_comp: float + mean_square_comp: float + sa_diff: float + mean_square_of_difference: float + ds_roughness_normalized_to_reference: float + ds_roughness_normalized_to_compared: float + ds_roughness_normalized_to_reference_and_compared: float diff --git a/packages/scratch-core/src/conversion/profile_correlator/statistics.py b/packages/scratch-core/src/conversion/profile_correlator/statistics.py new file mode 100644 index 00000000..66871a5c --- /dev/null +++ b/packages/scratch-core/src/conversion/profile_correlator/statistics.py @@ -0,0 +1,152 @@ +""" +Statistical metrics for profile comparison. + +This module provides functions for computing statistical metrics between +1D profiles, including: + +- compute_cross_correlation: NaN-aware normalized cross-correlation +- compute_roughness_sa: Arithmetic mean roughness +- compute_roughness_sq: Root mean square roughness +- compute_overlap_ratio: Overlap ratio relative to shorter profile +- compute_normalized_square_based_roughness_differences: Normalized signature difference metrics + +All length and height measurements are in meters (SI units). +""" + +import numpy as np + +from container_models.base import FloatArray1D +from conversion.profile_correlator.data_types import ( + RoughnessMetrics, + NormalizedSquareBasedRoughnessDifferences, +) + + +def compute_cross_correlation( + profile_1: FloatArray1D, + profile_2: FloatArray1D, +) -> float | None: + """ + Compute cross-correlation between two profiles. + + This function computes the Pearson correlation coefficient between two + 1D profiles, properly handling NaN values by excluding them from the + calculation. + + :param profile_1: First profile as a 1D array. May contain NaN values. + :param profile_2: Second profile as a 1D array. Must have the same length + as profile_1. May contain NaN values. + :returns: Correlation coefficient in the range [-1, 1]. Returns NaN if + there are fewer than 2 valid (non-NaN) overlapping samples or if + either profile has zero variance. + :raises ValueError: If profiles have different lengths. + """ + profile_1 = profile_1.ravel() + profile_2 = profile_2.ravel() + + if len(profile_1) != len(profile_2): + raise ValueError( + f"Profiles must have the same length. " + f"Got {len(profile_1)} and {len(profile_2)}." + ) + + valid_mask = ~(np.isnan(profile_1) | np.isnan(profile_2)) + + if np.sum(valid_mask) < 2: + return None + + return float(np.corrcoef(profile_1[valid_mask], profile_2[valid_mask])[0, 1]) + + +def compute_roughness_sa(profile: FloatArray1D) -> float: + """ + Compute arithmetic mean roughness (ISO 25178 Sa parameter) of a profile. + + Sa is the arithmetic mean of the absolute values of the profile heights, + calculated as: mean(|z|). The 'S' denotes a surface/areal parameter and + 'a' denotes arithmetical mean. + + :param profile: 1D profile array. May contain NaN values which are ignored. + :returns: Arithmetic mean roughness (Sa) in the same units as the input profile. + """ + return float(np.nanmean(np.abs(profile - np.nanmean(profile)))) + + +def compute_roughness_sq(profile: FloatArray1D) -> float: + """ + Compute root-mean-square roughness (ISO 25178 Sq parameter) of a profile. + + Sq is the root-mean-square of the profile heights, calculated as: + sqrt(mean(z^2)). The 'S' denotes a surface/areal parameter and 'q' + denotes quadratic mean (root-mean-square). + + :param profile: 1D profile array. May contain NaN values which are ignored. + :returns: Root-mean-square roughness (Sq) in the same units as the input profile. + """ + return float(np.sqrt(np.nanmean((profile - np.nanmean(profile)) ** 2))) + + +def compute_overlap_ratio( + overlap_length: float, + ref_length: float, + comp_length: float, +) -> float: + """ + Compute mean overlap ratio relative to the two profile lengths. + + The overlap ratio indicates what fraction of the profiles are + covered by the overlap region after alignment. + + :param overlap_length: Length of the overlap region in meters. + :param ref_length: Length of the reference profile in meters. + :param comp_length: Length of the comparison profile in meters. + :returns: Overlap ratio in range [0, 1]. Returns NaN if shorter length is 0 + or if overlap_length exceeds shorter_length (invalid input). + """ + shorter_length = min(ref_length, comp_length) + if np.isclose(shorter_length, 0.0): + return np.nan + if overlap_length > shorter_length: + return np.nan + return 0.5 * (overlap_length / ref_length + overlap_length / comp_length) + + +def compute_normalized_square_based_roughness_differences( + roughness: RoughnessMetrics, +) -> NormalizedSquareBasedRoughnessDifferences: + """ + Compute normalized signature difference metrics. + + These metrics quantify the difference between profiles normalized by + their roughness, providing dimensionless measures of dissimilarity. + + :param roughness: Container with quadratic mean roughness (Sq) values for + the reference profile, comparison profile, and difference profile. + :returns: NormalizedSquareBasedRoughnessDifferences containing normalized metrics. + Returns NaN for any metric where division by zero would occur. + """ + mean_square_ref = roughness.mean_square_ref + mean_square_comp = roughness.mean_square_comp + mean_square_of_difference = roughness.mean_square_of_difference + + with np.errstate(divide="ignore", invalid="ignore"): + roughness_normalized_to_reference = ( + (mean_square_of_difference / mean_square_ref) ** 2 + if mean_square_ref > 0 + else np.nan + ) + roughness_normalized_to_compared = ( + (mean_square_of_difference / mean_square_comp) ** 2 + if mean_square_comp > 0 + else np.nan + ) + roughness_normalized_to_reference_and_compared = ( + mean_square_of_difference**2 / (mean_square_ref * mean_square_comp) + if (mean_square_ref > 0 and mean_square_comp > 0) + else np.nan + ) + return NormalizedSquareBasedRoughnessDifferences( + roughness_normalized_to_reference=roughness_normalized_to_reference, + roughness_normalized_to_compared=roughness_normalized_to_compared, + roughness_normalized_to_reference_and_compared=roughness_normalized_to_reference_and_compared, + ) diff --git a/packages/scratch-core/src/conversion/profile_correlator/transforms.py b/packages/scratch-core/src/conversion/profile_correlator/transforms.py new file mode 100644 index 00000000..6a892bcb --- /dev/null +++ b/packages/scratch-core/src/conversion/profile_correlator/transforms.py @@ -0,0 +1,56 @@ +""" +Transform operations for profile alignment. + +This module provides scaling functions for profile alignment: + +- equalize_pixel_scale: Downsample profiles to matching pixel sizes +- apply_scaling: Apply scale transformation to a profile +""" + +import numpy as np +from conversion.profile_correlator.data_types import Profile +from conversion.resample import resample_array_1d + + +def equalize_pixel_scale( + profile_1: Profile, + profile_2: Profile, +) -> tuple[Profile, Profile]: + """ + Downsample the higher-resolution profile to match pixel sizes. + + This function compares the pixel sizes (sampling distances) of two profiles + and downsamples the one with smaller pixel size (higher resolution) to match + the one with larger pixel size (lower resolution). The downsampling uses + cubic spline interpolation. + + The profile with the larger pixel size is returned unchanged, while the + other profile is downsampled. + + :param profile_1: First profile with heights and pixel_size. + :param profile_2: Second profile with heights and pixel_size. + :returns: Tuple of (profile_1_out, profile_2_out) with equal pixel sizes. + The profile that was downsampled will have updated heights and pixel_size. + """ + pixel_1 = profile_1.pixel_size + pixel_2 = profile_2.pixel_size + + if np.isclose(pixel_1, pixel_2, atol=1e-12): + return profile_1, profile_2 + + # Downsample the higher-resolution profile to match the lower-resolution one + if pixel_1 > pixel_2: + to_downsample, target_pixel_size = profile_2, profile_1.pixel_size + else: + to_downsample, target_pixel_size = profile_1, profile_2.pixel_size + + factor = target_pixel_size / to_downsample.pixel_size + downsampled = Profile( + heights=resample_array_1d(to_downsample.heights, factor), + pixel_size=target_pixel_size, + ) + + if pixel_1 > pixel_2: + return profile_1, downsampled + else: + return downsampled, profile_2 diff --git a/packages/scratch-core/src/conversion/resample.py b/packages/scratch-core/src/conversion/resample.py index ebb12111..72f5d292 100644 --- a/packages/scratch-core/src/conversion/resample.py +++ b/packages/scratch-core/src/conversion/resample.py @@ -1,10 +1,12 @@ from typing import Optional, TypeVar import numpy as np +from scipy.signal import resample as signal_resample from skimage.transform import resize -from conversion.data_formats import Mark + +from container_models.base import BinaryMask, FloatArray1D, FloatArray2D from container_models.scan_image import ScanImage -from container_models.base import BinaryMask, FloatArray2D +from conversion.data_formats import Mark T = TypeVar("T", FloatArray2D, BinaryMask) @@ -42,7 +44,7 @@ def resample_scan_image_and_mask( return scan_image, mask image = _resample_scan_image(scan_image, factors=factors) if mask is not None: - mask = resample_image_array(mask, factors=factors) + mask = resample_array_2d(mask, factors=factors) return image, mask @@ -69,7 +71,7 @@ def _resample_scan_image(image: ScanImage, factors: tuple[float, float]) -> Scan :param factors: The multipliers for the scale of the X- and Y-axis. :returns: The resampled ScanImage. """ - image_array_resampled = resample_image_array(image.data, factors=factors) + image_array_resampled = resample_array_2d(image.data, factors=factors) return ScanImage( data=image_array_resampled, scale_x=image.scale_x * factors[0], @@ -77,12 +79,37 @@ def _resample_scan_image(image: ScanImage, factors: tuple[float, float]) -> Scan ) -def resample_image_array( +def resample_array_1d( + data: FloatArray1D, + factor: float, +) -> FloatArray1D: + """ + Resample a 1D array with anti-aliasing. + + Uses scipy.signal.resample which applies an anti-aliasing filter before + resampling, matching MATLAB's resample behavior. + + :param data: 1D input array. + :param factor: Scale factor for pixel size. factor > 1 means downsampling + (fewer output samples), factor < 1 means upsampling. + :returns: Resampled 1D array of length max(1, round(len(data) / factor)). + """ + n_in = len(data) + n_out = max(1, int(round(n_in / factor))) + + if n_out == n_in: + return data.copy() + + result = signal_resample(data, n_out) # type: ignore[assignment] + return result + + +def resample_array_2d( array: T, factors: tuple[float, float], ) -> T: """ - Resample an array using the specified resampling factors. + Resample a 2D array using the specified resampling factors. For example, if the scale factor is 0.5, then the image output shape will be scaled by 1 / 0.5 = 2. diff --git a/packages/scratch-core/tests/conversion/profile_correlator/__init__.py b/packages/scratch-core/tests/conversion/profile_correlator/__init__.py new file mode 100644 index 00000000..a3c520d4 --- /dev/null +++ b/packages/scratch-core/tests/conversion/profile_correlator/__init__.py @@ -0,0 +1 @@ +"""Tests for the profile_correlator module.""" diff --git a/packages/scratch-core/tests/conversion/profile_correlator/conftest.py b/packages/scratch-core/tests/conversion/profile_correlator/conftest.py new file mode 100644 index 00000000..e49dd178 --- /dev/null +++ b/packages/scratch-core/tests/conversion/profile_correlator/conftest.py @@ -0,0 +1,218 @@ +""" +Shared test fixtures for profile_correlator tests. + +This module provides fixtures for creating synthetic profiles with known +properties, useful for testing the profile correlation functions. +""" + +from pathlib import Path + +import numpy as np +import pytest + +from conversion.profile_correlator import Profile, AlignmentParameters + + +# Directory for test data files (MATLAB .mat files for validation) +DATA_DIR = Path(__file__).parent / "data" + + +@pytest.fixture +def pixel_size_05um() -> float: + """Standard pixel size of 0.5 micrometers in meters.""" + return 0.5e-6 + + +@pytest.fixture +def pixel_size_1um() -> float: + """Standard pixel size of 1.0 micrometer in meters.""" + return 1.0e-6 + + +@pytest.fixture +def simple_sine_profile(pixel_size_05um: float) -> Profile: + """ + Create a simple sinusoidal profile for basic tests. + + The profile contains 1000 samples of a sine wave with some added noise. + """ + np.random.seed(42) + x = np.linspace(0, 10 * np.pi, 1000) + data = np.sin(x) * 1e-6 # Heights in micrometers scale + data += np.random.normal(0, 0.01e-6, len(data)) # Add a small noise + + return Profile(heights=data, pixel_size=pixel_size_05um) + + +@pytest.fixture +def shifted_sine_profile(pixel_size_05um: float) -> Profile: + """ + Create a shifted version of the sine profile for alignment tests. + + The profile is shifted by approximately 20 samples. + """ + np.random.seed(43) # Different seed for different noise + x = np.linspace(0, 10 * np.pi, 1000) + shift = 0.2 # radians, approximately 20 samples + data = np.sin(x + shift) * 1e-6 + data += np.random.normal(0, 0.01e-6, len(data)) + + return Profile(heights=data, pixel_size=pixel_size_05um) + + +@pytest.fixture +def scaled_sine_profile(pixel_size_05um: float) -> Profile: + """ + Create a scaled (stretched) version of the sine profile. + + The profile is scaled by 1.02 (2% stretch). + """ + np.random.seed(44) + # Create profile with 1.02x scale (fewer periods in same length) + x = np.linspace(0, 10 * np.pi / 1.02, 1000) + data = np.sin(x) * 1e-6 + data += np.random.normal(0, 0.01e-6, len(data)) + + return Profile(heights=data, pixel_size=pixel_size_05um) + + +@pytest.fixture +def profile_with_nans(pixel_size_05um: float) -> Profile: + """Create a profile with some NaN values for NaN handling tests.""" + np.random.seed(45) + x = np.linspace(0, 10 * np.pi, 1000) + data = np.sin(x) * 1e-6 + data += np.random.normal(0, 0.01e-6, len(data)) + + # Insert some NaN values + data[100:110] = np.nan # Block of NaNs + data[500] = np.nan # Single NaN + data[700:750] = np.nan # Larger block + + return Profile(heights=data, pixel_size=pixel_size_05um) + + +@pytest.fixture +def partial_profile(pixel_size_05um: float) -> Profile: + """ + Create a partial (shorter) profile for partial matching tests. + + This profile is a subset of a longer reference, starting at index 300 + and having length 400. + """ + np.random.seed(46) + x = np.linspace(0, 10 * np.pi, 1000) + full_data = np.sin(x) * 1e-6 + + # Extract partial segment + partial_data = full_data[300:700].copy() + partial_data += np.random.normal(0, 0.01e-6, len(partial_data)) + + return Profile(heights=partial_data, pixel_size=pixel_size_05um) + + +@pytest.fixture +def different_resolution_profile(pixel_size_1um: float) -> Profile: + """Create a profile with different pixel size for resampling tests.""" + np.random.seed(48) + # Half the number of samples due to double pixel size + x = np.linspace(0, 10 * np.pi, 500) + data = np.sin(x) * 1e-6 + data += np.random.normal(0, 0.01e-6, len(data)) + + return Profile(heights=data, pixel_size=pixel_size_1um) + + +@pytest.fixture +def default_params() -> AlignmentParameters: + """Default alignment parameters for tests.""" + return AlignmentParameters() + + +@pytest.fixture +def fast_params() -> AlignmentParameters: + """Fast alignment parameters for quicker tests.""" + return AlignmentParameters( + max_scaling=0.05, + ) + + +def make_synthetic_striation_profile( + n_samples: int = 1000, + n_striations: int = 20, + amplitude_um: float = 0.5, + noise_level: float = 0.05, + pixel_size_m: float = 0.5e-6, + seed: int | None = None, +) -> Profile: + """ + Create a synthetic striation profile for testing. + + This function generates a profile that mimics the appearance of striation + marks with multiple ridges and valleys. + + :param n_samples: Number of samples in the profile. + :param n_striations: Number of striation features. + :param amplitude_um: Amplitude of striations in micrometers. + :param noise_level: Relative noise level (0 to 1). + :param pixel_size_m: Pixel size in meters. + :param seed: Random seed for reproducibility. + :returns: Profile with synthetic striation data. + """ + if seed is not None: + np.random.seed(seed) + + # Create base profile with multiple frequency components + x = np.linspace(0, n_striations * 2 * np.pi, n_samples) + + # Primary striation pattern + data = np.sin(x) * amplitude_um * 1e-6 + + # Add some harmonics for realism + data += np.sin(2 * x) * amplitude_um * 0.3 * 1e-6 + data += np.sin(0.5 * x) * amplitude_um * 0.5 * 1e-6 + + # Add noise + noise = np.random.normal(0, amplitude_um * noise_level * 1e-6, n_samples) + data += noise + + return Profile(heights=data, pixel_size=pixel_size_m) + + +def make_shifted_profile( + profile: Profile, + shift_samples: float, + scale_factor: float = 1.0, + seed: int | None = None, +) -> Profile: + """ + Create a shifted and optionally scaled version of a profile. + + :param profile: Original profile. + :param shift_samples: Number of samples to shift (can be fractional). + :param scale_factor: Scaling factor (1.0 = no scaling). + :param seed: Random seed for added noise. + :returns: New Profile with shifted/scaled data. + """ + from scipy.interpolate import interp1d + + if seed is not None: + np.random.seed(seed) + + data = profile.heights + n = len(data) + + # Create interpolator + x_orig = np.arange(n) + interpolator = interp1d( + x_orig, data, kind="linear", fill_value=0, bounds_error=False + ) + + # Create new coordinates with shift and scale + x_new = x_orig * scale_factor + shift_samples + new_data = interpolator(x_new) + + # Add a small amount of noise + new_data += np.random.normal(0, np.nanstd(data) * 0.01, n) + + return Profile(heights=new_data, pixel_size=profile.pixel_size) diff --git a/packages/scratch-core/tests/conversion/profile_correlator/test_correlator.py b/packages/scratch-core/tests/conversion/profile_correlator/test_correlator.py new file mode 100644 index 00000000..6e617988 --- /dev/null +++ b/packages/scratch-core/tests/conversion/profile_correlator/test_correlator.py @@ -0,0 +1,292 @@ +"""Tests for the correlator module.""" + +import numpy as np +import pytest +from numpy.testing import assert_allclose +from scipy.interpolate import interp1d + +from conversion.profile_correlator import ( + AlignmentParameters, + StriationComparisonResults, + Profile, + correlate_profiles, +) +from .conftest import make_synthetic_striation_profile, make_shifted_profile + +PIXEL_SIZE_M = 1.5e-6 # 1.5 μm + + +# --- Synthetic profile helpers --- + + +def create_base_profile(n_samples: int = 1000, seed: int = 42) -> np.ndarray: + """Generate a striation-like profile with multiple sine frequencies.""" + np.random.seed(seed) + x = np.linspace(0, 20 * np.pi, n_samples) + data = np.sin(x) * 0.5e-6 + data += np.sin(2.3 * x) * 0.2e-6 + data += np.sin(0.7 * x) * 0.15e-6 + data += np.random.normal(0, 0.01e-6, n_samples) + return data + + +def create_shifted_profiles( + base: np.ndarray, shift_samples: int +) -> tuple[Profile, Profile]: + """Create two profiles with a known shift.""" + n = len(base) + extended_length = n + shift_samples + x = np.linspace(0, 20 * np.pi * extended_length / n, extended_length) + extended = np.sin(x) * 0.5e-6 + np.sin(2.3 * x) * 0.2e-6 + np.sin(0.7 * x) * 0.15e-6 + ref_data = extended[:n].copy() + comp_data = extended[shift_samples : shift_samples + n].copy() + return ( + Profile(heights=ref_data, pixel_size=PIXEL_SIZE_M), + Profile(heights=comp_data, pixel_size=PIXEL_SIZE_M), + ) + + +def create_partial_profiles( + base: np.ndarray, partial_ratio: float +) -> tuple[Profile, Profile]: + """Create profiles where comparison is a subset of reference.""" + n = len(base) + partial_len = int(n * partial_ratio) + start = (n - partial_len) // 2 + return ( + Profile(heights=base.copy(), pixel_size=PIXEL_SIZE_M), + Profile( + heights=base[start : start + partial_len].copy(), pixel_size=PIXEL_SIZE_M + ), + ) + + +def create_scaled_profiles( + base: np.ndarray, scale_factor: float +) -> tuple[Profile, Profile]: + """Create profiles where comparison is stretched.""" + n = len(base) + x_orig = np.arange(n) + interp = interp1d(x_orig, base, kind="cubic", fill_value="extrapolate") # type: ignore[arg-type] + x_scaled = np.arange(n) / scale_factor + comp_data = interp(x_scaled) + return ( + Profile(heights=base.copy(), pixel_size=PIXEL_SIZE_M), + Profile(heights=comp_data, pixel_size=PIXEL_SIZE_M), + ) + + +# --- Basic functionality tests --- + + +class TestCorrelateProfilesBasic: + """Basic functionality tests for correlate_profiles.""" + + def test_returns_comparison_results(self): + """Should return a StriationComparisonResults object.""" + ref = make_synthetic_striation_profile(n_samples=1000, seed=42) + comp = make_shifted_profile(ref, 10.0, seed=43) + result = correlate_profiles(ref, comp, AlignmentParameters()) + assert isinstance(result, StriationComparisonResults) + + def test_correlation_coefficient_valid(self): + """Correlation coefficient should be computed and in valid range.""" + ref = make_synthetic_striation_profile(n_samples=1000, seed=42) + comp = make_shifted_profile(ref, 10.0, seed=43) + result = correlate_profiles(ref, comp, AlignmentParameters()) + assert result is not None + assert not np.isnan(result.correlation_coefficient) + assert -1 <= result.correlation_coefficient <= 1 + + def test_position_shift_computed(self): + """Position shift should be computed.""" + ref = make_synthetic_striation_profile(n_samples=1000, seed=42) + comp = make_shifted_profile(ref, 10.0, seed=43) + result = correlate_profiles(ref, comp, AlignmentParameters()) + assert result is not None + assert not np.isnan(result.position_shift) + + def test_scale_factor_computed(self): + """Scale factor should be computed.""" + ref = make_synthetic_striation_profile(n_samples=1000, seed=42) + comp = make_shifted_profile(ref, 0.0, 1.01, seed=43) + result = correlate_profiles(ref, comp, AlignmentParameters()) + assert result is not None + assert not np.isnan(result.scale_factor) + + def test_roughness_metrics_computed(self): + """Roughness metrics (Sa, Sq) should be computed.""" + ref = make_synthetic_striation_profile(n_samples=1000, seed=42) + comp = make_shifted_profile(ref, 10.0, seed=43) + result = correlate_profiles(ref, comp, AlignmentParameters()) + assert result is not None + assert not np.isnan(result.sa_ref) + assert not np.isnan(result.mean_square_ref) + assert result.sa_ref > 0 + assert result.mean_square_ref > 0 + + def test_overlap_metrics_computed(self): + """Overlap length and ratio should be computed.""" + ref = make_synthetic_striation_profile(n_samples=1000, seed=42) + comp = make_shifted_profile(ref, 10.0, seed=43) + result = correlate_profiles(ref, comp, AlignmentParameters()) + assert result is not None + assert not np.isnan(result.overlap_length) + assert result.overlap_length > 0 + + def test_pixel_sizes_recorded(self): + """Pixel sizes should be recorded in results.""" + ref = make_synthetic_striation_profile(n_samples=1000, seed=42) + comp = make_synthetic_striation_profile(n_samples=1000, seed=43) + result = correlate_profiles(ref, comp, AlignmentParameters()) + assert result is not None + assert_allclose(result.pixel_size, ref.pixel_size, atol=1e-16) + + def test_default_parameters_used(self): + """Should work with default parameters when none provided.""" + ref = make_synthetic_striation_profile(n_samples=1000, seed=42) + comp = make_shifted_profile(ref, 10.0, seed=43) + result = correlate_profiles(ref, comp) + assert result is not None + assert not np.isnan(result.correlation_coefficient) + + +# --- Synthetic profile alignment tests --- + + +class TestIdenticalProfiles: + """Tests for identical profiles.""" + + def test_identical_profiles_perfect_correlation(self): + """Identical profiles should have near-perfect correlation.""" + base = create_base_profile(n_samples=1000, seed=42) + ref = Profile(heights=base.copy(), pixel_size=PIXEL_SIZE_M) + comp = Profile(heights=base.copy(), pixel_size=PIXEL_SIZE_M) + result = correlate_profiles(ref, comp, AlignmentParameters()) + assert result is not None + assert result.correlation_coefficient > 0.999 + assert result.overlap_ratio > 0.99 + + +class TestShiftedProfiles: + """Tests for profiles with translation shifts.""" + + @pytest.mark.parametrize( + "shift_pct,min_corr", + [(3, 0.80), (5, 0.80), (10, 0.70), (20, 0.60), (30, 0.90), (50, 0.80)], + ) + def test_shifted_profiles(self, shift_pct: int, min_corr: float): + """Shifted profiles should align with good correlation.""" + base = create_base_profile(n_samples=1000, seed=42) + ref, comp = create_shifted_profiles(base, int(1000 * shift_pct / 100)) + result = correlate_profiles(ref, comp, AlignmentParameters()) + assert result is not None + assert result.correlation_coefficient >= min_corr + + @pytest.mark.parametrize( + "shift_pct,min_corr", + [(3, 0.75), (5, 0.80), (10, 0.70), (20, 0.60), (30, 0.90), (50, 0.80)], + ) + def test_shifted_profiles_flipped(self, shift_pct: int, min_corr: float): + """Shifted profiles work with swapped order.""" + base = create_base_profile(n_samples=1000, seed=42) + a, b = create_shifted_profiles(base, int(1000 * shift_pct / 100)) + result = correlate_profiles(b, a, AlignmentParameters()) + assert result is not None + assert result.correlation_coefficient >= min_corr + + +class TestPartialProfiles: + """Tests for partial profile matching.""" + + @pytest.mark.parametrize( + "length_pct, expected_overlap", + [ + (50, 3 / 4), + (30, 0.65), + ], + ) + def test_partial_profiles(self, length_pct: int, expected_overlap: float): + """Partial profiles should match with high correlation.""" + base = create_base_profile(n_samples=1000, seed=42) + ref, comp = create_partial_profiles(base, length_pct / 100.0) + result = correlate_profiles(ref, comp, AlignmentParameters()) + assert result is not None + assert result.correlation_coefficient == 1 + assert result.overlap_ratio == pytest.approx(expected_overlap, rel=1e-6) + + @pytest.mark.parametrize( + "length_pct, expected_overlap", + [ + (50, 3 / 4), + (30, 0.65), + ], + ) + def test_partial_profiles_flipped(self, length_pct: int, expected_overlap: float): + """Partial matching works with shorter profile as reference.""" + base = create_base_profile(n_samples=1000, seed=42) + long, short = create_partial_profiles(base, length_pct / 100.0) + result = correlate_profiles(short, long, AlignmentParameters()) + assert result is not None + assert result.correlation_coefficient == 1 + assert result.overlap_ratio == pytest.approx(expected_overlap, rel=1e-6) + + +class TestScaledProfiles: + """Tests for profiles with scaling differences.""" + + @pytest.mark.parametrize("scale_pct", [5, 10, 20]) + def test_scaled_profiles(self, scale_pct: int): + """Scaled profiles should be detected and aligned.""" + base = create_base_profile(n_samples=1000, seed=42) + scale = 1.0 + scale_pct / 100.0 + ref, comp = create_scaled_profiles(base, scale) + params = AlignmentParameters(max_scaling=scale_pct / 100.0) + result = correlate_profiles(ref, comp, params) + assert result is not None + assert result.correlation_coefficient >= 0.999 + assert abs(result.scale_factor - scale) == 0 + + @pytest.mark.parametrize("scale_pct", [5, 10]) + def test_scaled_profiles_flipped(self, scale_pct: int): + """Scaled profiles work with stretched as reference.""" + base = create_base_profile(n_samples=1000, seed=42) + scale = 1.0 + scale_pct / 100.0 + original, stretched = create_scaled_profiles(base, scale) + params = AlignmentParameters(max_scaling=scale_pct / 100.0) + result = correlate_profiles(stretched, original, params) + assert result is not None + assert result.correlation_coefficient >= 0.999 + assert abs(result.scale_factor - 1 / scale) == 0 + + +# --- Edge case tests --- + + +class TestEdgeCases: + """Edge case tests for correlate_profiles.""" + + def test_constant_profile(self): + """Constant profiles return None (no valid correlation possible).""" + ref = Profile(np.ones(500), pixel_size=0.5e-6) + comp = Profile(np.ones(500) * 2, pixel_size=0.5e-6) + result = correlate_profiles(ref, comp, AlignmentParameters()) + assert result is None + + def test_very_short_profiles(self): + """Very short profiles return None when below min_overlap_distance.""" + ref = Profile(np.random.randn(50), pixel_size=0.5e-6) + comp = Profile(np.random.randn(50), pixel_size=0.5e-6) + result = correlate_profiles(ref, comp, AlignmentParameters()) + assert result is None + + def test_different_length_profiles(self): + """Profiles with different lengths should correlate.""" + ref = make_synthetic_striation_profile(n_samples=1000, seed=42) + partial = Profile( + heights=ref.heights[100:900].copy(), + pixel_size=ref.pixel_size, + ) + result = correlate_profiles(ref, partial, AlignmentParameters()) + assert result is not None + assert not np.isnan(result.correlation_coefficient) diff --git a/packages/scratch-core/tests/conversion/profile_correlator/test_profile_pairs.py b/packages/scratch-core/tests/conversion/profile_correlator/test_profile_pairs.py new file mode 100644 index 00000000..b141ae9a --- /dev/null +++ b/packages/scratch-core/tests/conversion/profile_correlator/test_profile_pairs.py @@ -0,0 +1,80 @@ +"""Tests for profile correlation using real measurement data.""" + +from pathlib import Path + +import numpy as np +import pytest + +from conversion.profile_correlator import ( + AlignmentParameters, + Profile, + correlate_profiles, +) + +PROFILES_DIR = ( + Path(__file__).parent.parent.parent + / "resources" + / "profile_correlator" + / "profiles" +) + +PIXEL_SIZE_M = 1.5e-6 # Sample profiles use 1.5 μm pixel size + +# Expected correlation coefficients per profile pair +EXPECTED_CORRELATIONS = { + "sample_02": 0.8542, + "sample_03": 0.8036, + "sample_04": 0.5762, + "sample_05": 0.6500, + "sample_07": 0.5636, + "sample_08": 0.5961, + "sample_09": 0.5608, + "sample_11": 0.8213, + "sample_12": 0.8379, + "sample_14": 0.9809, + "sample_15": 0.9802, + "sample_16": 0.9756, + "sample_18": 0.9365, + "sample_19": 0.9412, + "sample_21": 0.9624, + "sample_22": 0.9371, + "sample_23": 0.9813, +} + + +def discover_profile_pairs() -> list[tuple[str, Path, Path]]: + """Discover all profile pairs in the profiles folder.""" + pairs = [] + for ref_path in sorted(PROFILES_DIR.glob("*_ref.npy")): + name = ref_path.stem.replace("_ref", "") + comp_path = ref_path.parent / f"{name}_comp.npy" + if comp_path.exists(): + pairs.append((name, ref_path, comp_path)) + return pairs + + +PROFILE_PAIRS = discover_profile_pairs() + + +class TestProfilePairs: + """Tests for real profile pairs.""" + + @pytest.mark.parametrize( + "name,ref_path,comp_path", + PROFILE_PAIRS, + ids=[p[0] for p in PROFILE_PAIRS], + ) + def test_correlation(self, name: str, ref_path: Path, comp_path: Path): + """Test correlation matches expected value.""" + ref_data = np.load(ref_path).astype(np.float64) + comp_data = np.load(comp_path).astype(np.float64) + + ref = Profile(heights=ref_data, pixel_size=PIXEL_SIZE_M) + comp = Profile(heights=comp_data, pixel_size=PIXEL_SIZE_M) + + result = correlate_profiles(ref, comp, AlignmentParameters()) + + expected = EXPECTED_CORRELATIONS[name] + assert abs(result.correlation_coefficient - expected) < 0.01, ( + f"Expected {expected:.4f}, got {result.correlation_coefficient:.4f}" + ) diff --git a/packages/scratch-core/tests/conversion/profile_correlator/test_statistics.py b/packages/scratch-core/tests/conversion/profile_correlator/test_statistics.py new file mode 100644 index 00000000..5b034399 --- /dev/null +++ b/packages/scratch-core/tests/conversion/profile_correlator/test_statistics.py @@ -0,0 +1,72 @@ +"""Tests for the similarity module.""" + +import numpy as np +import pytest +from numpy.testing import assert_allclose + +from conversion.profile_correlator import ( + compute_cross_correlation, +) + + +class TestComputeCrossCorrelation: + """Tests for compute_cross_correlation function.""" + + def test_identical_profiles_give_correlation_one(self): + """Identical profiles should have correlation coefficient of 1.0.""" + profile = np.array([1.0, 2.0, 3.0, 4.0, 5.0]) + result = compute_cross_correlation(profile, profile) + assert_allclose(result, 1.0, atol=1e-10) + + def test_negatively_correlated_profiles(self): + """Negatively correlated profiles should have correlation near -1.0.""" + p1 = np.array([1.0, 2.0, 3.0, 4.0, 5.0]) + p2 = np.array([5.0, 4.0, 3.0, 2.0, 1.0]) + result = compute_cross_correlation(p1, p2) + assert_allclose(result, -1.0, atol=1e-10) + + def test_uncorrelated_profiles(self): + """Uncorrelated profiles should have correlation near 0.""" + np.random.seed(42) + p1 = np.random.randn(1000) + p2 = np.random.randn(1000) + result = compute_cross_correlation(p1, p2) + assert abs(result) < 0.1 # Should be close to 0 for random data + + def test_handles_nan_values(self): + """NaN values should be excluded from correlation calculation.""" + p1 = np.array([1.0, np.nan, 3.0, 4.0, 5.0]) + p2 = np.array([1.1, 2.0, 2.9, 4.0, 5.1]) + result = compute_cross_correlation(p1, p2) + # Should still compute valid correlation using non-NaN values + assert not np.isnan(result) + assert result > 0.9 # Should be highly correlated + + def test_all_nan_returns_none(self): + """If all values are NaN, should return None.""" + p1 = np.array([np.nan, np.nan, np.nan]) + p2 = np.array([1.0, 2.0, 3.0]) + result = compute_cross_correlation(p1, p2) + assert result is None + + def test_different_lengths_raises_error(self): + """Profiles with different lengths should raise ValueError.""" + p1 = np.array([1.0, 2.0, 3.0]) + p2 = np.array([1.0, 2.0]) + with pytest.raises(ValueError): + compute_cross_correlation(p1, p2) + + def test_constant_profile_returns_nan(self): + """Constant profile (zero variance) should return NaN.""" + p1 = np.array([5.0, 5.0, 5.0, 5.0, 5.0]) + p2 = np.array([1.0, 2.0, 3.0, 4.0, 5.0]) + result = compute_cross_correlation(p1, p2) + assert np.isnan(result) + + def test_sinusoidal_profiles(self): + """Test correlation of sinusoidal profiles.""" + x = np.linspace(0, 2 * np.pi, 100) + p1 = np.sin(x) + p2 = np.sin(x + 0.1) # Slightly phase-shifted + result = compute_cross_correlation(p1, p2) + assert result > 0.95 # Should be very highly correlated diff --git a/packages/scratch-core/tests/conversion/profile_correlator/test_transforms.py b/packages/scratch-core/tests/conversion/profile_correlator/test_transforms.py new file mode 100644 index 00000000..aca01519 --- /dev/null +++ b/packages/scratch-core/tests/conversion/profile_correlator/test_transforms.py @@ -0,0 +1,41 @@ +"""Tests for the transforms module.""" + +import numpy as np +import pytest +from numpy.testing import assert_array_equal + +from conversion.profile_correlator import ( + Profile, + equalize_pixel_scale, +) + + +class TestEqualizePixelScale: + """Tests for equalize_pixel_scale function.""" + + def test_same_pixel_size_unchanged(self): + """Profiles with same pixel size should be returned unchanged.""" + pixel_size = 0.5e-6 + p1 = Profile(np.random.randn(100), pixel_size=pixel_size) + p2 = Profile(np.random.randn(100), pixel_size=pixel_size) + + p1_out, p2_out = equalize_pixel_scale(p1, p2) + + assert_array_equal(p1_out.heights, p1.heights) + assert_array_equal(p2_out.heights, p2.heights) + + def test_downsamples_higher_resolution_profile(self): + """Higher resolution profile should be downsampled to lower resolution.""" + p1 = Profile(np.random.randn(100), pixel_size=1.0e-6) # Lower resolution + p2 = Profile(np.random.randn(200), pixel_size=0.5e-6) # Higher resolution + + p1_out, p2_out = equalize_pixel_scale(p1, p2) + + # p1 should be unchanged (lower resolution) + assert_array_equal(p1_out.heights, p1.heights) + assert p1_out.pixel_size == p1.pixel_size + + # p2 should be downsampled to match p1's pixel size + assert p2_out.pixel_size == p1.pixel_size + # Length should be approximately halved + assert len(p2_out.heights) == pytest.approx(len(p2.heights) / 2, abs=2) diff --git a/packages/scratch-core/tests/conversion/test_resample.py b/packages/scratch-core/tests/conversion/test_resample.py index ad31678e..7ccc4ec6 100644 --- a/packages/scratch-core/tests/conversion/test_resample.py +++ b/packages/scratch-core/tests/conversion/test_resample.py @@ -8,7 +8,7 @@ _resample_scan_image, get_scaling_factors, _clip_factors, - resample_image_array, + resample_array_2d, resample_mark, ) @@ -44,7 +44,7 @@ def test_calculates_output_shape_correctly(self, mock_resize: MagicMock): array = np.zeros((100, 200)) mock_resize.return_value = np.zeros((50, 100)) - resample_image_array(array, factors=(2.0, 2.0)) + resample_array_2d(array, factors=(2.0, 2.0)) call_args = mock_resize.call_args[1] assert call_args["output_shape"] == (50.0, 100.0) @@ -54,7 +54,7 @@ def test_anti_aliasing_on_upsampling(self, mock_resize: MagicMock): array = np.zeros((100, 100)) mock_resize.return_value = np.zeros((200, 200)) - resample_image_array(array, factors=(0.5, 0.5)) + resample_array_2d(array, factors=(0.5, 0.5)) assert mock_resize.call_args[1]["anti_aliasing"] is False @@ -63,14 +63,14 @@ def test_no_anti_aliasing_on_downsampling(self, mock_resize: MagicMock): array = np.zeros((100, 100)) mock_resize.return_value = np.zeros((50, 50)) - resample_image_array(array, factors=(2.0, 2.0)) + resample_array_2d(array, factors=(2.0, 2.0)) assert mock_resize.call_args[1]["anti_aliasing"] is True class TestResampleScanImage: def test_updates_scales(self, scan_image_rectangular_with_nans: ScanImage): - with patch("conversion.resample.resample_image_array") as mock: + with patch("conversion.resample.resample_array_2d") as mock: mock.return_value = np.zeros((50, 50)) result = _resample_scan_image(scan_image_rectangular_with_nans, (2.0, 2.0)) @@ -138,7 +138,7 @@ def test_resamples_mask_when_provided( ): mask = np.ones((100, 100), dtype=np.bool_) - with patch("conversion.resample.resample_image_array") as mock: + with patch("conversion.resample.resample_array_2d") as mock: mock.return_value = np.zeros((50, 50)) _, result_mask = resample_scan_image_and_mask( diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_02_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_02_comp.npy new file mode 100644 index 00000000..5a1941a1 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_02_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_02_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_02_ref.npy new file mode 100644 index 00000000..0654faf8 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_02_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_03_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_03_comp.npy new file mode 100644 index 00000000..928f1611 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_03_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_03_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_03_ref.npy new file mode 100644 index 00000000..ee8d5d46 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_03_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_04_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_04_comp.npy new file mode 100644 index 00000000..2de904d4 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_04_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_04_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_04_ref.npy new file mode 100644 index 00000000..4af5aade Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_04_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_05_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_05_comp.npy new file mode 100644 index 00000000..1343932c Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_05_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_05_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_05_ref.npy new file mode 100644 index 00000000..60414bb3 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_05_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_07_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_07_comp.npy new file mode 100644 index 00000000..49215203 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_07_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_07_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_07_ref.npy new file mode 100644 index 00000000..89642e01 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_07_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_08_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_08_comp.npy new file mode 100644 index 00000000..642e5803 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_08_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_08_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_08_ref.npy new file mode 100644 index 00000000..956fbd30 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_08_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_09_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_09_comp.npy new file mode 100644 index 00000000..b34165b1 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_09_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_09_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_09_ref.npy new file mode 100644 index 00000000..839c57de Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_09_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_11_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_11_comp.npy new file mode 100644 index 00000000..a6ffbd43 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_11_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_11_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_11_ref.npy new file mode 100644 index 00000000..51a9e89d Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_11_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_12_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_12_comp.npy new file mode 100644 index 00000000..956fbd30 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_12_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_12_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_12_ref.npy new file mode 100644 index 00000000..6fb03b59 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_12_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_14_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_14_comp.npy new file mode 100644 index 00000000..bb561072 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_14_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_14_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_14_ref.npy new file mode 100644 index 00000000..feb28efc Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_14_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_15_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_15_comp.npy new file mode 100644 index 00000000..993bc346 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_15_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_15_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_15_ref.npy new file mode 100644 index 00000000..3e0436ac Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_15_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_16_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_16_comp.npy new file mode 100644 index 00000000..ba457828 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_16_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_16_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_16_ref.npy new file mode 100644 index 00000000..5d686f3b Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_16_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_18_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_18_comp.npy new file mode 100644 index 00000000..0da0ee89 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_18_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_18_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_18_ref.npy new file mode 100644 index 00000000..937c0cb7 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_18_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_19_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_19_comp.npy new file mode 100644 index 00000000..05b2de28 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_19_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_19_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_19_ref.npy new file mode 100644 index 00000000..429f16e9 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_19_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_21_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_21_comp.npy new file mode 100644 index 00000000..7cdec28e Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_21_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_21_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_21_ref.npy new file mode 100644 index 00000000..91ab1833 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_21_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_22_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_22_comp.npy new file mode 100644 index 00000000..502cada0 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_22_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_22_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_22_ref.npy new file mode 100644 index 00000000..4ae3f3d4 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_22_ref.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_23_comp.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_23_comp.npy new file mode 100644 index 00000000..5d2cde9e Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_23_comp.npy differ diff --git a/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_23_ref.npy b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_23_ref.npy new file mode 100644 index 00000000..75421272 Binary files /dev/null and b/packages/scratch-core/tests/resources/profile_correlator/profiles/sample_23_ref.npy differ