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/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/ViewSetupExplorerPanel.java b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/ViewSetupExplorerPanel.java index 67dabc81f..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 @@ -86,6 +86,7 @@ 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; @@ -694,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() ); 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..d1fce7938 --- /dev/null +++ b/src/main/java/net/preibisch/mvrecon/fiji/spimdata/explorer/popup/ComputeCrossCorrelationPopup.java @@ -0,0 +1,641 @@ +/*- + * #%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.ArrayList; +import java.util.Collections; +import java.util.Comparator; +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 java.util.stream.Collectors; + +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.Channel; +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.RealType; +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; +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; + + /** + * 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"); + 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(); + + // Get selected groups for potential use + List> selectedGroups = ((GroupedRowWindow) panel).selectedRowsViewIdGroups(); + + // Get available downsampling levels + List allViewIds = new ArrayList<>(spimData.getSequenceDescription().getViewDescriptions().keySet()); + if (allViewIds.isEmpty()) + { + IOFunctions.println(new Date(System.currentTimeMillis()) + + ": ERROR: No views found in dataset"); + return; + } + + 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.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()) + { + IOFunctions.println(new Date(System.currentTimeMillis()) + + ": Cross-correlation computation cancelled by user"); + return; + } + + 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); + + // Get the list of views to process + List viewIdsToProcess; + + if (selectedOnly) + { + // 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 + viewIdsToProcess = new ArrayList<>(); + for (List group : selectedGroups) + { + SpimData2.filterMissingViews(spimData, group); + if (!group.isEmpty()) + { + // Get first ViewId from each group + viewIdsToProcess.add(group.get(0)); + } + } + + if (viewIdsToProcess.size() < 2) + { + IOFunctions.println(new Date(System.currentTimeMillis()) + + ": ERROR: Selected groups contain fewer than 2 valid views"); + return; + } + + IOFunctions.println(" Processing " + viewIdsToProcess.size() + " selected views"); + } + else + { + // Get all views + viewIdsToProcess = new ArrayList<>(); + for (ViewId viewId : allViewIds) + { + if (!spimData.getSequenceDescription().getMissingViews().getMissingViews().contains(viewId)) + { + viewIdsToProcess.add(viewId); + } + } + + IOFunctions.println(" Found " + viewIdsToProcess.size() + " valid views"); + } + + // Find all overlapping pairs + List pairsToProcess = findOverlappingPairs( + spimData, viewIdsToProcess, chChoice); + + 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()); + 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"); + + List> futures = new ArrayList<>(); + + for (final ViewId[] pair : pairsToProcess) + { + Future future = executor.submit(new Runnable() + { + @Override + public void run() + { + try + { + CorrelationResult result = computeCorrelationForPair(spimData, pair[0], pair[1], downsampleFactors); + if (result != null) + { + results.add(result); + successCount.incrementAndGet(); + } + } + catch (Exception ex) + { +// // Temporarily log first error to diagnose issue +// if (successCount.get() == 0) +// { +// IOFunctions.println(new Date(System.currentTimeMillis()) + +// ": ERROR (first failure): " + ex.getMessage()); +// ex.printStackTrace(); +// } + } + } + }); + futures.add(future); + } + + // Wait for all tasks to complete + for (Future future : futures) + { + try + { + future.get(); + } + catch (Exception ex) + { + // Already handled in the task + } + } + + // 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) + { + IOFunctions.println(new Date(System.currentTimeMillis()) + + ": ERROR: Exception during correlation computation: " + ex.getMessage()); + ex.printStackTrace(); + } + } + }).start(); + } + } + + /** + * Compute cross-correlation for a specific pair of views + */ + private static CorrelationResult computeCorrelationForPair( + SpimData2 spimData, + ViewId viewId1, + ViewId viewId2, + long[] downsampleFactors + ) throws Exception + { + // 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); + + // 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(); + } + + // Use pixel validation to find accurate overlap regions that account for rotations + RealInterval[] localOverlaps = SimpleBoundingBoxOverlap.getLocalOverlapsUsingPixelValidation( + vs1.getSize(), vs2.getSize(), + m1, m2); + + if (localOverlaps == null) + { + throw new Exception("No overlap found"); + } + + RealInterval localOverlap1 = localOverlaps[0]; + + double[] localMin1 = new double[3]; + double[] localMax1 = new double[3]; + for (int d = 0; d < 3; 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]; + + 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])); + } + + // Calculate downsampled coordinates + long[] dsRasterMin1 = new long[3]; + long[] dsRasterMax1 = new long[3]; + long[] dsDims1 = new long[3]; + + for (int d = 0; d < 3; d++) + { + dsDims1[d] = dims1[d] / downsampleFactors[d]; + dsRasterMin1[d] = Math.max(0, rasterMin1[d] / downsampleFactors[d]); + dsRasterMax1[d] = Math.min(dsDims1[d] - 1, rasterMax1[d] / downsampleFactors[d]); + } + + Interval dsInterval1 = new FinalInterval(dsRasterMin1, dsRasterMax1); + + // Check for zero-sized interval + for (int d = 0; d < 3; d++) + { + if (dsInterval1.dimension(d) <= 0) + { + throw new Exception("Rastered overlap volume is zero"); + } + } + + // Load image data at specified downsampling level + final BasicImgLoader imgLoader = spimData.getSequenceDescription().getImgLoader(); + + net.imglib2.util.Pair opened1 = + DownsampleTools.openAndDownsample(imgLoader, viewId1, downsampleFactors, false); + net.imglib2.util.Pair opened2 = + DownsampleTools.openAndDownsample(imgLoader, viewId2, downsampleFactors, false); + + RandomAccessibleInterval img1Raw = opened1.getA(); + RandomAccessibleInterval img2Raw = opened2.getA(); + + // 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, + (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 region from view1 + RandomAccessibleInterval overlap1 = Views.zeroMin( + 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.raster(transformed2), dsInterval1)); + + // 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()); + 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); + } + + // 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); + } + + /** + * 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; + } + + /** + * 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; + } + + /** + * 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; + } +} 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..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; @@ -164,10 +166,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[ 3 ]; + final int[] maxInt = new int[ 3 ]; - for ( int d = 0; d < dims.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; @@ -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 ); + } } 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 ); + + } + + +}