From b891053fabe85a2ee1e9ed4b75245e2fc10d95ad Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Thu, 11 Dec 2025 11:40:54 -0500 Subject: [PATCH 01/16] refactored cross corr code to be in mvrecon instead of BS --- .../constellation/overlap/SimpleBoundingBoxOverlap.java | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/main/java/net/preibisch/mvrecon/process/interestpointregistration/pairwise/constellation/overlap/SimpleBoundingBoxOverlap.java b/src/main/java/net/preibisch/mvrecon/process/interestpointregistration/pairwise/constellation/overlap/SimpleBoundingBoxOverlap.java index 8c70e0dfc..9b7797eea 100644 --- a/src/main/java/net/preibisch/mvrecon/process/interestpointregistration/pairwise/constellation/overlap/SimpleBoundingBoxOverlap.java +++ b/src/main/java/net/preibisch/mvrecon/process/interestpointregistration/pairwise/constellation/overlap/SimpleBoundingBoxOverlap.java @@ -164,10 +164,10 @@ public static BoundingBox getBoundingBox( final Dimensions dims, final AffineTra { final RealInterval interval = getBoundingBoxReal( dims, transform ); - final int[] minInt = new int[ 3 ]; - final int[] maxInt = new int[ 3 ]; + final int[] minInt = new int[ interval.numDimensions() ]; + final int[] maxInt = new int[ interval.numDimensions() ]; - for ( int d = 0; d < dims.numDimensions(); ++d ) + for ( int d = 0; d < interval.numDimensions(); ++d ) { minInt[ d ] = (int)Math.round( interval.realMin( d ) ) - 1; maxInt[ d ] = (int)Math.round( interval.realMax( d ) ) + 1; From 65d7b2d1993967d6433ada065a016526597d0184 Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Thu, 11 Dec 2025 11:47:22 -0500 Subject: [PATCH 02/16] removing extraneous change --- .../constellation/overlap/SimpleBoundingBoxOverlap.java | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/main/java/net/preibisch/mvrecon/process/interestpointregistration/pairwise/constellation/overlap/SimpleBoundingBoxOverlap.java b/src/main/java/net/preibisch/mvrecon/process/interestpointregistration/pairwise/constellation/overlap/SimpleBoundingBoxOverlap.java index 9b7797eea..a0b8995da 100644 --- a/src/main/java/net/preibisch/mvrecon/process/interestpointregistration/pairwise/constellation/overlap/SimpleBoundingBoxOverlap.java +++ b/src/main/java/net/preibisch/mvrecon/process/interestpointregistration/pairwise/constellation/overlap/SimpleBoundingBoxOverlap.java @@ -164,10 +164,10 @@ public static BoundingBox getBoundingBox( final Dimensions dims, final AffineTra { final RealInterval interval = getBoundingBoxReal( dims, transform ); - final int[] minInt = new int[ interval.numDimensions() ]; - final int[] maxInt = new int[ interval.numDimensions() ]; + final int[] minInt = new int[ 3 ]; + final int[] maxInt = new int[ 3 ]; - for ( int d = 0; d < interval.numDimensions(); ++d ) + for ( int d = 0; d < dims.numDimensions(); ++d ) { minInt[ d ] = (int)Math.round( interval.realMin( d ) ) - 1; maxInt[ d ] = (int)Math.round( interval.realMax( d ) ) + 1; From 080e2ccf980b933234dfb20b1012fa18535033a0 Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Thu, 11 Dec 2025 17:23:37 -0500 Subject: [PATCH 03/16] first draft; appears buggy and slow --- .../mvrecon/fiji/spimdata/explorer/ViewSetupExplorerPanel.java | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/ViewSetupExplorerPanel.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/ViewSetupExplorerPanel.java index 67dabc81f..8eba44b38 100644 --- a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/ViewSetupExplorerPanel.java +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/ViewSetupExplorerPanel.java @@ -82,6 +82,7 @@ import net.preibisch.mvrecon.fiji.spimdata.XmlIoSpimData2; import net.preibisch.mvrecon.fiji.spimdata.explorer.bdv.ScrollableBrightnessDialog; import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.AnalyzeErrorPopup; +import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.AnalyzeOverlapCrossCorrelationPopup; import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.ApplyTransformationPopup; import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.BDVPopup; import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.BakeManualTransformationPopup; @@ -724,6 +725,7 @@ public ArrayList< ExplorerWindowSetable > initPopups() popups.add( new LabelPopUp( " Interest Points" ) ); popups.add( new InterestPointsExplorerPopup() ); popups.add( new AnalyzeErrorPopup() ); + popups.add( new AnalyzeOverlapCrossCorrelationPopup() ); popups.add( new RemoveDetectionsPopup() ); popups.add( new VisualizeDetectionsPopup() ); popups.add( new Separator() ); From e5a39b9b681f5067c9ff683843521603ec6f9f78 Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Thu, 11 Dec 2025 17:46:31 -0500 Subject: [PATCH 04/16] fixed parallelization race condition --- dependency-reduced-pom.xml | 171 ++++ .../AnalyzeOverlapCrossCorrelationPopup.java | 583 ++++++++++++ ...endedExtendedMirroredRandomAccesible2.java | 246 +++++ .../FourNeighborhoodExtrema.java | 318 +++++++ .../process/phasecorrelation/ImgLib2Util.java | 137 +++ .../phasecorrelation/PhaseCorrelation2.java | 313 +++++++ .../PhaseCorrelation2Util.java | 882 ++++++++++++++++++ .../PhaseCorrelationPeak2.java | 406 ++++++++ .../phasecorrelation/deprecated/Blending.java | 78 ++ .../deprecated/BlendingRealRandomAccess.java | 296 ++++++ .../FourNeighborhoodExtremaTest.java | 344 +++++++ .../PhaseCorrelationTest.java | 187 ++++ 12 files changed, 3961 insertions(+) create mode 100644 dependency-reduced-pom.xml create mode 100644 src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/AnalyzeOverlapCrossCorrelationPopup.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/phasecorrelation/BlendedExtendedMirroredRandomAccesible2.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/phasecorrelation/FourNeighborhoodExtrema.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/phasecorrelation/ImgLib2Util.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/phasecorrelation/PhaseCorrelation2.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/phasecorrelation/PhaseCorrelation2Util.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/phasecorrelation/PhaseCorrelationPeak2.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/phasecorrelation/deprecated/Blending.java create mode 100644 src/main/java/net/preibisch/mvrecon/process/phasecorrelation/deprecated/BlendingRealRandomAccess.java create mode 100644 src/test/java/net/preibisch/mvrecon/process/phasecorrelation/FourNeighborhoodExtremaTest.java create mode 100644 src/test/java/net/preibisch/mvrecon/process/phasecorrelation/PhaseCorrelationTest.java diff --git a/dependency-reduced-pom.xml b/dependency-reduced-pom.xml new file mode 100644 index 000000000..3b09f9b31 --- /dev/null +++ b/dependency-reduced-pom.xml @@ -0,0 +1,171 @@ + + + + pom-scijava + org.scijava + 43.0.0 + pom.xml + + 4.0.0 + net.preibisch + multiview-reconstruction + Multiview Reconstruction + 8.1.2-SNAPSHOT + Software for the reconstruction of multi-view microscopic acquisitions +like Selective Plane Illumination Microscopy (SPIM) Data. + https://github.com/JaneliaSciComp/multiview-reconstruction + + GitHub Issues + https://github.com/JaneliaSciComp/multiview-reconstruction/issues + + + GitHub Actions + https://github.com/JaneliaSciComp/multiview-reconstruction/actions + + 2012 + + + Image.sc Forum + https://forum.image.sc/tag/multiview-reconstruction + + + + + StephanPreibisch + Stephan Preibisch + https://imagej.net/people/StephanPreibisch + + founder + lead + developer + debugger + reviewer + support + maintainer + + + + hoerldavid + David Hoerl + https://imagej.net/people/hoerldavid + + developer + debugger + reviewer + support + maintainer + + + + + + Tobias Pietzsch + https://imagej.net/people/tpietzsch + + tpietzsch + + + + Curtis Rueden + https://imagej.net/people/ctrueden + + ctrueden + + + + + + GNU General Public License v2+ + https://www.gnu.org/licenses/old-licenses/gpl-2.0.en.html + repo + + + + scm:git:https://github.com/JaneliaSciComp/multiview-reconstruction + scm:git:git@github.com:JaneliaSciComp/multiview-reconstruction + https://github.com/JaneliaSciComp/multiview-reconstruction + + + Preibisch Lab + http://preibischlab.mdc-berlin.de + + + + fatjar + + + + maven-shade-plugin + + + package + + shade + + + + + + + META-INF/json/mpicbg.spim.data.generic.sequence.ImgLoaderIo + + + + + *:* + + META-INF/*.SF + META-INF/*.DSA + META-INF/*.RSA + + + + + + org.apache.commons.compress + org.janelia.saalfeldlab.org.apache.commons.compress + + + + + + + + + + + scijava.public + https://maven.scijava.org/content/groups/public + + + + + ch.qos.logback + logback-classic + 1.3.15 + test + + + activation + javax.activation + + + logback-core + ch.qos.logback + + + + + + 8.0.0 + net.preibisch.mvrecon + 1.0.0-beta-20 + Multiview Reconstruction developers. + 10.6.5 + 0.6.2 + 1.1.0 + gpl_v2 + sign,deploy-to-scijava + 1.1.0 + + diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/AnalyzeOverlapCrossCorrelationPopup.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/AnalyzeOverlapCrossCorrelationPopup.java new file mode 100644 index 000000000..8a68576c7 --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/AnalyzeOverlapCrossCorrelationPopup.java @@ -0,0 +1,583 @@ +/*- + * #%L + * Software for the reconstruction of multi-view microscopic acquisitions + * like Selective Plane Illumination Microscopy (SPIM) Data. + * %% + * Copyright (C) 2012 - 2025 Multiview Reconstruction developers. + * %% + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program. If not, see + * . + * #L% + */ +package net.preibisch.mvrecon.fiji.spimdata.explorer.popup; + +import java.awt.event.ActionEvent; +import java.awt.event.ActionListener; +import java.util.ArrayList; +import java.util.Collection; +import java.util.Collections; +import java.util.Comparator; +import java.util.HashMap; +import java.util.List; +import java.util.Map; +import java.util.concurrent.Callable; +import java.util.concurrent.ExecutorService; +import java.util.concurrent.Executors; +import java.util.concurrent.Future; + +import javax.swing.JMenuItem; + +import ij.gui.GenericDialog; +import mpicbg.spim.data.generic.sequence.AbstractSequenceDescription; +import mpicbg.spim.data.generic.sequence.BasicImgLoader; +import mpicbg.spim.data.generic.sequence.BasicViewDescription; +import mpicbg.spim.data.sequence.MultiResolutionImgLoader; +import mpicbg.spim.data.registration.ViewRegistration; +import mpicbg.spim.data.registration.ViewRegistrations; +import mpicbg.spim.data.sequence.ViewId; +import net.imglib2.Cursor; +import net.imglib2.FinalInterval; +import net.imglib2.Interval; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.realtransform.AffineTransform3D; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.view.Views; +import net.preibisch.legacy.io.IOFunctions; +import net.preibisch.mvrecon.fiji.spimdata.SpimData2; +import net.preibisch.mvrecon.fiji.spimdata.boundingbox.BoundingBox; +import net.preibisch.mvrecon.fiji.spimdata.explorer.ExplorerWindow; +import net.preibisch.mvrecon.fiji.spimdata.interestpoints.CorrespondingInterestPoints; +import net.preibisch.mvrecon.fiji.spimdata.interestpoints.ViewInterestPointLists; +import net.preibisch.mvrecon.process.boundingbox.BoundingBoxMaximalGroupOverlap; +import net.preibisch.mvrecon.process.downsampling.DownsampleTools; +import net.preibisch.mvrecon.process.fusion.transformed.FusedRandomAccessibleInterval; +import net.preibisch.mvrecon.process.fusion.transformed.TransformView; +import net.preibisch.mvrecon.process.fusion.transformed.TransformVirtual; +import net.preibisch.mvrecon.process.interestpointregistration.pairwise.constellation.grouping.Group; +import net.preibisch.mvrecon.process.phasecorrelation.PhaseCorrelation2Util; + +public class AnalyzeOverlapCrossCorrelationPopup extends JMenuItem implements ExplorerWindowSetable +{ + private static final long serialVersionUID = 1L; + + ExplorerWindow< ? > panel; + + public static int defaultDownsamplingChoiceIndex = 0; // First available resolution + public static int defaultMinCorrespondences = 10; + + @Override + public JMenuItem setExplorerWindow( ExplorerWindow< ? > panel ) + { + this.panel = panel; + return this; + } + + public AnalyzeOverlapCrossCorrelationPopup() + { + super( "Analyze Overlap Cross-Correlation ..." ); + this.addActionListener( new MyActionListener() ); + } + + public class MyActionListener implements ActionListener + { + @Override + public void actionPerformed( ActionEvent e ) + { + if ( panel == null ) + { + IOFunctions.println( "Panel not set for " + this.getClass().getSimpleName() ); + return; + } + + if ( !SpimData2.class.isInstance( panel.getSpimData() ) ) + { + IOFunctions.println( "Only supported for SpimData2 objects: " + this.getClass().getSimpleName() ); + return; + } + + final SpimData2 spimData = (SpimData2) panel.getSpimData(); + final List< ViewId > viewIds = ApplyTransformationPopup.getSelectedViews( panel ); + + if ( viewIds.isEmpty() ) + { + IOFunctions.println( "No views selected." ); + return; + } + + // Get available downsampling levels from dataset metadata + final String[] downsamplingChoices = DownsampleTools.availableDownsamplings( spimData, viewIds.get( 0 ) ); + + // Ask user for parameters + final GenericDialog gd = new GenericDialog( "Cross-Correlation Analysis Parameters" ); + gd.addChoice( "Downsampling", downsamplingChoices, downsamplingChoices[ Math.min( defaultDownsamplingChoiceIndex, downsamplingChoices.length - 1 ) ] ); + gd.addNumericField( "Minimum_correspondences (0 = analyze all overlaps)", defaultMinCorrespondences, 0 ); + gd.addMessage( "Note: Analysis will be performed on the maximal bounding box\nof the overlap between each pair of selected tiles." ); + gd.addMessage( "Set minimum correspondences to 0 to analyze all overlapping pairs\nwithout requiring pre-computed interest points." ); + + gd.showDialog(); + + if ( gd.wasCanceled() ) + return; + + defaultDownsamplingChoiceIndex = gd.getNextChoiceIndex(); + defaultMinCorrespondences = (int) gd.getNextNumber(); + + + new Thread( new Runnable() + { + @Override + public void run() + { + analyzeCrossCorrelations( spimData, viewIds, defaultDownsamplingChoiceIndex, defaultMinCorrespondences ); + } + } ).start(); + } + } + + public static class OverlapResult + { + public ViewId viewA; + public ViewId viewB; + public int numCorrespondences; + public double crossCorrelation; + public double avgIntensityA; + public double avgIntensityB; + public String label; + public int downsamplingIndex; + + public OverlapResult( ViewId viewA, ViewId viewB, String label, int numCorr, double cc, double avgA, double avgB, int dsIndex ) + { + this.viewA = viewA; + this.viewB = viewB; + this.label = label; + this.numCorrespondences = numCorr; + this.crossCorrelation = cc; + this.avgIntensityA = avgA; + this.avgIntensityB = avgB; + this.downsamplingIndex = dsIndex; + } + } + + public static void analyzeCrossCorrelations( + final SpimData2 spimData, + final List< ViewId > viewIds, + final int downsamplingIndex, + final int minCorrespondences ) + { + IOFunctions.println( "Starting Cross-Correlation Analysis..." ); + IOFunctions.println( "Downsampling index: " + downsamplingIndex ); + IOFunctions.println( "Minimum correspondences: " + minCorrespondences ); + IOFunctions.println( "Number of views: " + viewIds.size() ); + + final ExecutorService service = Executors.newFixedThreadPool( Runtime.getRuntime().availableProcessors() ); + final List< OverlapResult > results = new ArrayList<>(); + + final String selectedLabel; + + // Only load interest points if minCorrespondences > 0 + if ( minCorrespondences > 0 ) + { + // Load all interest points and correspondences + IOFunctions.println( "Loading interest points and correspondences..." ); + final Map< ViewId, ViewInterestPointLists > interestPoints = new HashMap<>(); + + for ( final ViewId viewId : viewIds ) + { + final ViewInterestPointLists vipl = spimData.getViewInterestPoints().getViewInterestPointLists( viewId ); + interestPoints.put( viewId, vipl ); + } + + // Get available labels (interest point types) + final Map< String, Integer > labelCounts = new HashMap<>(); + for ( final ViewId viewId : viewIds ) + { + final ViewInterestPointLists vipl = interestPoints.get( viewId ); + for ( final String label : vipl.getHashMap().keySet() ) + { + labelCounts.put( label, labelCounts.getOrDefault( label, 0 ) + 1 ); + } + } + + // Find the most common label + String tempLabel = null; + int maxCount = 0; + for ( final Map.Entry< String, Integer > entry : labelCounts.entrySet() ) + { + if ( entry.getValue() > maxCount ) + { + maxCount = entry.getValue(); + tempLabel = entry.getKey(); + } + } + + if ( tempLabel == null ) + { + IOFunctions.println( "ERROR: No interest points found in selected views." ); + service.shutdown(); + return; + } + + selectedLabel = tempLabel; + IOFunctions.println( "Using interest point label: " + selectedLabel ); + } + else + { + IOFunctions.println( "Minimum correspondences = 0, analyzing all overlapping pairs without checking interest points." ); + selectedLabel = null; + } + + // Analyze all pairs in parallel + final List< Future< OverlapResult > > futures = new ArrayList<>(); + + for ( int i = 0; i < viewIds.size() - 1; ++i ) + { + for ( int j = i + 1; j < viewIds.size(); ++j ) + { + final ViewId viewA = viewIds.get( i ); + final ViewId viewB = viewIds.get( j ); + + // Submit pair analysis as a parallel task + futures.add( service.submit( new Callable< OverlapResult >() + { + @Override + public OverlapResult call() + { + return analyzePair( spimData, viewA, viewB, selectedLabel, + downsamplingIndex, minCorrespondences, service ); + } + } ) ); + } + } + + // Collect results from all tasks + for ( Future< OverlapResult > future : futures ) + { + try + { + final OverlapResult result = future.get(); + if ( result != null ) + results.add( result ); + } + catch ( Exception e ) + { + IOFunctions.println( "Error processing pair: " + e.getMessage() ); + e.printStackTrace(); + } + } + + service.shutdown(); + + // Sort by cross-correlation (ascending) + Collections.sort( results, new Comparator< OverlapResult >() + { + @Override + public int compare( OverlapResult o1, OverlapResult o2 ) + { + return Double.compare( o1.crossCorrelation, o2.crossCorrelation ); + } + } ); + + // Print results + IOFunctions.println( "\n========================================" ); + IOFunctions.println( "Cross-Correlation Analysis Results" ); + IOFunctions.println( "========================================" ); + if ( minCorrespondences > 0 ) + { + IOFunctions.println( String.format( "%-30s %-15s %-15s %-15s %-15s %-15s", + "Tile Pair", "# Corresp.", "Cross-Corr", "Avg Int. A", "Avg Int. B", "Label" ) ); + } + else + { + IOFunctions.println( String.format( "%-30s %-15s %-15s %-15s", + "Tile Pair", "Cross-Corr", "Avg Int. A", "Avg Int. B" ) ); + } + IOFunctions.println( "----------------------------------------" ); + + for ( final OverlapResult result : results ) + { + final String pairName = getTileName( spimData, result.viewA ) + " <-> " + getTileName( spimData, result.viewB ); + if ( minCorrespondences > 0 ) + { + IOFunctions.println( String.format( "%-30s %-15d %-15.4f %-15.2f %-15.2f %-15s", + pairName, + result.numCorrespondences, + result.crossCorrelation, + result.avgIntensityA, + result.avgIntensityB, + result.label != null ? result.label : "N/A" ) ); + } + else + { + IOFunctions.println( String.format( "%-30s %-15.4f %-15.2f %-15.2f", + pairName, + result.crossCorrelation, + result.avgIntensityA, + result.avgIntensityB ) ); + } + } + + IOFunctions.println( "========================================\n" ); + } + + private static OverlapResult analyzePair( + final SpimData2 spimData, + final ViewId viewA, + final ViewId viewB, + final String label, + final int downsamplingIndex, + final int minCorrespondences, + final ExecutorService service ) + { + int numCorr = 0; + + // Only check correspondences if minCorrespondences > 0 + if ( minCorrespondences > 0 ) + { + // Get interest points + final ViewInterestPointLists viplA = spimData.getViewInterestPoints().getViewInterestPointLists( viewA ); + final ViewInterestPointLists viplB = spimData.getViewInterestPoints().getViewInterestPointLists( viewB ); + + if ( viplA.getInterestPointList( label ) == null || viplB.getInterestPointList( label ) == null ) + return null; + + // Count correspondences between A and B + final Collection< CorrespondingInterestPoints > corrsA = viplA.getInterestPointList( label ).getCorrespondingInterestPointsCopy(); + + for ( final CorrespondingInterestPoints cip : corrsA ) + { + if ( cip.getCorrespondingViewId().equals( viewB ) && cip.getCorrespodingLabel().equals( label ) ) + numCorr++; + } + + if ( numCorr < minCorrespondences ) + { + IOFunctions.println( getTileName( spimData, viewA ) + " <-> " + getTileName( spimData, viewB ) + + ": " + numCorr + " correspondences (< " + minCorrespondences + "), skipping." ); + return null; + } + + IOFunctions.println( "Processing: " + getTileName( spimData, viewA ) + " <-> " + getTileName( spimData, viewB ) + + " (" + numCorr + " correspondences)" ); + } + else + { + IOFunctions.println( "Processing: " + getTileName( spimData, viewA ) + " <-> " + getTileName( spimData, viewB ) ); + } + + // Compute overlap bounding box + final List< List< ViewId > > viewGroups = new ArrayList<>(); + final List< ViewId > groupA = new ArrayList<>(); + groupA.add( viewA ); + final List< ViewId > groupB = new ArrayList<>(); + groupB.add( viewB ); + viewGroups.add( groupA ); + viewGroups.add( groupB ); + + final BoundingBoxMaximalGroupOverlap< ViewId > bbDet = new BoundingBoxMaximalGroupOverlap<>( viewGroups, spimData ); + final BoundingBox bbOverlap = bbDet.estimate( "overlap" ); + + if ( bbOverlap == null ) + { + IOFunctions.println( " No overlap found, skipping." ); + return null; + } + + // Use downsampling index to open image at specified resolution level + // No need to specify downsampling factors - the index handles it + + try + { + final RandomAccessibleInterval< FloatType > imgA = loadImageInBB( + spimData.getSequenceDescription(), + spimData.getViewRegistrations(), + viewA, + bbOverlap, + downsamplingIndex ); + + final RandomAccessibleInterval< FloatType > imgB = loadImageInBB( + spimData.getSequenceDescription(), + spimData.getViewRegistrations(), + viewB, + bbOverlap, + downsamplingIndex ); + + // Compute downsampling factor for adaptive sampling + final double downsamplingFactor = Math.pow( 2, downsamplingIndex ); + + // Compute average intensities and cross-correlation + final double avgA = computeAverageIntensity( imgA, downsamplingFactor ); + final double avgB = computeAverageIntensity( imgB, downsamplingFactor ); + final double cc = PhaseCorrelation2Util.getCorrelation( imgA, imgB ); + + IOFunctions.println( " Cross-correlation: " + cc + ", Avg intensities: " + avgA + ", " + avgB ); + + return new OverlapResult( viewA, viewB, label, numCorr, cc, avgA, avgB, downsamplingIndex ); + } + catch ( Exception e ) + { + IOFunctions.println( " Error processing pair: " + e.getMessage() ); + e.printStackTrace(); + return null; + } + } + + private static RandomAccessibleInterval< FloatType > loadImageInBB( + final AbstractSequenceDescription< ?, ?, ? > sd, + final ViewRegistrations vrs, + final ViewId viewId, + final Interval boundingBox, + final int downsamplingIndex ) + { + final BasicImgLoader imgLoader = sd.getImgLoader(); + final ViewRegistration vr = vrs.getViewRegistration( viewId ); + + // Synchronize access to ViewRegistration to prevent race conditions when running in parallel + final AffineTransform3D model; + synchronized ( vr ) + { + vr.updateModel(); + model = vr.getModel().copy(); + } + + // Load image at the specific resolution level chosen by user + RandomAccessibleInterval inputImg; + + if ( imgLoader instanceof MultiResolutionImgLoader ) + { + final MultiResolutionImgLoader mrImgLoader = (MultiResolutionImgLoader) imgLoader; + final int setupId = viewId.getViewSetupId(); + final int timepointId = viewId.getTimePointId(); + + // Load image at user-specified level + inputImg = mrImgLoader.getSetupImgLoader( setupId ).getImage( timepointId, downsamplingIndex ); + + // Get and apply the mipmap transform for this level + final AffineTransform3D mipmapTransform = mrImgLoader.getSetupImgLoader( setupId ).getMipmapTransforms()[ downsamplingIndex ]; + model.concatenate( mipmapTransform ); + + // Use bounding box in world coordinates (do NOT scale it) + // The model transform (with mipmap concatenated) handles the coordinate mapping + // NOTE: The result image will have the same dimensions regardless of downsampling level + // because the bounding box is in world coordinates. The speedup from downsampling comes + // from faster access to the pre-downsampled input image, not from smaller output dimensions. + final Interval bb = new FinalInterval( boundingBox ); + + // Transform view to bounding box + return TransformView.transformView( inputImg, model, bb, 0, 1 ); + } + else + { + // Fallback for non-multiresolution loaders: load full resolution + inputImg = imgLoader.getSetupImgLoader( viewId.getViewSetupId() ).getImage( viewId.getTimePointId() ); + + // Compute downsampling factor: 2^index + final double downsamplingFactor = Math.pow( 2, downsamplingIndex ); + final double[] downsamplingFactors = new double[ 3 ]; + for ( int d = 0; d < 3; ++d ) + downsamplingFactors[ d ] = downsamplingFactor; + + // Scale bounding box + final Interval bbSc = TransformVirtual.scaleBoundingBox( new FinalInterval( boundingBox ), inverse( downsamplingFactors ) ); + + // Scale transform + TransformVirtual.scaleTransform( model, inverse( downsamplingFactors ) ); + + // Transform view to bounding box + return TransformView.transformView( inputImg, model, bbSc, 0, 1 ); + } + } + + private static double[] inverse( double[] in ) + { + final double[] res = new double[ in.length ]; + for ( int i = 0; i < in.length; i++ ) + res[ i ] = 1.0 / in[ i ]; + return res; + } + + private static double computeAverageIntensity( final RandomAccessibleInterval< FloatType > img, final double downsamplingFactor ) + { + // Sample every Nth pixel in each dimension for speed + // Scale sampling step with downsampling level: at 8x downsampling, sample 2x less frequently than at 4x + final int baseStep = 4; + final int step = Math.max( 1, (int) (baseStep * Math.sqrt( downsamplingFactor )) ); + + double sum = 0; + long count = 0; + + final long[] min = new long[ img.numDimensions() ]; + final long[] max = new long[ img.numDimensions() ]; + img.min( min ); + img.max( max ); + + final net.imglib2.RandomAccess< FloatType > ra = img.randomAccess(); + final long[] pos = new long[ img.numDimensions() ]; + + // Sample pixels at regular intervals + if ( img.numDimensions() == 3 ) + { + for ( long z = min[2]; z <= max[2]; z += step ) + { + pos[2] = z; + for ( long y = min[1]; y <= max[1]; y += step ) + { + pos[1] = y; + for ( long x = min[0]; x <= max[0]; x += step ) + { + pos[0] = x; + ra.setPosition( pos ); + sum += ra.get().get(); + count++; + } + } + } + } + else if ( img.numDimensions() == 2 ) + { + for ( long y = min[1]; y <= max[1]; y += step ) + { + pos[1] = y; + for ( long x = min[0]; x <= max[0]; x += step ) + { + pos[0] = x; + ra.setPosition( pos ); + sum += ra.get().get(); + count++; + } + } + } + else + { + // Fallback for other dimensions: use cursor (slower) + final Cursor< FloatType > cursor = Views.iterable( img ).cursor(); + while ( cursor.hasNext() ) + { + cursor.fwd(); + sum += cursor.get().get(); + count++; + } + } + + return count > 0 ? sum / count : 0; + } + + private static String getTileName( final SpimData2 spimData, final ViewId viewId ) + { + final BasicViewDescription< ? > vd = spimData.getSequenceDescription().getViewDescriptions().get( viewId ); + if ( vd != null ) + return Group.pvid( viewId ); + else + return viewId.toString(); + } +} diff --git a/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/BlendedExtendedMirroredRandomAccesible2.java b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/BlendedExtendedMirroredRandomAccesible2.java new file mode 100644 index 000000000..ec8bfebeb --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/BlendedExtendedMirroredRandomAccesible2.java @@ -0,0 +1,246 @@ +/*- + * #%L + * Multiview stitching of large datasets. + * %% + * Copyright (C) 2016 - 2025 Big Stitcher developers. + * %% + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program. If not, see + * . + * #L% + */ +package net.preibisch.mvrecon.process.phasecorrelation; + +import java.io.File; + +import ij.ImageJ; +import net.imglib2.Cursor; +import net.imglib2.FinalInterval; +import net.imglib2.Interval; +import net.imglib2.Localizable; +import net.imglib2.Point; +import net.imglib2.RandomAccess; +import net.imglib2.RandomAccessible; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.RealRandomAccess; +import net.imglib2.Sampler; +import net.preibisch.mvrecon.process.phasecorrelation.deprecated.Blending; +import net.imglib2.img.Img; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.img.display.imagej.ImageJFunctions; +import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.util.Intervals; +import net.imglib2.view.Views; +import net.preibisch.mvrecon.process.fusion.transformed.weights.BlendingRealRandomAccessible; + + +public class BlendedExtendedMirroredRandomAccesible2 >implements RandomAccessible { + + private RandomAccessibleInterval img; + BlendingRealRandomAccessible blending; + private int numDimensions; + private FinalInterval extDims; + + public BlendedExtendedMirroredRandomAccesible2(RandomAccessibleInterval img, int[] border) { + this.img = img; + this.numDimensions = img.numDimensions(); + + float[] blendingBorder = new float[numDimensions]; + float[] border2 = new float[numDimensions]; + + this.extDims = new FinalInterval(img); + for (int i = 0; i < numDimensions; i++) + { + extDims = Intervals.expand(extDims, border[i], i); + blendingBorder[i] = border[i]; + border2[i] = 0.0f; + } + + this.blending = new BlendingRealRandomAccessible(extDims, border2, blendingBorder); + } + + + @Override + public int numDimensions() { + return numDimensions; + } + + public FinalInterval getExtInterval() { + return extDims; + } + + @Override + public T getType() + { + return img.getType(); + } + + /** + * TODO: For efficiency reasons we should implement it as a RandomAccess that actually updates the underlying + * imgRA for every move. This way, the outofbounds can work very efficiently when it is iterated through Views.iterable().cursor() + */ + public class BlendedRandomAccess extends Point implements RandomAccess + { + public BlendedRandomAccess() { + super(img.numDimensions()); + } + + RandomAccess imgRA = Views.extendMirrorSingle(img).randomAccess(); + RealRandomAccess blendRA = blending.realRandomAccess(); + T val = imgRA.get().createVariable(); + + @Override + public T get() { + val.setReal(imgRA.get().getRealFloat() * blendRA.get().getRealFloat()); + return val; + } + + @Override + public void fwd(int d) { + super.fwd(d); + imgRA.fwd(d); + blendRA.fwd(d); + } + + @Override + public void bck(int d) { + super.bck(d); + imgRA.bck(d); + blendRA.bck(d); + } + + @Override + public void move(int distance, int d) { + super.move(distance, d); + imgRA.move(distance, d); + blendRA.move(distance, d); + } + + @Override + public void move(long distance, int d) { + super.move(distance, d); + imgRA.move(distance, d); + blendRA.move(distance, d); + } + + @Override + public void move(Localizable localizable) { + super.move(localizable); + imgRA.move(localizable); + blendRA.move(localizable); + } + + @Override + public void move(int[] distance) { + super.move(distance); + imgRA.move(distance); + blendRA.move(distance); + } + + @Override + public void move(long[] distance) { + super.move(distance); + imgRA.move(distance); + blendRA.move(distance); + } + + @Override + public void setPosition(Localizable localizable) { + super.setPosition(localizable); + imgRA.setPosition(localizable); + blendRA.setPosition(localizable); + } + + @Override + public void setPosition(int[] position) { + super.setPosition(position); + imgRA.setPosition(position); + blendRA.setPosition(position); + } + + @Override + public void setPosition(long[] position) { + super.setPosition(position); + imgRA.setPosition(position); + blendRA.setPosition(position); + } + + @Override + public void setPosition(int position, int d) { + super.setPosition(position, d); + imgRA.setPosition(position, d); + blendRA.setPosition(position, d); + } + + @Override + public void setPosition(long position, int d) { + super.setPosition(position, d); + imgRA.setPosition(position, d); + blendRA.setPosition(position, d); + } + + @Override + public RandomAccess copy() { + BlendedRandomAccess a = new BlendedRandomAccess(); + a.move(this); + return a; + } + + } + + @Override + public RandomAccess randomAccess() { + return new BlendedRandomAccess(); + } + + @Override + public RandomAccess randomAccess(Interval interval) { + return randomAccess(); + } + + + public static void main(String[] args) { + + new ImageJ(); + Img img1 = ImgLib2Util.openAs32Bit(new File("src/main/resources/img1.tif")); + long[] dims = new long[img1.numDimensions()]; + + BlendedExtendedMirroredRandomAccesible2 ext = new BlendedExtendedMirroredRandomAccesible2(img1, new int[]{100, 100, 10}); + + ext.getExtInterval().dimensions(dims); + Img img2 = ArrayImgs.floats(dims); + + // TODO: For efficiency reasons we should now also iterate the BlendedExtendedMirroredRandomAccesible and not the image + long start = System.currentTimeMillis(); + + Cursor c = img2.cursor(); + for (FloatType e : Views.iterable(Views.interval(ext, ext.getExtInterval()))) + { + c.fwd(); + c.get().set(e); + } + + long end = System.currentTimeMillis(); + System.out.println(end-start); + + ImageJFunctions.show(img2); + + //RandomAccessibleInterval img3 = Views.interval(ext, ext.getExtInterval()); + //ImageJFunctions.show(img3); + + + + + } +} diff --git a/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/FourNeighborhoodExtrema.java b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/FourNeighborhoodExtrema.java new file mode 100644 index 000000000..6c3083157 --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/FourNeighborhoodExtrema.java @@ -0,0 +1,318 @@ +/*- + * #%L + * Multiview stitching of large datasets. + * %% + * Copyright (C) 2016 - 2025 Big Stitcher developers. + * %% + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program. If not, see + * . + * #L% + */ +package net.preibisch.mvrecon.process.phasecorrelation; + +import java.util.List; +import java.util.ArrayList; +import java.util.Comparator; +import java.util.concurrent.Callable; +import java.util.concurrent.ExecutionException; +import java.util.concurrent.ExecutorService; +import java.util.concurrent.Future; + +import net.imglib2.Cursor; +import net.imglib2.FinalInterval; +import net.imglib2.Interval; +import net.imglib2.Localizable; +import net.imglib2.Point; +import net.imglib2.RandomAccess; +import net.imglib2.RandomAccessible; +import net.imglib2.type.numeric.RealType; +import net.imglib2.util.Pair; +import net.imglib2.util.Util; +import net.imglib2.util.ValuePair; +import net.imglib2.view.Views; + +public class FourNeighborhoodExtrema +{ + + /** + * merge the pre-sorted lists in lists, keep at most the maxN biggest values from any list in the result + * @param lists lists to merge + * @param maxN maximum size of output list + * @param compare comparator + * @param list content type + * @return merged list with size {@literal < maxN} + */ + public static ArrayList merge(List> lists, final int maxN, Comparator compare){ + ArrayList res = new ArrayList(); + int[] idxs = new int[lists.size()]; + boolean[] hasMore = new boolean[lists.size()]; + boolean allDone = true; + + for (int i = 0; i < lists.size(); i++){ + hasMore[i] = lists.get(i).size() > 0; + allDone &= !hasMore[i]; + } + + while(!allDone && res.size() < maxN ){ + + int activeLists = 0; + int maxList = 0; + for (int i = 0; i= 0){ + maxList = i; + } + } + + res.add(lists.get(maxList).get(idxs[maxList])); + idxs[maxList]++; + if (idxs[maxList] >= lists.get(maxList).size()){ + hasMore[maxList] = false; + activeLists--; + } + + allDone = activeLists == 0; + + } + + return res; + } + + /** + * split the given Interval into nSplits intervals along the largest dimension + * @param interval input interval + * @param nSplits how may splits + * @return list of intervals input was split into + */ + public static List splitAlongLargestDimension(Interval interval, long nSplits){ + + List res = new ArrayList(); + + long[] min = new long[interval.numDimensions()]; + long[] max = new long[interval.numDimensions()]; + interval.min(min); + interval.max(max); + + int splitDim = 0; + for (int i = 0; i< interval.numDimensions(); i++){ + if (interval.dimension(i) > interval.dimension(splitDim)) splitDim = i; + } + + // there could be more splits than actual dimension entries + nSplits = Math.min( nSplits, interval.dimension(splitDim) ); + + long chunkSize = interval.dimension(splitDim) / nSplits; + long maxSplitDim = max[splitDim]; + + for (int i = 0; i > ArrayList< Pair< Localizable, Double > > findMaxMT( final RandomAccessible< T > img, final Interval region, final int maxN , ExecutorService service){ + + + int nTasks = Runtime.getRuntime().availableProcessors() * 4; + List intervals = splitAlongLargestDimension(region, nTasks); + List >>> futures = new ArrayList>>>(); + + for (final Interval i : intervals){ + futures.add(service.submit(new Callable >>() { + + @Override + public ArrayList> call() throws Exception { + return findMax(img, i, maxN); + } + })); + } + + List >> toMerge = new ArrayList>>(); + + for (Future >> f : futures){ + try { + toMerge.add(f.get()); + } catch (InterruptedException e) { + // TODO Auto-generated catch block + e.printStackTrace(); + } catch (ExecutionException e) { + // TODO Auto-generated catch block + e.printStackTrace(); + } + } + + ArrayList< Pair< Localizable, Double > > res = merge(toMerge, maxN, new Comparator>() { + + @Override + public int compare(Pair o1, Pair o2) { + return (int) Math.signum(o1.getB() - o2.getB()); + } + }); + + return res; + } + + public static < T extends RealType< T > > ArrayList< Pair< Localizable, Double > > findMax( final RandomAccessible< T > img, final Interval region, final int maxN ) + { + final Cursor< T > c = Views.iterable( Views.interval( img, region ) ).localizingCursor(); + final RandomAccess< T > r = img.randomAccess(); + final int n = img.numDimensions(); + + final ArrayList< Pair< Localizable, Double > > list = new ArrayList< Pair< Localizable, Double > >(); + + for ( int i = 0; i < maxN; ++i ) + list.add( new ValuePair< Localizable, Double >( null, -Double.MAX_VALUE ) ); + +A: while ( c.hasNext() ) + { + final double type = c.next().getRealDouble(); + r.setPosition( c ); + + for ( int d = 0; d < n; ++d ) + { + r.fwd( d ); + if ( type < r.get().getRealDouble() ) + continue A; + + r.bck( d ); + r.bck( d ); + + if ( type < r.get().getRealDouble() ) + continue A; + + r.fwd( d ); + } + + + for ( int i = maxN - 1; i >= 0; --i ) + { + if ( type < list.get( i ).getB() ) + { + if ( i == maxN - 1 ) + { + continue A; + } + else + { + list.add( i + 1, new ValuePair< Localizable, Double >( new Point( c ), type ) ); + list.remove( maxN ); + continue A; + } + } + } + + list.add( 0, new ValuePair< Localizable, Double >( new Point( c ), type ) ); + list.remove( maxN ); + } + + // remove all null elements + for ( int i = maxN -1; i >= 0; --i ) + if ( list.get( i ).getA() == null ) + list.remove( i ); + + return list; + } + + public static void main( String[] args ) + { + int maxN = 3; + ArrayList< Double > list = new ArrayList< Double >(); + + list.add( 5.0 ); + list.add( 2.0 ); + list.add( 0.0 ); + + /* + double type = 10; + + for ( int i = maxN - 1; i >= 0; --i ) + { + if ( type < list.get( i ) ) + { + if ( i == maxN - 1 ) + { + printList( list ); + System.exit( 0 ); + } + else + { + list.add( i + 1, type ); + list.remove( maxN ); + printList( list ); + System.exit( 0 ); + } + } + } + + list.add( 0, type ); + list.remove( maxN ); + printList( list ); + + */ + ArrayList< Double > list1 = new ArrayList< Double >(); + list1.add( 5.0 ); + list1.add( 2.0 ); + list1.add( 0.0 ); + ArrayList< Double > list2 = new ArrayList< Double >(); + list2.add( 8.0 ); + list2.add( 3.0 ); + list2.add( 1.0 ); + + ArrayList> both = new ArrayList>(); + both.add(list1); + both.add(list2); + + ArrayList res = merge(both, 12, new Comparator() { + + @Override + public int compare(Double o1, Double o2) { + return Double.compare(o1, o2); + } + }); + + printList(res); + + + Interval toSplit = new FinalInterval(new long[] {20, 40, 70}); + + List splits = splitAlongLargestDimension(toSplit, 5); + + + + + + } + + public static void printList( ArrayList< Double > list ) + { + for ( double d : list ) + System.out.println( d ); + } +} diff --git a/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/ImgLib2Util.java b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/ImgLib2Util.java new file mode 100644 index 000000000..494cff997 --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/ImgLib2Util.java @@ -0,0 +1,137 @@ +/*- + * #%L + * Multiview stitching of large datasets. + * %% + * Copyright (C) 2016 - 2025 Big Stitcher developers. + * %% + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program. If not, see + * . + * #L% + */ +package net.preibisch.mvrecon.process.phasecorrelation; + + +import ij.ImageJ; +import ij.ImagePlus; +import ij.io.Opener; +import ij.process.ImageProcessor; + +import java.io.File; +import java.util.ArrayList; + +import edu.mines.jtk.dsp.FftComplex; +import net.imglib2.Cursor; +import net.imglib2.IterableInterval; +import net.imglib2.RandomAccess; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.algorithm.fft2.FFTMethods; +import net.imglib2.img.Img; +import net.imglib2.img.ImgFactory; +import net.imglib2.img.array.ArrayImgFactory; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.img.display.imagej.ImageJFunctions; +import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.view.Views; + +public class ImgLib2Util +{ + public static Img< FloatType > openAs32Bit( final File file ) + { + return openAs32Bit( file, new ArrayImgFactory< FloatType >() ); + } + + public static Img< FloatType > openAs32Bit( final File file, final ImgFactory< FloatType > factory ) + { + if ( !file.exists() ) + throw new RuntimeException( "File '" + file.getAbsolutePath() + "' does not exisit." ); + + final ImagePlus imp = new Opener().openImage( file.getAbsolutePath() ); + + if ( imp == null ) + throw new RuntimeException( "File '" + file.getAbsolutePath() + "' coult not be opened." ); + + final Img< FloatType > img; + + if ( imp.getStack().getSize() == 1 ) + { + // 2d + img = factory.create( new int[]{ imp.getWidth(), imp.getHeight() }, new FloatType() ); + final ImageProcessor ip = imp.getProcessor(); + + final Cursor< FloatType > c = img.localizingCursor(); + + while ( c.hasNext() ) + { + c.fwd(); + + final int x = c.getIntPosition( 0 ); + final int y = c.getIntPosition( 1 ); + + c.get().set( ip.getf( x, y ) ); + } + + } + else + { + // >2d + img = factory.create( new int[]{ imp.getWidth(), imp.getHeight(), imp.getStack().getSize() }, new FloatType() ); + + final Cursor< FloatType > c = img.localizingCursor(); + + // for efficiency reasons + final ArrayList< ImageProcessor > ips = new ArrayList< ImageProcessor >(); + + for ( int z = 0; z < imp.getStack().getSize(); ++z ) + ips.add( imp.getStack().getProcessor( z + 1 ) ); + + while ( c.hasNext() ) + { + c.fwd(); + + final int x = c.getIntPosition( 0 ); + final int y = c.getIntPosition( 1 ); + final int z = c.getIntPosition( 2 ); + + c.get().set( ips.get( z ).getf( x, y ) ); + } + } + + return img; + } + + + public static , S extends RealType> void copyRealImage(IterableInterval source, RandomAccessibleInterval dest) { + RandomAccess destRA = dest.randomAccess(); + Cursor srcC = source.cursor(); + + + while (srcC.hasNext()){ + srcC.fwd(); + destRA.setPosition(srcC); + destRA.get().setReal(srcC.get().getRealDouble()); + } + } + + + public static void main( String[] args ) + { + new ImageJ(); + + final Img< FloatType > img = openAs32Bit( new File( "src/main/resources/mri-stack.tif" ) ); + //final Img< FloatType > img = openAs32Bit( new File( "src/main/resources/bridge.png" ) ); + + ImageJFunctions.show( img ); + } +} diff --git a/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/PhaseCorrelation2.java b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/PhaseCorrelation2.java new file mode 100644 index 000000000..bb1b81123 --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/PhaseCorrelation2.java @@ -0,0 +1,313 @@ +/*- + * #%L + * Multiview stitching of large datasets. + * %% + * Copyright (C) 2016 - 2025 Big Stitcher developers. + * %% + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program. If not, see + * . + * #L% + */ +package net.preibisch.mvrecon.process.phasecorrelation; + +import java.io.File; +import java.util.Arrays; +import java.util.Collections; +import java.util.List; +import java.util.concurrent.ExecutorService; +import java.util.concurrent.Executors; + +import ij.ImageJ; +import net.imglib2.Dimensions; +import net.imglib2.FinalInterval; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.algorithm.fft2.FFT; +import net.imglib2.algorithm.fft2.FFTMethods; +import net.imglib2.exception.IncompatibleTypeException; +import net.imglib2.img.Img; +import net.imglib2.img.ImgFactory; +import net.imglib2.img.array.ArrayImgFactory; +import net.imglib2.img.display.imagej.ImageJFunctions; +import net.imglib2.type.NativeType; +import net.imglib2.type.numeric.ComplexType; +import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.complex.ComplexFloatType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.view.Views; + +public class PhaseCorrelation2 { + + public static boolean debug = false; + + /* + * calculate the phase correlation of fft1 and fft2, save result to res + * fft1 and fft2 will be altered by the function + * @param fft1 fft of first image + * @param fft2 fft of second image + * @param pcm + * @param service + * @param + * @param + * @param + */ + public static , S extends ComplexType, R extends RealType> void calculatePCMInPlace( + RandomAccessibleInterval fft1, RandomAccessibleInterval fft2, RandomAccessibleInterval pcm, ExecutorService service) + { + calculatePCM( fft1, fft1, fft2, fft2, pcm , service); + } + + /* + * calculate the phase correlation of fft1 and fft2, save result to res + * fft1 and fft2 will NOT be altered by the function + * @param fft1 + * @param fft1Copy - a temporary image same size as fft1 and fft2 + * @param fft2 + * @param fft2Copy - a temporary image same size as fft1 and fft2 + * @param pcm + */ + public static , S extends ComplexType, R extends RealType> void calculatePCM( + RandomAccessibleInterval fft1, RandomAccessibleInterval fft1Copy, RandomAccessibleInterval fft2, RandomAccessibleInterval fft2Copy, RandomAccessibleInterval pcm, + ExecutorService service) + { + // TODO: multithreaded & check for cursor vs randomaccess + + // normalize, save to copies + PhaseCorrelation2Util.normalizeInterval(fft1, fft1Copy, service); + PhaseCorrelation2Util.normalizeInterval(fft2, fft2Copy, service); + // conjugate + PhaseCorrelation2Util.complexConjInterval(fft2Copy, fft2Copy, service); + // in-place multiplication + PhaseCorrelation2Util.multiplyComplexIntervals(fft1Copy, fft2Copy, fft1Copy, service); + FFT.complexToReal(fft1Copy, pcm, service); + } + + + /* + * calculate phase correlation of fft1 and fft2, return result in a new Img + * @param fft1 + * @param fft2 + * @return + */ + public static , S extends ComplexType, R extends RealType> RandomAccessibleInterval calculatePCMInPlace( + RandomAccessibleInterval fft1, RandomAccessibleInterval fft2, ImgFactory factory, R type, ExecutorService service){ + + long[] paddedDimensions = new long[fft1.numDimensions()]; + long[] realSize = new long[fft1.numDimensions()]; + + FFTMethods.dimensionsComplexToRealFast(fft1, paddedDimensions, realSize); + RandomAccessibleInterval res = factory.create(realSize, type); + + calculatePCMInPlace(fft1, fft2, res, service); + + return res; + } + + /* + * calculate phase correlation of fft1 and fft2, doing the calculations in copies of the ffts, return result in a new Img + * @param fft1 + * @param fft2 + * @return + */ + public static & NativeType, S extends ComplexType & NativeType , R extends RealType> RandomAccessibleInterval calculatePCM( + RandomAccessibleInterval fft1, RandomAccessibleInterval fft2, ImgFactory factory, R type, ExecutorService service){ + + long[] paddedDimensions = new long[fft1.numDimensions()]; + long[] realSize = new long[fft1.numDimensions()]; + + FFTMethods.dimensionsComplexToRealFast(fft1, paddedDimensions, realSize); + RandomAccessibleInterval res = factory.create(realSize, type); + + final T typeT = Views.iterable(fft1).firstElement().createVariable(); + final S typeS = Views.iterable(fft2).firstElement().createVariable(); + RandomAccessibleInterval< T > fft1Copy; + RandomAccessibleInterval< S > fft2Copy; + + try + { + fft1Copy = factory.imgFactory( typeT ).create(fft1, typeT ); + fft2Copy = factory.imgFactory( typeS ).create(fft2, typeS ); + } + catch ( IncompatibleTypeException e ) + { + throw new RuntimeException( "Cannot instantiate Img for type " + typeS.getClass().getSimpleName() + " or " + typeT.getClass().getSimpleName() ); + } + + + calculatePCM(fft1, fft1Copy, fft2, fft2Copy, res, service); + + return res; + } + + /* + * calculate and return the phase correlation matrix of two images + * @param img1 + * @param img2 + * @return + */ + public static , S extends RealType, R extends RealType, C extends ComplexType> RandomAccessibleInterval calculatePCM( + RandomAccessibleInterval img1, RandomAccessibleInterval img2, int[] extension, + ImgFactory factory, R type, ImgFactory fftFactory, C fftType, ExecutorService service){ + + + // TODO: Extension absolute per dimension in pixels, i.e. int[] extension + // TODO: not bigger than the image dimension because the second mirroring is identical to the image + + Dimensions extSize = PhaseCorrelation2Util.getExtendedSize(img1, img2, extension); + long[] paddedDimensions = new long[extSize.numDimensions()]; + long[] fftSize = new long[extSize.numDimensions()]; + FFTMethods.dimensionsRealToComplexFast(extSize, paddedDimensions, fftSize); + + RandomAccessibleInterval fft1 = fftFactory.create(fftSize, fftType); + RandomAccessibleInterval fft2 = fftFactory.create(fftSize, fftType); + + FFT.realToComplex(Views.interval(PhaseCorrelation2Util.extendImageByFactor(img1, extension), + FFTMethods.paddingIntervalCentered(img1, new FinalInterval(paddedDimensions))), fft1, service); + FFT.realToComplex(Views.interval(PhaseCorrelation2Util.extendImageByFactor(img2, extension), + FFTMethods.paddingIntervalCentered(img2, new FinalInterval(paddedDimensions))), fft2, service); + + RandomAccessibleInterval pcm = calculatePCMInPlace(fft1, fft2, factory, type, service); + return pcm; + + } + + /* + * calculate PCM with default extension + * @param img1 + * @param img2 + * @return + */ + public static , S extends RealType, R extends RealType, C extends ComplexType> RandomAccessibleInterval calculatePCM( + RandomAccessibleInterval img1, RandomAccessibleInterval img2, ImgFactory factory, R type, + ImgFactory fftFactory, C fftType, ExecutorService service) { + + int [] extension = new int[img1.numDimensions()]; + Arrays.fill(extension, 10); + return calculatePCM(img1, img2, extension, factory, type, fftFactory, fftType, service); + } + + /** + * calculate the shift between two images from the phase correlation matrix + * @param pcm the phase correlation matrix of img1 and img2 + * @param img1 source image 1 + * @param img2 source image 2 + * @param nHighestPeaks the number of peaks in pcm to check via cross. corr. + * @param minOverlap minimal overlap (in pixels) + * @param subpixelAccuracy whether to do subpixel shift peak localization or not + * @param interpolateSubpixel whether to interpolate the subpixel shift in cross. corr. + * @param service thread pool + * @param PCM pixel type + * @param image 1 pixel type + * @param image 2 pixel type + * @return best (highest c.c.) shift peak + */ + public static , S extends RealType, R extends RealType> PhaseCorrelationPeak2 getShift( + RandomAccessibleInterval pcm, RandomAccessibleInterval img1, RandomAccessibleInterval img2, int nHighestPeaks, + long minOverlap, boolean subpixelAccuracy, boolean interpolateSubpixel, ExecutorService service) + { + if ( debug ) + System.out.println( "PCM" ); + + List peaks = PhaseCorrelation2Util.getPCMMaxima(pcm, service, nHighestPeaks, subpixelAccuracy); + //peaks = PhaseCorrelation2Util.getHighestPCMMaxima(peaks, nHighestPeaks); + + if ( debug ) + System.out.println( "expand" ); + + PhaseCorrelation2Util.expandPeakListToPossibleShifts(peaks, pcm, img1, img2); + + if ( debug ) + System.out.print( "cross " ); + + long t = System.currentTimeMillis(); + + PhaseCorrelation2Util.calculateCrossCorrParallel(peaks, img1, img2, minOverlap, service, interpolateSubpixel); + + if ( debug ) + { + System.out.println( (System.currentTimeMillis() - t) ); + System.out.println( "sort" ); + } + + Collections.sort(peaks, Collections.reverseOrder(new PhaseCorrelationPeak2.ComparatorByCrossCorrelation())); + + if ( debug ) + System.out.println( "done" ); + + if (peaks.size() > 0) + return peaks.get(0); + else + return null; + } + + /** + * get shift, do not interpolate subpixel offset for cross correlation + * @param pcm the phase correlation matrix of img1 and img2 + * @param img1 source image 1 + * @param img2 source image 2 + * @param nHighestPeaks the number of peaks in pcm to check via cross. corr. + * @param minOverlap minimal overlap (in pixels) + * @param subpixelAccuracy whether to do subpixel shift peak localization or not + * @param service thread pool + * @param PCM pixel type + * @param image 1 pixel type + * @param image 2 pixel type + * @return best (highest c.c.) shift peak + */ + public static , S extends RealType, R extends RealType> PhaseCorrelationPeak2 getShift( + RandomAccessibleInterval pcm, RandomAccessibleInterval img1, RandomAccessibleInterval img2, int nHighestPeaks, + long minOverlap, boolean subpixelAccuracy, ExecutorService service) + { + return getShift( pcm, img1, img2, nHighestPeaks, minOverlap, subpixelAccuracy, false, service ); + } + + /** + * calculate the sift with default parameters (5 highest pcm peaks are considered, no minimum overlap, temporary thread pool, + * no subpixel interpolation) + * @param pcm the phase correlation matrix of img1 and img2 + * @param img1 source image 1 + * @param img2 source image 2 + * @param PCM pixel type + * @param image 1 pixel type + * @param image 2 pixel type + * @return best (highest c.c.) shift peak + */ + public static , S extends RealType, R extends RealType> PhaseCorrelationPeak2 getShift( + RandomAccessibleInterval pcm, RandomAccessibleInterval img1, RandomAccessibleInterval img2) + { + ExecutorService service = Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors()); + PhaseCorrelationPeak2 res = getShift(pcm, img1, img2, 5, 0, true, false, service); + service.shutdown(); + return res; + } + + public static void main(String[] args) { + + new ImageJ(); + + Img img1 = ImgLib2Util.openAs32Bit(new File("src/main/resources/img1singleplane.tif")); + Img img2 = ImgLib2Util.openAs32Bit(new File("src/main/resources/img2singleplane.tif")); + + ExecutorService service = Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors()); + + RandomAccessibleInterval pcm = calculatePCM(img1, img2, new ArrayImgFactory(), new FloatType(), + new ArrayImgFactory(), new ComplexFloatType(), service ); + + PhaseCorrelationPeak2 shiftPeak = getShift(pcm, img1, img2); + + RandomAccessibleInterval res = PhaseCorrelation2Util.dummyFuse(img1, img2, shiftPeak,service); + + ImageJFunctions.show(res); + } + +} diff --git a/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/PhaseCorrelation2Util.java b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/PhaseCorrelation2Util.java new file mode 100644 index 000000000..8a6c62e4c --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/PhaseCorrelation2Util.java @@ -0,0 +1,882 @@ +/*- + * #%L + * Multiview stitching of large datasets. + * %% + * Copyright (C) 2016 - 2025 Big Stitcher developers. + * %% + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program. If not, see + * . + * #L% + */ +package net.preibisch.mvrecon.process.phasecorrelation; + +import java.io.File; +import java.util.ArrayList; +import java.util.Collections; +import java.util.List; +import java.util.Vector; +import java.util.concurrent.ExecutionException; +import java.util.concurrent.ExecutorService; +import java.util.concurrent.Executors; +import java.util.concurrent.Future; +import java.util.concurrent.atomic.AtomicInteger; + +import net.imglib2.Cursor; +import net.imglib2.Dimensions; +import net.imglib2.FinalDimensions; +import net.imglib2.FinalInterval; +import net.imglib2.Interval; +import net.imglib2.IterableInterval; +import net.imglib2.Localizable; +import net.imglib2.Point; +import net.imglib2.RandomAccess; +import net.imglib2.RandomAccessible; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.RealPoint; +import net.imglib2.img.Img; +import net.imglib2.img.array.ArrayImgFactory; +import net.imglib2.realtransform.AffineTransform3D; +import net.imglib2.realtransform.Translation3D; +import net.imglib2.type.numeric.ComplexType; +import net.imglib2.type.numeric.RealType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.util.BenchmarkHelper; +import net.imglib2.util.Pair; +import net.imglib2.util.Util; +import net.imglib2.util.ValuePair; +import net.imglib2.view.Views; +import net.preibisch.mvrecon.process.fusion.FusionTools; +import net.preibisch.mvrecon.process.fusion.ImagePortion; + + + +public class PhaseCorrelation2Util { + + /* + * copy source to dest. they do not have to be of the same size, but source must fit in dest + * @param source + * @param dest + */ + public static , S extends RealType> void copyRealImage(final IterableInterval source, final RandomAccessibleInterval dest, ExecutorService service) { + + final Vector portions = FusionTools.divideIntoPortions( source.size() ); + ArrayList> futures = new ArrayList>(); + final AtomicInteger ai = new AtomicInteger(-1); + + for (int i = 0; i < portions.size(); i++){ + futures.add(service.submit(new Runnable() { + @Override + public void run() { + RandomAccess destRA = dest.randomAccess(); + Cursor srcC = source.localizingCursor(); + + ImagePortion ip = portions.get(ai.incrementAndGet()); + + long loopSize = ip.getLoopSize(); + + srcC.jumpFwd(ip.getStartPosition()); + + for (long l = 0; l < loopSize; l++){ + srcC.fwd(); + destRA.setPosition(srcC); + destRA.get().setReal(srcC.get().getRealDouble()); + } + + } + })); + } + + for (Future f : futures){ + try { + f.get(); + } catch (InterruptedException e) { + // TODO Auto-generated catch block + e.printStackTrace(); + } catch (ExecutionException e) { + // TODO Auto-generated catch block + e.printStackTrace(); + } + } + + } + + /** + * calculate the size difference of two Dimensions objects (dim2-dim1) + * @param dim1 first Dimensions + * @param dim2 second Dimensions + * @return int array of difference + */ + public static int[] getSizeDifference(Dimensions dim1, Dimensions dim2) { + int[] diff = new int[dim1.numDimensions()]; + for (int i = 0; i < dim1.numDimensions(); i++){ + diff[i] = (int) (dim2.dimension(i) - dim1.dimension(i)); + } + return diff; + } + + /** + * calculate the size of an extended image big enough to hold dim1 and dim2 + * with each dimension also enlarged by extension pixels on each side (but at most by the original image size) + * @param dim1 first Dimensions + * @param dim2 second Dimensions + * @param extension: number of pixels to add at each side in each dimension + * @return extended dimensions + */ + public static FinalDimensions getExtendedSize(Dimensions dim1, Dimensions dim2, int [] extension) { + long[] extDims = new long[dim1.numDimensions()]; + for (int i = 0; i dim2.dimension(i) ? dim1.dimension(i) : dim2.dimension(i); + long extBothSides = extDims[i] < extension[i] ? extDims[i] * 2 : extension[i] * 2; + extDims[i] += extBothSides; + } + return new FinalDimensions(extDims); + } + + /* + * return a BlendedExtendedMirroredRandomAccesible of img extended extension pixels on each side (but at most by the original image size) + * @param img + * @param extension: number of blending pixels to add at each side in each dimension + * @return + */ + public static > RandomAccessible extendImageByFactor(RandomAccessibleInterval img, int [] extension) + { + int[] extEachSide = new int[img.numDimensions()]; + for (int i = 0; i (img, extEachSide); + } + + /* + * returns the extension at each side if an image is enlarged by a factor of extensionFactor at each side + * @param dims + * @param extensionFactor + * @return + */ + public static int[] extensionByFactor(Dimensions dims, double extensionFactor){ + int[] res = new int[dims.numDimensions()]; + for (int i = 0; i< dims.numDimensions(); i++){ + res[i] = (int) (dims.dimension(i)*extensionFactor); + } + return res; + } + + + /* + * return a BlendedExtendedMirroredRandomAccesible of img extended to extDims + * @param img + * @param extDims + * @return + */ + public static > RandomAccessible extendImageToSize(RandomAccessibleInterval img, Dimensions extDims) + { + int[] extEachSide = getSizeDifference(img, extDims); + for (int i = 0; i< img.numDimensions(); i++){ + extEachSide[i] /= 2; + } + return new BlendedExtendedMirroredRandomAccesible2(img, extEachSide); + } + + + public static , S extends RealType> void calculateCrossCorrParallel( + List peaks, final RandomAccessibleInterval img1, final RandomAccessibleInterval img2, + final long minOverlapPx, ExecutorService service) + { + calculateCrossCorrParallel( peaks, img1, img2, minOverlapPx, service, false ); + } + + /* + * calculate the crosscorrelation of img1 and img2 for all shifts represented by a PhasecorrelationPeak List in parallel using a specified + * ExecutorService. service remains functional after the call + * @param peaks + * @param img1 + * @param img2 + * @param minOverlapPx minimal number of overlapping pixels in each Dimension, may be null to indicate no minimum + * @param service + */ + public static , S extends RealType> void calculateCrossCorrParallel( + List peaks, final RandomAccessibleInterval img1, final RandomAccessibleInterval img2, + final long minOverlapPx, ExecutorService service, boolean interpolateSubpixel) + { + List> futures = new ArrayList>(); + + for (final PhaseCorrelationPeak2 p : peaks){ + futures.add(service.submit(new Runnable() { + @Override + public void run() { + p.calculateCrossCorr(img1, img2, minOverlapPx, interpolateSubpixel); + } + })); + } + + for (Future f: futures){ + try { + f.get(); + } catch (InterruptedException e) { + // TODO Auto-generated catch block + e.printStackTrace(); + } catch (ExecutionException e) { + // TODO Auto-generated catch block + e.printStackTrace(); + } + } + } + + /* + * find local maxima in PCM + * @param pcm + * @param service + * @param maxN + * @return + */ + public static > List getPCMMaxima(RandomAccessibleInterval pcm, ExecutorService service, int maxN, boolean subpixelAccuracy){ + + List res = new ArrayList(); + + ArrayList> maxima = FourNeighborhoodExtrema.findMaxMT(Views.extendPeriodic(pcm), pcm, maxN, service); + //ArrayList> maxima = FourNeighborhoodExtrema.findMax(Views.extendPeriodic(pcm), pcm, maxN); + + + + for (Pair p: maxima){ + PhaseCorrelationPeak2 pcp = new PhaseCorrelationPeak2(p.getA(), p.getB()); + if (subpixelAccuracy) + pcp.calculateSubpixelLocalization(pcm); + + res.add(pcp); + } + return res; + } + + /* + * find maxima in PCM, use a temporary thread pool for calculation + * @param pcm + * @param nMax + * @return + */ + public static > List getPCMMaxima(RandomAccessibleInterval pcm, int nMax, boolean subpixelAccuracy){ + ExecutorService tExecService = Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors()); + List res = getPCMMaxima(pcm, tExecService, nMax, subpixelAccuracy); + tExecService.shutdown(); + return res; + } + + /* + * sort PCM Peaks by phaseCorrelation and return a new list containing just the nToKeep highest peaks + * @param rawPeaks + * @param nToKeep + * @return + */ + @Deprecated + public static List getHighestPCMMaxima(List rawPeaks, long nToKeep){ + Collections.sort(rawPeaks, Collections.reverseOrder( new PhaseCorrelationPeak2.ComparatorByPhaseCorrelation())); + List res = new ArrayList(); + for (int i = 0; i < nToKeep; i++){ + res.add(new PhaseCorrelationPeak2(rawPeaks.get(i))); + } + return res; + } + /* + * expand a list of PCM maxima to to a list containing all possible shifts corresponding to these maxima + * @param peaks + * @param pcmDims + * @param img1Dims + * @param img2Dims + */ + public static void expandPeakListToPossibleShifts(List peaks, + Dimensions pcmDims, Dimensions img1Dims, Dimensions img2Dims) + { + List res = new ArrayList(); + for (PhaseCorrelationPeak2 p : peaks){ + res.addAll(expandPeakToPossibleShifts(p, pcmDims, img1Dims, img2Dims)); + } + peaks.clear(); + peaks.addAll(res); + } + + /* + * expand a single maximum in the PCM to a list of possible shifts corresponding to that peak + * an offset due to different images sizes is accounted for + * @param peak + * @param pcmDims + * @param img1Dims + * @param img2Dims + * @return + */ + public static List expandPeakToPossibleShifts( + PhaseCorrelationPeak2 peak, Dimensions pcmDims, Dimensions img1Dims, Dimensions img2Dims) + { + int n = pcmDims.numDimensions(); + double[] subpixelDiff = new double[n]; + + if (peak.getSubpixelPcmLocation() != null) + for (int i = 0; i < n; i++) + subpixelDiff[i] = peak.getSubpixelPcmLocation().getDoublePosition( i ) - peak.getPcmLocation().getIntPosition( i ); + + int[] originalPCMPeakWithOffset = new int[n]; + peak.getPcmLocation().localize(originalPCMPeakWithOffset); + + int[] extensionImg1 = getSizeDifference(img1Dims, pcmDims); + int[] extensionImg2 = getSizeDifference(img2Dims, pcmDims); + int[] offset = new int[pcmDims.numDimensions()]; + for(int i = 0; i < offset.length; i++){ + offset[i] = (extensionImg2[i] - extensionImg1[i] ) / 2; + originalPCMPeakWithOffset[i] += offset[i]; + originalPCMPeakWithOffset[i] %= pcmDims.dimension(i); + } + + List shiftedPeaks = new ArrayList(); + for (int i = 0; i < Math.pow(2, pcmDims.numDimensions()); i++){ + int[] possibleShift = originalPCMPeakWithOffset.clone(); + PhaseCorrelationPeak2 peakWithShift = new PhaseCorrelationPeak2(peak); + for (int d = 0; d < pcmDims.numDimensions(); d++){ + /* + * mirror the shift around the origin in dimension d if (i / 2^d) is even + * --> all possible shifts + */ + if ((i / (int) Math.pow(2, d) % 2) == 0){ + possibleShift[d] = possibleShift[d] < 0 ? possibleShift[d] + (int) pcmDims.dimension(d) : possibleShift[d] - (int) pcmDims.dimension(d); + } + } + peakWithShift.setShift(new Point(possibleShift)); + + if (peakWithShift.getSubpixelPcmLocation() != null) + { + double[] subpixelShift = new double[n]; + for (int j =0; j getOverlapIntervals(Dimensions img1, Dimensions img2, Localizable shift){ + + final int numDimensions = img1.numDimensions(); + final long[] offsetImage1 = new long[ numDimensions ]; + final long[] offsetImage2 = new long[ numDimensions ]; + final long[] maxImage1 = new long[ numDimensions ]; + final long[] maxImage2 = new long[ numDimensions ]; + + long overlapSize; + + for ( int d = 0; d < numDimensions; ++d ) + { + if ( shift.getLongPosition(d) >= 0 ) + { + // two possiblities + // + // shift=start end + // | | + // A: Image 1 ------------------------------ + // Image 2 ---------------------------------- + // + // shift=start end + // | | + // B: Image 1 ------------------------------ + // Image 2 ------------------- + + // they are not overlapping ( this might happen due to fft zeropadding and extension ) + if ( shift.getLongPosition(d) >= img1.dimension( d ) ) + { + return null; + } + + offsetImage1[ d ] = shift.getLongPosition(d); + offsetImage2[ d ] = 0; + overlapSize = Math.min( img1.dimension( d ) - shift.getLongPosition(d), img2.dimension( d ) ); + maxImage1[ d ] = offsetImage1[d] + overlapSize -1; + maxImage2[ d ] = offsetImage2[d] + overlapSize -1; + } + else + { + // two possiblities + // + // shift start end + // | | ` | + // A: Image 1 ------------------------------ + // Image 2 ------------------------------ + // + // shift start end + // | | | + // B: Image 1 ------------ + // Image 2 ------------------- + + // they are not overlapping ( this might happen due to fft zeropadding and extension + if ( shift.getLongPosition(d) <= -img2.dimension( d ) ) + { + return null; + } + + offsetImage1[ d ] = 0; + offsetImage2[ d ] = -shift.getLongPosition(d); + overlapSize = Math.min( img2.dimension( d ) + shift.getLongPosition(d), img1.dimension( d ) ); + maxImage1[ d ] = offsetImage1[d] + overlapSize -1; + maxImage2[ d ] = offsetImage2[d] + overlapSize -1; + } + + } + + FinalInterval img1Interval = new FinalInterval(offsetImage1, maxImage1); + FinalInterval img2Interval = new FinalInterval(offsetImage2, maxImage2); + + Pair res = new ValuePair(img1Interval, img2Interval); + return res; + } + + /* + * multiply complex numbers c1 and c2, set res to the result of multiplication + * @param c1 + * @param c2 + * @param res + */ + public static , S extends ComplexType, T extends ComplexType> void multiplyComplex( + R c1, S c2, T res) + { + double a = c1.getRealDouble(); + double b = c1.getImaginaryDouble(); + double c = c2.getRealDouble(); + double d = c2.getImaginaryDouble(); + res.setReal(a*c - b*d); + res.setImaginary(a*d + b*c); + } + + /* + * pixel-wise multiplication of img1 and img2 + * res is overwritten by the result + * @param img1 + * @param img2 + * @param res + */ + public static , S extends ComplexType, T extends ComplexType> void multiplyComplexIntervals( + final RandomAccessibleInterval img1, final RandomAccessibleInterval img2, final RandomAccessibleInterval res, ExecutorService service) + { + + + final Vector portions = FusionTools.divideIntoPortions (Views.iterable(img1).size() ); + List> futures = new ArrayList>(); + final AtomicInteger ai = new AtomicInteger(-1); + + for (int i = 0; i < portions.size(); i++){ + futures.add(service.submit(new Runnable() { + + @Override + public void run() { + + ImagePortion ip = portions.get(ai.incrementAndGet()); + long loopSize = ip.getLoopSize(); + + if (Views.iterable(img1).iterationOrder().equals(Views.iterable(img2).iterationOrder()) && + Views.iterable(img1).iterationOrder().equals(Views.iterable(res).iterationOrder())){ + Cursor cRes = Views.iterable(res).cursor(); + Cursor cSrc1 = Views.iterable(img1).cursor(); + Cursor cSrc2 = Views.iterable(img2).cursor(); + + cSrc1.jumpFwd(ip.getStartPosition()); + cSrc2.jumpFwd(ip.getStartPosition()); + cRes.jumpFwd(ip.getStartPosition()); + + for (long l = 0; l < loopSize; l++){ + cRes.fwd(); + cSrc1.fwd(); + cSrc2.fwd(); + multiplyComplex(cSrc1.get(), cSrc2.get(), cRes.get()); + } + } + + else { + final RandomAccess ra1 = img1.randomAccess(); + final RandomAccess ra2 = img2.randomAccess(); + final Cursor cRes = Views.iterable(res).localizingCursor(); + + cRes.jumpFwd(ip.getStartPosition()); + + for (long l = 0; l < loopSize; l++){ + cRes.fwd(); + ra1.setPosition(cRes); + ra2.setPosition(cRes); + multiplyComplex(ra1.get(), ra2.get(), cRes.get()); + } + } + + } + })); + + for (Future f : futures){ + try { + f.get(); + } catch (InterruptedException e) { + // TODO Auto-generated catch block + e.printStackTrace(); + } catch (ExecutionException e) { + // TODO Auto-generated catch block + e.printStackTrace(); + } + } + } + + } + + /* + * calculate complex conjugate of c, save result to res + * @param c + * @param res + */ + public static , S extends ComplexType> void complexConj(R c, S res) { + res.setComplexNumber(c.getRealDouble(), - c.getImaginaryDouble()); + } + + /* + * calculate element-wise complex conjugate of img, save result to res + * @param img + * @param res + */ + public static , S extends ComplexType> void complexConjInterval( + final RandomAccessibleInterval img, final RandomAccessibleInterval res, ExecutorService service) + { + + final Vector portions = FusionTools.divideIntoPortions( Views.iterable(img).size() ); + List> futures = new ArrayList>(); + final AtomicInteger ai = new AtomicInteger(-1); + + for (int i = 0; i < portions.size(); i++){ + futures.add(service.submit(new Runnable() { + + @Override + public void run() { + + ImagePortion ip = portions.get(ai.incrementAndGet()); + long loopSize = ip.getLoopSize(); + + if (Views.iterable(img).iterationOrder().equals(Views.iterable(res).iterationOrder())){ + Cursor cRes = Views.iterable(res).cursor(); + Cursor cSrc = Views.iterable(img).cursor(); + + cSrc.jumpFwd(ip.getStartPosition()); + cRes.jumpFwd(ip.getStartPosition()); + + for (long l = 0; l < loopSize; l++){ + cRes.fwd(); + cSrc.fwd(); + complexConj(cSrc.get(), cRes.get()); + } + } + + else { + Cursor cRes = Views.iterable(res).localizingCursor(); + RandomAccess raImg = img.randomAccess(); + + cRes.jumpFwd(ip.getStartPosition()); + + for (long l = 0; l < loopSize; l++){ + cRes.fwd(); + raImg.setPosition(cRes); + complexConj(raImg.get(), cRes.get()); + } + } + + } + })); + + for (Future f : futures){ + try { + f.get(); + } catch (InterruptedException e) { + // TODO Auto-generated catch block + e.printStackTrace(); + } catch (ExecutionException e) { + // TODO Auto-generated catch block + e.printStackTrace(); + } + } + } + + } + + /* + * normalize complex number c1 to length 1, save result to res + * if the length of c1 is less than normalizationThreshold, set res to 0 + * @param c1 + * @param res + * @param normalizationThreshold + */ + public static , S extends ComplexType> void normalize( R c1, S res, double normalizationThreshold) + { + double len = c1.getPowerDouble(); + if (len > normalizationThreshold){ + res.setReal(c1.getRealDouble()/len); + res.setImaginary(c1.getImaginaryDouble()/len); + } else { + res.setComplexNumber(0, 0); + } + + } + + /* + * normalization with default threshold + * @param c1 + * @param res + */ + public static , S extends ComplexType> void normalize( R c1, S res){ + normalize(c1, res, 1e-5); + } + + + /* + * normalize complex valued img to length 1, pixel-wise, saving result to res + * if the length of a pixel is less than normalizationThreshold, set res to 0 + * @param img + * @param res + * @param normalizationThreshold + */ + public static , S extends ComplexType>void normalizeInterval( + final RandomAccessibleInterval img, final RandomAccessibleInterval res, final double normalizationThreshold, ExecutorService service) + { + + final Vector portions = FusionTools.divideIntoPortions( Views.iterable(img).size() ); + List> futures = new ArrayList>(); + final AtomicInteger ai = new AtomicInteger(-1); + + for (int i = 0; i < portions.size(); i++){ + futures.add(service.submit(new Runnable() { + + @Override + public void run() { + + ImagePortion ip = portions.get(ai.incrementAndGet()); + long loopSize = ip.getLoopSize(); + + if (Views.iterable(img).iterationOrder().equals(Views.iterable(res).iterationOrder())){ + Cursor cRes = Views.iterable(res).cursor(); + Cursor cSrc = Views.iterable(img).cursor(); + + cSrc.jumpFwd(ip.getStartPosition()); + cRes.jumpFwd(ip.getStartPosition()); + + for (long l = 0; l < loopSize; l++){ + cRes.fwd(); + cSrc.fwd(); + normalize(cSrc.get(), cRes.get(), normalizationThreshold); + } + } + + else { + Cursor cRes = Views.iterable(res).localizingCursor(); + RandomAccess raImg = img.randomAccess(); + + cRes.jumpFwd(ip.getStartPosition()); + + for (long l = 0; l < loopSize; l++){ + cRes.fwd(); + raImg.setPosition(cRes); + normalize(raImg.get(), cRes.get(), normalizationThreshold); + } + } + + } + })); + + for (Future f : futures){ + try { + f.get(); + } catch (InterruptedException e) { + // TODO Auto-generated catch block + e.printStackTrace(); + } catch (ExecutionException e) { + // TODO Auto-generated catch block + e.printStackTrace(); + } + } + } + + + + + } + /* + * normalization with default threshold + * @param img + * @param res + */ + public static , S extends ComplexType>void normalizeInterval( + RandomAccessibleInterval img, RandomAccessibleInterval res, ExecutorService service){ + normalizeInterval(img, res, 1E-5, service); + } + + /* + * get the mean pixel intensity of an img + * @param img + * @return + */ + public static > double getMean(RandomAccessibleInterval img) + { + // TODO: if #pixels > ???? else RealSum + // TODO: integral image? + double sum = 0.0; + long n = 0; + for (T pix: Views.iterable(img)){ + sum += pix.getRealDouble(); + n++; + } + return sum/n; + } + + /* + * get pixel-value correlation of two RandomAccessibleIntervals + * @param img1 + * @param img2 + * @return + */ + public static , S extends RealType> double getCorrelation ( + final RandomAccessibleInterval img1, final RandomAccessibleInterval img2) + { + final double m1 = getMean(img1); + final double m2 = getMean(img2); + + // square sums + double sum11 = 0.0, sum22 = 0.0, sum12 = 0.0; + + final Cursor c1 = Views.iterable(img1).cursor(); + + if (Views.iterable( img1 ).iterationOrder().equals( Views.iterable( img2 ).iterationOrder() )) + { + final Cursor< S > c2 = Views.iterable( img2 ).cursor(); + while (c1.hasNext()){ + final double c = c1.next().getRealDouble(); + final double r = c2.next().getRealDouble(); + + sum11 += (c - m1) * (c - m1); + sum22 += (r - m2) * (r - m2); + sum12 += (c - m1) * (r - m2); + } + } + else + { + final RandomAccess r2 = img2.randomAccess(); + while (c1.hasNext()){ + final double c = c1.next().getRealDouble(); + r2.setPosition(c1); + final double r = r2.get().getRealDouble(); + + sum11 += (c - m1) * (c - m1); + sum22 += (r - m2) * (r - m2); + sum12 += (c - m1) * (r - m2); + } + } + + // all pixels had the same color.... + if (sum11 == 0 || sum22 == 0) + { + // having the same means and same sums means the overlapping area was simply identically the same color + // this is most likely an artifact and we return 0 + /* if ( sum11 == sum22 && m1 == m2 ) + return 1; + else */ + return 0; + } + + return sum12 / Math.sqrt(sum11 * sum22); + } + + + /* + * test stitching, create new image with img2 copied over img1 at the specified shift + * @param img1 + * @param img2 + * @param shiftPeak + * @return + */ + public static , S extends RealType> RandomAccessibleInterval dummyFuse(RandomAccessibleInterval img1, RandomAccessibleInterval img2, PhaseCorrelationPeak2 shiftPeak, ExecutorService service) + { + long[] shift = new long[img1.numDimensions()]; + shiftPeak.getShift().localize(shift); + long[] minImg1 = new long[img1.numDimensions()]; + long[] minImg2 = new long[img1.numDimensions()]; + long[] maxImg1 = new long[img1.numDimensions()]; + long[] maxImg2 = new long[img1.numDimensions()]; + long[] min = new long[img1.numDimensions()]; + long[] max = new long[img1.numDimensions()]; + + for (int i = 0; i < img1.numDimensions(); i++){ + minImg1[i] = 0; + maxImg1[i] = img1.dimension(i) -1; + minImg2[i] = shiftPeak.getShift().getLongPosition(i); + maxImg2[i] = img2.dimension(i) + minImg2[i] - 1; + + min[i] = Math.min(minImg1[i], minImg2[i]); + max[i] = Math.max(maxImg1[i], maxImg2[i]); + } + + + RandomAccessibleInterval res = new ArrayImgFactory().create(new FinalInterval(min, max), new FloatType()); + copyRealImage(Views.iterable(img1), Views.translate(res, min), service); + copyRealImage(Views.iterable(Views.translate(img2, shift)), Views.translate(res, min), service); + return res; + + } + + public static void main( String[] args ) + { + final Point p = new Point( 90, 90 ); // identical to (-10,-10), so subpixel localization can move on periodic condition outofbounds + PhaseCorrelationPeak2 pcp = new PhaseCorrelationPeak2( p, 5 ); + + Dimensions pcmDims = new FinalDimensions( 100, 100 ); + Dimensions p1 = new FinalDimensions( 80, 81 ); + Dimensions p2 = new FinalDimensions( 91, 90 ); + + final List peaks = expandPeakToPossibleShifts( pcp, pcmDims, p1, p2 ); + + for ( final PhaseCorrelationPeak2 pc : peaks ) + System.out.println( Util.printCoordinates( pc.getShift() ) ); + + Img< FloatType > a = ImgLib2Util.openAs32Bit( new File( "73.tif.zip" ) ); + Img< FloatType > b = ImgLib2Util.openAs32Bit( new File( "74.tif.zip" ) ); + +// BenchmarkHelper.benchmarkAndPrint( 10, true, new Runnable() +// { +// @Override +// public void run() +// { +// System.out.println( getCorrelation ( a, b ) ); +// } +// } ); + + // Commented out to avoid circular dependency with BigStitcher +// BenchmarkHelper.benchmarkAndPrint( 10, true, new Runnable() +// { +// @Override +// public void run() +// { +// PairwiseStitching.getShift( a, b, new Translation3D(), new Translation3D(), +// new PairwiseStitchingParameters(), Executors.newFixedThreadPool( Runtime.getRuntime().availableProcessors() ) ); +// } +// } ); + +// System.out.println( getCorrelation ( a, b ) ); +// System.out.println( getCorrelation ( a, c ) ); +// System.out.println( getCorrelation ( b, c ) ); + } +} diff --git a/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/PhaseCorrelationPeak2.java b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/PhaseCorrelationPeak2.java new file mode 100644 index 000000000..69385edb3 --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/PhaseCorrelationPeak2.java @@ -0,0 +1,406 @@ +/*- + * #%L + * Multiview stitching of large datasets. + * %% + * Copyright (C) 2016 - 2025 Big Stitcher developers. + * %% + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program. If not, see + * . + * #L% + */ +package net.preibisch.mvrecon.process.phasecorrelation; + +import java.util.ArrayList; +import java.util.Comparator; +import java.util.List; + +import net.imglib2.Dimensions; +import net.imglib2.FinalDimensions; +import net.imglib2.Interval; +import net.imglib2.Localizable; +import net.imglib2.Point; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.RealLocalizable; +import net.imglib2.RealPoint; +import net.imglib2.RealRandomAccessible; +import net.imglib2.algorithm.localextrema.RefinedPeak; +import net.imglib2.algorithm.localextrema.SubpixelLocalization; +import net.imglib2.interpolation.randomaccess.NLinearInterpolatorFactory; +import net.imglib2.realtransform.InvertibleRealTransform; +import net.imglib2.realtransform.RealViews; +import net.imglib2.realtransform.Translation2D; +import net.imglib2.realtransform.Translation3D; +import net.imglib2.type.numeric.RealType; +import net.imglib2.util.Pair; +import net.imglib2.view.Views; + +public class PhaseCorrelationPeak2 { + + Localizable pcmLocation; // location in the raw PCM + RealLocalizable subpixelPcmLocation; // subpixel localized PCM location + + Localizable shift; // corresponding shift between images + RealLocalizable subpixelShift; // corresponding subpixel accurate shift between images + + double phaseCorr; // value in PCM + double crossCorr; // crosscorrelation bewtween imgs + long nPixel; // number of overlapping pixels + + public Localizable getPcmLocation() { + return pcmLocation; + } + + public void setPcmLocation(Localizable pcmLocation) { + this.pcmLocation = pcmLocation; + } + + public Localizable getShift() { + return shift; + } + + public void setShift(Localizable shift) { + this.shift = shift; + } + + public RealLocalizable getSubpixelPcmLocation() { + return subpixelPcmLocation; + } + + public void setSubpixelPcmLocation(RealLocalizable subpixelPcmLocation) { + this.subpixelPcmLocation = subpixelPcmLocation; + } + + public RealLocalizable getSubpixelShift() { + return subpixelShift; + } + + public void setSubpixelShift(RealLocalizable subpixelShift) { + this.subpixelShift = subpixelShift; + } + + public double getPhaseCorr() { + return phaseCorr; + } + + public void setPhaseCorr(double phaseCorr) { + this.phaseCorr = phaseCorr; + } + + public double getCrossCorr() { + return crossCorr; + } + + public void setCrossCorr(double crossCorr) { + this.crossCorr = crossCorr; + } + + public long getnPixel() { + return nPixel; + } + + public void setnPixel(long nPixel) { + this.nPixel = nPixel; + } + + /* + * constructor with just raw PCM location and value + * @param pcmPosition + * @param phaseCorr + */ + public PhaseCorrelationPeak2(Localizable pcmPosition, double phaseCorr) + { + this.pcmLocation = new Point( pcmPosition ); + this.shift = null; + this.subpixelPcmLocation = null; + this.subpixelShift = null; + this.phaseCorr = phaseCorr; + this.crossCorr = 0.0; + this.nPixel = 0; + } + + /* + * copy contructor + * @param src + */ + public PhaseCorrelationPeak2(PhaseCorrelationPeak2 src){ + this.pcmLocation = new Point(src.pcmLocation); + this.shift = src.shift == null? null : new Point(src.shift); + this.subpixelPcmLocation = src.subpixelPcmLocation == null ? null : new RealPoint( src.subpixelPcmLocation ); + this.subpixelShift = src.subpixelShift == null ? null : new RealPoint( src.subpixelShift ); + this.phaseCorr = src.phaseCorr; + this.crossCorr = src.crossCorr; + this.nPixel = src.nPixel; + } + + public static class ComparatorByPhaseCorrelation implements Comparator { + @Override + public int compare(PhaseCorrelationPeak2 o1, PhaseCorrelationPeak2 o2) { + return Double.compare(o1.phaseCorr, o2.phaseCorr); + } + } + + public static class ComparatorByCrossCorrelation implements Comparator{ + @Override + public int compare(PhaseCorrelationPeak2 o1, PhaseCorrelationPeak2 o2) { + int ccCompare = Double.compare(o1.crossCorr, o2.crossCorr); + if (ccCompare != 0){ + return ccCompare; + } else { + return (int)(o1.nPixel - o2.nPixel); + } + } + } + + public , S extends RealType> void calculateCrossCorr(RandomAccessibleInterval img1, RandomAccessibleInterval img2, + long minOverlapPx) + { + calculateCrossCorr( img1, img2, minOverlapPx, false ); + } + + /* + * checks the cross correlation of two images shifted as indicated by this phaseCorrelationPeak, + * update the values of crossCor and nPixels accordingly + * if the images do not overlap (or overlap less than a specified minimal overlap), + * crossCorr and nPixel will be set to 0 + * @param img1 + * @param img2 + * @param minOverlapPx + */ + + public , S extends RealType> void calculateCrossCorr(RandomAccessibleInterval img1, RandomAccessibleInterval img2, + long minOverlapPx, boolean interpolateSubpixel) + { + Pair intervals = PhaseCorrelation2Util.getOverlapIntervals(img1, img2, shift); + + // no overlap found + if (intervals == null) { + crossCorr = Double.NEGATIVE_INFINITY; + nPixel = 0; + return; + } + + nPixel = 1; + for (int i = 0; i< intervals.getA().numDimensions(); i++){ + nPixel *= intervals.getA().dimension(i); + } + + if (nPixel < minOverlapPx){ + crossCorr = Double.NEGATIVE_INFINITY; + nPixel = 0; + return; + } + + // for subpixel move the underlying Img2 by the subpixel offset + if ( subpixelShift != null && interpolateSubpixel ) + { + RealRandomAccessible< S > rra = Views.interpolate( Views.extendMirrorSingle( img2 ), new NLinearInterpolatorFactory< S >() ); + + InvertibleRealTransform transform = null; + + // e.g. subpixel = (-0.4, 0.1, -0.145) + final double tx = subpixelShift.getDoublePosition( 0 ) - shift.getDoublePosition( 0 ); + final double ty = subpixelShift.getDoublePosition( 1 ) - shift.getDoublePosition( 1 ); + + if ( rra.numDimensions() == 2 ) + transform = new Translation2D( -tx, -ty ); // -relative subpixel shift only + else if ( rra.numDimensions() == 3 ) + transform = new Translation3D( -tx, -ty, shift.getDoublePosition( 2 ) - subpixelShift.getDoublePosition( 2 ) ); // -relative subpixel shift only + + img2 = Views.interval( Views.raster( RealViews.transform( rra, transform ) ), img2 ); + } + + // calculate cross correlation. + // note that the overlap we calculate assumes zero-min input + crossCorr = PhaseCorrelation2Util.getCorrelation( + Views.zeroMin( Views.interval(Views.zeroMin(img1), intervals.getA())), + Views.zeroMin( Views.interval(Views.zeroMin(img2), intervals.getB())) + ); + + } + + /* + * calculate cross correlation of two images with no minimal overlap size + * @param img1 + * @param img2 + */ + public , S extends RealType> void calculateCrossCorr(RandomAccessibleInterval img1, RandomAccessibleInterval img2){ + calculateCrossCorr(img1, img2, 0); + } + + /* + * refine the shift using subpixel localization in the original PCM + * @param pcm + */ + public > void calculateSubpixelLocalization(RandomAccessibleInterval pcm){ + + List peaks = new ArrayList(); + peaks.add(new Point(pcmLocation)); + + final int n = pcm.numDimensions(); + + boolean[] allowedToMoveInDim = new boolean[ n ]; + for (int i = 0; i< allowedToMoveInDim.length; i++){ + allowedToMoveInDim[i]=false; + } + + // TODO: It doesnt look like this does anything? Subpixel peaks are just regular peaks as RealPoint - with maxNumMoves == 1 it should now :) + // subpixel localization can move on periodic condition outofbounds + List> res = SubpixelLocalization.refinePeaks(peaks, Views.extendPeriodic( pcm ), null, true, + 1, false, 0.0f, allowedToMoveInDim); + + final RefinedPeak< Point > rp = res.get( 0 ); + this.subpixelPcmLocation = new RealPoint( rp ); + /* + final long[] range = Util.getArrayFromValue( 4l, pcm.numDimensions() ); + final Gradient derivative = new GradientOnDemand< T >( Views.extendPeriodic( pcm ) ); + final ArrayList< long[] > peaksRS = new ArrayList<>(); + + final long[] loc = new long[ pcm.numDimensions() ]; + pcmLocation.localize( loc ); + peaksRS.add( loc ); + + final ArrayList< Spot > spots = Spot.extractSpots( peaksRS, derivative, range ); + try + { + Spot.fitCandidates( spots ); + System.out.println( "rs: " + Util.printCoordinates( spots.get( 0 ) ) ); + //Spot.ransac( spots, 1000, 0.1, 1.0/100.0 ); + //System.out.println( "rsransac: " + Util.printCoordinates( spots.get( 0 ) ) ); + + } catch ( Exception e ) + { + // TODO Auto-generated catch block + e.printStackTrace(); + } + + final long[] minSF = new long[ pcm.numDimensions() ]; + final long[] maxSF = new long[ pcm.numDimensions() ]; + + pcmLocation.localize( minSF ); + pcmLocation.localize( maxSF ); + + for ( int d = 0; d < pcm.numDimensions(); ++d ) + { + minSF[ d ] -= 2; + maxSF[ d ] += 2; + } + + final Cursor< T > cursor = Views.iterable( Views.interval( Views.extendPeriodic( pcm ), minSF, maxSF ) ).localizingCursor(); + + double[] sumW = new double[ pcm.numDimensions() ]; + double[] sum = new double[ pcm.numDimensions() ]; + + while ( cursor.hasNext() ) + { + double w = cursor.next().getRealDouble(); + + for ( int d = 0; d < pcm.numDimensions(); ++d ) + { + sumW[ d ] += w; + sum[ d ] += cursor.getLongPosition( d ) * w; + } + + //System.out.println( Util.printCoordinates( cursor ) + ": " + w ); + } + + double[] simpleFit = new double[ pcm.numDimensions() ]; + for ( int d = 0; d < pcm.numDimensions(); ++d ) + simpleFit[ d ] = sum[ d ] / sumW[ d ]; + + long size = 41; + Img< FloatType > img = ArrayImgs.floats( Util.getArrayFromValue( size, pcm.numDimensions() ) ); + Cursor< FloatType > c = img.localizingCursor(); + RandomAccess< T > r = Views.extendPeriodic( pcm ).randomAccess(); + long[] pos = new long[ img.numDimensions() ]; + + while ( c.hasNext() ) + { + c.fwd(); + c.localize( pos ); + for ( int d = 0; d < pcm.numDimensions(); ++d ) + { + // relative to center + pos[ d ] = size / 2 - pos[ d ]; + pos[ d ] += pcmLocation.getLongPosition( d ); + } + r.setPosition( pos ); + c.get().set( r.get().getRealFloat() ); + } + + Collection peaksG = new ArrayList(); + peaksG.add( new Point( Util.getArrayFromValue( size/2, pcm.numDimensions() ) ) ); + + PeakFitter pf = new PeakFitter<>(img, peaksG, + new LevenbergMarquardtSolver(), new EllipticGaussianOrtho(), + new MLEllipticGaussianEstimator( Util.getArrayFromValue( 0.9, pcm.numDimensions() ) ) ); + + //ImageJFunctions.show( img ); + if ( !pf.process() ) + System.out.println( pf.getErrorMessage() ); + + final double[] peakGauss = pf.getResult().values().iterator().next(); + + for ( int d = 0; d < pcm.numDimensions(); ++d ) + peakGauss[ d ] = peakGauss[ d ] - size/2 + pcmLocation.getLongPosition( d ); + + System.out.println( "gf: " + Util.printCoordinates( peakGauss ) ); + System.out.println( "sf: " + Util.printCoordinates( simpleFit ) ); + System.out.println( "qf: " + Util.printCoordinates( subpixelPcmLocation ) ); + System.out.println(); + + //this.subpixelPcmLocation = new RealPoint( spots.get( 0 ) ); + */ + + double maxDist = 0.0; + + for ( int d = 0; d < n; ++d ) + maxDist = Math.max( maxDist, Math.abs( rp.getOriginalPeak().getDoublePosition( d ) - this.subpixelPcmLocation.getDoublePosition( d ) ) ); + + //SimpleMultiThreading.threadHaltUnClean(); + // not a stable peak + if ( maxDist > 0.75 ) + { + this.subpixelPcmLocation = null; + } + else + { + this.phaseCorr = rp.getValue(); + } + } + + public static void main(String[] args) { + + double o1 = 6; + double o2 = Double.NEGATIVE_INFINITY; + int np1 = 30; + int np2 = 20; + + System.out.println( Double.isInfinite( o2 )); + + int ccCompare = Double.compare(o1, o2); + if (ccCompare != 0) + System.out.println( ccCompare ); + else + System.out.println( (int)(np1 - np2) ); + + System.exit( 0 ); + PhaseCorrelationPeak2 peaks = new PhaseCorrelationPeak2(new Point(new int[] {10, 10}), 1.0); + Dimensions pcmDims = new FinalDimensions(new int[] {50, 50}); + Dimensions imgDims = new FinalDimensions(new int[] {30, 30}); + PhaseCorrelation2Util.expandPeakToPossibleShifts(peaks, pcmDims, imgDims, imgDims); + + } + + +} diff --git a/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/deprecated/Blending.java b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/deprecated/Blending.java new file mode 100644 index 000000000..d4843a8a3 --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/deprecated/Blending.java @@ -0,0 +1,78 @@ +/*- + * #%L + * Multiview stitching of large datasets. + * %% + * Copyright (C) 2016 - 2025 Big Stitcher developers. + * %% + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program. If not, see + * . + * #L% + */ +package net.preibisch.mvrecon.process.phasecorrelation.deprecated; + +import net.imglib2.FinalInterval; +import net.imglib2.Interval; +import net.imglib2.RealInterval; +import net.imglib2.RealRandomAccess; +import net.imglib2.RealRandomAccessible; +import net.imglib2.type.numeric.real.FloatType; + +/** + * + * RealRandomAccess that computed cosine-blending for a certain interval + * + * @author Stephan Preibisch (stephan.preibisch@gmx.de) + * @deprecated + */ +public class Blending implements RealRandomAccessible< FloatType > +{ + final Interval interval; + final float[] border, blending; + + /** + * RealRandomAccess that computes a blending function for a certain {@link Interval} + * + * @param interval - the interval it is defined on (return zero outside of it) + * @param border - how many pixels to skip before starting blending (on each side of each dimension) + * @param blending - how many pixels to compute the blending function on (on each side of each dimension) + */ + public Blending( final Interval interval, final float[] border, final float[] blending ) + { + // in case the interval is actually image data re-instantiate just a simple FinalInterval + this.interval = new FinalInterval( interval ); + this.border = border; + this.blending = blending; + } + + @Override + public int numDimensions() { return interval.numDimensions(); } + + @Override + public RealRandomAccess realRandomAccess() + { + return new BlendingRealRandomAccess( interval, border, blending ); + } + + @Override + public RealRandomAccess realRandomAccess( final RealInterval interval ) + { + return realRandomAccess(); + } + + @Override + public FloatType getType() + { + return new FloatType(); + } +} diff --git a/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/deprecated/BlendingRealRandomAccess.java b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/deprecated/BlendingRealRandomAccess.java new file mode 100644 index 000000000..f1aec6cff --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/process/phasecorrelation/deprecated/BlendingRealRandomAccess.java @@ -0,0 +1,296 @@ +/*- + * #%L + * Multiview stitching of large datasets. + * %% + * Copyright (C) 2016 - 2025 Big Stitcher developers. + * %% + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program. If not, see + * . + * #L% + */ +package net.preibisch.mvrecon.process.phasecorrelation.deprecated; + +import ij.ImageJ; +import net.imglib2.Cursor; +import net.imglib2.Interval; +import net.imglib2.Localizable; +import net.imglib2.RealLocalizable; +import net.imglib2.RealRandomAccess; +import net.imglib2.Sampler; +import net.imglib2.img.Img; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.img.display.imagej.ImageJFunctions; +import net.imglib2.type.numeric.real.FloatType; + +/** + * + * @author preibisch + * + * @deprecated + */ +public class BlendingRealRandomAccess implements RealRandomAccess< FloatType > +{ + final Interval interval; + final int[] min, dimMinus1; + final float[] l, border, blending; + final int n; + final FloatType v; + + // static lookup table for the blending function + final static private double[] lookUp; + private static final int indexFor( final double d ) { return (int)Math.round( d * 1000.0 ); } + + static + { + lookUp = new double[ 1001 ]; + + for ( double d = 0; d <= 1.0001; d = d + 0.001 ) + lookUp[ indexFor( d ) ] = ( Math.cos( ( 1 - d ) * Math.PI ) + 1 ) / 2; + } + + /** + * RealRandomAccess that computes a blending function for a certain {@link Interval} + * + * @param interval - the interval it is defined on (return zero outside of it) + * @param border - how many pixels to skip before starting blending (on each side of each dimension) + * @param blending - how many pixels to compute the blending function on (on each side of each dimension) + */ + public BlendingRealRandomAccess( + final Interval interval, + final float[] border, + final float[] blending ) + { + this.interval = interval; + this.n = interval.numDimensions(); + this.l = new float[ n ]; + this.border = border; + this.blending = blending; + this.v = new FloatType(); + + this.min = new int[ n ]; + this.dimMinus1 = new int[ n ]; + + for ( int d = 0; d < n; ++d ) + { + this.min[ d ] = (int)interval.min( d ); + this.dimMinus1[ d ] = (int)interval.max( d ) - min[ d ]; + } + } + + @Override + public FloatType get() + { + v.set( computeWeight( l, min, dimMinus1, border, blending, n ) ); + return v; + } + + final private static float computeWeight( + final float[] location, + final int[] min, + final int[] dimMinus1, + final float[] border, + final float[] blending, + final int n ) + { + // compute multiplicative distance to the respective borders [0...1] + float minDistance = 1; + + for ( int d = 0; d < n; ++d ) + { + // the position in the image relative to the boundaries and the border + final float l = ( location[ d ] - min[ d ] ); + + // the distance to the border that is closer + final float dist = Math.max( 0, Math.min( l - border[ d ], dimMinus1[ d ] - l - border[ d ] ) ); + + // if this is 0, the total result will be 0, independent of the number of dimensions + if ( dist == 0 ) + return 0; + + final float relDist = dist / blending[ d ]; + + if ( relDist < 1 ) + minDistance *= lookUp[ indexFor( relDist ) ]; //( Math.cos( ( 1 - relDist ) * Math.PI ) + 1 ) / 2; + } + + return minDistance; + } + + @Override + public void localize( final float[] position ) + { + for ( int d = 0; d < n; ++d ) + position[ d ] = l[ d ]; + } + + @Override + public void localize( final double[] position ) + { + for ( int d = 0; d < n; ++d ) + position[ d ] = l[ d ]; + } + + @Override + public float getFloatPosition( final int d ){ return l[ d ]; } + + @Override + public double getDoublePosition( final int d ) { return l[ d ]; } + + @Override + public int numDimensions() { return n; } + + @Override + public void move( final float distance, final int d ) { l[ d ] += distance; } + + @Override + public void move( final double distance, final int d ) { l[ d ] += distance; } + + @Override + public void move( final RealLocalizable localizable ) + { + for ( int d = 0; d < n; ++d ) + l[ d ] += localizable.getFloatPosition( d ); + } + + @Override + public void move( final float[] distance ) + { + for ( int d = 0; d < n; ++d ) + l[ d ] += distance[ d ]; + } + + @Override + public void move( final double[] distance ) + { + for ( int d = 0; d < n; ++d ) + l[ d ] += distance[ d ]; + } + + @Override + public void setPosition( final RealLocalizable localizable ) + { + for ( int d = 0; d < n; ++d ) + l[ d ] = localizable.getFloatPosition( d ); + } + + @Override + public void setPosition( final float[] position ) + { + for ( int d = 0; d < n; ++d ) + l[ d ] = position[ d ]; + } + + @Override + public void setPosition( final double[] position ) + { + for ( int d = 0; d < n; ++d ) + l[ d ] =(float)position[ d ]; + } + + @Override + public void setPosition( final float position, final int d ) { l[ d ] = position; } + + @Override + public void setPosition( final double position, final int d ) { l[ d ] = (float)position; } + + @Override + public void fwd( final int d ) { ++l[ d ]; } + + @Override + public void bck( final int d ) { --l[ d ]; } + + @Override + public void move( final int distance, final int d ) { l[ d ] += distance; } + + @Override + public void move( final long distance, final int d ) { l[ d ] += distance; } + + @Override + public void move( final Localizable localizable ) + { + for ( int d = 0; d < n; ++d ) + l[ d ] += localizable.getFloatPosition( d ); + } + + @Override + public void move( final int[] distance ) + { + for ( int d = 0; d < n; ++d ) + l[ d ] += distance[ d ]; + } + + @Override + public void move( final long[] distance ) + { + for ( int d = 0; d < n; ++d ) + l[ d ] += distance[ d ]; + } + + @Override + public void setPosition( final Localizable localizable ) + { + for ( int d = 0; d < n; ++d ) + l[ d ] = localizable.getFloatPosition( d ); + } + + @Override + public void setPosition( final int[] position ) + { + for ( int d = 0; d < n; ++d ) + l[ d ] = position[ d ]; + } + + @Override + public void setPosition( final long[] position ) + { + for ( int d = 0; d < n; ++d ) + l[ d ] = position[ d ]; + } + + @Override + public void setPosition( final int position, final int d ) { l[ d ] = position; } + + @Override + public void setPosition( final long position, final int d ) { l[ d ] = position; } + + @Override + public RealRandomAccess copy() + { + final BlendingRealRandomAccess r = new BlendingRealRandomAccess( interval, border, blending ); + r.setPosition( this ); + return r; + } + + public static void main( String[] args ) + { + new ImageJ(); + + Img< FloatType > img = ArrayImgs.floats( 500, 500 ); + BlendingRealRandomAccess blend = new BlendingRealRandomAccess( + img, + new float[]{ 100, 0 }, + new float[]{ 12, 150 } ); + + Cursor< FloatType > c = img.localizingCursor(); + + while ( c.hasNext() ) + { + c.fwd(); + blend.setPosition( c ); + c.get().setReal( blend.get().getRealFloat() ); + } + + ImageJFunctions.show( img ); + } +} diff --git a/src/test/java/net/preibisch/mvrecon/process/phasecorrelation/FourNeighborhoodExtremaTest.java b/src/test/java/net/preibisch/mvrecon/process/phasecorrelation/FourNeighborhoodExtremaTest.java new file mode 100644 index 000000000..d4e3e0351 --- /dev/null +++ b/src/test/java/net/preibisch/mvrecon/process/phasecorrelation/FourNeighborhoodExtremaTest.java @@ -0,0 +1,344 @@ +/*- + * #%L + * Multiview stitching of large datasets. + * %% + * Copyright (C) 2016 - 2025 Big Stitcher developers. + * %% + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program. If not, see + * . + * #L% + */ +package net.preibisch.mvrecon.process.phasecorrelation; + +import static org.junit.Assert.*; + +import java.util.ArrayList; +import java.util.Collections; +import java.util.Comparator; +import java.util.Random; +import java.util.concurrent.Executors; + +import net.imglib2.Localizable; +import net.imglib2.Point; +import net.imglib2.RandomAccess; +import net.imglib2.img.Img; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.img.display.imagej.ImageJFunctions; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.util.Pair; +import net.imglib2.util.ValuePair; +import net.imglib2.view.Views; + +import org.junit.Test; + +public class FourNeighborhoodExtremaTest +{ + + public static long seed = 4353; + + @Test + public void testRandomPeaks() + { + Img< FloatType > img = ArrayImgs.floats( 50, 50 ); + Random rnd = new Random( seed ); + + for( FloatType t : img ) + t.set( rnd.nextFloat() ); + + ArrayList< Pair< Localizable, Double > > correct = new ArrayList< Pair >(); + + RandomAccess ra = img.randomAccess(); + + for ( int i = 0; i < 10; ++i ){ + for (int d = 0; d< img.numDimensions(); d++){ + ra.setPosition((int) (rnd.nextDouble() * 50), d); + } + ra.get().set((float) (rnd.nextDouble() + 1.0)); + + correct.add(new ValuePair(new Point(ra), new Double(ra.get().get()))); + } + + // sort the peaks in descending order + Collections.sort(correct, new Comparator>() { + @Override + public int compare(Pair o1, Pair o2) { + return Double.compare(o2.getB(), o1.getB()); + } + }); + + int nMax = 5; + + ArrayList< Pair< Localizable, Double > > found = FourNeighborhoodExtrema.findMax(Views.extendPeriodic(img), img, nMax); + + + assertEquals(nMax, found.size()); + + + long[] posCorrect = new long[img.numDimensions()]; + long[] posFound = new long[img.numDimensions()]; + + for (int i = 0; i img = ArrayImgs.floats( 50, 50 ); + Random rnd = new Random( seed ); + + for( FloatType t : img ) + t.set( rnd.nextFloat() ); + + ArrayList< Pair< Localizable, Double > > correct = new ArrayList< Pair >(); + + RandomAccess ra = img.randomAccess(); + + for ( int i = 0; i < 10; ++i ){ + for (int d = 0; d< img.numDimensions(); d++){ + ra.setPosition((int) (rnd.nextDouble() * 50), d); + } + ra.get().set((float) (rnd.nextDouble() + 1.0)); + + correct.add(new ValuePair(new Point(ra), new Double(ra.get().get()))); + } + + // sort the peaks in descending order + Collections.sort(correct, new Comparator>() { + @Override + public int compare(Pair o1, Pair o2) { + return Double.compare(o2.getB(), o1.getB()); + } + }); + + int nMax = 5; + ArrayList< Pair< Localizable, Double > > found = FourNeighborhoodExtrema.findMaxMT(Views.extendPeriodic(img), img, nMax, Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors())); + assertEquals(nMax, found.size()); + + + long[] posCorrect = new long[img.numDimensions()]; + long[] posFound = new long[img.numDimensions()]; + + for (int i = 0; i img = ArrayImgs.floats( 50, 50 ); + Random rnd = new Random( seed ); + + for( FloatType t : img ) + t.set( rnd.nextFloat() ); + + ArrayList< Pair< Localizable, Double > > correct = new ArrayList< Pair >(); + + RandomAccess ra = img.randomAccess(); + + ra.setPosition(new long[]{0,0}); + ra.get().set((float) (rnd.nextDouble() + 1.0)); + correct.add(new ValuePair(new Point(ra), new Double(ra.get().get()))); + + ra.setPosition(new long[]{0,49}); + ra.get().set((float) (rnd.nextDouble() + 1.0)); + correct.add(new ValuePair(new Point(ra), new Double(ra.get().get()))); + + ra.setPosition(new long[]{49,0}); + ra.get().set((float) (rnd.nextDouble() + 1.0)); + correct.add(new ValuePair(new Point(ra), new Double(ra.get().get()))); + + ra.setPosition(new long[]{49,49}); + ra.get().set((float) (rnd.nextDouble() + 1.0)); + correct.add(new ValuePair(new Point(ra), new Double(ra.get().get()))); + + // sort the peaks in descending order + Collections.sort(correct, new Comparator>() { + @Override + public int compare(Pair o1, Pair o2) { + return Double.compare(o2.getB(), o1.getB()); + } + }); + + int nMax = 4; + + ArrayList< Pair< Localizable, Double > > found = FourNeighborhoodExtrema.findMaxMT(Views.extendPeriodic(img), img, nMax, Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors())); + + + assertEquals(nMax, found.size()); + + + long[] posCorrect = new long[img.numDimensions()]; + long[] posFound = new long[img.numDimensions()]; + + for (int i = 0; i img = ArrayImgs.floats( 50, 50 ); + Random rnd = new Random( seed ); + + for( FloatType t : img ) + t.set( rnd.nextFloat() ); + + ArrayList< Pair< Localizable, Double > > correct = new ArrayList< Pair >(); + + RandomAccess ra = img.randomAccess(); + + ra.setPosition(new long[]{0,0}); + ra.get().set((float) (rnd.nextDouble() + 1.0)); + correct.add(new ValuePair(new Point(ra), new Double(ra.get().get()))); + + ra.setPosition(new long[]{0,49}); + ra.get().set((float) (rnd.nextDouble() + 1.0)); + correct.add(new ValuePair(new Point(ra), new Double(ra.get().get()))); + + ra.setPosition(new long[]{49,0}); + ra.get().set((float) (rnd.nextDouble() + 1.0)); + correct.add(new ValuePair(new Point(ra), new Double(ra.get().get()))); + + ra.setPosition(new long[]{49,49}); + ra.get().set((float) (rnd.nextDouble() + 1.0)); + correct.add(new ValuePair(new Point(ra), new Double(ra.get().get()))); + + // sort the peaks in descending order + Collections.sort(correct, new Comparator>() { + @Override + public int compare(Pair o1, Pair o2) { + return Double.compare(o2.getB(), o1.getB()); + } + }); + + int nMax = 4; + + ArrayList< Pair< Localizable, Double > > found = FourNeighborhoodExtrema.findMaxMT(Views.extendMirrorSingle(img), img, nMax, Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors())); + + + assertEquals(nMax, found.size()); + + + long[] posCorrect = new long[img.numDimensions()]; + long[] posFound = new long[img.numDimensions()]; + + for (int i = 0; i img = ArrayImgs.floats( 50, 50 ); + Random rnd = new Random( seed ); + + for( FloatType t : img ) + t.set( rnd.nextFloat() ); + + ArrayList< Pair< Localizable, Double > > correct = new ArrayList< Pair >(); + + RandomAccess ra = img.randomAccess(); + + ra.setPosition(new long[]{2,0}); + ra.get().set((float) 2.0); + correct.add(new ValuePair(new Point(ra), new Double(ra.get().get()))); + + ra.setPosition(new long[]{1,0}); + ra.get().set((float) 2.0); + correct.add(new ValuePair(new Point(ra), new Double(ra.get().get()))); + + ra.setPosition(new long[]{3,0}); + ra.get().set((float) 2.0); + correct.add(new ValuePair(new Point(ra), new Double(ra.get().get()))); + + ra.setPosition(new long[]{4,0}); + ra.get().set((float) 2.0); + correct.add(new ValuePair(new Point(ra), new Double(ra.get().get()))); + + // sort the peaks in descending order + Collections.sort(correct, new Comparator>() { + @Override + public int compare(Pair o1, Pair o2) { + return Double.compare(o2.getB(), o1.getB()); + } + }); + + int nMax = 4; + + ArrayList< Pair< Localizable, Double > > found = FourNeighborhoodExtrema.findMaxMT(Views.extendPeriodic(img), img, nMax, Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors())); + + + assertEquals(nMax, found.size()); + + + long[] posCorrect = new long[img.numDimensions()]; + long[] posFound = new long[img.numDimensions()]; + + for (int i = 0; i. + * #L% + */ +package net.preibisch.mvrecon.process.phasecorrelation; + +import static org.junit.Assert.assertArrayEquals; +import static org.junit.Assert.assertTrue; + +import java.util.Arrays; +import java.util.Random; +import java.util.concurrent.Executors; + +import org.junit.Test; + +import net.imglib2.FinalInterval; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.img.Img; +import net.imglib2.img.array.ArrayImgFactory; +import net.imglib2.img.array.ArrayImgs; +import net.imglib2.interpolation.randomaccess.NLinearInterpolatorFactory; +import net.imglib2.realtransform.AffineGet; +import net.imglib2.realtransform.AffineRandomAccessible; +import net.imglib2.realtransform.RealViews; +import net.imglib2.realtransform.Translation2D; +import net.imglib2.type.numeric.complex.ComplexFloatType; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.util.Intervals; +import net.imglib2.util.Util; +import net.imglib2.view.IntervalView; +import net.imglib2.view.Views; + +public class PhaseCorrelationTest { + + public static long seed = 4353; + + @Test + public void testPC() { + + // TODO: very large shifts (nearly no overlap) lead to incorrect shift determination (as expected) + // maybe we can optimize behaviour in this situation + Img< FloatType > img = ArrayImgs.floats( 200, 200 ); + Random rnd = new Random( seed ); + + for( FloatType t : img ) + t.set( rnd.nextFloat()); + + long shiftX = 28; + long shiftY = 0; + + FinalInterval interval1 = new FinalInterval(new long[] {50, 50}); + FinalInterval interval2 = Intervals.translate(interval1, shiftX, 0); + interval2 = Intervals.translate(interval2, shiftY, 1); + + + int [] extension = new int[img.numDimensions()]; + Arrays.fill(extension, 10); + + RandomAccessibleInterval pcm = PhaseCorrelation2.calculatePCM(Views.interval(img, interval1), Views.zeroMin(Views.interval(img, interval2)), extension, new ArrayImgFactory(), + new FloatType(), new ArrayImgFactory(), new ComplexFloatType(), Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors())); + + PhaseCorrelationPeak2 shiftPeak = PhaseCorrelation2.getShift(pcm, Views.interval(img, interval1), Views.zeroMin(Views.interval(img, interval2)), 20, 0, false, Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors())); + + long[] expected = new long[]{shiftX, shiftY}; + long[] found = new long[img.numDimensions()]; + + + + shiftPeak.getShift().localize(found); + + assertArrayEquals(expected, found); + + } + + + @Test + public void testPCNegativeShift() { + + // TODO: very large shifts (nearly no overlap) lead to incorrect shift determination (as expected) + // maybe we can optimize behaviour in this situation + Img< FloatType > img = ArrayImgs.floats( 200, 200 ); + Random rnd = new Random( seed ); + + for( FloatType t : img ) + t.set( rnd.nextFloat()); + + long shiftX = -2; + long shiftY = -2; + + FinalInterval interval1 = new FinalInterval(new long[] {50, 50}); + FinalInterval interval2 = Intervals.translate(interval1, shiftX, 0); + interval2 = Intervals.translate(interval2, shiftY, 1); + + int [] extension = new int[img.numDimensions()]; + Arrays.fill(extension, 10); + + RandomAccessibleInterval pcm = PhaseCorrelation2.calculatePCM( + Views.zeroMin(Views.interval(Views.extendZero( img ), interval1)), + Views.zeroMin(Views.interval(Views.extendZero( img ), interval2)), + extension, new ArrayImgFactory(), + new FloatType(), new ArrayImgFactory(), new ComplexFloatType(), Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors())); + + PhaseCorrelationPeak2 shiftPeak = PhaseCorrelation2.getShift(pcm, + Views.zeroMin(Views.interval(Views.extendZero( img ), interval1)), + Views.zeroMin(Views.interval(Views.extendZero( img ), interval2)), + 20, 0, false, Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors())); + + long[] expected = new long[]{shiftX, shiftY}; + long[] found = new long[img.numDimensions()]; + + + + shiftPeak.getShift().localize(found); + System.out.println( Util.printCoordinates( found ) ); + + assertArrayEquals(expected, found); + + } + + @Test + public void testPCRealShift() { + + // TODO: very large shifts (nearly no overlap) lead to incorrect shift determination (as expected) + // maybe we can optimize behaviour in this situation + Img< FloatType > img = ArrayImgs.floats( 200, 200 ); + Random rnd = new Random( seed ); + + for( FloatType t : img ) + t.set( rnd.nextFloat()); + + double shiftX = -20.9; + double shiftY = 1.9; + + // to test < 0.5 px off + final double eps = 0.5; + + FinalInterval interval2 = new FinalInterval(new long[] {50, 50}); + + + + AffineRandomAccessible< FloatType, AffineGet > imgTr = RealViews.affine( Views.interpolate( Views.extendZero( img ), new NLinearInterpolatorFactory<>() ), new Translation2D( shiftX, shiftY )); + IntervalView< FloatType > img2 = Views.interval( Views.raster( imgTr ), interval2); + + int [] extension = new int[img.numDimensions()]; + Arrays.fill(extension, 10); + + RandomAccessibleInterval pcm = PhaseCorrelation2.calculatePCM(Views.zeroMin(img2), Views.zeroMin(Views.interval(img, interval2)), extension, new ArrayImgFactory(), + new FloatType(), new ArrayImgFactory(), new ComplexFloatType(), Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors())); + + PhaseCorrelationPeak2 shiftPeak = PhaseCorrelation2.getShift(pcm, Views.zeroMin(img2), Views.zeroMin(Views.interval(img, interval2)), 20, 0, true, Executors.newFixedThreadPool(Runtime.getRuntime().availableProcessors())); + + + double[] expected = new double[]{shiftX, shiftY}; + double[] found = new double[img.numDimensions()]; + + + + + shiftPeak.getSubpixelShift().localize(found); + + System.out.println( Util.printCoordinates( found ) ); + + + for (int d = 0; d < expected.length; d++) + assertTrue( Math.abs( expected[d] - found[d] ) < eps ); + + } + + +} From 4913138242fe4ba1191cfe9092ef1162c6376d1f Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Thu, 8 Jan 2026 16:46:19 -0500 Subject: [PATCH 05/16] scaled-back version that computes on 2 tiles of your choosing --- .../explorer/ViewSetupExplorerPanel.java | 4 +- .../AnalyzeOverlapCrossCorrelationPopup.java | 583 ------------------ 2 files changed, 2 insertions(+), 585 deletions(-) delete mode 100644 src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/AnalyzeOverlapCrossCorrelationPopup.java diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/ViewSetupExplorerPanel.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/ViewSetupExplorerPanel.java index 8eba44b38..b4c83708c 100644 --- a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/ViewSetupExplorerPanel.java +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/ViewSetupExplorerPanel.java @@ -82,11 +82,11 @@ import net.preibisch.mvrecon.fiji.spimdata.XmlIoSpimData2; import net.preibisch.mvrecon.fiji.spimdata.explorer.bdv.ScrollableBrightnessDialog; import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.AnalyzeErrorPopup; -import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.AnalyzeOverlapCrossCorrelationPopup; import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.ApplyTransformationPopup; import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.BDVPopup; import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.BakeManualTransformationPopup; import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.BoundingBoxPopup; +import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.ComputeCrossCorrelationPopup; import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.DeconvolutionPopup; import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.DetectInterestPointsPopup; import net.preibisch.mvrecon.fiji.spimdata.explorer.popup.DisplayFusedImagesPopup; @@ -695,6 +695,7 @@ public ArrayList< ExplorerWindowSetable > initPopups() popups.add( new BDVPopup() ); popups.add( new DisplayRawImagesPopup() ); popups.add( new DisplayFusedImagesPopup() ); + popups.add( new ComputeCrossCorrelationPopup() ); popups.add( new VisualizeNonRigid() ); popups.add( new MaxProjectPopup() ); popups.add( new Separator() ); @@ -725,7 +726,6 @@ public ArrayList< ExplorerWindowSetable > initPopups() popups.add( new LabelPopUp( " Interest Points" ) ); popups.add( new InterestPointsExplorerPopup() ); popups.add( new AnalyzeErrorPopup() ); - popups.add( new AnalyzeOverlapCrossCorrelationPopup() ); popups.add( new RemoveDetectionsPopup() ); popups.add( new VisualizeDetectionsPopup() ); popups.add( new Separator() ); diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/AnalyzeOverlapCrossCorrelationPopup.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/AnalyzeOverlapCrossCorrelationPopup.java deleted file mode 100644 index 8a68576c7..000000000 --- a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/AnalyzeOverlapCrossCorrelationPopup.java +++ /dev/null @@ -1,583 +0,0 @@ -/*- - * #%L - * Software for the reconstruction of multi-view microscopic acquisitions - * like Selective Plane Illumination Microscopy (SPIM) Data. - * %% - * Copyright (C) 2012 - 2025 Multiview Reconstruction developers. - * %% - * This program is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as - * published by the Free Software Foundation, either version 2 of the - * License, or (at your option) any later version. - * - * This program is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public - * License along with this program. If not, see - * . - * #L% - */ -package net.preibisch.mvrecon.fiji.spimdata.explorer.popup; - -import java.awt.event.ActionEvent; -import java.awt.event.ActionListener; -import java.util.ArrayList; -import java.util.Collection; -import java.util.Collections; -import java.util.Comparator; -import java.util.HashMap; -import java.util.List; -import java.util.Map; -import java.util.concurrent.Callable; -import java.util.concurrent.ExecutorService; -import java.util.concurrent.Executors; -import java.util.concurrent.Future; - -import javax.swing.JMenuItem; - -import ij.gui.GenericDialog; -import mpicbg.spim.data.generic.sequence.AbstractSequenceDescription; -import mpicbg.spim.data.generic.sequence.BasicImgLoader; -import mpicbg.spim.data.generic.sequence.BasicViewDescription; -import mpicbg.spim.data.sequence.MultiResolutionImgLoader; -import mpicbg.spim.data.registration.ViewRegistration; -import mpicbg.spim.data.registration.ViewRegistrations; -import mpicbg.spim.data.sequence.ViewId; -import net.imglib2.Cursor; -import net.imglib2.FinalInterval; -import net.imglib2.Interval; -import net.imglib2.RandomAccessibleInterval; -import net.imglib2.realtransform.AffineTransform3D; -import net.imglib2.type.numeric.real.FloatType; -import net.imglib2.view.Views; -import net.preibisch.legacy.io.IOFunctions; -import net.preibisch.mvrecon.fiji.spimdata.SpimData2; -import net.preibisch.mvrecon.fiji.spimdata.boundingbox.BoundingBox; -import net.preibisch.mvrecon.fiji.spimdata.explorer.ExplorerWindow; -import net.preibisch.mvrecon.fiji.spimdata.interestpoints.CorrespondingInterestPoints; -import net.preibisch.mvrecon.fiji.spimdata.interestpoints.ViewInterestPointLists; -import net.preibisch.mvrecon.process.boundingbox.BoundingBoxMaximalGroupOverlap; -import net.preibisch.mvrecon.process.downsampling.DownsampleTools; -import net.preibisch.mvrecon.process.fusion.transformed.FusedRandomAccessibleInterval; -import net.preibisch.mvrecon.process.fusion.transformed.TransformView; -import net.preibisch.mvrecon.process.fusion.transformed.TransformVirtual; -import net.preibisch.mvrecon.process.interestpointregistration.pairwise.constellation.grouping.Group; -import net.preibisch.mvrecon.process.phasecorrelation.PhaseCorrelation2Util; - -public class AnalyzeOverlapCrossCorrelationPopup extends JMenuItem implements ExplorerWindowSetable -{ - private static final long serialVersionUID = 1L; - - ExplorerWindow< ? > panel; - - public static int defaultDownsamplingChoiceIndex = 0; // First available resolution - public static int defaultMinCorrespondences = 10; - - @Override - public JMenuItem setExplorerWindow( ExplorerWindow< ? > panel ) - { - this.panel = panel; - return this; - } - - public AnalyzeOverlapCrossCorrelationPopup() - { - super( "Analyze Overlap Cross-Correlation ..." ); - this.addActionListener( new MyActionListener() ); - } - - public class MyActionListener implements ActionListener - { - @Override - public void actionPerformed( ActionEvent e ) - { - if ( panel == null ) - { - IOFunctions.println( "Panel not set for " + this.getClass().getSimpleName() ); - return; - } - - if ( !SpimData2.class.isInstance( panel.getSpimData() ) ) - { - IOFunctions.println( "Only supported for SpimData2 objects: " + this.getClass().getSimpleName() ); - return; - } - - final SpimData2 spimData = (SpimData2) panel.getSpimData(); - final List< ViewId > viewIds = ApplyTransformationPopup.getSelectedViews( panel ); - - if ( viewIds.isEmpty() ) - { - IOFunctions.println( "No views selected." ); - return; - } - - // Get available downsampling levels from dataset metadata - final String[] downsamplingChoices = DownsampleTools.availableDownsamplings( spimData, viewIds.get( 0 ) ); - - // Ask user for parameters - final GenericDialog gd = new GenericDialog( "Cross-Correlation Analysis Parameters" ); - gd.addChoice( "Downsampling", downsamplingChoices, downsamplingChoices[ Math.min( defaultDownsamplingChoiceIndex, downsamplingChoices.length - 1 ) ] ); - gd.addNumericField( "Minimum_correspondences (0 = analyze all overlaps)", defaultMinCorrespondences, 0 ); - gd.addMessage( "Note: Analysis will be performed on the maximal bounding box\nof the overlap between each pair of selected tiles." ); - gd.addMessage( "Set minimum correspondences to 0 to analyze all overlapping pairs\nwithout requiring pre-computed interest points." ); - - gd.showDialog(); - - if ( gd.wasCanceled() ) - return; - - defaultDownsamplingChoiceIndex = gd.getNextChoiceIndex(); - defaultMinCorrespondences = (int) gd.getNextNumber(); - - - new Thread( new Runnable() - { - @Override - public void run() - { - analyzeCrossCorrelations( spimData, viewIds, defaultDownsamplingChoiceIndex, defaultMinCorrespondences ); - } - } ).start(); - } - } - - public static class OverlapResult - { - public ViewId viewA; - public ViewId viewB; - public int numCorrespondences; - public double crossCorrelation; - public double avgIntensityA; - public double avgIntensityB; - public String label; - public int downsamplingIndex; - - public OverlapResult( ViewId viewA, ViewId viewB, String label, int numCorr, double cc, double avgA, double avgB, int dsIndex ) - { - this.viewA = viewA; - this.viewB = viewB; - this.label = label; - this.numCorrespondences = numCorr; - this.crossCorrelation = cc; - this.avgIntensityA = avgA; - this.avgIntensityB = avgB; - this.downsamplingIndex = dsIndex; - } - } - - public static void analyzeCrossCorrelations( - final SpimData2 spimData, - final List< ViewId > viewIds, - final int downsamplingIndex, - final int minCorrespondences ) - { - IOFunctions.println( "Starting Cross-Correlation Analysis..." ); - IOFunctions.println( "Downsampling index: " + downsamplingIndex ); - IOFunctions.println( "Minimum correspondences: " + minCorrespondences ); - IOFunctions.println( "Number of views: " + viewIds.size() ); - - final ExecutorService service = Executors.newFixedThreadPool( Runtime.getRuntime().availableProcessors() ); - final List< OverlapResult > results = new ArrayList<>(); - - final String selectedLabel; - - // Only load interest points if minCorrespondences > 0 - if ( minCorrespondences > 0 ) - { - // Load all interest points and correspondences - IOFunctions.println( "Loading interest points and correspondences..." ); - final Map< ViewId, ViewInterestPointLists > interestPoints = new HashMap<>(); - - for ( final ViewId viewId : viewIds ) - { - final ViewInterestPointLists vipl = spimData.getViewInterestPoints().getViewInterestPointLists( viewId ); - interestPoints.put( viewId, vipl ); - } - - // Get available labels (interest point types) - final Map< String, Integer > labelCounts = new HashMap<>(); - for ( final ViewId viewId : viewIds ) - { - final ViewInterestPointLists vipl = interestPoints.get( viewId ); - for ( final String label : vipl.getHashMap().keySet() ) - { - labelCounts.put( label, labelCounts.getOrDefault( label, 0 ) + 1 ); - } - } - - // Find the most common label - String tempLabel = null; - int maxCount = 0; - for ( final Map.Entry< String, Integer > entry : labelCounts.entrySet() ) - { - if ( entry.getValue() > maxCount ) - { - maxCount = entry.getValue(); - tempLabel = entry.getKey(); - } - } - - if ( tempLabel == null ) - { - IOFunctions.println( "ERROR: No interest points found in selected views." ); - service.shutdown(); - return; - } - - selectedLabel = tempLabel; - IOFunctions.println( "Using interest point label: " + selectedLabel ); - } - else - { - IOFunctions.println( "Minimum correspondences = 0, analyzing all overlapping pairs without checking interest points." ); - selectedLabel = null; - } - - // Analyze all pairs in parallel - final List< Future< OverlapResult > > futures = new ArrayList<>(); - - for ( int i = 0; i < viewIds.size() - 1; ++i ) - { - for ( int j = i + 1; j < viewIds.size(); ++j ) - { - final ViewId viewA = viewIds.get( i ); - final ViewId viewB = viewIds.get( j ); - - // Submit pair analysis as a parallel task - futures.add( service.submit( new Callable< OverlapResult >() - { - @Override - public OverlapResult call() - { - return analyzePair( spimData, viewA, viewB, selectedLabel, - downsamplingIndex, minCorrespondences, service ); - } - } ) ); - } - } - - // Collect results from all tasks - for ( Future< OverlapResult > future : futures ) - { - try - { - final OverlapResult result = future.get(); - if ( result != null ) - results.add( result ); - } - catch ( Exception e ) - { - IOFunctions.println( "Error processing pair: " + e.getMessage() ); - e.printStackTrace(); - } - } - - service.shutdown(); - - // Sort by cross-correlation (ascending) - Collections.sort( results, new Comparator< OverlapResult >() - { - @Override - public int compare( OverlapResult o1, OverlapResult o2 ) - { - return Double.compare( o1.crossCorrelation, o2.crossCorrelation ); - } - } ); - - // Print results - IOFunctions.println( "\n========================================" ); - IOFunctions.println( "Cross-Correlation Analysis Results" ); - IOFunctions.println( "========================================" ); - if ( minCorrespondences > 0 ) - { - IOFunctions.println( String.format( "%-30s %-15s %-15s %-15s %-15s %-15s", - "Tile Pair", "# Corresp.", "Cross-Corr", "Avg Int. A", "Avg Int. B", "Label" ) ); - } - else - { - IOFunctions.println( String.format( "%-30s %-15s %-15s %-15s", - "Tile Pair", "Cross-Corr", "Avg Int. A", "Avg Int. B" ) ); - } - IOFunctions.println( "----------------------------------------" ); - - for ( final OverlapResult result : results ) - { - final String pairName = getTileName( spimData, result.viewA ) + " <-> " + getTileName( spimData, result.viewB ); - if ( minCorrespondences > 0 ) - { - IOFunctions.println( String.format( "%-30s %-15d %-15.4f %-15.2f %-15.2f %-15s", - pairName, - result.numCorrespondences, - result.crossCorrelation, - result.avgIntensityA, - result.avgIntensityB, - result.label != null ? result.label : "N/A" ) ); - } - else - { - IOFunctions.println( String.format( "%-30s %-15.4f %-15.2f %-15.2f", - pairName, - result.crossCorrelation, - result.avgIntensityA, - result.avgIntensityB ) ); - } - } - - IOFunctions.println( "========================================\n" ); - } - - private static OverlapResult analyzePair( - final SpimData2 spimData, - final ViewId viewA, - final ViewId viewB, - final String label, - final int downsamplingIndex, - final int minCorrespondences, - final ExecutorService service ) - { - int numCorr = 0; - - // Only check correspondences if minCorrespondences > 0 - if ( minCorrespondences > 0 ) - { - // Get interest points - final ViewInterestPointLists viplA = spimData.getViewInterestPoints().getViewInterestPointLists( viewA ); - final ViewInterestPointLists viplB = spimData.getViewInterestPoints().getViewInterestPointLists( viewB ); - - if ( viplA.getInterestPointList( label ) == null || viplB.getInterestPointList( label ) == null ) - return null; - - // Count correspondences between A and B - final Collection< CorrespondingInterestPoints > corrsA = viplA.getInterestPointList( label ).getCorrespondingInterestPointsCopy(); - - for ( final CorrespondingInterestPoints cip : corrsA ) - { - if ( cip.getCorrespondingViewId().equals( viewB ) && cip.getCorrespodingLabel().equals( label ) ) - numCorr++; - } - - if ( numCorr < minCorrespondences ) - { - IOFunctions.println( getTileName( spimData, viewA ) + " <-> " + getTileName( spimData, viewB ) + - ": " + numCorr + " correspondences (< " + minCorrespondences + "), skipping." ); - return null; - } - - IOFunctions.println( "Processing: " + getTileName( spimData, viewA ) + " <-> " + getTileName( spimData, viewB ) + - " (" + numCorr + " correspondences)" ); - } - else - { - IOFunctions.println( "Processing: " + getTileName( spimData, viewA ) + " <-> " + getTileName( spimData, viewB ) ); - } - - // Compute overlap bounding box - final List< List< ViewId > > viewGroups = new ArrayList<>(); - final List< ViewId > groupA = new ArrayList<>(); - groupA.add( viewA ); - final List< ViewId > groupB = new ArrayList<>(); - groupB.add( viewB ); - viewGroups.add( groupA ); - viewGroups.add( groupB ); - - final BoundingBoxMaximalGroupOverlap< ViewId > bbDet = new BoundingBoxMaximalGroupOverlap<>( viewGroups, spimData ); - final BoundingBox bbOverlap = bbDet.estimate( "overlap" ); - - if ( bbOverlap == null ) - { - IOFunctions.println( " No overlap found, skipping." ); - return null; - } - - // Use downsampling index to open image at specified resolution level - // No need to specify downsampling factors - the index handles it - - try - { - final RandomAccessibleInterval< FloatType > imgA = loadImageInBB( - spimData.getSequenceDescription(), - spimData.getViewRegistrations(), - viewA, - bbOverlap, - downsamplingIndex ); - - final RandomAccessibleInterval< FloatType > imgB = loadImageInBB( - spimData.getSequenceDescription(), - spimData.getViewRegistrations(), - viewB, - bbOverlap, - downsamplingIndex ); - - // Compute downsampling factor for adaptive sampling - final double downsamplingFactor = Math.pow( 2, downsamplingIndex ); - - // Compute average intensities and cross-correlation - final double avgA = computeAverageIntensity( imgA, downsamplingFactor ); - final double avgB = computeAverageIntensity( imgB, downsamplingFactor ); - final double cc = PhaseCorrelation2Util.getCorrelation( imgA, imgB ); - - IOFunctions.println( " Cross-correlation: " + cc + ", Avg intensities: " + avgA + ", " + avgB ); - - return new OverlapResult( viewA, viewB, label, numCorr, cc, avgA, avgB, downsamplingIndex ); - } - catch ( Exception e ) - { - IOFunctions.println( " Error processing pair: " + e.getMessage() ); - e.printStackTrace(); - return null; - } - } - - private static RandomAccessibleInterval< FloatType > loadImageInBB( - final AbstractSequenceDescription< ?, ?, ? > sd, - final ViewRegistrations vrs, - final ViewId viewId, - final Interval boundingBox, - final int downsamplingIndex ) - { - final BasicImgLoader imgLoader = sd.getImgLoader(); - final ViewRegistration vr = vrs.getViewRegistration( viewId ); - - // Synchronize access to ViewRegistration to prevent race conditions when running in parallel - final AffineTransform3D model; - synchronized ( vr ) - { - vr.updateModel(); - model = vr.getModel().copy(); - } - - // Load image at the specific resolution level chosen by user - RandomAccessibleInterval inputImg; - - if ( imgLoader instanceof MultiResolutionImgLoader ) - { - final MultiResolutionImgLoader mrImgLoader = (MultiResolutionImgLoader) imgLoader; - final int setupId = viewId.getViewSetupId(); - final int timepointId = viewId.getTimePointId(); - - // Load image at user-specified level - inputImg = mrImgLoader.getSetupImgLoader( setupId ).getImage( timepointId, downsamplingIndex ); - - // Get and apply the mipmap transform for this level - final AffineTransform3D mipmapTransform = mrImgLoader.getSetupImgLoader( setupId ).getMipmapTransforms()[ downsamplingIndex ]; - model.concatenate( mipmapTransform ); - - // Use bounding box in world coordinates (do NOT scale it) - // The model transform (with mipmap concatenated) handles the coordinate mapping - // NOTE: The result image will have the same dimensions regardless of downsampling level - // because the bounding box is in world coordinates. The speedup from downsampling comes - // from faster access to the pre-downsampled input image, not from smaller output dimensions. - final Interval bb = new FinalInterval( boundingBox ); - - // Transform view to bounding box - return TransformView.transformView( inputImg, model, bb, 0, 1 ); - } - else - { - // Fallback for non-multiresolution loaders: load full resolution - inputImg = imgLoader.getSetupImgLoader( viewId.getViewSetupId() ).getImage( viewId.getTimePointId() ); - - // Compute downsampling factor: 2^index - final double downsamplingFactor = Math.pow( 2, downsamplingIndex ); - final double[] downsamplingFactors = new double[ 3 ]; - for ( int d = 0; d < 3; ++d ) - downsamplingFactors[ d ] = downsamplingFactor; - - // Scale bounding box - final Interval bbSc = TransformVirtual.scaleBoundingBox( new FinalInterval( boundingBox ), inverse( downsamplingFactors ) ); - - // Scale transform - TransformVirtual.scaleTransform( model, inverse( downsamplingFactors ) ); - - // Transform view to bounding box - return TransformView.transformView( inputImg, model, bbSc, 0, 1 ); - } - } - - private static double[] inverse( double[] in ) - { - final double[] res = new double[ in.length ]; - for ( int i = 0; i < in.length; i++ ) - res[ i ] = 1.0 / in[ i ]; - return res; - } - - private static double computeAverageIntensity( final RandomAccessibleInterval< FloatType > img, final double downsamplingFactor ) - { - // Sample every Nth pixel in each dimension for speed - // Scale sampling step with downsampling level: at 8x downsampling, sample 2x less frequently than at 4x - final int baseStep = 4; - final int step = Math.max( 1, (int) (baseStep * Math.sqrt( downsamplingFactor )) ); - - double sum = 0; - long count = 0; - - final long[] min = new long[ img.numDimensions() ]; - final long[] max = new long[ img.numDimensions() ]; - img.min( min ); - img.max( max ); - - final net.imglib2.RandomAccess< FloatType > ra = img.randomAccess(); - final long[] pos = new long[ img.numDimensions() ]; - - // Sample pixels at regular intervals - if ( img.numDimensions() == 3 ) - { - for ( long z = min[2]; z <= max[2]; z += step ) - { - pos[2] = z; - for ( long y = min[1]; y <= max[1]; y += step ) - { - pos[1] = y; - for ( long x = min[0]; x <= max[0]; x += step ) - { - pos[0] = x; - ra.setPosition( pos ); - sum += ra.get().get(); - count++; - } - } - } - } - else if ( img.numDimensions() == 2 ) - { - for ( long y = min[1]; y <= max[1]; y += step ) - { - pos[1] = y; - for ( long x = min[0]; x <= max[0]; x += step ) - { - pos[0] = x; - ra.setPosition( pos ); - sum += ra.get().get(); - count++; - } - } - } - else - { - // Fallback for other dimensions: use cursor (slower) - final Cursor< FloatType > cursor = Views.iterable( img ).cursor(); - while ( cursor.hasNext() ) - { - cursor.fwd(); - sum += cursor.get().get(); - count++; - } - } - - return count > 0 ? sum / count : 0; - } - - private static String getTileName( final SpimData2 spimData, final ViewId viewId ) - { - final BasicViewDescription< ? > vd = spimData.getSequenceDescription().getViewDescriptions().get( viewId ); - if ( vd != null ) - return Group.pvid( viewId ); - else - return viewId.toString(); - } -} From a7696d3480f5f18ff8d0b20cbc31750ce377d4d6 Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Fri, 9 Jan 2026 11:07:10 -0500 Subject: [PATCH 06/16] added downsampling support --- .../popup/ComputeCrossCorrelationPopup.java | 292 ++++++++++++++++++ 1 file changed, 292 insertions(+) create mode 100644 src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java new file mode 100644 index 000000000..6790a95a0 --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java @@ -0,0 +1,292 @@ +/*- + * #%L + * Software for the reconstruction of multi-view microscopic acquisitions + * like Selective Plane Illumination Microscopy (SPIM) Data. + * %% + * Copyright (C) 2012 - 2026 Multiview Reconstruction developers. + * %% + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as + * published by the Free Software Foundation, either version 2 of the + * License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public + * License along with this program. If not, see + * . + * #L% + */ +package net.preibisch.mvrecon.fiji.spimdata.explorer.popup; + +import java.awt.event.ActionEvent; +import java.awt.event.ActionListener; +import java.util.Date; +import java.util.List; + +import javax.swing.JComponent; +import javax.swing.JMenuItem; + +import ij.gui.GenericDialog; + +import mpicbg.spim.data.generic.sequence.BasicImgLoader; +import mpicbg.spim.data.generic.sequence.BasicViewSetup; +import mpicbg.spim.data.registration.ViewRegistration; +import mpicbg.spim.data.sequence.ViewId; +import net.imglib2.FinalInterval; +import net.imglib2.Interval; +import net.imglib2.Point; +import net.imglib2.RandomAccessibleInterval; +import net.imglib2.RealInterval; +import net.imglib2.realtransform.AffineTransform3D; +import net.imglib2.type.numeric.real.FloatType; +import net.imglib2.util.Util; +import net.imglib2.view.Views; +import net.preibisch.legacy.io.IOFunctions; +import net.preibisch.mvrecon.fiji.spimdata.SpimData2; +import net.preibisch.mvrecon.fiji.spimdata.explorer.ExplorerWindow; +import net.preibisch.mvrecon.fiji.spimdata.explorer.GroupedRowWindow; +import net.preibisch.mvrecon.process.downsampling.DownsampleTools; +import net.preibisch.mvrecon.process.interestpointregistration.pairwise.constellation.overlap.SimpleBoundingBoxOverlap; +import net.preibisch.mvrecon.process.phasecorrelation.PhaseCorrelationPeak2; + +public class ComputeCrossCorrelationPopup extends JMenuItem implements ExplorerWindowSetable +{ + private static final long serialVersionUID = 1L; + ExplorerWindow panel; + + public ComputeCrossCorrelationPopup() + { + super("Compute Cross-Correlation for Selected Tiles"); + this.addActionListener(new MyActionListener()); + } + + @Override + public JComponent setExplorerWindow(ExplorerWindow panel) + { + this.panel = panel; + return this; + } + + private class MyActionListener implements ActionListener + { + @Override + public void actionPerformed(ActionEvent e) + { + if (panel == null) + { + IOFunctions.println("Panel not set for " + this.getClass().getSimpleName()); + return; + } + + new Thread(new Runnable() + { + @Override + public void run() + { + try + { + // Get SpimData + final SpimData2 spimData = panel.getSpimData(); + + // Step 1: Validate selection - exactly 2 tiles + List> selectedGroups = ((GroupedRowWindow) panel).selectedRowsViewIdGroups(); + + if (selectedGroups.size() != 2) + { + IOFunctions.println(new Date(System.currentTimeMillis()) + + ": ERROR: Please select exactly 2 tiles. Currently selected: " + + selectedGroups.size()); + return; + } + + // Filter missing views + for (List group : selectedGroups) + { + SpimData2.filterMissingViews(spimData, group); + } + + // Check that each group has at least one view + if (selectedGroups.get(0).isEmpty() || selectedGroups.get(1).isEmpty()) + { + IOFunctions.println(new Date(System.currentTimeMillis()) + + ": ERROR: Selected groups contain no valid views"); + return; + } + + // Get first ViewId from each group + ViewId viewId1 = selectedGroups.get(0).get(0); + ViewId viewId2 = selectedGroups.get(1).get(0); + +// IOFunctions.println(new Date(System.currentTimeMillis()) + +// ": Computing cross-correlation between tiles..."); +// IOFunctions.println(" Tile 1: ViewId(" + viewId1.getTimePointId() + ", " + +// viewId1.getViewSetupId() + ")"); +// IOFunctions.println(" Tile 2: ViewId(" + viewId2.getTimePointId() + ", " + +// viewId2.getViewSetupId() + ")"); + + // Step 2: Calculate overlap using SimpleBoundingBoxOverlap + SimpleBoundingBoxOverlap overlapDetection = new SimpleBoundingBoxOverlap<>(spimData); + + RealInterval globalOverlap = overlapDetection.getOverlapInterval(viewId1, viewId2); + + if (globalOverlap == null) + { + IOFunctions.println(new Date(System.currentTimeMillis()) + + ": ERROR: Selected tiles do not overlap"); + return; + } + + // excess debug +// // Print global overlap region +// StringBuilder sb = new StringBuilder("("); +// for (int d = 0; d < globalOverlap.numDimensions(); d++) +// { +// if (d > 0) sb.append(", "); +// sb.append(String.format("%.1f", globalOverlap.realMin(d))); +// sb.append(" - "); +// sb.append(String.format("%.1f", globalOverlap.realMax(d))); +// } +// sb.append(")"); +// IOFunctions.println(" Global overlap region: " + sb); + + // Step 3: Get bounding boxes and convert overlap to local coordinates + BasicViewSetup vs1 = spimData.getSequenceDescription().getViewSetups().get(viewId1.getViewSetupId()); + BasicViewSetup vs2 = spimData.getSequenceDescription().getViewSetups().get(viewId2.getViewSetupId()); + ViewRegistration vr1 = spimData.getViewRegistrations().getViewRegistration(viewId1); + ViewRegistration vr2 = spimData.getViewRegistrations().getViewRegistration(viewId2); + + vr1.updateModel(); + vr2.updateModel(); + + RealInterval bbox1 = SimpleBoundingBoxOverlap.getBoundingBoxReal(vs1, vr1); + RealInterval bbox2 = SimpleBoundingBoxOverlap.getBoundingBoxReal(vs2, vr2); + + // Convert global overlap to local coordinates + double[] localMin1 = new double[globalOverlap.numDimensions()]; + double[] localMax1 = new double[globalOverlap.numDimensions()]; + double[] localMin2 = new double[globalOverlap.numDimensions()]; + double[] localMax2 = new double[globalOverlap.numDimensions()]; + + for (int d = 0; d < globalOverlap.numDimensions(); d++) + { + localMin1[d] = globalOverlap.realMin(d) - bbox1.realMin(d); + localMax1[d] = globalOverlap.realMax(d) - bbox1.realMin(d); + localMin2[d] = globalOverlap.realMin(d) - bbox2.realMin(d); + localMax2[d] = globalOverlap.realMax(d) - bbox2.realMin(d); + } + + // Convert to raster coordinates with bounds checking + long[] dims1 = new long[globalOverlap.numDimensions()]; + long[] dims2 = new long[globalOverlap.numDimensions()]; + vs1.getSize().dimensions(dims1); + vs2.getSize().dimensions(dims2); + + long[] rasterMin1 = new long[globalOverlap.numDimensions()]; + long[] rasterMax1 = new long[globalOverlap.numDimensions()]; + long[] rasterMin2 = new long[globalOverlap.numDimensions()]; + long[] rasterMax2 = new long[globalOverlap.numDimensions()]; + + for (int d = 0; d < globalOverlap.numDimensions(); d++) + { + rasterMin1[d] = Math.max(0, (long) Math.ceil(localMin1[d])); + rasterMax1[d] = Math.min(dims1[d] - 1, (long) Math.floor(localMax1[d])); + rasterMin2[d] = Math.max(0, (long) Math.ceil(localMin2[d])); + rasterMax2[d] = Math.min(dims2[d] - 1, (long) Math.floor(localMax2[d])); + } + + // Step 4: Ask user for downsampling level + // Get available downsampling levels from the data pyramid + String[] availableDownsamplings = DownsampleTools.availableDownsamplings(spimData, viewId1); + + final GenericDialog gd = new GenericDialog("Cross-Correlation Parameters"); + gd.addChoice("Downsampling", availableDownsamplings, availableDownsamplings[0]); + gd.addMessage("Choose downsampling level:"); + gd.addMessage("1, 1, 1 = full resolution (slower, more accurate)"); + gd.addMessage("Higher values = faster computation (less accurate)"); + gd.showDialog(); + + if (gd.wasCanceled()) + { + IOFunctions.println(new Date(System.currentTimeMillis()) + + ": Cross-correlation computation cancelled by user"); + return; + } + + final String downsamplingChoice = gd.getNextChoice(); + final long[] downsampleFactors = DownsampleTools.parseDownsampleChoice(downsamplingChoice); + + IOFunctions.println(" Using downsampling: " + downsamplingChoice); + + // Adjust overlap intervals for downsampling + long[] dsRasterMin1 = new long[rasterMin1.length]; + long[] dsRasterMax1 = new long[rasterMax1.length]; + long[] dsRasterMin2 = new long[rasterMin2.length]; + long[] dsRasterMax2 = new long[rasterMax2.length]; + + for (int d = 0; d < rasterMin1.length; d++) + { + dsRasterMin1[d] = rasterMin1[d] / downsampleFactors[d]; + dsRasterMax1[d] = rasterMax1[d] / downsampleFactors[d]; + dsRasterMin2[d] = rasterMin2[d] / downsampleFactors[d]; + dsRasterMax2[d] = rasterMax2[d] / downsampleFactors[d]; + } + + Interval dsInterval1 = new FinalInterval(dsRasterMin1, dsRasterMax1); + Interval dsInterval2 = new FinalInterval(dsRasterMin2, dsRasterMax2); + + // Step 5: Load image data at specified downsampling level using pyramid + final BasicImgLoader imgLoader = spimData.getSequenceDescription().getImgLoader(); + + // Open images at the specified downsampling level using pre-computed pyramid + net.imglib2.util.Pair opened1 = + DownsampleTools.openAndDownsample(imgLoader, viewId1, downsampleFactors, false); + net.imglib2.util.Pair opened2 = + DownsampleTools.openAndDownsample(imgLoader, viewId2, downsampleFactors, false); + + RandomAccessibleInterval img1 = opened1.getA(); + RandomAccessibleInterval img2 = opened2.getA(); + + // Extract overlap regions and convert to FloatType + RandomAccessibleInterval overlap1 = Views.zeroMin( + Views.interval(Views.zeroMin(img1), dsInterval1)); + RandomAccessibleInterval overlap2 = Views.zeroMin( + Views.interval(Views.zeroMin(img2), dsInterval2)); + + // Step 6: Compute cross-correlation directly + IOFunctions.println(" Computing cross-correlation..."); + + // Create a peak with zero shift to compute correlation at current alignment + Point zeroShift = new Point(overlap1.numDimensions()); + PhaseCorrelationPeak2 peak = new PhaseCorrelationPeak2(zeroShift, 0.0); + + // Set the shift field (required by calculateCrossCorr) + peak.setShift(zeroShift); + + // Calculate cross-correlation at zero shift (current alignment) + peak.calculateCrossCorr(overlap1, overlap2); + + // Step 7: Log results + double correlationCoefficient = peak.getCrossCorr(); + long nPixels = peak.getnPixel(); + + IOFunctions.println(new Date(System.currentTimeMillis()) + + String.format(" %d-%d <> %d-%d", viewId1.getTimePointId(), viewId1.getViewSetupId(), + viewId2.getTimePointId(), viewId2.getViewSetupId()) + + String.format(": r=%.4f", correlationCoefficient)); + + } + catch (Exception ex) + { + IOFunctions.println(new Date(System.currentTimeMillis()) + + ": ERROR: Exception during correlation computation: " + ex.getMessage()); + ex.printStackTrace(); + } + } + }).start(); + } + } +} From 29623ed1a33fcdfdd4cc1929ded13f406bdd2913 Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Fri, 9 Jan 2026 11:31:52 -0500 Subject: [PATCH 07/16] extended to all overlapping pairs with option to only calculate on selected tiles --- .../popup/ComputeCrossCorrelationPopup.java | 425 +++++++++++------- 1 file changed, 272 insertions(+), 153 deletions(-) diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java index 6790a95a0..9f08bbf2c 100644 --- a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java @@ -24,6 +24,7 @@ import java.awt.event.ActionEvent; import java.awt.event.ActionListener; +import java.util.ArrayList; import java.util.Date; import java.util.List; @@ -60,7 +61,7 @@ public class ComputeCrossCorrelationPopup extends JMenuItem implements ExplorerW public ComputeCrossCorrelationPopup() { - super("Compute Cross-Correlation for Selected Tiles"); + super("Compute Cross-Correlation for Overlapping Tiles"); this.addActionListener(new MyActionListener()); } @@ -92,201 +93,319 @@ public void run() // Get SpimData final SpimData2 spimData = panel.getSpimData(); - // Step 1: Validate selection - exactly 2 tiles + // Get selected groups for potential use List> selectedGroups = ((GroupedRowWindow) panel).selectedRowsViewIdGroups(); - if (selectedGroups.size() != 2) + // Get available downsampling levels + List allViewIds = new ArrayList<>(spimData.getSequenceDescription().getViewDescriptions().keySet()); + if (allViewIds.isEmpty()) { IOFunctions.println(new Date(System.currentTimeMillis()) + - ": ERROR: Please select exactly 2 tiles. Currently selected: " + - selectedGroups.size()); + ": ERROR: No views found in dataset"); return; } - // Filter missing views - for (List group : selectedGroups) - { - SpimData2.filterMissingViews(spimData, group); - } + ViewId sampleViewId = allViewIds.get(0); + String[] availableDownsamplings = DownsampleTools.availableDownsamplings(spimData, sampleViewId); + + // Show dialog with checkbox and downsampling choice + final GenericDialog gd = new GenericDialog("Cross-Correlation Parameters"); + gd.addCheckbox("Only compute for selected tiles (requires at least 2 selected)", false); + gd.addChoice("Downsampling", availableDownsamplings, availableDownsamplings[0]); + gd.addMessage("Choose downsampling level:"); + gd.addMessage("1, 1, 1 = full resolution (slower, more accurate)"); + gd.addMessage("Higher values = faster computation (less accurate)"); + gd.showDialog(); - // Check that each group has at least one view - if (selectedGroups.get(0).isEmpty() || selectedGroups.get(1).isEmpty()) + if (gd.wasCanceled()) { IOFunctions.println(new Date(System.currentTimeMillis()) + - ": ERROR: Selected groups contain no valid views"); + ": Cross-correlation computation cancelled by user"); return; } - // Get first ViewId from each group - ViewId viewId1 = selectedGroups.get(0).get(0); - ViewId viewId2 = selectedGroups.get(1).get(0); - -// IOFunctions.println(new Date(System.currentTimeMillis()) + -// ": Computing cross-correlation between tiles..."); -// IOFunctions.println(" Tile 1: ViewId(" + viewId1.getTimePointId() + ", " + -// viewId1.getViewSetupId() + ")"); -// IOFunctions.println(" Tile 2: ViewId(" + viewId2.getTimePointId() + ", " + -// viewId2.getViewSetupId() + ")"); + final boolean selectedOnly = gd.getNextBoolean(); + final String downsamplingChoice = gd.getNextChoice(); + final long[] downsampleFactors = DownsampleTools.parseDownsampleChoice(downsamplingChoice); - // Step 2: Calculate overlap using SimpleBoundingBoxOverlap - SimpleBoundingBoxOverlap overlapDetection = new SimpleBoundingBoxOverlap<>(spimData); + IOFunctions.println(new Date(System.currentTimeMillis()) + + ": Computing cross-correlations..."); + IOFunctions.println(" Using downsampling: " + downsamplingChoice); - RealInterval globalOverlap = overlapDetection.getOverlapInterval(viewId1, viewId2); + List pairsToProcess = new ArrayList<>(); - if (globalOverlap == null) + if (selectedOnly) { - IOFunctions.println(new Date(System.currentTimeMillis()) + - ": ERROR: Selected tiles do not overlap"); - return; + // Validate at least 2 tiles selected + if (selectedGroups.size() < 2) + { + IOFunctions.println(new Date(System.currentTimeMillis()) + + ": ERROR: Please select at least 2 tiles. Currently selected: " + + selectedGroups.size()); + return; + } + + // Filter missing views and collect valid ViewIds from selected groups + List selectedViewIds = new ArrayList<>(); + for (List group : selectedGroups) + { + SpimData2.filterMissingViews(spimData, group); + if (!group.isEmpty()) + { + // Get first ViewId from each group + selectedViewIds.add(group.get(0)); + } + } + + if (selectedViewIds.size() < 2) + { + IOFunctions.println(new Date(System.currentTimeMillis()) + + ": ERROR: Selected groups contain fewer than 2 valid views"); + return; + } + + IOFunctions.println(" Processing " + selectedViewIds.size() + " selected views"); + + SimpleBoundingBoxOverlap overlapDetection = new SimpleBoundingBoxOverlap<>(spimData); + + // Find all overlapping pairs among selected views + for (int i = 0; i < selectedViewIds.size(); i++) + { + for (int j = i + 1; j < selectedViewIds.size(); j++) + { + ViewId viewId1 = selectedViewIds.get(i); + ViewId viewId2 = selectedViewIds.get(j); + + RealInterval overlap = overlapDetection.getOverlapInterval(viewId1, viewId2); + if (overlap != null) + { + pairsToProcess.add(new ViewId[] { viewId1, viewId2 }); + } + } + } + + IOFunctions.println(" Found " + pairsToProcess.size() + " overlapping pairs among selected views"); + } + else + { + // Get all views and find overlapping pairs + List viewIds = new ArrayList<>(); + for (ViewId viewId : allViewIds) + { + if (!spimData.getSequenceDescription().getMissingViews().getMissingViews().contains(viewId)) + { + viewIds.add(viewId); + } + } + + IOFunctions.println(" Found " + viewIds.size() + " valid views"); + + SimpleBoundingBoxOverlap overlapDetection = new SimpleBoundingBoxOverlap<>(spimData); + + // Find all overlapping pairs + for (int i = 0; i < viewIds.size(); i++) + { + for (int j = i + 1; j < viewIds.size(); j++) + { + ViewId viewId1 = viewIds.get(i); + ViewId viewId2 = viewIds.get(j); + + RealInterval overlap = overlapDetection.getOverlapInterval(viewId1, viewId2); + if (overlap != null) + { + pairsToProcess.add(new ViewId[] { viewId1, viewId2 }); + } + } + } + + IOFunctions.println(" Found " + pairsToProcess.size() + " overlapping pairs"); } - // excess debug -// // Print global overlap region -// StringBuilder sb = new StringBuilder("("); -// for (int d = 0; d < globalOverlap.numDimensions(); d++) -// { -// if (d > 0) sb.append(", "); -// sb.append(String.format("%.1f", globalOverlap.realMin(d))); -// sb.append(" - "); -// sb.append(String.format("%.1f", globalOverlap.realMax(d))); -// } -// sb.append(")"); -// IOFunctions.println(" Global overlap region: " + sb); - - // Step 3: Get bounding boxes and convert overlap to local coordinates - BasicViewSetup vs1 = spimData.getSequenceDescription().getViewSetups().get(viewId1.getViewSetupId()); - BasicViewSetup vs2 = spimData.getSequenceDescription().getViewSetups().get(viewId2.getViewSetupId()); - ViewRegistration vr1 = spimData.getViewRegistrations().getViewRegistration(viewId1); - ViewRegistration vr2 = spimData.getViewRegistrations().getViewRegistration(viewId2); - - vr1.updateModel(); - vr2.updateModel(); - - RealInterval bbox1 = SimpleBoundingBoxOverlap.getBoundingBoxReal(vs1, vr1); - RealInterval bbox2 = SimpleBoundingBoxOverlap.getBoundingBoxReal(vs2, vr2); - - // Convert global overlap to local coordinates - double[] localMin1 = new double[globalOverlap.numDimensions()]; - double[] localMax1 = new double[globalOverlap.numDimensions()]; - double[] localMin2 = new double[globalOverlap.numDimensions()]; - double[] localMax2 = new double[globalOverlap.numDimensions()]; - - for (int d = 0; d < globalOverlap.numDimensions(); d++) + // Process all pairs + int successCount = 0; + for (ViewId[] pair : pairsToProcess) + { + try + { + computeCorrelationForPair(spimData, pair[0], pair[1], downsampleFactors, downsamplingChoice); + successCount++; + } + catch (Exception ex) { - localMin1[d] = globalOverlap.realMin(d) - bbox1.realMin(d); - localMax1[d] = globalOverlap.realMax(d) - bbox1.realMin(d); - localMin2[d] = globalOverlap.realMin(d) - bbox2.realMin(d); - localMax2[d] = globalOverlap.realMax(d) - bbox2.realMin(d); + // Silently skip pairs with errors (e.g., no overlap after downsampling) + } } - // Convert to raster coordinates with bounds checking - long[] dims1 = new long[globalOverlap.numDimensions()]; - long[] dims2 = new long[globalOverlap.numDimensions()]; - vs1.getSize().dimensions(dims1); - vs2.getSize().dimensions(dims2); + IOFunctions.println(new Date(System.currentTimeMillis()) + + ": Completed " + successCount + " of " + pairsToProcess.size() + " correlations"); + } + catch (Exception ex) + { + IOFunctions.println(new Date(System.currentTimeMillis()) + + ": ERROR: Exception during correlation computation: " + ex.getMessage()); + ex.printStackTrace(); + } + } + }).start(); + } + } - long[] rasterMin1 = new long[globalOverlap.numDimensions()]; - long[] rasterMax1 = new long[globalOverlap.numDimensions()]; - long[] rasterMin2 = new long[globalOverlap.numDimensions()]; - long[] rasterMax2 = new long[globalOverlap.numDimensions()]; + /** + * Compute cross-correlation for a specific pair of views + */ + private static void computeCorrelationForPair( + SpimData2 spimData, + ViewId viewId1, + ViewId viewId2, + long[] downsampleFactors, + String downsamplingChoice) throws Exception + { + // Calculate overlap using SimpleBoundingBoxOverlap + SimpleBoundingBoxOverlap overlapDetection = new SimpleBoundingBoxOverlap<>(spimData); + RealInterval globalOverlap = overlapDetection.getOverlapInterval(viewId1, viewId2); - for (int d = 0; d < globalOverlap.numDimensions(); d++) - { - rasterMin1[d] = Math.max(0, (long) Math.ceil(localMin1[d])); - rasterMax1[d] = Math.min(dims1[d] - 1, (long) Math.floor(localMax1[d])); - rasterMin2[d] = Math.max(0, (long) Math.ceil(localMin2[d])); - rasterMax2[d] = Math.min(dims2[d] - 1, (long) Math.floor(localMax2[d])); - } + if (globalOverlap == null) + { + throw new Exception("No overlap found"); + } - // Step 4: Ask user for downsampling level - // Get available downsampling levels from the data pyramid - String[] availableDownsamplings = DownsampleTools.availableDownsamplings(spimData, viewId1); + // Get bounding boxes and convert overlap to local coordinates + BasicViewSetup vs1 = spimData.getSequenceDescription().getViewSetups().get(viewId1.getViewSetupId()); + BasicViewSetup vs2 = spimData.getSequenceDescription().getViewSetups().get(viewId2.getViewSetupId()); + ViewRegistration vr1 = spimData.getViewRegistrations().getViewRegistration(viewId1); + ViewRegistration vr2 = spimData.getViewRegistrations().getViewRegistration(viewId2); - final GenericDialog gd = new GenericDialog("Cross-Correlation Parameters"); - gd.addChoice("Downsampling", availableDownsamplings, availableDownsamplings[0]); - gd.addMessage("Choose downsampling level:"); - gd.addMessage("1, 1, 1 = full resolution (slower, more accurate)"); - gd.addMessage("Higher values = faster computation (less accurate)"); - gd.showDialog(); + vr1.updateModel(); + vr2.updateModel(); - if (gd.wasCanceled()) - { - IOFunctions.println(new Date(System.currentTimeMillis()) + - ": Cross-correlation computation cancelled by user"); - return; - } + RealInterval bbox1 = SimpleBoundingBoxOverlap.getBoundingBoxReal(vs1, vr1); + RealInterval bbox2 = SimpleBoundingBoxOverlap.getBoundingBoxReal(vs2, vr2); - final String downsamplingChoice = gd.getNextChoice(); - final long[] downsampleFactors = DownsampleTools.parseDownsampleChoice(downsamplingChoice); + // Convert global overlap to local coordinates + double[] localMin1 = new double[globalOverlap.numDimensions()]; + double[] localMax1 = new double[globalOverlap.numDimensions()]; + double[] localMin2 = new double[globalOverlap.numDimensions()]; + double[] localMax2 = new double[globalOverlap.numDimensions()]; - IOFunctions.println(" Using downsampling: " + downsamplingChoice); + for (int d = 0; d < globalOverlap.numDimensions(); d++) + { + localMin1[d] = globalOverlap.realMin(d) - bbox1.realMin(d); + localMax1[d] = globalOverlap.realMax(d) - bbox1.realMin(d); + localMin2[d] = globalOverlap.realMin(d) - bbox2.realMin(d); + localMax2[d] = globalOverlap.realMax(d) - bbox2.realMin(d); + } - // Adjust overlap intervals for downsampling - long[] dsRasterMin1 = new long[rasterMin1.length]; - long[] dsRasterMax1 = new long[rasterMax1.length]; - long[] dsRasterMin2 = new long[rasterMin2.length]; - long[] dsRasterMax2 = new long[rasterMax2.length]; + // Convert to raster coordinates with bounds checking + long[] dims1 = new long[globalOverlap.numDimensions()]; + long[] dims2 = new long[globalOverlap.numDimensions()]; + vs1.getSize().dimensions(dims1); + vs2.getSize().dimensions(dims2); - for (int d = 0; d < rasterMin1.length; d++) - { - dsRasterMin1[d] = rasterMin1[d] / downsampleFactors[d]; - dsRasterMax1[d] = rasterMax1[d] / downsampleFactors[d]; - dsRasterMin2[d] = rasterMin2[d] / downsampleFactors[d]; - dsRasterMax2[d] = rasterMax2[d] / downsampleFactors[d]; - } + long[] rasterMin1 = new long[globalOverlap.numDimensions()]; + long[] rasterMax1 = new long[globalOverlap.numDimensions()]; + long[] rasterMin2 = new long[globalOverlap.numDimensions()]; + long[] rasterMax2 = new long[globalOverlap.numDimensions()]; - Interval dsInterval1 = new FinalInterval(dsRasterMin1, dsRasterMax1); - Interval dsInterval2 = new FinalInterval(dsRasterMin2, dsRasterMax2); + for (int d = 0; d < globalOverlap.numDimensions(); d++) + { + rasterMin1[d] = Math.max(0, (long) Math.ceil(localMin1[d])); + rasterMax1[d] = Math.min(dims1[d] - 1, (long) Math.floor(localMax1[d])); + rasterMin2[d] = Math.max(0, (long) Math.ceil(localMin2[d])); + rasterMax2[d] = Math.min(dims2[d] - 1, (long) Math.floor(localMax2[d])); + } - // Step 5: Load image data at specified downsampling level using pyramid - final BasicImgLoader imgLoader = spimData.getSequenceDescription().getImgLoader(); + // Adjust overlap intervals for downsampling + long[] dsRasterMin1 = new long[rasterMin1.length]; + long[] dsRasterMax1 = new long[rasterMax1.length]; + long[] dsRasterMin2 = new long[rasterMin2.length]; + long[] dsRasterMax2 = new long[rasterMax2.length]; - // Open images at the specified downsampling level using pre-computed pyramid - net.imglib2.util.Pair opened1 = - DownsampleTools.openAndDownsample(imgLoader, viewId1, downsampleFactors, false); - net.imglib2.util.Pair opened2 = - DownsampleTools.openAndDownsample(imgLoader, viewId2, downsampleFactors, false); + for (int d = 0; d < rasterMin1.length; d++) + { + dsRasterMin1[d] = rasterMin1[d] / downsampleFactors[d]; + dsRasterMax1[d] = rasterMax1[d] / downsampleFactors[d]; + dsRasterMin2[d] = rasterMin2[d] / downsampleFactors[d]; + dsRasterMax2[d] = rasterMax2[d] / downsampleFactors[d]; + } - RandomAccessibleInterval img1 = opened1.getA(); - RandomAccessibleInterval img2 = opened2.getA(); + Interval dsInterval1 = new FinalInterval(dsRasterMin1, dsRasterMax1); + Interval dsInterval2 = new FinalInterval(dsRasterMin2, dsRasterMax2); - // Extract overlap regions and convert to FloatType - RandomAccessibleInterval overlap1 = Views.zeroMin( - Views.interval(Views.zeroMin(img1), dsInterval1)); - RandomAccessibleInterval overlap2 = Views.zeroMin( - Views.interval(Views.zeroMin(img2), dsInterval2)); - // Step 6: Compute cross-correlation directly - IOFunctions.println(" Computing cross-correlation..."); + // Validate that intervals are non-empty (have positive size in all dimensions) + for (int d = 0; d < dsInterval1.numDimensions(); d++) + { + if (dsInterval1.min(d) > dsInterval1.max(d) || dsInterval2.min(d) > dsInterval2.max(d)) + { + throw new Exception("Empty overlap interval after downsampling"); + } + } - // Create a peak with zero shift to compute correlation at current alignment - Point zeroShift = new Point(overlap1.numDimensions()); - PhaseCorrelationPeak2 peak = new PhaseCorrelationPeak2(zeroShift, 0.0); + // Calculate number of overlapping pixels to validate + long numPixels = 1; + for (int d = 0; d < dsInterval1.numDimensions(); d++) + { + numPixels *= (dsInterval1.max(d) - dsInterval1.min(d) + 1); + } - // Set the shift field (required by calculateCrossCorr) - peak.setShift(zeroShift); + if (numPixels == 0) + { + throw new Exception("No overlapping pixels after downsampling"); + } - // Calculate cross-correlation at zero shift (current alignment) - peak.calculateCrossCorr(overlap1, overlap2); + // Load image data at specified downsampling level using pyramid + final BasicImgLoader imgLoader = spimData.getSequenceDescription().getImgLoader(); - // Step 7: Log results - double correlationCoefficient = peak.getCrossCorr(); - long nPixels = peak.getnPixel(); + // Open images at the specified downsampling level using pre-computed pyramid + net.imglib2.util.Pair opened1 = + DownsampleTools.openAndDownsample(imgLoader, viewId1, downsampleFactors, false); + net.imglib2.util.Pair opened2 = + DownsampleTools.openAndDownsample(imgLoader, viewId2, downsampleFactors, false); - IOFunctions.println(new Date(System.currentTimeMillis()) + - String.format(" %d-%d <> %d-%d", viewId1.getTimePointId(), viewId1.getViewSetupId(), - viewId2.getTimePointId(), viewId2.getViewSetupId()) - + String.format(": r=%.4f", correlationCoefficient)); + RandomAccessibleInterval img1 = opened1.getA(); + RandomAccessibleInterval img2 = opened2.getA(); - } - catch (Exception ex) - { - IOFunctions.println(new Date(System.currentTimeMillis()) + - ": ERROR: Exception during correlation computation: " + ex.getMessage()); - ex.printStackTrace(); - } - } - }).start(); + // Extract overlap regions and convert to FloatType + RandomAccessibleInterval overlap1 = Views.zeroMin( + Views.interval( + Views.zeroMin((RandomAccessibleInterval) img1), + dsInterval1)); + RandomAccessibleInterval overlap2 = Views.zeroMin( + Views.interval( + Views.zeroMin((RandomAccessibleInterval) img2), + dsInterval2)); + + // Compute cross-correlation directly + // Create a peak with zero shift to compute correlation at current alignment + Point zeroShift = new Point(overlap1.numDimensions()); + PhaseCorrelationPeak2 peak = new PhaseCorrelationPeak2(zeroShift, 0.0); + + // Set the shift field (required by calculateCrossCorr) + peak.setShift(zeroShift); + + // Calculate cross-correlation at zero shift (current alignment) + peak.calculateCrossCorr(overlap1, overlap2); + + // Log results + double correlationCoefficient = peak.getCrossCorr(); + long nPixels = peak.getnPixel(); + + + // Validate results before logging + if (nPixels == 0) + { + throw new Exception("No overlapping pixels found"); } + + if (!Double.isFinite(correlationCoefficient)) + { + throw new Exception("Invalid correlation coefficient: " + correlationCoefficient); + } + + IOFunctions.println(new Date(System.currentTimeMillis()) + + String.format(" %d-%d <> %d-%d", + viewId1.getTimePointId(), viewId1.getViewSetupId(), + viewId2.getTimePointId(), viewId2.getViewSetupId()) + + String.format(": r=%.4f (n=%d)", correlationCoefficient, nPixels)); } } From 2f0ec68bba62155baf7290b8ecc3cc81b665f966 Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Fri, 9 Jan 2026 11:48:18 -0500 Subject: [PATCH 08/16] parallelized tasks --- .../popup/ComputeCrossCorrelationPopup.java | 54 +++++++++++++++---- 1 file changed, 45 insertions(+), 9 deletions(-) diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java index 9f08bbf2c..35426ce9c 100644 --- a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java @@ -27,6 +27,10 @@ import java.util.ArrayList; import java.util.Date; import java.util.List; +import java.util.concurrent.ExecutorService; +import java.util.concurrent.Executors; +import java.util.concurrent.Future; +import java.util.concurrent.atomic.AtomicInteger; import javax.swing.JComponent; import javax.swing.JMenuItem; @@ -221,23 +225,55 @@ public void run() IOFunctions.println(" Found " + pairsToProcess.size() + " overlapping pairs"); } - // Process all pairs - int successCount = 0; - for (ViewId[] pair : pairsToProcess) + // Process all pairs in parallel + final int numThreads = Math.min(pairsToProcess.size(), Runtime.getRuntime().availableProcessors()); + final ExecutorService executor = Executors.newFixedThreadPool(numThreads); + final AtomicInteger successCount = new AtomicInteger(0); + + IOFunctions.println(" Using " + numThreads + " threads for parallel processing"); + + List> futures = new ArrayList<>(); + + for (final ViewId[] pair : pairsToProcess) + { + Future future = executor.submit(new Runnable() + { + @Override + public void run() + { + try + { + computeCorrelationForPair(spimData, pair[0], pair[1], downsampleFactors, downsamplingChoice); + successCount.incrementAndGet(); + } + catch (Exception ex) + { + // Silently skip pairs with errors (e.g., no overlap after downsampling) + } + } + }); + futures.add(future); + } + + // Wait for all tasks to complete + for (Future future : futures) { try { - computeCorrelationForPair(spimData, pair[0], pair[1], downsampleFactors, downsamplingChoice); - successCount++; + future.get(); } catch (Exception ex) - { - // Silently skip pairs with errors (e.g., no overlap after downsampling) - } + { + // Already handled in the task + } } + // Shut down executor + executor.shutdown(); + + IOFunctions.println(new Date(System.currentTimeMillis()) + - ": Completed " + successCount + " of " + pairsToProcess.size() + " correlations"); + ": Completed " + successCount.get() + " of " + pairsToProcess.size() + " correlations"); } catch (Exception ex) { From 421bb8318a5b7f78dc56ee6c0235536e79ea61d2 Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Fri, 9 Jan 2026 12:25:02 -0500 Subject: [PATCH 09/16] added sampled avg intensity calculation --- .../popup/ComputeCrossCorrelationPopup.java | 105 +++++++++++++++--- 1 file changed, 90 insertions(+), 15 deletions(-) diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java index 35426ce9c..fde511da6 100644 --- a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java @@ -47,9 +47,10 @@ import net.imglib2.RandomAccessibleInterval; import net.imglib2.RealInterval; import net.imglib2.realtransform.AffineTransform3D; +import net.imglib2.type.numeric.RealType; import net.imglib2.type.numeric.real.FloatType; -import net.imglib2.util.Util; import net.imglib2.view.Views; +import net.imglib2.converter.Converters; import net.preibisch.legacy.io.IOFunctions; import net.preibisch.mvrecon.fiji.spimdata.SpimData2; import net.preibisch.mvrecon.fiji.spimdata.explorer.ExplorerWindow; @@ -180,7 +181,11 @@ public void run() ViewId viewId1 = selectedViewIds.get(i); ViewId viewId2 = selectedViewIds.get(j); - RealInterval overlap = overlapDetection.getOverlapInterval(viewId1, viewId2); + // Skip comparisons between different timepoints + if (viewId1.getTimePointId() != viewId2.getTimePointId()) + continue; + + RealInterval overlap = overlapDetection.getOverlapInterval(viewId1, viewId2); if (overlap != null) { pairsToProcess.add(new ViewId[] { viewId1, viewId2 }); @@ -214,11 +219,15 @@ public void run() ViewId viewId1 = viewIds.get(i); ViewId viewId2 = viewIds.get(j); - RealInterval overlap = overlapDetection.getOverlapInterval(viewId1, viewId2); - if (overlap != null) - { - pairsToProcess.add(new ViewId[] { viewId1, viewId2 }); - } + // Skip comparisons between different timepoints + if (viewId1.getTimePointId() != viewId2.getTimePointId()) + continue; + + RealInterval overlap = overlapDetection.getOverlapInterval(viewId1, viewId2); + if (overlap != null) + { + pairsToProcess.add(new ViewId[] { viewId1, viewId2 }); + } } } @@ -248,7 +257,13 @@ public void run() } catch (Exception ex) { - // Silently skip pairs with errors (e.g., no overlap after downsampling) +// // Temporarily log first error to diagnose issue +// if (successCount.get() == 0) +// { +// IOFunctions.println(new Date(System.currentTimeMillis()) + +// ": ERROR (first failure): " + ex.getMessage()); +// ex.printStackTrace(); +// } } } }); @@ -398,19 +413,35 @@ private static void computeCorrelationForPair( net.imglib2.util.Pair opened2 = DownsampleTools.openAndDownsample(imgLoader, viewId2, downsampleFactors, false); - RandomAccessibleInterval img1 = opened1.getA(); - RandomAccessibleInterval img2 = opened2.getA(); - - // Extract overlap regions and convert to FloatType + RandomAccessibleInterval img1Raw = opened1.getA(); + RandomAccessibleInterval img2Raw = opened2.getA(); + + // Convert to FloatType (handles any numeric type) + @SuppressWarnings("unchecked") + RandomAccessibleInterval img1 = Converters.convert( + (RandomAccessibleInterval>) img1Raw, + (in, out) -> out.setReal(in.getRealDouble()), + new FloatType()); + @SuppressWarnings("unchecked") + RandomAccessibleInterval img2 = Converters.convert( + (RandomAccessibleInterval>) img2Raw, + (in, out) -> out.setReal(in.getRealDouble()), + new FloatType()); + + // Extract overlap regions RandomAccessibleInterval overlap1 = Views.zeroMin( Views.interval( - Views.zeroMin((RandomAccessibleInterval) img1), + Views.zeroMin(img1), dsInterval1)); RandomAccessibleInterval overlap2 = Views.zeroMin( Views.interval( - Views.zeroMin((RandomAccessibleInterval) img2), + Views.zeroMin(img2), dsInterval2)); + // Compute average intensities in the overlapping regions (with sampling for speed) + double avgIntensity1 = computeAverageIntensity(overlap1, 10); + double avgIntensity2 = computeAverageIntensity(overlap2, 10); + // Compute cross-correlation directly // Create a peak with zero shift to compute correlation at current alignment Point zeroShift = new Point(overlap1.numDimensions()); @@ -442,6 +473,50 @@ private static void computeCorrelationForPair( String.format(" %d-%d <> %d-%d", viewId1.getTimePointId(), viewId1.getViewSetupId(), viewId2.getTimePointId(), viewId2.getViewSetupId()) + - String.format(": r=%.4f (n=%d)", correlationCoefficient, nPixels)); + String.format(": r=%.4f (n=%d) avg int=[%.1f, %.1f]", + correlationCoefficient, nPixels, avgIntensity1, avgIntensity2)); + } + + /** + * Compute average intensity of an image using sampling for efficiency. + * Samples every stepSize-th pixel in each dimension. + * + * @param img The image to sample + * @param stepSize Sampling step (e.g., 10 means sample every 10th pixel) + * @return Average intensity of sampled pixels + */ + private static double computeAverageIntensity(RandomAccessibleInterval img, int stepSize) + { + double sum = 0; + long count = 0; + + // Create a cursor that samples at regular intervals + final net.imglib2.Cursor cursor = Views.iterable(img).cursor(); + final long[] position = new long[img.numDimensions()]; + + while (cursor.hasNext()) + { + cursor.fwd(); + cursor.localize(position); + + // Sample only at regular intervals + boolean shouldSample = true; + for (int d = 0; d < position.length; d++) + { + if (position[d] % stepSize != 0) + { + shouldSample = false; + break; + } + } + + if (shouldSample) + { + sum += cursor.get().get(); + count++; + } + } + + return count > 0 ? sum / count : 0.0; } } From 789d7d21a4686dc18a471fbfd841b9dd63dcf6a2 Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Fri, 9 Jan 2026 14:13:38 -0500 Subject: [PATCH 10/16] added channel choice --- .../popup/ComputeCrossCorrelationPopup.java | 97 ++++++++++++++----- 1 file changed, 71 insertions(+), 26 deletions(-) diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java index fde511da6..affba7db2 100644 --- a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java @@ -31,6 +31,7 @@ import java.util.concurrent.Executors; import java.util.concurrent.Future; import java.util.concurrent.atomic.AtomicInteger; +import java.util.stream.Collectors; import javax.swing.JComponent; import javax.swing.JMenuItem; @@ -40,6 +41,7 @@ import mpicbg.spim.data.generic.sequence.BasicImgLoader; import mpicbg.spim.data.generic.sequence.BasicViewSetup; import mpicbg.spim.data.registration.ViewRegistration; +import mpicbg.spim.data.sequence.Channel; import mpicbg.spim.data.sequence.ViewId; import net.imglib2.FinalInterval; import net.imglib2.Interval; @@ -112,14 +114,22 @@ public void run() ViewId sampleViewId = allViewIds.get(0); String[] availableDownsamplings = DownsampleTools.availableDownsamplings(spimData, sampleViewId); + List allChannels = new ArrayList<>(spimData.getSequenceDescription().getAllChannels().keySet()); + List allChannelsStr = allChannels.stream() + .map(i -> String.format("Channel %d", i)) + .collect(Collectors.toList()); + allChannelsStr.add("Compare across all channels"); + String[] allChannelsArr = allChannelsStr.toArray(new String[0]); + // Show dialog with checkbox and downsampling choice final GenericDialog gd = new GenericDialog("Cross-Correlation Parameters"); gd.addCheckbox("Only compute for selected tiles (requires at least 2 selected)", false); - gd.addChoice("Downsampling", availableDownsamplings, availableDownsamplings[0]); gd.addMessage("Choose downsampling level:"); gd.addMessage("1, 1, 1 = full resolution (slower, more accurate)"); gd.addMessage("Higher values = faster computation (less accurate)"); + gd.addChoice("Downsampling", availableDownsamplings, availableDownsamplings[0]); + gd.addChoice("Channel", allChannelsArr, allChannelsArr[0]); gd.showDialog(); if (gd.wasCanceled()) @@ -131,11 +141,22 @@ public void run() final boolean selectedOnly = gd.getNextBoolean(); final String downsamplingChoice = gd.getNextChoice(); + final String chChoiceStr = gd.getNextChoice(); final long[] downsampleFactors = DownsampleTools.parseDownsampleChoice(downsamplingChoice); + final int chChoice; + if (chChoiceStr.equals("Compare across all channels")) { + chChoice = -1; + } + else { + String[] parts = chChoiceStr.split(" "); + chChoice = Integer.parseInt(parts[parts.length - 1]); + } + IOFunctions.println(new Date(System.currentTimeMillis()) + ": Computing cross-correlations..."); IOFunctions.println(" Using downsampling: " + downsamplingChoice); + IOFunctions.println(" Using channel: " + chChoiceStr); List pairsToProcess = new ArrayList<>(); @@ -185,6 +206,15 @@ public void run() if (viewId1.getTimePointId() != viewId2.getTimePointId()) continue; + // Skip cross-channel comparisons if specific channel selected + if (chChoice != -1) + { + int ch1 = getChannelId(spimData, viewId1); + int ch2 = getChannelId(spimData, viewId2); + if (ch1 != chChoice || ch2 != chChoice) + continue; + } + RealInterval overlap = overlapDetection.getOverlapInterval(viewId1, viewId2); if (overlap != null) { @@ -223,6 +253,15 @@ public void run() if (viewId1.getTimePointId() != viewId2.getTimePointId()) continue; + // Skip cross-channel comparisons if specific channel selected + if (chChoice != -1) + { + int ch1 = getChannelId(spimData, viewId1); + int ch2 = getChannelId(spimData, viewId2); + if (ch1 != chChoice || ch2 != chChoice) + continue; + } + RealInterval overlap = overlapDetection.getOverlapInterval(viewId1, viewId2); if (overlap != null) { @@ -252,7 +291,7 @@ public void run() { try { - computeCorrelationForPair(spimData, pair[0], pair[1], downsampleFactors, downsamplingChoice); + computeCorrelationForPair(spimData, pair[0], pair[1], downsampleFactors); successCount.incrementAndGet(); } catch (Exception ex) @@ -308,8 +347,8 @@ private static void computeCorrelationForPair( SpimData2 spimData, ViewId viewId1, ViewId viewId2, - long[] downsampleFactors, - String downsamplingChoice) throws Exception + long[] downsampleFactors + ) throws Exception { // Calculate overlap using SimpleBoundingBoxOverlap SimpleBoundingBoxOverlap overlapDetection = new SimpleBoundingBoxOverlap<>(spimData); @@ -382,27 +421,16 @@ private static void computeCorrelationForPair( Interval dsInterval1 = new FinalInterval(dsRasterMin1, dsRasterMax1); Interval dsInterval2 = new FinalInterval(dsRasterMin2, dsRasterMax2); + // check whether we have 0-sized (or negative sized) + // ignore this pair in that case + for ( int d = 0; d < dsInterval1.numDimensions(); ++d ) + { + if ( dsInterval1.dimension( d ) <= 0 || dsInterval2.dimension( d ) <= 0) + { + throw new Exception("Rastered overlap volume is zero, skipping." ); + } + } - // Validate that intervals are non-empty (have positive size in all dimensions) - for (int d = 0; d < dsInterval1.numDimensions(); d++) - { - if (dsInterval1.min(d) > dsInterval1.max(d) || dsInterval2.min(d) > dsInterval2.max(d)) - { - throw new Exception("Empty overlap interval after downsampling"); - } - } - - // Calculate number of overlapping pixels to validate - long numPixels = 1; - for (int d = 0; d < dsInterval1.numDimensions(); d++) - { - numPixels *= (dsInterval1.max(d) - dsInterval1.min(d) + 1); - } - - if (numPixels == 0) - { - throw new Exception("No overlapping pixels after downsampling"); - } // Load image data at specified downsampling level using pyramid final BasicImgLoader imgLoader = spimData.getSequenceDescription().getImgLoader(); @@ -473,8 +501,25 @@ private static void computeCorrelationForPair( String.format(" %d-%d <> %d-%d", viewId1.getTimePointId(), viewId1.getViewSetupId(), viewId2.getTimePointId(), viewId2.getViewSetupId()) + - String.format(": r=%.4f (n=%d) avg int=[%.1f, %.1f]", - correlationCoefficient, nPixels, avgIntensity1, avgIntensity2)); + String.format(": r=%.4f avg int=[%.1f, %.1f]", + correlationCoefficient, avgIntensity1, avgIntensity2)); + } + + /** + * Get the channel ID for a given ViewId. + * + * @param spimData The SpimData2 object + * @param viewId The ViewId to query + * @return The channel ID, or -1 if channel information is not available + */ + private static int getChannelId(SpimData2 spimData, ViewId viewId) + { + mpicbg.spim.data.sequence.ViewDescription vd = spimData.getSequenceDescription().getViewDescriptions().get(viewId); + if (vd != null && vd.getViewSetup().getChannel() != null) + { + return vd.getViewSetup().getChannel().getId(); + } + return -1; } /** From 54f32b2763670ff22aa6a9f9e5dc72ce07e7555b Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Fri, 9 Jan 2026 14:21:02 -0500 Subject: [PATCH 11/16] added summary class at end --- .../popup/ComputeCrossCorrelationPopup.java | 73 ++++++++++++++++++- 1 file changed, 69 insertions(+), 4 deletions(-) diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java index affba7db2..2c417807c 100644 --- a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java @@ -25,6 +25,8 @@ import java.awt.event.ActionEvent; import java.awt.event.ActionListener; import java.util.ArrayList; +import java.util.Collections; +import java.util.Comparator; import java.util.Date; import java.util.List; import java.util.concurrent.ExecutorService; @@ -66,6 +68,27 @@ public class ComputeCrossCorrelationPopup extends JMenuItem implements ExplorerW private static final long serialVersionUID = 1L; ExplorerWindow panel; + /** + * Class to store correlation result for a pair of views + */ + private static class CorrelationResult + { + final ViewId viewId1; + final ViewId viewId2; + final double correlation; + final double avgIntensity1; + final double avgIntensity2; + + CorrelationResult(ViewId viewId1, ViewId viewId2, double correlation, double avgIntensity1, double avgIntensity2) + { + this.viewId1 = viewId1; + this.viewId2 = viewId2; + this.correlation = correlation; + this.avgIntensity1 = avgIntensity1; + this.avgIntensity2 = avgIntensity2; + } + } + public ComputeCrossCorrelationPopup() { super("Compute Cross-Correlation for Overlapping Tiles"); @@ -277,6 +300,7 @@ public void run() final int numThreads = Math.min(pairsToProcess.size(), Runtime.getRuntime().availableProcessors()); final ExecutorService executor = Executors.newFixedThreadPool(numThreads); final AtomicInteger successCount = new AtomicInteger(0); + final List results = Collections.synchronizedList(new ArrayList<>()); IOFunctions.println(" Using " + numThreads + " threads for parallel processing"); @@ -291,8 +315,12 @@ public void run() { try { - computeCorrelationForPair(spimData, pair[0], pair[1], downsampleFactors); - successCount.incrementAndGet(); + CorrelationResult result = computeCorrelationForPair(spimData, pair[0], pair[1], downsampleFactors); + if (result != null) + { + results.add(result); + successCount.incrementAndGet(); + } } catch (Exception ex) { @@ -325,9 +353,42 @@ public void run() // Shut down executor executor.shutdown(); - IOFunctions.println(new Date(System.currentTimeMillis()) + ": Completed " + successCount.get() + " of " + pairsToProcess.size() + " correlations"); + + // Print summary table sorted by correlation (ascending) + if (!results.isEmpty()) + { + IOFunctions.println("\n" + new Date(System.currentTimeMillis()) + + ": Summary of correlations (sorted by correlation):"); + IOFunctions.println("================================================================================"); + IOFunctions.println(String.format("%-20s %-20s %12s %12s %12s", + "ViewId 1", "ViewId 2", "Correlation", "Avg Int 1", "Avg Int 2")); + IOFunctions.println("================================================================================"); + + // Sort by correlation (ascending) + Collections.sort(results, new Comparator() + { + @Override + public int compare(CorrelationResult r1, CorrelationResult r2) + { + return Double.compare(r1.correlation, r2.correlation); + } + }); + + // Print each result + for (CorrelationResult result : results) + { + String viewId1Str = String.format("%d-%d", + result.viewId1.getTimePointId(), result.viewId1.getViewSetupId()); + String viewId2Str = String.format("%d-%d", + result.viewId2.getTimePointId(), result.viewId2.getViewSetupId()); + IOFunctions.println(String.format("%-20s %-20s %12.4f %12.1f %12.1f", + viewId1Str, viewId2Str, result.correlation, + result.avgIntensity1, result.avgIntensity2)); + } + IOFunctions.println("================================================================================"); + } } catch (Exception ex) { @@ -343,7 +404,7 @@ public void run() /** * Compute cross-correlation for a specific pair of views */ - private static void computeCorrelationForPair( + private static CorrelationResult computeCorrelationForPair( SpimData2 spimData, ViewId viewId1, ViewId viewId2, @@ -497,12 +558,16 @@ private static void computeCorrelationForPair( throw new Exception("Invalid correlation coefficient: " + correlationCoefficient); } + // Log individual result IOFunctions.println(new Date(System.currentTimeMillis()) + String.format(" %d-%d <> %d-%d", viewId1.getTimePointId(), viewId1.getViewSetupId(), viewId2.getTimePointId(), viewId2.getViewSetupId()) + String.format(": r=%.4f avg int=[%.1f, %.1f]", correlationCoefficient, avgIntensity1, avgIntensity2)); + + // Return result for summary table + return new CorrelationResult(viewId1, viewId2, correlationCoefficient, avgIntensity1, avgIntensity2); } /** From 891e92fc7c4492c252b06ee462e7164b0fdd1587 Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Mon, 12 Jan 2026 11:20:32 -0500 Subject: [PATCH 12/16] bug fix to solve 1st bbox transformation problem --- .../popup/ComputeCrossCorrelationPopup.java | 218 +++++++++++------- 1 file changed, 135 insertions(+), 83 deletions(-) diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java index 2c417807c..c49733a7f 100644 --- a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java @@ -181,7 +181,8 @@ public void run() IOFunctions.println(" Using downsampling: " + downsamplingChoice); IOFunctions.println(" Using channel: " + chChoiceStr); - List pairsToProcess = new ArrayList<>(); + // Get the list of views to process + List viewIdsToProcess; if (selectedOnly) { @@ -195,106 +196,47 @@ public void run() } // Filter missing views and collect valid ViewIds from selected groups - List selectedViewIds = new ArrayList<>(); + viewIdsToProcess = new ArrayList<>(); for (List group : selectedGroups) { SpimData2.filterMissingViews(spimData, group); if (!group.isEmpty()) { // Get first ViewId from each group - selectedViewIds.add(group.get(0)); + viewIdsToProcess.add(group.get(0)); } } - if (selectedViewIds.size() < 2) + if (viewIdsToProcess.size() < 2) { IOFunctions.println(new Date(System.currentTimeMillis()) + ": ERROR: Selected groups contain fewer than 2 valid views"); return; } - IOFunctions.println(" Processing " + selectedViewIds.size() + " selected views"); - - SimpleBoundingBoxOverlap overlapDetection = new SimpleBoundingBoxOverlap<>(spimData); - - // Find all overlapping pairs among selected views - for (int i = 0; i < selectedViewIds.size(); i++) - { - for (int j = i + 1; j < selectedViewIds.size(); j++) - { - ViewId viewId1 = selectedViewIds.get(i); - ViewId viewId2 = selectedViewIds.get(j); - - // Skip comparisons between different timepoints - if (viewId1.getTimePointId() != viewId2.getTimePointId()) - continue; - - // Skip cross-channel comparisons if specific channel selected - if (chChoice != -1) - { - int ch1 = getChannelId(spimData, viewId1); - int ch2 = getChannelId(spimData, viewId2); - if (ch1 != chChoice || ch2 != chChoice) - continue; - } - - RealInterval overlap = overlapDetection.getOverlapInterval(viewId1, viewId2); - if (overlap != null) - { - pairsToProcess.add(new ViewId[] { viewId1, viewId2 }); - } - } - } - - IOFunctions.println(" Found " + pairsToProcess.size() + " overlapping pairs among selected views"); + IOFunctions.println(" Processing " + viewIdsToProcess.size() + " selected views"); } else { - // Get all views and find overlapping pairs - List viewIds = new ArrayList<>(); + // Get all views + viewIdsToProcess = new ArrayList<>(); for (ViewId viewId : allViewIds) { if (!spimData.getSequenceDescription().getMissingViews().getMissingViews().contains(viewId)) { - viewIds.add(viewId); + viewIdsToProcess.add(viewId); } } - IOFunctions.println(" Found " + viewIds.size() + " valid views"); - - SimpleBoundingBoxOverlap overlapDetection = new SimpleBoundingBoxOverlap<>(spimData); - - // Find all overlapping pairs - for (int i = 0; i < viewIds.size(); i++) - { - for (int j = i + 1; j < viewIds.size(); j++) - { - ViewId viewId1 = viewIds.get(i); - ViewId viewId2 = viewIds.get(j); - - // Skip comparisons between different timepoints - if (viewId1.getTimePointId() != viewId2.getTimePointId()) - continue; - - // Skip cross-channel comparisons if specific channel selected - if (chChoice != -1) - { - int ch1 = getChannelId(spimData, viewId1); - int ch2 = getChannelId(spimData, viewId2); - if (ch1 != chChoice || ch2 != chChoice) - continue; - } + IOFunctions.println(" Found " + viewIdsToProcess.size() + " valid views"); + } - RealInterval overlap = overlapDetection.getOverlapInterval(viewId1, viewId2); - if (overlap != null) - { - pairsToProcess.add(new ViewId[] { viewId1, viewId2 }); - } - } - } + // Find all overlapping pairs + List pairsToProcess = findOverlappingPairs( + spimData, viewIdsToProcess, chChoice); - IOFunctions.println(" Found " + pairsToProcess.size() + " overlapping pairs"); - } + IOFunctions.println(" Found " + pairsToProcess.size() + " overlapping pairs" + + (selectedOnly ? " among selected views" : "")); // Process all pairs in parallel final int numThreads = Math.min(pairsToProcess.size(), Runtime.getRuntime().availableProcessors()); @@ -432,7 +374,34 @@ private static CorrelationResult computeCorrelationForPair( RealInterval bbox1 = SimpleBoundingBoxOverlap.getBoundingBoxReal(vs1, vr1); RealInterval bbox2 = SimpleBoundingBoxOverlap.getBoundingBoxReal(vs2, vr2); - // Convert global overlap to local coordinates + // Get image dimensions + long[] dims1 = new long[globalOverlap.numDimensions()]; + long[] dims2 = new long[globalOverlap.numDimensions()]; + vs1.getSize().dimensions(dims1); + vs2.getSize().dimensions(dims2); + + // Debug: Print bounding boxes and global overlap + IOFunctions.println(String.format("DEBUG: Analyzing %d-%d <> %d-%d", + viewId1.getTimePointId(), viewId1.getViewSetupId(), + viewId2.getTimePointId(), viewId2.getViewSetupId())); + IOFunctions.println(String.format(" Image sizes: view1=[%d,%d,%d] view2=[%d,%d,%d]", + dims1[0], dims1[1], dims1[2], dims2[0], dims2[1], dims2[2])); + for (int d = 0; d < globalOverlap.numDimensions(); d++) + { + IOFunctions.println(String.format(" Dim %d: bbox1=[%.1f,%.1f] bbox2=[%.1f,%.1f] overlap=[%.1f,%.1f]", + d, bbox1.realMin(d), bbox1.realMax(d), + bbox2.realMin(d), bbox2.realMax(d), + globalOverlap.realMin(d), globalOverlap.realMax(d))); + } + + // Convert global overlap to local coordinates using inverse transforms + AffineTransform3D transform1 = vr1.getModel(); + AffineTransform3D transform2 = vr2.getModel(); + AffineTransform3D invTransform1 = transform1.inverse(); + AffineTransform3D invTransform2 = transform2.inverse(); + + double[] globalMin = new double[globalOverlap.numDimensions()]; + double[] globalMax = new double[globalOverlap.numDimensions()]; double[] localMin1 = new double[globalOverlap.numDimensions()]; double[] localMax1 = new double[globalOverlap.numDimensions()]; double[] localMin2 = new double[globalOverlap.numDimensions()]; @@ -440,17 +409,41 @@ private static CorrelationResult computeCorrelationForPair( for (int d = 0; d < globalOverlap.numDimensions(); d++) { - localMin1[d] = globalOverlap.realMin(d) - bbox1.realMin(d); - localMax1[d] = globalOverlap.realMax(d) - bbox1.realMin(d); - localMin2[d] = globalOverlap.realMin(d) - bbox2.realMin(d); - localMax2[d] = globalOverlap.realMax(d) - bbox2.realMin(d); + globalMin[d] = globalOverlap.realMin(d); + globalMax[d] = globalOverlap.realMax(d); + } + + // Apply inverse transforms to convert world coords to local pixel coords + invTransform1.apply(globalMin, localMin1); + invTransform1.apply(globalMax, localMax1); + invTransform2.apply(globalMin, localMin2); + invTransform2.apply(globalMax, localMax2); + + // Ensure min < max (inverse transform might swap them) + for (int d = 0; d < globalOverlap.numDimensions(); d++) + { + if (localMin1[d] > localMax1[d]) + { + double tmp = localMin1[d]; + localMin1[d] = localMax1[d]; + localMax1[d] = tmp; + } + if (localMin2[d] > localMax2[d]) + { + double tmp = localMin2[d]; + localMin2[d] = localMax2[d]; + localMax2[d] = tmp; + } + } + + // Debug: Print local coordinates + for (int d = 0; d < globalOverlap.numDimensions(); d++) + { + IOFunctions.println(String.format(" Dim %d local: view1=[%.1f,%.1f] view2=[%.1f,%.1f]", + d, localMin1[d], localMax1[d], localMin2[d], localMax2[d])); } // Convert to raster coordinates with bounds checking - long[] dims1 = new long[globalOverlap.numDimensions()]; - long[] dims2 = new long[globalOverlap.numDimensions()]; - vs1.getSize().dimensions(dims1); - vs2.getSize().dimensions(dims2); long[] rasterMin1 = new long[globalOverlap.numDimensions()]; long[] rasterMax1 = new long[globalOverlap.numDimensions()]; @@ -488,6 +481,18 @@ private static CorrelationResult computeCorrelationForPair( { if ( dsInterval1.dimension( d ) <= 0 || dsInterval2.dimension( d ) <= 0) { + // Debug logging + IOFunctions.println(String.format("DEBUG: Zero overlap for %d-%d <> %d-%d", + viewId1.getTimePointId(), viewId1.getViewSetupId(), + viewId2.getTimePointId(), viewId2.getViewSetupId())); + IOFunctions.println(String.format(" Dimension %d: view1=[%d,%d] (%d px), view2=[%d,%d] (%d px)", + d, dsRasterMin1[d], dsRasterMax1[d], dsInterval1.dimension(d), + dsRasterMin2[d], dsRasterMax2[d], dsInterval2.dimension(d))); + IOFunctions.println(String.format(" Original raster: view1=[%d,%d], view2=[%d,%d]", + rasterMin1[d], rasterMax1[d], rasterMin2[d], rasterMax2[d])); + IOFunctions.println(String.format(" Downsampling factor: %d", downsampleFactors[d])); + IOFunctions.println(String.format(" Global overlap: [%.1f, %.1f]", + globalOverlap.realMin(d), globalOverlap.realMax(d))); throw new Exception("Rastered overlap volume is zero, skipping." ); } } @@ -629,4 +634,51 @@ private static double computeAverageIntensity(RandomAccessibleInterval 0 ? sum / count : 0.0; } + + /** + * Find all overlapping pairs among a list of views, with optional channel filtering. + * + * @param spimData The SpimData2 object + * @param viewIds List of views to check for overlaps + * @param chChoice Channel to filter by, or -1 for all channels + * @return List of overlapping view pairs + */ + private static List findOverlappingPairs( + SpimData2 spimData, + List viewIds, + int chChoice) + { + List pairsToProcess = new ArrayList<>(); + SimpleBoundingBoxOverlap overlapDetection = new SimpleBoundingBoxOverlap<>(spimData); + + for (int i = 0; i < viewIds.size(); i++) + { + for (int j = i + 1; j < viewIds.size(); j++) + { + ViewId viewId1 = viewIds.get(i); + ViewId viewId2 = viewIds.get(j); + + // Skip comparisons between different timepoints + if (viewId1.getTimePointId() != viewId2.getTimePointId()) + continue; + + // Skip cross-channel comparisons if specific channel selected + if (chChoice != -1) + { + int ch1 = getChannelId(spimData, viewId1); + int ch2 = getChannelId(spimData, viewId2); + if (ch1 != chChoice || ch2 != chChoice) + continue; + } + + RealInterval overlap = overlapDetection.getOverlapInterval(viewId1, viewId2); + if (overlap != null) + { + pairsToProcess.add(new ViewId[] { viewId1, viewId2 }); + } + } + } + + return pairsToProcess; + } } From 843bfa7efc0b66f452f37f11adf9e103bb8495a7 Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Mon, 12 Jan 2026 18:50:00 -0500 Subject: [PATCH 13/16] changed logic to sampling --- .../popup/ComputeCrossCorrelationPopup.java | 132 ++-- .../overlap/SimpleBoundingBoxOverlap.java | 700 ++++++++++++++++++ 2 files changed, 751 insertions(+), 81 deletions(-) diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java index c49733a7f..e3388bbf5 100644 --- a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java @@ -353,109 +353,81 @@ private static CorrelationResult computeCorrelationForPair( long[] downsampleFactors ) throws Exception { - // Calculate overlap using SimpleBoundingBoxOverlap - SimpleBoundingBoxOverlap overlapDetection = new SimpleBoundingBoxOverlap<>(spimData); - RealInterval globalOverlap = overlapDetection.getOverlapInterval(viewId1, viewId2); - - if (globalOverlap == null) - { - throw new Exception("No overlap found"); - } - - // Get bounding boxes and convert overlap to local coordinates + // Get view information BasicViewSetup vs1 = spimData.getSequenceDescription().getViewSetups().get(viewId1.getViewSetupId()); BasicViewSetup vs2 = spimData.getSequenceDescription().getViewSetups().get(viewId2.getViewSetupId()); ViewRegistration vr1 = spimData.getViewRegistrations().getViewRegistration(viewId1); ViewRegistration vr2 = spimData.getViewRegistrations().getViewRegistration(viewId2); - vr1.updateModel(); - vr2.updateModel(); + // thread-safe update of the model and creation of a copy for processing + AffineTransform3D m1, m2; + + synchronized (vr1) + { + vr1.updateModel(); + m1 = vr1.getModel().copy(); + } + + synchronized (vr2) + { + vr2.updateModel(); + m2 = vr2.getModel().copy(); + } + + // Calculate local overlaps using pixel validation that handles rotations correctly + RealInterval[] localOverlaps = SimpleBoundingBoxOverlap.getLocalOverlapsUsingPixelValidation( + vs1.getSize(), vs2.getSize(), + m1, m2); + + if (localOverlaps == null) + { + throw new Exception("No overlap found"); + } - RealInterval bbox1 = SimpleBoundingBoxOverlap.getBoundingBoxReal(vs1, vr1); - RealInterval bbox2 = SimpleBoundingBoxOverlap.getBoundingBoxReal(vs2, vr2); + RealInterval localOverlap1 = localOverlaps[0]; + RealInterval localOverlap2 = localOverlaps[1]; // Get image dimensions - long[] dims1 = new long[globalOverlap.numDimensions()]; - long[] dims2 = new long[globalOverlap.numDimensions()]; + long[] dims1 = new long[3]; + long[] dims2 = new long[3]; vs1.getSize().dimensions(dims1); vs2.getSize().dimensions(dims2); - // Debug: Print bounding boxes and global overlap + // Debug: Print local overlaps IOFunctions.println(String.format("DEBUG: Analyzing %d-%d <> %d-%d", viewId1.getTimePointId(), viewId1.getViewSetupId(), viewId2.getTimePointId(), viewId2.getViewSetupId())); IOFunctions.println(String.format(" Image sizes: view1=[%d,%d,%d] view2=[%d,%d,%d]", dims1[0], dims1[1], dims1[2], dims2[0], dims2[1], dims2[2])); - for (int d = 0; d < globalOverlap.numDimensions(); d++) + for (int d = 0; d < 3; d++) { - IOFunctions.println(String.format(" Dim %d: bbox1=[%.1f,%.1f] bbox2=[%.1f,%.1f] overlap=[%.1f,%.1f]", - d, bbox1.realMin(d), bbox1.realMax(d), - bbox2.realMin(d), bbox2.realMax(d), - globalOverlap.realMin(d), globalOverlap.realMax(d))); + IOFunctions.println(String.format(" Dim %d local: view1=[%.1f,%.1f] (size=%.1f) view2=[%.1f,%.1f] (size=%.1f)", + d, localOverlap1.realMin(d), localOverlap1.realMax(d), + localOverlap1.realMax(d) - localOverlap1.realMin(d), + localOverlap2.realMin(d), localOverlap2.realMax(d), + localOverlap2.realMax(d) - localOverlap2.realMin(d))); } - // Convert global overlap to local coordinates using inverse transforms - AffineTransform3D transform1 = vr1.getModel(); - AffineTransform3D transform2 = vr2.getModel(); - AffineTransform3D invTransform1 = transform1.inverse(); - AffineTransform3D invTransform2 = transform2.inverse(); - - double[] globalMin = new double[globalOverlap.numDimensions()]; - double[] globalMax = new double[globalOverlap.numDimensions()]; - double[] localMin1 = new double[globalOverlap.numDimensions()]; - double[] localMax1 = new double[globalOverlap.numDimensions()]; - double[] localMin2 = new double[globalOverlap.numDimensions()]; - double[] localMax2 = new double[globalOverlap.numDimensions()]; - - for (int d = 0; d < globalOverlap.numDimensions(); d++) - { - globalMin[d] = globalOverlap.realMin(d); - globalMax[d] = globalOverlap.realMax(d); - } - - // Apply inverse transforms to convert world coords to local pixel coords - invTransform1.apply(globalMin, localMin1); - invTransform1.apply(globalMax, localMax1); - invTransform2.apply(globalMin, localMin2); - invTransform2.apply(globalMax, localMax2); - - // Ensure min < max (inverse transform might swap them) - for (int d = 0; d < globalOverlap.numDimensions(); d++) - { - if (localMin1[d] > localMax1[d]) - { - double tmp = localMin1[d]; - localMin1[d] = localMax1[d]; - localMax1[d] = tmp; - } - if (localMin2[d] > localMax2[d]) - { - double tmp = localMin2[d]; - localMin2[d] = localMax2[d]; - localMax2[d] = tmp; - } - } + // Convert to raster coordinates with bounds checking + long[] rasterMin1 = new long[3]; + long[] rasterMax1 = new long[3]; + long[] rasterMin2 = new long[3]; + long[] rasterMax2 = new long[3]; - // Debug: Print local coordinates - for (int d = 0; d < globalOverlap.numDimensions(); d++) + for (int d = 0; d < 3; d++) { - IOFunctions.println(String.format(" Dim %d local: view1=[%.1f,%.1f] view2=[%.1f,%.1f]", - d, localMin1[d], localMax1[d], localMin2[d], localMax2[d])); + rasterMin1[d] = Math.max(0, (long) Math.ceil(localOverlap1.realMin(d))); + rasterMax1[d] = Math.min(dims1[d] - 1, (long) Math.floor(localOverlap1.realMax(d))); + rasterMin2[d] = Math.max(0, (long) Math.ceil(localOverlap2.realMin(d))); + rasterMax2[d] = Math.min(dims2[d] - 1, (long) Math.floor(localOverlap2.realMax(d))); } - // Convert to raster coordinates with bounds checking - - long[] rasterMin1 = new long[globalOverlap.numDimensions()]; - long[] rasterMax1 = new long[globalOverlap.numDimensions()]; - long[] rasterMin2 = new long[globalOverlap.numDimensions()]; - long[] rasterMax2 = new long[globalOverlap.numDimensions()]; - - for (int d = 0; d < globalOverlap.numDimensions(); d++) + // Debug: Print raster coordinates + for (int d = 0; d < 3; d++) { - rasterMin1[d] = Math.max(0, (long) Math.ceil(localMin1[d])); - rasterMax1[d] = Math.min(dims1[d] - 1, (long) Math.floor(localMax1[d])); - rasterMin2[d] = Math.max(0, (long) Math.ceil(localMin2[d])); - rasterMax2[d] = Math.min(dims2[d] - 1, (long) Math.floor(localMax2[d])); + IOFunctions.println(String.format(" Dim %d raster: view1=[%d,%d] (size=%d) view2=[%d,%d] (size=%d)", + d, rasterMin1[d], rasterMax1[d], rasterMax1[d] - rasterMin1[d] + 1, + rasterMin2[d], rasterMax2[d], rasterMax2[d] - rasterMin2[d] + 1)); } // Adjust overlap intervals for downsampling @@ -491,8 +463,6 @@ private static CorrelationResult computeCorrelationForPair( IOFunctions.println(String.format(" Original raster: view1=[%d,%d], view2=[%d,%d]", rasterMin1[d], rasterMax1[d], rasterMin2[d], rasterMax2[d])); IOFunctions.println(String.format(" Downsampling factor: %d", downsampleFactors[d])); - IOFunctions.println(String.format(" Global overlap: [%.1f, %.1f]", - globalOverlap.realMin(d), globalOverlap.realMax(d))); throw new Exception("Rastered overlap volume is zero, skipping." ); } } diff --git a/src/main/java/net/preibisch/mvrecon/process/interestpointregistration/pairwise/constellation/overlap/SimpleBoundingBoxOverlap.java b/src/main/java/net/preibisch/mvrecon/process/interestpointregistration/pairwise/constellation/overlap/SimpleBoundingBoxOverlap.java index a0b8995da..962ecc09e 100644 --- a/src/main/java/net/preibisch/mvrecon/process/interestpointregistration/pairwise/constellation/overlap/SimpleBoundingBoxOverlap.java +++ b/src/main/java/net/preibisch/mvrecon/process/interestpointregistration/pairwise/constellation/overlap/SimpleBoundingBoxOverlap.java @@ -32,7 +32,9 @@ import mpicbg.spim.data.sequence.SequenceDescription; import mpicbg.spim.data.sequence.ViewId; import net.imglib2.Dimensions; +import net.imglib2.FinalInterval; import net.imglib2.FinalRealInterval; +import net.imglib2.Interval; import net.imglib2.RealInterval; import net.imglib2.realtransform.AffineTransform3D; import net.preibisch.mvrecon.fiji.spimdata.boundingbox.BoundingBox; @@ -198,4 +200,702 @@ public static RealInterval getBoundingBoxReal( final Dimensions dims, final Affi return transform.estimateBounds( new FinalRealInterval( min, max ) ); } + + /** + * Compute overlap between two views by sampling points and validating they fall + * within both views. Returns two local intervals - one for each view. + * + * This correctly handles rotations by validating pixel-by-pixel. + * + * @param dims1 Dimensions of view 1 + * @param dims2 Dimensions of view 2 + * @param transform1 Local-to-global transform for view 1 + * @param transform2 Local-to-global transform for view 2 + * @return Array of [localOverlap1, localOverlap2], or null if no overlap + */ + public static RealInterval[] getLocalOverlapsUsingPixelValidation( + final Dimensions dims1, + final Dimensions dims2, + final AffineTransform3D transform1, + final AffineTransform3D transform2 ) + { + final AffineTransform3D invTransform1 = transform1.inverse(); + final AffineTransform3D invTransform2 = transform2.inverse(); + + // Sample in view 1's local space and validate pixels are in view 2 + final RealInterval localOverlap1 = sampleLocalSpaceAndValidate( + dims1, dims2, transform1, invTransform2 ); + + // Sample in view 2's local space and validate pixels are in view 1 + final RealInterval localOverlap2 = sampleLocalSpaceAndValidate( + dims2, dims1, transform2, invTransform1 ); + + if ( localOverlap1 == null || localOverlap2 == null ) + return null; + + return new RealInterval[] { localOverlap1, localOverlap2 }; + } + + /** + * Sample pixels in a view's local space and validate they map to valid pixels in the other view. + * Builds bounding box directly in local space from valid pixels. + * Includes refinement to find tighter boundaries. + * + * @param dimsLocal Dimensions of this view + * @param dimsOther Dimensions of other view + * @param transformLocal Local-to-global transform for this view + * @param invTransformOther Global-to-local transform for other view + * @return Bounding box of valid pixels in local space, or null if none found + */ + private static RealInterval sampleLocalSpaceAndValidate( + final Dimensions dimsLocal, + final Dimensions dimsOther, + final AffineTransform3D transformLocal, + final AffineTransform3D invTransformOther ) + { + final int n = dimsLocal.numDimensions(); + final double[] minValid = new double[ n ]; + final double[] maxValid = new double[ n ]; + + // Store the points that generated the min/max values to use as seeds for refinement + final double[][] minPoints = new double[ n ][ n ]; + final double[][] maxPoints = new double[ n ][ n ]; + + // Initialize with invalid values + for ( int d = 0; d < n; d++ ) + { + minValid[ d ] = Double.MAX_VALUE; + maxValid[ d ] = -Double.MAX_VALUE; + } + + final int stride = 10; + final double[] localPoint = new double[ n ]; + final double[] globalPoint = new double[ n ]; + final double[] otherLocalPoint = new double[ n ]; + + boolean foundValid = false; + + // Sample in this view's local space + for ( long z = 0; z < dimsLocal.dimension( 2 ); z += stride ) + { + for ( long y = 0; y < dimsLocal.dimension( 1 ); y += stride ) + { + for ( long x = 0; x < dimsLocal.dimension( 0 ); x += stride ) + { + localPoint[ 0 ] = x; + localPoint[ 1 ] = y; + localPoint[ 2 ] = z; + + if ( isValid( localPoint, dimsOther, transformLocal, invTransformOther, globalPoint, otherLocalPoint ) ) + { + foundValid = true; + for ( int d = 0; d < n; d++ ) + { + if ( localPoint[ d ] < minValid[ d ] ) + { + minValid[ d ] = localPoint[ d ]; + System.arraycopy( localPoint, 0, minPoints[ d ], 0, n ); + } + if ( localPoint[ d ] > maxValid[ d ] ) + { + maxValid[ d ] = localPoint[ d ]; + System.arraycopy( localPoint, 0, maxPoints[ d ], 0, n ); + } + } + } + } + } + } + + if ( !foundValid ) + return null; + + // Refine boundaries + // For each dimension, search outwards from the min/max points found + // First with step 1, then step 0.1 + refineBoundaries( minValid, maxValid, minPoints, maxPoints, dimsOther, transformLocal, invTransformOther, stride ); + + return new FinalRealInterval( minValid, maxValid ); + } + + private static void refineBoundaries( + final double[] minValid, + final double[] maxValid, + final double[][] minPoints, + final double[][] maxPoints, + final Dimensions dimsOther, + final AffineTransform3D transformLocal, + final AffineTransform3D invTransformOther, + final int initialStride ) + { + final int n = minValid.length; + final double[] globalPoint = new double[ n ]; + final double[] otherLocalPoint = new double[ n ]; + final double[] testPoint = new double[ n ]; + + // Refinement steps + final double[] steps = { 1.0, 0.1 }; + + for ( int d = 0; d < n; d++ ) + { + // Refine Min + System.arraycopy( minPoints[ d ], 0, testPoint, 0, n ); + // Search backwards from current min + // We search up to 'initialStride' distance because that's the max error from coarse scan + // But we do it in two passes (step 1, then step 0.1) + + double currentMin = minValid[ d ]; + double range = initialStride; + + for ( double step : steps ) + { + boolean improved = false; + // Search backwards + for ( double val = currentMin - step; val >= currentMin - range; val -= step ) + { + testPoint[ d ] = val; + if ( isValid( testPoint, dimsOther, transformLocal, invTransformOther, globalPoint, otherLocalPoint ) ) + { + minValid[ d ] = val; + // Update the point to this new valid location for next pass + // (though for 1D search on this dim, just keeping the coordinate is enough, + // but keeping other coords fixed is correct behavior for local refinement) + improved = true; + } + else + { + // Once we hit invalid, we stop this pass + // The next finer pass will start from the last valid point (minValid[d]) + break; + } + } + // For the next finer pass, we only need to search a small range around the new min + // specifically, we just searched with 'step', so error is at most 'step'. + // So set range = step for next pass. + // Also update currentMin to the new best found. + currentMin = minValid[ d ]; + range = step; + } + + + // Refine Max + System.arraycopy( maxPoints[ d ], 0, testPoint, 0, n ); + double currentMax = maxValid[ d ]; + range = initialStride; + + for ( double step : steps ) + { + boolean improved = false; + // Search forwards + for ( double val = currentMax + step; val <= currentMax + range; val += step ) + { + testPoint[ d ] = val; + if ( isValid( testPoint, dimsOther, transformLocal, invTransformOther, globalPoint, otherLocalPoint ) ) + { + maxValid[ d ] = val; + improved = true; + } + else + { + break; + } + } + currentMax = maxValid[ d ]; + range = step; + } + } + } + + private static boolean isValid( + final double[] localPoint, + final Dimensions dimsOther, + final AffineTransform3D transformLocal, + final AffineTransform3D invTransformOther, + final double[] globalPoint, // temp buffer + final double[] otherLocalPoint // temp buffer + ) + { + // Transform to global, then to other view's local space + transformLocal.apply( localPoint, globalPoint ); + invTransformOther.apply( globalPoint, otherLocalPoint ); + + // Check if it falls within other view's valid bounds + for ( int d = 0; d < localPoint.length; d++ ) + { + if ( otherLocalPoint[ d ] < 0 || otherLocalPoint[ d ] >= dimsOther.dimension( d ) ) + { + return false; + } + } + return true; + } + + /** + * Sample points in global overlap region and validate they map to valid pixels in both views. + * Returns refined global overlap bounding box containing only valid pixels. + * + * @param globalOverlap Candidate global overlap region + * @param dims1 Dimensions of view 1 + * @param dims2 Dimensions of view 2 + * @param invTransform1 Global-to-local transform for view 1 + * @param invTransform2 Global-to-local transform for view 2 + * @return Refined global overlap, or null if no valid pixels + */ + private static RealInterval sampleAndRefineGlobalOverlap( + final RealInterval globalOverlap, + final Dimensions dims1, + final Dimensions dims2, + final AffineTransform3D invTransform1, + final AffineTransform3D invTransform2 ) + { + final int n = globalOverlap.numDimensions(); + final double[] minValid = new double[ n ]; + final double[] maxValid = new double[ n ]; + + // Initialize with invalid values + for ( int d = 0; d < n; d++ ) + { + minValid[ d ] = Double.MAX_VALUE; + maxValid[ d ] = -Double.MAX_VALUE; + } + + // Sample stride + final int stride = 10; + + // Get sampling bounds in global space + final long[] min = new long[ n ]; + final long[] max = new long[ n ]; + for ( int d = 0; d < n; d++ ) + { + min[ d ] = (long) Math.ceil( globalOverlap.realMin( d ) ); + max[ d ] = (long) Math.floor( globalOverlap.realMax( d ) ); + } + + // Sample and validate points + final double[] globalPoint = new double[ n ]; + final double[] local1Point = new double[ n ]; + final double[] local2Point = new double[ n ]; + + boolean foundValid = false; + + // Sample in 3D grid in global space + for ( long z = min[ 2 ]; z <= max[ 2 ]; z += stride ) + { + for ( long y = min[ 1 ]; y <= max[ 1 ]; y += stride ) + { + for ( long x = min[ 0 ]; x <= max[ 0 ]; x += stride ) + { + globalPoint[ 0 ] = x; + globalPoint[ 1 ] = y; + globalPoint[ 2 ] = z; + + // Transform to both views' local spaces + invTransform1.apply( globalPoint, local1Point ); + invTransform2.apply( globalPoint, local2Point ); + + // Check if it falls within both views' valid bounds + boolean validInView1 = true; + boolean validInView2 = true; + for ( int d = 0; d < n; d++ ) + { + if ( local1Point[ d ] < 0 || local1Point[ d ] >= dims1.dimension( d ) ) + validInView1 = false; + if ( local2Point[ d ] < 0 || local2Point[ d ] >= dims2.dimension( d ) ) + validInView2 = false; + } + + if ( validInView1 && validInView2 ) + { + foundValid = true; + for ( int d = 0; d < n; d++ ) + { + minValid[ d ] = Math.min( minValid[ d ], globalPoint[ d ] ); + maxValid[ d ] = Math.max( maxValid[ d ], globalPoint[ d ] ); + } + } + } + } + } + + if ( !foundValid ) + return null; + + return new FinalRealInterval( minValid, maxValid ); + } + + /** + * Compute overlap between two views by transforming each view's bounding box corners + * into the other view's local space, computing overlaps there, transforming back to + * global space, and intersecting. + * + * This correctly handles rotations by working in each view's local coordinate system. + * + * @param dims1 Dimensions of view 1 + * @param dims2 Dimensions of view 2 + * @param transform1 Local-to-global transform for view 1 + * @param transform2 Local-to-global transform for view 2 + * @return The overlap interval in global coordinates, or null if no overlap + */ + public static RealInterval getOverlapIntervalUsingCorners( + final Dimensions dims1, + final Dimensions dims2, + final AffineTransform3D transform1, + final AffineTransform3D transform2 ) + { + // Get bounding boxes in global space + final RealInterval bbox1Global = getBoundingBoxReal( dims1, transform1 ); + final RealInterval bbox2Global = getBoundingBoxReal( dims2, transform2 ); + + // Generate corners of each view's global bounding box + final double[][] corners1Global = generateCornersFromInterval( bbox1Global ); + final double[][] corners2Global = generateCornersFromInterval( bbox2Global ); + + // Transform bbox1's corners to view 2's local space + final AffineTransform3D invTransform2 = transform2.inverse(); + final double[][] corners1InLocal2 = transformCorners( corners1Global, invTransform2 ); + + // Compute bounding box of transformed corners in view 2's local space + final RealInterval bbox1InLocal2 = getBoundingBoxFromCorners( corners1InLocal2 ); + + // Intersect with view 2's valid image bounds [0, dims2-1] + final RealInterval overlapInLocal2 = intersectWithImageBounds( bbox1InLocal2, dims2 ); + if ( overlapInLocal2 == null ) + return null; + + // Transform just the min/max points of this overlap back to global space + final RealInterval overlap2Global = transformMinMaxToGlobal( overlapInLocal2, transform2 ); + + // Do the same for view 2 → view 1 + final AffineTransform3D invTransform1 = transform1.inverse(); + final double[][] corners2InLocal1 = transformCorners( corners2Global, invTransform1 ); + final RealInterval bbox2InLocal1 = getBoundingBoxFromCorners( corners2InLocal1 ); + final RealInterval overlapInLocal1 = intersectWithImageBounds( bbox2InLocal1, dims1 ); + if ( overlapInLocal1 == null ) + return null; + + // Transform just the min/max points of this overlap back to global space + final RealInterval overlap1Global = transformMinMaxToGlobal( overlapInLocal1, transform1 ); + + // Intersect the two global overlaps + return intersectIntervals( overlap1Global, overlap2Global ); + } + + /** + * Generate corner points from an existing RealInterval. + * + * @param interval The interval to generate corners from + * @return Array of corner points [numCorners][numDimensions] + */ + private static double[][] generateCornersFromInterval( final RealInterval interval ) + { + final int n = interval.numDimensions(); + final int numCorners = (int) Math.pow( 2, n ); + final double[][] corners = new double[ numCorners ][ n ]; + + final double[] min = new double[ n ]; + final double[] max = new double[ n ]; + interval.realMin( min ); + interval.realMax( max ); + + // Generate all 2^n corner combinations using bit pattern + for ( int i = 0; i < numCorners; i++ ) + { + int j = i; + for ( int d = 0; d < n; d++ ) + { + corners[ i ][ d ] = ( j % 2 == 0 ) ? min[ d ] : max[ d ]; + j /= 2; + } + } + + return corners; + } + + /** + * Intersect an interval with valid image bounds [0, dims[d]-1] for each dimension. + * + * @param interval The interval to intersect + * @param dims The image dimensions defining the bounds + * @return The intersected interval, or null if no overlap + */ + private static RealInterval intersectWithImageBounds( final RealInterval interval, final Dimensions dims ) + { + final int n = interval.numDimensions(); + final double[] min = new double[ n ]; + final double[] max = new double[ n ]; + + for ( int d = 0; d < n; d++ ) + { + min[ d ] = Math.max( 0, interval.realMin( d ) ); + max[ d ] = Math.min( dims.dimension( d ) - 1, interval.realMax( d ) ); + + // Check for no overlap + if ( max[ d ] < min[ d ] ) + return null; + } + + return new FinalRealInterval( min, max ); + } + + /** + * Transform just the min and max corner points of an interval to local space. + * This avoids the expansion that would occur if all corners were transformed. + * + * @param globalInterval The interval in global coordinates + * @param invTransform The global-to-local transform + * @return The interval in local coordinates + */ + private static RealInterval transformMinMaxToLocal( final RealInterval globalInterval, final AffineTransform3D invTransform ) + { + final int n = globalInterval.numDimensions(); + final double[] globalMin = new double[ n ]; + final double[] globalMax = new double[ n ]; + final double[] localMin = new double[ n ]; + final double[] localMax = new double[ n ]; + + globalInterval.realMin( globalMin ); + globalInterval.realMax( globalMax ); + + // Transform min and max points + invTransform.apply( globalMin, localMin ); + invTransform.apply( globalMax, localMax ); + + // Ensure min < max (transform might swap them) + final double[] resultMin = new double[ n ]; + final double[] resultMax = new double[ n ]; + for ( int d = 0; d < n; d++ ) + { + resultMin[ d ] = Math.min( localMin[ d ], localMax[ d ] ); + resultMax[ d ] = Math.max( localMin[ d ], localMax[ d ] ); + } + + return new FinalRealInterval( resultMin, resultMax ); + } + + /** + * Transform just the min and max corner points of an interval to global space. + * This avoids the expansion that would occur if all corners were transformed. + * + * @param localInterval The interval in local coordinates + * @param transform The local-to-global transform + * @return The interval in global coordinates + */ + private static RealInterval transformMinMaxToGlobal( final RealInterval localInterval, final AffineTransform3D transform ) + { + final int n = localInterval.numDimensions(); + final double[] localMin = new double[ n ]; + final double[] localMax = new double[ n ]; + final double[] globalMin = new double[ n ]; + final double[] globalMax = new double[ n ]; + + localInterval.realMin( localMin ); + localInterval.realMax( localMax ); + + // Transform min and max points + transform.apply( localMin, globalMin ); + transform.apply( localMax, globalMax ); + + // Ensure min < max (transform might swap them) + final double[] resultMin = new double[ n ]; + final double[] resultMax = new double[ n ]; + for ( int d = 0; d < n; d++ ) + { + resultMin[ d ] = Math.min( globalMin[ d ], globalMax[ d ] ); + resultMax[ d ] = Math.max( globalMin[ d ], globalMax[ d ] ); + } + + return new FinalRealInterval( resultMin, resultMax ); + } + + /** + * Intersect two intervals by taking max of mins and min of maxs. + * + * @param interval1 First interval + * @param interval2 Second interval + * @return The intersection, or null if no overlap + */ + private static RealInterval intersectIntervals( final RealInterval interval1, final RealInterval interval2 ) + { + final int n = interval1.numDimensions(); + final double[] min = new double[ n ]; + final double[] max = new double[ n ]; + + for ( int d = 0; d < n; d++ ) + { + min[ d ] = Math.max( interval1.realMin( d ), interval2.realMin( d ) ); + max[ d ] = Math.min( interval1.realMax( d ), interval2.realMax( d ) ); + + // Check for no overlap + if ( max[ d ] < min[ d ] ) + return null; + } + + return new FinalRealInterval( min, max ); + } + + /** + * Generate all 2^n corner points of an n-dimensional bounding box in local coordinates. + * For a 3D box with dimensions (width, height, depth), generates 8 corners from + * (0,0,0) to (width-1, height-1, depth-1). + * + * @param dims The dimensions of the view in local coordinates + * @return Array of corner points [numCorners][numDimensions] + */ + private static double[][] generateCorners( final Dimensions dims ) + { + final int n = dims.numDimensions(); + final int numCorners = (int) Math.pow( 2, n ); + final double[][] corners = new double[ numCorners ][ n ]; + + final double[] min = new double[ n ]; + final double[] max = new double[ n ]; + + // min is always [0, 0, 0, ...] in local coordinates + // max is [dim0-1, dim1-1, dim2-1, ...] + for ( int d = 0; d < n; d++ ) + { + min[ d ] = 0; + max[ d ] = dims.dimension( d ) - 1; + } + + // Generate all 2^n corner combinations using bit pattern + for ( int i = 0; i < numCorners; i++ ) + { + int j = i; + for ( int d = 0; d < n; d++ ) + { + corners[ i ][ d ] = ( j % 2 == 0 ) ? min[ d ] : max[ d ]; + j /= 2; + } + } + + return corners; + } + + /** + * Transform an array of corner points using the given affine transformation. + * + * @param corners Array of corner points to transform [numCorners][numDimensions] + * @param transform The affine transformation to apply + * @return Array of transformed corner points [numCorners][numDimensions] + */ + private static double[][] transformCorners( final double[][] corners, final AffineTransform3D transform ) + { + final int numCorners = corners.length; + final int n = corners[ 0 ].length; + final double[][] transformed = new double[ numCorners ][ n ]; + + for ( int i = 0; i < numCorners; i++ ) + { + transform.apply( corners[ i ], transformed[ i ] ); + } + + return transformed; + } + + /** + * Calculate the axis-aligned bounding box that encompasses all given corner points. + * Finds the minimum and maximum coordinate values across all corners for each dimension. + * + * @param transformedCorners Array of corner points [numCorners][numDimensions] + * @return RealInterval representing the axis-aligned bounding box + */ + private static RealInterval getBoundingBoxFromCorners( final double[][] transformedCorners ) + { + final int n = transformedCorners[ 0 ].length; + final double[] min = new double[ n ]; + final double[] max = new double[ n ]; + + // Initialize with extreme values + for ( int d = 0; d < n; d++ ) + { + min[ d ] = Double.MAX_VALUE; + max[ d ] = -Double.MAX_VALUE; + } + + // Find min and max across all corners + for ( double[] corner : transformedCorners ) + { + for ( int d = 0; d < n; d++ ) + { + min[ d ] = Math.min( min[ d ], corner[ d ] ); + max[ d ] = Math.max( max[ d ], corner[ d ] ); + } + } + + return new FinalRealInterval( min, max ); + } + + /** + * Transform a global space overlap interval to local coordinates by transforming + * all corners of the overlap box and finding the bounding box of transformed corners. + * This correctly handles non-axis-aligned transformations like rotations. + * + * @param globalOverlap The overlap interval in global coordinates + * @param invTransform The inverse transform (global-to-local) + * @return RealInterval in local coordinates + */ + public static RealInterval transformOverlapToLocal( + final RealInterval globalOverlap, + final AffineTransform3D invTransform ) + { + // Generate corners of the global overlap box + final int n = globalOverlap.numDimensions(); + final int numCorners = (int) Math.pow( 2, n ); + final double[][] globalCorners = new double[ numCorners ][ n ]; + + final double[] gMin = new double[ n ]; + final double[] gMax = new double[ n ]; + globalOverlap.realMin( gMin ); + globalOverlap.realMax( gMax ); + + // Generate all corners of global overlap box + for ( int i = 0; i < numCorners; i++ ) + { + int j = i; + for ( int d = 0; d < n; d++ ) + { + globalCorners[ i ][ d ] = ( j % 2 == 0 ) ? gMin[ d ] : gMax[ d ]; + j /= 2; + } + } + + // Transform all corners to local space + final double[][] localCorners = transformCorners( globalCorners, invTransform ); + + // Find bounding box of transformed corners + return getBoundingBoxFromCorners( localCorners ); + } + + /** + * Convert real-valued local overlap coordinates to integer raster coordinates + * with proper bounds checking. Clamps to [0, dims[d]-1] for each dimension. + * + * @param localOverlap The overlap in local real coordinates + * @param dims The dimensions of the view (for bounds checking) + * @return FinalInterval with integer raster coordinates, or null if no valid overlap + */ + public static FinalInterval getRasterOverlap( + final RealInterval localOverlap, + final Dimensions dims ) + { + final int n = localOverlap.numDimensions(); + final long[] min = new long[ n ]; + final long[] max = new long[ n ]; + + for ( int d = 0; d < n; d++ ) + { + // Use ceiling for min, floor for max to be conservative + min[ d ] = Math.max( 0, (long) Math.ceil( localOverlap.realMin( d ) ) ); + max[ d ] = Math.min( dims.dimension( d ) - 1, (long) Math.floor( localOverlap.realMax( d ) ) ); + + // Validate that we have positive size + if ( max[ d ] < min[ d ] ) + { + return null; // No valid overlap + } + } + + return new FinalInterval( min, max ); + } } From e370c75ab23b8a1bca10f554eab863326203742f Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Tue, 13 Jan 2026 10:23:46 -0500 Subject: [PATCH 14/16] downsample clamp fix --- .../popup/ComputeCrossCorrelationPopup.java | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java index e3388bbf5..e2f122b9f 100644 --- a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java @@ -430,7 +430,17 @@ private static CorrelationResult computeCorrelationForPair( rasterMin2[d], rasterMax2[d], rasterMax2[d] - rasterMin2[d] + 1)); } + // Calculate downsampled image dimensions + long[] dsDims1 = new long[3]; + long[] dsDims2 = new long[3]; + for (int d = 0; d < 3; d++) + { + dsDims1[d] = dims1[d] / downsampleFactors[d]; + dsDims2[d] = dims2[d] / downsampleFactors[d]; + } + // Adjust overlap intervals for downsampling + // Apply downsampling to real coordinates first to maintain alignment long[] dsRasterMin1 = new long[rasterMin1.length]; long[] dsRasterMax1 = new long[rasterMax1.length]; long[] dsRasterMin2 = new long[rasterMin2.length]; @@ -438,10 +448,11 @@ private static CorrelationResult computeCorrelationForPair( for (int d = 0; d < rasterMin1.length; d++) { - dsRasterMin1[d] = rasterMin1[d] / downsampleFactors[d]; - dsRasterMax1[d] = rasterMax1[d] / downsampleFactors[d]; - dsRasterMin2[d] = rasterMin2[d] / downsampleFactors[d]; - dsRasterMax2[d] = rasterMax2[d] / downsampleFactors[d]; + // Calculate directly from real coordinates to preserve alignment + dsRasterMin1[d] = Math.max(0, (long) Math.ceil(localOverlap1.realMin(d) / downsampleFactors[d])); + dsRasterMax1[d] = Math.min(dsDims1[d] - 1, (long) Math.floor(localOverlap1.realMax(d) / downsampleFactors[d])); + dsRasterMin2[d] = Math.max(0, (long) Math.ceil(localOverlap2.realMin(d) / downsampleFactors[d])); + dsRasterMax2[d] = Math.min(dsDims2[d] - 1, (long) Math.floor(localOverlap2.realMax(d) / downsampleFactors[d])); } Interval dsInterval1 = new FinalInterval(dsRasterMin1, dsRasterMax1); From e0b93f6ebbd0d5039690de83c307eb49731f6e41 Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Tue, 13 Jan 2026 11:12:40 -0500 Subject: [PATCH 15/16] changed extension strategy from zero to mirror --- .../popup/ComputeCrossCorrelationPopup.java | 153 ++++++++---------- 1 file changed, 71 insertions(+), 82 deletions(-) diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java index e2f122b9f..c65ed2b48 100644 --- a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java @@ -55,6 +55,8 @@ import net.imglib2.type.numeric.real.FloatType; import net.imglib2.view.Views; import net.imglib2.converter.Converters; +import net.imglib2.interpolation.randomaccess.NLinearInterpolatorFactory; +import net.imglib2.realtransform.RealViews; import net.preibisch.legacy.io.IOFunctions; import net.preibisch.mvrecon.fiji.spimdata.SpimData2; import net.preibisch.mvrecon.fiji.spimdata.explorer.ExplorerWindow; @@ -374,7 +376,7 @@ private static CorrelationResult computeCorrelationForPair( m2 = vr2.getModel().copy(); } - // Calculate local overlaps using pixel validation that handles rotations correctly + // Use pixel validation to find accurate overlap regions that account for rotations RealInterval[] localOverlaps = SimpleBoundingBoxOverlap.getLocalOverlapsUsingPixelValidation( vs1.getSize(), vs2.getSize(), m1, m2); @@ -385,104 +387,68 @@ private static CorrelationResult computeCorrelationForPair( } RealInterval localOverlap1 = localOverlaps[0]; - RealInterval localOverlap2 = localOverlaps[1]; - // Get image dimensions - long[] dims1 = new long[3]; - long[] dims2 = new long[3]; - vs1.getSize().dimensions(dims1); - vs2.getSize().dimensions(dims2); - - // Debug: Print local overlaps - IOFunctions.println(String.format("DEBUG: Analyzing %d-%d <> %d-%d", - viewId1.getTimePointId(), viewId1.getViewSetupId(), - viewId2.getTimePointId(), viewId2.getViewSetupId())); - IOFunctions.println(String.format(" Image sizes: view1=[%d,%d,%d] view2=[%d,%d,%d]", - dims1[0], dims1[1], dims1[2], dims2[0], dims2[1], dims2[2])); + double[] localMin1 = new double[3]; + double[] localMax1 = new double[3]; for (int d = 0; d < 3; d++) { - IOFunctions.println(String.format(" Dim %d local: view1=[%.1f,%.1f] (size=%.1f) view2=[%.1f,%.1f] (size=%.1f)", - d, localOverlap1.realMin(d), localOverlap1.realMax(d), - localOverlap1.realMax(d) - localOverlap1.realMin(d), - localOverlap2.realMin(d), localOverlap2.realMax(d), - localOverlap2.realMax(d) - localOverlap2.realMin(d))); + localMin1[d] = localOverlap1.realMin(d); + localMax1[d] = localOverlap1.realMax(d); } + // Get image dimensions + long[] dims1 = new long[3]; + vs1.getSize().dimensions(dims1); + // Convert to raster coordinates with bounds checking long[] rasterMin1 = new long[3]; long[] rasterMax1 = new long[3]; - long[] rasterMin2 = new long[3]; - long[] rasterMax2 = new long[3]; - - for (int d = 0; d < 3; d++) - { - rasterMin1[d] = Math.max(0, (long) Math.ceil(localOverlap1.realMin(d))); - rasterMax1[d] = Math.min(dims1[d] - 1, (long) Math.floor(localOverlap1.realMax(d))); - rasterMin2[d] = Math.max(0, (long) Math.ceil(localOverlap2.realMin(d))); - rasterMax2[d] = Math.min(dims2[d] - 1, (long) Math.floor(localOverlap2.realMax(d))); - } - // Debug: Print raster coordinates for (int d = 0; d < 3; d++) { - IOFunctions.println(String.format(" Dim %d raster: view1=[%d,%d] (size=%d) view2=[%d,%d] (size=%d)", - d, rasterMin1[d], rasterMax1[d], rasterMax1[d] - rasterMin1[d] + 1, - rasterMin2[d], rasterMax2[d], rasterMax2[d] - rasterMin2[d] + 1)); + rasterMin1[d] = Math.max(0, (long) Math.ceil(localMin1[d])); + rasterMax1[d] = Math.min(dims1[d] - 1, (long) Math.floor(localMax1[d])); } - // Calculate downsampled image dimensions + // Calculate downsampled coordinates + long[] dsRasterMin1 = new long[3]; + long[] dsRasterMax1 = new long[3]; long[] dsDims1 = new long[3]; - long[] dsDims2 = new long[3]; + for (int d = 0; d < 3; d++) { dsDims1[d] = dims1[d] / downsampleFactors[d]; - dsDims2[d] = dims2[d] / downsampleFactors[d]; + dsRasterMin1[d] = Math.max(0, rasterMin1[d] / downsampleFactors[d]); + dsRasterMax1[d] = Math.min(dsDims1[d] - 1, rasterMax1[d] / downsampleFactors[d]); } - // Adjust overlap intervals for downsampling - // Apply downsampling to real coordinates first to maintain alignment - long[] dsRasterMin1 = new long[rasterMin1.length]; - long[] dsRasterMax1 = new long[rasterMax1.length]; - long[] dsRasterMin2 = new long[rasterMin2.length]; - long[] dsRasterMax2 = new long[rasterMax2.length]; + Interval dsInterval1 = new FinalInterval(dsRasterMin1, dsRasterMax1); - for (int d = 0; d < rasterMin1.length; d++) + // Check for zero-sized interval + for (int d = 0; d < 3; d++) { - // Calculate directly from real coordinates to preserve alignment - dsRasterMin1[d] = Math.max(0, (long) Math.ceil(localOverlap1.realMin(d) / downsampleFactors[d])); - dsRasterMax1[d] = Math.min(dsDims1[d] - 1, (long) Math.floor(localOverlap1.realMax(d) / downsampleFactors[d])); - dsRasterMin2[d] = Math.max(0, (long) Math.ceil(localOverlap2.realMin(d) / downsampleFactors[d])); - dsRasterMax2[d] = Math.min(dsDims2[d] - 1, (long) Math.floor(localOverlap2.realMax(d) / downsampleFactors[d])); + if (dsInterval1.dimension(d) <= 0) + { + throw new Exception("Rastered overlap volume is zero"); + } } - Interval dsInterval1 = new FinalInterval(dsRasterMin1, dsRasterMax1); - Interval dsInterval2 = new FinalInterval(dsRasterMin2, dsRasterMax2); - - // check whether we have 0-sized (or negative sized) - // ignore this pair in that case - for ( int d = 0; d < dsInterval1.numDimensions(); ++d ) - { - if ( dsInterval1.dimension( d ) <= 0 || dsInterval2.dimension( d ) <= 0) - { - // Debug logging - IOFunctions.println(String.format("DEBUG: Zero overlap for %d-%d <> %d-%d", - viewId1.getTimePointId(), viewId1.getViewSetupId(), - viewId2.getTimePointId(), viewId2.getViewSetupId())); - IOFunctions.println(String.format(" Dimension %d: view1=[%d,%d] (%d px), view2=[%d,%d] (%d px)", - d, dsRasterMin1[d], dsRasterMax1[d], dsInterval1.dimension(d), - dsRasterMin2[d], dsRasterMax2[d], dsInterval2.dimension(d))); - IOFunctions.println(String.format(" Original raster: view1=[%d,%d], view2=[%d,%d]", - rasterMin1[d], rasterMax1[d], rasterMin2[d], rasterMax2[d])); - IOFunctions.println(String.format(" Downsampling factor: %d", downsampleFactors[d])); - throw new Exception("Rastered overlap volume is zero, skipping." ); - } - } - - - // Load image data at specified downsampling level using pyramid + // Debug output + IOFunctions.println(String.format("DEBUG: Analyzing %d-%d <> %d-%d", + viewId1.getTimePointId(), viewId1.getViewSetupId(), + viewId2.getTimePointId(), viewId2.getViewSetupId())); + IOFunctions.println(String.format(" View1 overlap local: [%.1f,%.1f] [%.1f,%.1f] [%.1f,%.1f]", + localMin1[0], localMax1[0], localMin1[1], localMax1[1], localMin1[2], localMax1[2])); + IOFunctions.println(String.format(" View1 raster: [%d,%d] [%d,%d] [%d,%d] (size: %dx%dx%d)", + rasterMin1[0], rasterMax1[0], rasterMin1[1], rasterMax1[1], rasterMin1[2], rasterMax1[2], + rasterMax1[0]-rasterMin1[0]+1, rasterMax1[1]-rasterMin1[1]+1, rasterMax1[2]-rasterMin1[2]+1)); + IOFunctions.println(String.format(" Downsampled interval: [%d,%d] [%d,%d] [%d,%d] (size: %dx%dx%d)", + dsRasterMin1[0], dsRasterMax1[0], dsRasterMin1[1], dsRasterMax1[1], dsRasterMin1[2], dsRasterMax1[2], + dsInterval1.dimension(0), dsInterval1.dimension(1), dsInterval1.dimension(2))); + + // Load image data at specified downsampling level final BasicImgLoader imgLoader = spimData.getSequenceDescription().getImgLoader(); - // Open images at the specified downsampling level using pre-computed pyramid net.imglib2.util.Pair opened1 = DownsampleTools.openAndDownsample(imgLoader, viewId1, downsampleFactors, false); net.imglib2.util.Pair opened2 = @@ -491,7 +457,17 @@ private static CorrelationResult computeCorrelationForPair( RandomAccessibleInterval img1Raw = opened1.getA(); RandomAccessibleInterval img2Raw = opened2.getA(); - // Convert to FloatType (handles any numeric type) + // Get the downsampled transforms + AffineTransform3D dsTransform1 = opened1.getB(); + AffineTransform3D dsTransform2 = opened2.getB(); + + // Update transforms to account for downsampling + AffineTransform3D dsM1 = m1.copy(); + dsM1.concatenate(dsTransform1); + AffineTransform3D dsM2 = m2.copy(); + dsM2.concatenate(dsTransform2); + + // Convert to FloatType @SuppressWarnings("unchecked") RandomAccessibleInterval img1 = Converters.convert( (RandomAccessibleInterval>) img1Raw, @@ -503,15 +479,28 @@ private static CorrelationResult computeCorrelationForPair( (in, out) -> out.setReal(in.getRealDouble()), new FloatType()); - // Extract overlap regions + // Extract overlap region from view1 RandomAccessibleInterval overlap1 = Views.zeroMin( - Views.interval( - Views.zeroMin(img1), - dsInterval1)); + Views.interval(Views.zeroMin(img1), dsInterval1)); + + // Transform view2 into view1's coordinate system + // The transform is: view2_local -> view2_global -> view1_local + // Which is: invDsM1 * dsM2 + AffineTransform3D view2ToView1 = dsM1.inverse(); + view2ToView1.concatenate(dsM2); + + // Create interpolated view of img2 transformed into view1's space + // Use mirror extension to handle edge cases more gracefully than zero-padding + net.imglib2.RealRandomAccessible interpolated2 = Views.interpolate( + Views.extendMirrorSingle(img2), + new NLinearInterpolatorFactory<>()); + + net.imglib2.RealRandomAccessible transformed2 = + RealViews.transform(interpolated2, view2ToView1); + + // Rasterize the transformed view2 over the same interval as overlap1 RandomAccessibleInterval overlap2 = Views.zeroMin( - Views.interval( - Views.zeroMin(img2), - dsInterval2)); + Views.interval(Views.raster(transformed2), dsInterval1)); // Compute average intensities in the overlapping regions (with sampling for speed) double avgIntensity1 = computeAverageIntensity(overlap1, 10); From c039495ce39428bc4929aa3d63b6762f19c7dad9 Mon Sep 17 00:00:00 2001 From: indecisiveuser Date: Tue, 13 Jan 2026 11:57:10 -0500 Subject: [PATCH 16/16] cleaned up debug statements, updated claude.md --- claude.md | 66 +++++++++++++++++++ .../popup/ComputeCrossCorrelationPopup.java | 13 ---- 2 files changed, 66 insertions(+), 13 deletions(-) diff --git a/claude.md b/claude.md index dabd2b0a3..34dff2b5f 100644 --- a/claude.md +++ b/claude.md @@ -331,6 +331,72 @@ A complete help system was added to the InterestPointExplorer: - `8994af77` - Initial help window implementation - `34bff636` - Fixed focus issue for immediate F1 response +## Cross-Correlation Computation for Rotated Views + +### Problem Statement +Computing accurate cross-correlation values between views with arbitrary affine transformations (especially rotations) requires special handling. Naive approaches that compare axis-aligned bounding boxes in local coordinate systems fail because: + +1. **Simple bounding box intersection** in global space can create overlap regions that extend outside actual image bounds when views are rotated, leading to excessive zero-padding +2. **Independent local overlaps** (sampling in each view separately) find regions with different physical extents that don't actually correspond pixel-by-pixel + +### Solution: Hybrid Approach with Resampling + +The implementation in `ComputeCrossCorrelationPopup.java` uses a hybrid approach combining: + +1. **Accurate Overlap Detection** (via `SimpleBoundingBoxOverlap.getLocalOverlapsUsingPixelValidation`): + - Samples pixels in each view's local space + - Validates they map to valid pixels in the other view + - Accounts for rotations and complex transformations + - Returns tight bounding boxes in each view's local coordinate system + +2. **Proper Data Alignment** (via resampling with interpolation): + - Uses View1's overlap region as the reference coordinate system + - Transforms View2 into View1's coordinate system using the full transform chain + - Applies N-linear interpolation for sub-pixel accuracy + - Uses mirror extension (`extendMirrorSingle`) to handle edge pixels gracefully + +### Implementation Details + +**Key code section** (lines ~379-520 in ComputeCrossCorrelationPopup.java): + +```java +// 1. Find accurate overlaps using pixel validation +RealInterval[] localOverlaps = SimpleBoundingBoxOverlap.getLocalOverlapsUsingPixelValidation( + vs1.getSize(), vs2.getSize(), m1, m2); +RealInterval localOverlap1 = localOverlaps[0]; + +// 2. Extract overlap region from View1 +RandomAccessibleInterval overlap1 = Views.interval(img1, dsInterval1); + +// 3. Transform View2 into View1's coordinate system +AffineTransform3D view2ToView1 = dsM1.inverse(); +view2ToView1.concatenate(dsM2); + +// 4. Create interpolated, transformed view of img2 +RealRandomAccessible interpolated2 = Views.interpolate( + Views.extendMirrorSingle(img2), + new NLinearInterpolatorFactory<>()); +RealRandomAccessible transformed2 = RealViews.transform(interpolated2, view2ToView1); + +// 5. Rasterize over the SAME interval as View1 +RandomAccessibleInterval overlap2 = Views.interval( + Views.raster(transformed2), dsInterval1); + +// 6. Compute correlation on perfectly aligned regions +peak.calculateCrossCorr(overlap1, overlap2); +``` + +### Why This Works + +- **Pixel-by-pixel correspondence**: Both overlap regions use the same interval (from View1's coordinate system), ensuring perfect alignment +- **Rotation handling**: The transformation and interpolation properly handle arbitrary rotations without introducing misalignment +- **Edge handling**: Mirror extension prevents low correlation due to zero-padding at boundaries +- **Efficient**: Only processes the actual overlapping region, not entire images + +### Results + +This approach produces accurate correlation values (typically >0.9 for good alignments) across all view pairs, regardless of relative rotation or transformation complexity. + ## Important Notes - **Never commit without explicit user consent** diff --git a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java index c65ed2b48..d1fce7938 100644 --- a/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java @@ -433,19 +433,6 @@ private static CorrelationResult computeCorrelationForPair( } } - // Debug output - IOFunctions.println(String.format("DEBUG: Analyzing %d-%d <> %d-%d", - viewId1.getTimePointId(), viewId1.getViewSetupId(), - viewId2.getTimePointId(), viewId2.getViewSetupId())); - IOFunctions.println(String.format(" View1 overlap local: [%.1f,%.1f] [%.1f,%.1f] [%.1f,%.1f]", - localMin1[0], localMax1[0], localMin1[1], localMax1[1], localMin1[2], localMax1[2])); - IOFunctions.println(String.format(" View1 raster: [%d,%d] [%d,%d] [%d,%d] (size: %dx%dx%d)", - rasterMin1[0], rasterMax1[0], rasterMin1[1], rasterMax1[1], rasterMin1[2], rasterMax1[2], - rasterMax1[0]-rasterMin1[0]+1, rasterMax1[1]-rasterMin1[1]+1, rasterMax1[2]-rasterMin1[2]+1)); - IOFunctions.println(String.format(" Downsampled interval: [%d,%d] [%d,%d] [%d,%d] (size: %dx%dx%d)", - dsRasterMin1[0], dsRasterMax1[0], dsRasterMin1[1], dsRasterMax1[1], dsRasterMin1[2], dsRasterMax1[2], - dsInterval1.dimension(0), dsInterval1.dimension(1), dsInterval1.dimension(2))); - // Load image data at specified downsampling level final BasicImgLoader imgLoader = spimData.getSequenceDescription().getImgLoader();