Correlate two profiles of striation marks #98
+1,470
−14
There are no files selected for viewing
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,66 @@ | ||
| """ | ||
laurensWe marked this conversation as resolved.
Show resolved
Hide resolved
|
||
| 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", | ||
| ] | ||
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change | ||||
|---|---|---|---|---|---|---|
| @@ -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))) | ||||||
|
Collaborator
There was a problem hiding this comment. Comments in code usually suggests you should perhaps use more functions: def _get_scale_factors(max_scaling: float, n_steps: int) -> NDArray:
"""Generate symmetric scale factors."""
scales = np.linspace(1.0 - max_scaling, 1.0 + max_scaling, n_steps)
scales = np.unique(np.concatenate((scale_factors, 1 / scale_factors)))
return scales |
||||||
|
|
||||||
| # 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 | ||||||
|
Collaborator
There was a problem hiding this comment.
Suggested change
|
||||||
| """ | ||||||
|
|
||||||
| 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 | ||||||
|
Collaborator
There was a problem hiding this comment. I understand it was like this in MATLAB, so this is not a suggestion for this PR, but: shouldn't the overlap length be a variable as well (ranging from |
||||||
|
|
||||||
|
|
||||||
| 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 | ||||||
|
|
||||||
|
Collaborator
There was a problem hiding this comment. NTH: You could consider splitting the function at the point into two separate functions by returning the parameters found here. Now your function is doing more than one thing (which is generally undesirable). I think for now it's okay, since the intended meaning is quite clear, but in general the functions become very "script"-like if they do too many things |
||||||
| # 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, | ||||||
| ) | ||||||
Oops, something went wrong.
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
Uh oh!
There was an error while loading. Please reload this page.