diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000000000000000000000000000000000000..a83e5ab37d01418d3c7eae627b384e347b13b4ee
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,7 @@
+.idea/
+.settings/
+*.iml
+.project
+.classpath
+target/
+build/
\ No newline at end of file
diff --git a/pom.xml b/pom.xml
new file mode 100644
index 0000000000000000000000000000000000000000..9db049ec44c38f98dae27d8b7ec026c929b0662f
--- /dev/null
+++ b/pom.xml
@@ -0,0 +1,96 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<project xmlns="http://maven.apache.org/POM/4.0.0"
+ xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
+ xsi:schemaLocation="http://maven.apache.org/POM/4.0.0 http://maven.apache.org/xsd/maven-4.0.0.xsd">
+ <modelVersion>4.0.0</modelVersion>
+
+ <!-- Inherited Icy Parent POM -->
+ <parent>
+ <groupId>org.bioimageanalysis.icy</groupId>
+ <artifactId>parent-pom-plugin</artifactId>
+ <version>1.0.3</version>
+ </parent>
+
+ <!-- Project Information -->
+ <artifactId>isotropic-wavelet-transform</artifactId>
+ <version>1.1.0</version>
+
+ <packaging>jar</packaging>
+
+ <name>Isotropic Wavelet Trabsform</name>
+ <description>
+ Isotropic wavelets transform for 2D images. The plugin provides a multiscale representation of 2D images without directional bias.
+ The transform defines a 'tight frame' so that images are easily reconstructed from wavelet coefficients. See: Chenouard, N.; Unser,
+ M., "3D Steerable Wavelets in Practice," Image Processing, IEEE Transactions on , vol.21, no.11, pp.4522,4533, Nov. 2012, for a
+ complete description in 3D.
+ Chenouard, N.; Unser, M., "3D Steerable Wavelets in Practice," Image Processing, IEEE Transactions on , vol.21, no.11, pp.4522,4533,
+ Nov. 2012: http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=6226458&tag=1
+ </description>
+ <url>http://icy.bioimageanalysis.org/plugin/isotropic-wavelet-transform/</url>
+ <inceptionYear>2020</inceptionYear>
+
+ <organization>
+ <name>Institut Pasteur</name>
+ <url>https://pasteur.fr</url>
+ </organization>
+
+ <licenses>
+ <license>
+ <name>GNU GPLv3</name>
+ <url>https://www.gnu.org/licenses/gpl-3.0.en.html</url>
+ <distribution>repo</distribution>
+ </license>
+ </licenses>
+
+ <developers>
+ <developer>
+ <id>sdallongeville</id>
+ <name>Stéphane Dallongeville</name>
+ <url>https://research.pasteur.fr/fr/member/stephane-dallongeville/</url>
+ <roles>
+ <role>founder</role>
+ <role>lead</role>
+ <role>architect</role>
+ <role>developer</role>
+ <role>debugger</role>
+ <role>tester</role>
+ <role>maintainer</role>
+ <role>support</role>
+ </roles>
+ </developer>
+ </developers>
+
+ <!-- Project properties -->
+ <properties>
+
+ </properties>
+
+ <!-- Project build configuration -->
+ <build>
+
+ </build>
+
+ <!-- List of project's dependencies -->
+ <dependencies>
+ <!-- The core of Icy -->
+ <dependency>
+ <groupId>org.bioimageanalysis.icy</groupId>
+ <artifactId>icy-kernel</artifactId>
+ </dependency>
+
+ <dependency>
+ <groupId>edu.emory.mathcs</groupId>
+ <artifactId>JTransforms</artifactId>
+ <version>2.4</version>
+ </dependency>
+ </dependencies>
+
+ <!-- Icy Maven repository (to find parent POM) -->
+ <repositories>
+ <repository>
+ <id>icy</id>
+ <name>Icy's Nexus</name>
+ <url>https://icy-nexus.pasteur.fr/repository/Icy/</url>
+ </repository>
+ </repositories>
+</project>
\ No newline at end of file
diff --git a/src/main/java/plugins/nchenouard/isotropicwavelets/CoefficientSelectionPanel.java b/src/main/java/plugins/nchenouard/isotropicwavelets/CoefficientSelectionPanel.java
new file mode 100644
index 0000000000000000000000000000000000000000..784f9141c0094ad2f9423e00fba148c5552dd5cb
--- /dev/null
+++ b/src/main/java/plugins/nchenouard/isotropicwavelets/CoefficientSelectionPanel.java
@@ -0,0 +1,116 @@
+package plugins.nchenouard.isotropicwavelets;
+
+import icy.gui.component.pool.SwimmingObjectChooser;
+import icy.swimmingPool.SwimmingObject;
+
+import java.awt.BorderLayout;
+import java.awt.event.ItemEvent;
+import java.awt.event.ItemListener;
+
+import javax.swing.JPanel;
+import javax.swing.JScrollPane;
+import javax.swing.JTable;
+import javax.swing.table.AbstractTableModel;
+
+/**
+ *
+ * GUI for wavelet coefficients manual selection
+ *
+ * @author Nicolas Chenouard (nicolas.chenouard.dev@gmail.com)
+ * @version 1.0
+ * 2014-05-22
+ * gpl v3.0
+ */
+
+
+public class CoefficientSelectionPanel extends JPanel
+{
+ private static final long serialVersionUID = -4324385449999084787L;
+ CoefficientsTableModel coeffTableModel;
+ JTable coefficientTable;
+ SwimmingObjectChooser resultsBox;
+
+ public CoefficientSelectionPanel()
+ {
+ this.setLayout(new BorderLayout());
+ resultsBox = new SwimmingObjectChooser(SequenceAnalysisResults.class);
+ this.add(resultsBox, BorderLayout.NORTH);
+ resultsBox.addItemListener(new ItemListener() {
+ @Override
+ public void itemStateChanged(ItemEvent arg0) {
+ updateResultsTable();
+ }
+ });
+
+ coeffTableModel = new CoefficientsTableModel();
+ coefficientTable = new JTable(coeffTableModel);
+ JScrollPane scrollTablePane = new JScrollPane(coefficientTable);
+ this.add(scrollTablePane, BorderLayout.CENTER);
+ }
+
+ void updateResultsTable()
+ {
+ Object selectedItem = resultsBox.getSelectedItem();
+ if (selectedItem != null)
+ {
+ SequenceAnalysisResults results = (SequenceAnalysisResults) ((SwimmingObject) selectedItem).getObject();
+ coeffTableModel.setResults(results);
+ }
+ else
+ coeffTableModel.setResults(null);
+ }
+
+ /**Table model for wavelet coefficients*/
+ public class CoefficientsTableModel extends AbstractTableModel
+ {
+ private static final long serialVersionUID = -4810024326298005026L;
+ SequenceAnalysisResults results = null;
+ String[] columnNames = new String[]{"Sequence name", "Frames", "Scales", "Wavelet profile", "Filter high frequencies", "Isotropic padding", "Padding X", "Padding Y"};
+ String[] valueList = new String[]{"", "", "", "","","","",""};
+
+ @Override
+ public String getColumnName(int column)
+ {
+ return columnNames[column];
+ };
+
+ public void setResults(SequenceAnalysisResults r)
+ {
+ results = r ;
+ if (results == null)
+ {
+ for (int i = 0; i < valueList.length; i++)
+ valueList[i] = "";
+ }
+ else
+ {
+ valueList[0] = results.getSequenceName();
+ valueList[1] = Integer.toString(results.getAllAnalyzedTimesResults().size());
+ WaveletAnalysisResults wr = results.getAllResults().get(0);
+ WaveletFilterSet filters = wr.getWaveletFilters();
+ valueList[2] = Integer.toString(filters.getNumScales());
+ valueList[3] = filters.getWaveletType().toString();
+ valueList[4] = Boolean.toString(filters.isPrefilter());
+ valueList[5] = Boolean.toString(filters.isIsotropicPadding());
+ valueList[6] = Integer.toString(wr.getPadX());
+ valueList[7] = Integer.toString(wr.getPadY());
+ }
+ this.fireTableDataChanged();
+ }
+
+ @Override
+ public int getColumnCount() {
+ return columnNames.length;
+ }
+
+ @Override
+ public int getRowCount() {
+ return 1;
+ }
+
+ @Override
+ public Object getValueAt(int arg0, int arg1) {
+ return valueList[arg1];
+ }
+ }
+}
\ No newline at end of file
diff --git a/src/main/java/plugins/nchenouard/isotropicwavelets/IsotropicWaveletTransform.java b/src/main/java/plugins/nchenouard/isotropicwavelets/IsotropicWaveletTransform.java
new file mode 100644
index 0000000000000000000000000000000000000000..fe4c7396af9bdb4306a4ec06b9583fdaaf575aa4
--- /dev/null
+++ b/src/main/java/plugins/nchenouard/isotropicwavelets/IsotropicWaveletTransform.java
@@ -0,0 +1,410 @@
+package plugins.nchenouard.isotropicwavelets;
+
+import icy.image.IcyBufferedImage;
+import icy.main.Icy;
+import icy.sequence.Sequence;
+import icy.type.collection.array.ArrayUtil;
+import edu.emory.mathcs.jtransforms.fft.DoubleFFT_2D;
+
+/**
+ *
+ * Isotropic wavelet transform for 2D images
+ *
+ * 2D wavelet decomposition and reconstruction algorithms are provided.
+ * Transforms are as decribed in:
+ *
+ * Chenouard, N.; Unser, M., "3D Steerable Wavelets in Practice"
+ * IEEE Transactions on Image Processing, vol.21, no.11, pp.4522,4533, Nov. 2012
+ * doi: 10.1109/TIP.2012.2206044
+ *
+ * except being implemented in 2D instead of 3D.
+ *
+ * Please cite the above reference in scientific communications upon use of these tools.
+ *
+ *
+ * @author Nicolas Chenouard (nicolas.chenouard.dev@gmail.com)
+ * @version 1.0
+ * 2014-05-22
+ * gpl v3.0
+ */
+
+public class IsotropicWaveletTransform
+{
+ /**
+ * @return Compute x and y direction padding such that the original image is centered in the padded image
+ *
+ * @param imageWidth the width of the original image
+ * @param imageHeight the height of the original image
+ * @param width the width of the image after padding
+ * @param height the height of the image after padding
+ *
+ * */
+ public static int[] computePadding(int imageWidth, int imageHeight, int width, int height)
+ {
+ int padY = (int) (height - imageHeight)/2;
+ int padX = (int) (width - imageWidth)/2;
+ return new int[]{padX, padY};
+ }
+
+ /**
+ * @return Pad an image with zero values
+ *
+ * @param image the 2D image to pad as a 1D array
+ * @param imageWidth the width of the original image
+ * @param imageHeight the height of the original image
+ * @param width the width of the image after padding
+ * @param height the height of the image after padding
+ * @param padX the size of the pad in x direction, on the left of the image
+ * @param padY the size of the pad in x direction, on the top of the image
+ * */
+ public static double[] padImage(double[] image, int imageWidth, int imageHeight, int width, int height, int padX, int padY)
+ {
+ if (imageWidth == width && imageHeight == height)
+ return image;
+ double[] paddedImage = new double[width*height];
+ for (int y = padY; y < padY + imageHeight; y++)
+ {
+ System.arraycopy(image, (y - padY)*imageWidth, paddedImage, y*width + padX, imageWidth);
+ }
+ return paddedImage;
+ }
+
+ /**
+ * @return rop a padded image to its original dimensions
+ *
+ * @param image the padded 2D image as a 1D array
+ * @param croppedImageWidth the width of the image after cropping
+ * @param croppedImageHeight the height of the image after cropping
+ * @param width the width of the padded image
+ * @param height the height of the padded image
+ * @param padX the size of the pad in x direction, on the left of the image
+ * @param padY the size of the pad in x direction, on the top of the image
+ * */
+ public static double[] unpadImage(double[] image, int croppedImageWidth, int croppedImageHeight, int width, int height, int padX, int padY)
+ {
+ if (croppedImageWidth == width && croppedImageHeight == height)
+ return image;
+ double[] croppedImage = new double[croppedImageWidth*croppedImageHeight];
+ for (int y = padY; y < padY + croppedImageHeight; y++)
+ {
+ System.arraycopy(image, y*width + padX, croppedImage, (y - padY)*croppedImageWidth, croppedImageWidth);
+ }
+ return croppedImage;
+ }
+
+ /**
+ * @return Downsample the Fourier representation of a 2D image by a factor 2
+ * @param dataFFT Fourier coefficients of the image
+ * @param scaleWidth width of the image
+ * @param scaleHeight height of the image
+ * */
+ static public double[] downSampleFourierDomain(double[] dataFFT, int scaleWidth, int scaleHeight)
+ {
+ // %downSampling in spatial domain = fold the frequency domain
+ // c2 = size(Im,2)/2;
+ // c1 = size(Im,1)/2;
+ // Im = 0.25*(Im(1:c1, 1:c2) + Im((1:c1)+c1, 1:c2) + Im((1:c1) + c1, (1:c2) +c2) + Im(1:c1, (1:c2) +c2));
+
+ int halfWidth = (int)(scaleWidth/2);
+ int halfHeight = (int)(scaleHeight/2);
+ double[] downSampledData = new double[2*halfWidth*halfHeight];
+ for (int y = 0; y < halfHeight; y++)
+ for (int x = 0; x < halfWidth; x++)
+ {
+ downSampledData[2*(y*halfWidth + x)] = 0.25f*(dataFFT[2*(y* scaleWidth + x)] + dataFFT[2*(y* scaleWidth + x + halfWidth)] + dataFFT[2*((y + halfHeight)* scaleWidth + x)] + dataFFT[2*((y + halfHeight) * scaleWidth + x + halfWidth)]);
+ downSampledData[2*(y*halfWidth + x) + 1] = 0.25f*(dataFFT[2*(y* scaleWidth + x) + 1] + dataFFT[2*(y* scaleWidth + x + halfWidth) + 1] + dataFFT[2*((y + halfHeight)* scaleWidth + x) + 1] + dataFFT[2*((y + halfHeight) * scaleWidth + x + halfWidth) + 1]);
+ }
+ return downSampledData;
+ }
+
+ /**
+ * Convolve two real images using their Fourier representation and write results at the first image location.
+ * @param dataFFT Fourier coefficients of the first image. Its values are overrided by the convolution results.
+ * @param prefilterRealFFT Fourirer of the second image
+ * */
+ public static void convolveRealFFTFilterInPlace(double[] dataFFT, double[] prefilterRealFFT) {
+ for (int i = 0; i < prefilterRealFFT.length; i++)
+ {
+ dataFFT[2*i] = dataFFT[2*i]*prefilterRealFFT[i];
+ dataFFT[2*i + 1] = dataFFT[2*i + 1]*prefilterRealFFT[i];
+ }
+ }
+
+ /**
+ * Convolve two real images using their Fourier representation.
+ * @param dataFFT Fourier coefficients of the first image
+ * @param prefilterRealFFT Fourirer of the second image
+ * @return Fourier coefficients of the convolution results
+ * */
+ public static double[] convolveRealFFTFilter(double[] dataFFT, double[] prefilterRealFFT) {
+ double[] dataCopy = new double[dataFFT.length];
+ System.arraycopy(dataFFT, 0, dataCopy, 0, dataCopy.length);
+ for (int i = 0; i < prefilterRealFFT.length; i++)
+ {
+ dataCopy[2*i] = dataFFT[2*i]*prefilterRealFFT[i];
+ dataCopy[2*i + 1] = dataFFT[2*i + 1]*prefilterRealFFT[i];
+ }
+ return dataCopy;
+ }
+
+ /**
+ * Reconstruct a 2D image from a set of wavelet coefficients
+ * @param waveletResults wavelet coefficients
+ * @return array corresponding to the reconstructed image
+ * */
+ public static double[] isotropicBandlimitedSynthesis(WaveletAnalysisResults waveletResults)
+ {
+ WaveletFilterSet filterSet = waveletResults.getWaveletFilters();
+ int scaleWidth = filterSet.getLPWidth();
+ int scaleHeight = filterSet.getLPHeight();
+
+ double[] lpFFT;
+ if (waveletResults.storeCoefficientsInFourierDomain)
+ {
+ lpFFT = waveletResults.lpResidual;
+ }
+ else
+ {
+ int numCoefficients = filterSet.getLPHeight()*filterSet.getLPWidth();
+ lpFFT = new double[2*numCoefficients];
+ for (int i = 0; i < numCoefficients; i++)
+ lpFFT[2*i] = waveletResults.lpResidual[i];
+ DoubleFFT_2D fftScale = new DoubleFFT_2D(scaleHeight, scaleWidth);
+ fftScale.complexForward(lpFFT);
+ }
+ for (int iter = filterSet.getNumScales() - 1; iter >= 0; iter--)
+ {
+ // upsampling of the lowpass image in Fourier domain
+ lpFFT = upsampleFFT(lpFFT, scaleWidth, scaleHeight);
+ // compute FFT for the wavelet band
+ scaleWidth = filterSet.getScaleWidth(iter);
+ scaleHeight = filterSet.getScaleHeight(iter);
+ double[] hpFFT;
+ if (waveletResults.storeCoefficientsInFourierDomain)
+ {
+ hpFFT = waveletResults.getWaveletBand(iter);
+ }
+ else
+ {
+ double[] scale = waveletResults.getWaveletBand(iter);
+ int numCoefficients = filterSet.getScaleWidth(iter)*filterSet.getScaleHeight(iter);
+ hpFFT = new double[2*numCoefficients];
+ for (int i = 0; i < numCoefficients; i++)
+ hpFFT[2*i] = scale[i];
+ DoubleFFT_2D fftScale = new DoubleFFT_2D(scaleHeight, scaleWidth);
+ fftScale.complexForward(hpFFT);
+ }
+ // retrieve the filters in fourier domain
+ double[][] filters = WaveletFunction.getFilterFFT(filterSet.getWaveletType(), scaleWidth, scaleHeight);
+ // combine low pass and high pass images
+ for(int i = 0; i < (int)(hpFFT.length/2); i++)
+ {
+ lpFFT[2*i] = 4*lpFFT[2*i]*filters[0][i] + hpFFT[2*i]*filters[1][i];
+ lpFFT[2*i + 1] = 4*lpFFT[2*i + 1]*filters[0][i] + hpFFT[2*i + 1]*filters[1][i];
+ }
+ }
+ if (filterSet.isPrefilter())
+ {
+ double[][] prefiltersFFT = WaveletFunction.getPreFiltersFFT(filterSet.waveletType, filterSet.getHPWidth(), filterSet.getHPHeight());
+ double[] hpFFT;
+ if (waveletResults.storeCoefficientsInFourierDomain)
+ {
+ //hpFFT = waveletResults.getHPResidual();
+ hpFFT = waveletResults.getHPResidual().clone();
+ }
+ else
+ {
+ int numCoefficients = filterSet.getHPWidth()*filterSet.getHPHeight();
+ hpFFT = new double[2*numCoefficients];
+ for (int i = 0; i < numCoefficients; i++)
+ hpFFT[2*i] = waveletResults.hpResidual[i];
+ DoubleFFT_2D fftScale = new DoubleFFT_2D(filterSet.getHPHeight(), filterSet.getHPWidth());
+ fftScale.complexForward(hpFFT);
+ }
+ convolveRealFFTFilterInPlace(hpFFT, prefiltersFFT[1]);
+ convolveRealFFTFilterInPlace(lpFFT, prefiltersFFT[0]);
+ for (int i = 0; i < lpFFT.length; i++)
+ lpFFT[i] = lpFFT[i] + hpFFT[i];
+ }
+ // project back to spatial domain
+ DoubleFFT_2D fftScale = new DoubleFFT_2D(filterSet.getHeight(), filterSet.getWidth());
+ fftScale.complexInverse(lpFFT, true);
+ double[] reconstruction = new double[lpFFT.length/2];
+ for (int i = 0; i < reconstruction.length; i++)
+ reconstruction[i] = lpFFT[2*i];
+ return reconstruction;
+ }
+
+ /**
+ * Up sample a 2D image by duplicating its Fourier coefficients
+ * @param dataFFT Fourier coefficients of the image to upsample
+ * @param width width of the image
+ * @param height height of the image
+ * @return an array containing the Fourier coefficients of the upsampled image
+ *
+ * */
+ static public double[] upsampleFFT(double[] dataFFT, int width, int height) {
+ // upsampling by a factor two corresponds to duplicating the fourier plane
+ double[] upsampledFFT = new double[4*dataFFT.length];
+ for (int y = 0; y < height; y++)
+ {
+ System.arraycopy(dataFFT, 2*y*width, upsampledFFT, 2*y*(2*width), 2*width);
+ System.arraycopy(dataFFT, 2*y*width, upsampledFFT, 2*(y*(2*width) + width), 2*width);
+ System.arraycopy(dataFFT, 2*y*width, upsampledFFT, 2*(y + height)*(2*width), 2*width);
+ System.arraycopy(dataFFT, 2*y*width, upsampledFFT, 2*((y + height)*(2*width) + width), 2*width);
+ }
+ return upsampledFFT;
+ }
+
+ /**
+ * @return Perform wavelet decomposition of a 2D image
+ * @param image the 2D image to decompose as a 1D array pixel value at (x, y) = image[x +y*width]
+ * @param imageWidth width of the 2D image
+ * @param imageHeight height of the 2D image
+ * @param waveletFilters set of wavelet filters to use for the decomposition
+ * */
+ public static WaveletAnalysisResults isotropicBandlimitedAnalysis(double[] image, int imageWidth, int imageHeight, WaveletFilterSet waveletFilters)
+ {
+ return isotropicBandlimitedAnalysis(image, imageWidth, imageHeight, waveletFilters, false);
+ }
+
+ /**
+ * @return Perform wavelet decomposition of a 2D image
+ * @param image the 2D image to decompose as a 1D array pixel value at (x, y) = image[x +y*width]
+ * @param imageWidth width of the 2D image
+ * @param imageHeight height of the 2D image
+ * @param waveletFilters set of wavelet filters to use for the decomposition
+ * @param coefficientsInFourierDomain true if wavelet coefficients are to be stored in the Fourier domain
+ * */
+ public static WaveletAnalysisResults isotropicBandlimitedAnalysis(double[] image, int imageWidth, int imageHeight, WaveletFilterSet waveletFilters, boolean coefficientsInFourierDomain)
+ {
+ WaveletAnalysisResults results = new WaveletAnalysisResults();
+ results.filterSet = waveletFilters;
+ results.storeCoefficientsInFourierDomain = coefficientsInFourierDomain;
+ results.bands = new double[waveletFilters.getNumScales()][];
+
+ int width = waveletFilters.width;
+ int height = waveletFilters.height;
+
+ // pad image to fit filter size
+ int [] pad = computePadding(imageWidth, imageHeight, width, height);
+ results.padX = pad[0];
+ results.padY = pad[1];
+ results.fullWidth = imageWidth;
+ results.fullHeight = imageHeight;
+ image = padImage(image, imageWidth, imageHeight, width, height, pad[0], pad[1]);
+ //Sequence padSequence = new Sequence("pad");
+ //padSequence.addImage(new IcyBufferedImage(width, height, image));
+ //Icy.getMainInterface().addSequence(padSequence);
+
+ double[] dataFFT_t = new double[image.length*2];
+ for (int i = 0; i < image.length; i++)
+ dataFFT_t[2*i] = image[i];
+ DoubleFFT_2D fft = new DoubleFFT_2D(height, width);
+ fft.complexForward(dataFFT_t);
+ if (waveletFilters.prefilter)
+ {
+ double[][] prefiltersFFT = waveletFilters.prefilterFFT;
+ double[] hpFFT = convolveRealFFTFilter(dataFFT_t, prefiltersFFT[1]);
+ convolveRealFFTFilterInPlace(dataFFT_t, prefiltersFFT[0]);
+ //results.hpSize = new int[]{width, height};
+ if (coefficientsInFourierDomain)
+ {
+ results.hpResidual = hpFFT;
+ }
+ else
+ {
+ DoubleFFT_2D fftScale = new DoubleFFT_2D(height, width);
+ fftScale.complexInverse(hpFFT, true);
+ results.hpResidual = new double[height*width];
+ for (int i = 0; i < results.hpResidual.length; i++)
+ results.hpResidual[i] = hpFFT[2*i];
+ }
+ }
+ for (int scale = 0; scale < waveletFilters.numScales; scale++)
+ {
+ int scaleWidth = waveletFilters.getScaleWidth(scale);
+ int scaleHeight = waveletFilters.getScaleHeight(scale);
+ double[][] filters = waveletFilters.waveletFiltersFFT.get(scale);
+ double[] fftHP = convolveRealFFTFilter(dataFFT_t, filters[1]);
+ convolveRealFFTFilterInPlace(dataFFT_t, filters[0]);
+ dataFFT_t = downSampleFourierDomain(dataFFT_t, scaleWidth, scaleHeight);
+ if (coefficientsInFourierDomain)
+ {
+ results.bands[scale] = fftHP;
+ }
+ else
+ {
+ DoubleFFT_2D fftScale = new DoubleFFT_2D(scaleHeight, scaleWidth);
+ fftScale.complexInverse(fftHP, true);
+ double[] hp = new double[scaleWidth*scaleHeight];
+ for (int i = 0; i < hp.length; i++)
+ hp[i] = fftHP[2*i];
+ results.bands[scale] = hp;
+ }
+ }
+ // project back the low pass residual
+ if (coefficientsInFourierDomain)
+ {
+ results.lpResidual = dataFFT_t;
+ }
+ else
+ {
+ int scaleHeight = waveletFilters.getLPHeight();
+ int scaleWidth = waveletFilters.getLPWidth();
+ DoubleFFT_2D fftScale = new DoubleFFT_2D(scaleHeight, scaleWidth);
+ fftScale.complexInverse(dataFFT_t, true);
+ results.lpResidual = new double[scaleHeight*scaleWidth];
+ for (int i = 0; i < results.lpResidual.length; i++)
+ results.lpResidual[i] = dataFFT_t[2*i];
+ }
+ return results;
+ }
+
+ /**
+ * Test function of Fourier domain decomposition and reconstruction.
+ * */
+ protected static void testFFT()
+ {
+ Sequence seq = Icy.getMainInterface().getActiveSequence();
+ if (seq == null)
+ return;
+ int width = seq.getSizeX();
+ int height = seq.getSizeY();
+ // just show FFT for the first image of the sequence
+ DoubleFFT_2D fft = new DoubleFFT_2D(height, width);
+ IcyBufferedImage image = seq.getFirstImage();
+ double[] data = (double[])ArrayUtil.arrayToDoubleArray(image.getDataXY(0), image.isSignedDataType() );
+ double[] dataFFT = new double[data.length*2];
+ for (int i = 0; i < data.length; i++)
+ dataFFT[2*i] = data[i];
+ fft.complexForward(dataFFT);
+ Sequence resSeq = new Sequence();
+ resSeq.setName("FFT");
+ double[] dataReal = new double[dataFFT.length/2];
+ double[] dataImaginary = new double[dataFFT.length/2];
+ for (int i = 0; i < dataReal.length; i++)
+ {
+ dataReal[i] = dataFFT[2*i];
+ dataImaginary[i] = dataFFT[2*i +1];
+ }
+ IcyBufferedImage realImage = new IcyBufferedImage(width, height, dataReal);
+ IcyBufferedImage complexImage = new IcyBufferedImage(width, height, dataImaginary);
+ resSeq.addImage(0, realImage);
+ resSeq.addImage(1, complexImage);
+ Icy.getMainInterface().addSequence(resSeq);
+
+ double[] reverse = new double[dataFFT.length];
+ System.arraycopy(dataFFT, 0, reverse, 0, reverse.length);
+ DoubleFFT_2D fftRec = new DoubleFFT_2D(height, width);
+ fftRec.complexInverse(reverse, true);
+ double[] revImage = new double[width*height];
+ for (int i = 0; i < revImage.length; i++)
+ revImage[i] = reverse[2*i];
+ Sequence reconstructSeq = new Sequence();
+ reconstructSeq.setName("reconstruction");
+ IcyBufferedImage inverseImage = new IcyBufferedImage(width, height, revImage);
+ reconstructSeq.addImage(0, inverseImage);
+ Icy.getMainInterface().addSequence(reconstructSeq);
+ }
+}
diff --git a/src/main/java/plugins/nchenouard/isotropicwavelets/IsotropicWaveletType.java b/src/main/java/plugins/nchenouard/isotropicwavelets/IsotropicWaveletType.java
new file mode 100644
index 0000000000000000000000000000000000000000..be21c5b1144b4d567b34a0bbc4d64d809d68dbf9
--- /dev/null
+++ b/src/main/java/plugins/nchenouard/isotropicwavelets/IsotropicWaveletType.java
@@ -0,0 +1,18 @@
+package plugins.nchenouard.isotropicwavelets;
+
+/**
+ * Available radial profiles for the isotropic wavelet functions
+ *
+ *
+ * @author Nicolas Chenouard (nicolas.chenouard.dev@gmail.com)
+ * @version 1.0
+ * 2014-05-22
+ * gpl v3.0
+ */
+
+public enum IsotropicWaveletType {
+ Simoncelli,
+ Shannon,
+ Papadakis,
+ Meyer,
+}
diff --git a/src/main/java/plugins/nchenouard/isotropicwavelets/IsotropicWavelets.java b/src/main/java/plugins/nchenouard/isotropicwavelets/IsotropicWavelets.java
new file mode 100644
index 0000000000000000000000000000000000000000..a98c4bf8a808ec0c0c138f914cb99f594a4f63fe
--- /dev/null
+++ b/src/main/java/plugins/nchenouard/isotropicwavelets/IsotropicWavelets.java
@@ -0,0 +1,432 @@
+package plugins.nchenouard.isotropicwavelets;
+
+import icy.gui.frame.IcyFrame;
+import icy.gui.frame.progress.AnnounceFrame;
+import icy.image.IcyBufferedImage;
+import icy.main.Icy;
+import icy.plugin.abstract_.PluginActionable;
+import icy.sequence.Sequence;
+import icy.swimmingPool.SwimmingObject;
+import icy.type.collection.array.ArrayUtil;
+
+import java.awt.BorderLayout;
+import java.awt.Dimension;
+import java.awt.GridBagConstraints;
+import java.awt.GridBagLayout;
+import java.awt.GridLayout;
+import java.awt.Insets;
+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.HashMap;
+
+import javax.swing.JButton;
+import javax.swing.JCheckBox;
+import javax.swing.JPanel;
+import javax.swing.JTabbedPane;
+import javax.swing.SwingUtilities;
+
+
+/**
+ *
+ * Isotropic wavelet transform plugin for 2D images
+ *
+ * The plugin provides tools for 2D wavelet decomposition and reconstruction of 2D grayscale images.
+ * Transforms are as decribed in:
+ *
+ * Chenouard, N.; Unser, M., "3D Steerable Wavelets in Practice"
+ * IEEE Transactions on Image Processing, vol.21, no.11, pp.4522,4533, Nov. 2012
+ * doi: 10.1109/TIP.2012.2206044 *
+ *
+ * except being implemented in 2D instead of 3D.
+ *
+ * Please cite the above reference in scientific communications upon use of these tools.
+ *
+ *
+ * @author Nicolas Chenouard (nicolas.chenouard.dev@gmail.com)
+ * @version 1.0
+ * 2014-05-22
+ * gpl v3.0
+ */
+
+/**
+ * Wavelet tranform GUI
+ * */
+public class IsotropicWavelets extends PluginActionable
+{
+ IcyFrame mainFrame;
+
+ JTabbedPane mainPanel;
+ WaveletConfigPanel optionPanel;
+ JButton processSequenceButton;
+ JButton processImageButton;
+ CoefficientSelectionPanel waveletResultsPanel;
+ JButton displayCoefficientsButton;
+ JButton reconstructButton;
+ JCheckBox displayCoefficientsBox;
+
+ Thread runningThread;
+ AnnounceFrame announceFrame;
+
+ @Override
+ public void run()
+ {
+ generateGUI();
+ }
+
+ /**
+ * Initialize the GUI of the plugin
+ * */
+ private void generateGUI()
+ {
+ mainFrame = new IcyFrame("Isotropic wavelet transform", true, true, false, true);
+
+ JTabbedPane mainPanel = new JTabbedPane();
+ mainFrame.setContentPane(mainPanel);
+
+ JPanel analysisPanel = new JPanel(new BorderLayout());
+ mainPanel.addTab("Image Analysis", analysisPanel);
+
+ optionPanel = new WaveletConfigPanel();
+ analysisPanel.add(optionPanel, BorderLayout.CENTER);
+
+ GridBagConstraints c = new GridBagConstraints();
+ c.gridheight = 1;
+ c.gridwidth = 1;
+ c.gridx = 0;
+ c.gridy = 0;
+ c.fill = GridBagConstraints.HORIZONTAL;
+ c.weightx = 1;
+ c.weighty = 1;
+ c.insets = new Insets(2, 2, 2, 2);
+
+ JPanel actionPanel = new JPanel(new GridLayout(3, 1));
+ analysisPanel.add(actionPanel, BorderLayout.SOUTH);
+
+ displayCoefficientsBox = new JCheckBox("Display coefficients", true);
+ displayCoefficientsBox.setSelected(true);
+ actionPanel.add(displayCoefficientsBox);
+ c.gridy++;
+
+ processImageButton = new JButton("Process the focused image");
+ processImageButton.addActionListener(new ActionListener(){
+ @Override
+ public void actionPerformed(ActionEvent e) {
+ processImage();
+ }});
+ actionPanel.add(processImageButton);
+
+ processSequenceButton = new JButton("Process the whole sequence");
+ actionPanel.add(processSequenceButton);
+ processSequenceButton.addActionListener(new ActionListener(){
+ @Override
+ public void actionPerformed(ActionEvent e) {
+ processSequence();
+ }});
+ c.gridy++;
+
+
+ JPanel reconstructionPanel = new JPanel(new BorderLayout());
+ mainPanel.addTab("Image reconstruction", reconstructionPanel);
+
+ waveletResultsPanel = new CoefficientSelectionPanel();
+ reconstructionPanel.add(waveletResultsPanel, BorderLayout.CENTER);
+ waveletResultsPanel.updateResultsTable();
+ c.gridy++;
+
+ JPanel coefficientsActionPanel = new JPanel(new GridBagLayout());
+ reconstructionPanel.add(coefficientsActionPanel, BorderLayout.SOUTH);
+ c.gridy = 0;
+
+ displayCoefficientsButton = new JButton("Display coefficients");
+ displayCoefficientsButton.addActionListener(new ActionListener() {
+ @Override
+ public void actionPerformed(ActionEvent arg0) {
+ displayCurrentCoefficients();
+ }
+ });
+ coefficientsActionPanel.add(displayCoefficientsButton, c);
+ c.gridy++;
+
+ reconstructButton = new JButton("Reconstruct sequence from coefficients");
+ reconstructButton.addActionListener(new ActionListener() {
+ @Override
+ public void actionPerformed(ActionEvent e)
+ {
+ reconstruct();
+ }
+ });
+ coefficientsActionPanel.add(reconstructButton, c);
+
+ mainFrame.setPreferredSize(new Dimension(400, 300));
+ mainFrame.pack();
+ addIcyFrame(mainFrame);
+ mainFrame.setVisible(true);
+ }
+
+ /**
+ * Display selected wavelet coefficients in the SwimmingPool as ICY images.
+ * */
+ private void displayCurrentCoefficients()
+ {
+
+ Object selectedItem = waveletResultsPanel.resultsBox.getSelectedItem();
+ if (selectedItem != null)
+ {
+ final SequenceAnalysisResults results = (SequenceAnalysisResults) ((SwimmingObject) selectedItem).getObject();
+ computationStarted();
+
+ runningThread = new Thread()
+ {
+ @Override
+ public void run()
+ {
+ displayCoefficients(results);
+ SwingUtilities.invokeLater(new Runnable() {
+ @Override
+ public void run() {computationEnded();}
+ });
+ }
+ };
+ runningThread.start();
+ }
+ }
+
+ /**
+ * Reconstruct 2D images from selected wavelet coefficients in the SwimmingPool
+ * */
+ private void reconstruct()
+ {
+ final Object selectedItem = waveletResultsPanel.resultsBox.getSelectedItem();
+ if (selectedItem != null)
+ {
+ final SequenceAnalysisResults results = (SequenceAnalysisResults) ((SwimmingObject) selectedItem).getObject();
+
+ computationStarted();
+
+ runningThread = new Thread()
+ {
+ @Override
+ public void run()
+ {
+ // sort results
+ ArrayList<WaveletAnalysisResults> coefficients = results.getAllResults();
+ ArrayList<Integer> timeList = results.getAllAnalyzedTimesResults();
+ final HashMap<WaveletAnalysisResults, Integer> mapResTime = new HashMap<WaveletAnalysisResults, Integer>();
+ for (int cnt = 0; cnt<coefficients.size(); cnt++)
+ mapResTime.put(coefficients.get(cnt), timeList.get(cnt));
+ // create sequences
+ Sequence reconstructionSequence = new Sequence();
+ reconstructionSequence.setName("Reconstruction");
+ Icy.getMainInterface().addSequence(reconstructionSequence);
+
+ Sequence cropppedReconstructionSequence = null;
+ if (coefficients.get(0).getPadX() > 0 || coefficients.get(0).getPadX() > 0)
+ {
+ cropppedReconstructionSequence = new Sequence();
+ cropppedReconstructionSequence.setName("Reconstruction (cropped)");
+ Icy.getMainInterface().addSequence(cropppedReconstructionSequence);
+ }
+ for (int t = 0; t < coefficients.size(); t++)
+ {
+ WaveletAnalysisResults result = coefficients.get(t);
+ WaveletFilterSet filterSet = result.getWaveletFilters();
+ double[] reconstruction = IsotropicWaveletTransform.isotropicBandlimitedSynthesis(result);
+ reconstructionSequence.addImage(t, new IcyBufferedImage(filterSet.getWidth(), filterSet.getHeight(), reconstruction));
+ if (cropppedReconstructionSequence != null)
+ {
+ cropppedReconstructionSequence.addImage(t, new IcyBufferedImage(result.getFullWidth(), result.getFullHeight(), IsotropicWaveletTransform.unpadImage(reconstruction, result.getFullWidth(), result.getFullHeight(), filterSet.getWidth(), filterSet.getHeight(), result.getPadX(), result.getPadY())));
+ }
+ }
+ SwingUtilities.invokeLater(new Runnable() {
+ @Override
+ public void run() {computationEnded();}
+ });
+ }
+ };
+ runningThread.start();
+ }
+ }
+
+ /**
+ * Compute wavelet coefficients for all the frames of the active ICY sequence
+ * */
+ private void processSequence()
+ {
+ final Sequence seq = getActiveSequence();
+ if ( seq != null)
+ {
+ final int width = seq.getSizeX();
+ final int height = seq.getSizeY();
+ final int numScales = optionPanel.numScaleModel.getNumber().intValue();
+ final boolean prefilter = optionPanel.filterHighPassResidual.isSelected();
+ final IsotropicWaveletType waveletType = (IsotropicWaveletType) optionPanel.waveletTypeBox.getSelectedItem();
+ final boolean isotropicPadding = optionPanel.padForIsotropyBox.isSelected();
+
+ computationStarted();
+
+ runningThread = new Thread()
+ {
+ @Override
+ public void run()
+ {
+ WaveletFilterSet filterSet = new WaveletFilterSet(waveletType, numScales, prefilter, width, height, isotropicPadding);
+ SequenceAnalysisResults results = new SequenceAnalysisResults(seq.getName());
+ for (int t = 0; t < seq.getSizeT(); t++)
+ {
+ IcyBufferedImage image = seq.getImage(t, 0);
+ double[] data_t = (double[])ArrayUtil.arrayToDoubleArray(image.getDataXY(0), image.isSignedDataType() );
+ WaveletAnalysisResults result = IsotropicWaveletTransform.isotropicBandlimitedAnalysis(data_t, seq.getSizeX(), seq.getSizeY(), filterSet);
+ results.setResult(t, result);
+ }
+ Icy.getMainInterface().getSwimmingPool().add(new SwimmingObject(results, "Wavelet coefficients for "+results.getSequenceName()));
+
+ if (displayCoefficientsBox.isSelected())
+ {
+ displayCoefficients(results);
+ }
+
+ SwingUtilities.invokeLater(new Runnable() {
+ @Override
+ public void run() {computationEnded();}
+ });
+ }
+ };
+ runningThread.start();
+ }
+
+ }
+
+ /**
+ * Disable the GUI during coefficient or image computation
+ * */
+ private void computationStarted() {
+ processImageButton.setEnabled(false);
+ processSequenceButton.setEnabled(false);
+ displayCoefficientsButton.setEnabled(false);
+ reconstructButton.setEnabled(false);
+ announceFrame = new AnnounceFrame("Computation started");
+ }
+
+ /**
+ * Re-enable the GUI after mathematical computation has ended
+ * */
+ private void computationEnded() {
+ processImageButton.setEnabled(true);
+ processSequenceButton.setEnabled(true);
+ displayCoefficientsButton.setEnabled(true);
+ reconstructButton.setEnabled(true);
+ announceFrame.close();
+ }
+
+ /**
+ * Compute wavelet coefficients for the active image only
+ * */
+ private void processImage()
+ {
+ final Sequence seq = getActiveSequence();
+ if ( seq != null)
+ {
+ final IcyBufferedImage image = seq.getFirstViewer().getCurrentImage();
+ final int t = seq.getFirstViewer().getPositionT();
+
+ final int width = seq.getSizeX();
+ final int height = seq.getSizeY();
+ final int numScales = optionPanel.numScaleModel.getNumber().intValue();
+ final boolean prefilter = optionPanel.filterHighPassResidual.isSelected();
+ final IsotropicWaveletType waveletType = (IsotropicWaveletType) optionPanel.waveletTypeBox.getSelectedItem();
+ final boolean isotropicPadding = optionPanel.padForIsotropyBox.isSelected();
+
+ computationStarted();
+
+ runningThread = new Thread()
+ {
+ @Override
+ public void run()
+ {
+ WaveletFilterSet filterSet = new WaveletFilterSet(waveletType, numScales, prefilter, width, height, isotropicPadding);
+
+ double[] data_t = (double[])ArrayUtil.arrayToDoubleArray(image.getDataXY(0), image.isSignedDataType() );
+ WaveletAnalysisResults result = IsotropicWaveletTransform.isotropicBandlimitedAnalysis(data_t, seq.getSizeX(), seq.getSizeY(), filterSet);
+ SequenceAnalysisResults results = new SequenceAnalysisResults(seq.getName());
+ results.setResult(t, result);
+ Icy.getMainInterface().getSwimmingPool().add(new SwimmingObject(results, "Wavelet coefficients for frame of "+results.getSequenceName()));
+
+ if (displayCoefficientsBox.isSelected())
+ {
+ displayCoefficients(results);
+ }
+ computationEnded();
+ }
+ };
+ runningThread.start();
+ }
+ }
+
+ /**
+ * Display some wavelet coefficients as ICY images.
+ * @param results the wavelet coefficients to display
+ * */
+ public static void displayCoefficients(SequenceAnalysisResults results)
+ {
+ ArrayList<WaveletAnalysisResults> coefficients = results.getAllResults();
+ ArrayList<Integer> timeList = results.getAllAnalyzedTimesResults();
+
+ if (coefficients.isEmpty())
+ return;
+ int numScales = coefficients.get(0).getWaveletFilters().getNumScales();
+
+ // sort results
+ final HashMap<WaveletAnalysisResults, Integer> mapResTime = new HashMap<WaveletAnalysisResults, Integer>();
+ for (int cnt = 0; cnt<coefficients.size(); cnt++)
+ mapResTime.put(coefficients.get(cnt), timeList.get(cnt));
+ Collections.sort(coefficients, new Comparator<WaveletAnalysisResults>() {
+ @Override
+ public int compare(WaveletAnalysisResults o1,
+ WaveletAnalysisResults o2) {
+ return mapResTime.get(o1).compareTo(mapResTime.get(o2));
+ }
+ });
+
+ // create output sequences
+ ArrayList<Sequence> waveletBandsSequenceList = new ArrayList<Sequence>(numScales);
+ for (int scale = 0; scale < numScales; scale++)
+ {
+ Sequence scaleSeq = new Sequence();
+ scaleSeq.setName("scale "+ scale);
+ waveletBandsSequenceList.add(scaleSeq);
+ Icy.getMainInterface().addSequence(scaleSeq);
+ }
+ Sequence lpSequence = new Sequence();
+ lpSequence.setName("Low frequency residual");
+ Icy.getMainInterface().addSequence(lpSequence);
+
+ Sequence hpSequence = null;
+ if (coefficients.get(0).getHPResidual() != null)
+ {
+ hpSequence = new Sequence();
+ hpSequence.setName("High frequency residual");
+ Icy.getMainInterface().addSequence(hpSequence);
+ }
+ for (int t = 0; t < coefficients.size(); t++)
+ {
+ WaveletAnalysisResults result = coefficients.get(t);
+ WaveletFilterSet filterSet = result.getWaveletFilters();
+ for (int scale = 0; scale < numScales; scale++)
+ {
+ IcyBufferedImage scaleImage = new IcyBufferedImage(filterSet.getScaleWidth(scale), result.getWaveletFilters().getScaleHeight(scale), result.getWaveletBand(scale));
+ Sequence scaleSeq = waveletBandsSequenceList.get(scale);
+ scaleSeq.addImage(t, scaleImage);
+ }
+ IcyBufferedImage lpImage = new IcyBufferedImage(filterSet.getLPWidth(), filterSet.getLPHeight(), result.getLPResidual());
+ lpSequence.addImage(t, lpImage);
+ if (hpSequence != null)
+ {
+ IcyBufferedImage hpImage = new IcyBufferedImage(filterSet.getHPWidth(), filterSet.getHPHeight(), result.getHPResidual());
+ hpSequence.addImage(t, hpImage);
+ }
+ }
+ }
+}
diff --git a/src/main/java/plugins/nchenouard/isotropicwavelets/SequenceAnalysisResults.java b/src/main/java/plugins/nchenouard/isotropicwavelets/SequenceAnalysisResults.java
new file mode 100644
index 0000000000000000000000000000000000000000..ce3cb08085bf577f22463ae40ea80872f4501197
--- /dev/null
+++ b/src/main/java/plugins/nchenouard/isotropicwavelets/SequenceAnalysisResults.java
@@ -0,0 +1,113 @@
+package plugins.nchenouard.isotropicwavelets;
+
+import java.util.ArrayList;
+
+/**
+ * Set of wavelet coefficients for a temporal sequence of 2D images
+ *
+ * @author Nicolas Chenouard (nicolas.chenouard.dev@gmail.com)
+ * @version 1.0
+ * 2014-05-22
+ * gpl v3.0
+ */
+
+
+public class SequenceAnalysisResults
+{
+ protected String sequenceName = "";// name of the analyzed sequence
+ protected ArrayList<Integer> analyzedTimes = new ArrayList<Integer>();// set of frames that have been analyzed with the wavelet transform
+ protected ArrayList<WaveletAnalysisResults> results = new ArrayList<WaveletAnalysisResults>();// set of wavelet coefficients corresponding to the analyzed frames
+
+ /**
+ * Initialize the results objects using a defined sequence name
+ * @param sequenceName the name of the analyzed sequence
+ * */
+ public SequenceAnalysisResults(String sequenceName)
+ {
+ this.sequenceName = sequenceName;
+ }
+ /**
+ * Get the name of the analyzed sequence
+ * @return name of the analyzed sequence
+ * */
+ public String getSequenceName()
+ {
+ return sequenceName;
+ }
+
+ /**
+ * Get the set of wavelet coefficients for all the analyzed frames with the wavelet transform
+ * @return array of wavelet coefficients corresponding to analyzed frames
+ * */
+ public ArrayList<WaveletAnalysisResults> getAllResults()
+ {
+ return new ArrayList<WaveletAnalysisResults>(results);
+ }
+
+ /**
+ * Get the time index of analyzed frames with the wavelet transform
+ * @return array of time index of analyzed frames
+ * */
+ public ArrayList<Integer> getAllAnalyzedTimesResults()
+ {
+ return new ArrayList<Integer>(analyzedTimes);
+ }
+
+ /**
+ * Store a set of wavelet coefficients for a determined time frame
+ * @param t the time index of the analyzed 2D image in the sequence
+ * @param result a set of wavelet coefficients
+ * */
+ public void setResult(int t, WaveletAnalysisResults result)
+ {
+ int idx = -1;
+ int cnt = 0;
+ for (int t2:analyzedTimes)
+ {
+ if (t2 == t)
+ {
+ idx = cnt;
+ break;
+ }
+ cnt++;
+ }
+ if (idx == -1)
+ {
+ analyzedTimes.add(t);
+ results.add(result);
+ }
+ else
+ {
+ results.set(idx, result);
+ }
+ }
+
+ /**
+ * Get the computed wavelet coefficients for a given time index in the sequence of images
+ * @param t time index
+ * @return a set of wavelet coefficient for the frame at time t
+ * */
+ public WaveletAnalysisResults getResult(int t)
+ {
+ int cnt = 0;
+ for (int t2:analyzedTimes)
+ {
+ if (t2 == t)
+ return results.get(cnt);
+ cnt++;
+ }
+ return null;
+ }
+ /**
+ * Reset the coefficient storage structure with a given set of wavelet coefficients
+ * @param results the set of wavelet coefficients to reset the storage structure with. Frames are assumed to start at 0 and end at results.size()-1
+ * */
+ public void setAllResults(ArrayList<WaveletAnalysisResults> results)
+ {
+ this.results.clear();
+ this.analyzedTimes.clear();
+ this.results.addAll(results);
+ for (int t = 0; t < results.size(); t++)
+ this.analyzedTimes.add(t);
+ }
+}
diff --git a/src/main/java/plugins/nchenouard/isotropicwavelets/WaveletAnalysisResults.java b/src/main/java/plugins/nchenouard/isotropicwavelets/WaveletAnalysisResults.java
new file mode 100644
index 0000000000000000000000000000000000000000..2f46f27593dfe2b3bd873731e85e4d88315c8e5e
--- /dev/null
+++ b/src/main/java/plugins/nchenouard/isotropicwavelets/WaveletAnalysisResults.java
@@ -0,0 +1,141 @@
+package plugins.nchenouard.isotropicwavelets;
+
+/**
+ *
+ * Coefficients and parameters corresponding to the wavelet decomposition of an 2D image*
+ *
+ * 2D wavelet decomposition and reconstruction algorithms are provided.
+ * Transforms are as decribed in:
+ *
+ * Chenouard, N.; Unser, M., "3D Steerable Wavelets in Practice"
+ * IEEE Transactions on Image Processing, vol.21, no.11, pp.4522,4533, Nov. 2012
+ * doi: 10.1109/TIP.2012.2206044 *
+ *
+ * except being implemented in 2D instead of 3D.
+ *
+ * Please cite the above reference in scientific communications upon use of these tools.
+ *
+ *
+ * @author Nicolas Chenouard (nicolas.chenouard.dev@gmail.com)
+ * @version 1.0
+ * 2014-05-22
+ * gpl v3.0
+ */
+
+public class WaveletAnalysisResults
+{
+ double[] hpResidual; // high frequency residual that cannot analyzed with an isotropic filter (i.e. 'corners' of the Fourier domain)
+ double[] lpResidual; // low frequency residual (i.e., center disk in the Fourier domain)
+ double[][] bands; // Wavelet coefficients. One double array per scale.
+ WaveletFilterSet filterSet; // set of filters used for image decomposition
+ int padX = 0; // left zero-padding of the image before wavelet decomposition
+ int padY = 0; // top zero-padding of the image before wavelet decomposition
+ int fullWidth; // width of the image after zero-padding
+ int fullHeight; // height of the image after zero-padding
+
+ boolean storeCoefficientsInFourierDomain = false; // true if coefficients are stored in the Fourier domain
+
+ /**
+ * Empty constructor
+ * */
+ protected WaveletAnalysisResults(){}
+
+ /**
+ * Initialize the wavelet analysis result object.
+ * @param bands wavelet coefficients. One array of double per scale.
+ * @param lpResidual low frequency residual
+ * @param hpResidual high frequency residual
+ * @param filterSet set of filters used for image decomposition
+ * @param coefficientsInFourierDomain true if coefficients are stored in the Fourier domain
+ *
+ * */
+ public WaveletAnalysisResults(double[][] bands, double[] lpResidual, double[] hpResidual, WaveletFilterSet filterSet, boolean coefficientsInFourierDomain)
+ {
+ this.bands = bands;
+ this.lpResidual = lpResidual;
+ this.hpResidual = hpResidual;
+ this.filterSet = filterSet;
+ this.storeCoefficientsInFourierDomain = coefficientsInFourierDomain;
+ }
+
+ /**
+ * Get the wavelet filters used for image decomposition
+ * @return the set of wavelet filters used for image decomposition
+ * */
+ public WaveletFilterSet getWaveletFilters()
+ {
+ return filterSet;
+ }
+
+ /**
+ * Get the set of wavelet coefficients at a given scale as a 1D array of double
+ * @param scale int
+ * @return wavelet coefficients
+ * */
+ public double[] getWaveletBand(int scale)
+ {
+ return bands[scale];
+ }
+
+ /**
+ * Get the high frequency residual as a 1D array
+ * @return high frequency residual
+ * */
+ public double[] getHPResidual() {
+ return hpResidual;
+ }
+
+ /**
+ * Get the low frequency residual as a 1D array
+ * @return high frequency residual
+ * */
+ public double[] getLPResidual() {
+ return lpResidual;
+ }
+
+ /**
+ * Get the number of scales used for wavelet decomposition
+ * @return number of wavelet scales
+ * */
+ public int getNumScales()
+ {
+ if (filterSet != null)
+ return filterSet.getNumScales();
+ else
+ return bands.length;
+ }
+
+ /**
+ * Get the padding size at the left of the image
+ * @return x direction padding size
+ * */
+ public int getPadX()
+ {
+ return padX;
+ }
+
+ /**
+ * Get the padding size at the top of the image
+ * @return y direction padding size
+ * */
+ public int getPadY()
+ {
+ return padY;
+ }
+
+ /**
+ * Get the width of the analyzed image after padding
+ * @return width after padding
+ * */
+ public int getFullWidth() {
+ return fullWidth;
+ }
+
+ /**
+ * Get the height of the analyzed image after padding
+ * @return height after padding
+ * */
+ public int getFullHeight() {
+ return fullHeight;
+ }
+}
\ No newline at end of file
diff --git a/src/main/java/plugins/nchenouard/isotropicwavelets/WaveletConfigPanel.java b/src/main/java/plugins/nchenouard/isotropicwavelets/WaveletConfigPanel.java
new file mode 100644
index 0000000000000000000000000000000000000000..fd7c701c41162db24c29321cf3ae7330b2da4538
--- /dev/null
+++ b/src/main/java/plugins/nchenouard/isotropicwavelets/WaveletConfigPanel.java
@@ -0,0 +1,105 @@
+package plugins.nchenouard.isotropicwavelets;
+
+import java.awt.GridBagConstraints;
+import java.awt.GridBagLayout;
+import java.awt.Insets;
+
+import javax.swing.JCheckBox;
+import javax.swing.JComboBox;
+import javax.swing.JLabel;
+import javax.swing.JPanel;
+import javax.swing.JSpinner;
+import javax.swing.SpinnerNumberModel;
+
+public class WaveletConfigPanel extends JPanel
+{
+ /**
+ * GUI for manually setting wavelet decomposition parameters
+ *
+ * @author Nicolas Chenouard (nicolas.chenouard.dev@gmail.com)
+ * @version 1.0
+ * 2014-05-22
+ * gpl v3.0
+ */
+
+ private static final long serialVersionUID = -7253036809006307841L;
+ JComboBox<IsotropicWaveletType> waveletTypeBox;
+ JSpinner numScaleSpinner;
+ SpinnerNumberModel numScaleModel;
+ JCheckBox filterHighPassResidual;
+ JCheckBox padForIsotropyBox;
+
+ public WaveletConfigPanel()
+ {
+ this.setLayout(new GridBagLayout());
+
+ GridBagConstraints c = new GridBagConstraints();
+ c.gridheight = 1;
+ c.gridwidth = 1;
+ c.gridx = 0;
+ c.gridy = 0;
+ c.fill = GridBagConstraints.HORIZONTAL;
+ c.weightx = 1;
+ c.weighty = 1;
+ c.insets = new Insets(2, 2, 2, 2);
+
+ this.add(new JLabel("Number of scales:"), c);
+ c.gridy++;
+
+ numScaleModel = new SpinnerNumberModel(3, 1, 100, 1);
+ numScaleSpinner = new JSpinner(numScaleModel);
+ this.add(numScaleSpinner, c);
+ c.gridy++;
+
+ this.add(new JLabel("Wavelet profile:"), c);
+ c.gridy++;
+
+ waveletTypeBox = new JComboBox<IsotropicWaveletType>(IsotropicWaveletType.values());
+ this.add(waveletTypeBox, c);
+ c.gridy++;
+
+ filterHighPassResidual = new JCheckBox("Filter out high frequencies", true);
+ this.add(filterHighPassResidual, c);
+ c.gridy++;
+
+ padForIsotropyBox = new JCheckBox("Isotropic zero-padding of image", false);
+ this.add(padForIsotropyBox, c);
+ c.gridy++;
+ }
+
+ /**
+ * Get the number of scales to use for wavelet decomposition
+ * @return number of wavelet scales
+ * */
+ public int getNumScales()
+ {
+ return numScaleModel.getNumber().intValue();
+ }
+
+ /**
+ * Determine whether a prefilter needs to be used to eliminate the high frequency residual
+ * @return true if prefiltering of high frequencies has to be performed
+ * */
+ public boolean isPrefilter()
+ {
+ return filterHighPassResidual.isSelected();
+ }
+
+ /**
+ * Get the type of radial function to use for wavelet decomposition
+ * @return type of radial wavelet function
+ * */
+ public IsotropicWaveletType getWaveletType()
+ {
+ return (IsotropicWaveletType) waveletTypeBox.getSelectedItem();
+ }
+
+ /**
+ * Determine whether zero-padding is to be used to make the analyzed image isotropic
+ * @return true if zero padding is required
+ * */
+ public boolean isIsotropicPadding()
+ {
+ return padForIsotropyBox.isSelected();
+ }
+}
\ No newline at end of file
diff --git a/src/main/java/plugins/nchenouard/isotropicwavelets/WaveletFilterSet.java b/src/main/java/plugins/nchenouard/isotropicwavelets/WaveletFilterSet.java
new file mode 100644
index 0000000000000000000000000000000000000000..d3ff48da1f7ef4ce4ed65738267dfeb501ff1ee7
--- /dev/null
+++ b/src/main/java/plugins/nchenouard/isotropicwavelets/WaveletFilterSet.java
@@ -0,0 +1,244 @@
+package plugins.nchenouard.isotropicwavelets;
+
+/**
+ *
+ * Set of wavelet filters for image decomposition and reconstruction
+ *
+ * 2D wavelet decomposition and reconstruction algorithms are provided.
+ * Transforms are as decribed in:
+ *
+ * Chenouard, N.; Unser, M., "3D Steerable Wavelets in Practice"
+ * IEEE Transactions on Image Processing, vol.21, no.11, pp.4522,4533, Nov. 2012
+ * doi: 10.1109/TIP.2012.2206044
+ *
+ * except being implemented in 2D instead of 3D.
+ *
+ * Please cite the above reference in scientific communications upon use of these tools.
+ *
+ *
+ * @author Nicolas Chenouard (nicolas.chenouard.dev@gmail.com)
+ * @version 1.0
+ * 2014-05-22
+ * gpl v3.0
+ */
+
+import java.util.ArrayList;
+
+public class WaveletFilterSet {
+ ArrayList<double[][]> waveletFiltersFFT;// set of wavelet filters in the Fourier domain. One filter as a double[][] per scale.
+ double[][] prefilterFFT; // set of filters in the Fourier domain for high frequency band filtering out
+ int numScales; // number of scales for wavelet decomposition
+ IsotropicWaveletType waveletType; // type of radial profile of the wavelet function
+ boolean prefilter; // true if high frequencies (corners of the Fourier domain) need to be filtered out
+ int width; // width of the 2D images to process
+ int height; // height of the 2D images to process
+ int[][] scaleDims; // dimensions of the wavelet bands
+ int[] lpSize; // dimensions of the low frequency residual
+ int[] hpSize; // dimensions of the high frequency residual
+ boolean isotropicPadding; // true if zero-padding is used to obtain an isotropic image
+
+ /**
+ * Initialize the set of filters without using zero-padding
+ *
+ * @param waveletType type of radial profile of the wavelet function
+ * @param numScales number of scales for wavelet decomposition
+ * @param prefilter true if high frequencies (corners of the Fourier domain) need to be filtered out
+ * @param width width of the 2D images to process
+ * @param height height of the 2D images to process
+ * */
+ public WaveletFilterSet(IsotropicWaveletType waveletType, int numScales, boolean prefilter, int width, int height)
+ {
+ this(waveletType, numScales, prefilter, width, height, false);
+ }
+
+ /**
+ * Initialize the set of filters
+ *
+ * @param waveletType type of radial profile of the wavelet function
+ * @param numScales number of scales for wavelet decomposition
+ * @param prefilter true if high frequencies (corners of the Fourier domain) need to be filtered out
+ * @param width width of the 2D images to process
+ * @param height height of the 2D images to process
+ * @param isotropicPadding true if zero-padding is used to obtain an isotropic image
+ * */
+ public WaveletFilterSet(IsotropicWaveletType waveletType, int numScales, boolean prefilter, int width, int height, boolean isotropicPadding)
+ {
+ this.waveletType = waveletType;
+ this.numScales = numScales;
+ this.isotropicPadding = isotropicPadding;
+ // check image dimensions
+ if (areFeasibleDimensions(width, height, numScales))
+ {
+ this.width = width;
+ this.height = height;
+ }
+ else
+ {
+ this.width = getFeasibleSize(width, numScales);
+ this.height = getFeasibleSize(height, numScales);
+ if (isotropicPadding)
+ {
+ this.width = Math.max(this.width, this.height);
+ this.height = this.width;
+ }
+ }
+ this.prefilter = prefilter;
+ if (prefilter)
+ this.prefilterFFT = WaveletFunction.getPreFiltersFFT(waveletType, this.width, this.height);
+ hpSize = new int[]{this.width, this.height};
+ waveletFiltersFFT = new ArrayList<double[][]>(numScales);
+ int scaleWidth = this.width;
+ int scaleHeight = this.height;
+ scaleDims = new int[numScales][2];
+ for (int i = 0; i < numScales; i++)
+ {
+ double[][] filters = WaveletFunction.getFilterFFT(waveletType, scaleWidth, scaleHeight);
+ waveletFiltersFFT.add(filters);
+ scaleDims[i][0] = scaleWidth;
+ scaleDims[i][1] = scaleHeight;
+ scaleWidth = (int)scaleWidth/2;
+ scaleHeight = (int) scaleHeight/2;
+ }
+ lpSize = new int[]{scaleWidth, scaleHeight};
+ }
+
+ /**
+ * Return the next image size (superior to the original image size) such that dyadic subsampling is feasible for a number of scales
+ * @param size size of the image along a dimension
+ * @param numScales number of wavelet scales
+ * @return the minimum size the image should be to make dyadic subsampling feasible
+ * */
+ private int getFeasibleSize(int size, int numScales) {
+ int feasibleSize = size;
+ if (feasibleSize % 2 != 0)
+ feasibleSize ++;
+ if (Math.floor(feasibleSize/Math.pow(2, numScales)) == feasibleSize/Math.pow(2, numScales))
+ return feasibleSize;
+ int minSize = (int) Math.pow(2, numScales);
+ feasibleSize = (int) Math.ceil(feasibleSize/(double)minSize)*minSize;
+ return feasibleSize;
+ }
+
+ /**
+ * Test whether the width and height of the image comply with the subsampling requirements of wavelet decomposition
+ * @param width width of the image to decompose
+ * @param height height of the image to decompose
+ * @param numScales number of wavelet scales to use
+ * @return true if the height and width of the image comply with the number of wavelet scales, false otherwise
+ * */
+ private boolean areFeasibleDimensions(int width, int height, int numScales) {
+ if (width%2 != 0 || height != 0)
+ return false;
+ if (Math.floor(width/Math.pow(2, numScales)) != (width/Math.pow(2, numScales)) ||Math.floor(height/Math.pow(2, numScales)) != (width/Math.pow(2, numScales)))
+ return false;
+ return true;
+ }
+
+ /**
+ * Get the width of the analyzed image
+ * @return width of the analyzed image
+ * */
+ public int getWidth()
+ {
+ return width;
+ }
+
+ /**
+ * Get the height of the analyzed image
+ * @return height of the analyzed image
+ * */
+ public int getHeight()
+ {
+ return height;
+ }
+
+ /**
+ * Get the number of scales for the wavelet decomposition
+ * @return number of wavelet scales
+ * */
+ public int getNumScales()
+ {
+ return numScales;
+ }
+
+ /**
+ * Determines whether high frequencies are filtered out
+ * @return true if the high frequency component is prefiltered, false else
+ * */
+ public boolean isPrefilter()
+ {
+ return prefilter;
+ }
+
+ /**
+ * Get the width of a given wavelet scale
+ * @param scale scale of the wavelet transform
+ * @return width of the wavelet band
+ * */
+ public int getScaleWidth(int scale)
+ {
+ return scaleDims[scale][0];
+ }
+
+ /**
+ * Get the height of a given wavelet scale
+ * @param scale scale of the wavelet transform
+ * @return height of the wavelet band
+ * */
+ public int getScaleHeight(int scale)
+ {
+ return scaleDims[scale][1];
+ }
+
+ /**
+ * Get the width of the low frequency residual
+ * @return width of the low frequency residual
+ * */
+ public int getLPWidth()
+ {
+ return lpSize[0];
+ }
+
+ /**
+ * Get the height of the low frequency residual
+ * @return height of the low frequency residual
+ * */
+ public int getLPHeight()
+ {
+ return lpSize[1];
+ }
+
+ /**
+ * Get the width of the high frequency residual
+ * @return width of the high frequency residual
+ * */
+ public int getHPWidth()
+ {
+ return hpSize[0];
+ }
+
+ /**
+ * Get the width of the height frequency residual
+ * @return width of the height frequency residual
+ * */
+ public int getHPHeight()
+ {
+ return hpSize[1];
+ }
+
+ /**
+ * get the type of radial profile of the wavelet function
+ * @return type of radial profile
+ * */
+ public IsotropicWaveletType getWaveletType() {
+ return this.waveletType;
+ }
+
+ /**
+ * Determines whether zero-padding is used to obtain an isotropic image
+ * @return true if zero-padding is used to obtain an isotropic image
+ * */
+ public boolean isIsotropicPadding() {
+ return isotropicPadding;
+ }
+}
\ No newline at end of file
diff --git a/src/main/java/plugins/nchenouard/isotropicwavelets/WaveletFunction.java b/src/main/java/plugins/nchenouard/isotropicwavelets/WaveletFunction.java
new file mode 100644
index 0000000000000000000000000000000000000000..0858b7b25a3d0bd3c0002be7c8c8560bf4303305
--- /dev/null
+++ b/src/main/java/plugins/nchenouard/isotropicwavelets/WaveletFunction.java
@@ -0,0 +1,320 @@
+package plugins.nchenouard.isotropicwavelets;
+
+import icy.image.IcyBufferedImage;
+import icy.main.Icy;
+import icy.sequence.Sequence;
+
+/**
+ *
+ * Filters for wavelet transform and low-pass prefiltering
+ *
+ * 2D wavelet decomposition and reconstruction algorithms are provided.
+ * Transforms are as decribed in:
+
+ * Chenouard, N.; Unser, M., "3D Steerable Wavelets in Practice"
+ * IEEE Transactions on Image Processing, vol.21, no.11, pp.4522,4533, Nov. 2012
+ * doi: 10.1109/TIP.2012.2206044 *
+ *
+ * except being implemented in 2D instead of 3D.
+ *
+ * Please cite the above reference in scientific communications upon use of these tools.
+ *
+ *
+ * @author Nicolas Chenouard (nicolas.chenouard.dev@gmail.com)
+ * @version 1.0
+ * 2014-05-22
+ * gpl v3.0
+ */
+
+public class WaveletFunction
+{
+ /**
+ * nu function for Meyer wavelets
+ * @param t
+ * @return nu(t)
+ * */
+ private static double nu(double t)
+ {
+ if (t <= 0 || t >= 1 )
+ return 0;
+ else return t*t*t*t*(35 - 84*t + 70*t*t - 20*t*t*t);
+ }
+
+ /**
+ * Get wavelet filter analysis filters in Fourier domain
+ * @param waveletType type of radial function for the wavelet
+ * @param width width of the filter
+ * @param height height of the filter
+ * @return 2D wavelet filter in the Fourier domain
+ * */
+ public static double[][] getFilterFFT(IsotropicWaveletType waveletType, int width, int height) {
+ double[] distMap = new double[width*height];
+ int cX = (int) Math.ceil((width + 1)/2);
+ int cY = (int) Math.ceil((height + 1)/2);
+ double normX = (double)(cX*cX);
+ double normY = (double)(cY*cY);
+
+ for (int i = 0; i <= cX; i++)
+ {
+ for (int j = 0; j <= cY; j++)
+ {
+ distMap[i + j*width] = (double) Math.sqrt(i*i/normX + j*j/normY);
+ }
+ for (int j = cY + 1; j < height; j++)
+ {
+ distMap[i + j*width] = (double) Math.sqrt(i*i/normX + (j - 2*cY)*(j - 2*cY)/normY);
+ }
+ }
+ for (int i = cX + 1; i < width; i++)
+ {
+ for (int j = 0; j <= cY; j++)
+ {
+ distMap[i + j*width] = (double) Math.sqrt((i - 2*cX)*(i - 2*cX)/normX + j*j/normY);
+ }
+ for (int j = cY + 1; j < height; j++)
+ {
+ distMap[i + j*width] = (double) Math.sqrt((i - 2*cX)*(i - 2*cX)/normX + (j - 2*cY)*(j - 2*cY)/normY);
+ }
+ }
+ double[] lp = new double[distMap.length];
+ double[] hp = new double[distMap.length];
+ switch(waveletType)
+ {
+ case Meyer:
+ {
+ // maskLP = mask4; < pi/4
+ // maskLP = maskLP + cos(0.5*pi*nu(4*dist-1)).*mask2.*(1-mask4);
+ // maskHP = sqrt(1-maskLP.^2);
+ for (int i = 0; i < distMap.length; i++)
+ {
+ double dist = distMap[i];
+ if (dist < 1d/4)
+ {
+ lp[i] = 1;
+ }
+ else if (dist < 1d/2)
+ {
+ lp[i] = Math.cos(0.5*Math.PI*nu(4*dist - 1));
+ }
+ hp[i] = Math.sqrt(1 - lp[i]*lp[i]);
+ }
+ break;
+ }
+ case Papadakis:
+ {
+ // %create mask selecting r<=3*pi/10
+ // maskLP = sqrt((1+cos(5*dist*pi - 3*pi/2))/2).*mask2.*(1-mask310);
+ // maskLP(1) = 0;
+ // maskLP = maskLP + mask310;
+ //
+ // maskHP = sqrt(1-maskLP.^2);
+
+ for (int i = 0; i < distMap.length; i++)
+ {
+ double dist = distMap[i];
+ if (dist <= 3d/10)
+ {
+ lp[i] = 1;
+ }
+ else if (dist <= 1d/2)
+ {
+ lp[i] = Math.sqrt((1 + Math.cos(5*dist*Math.PI - 3*Math.PI/2))/2);
+ hp[i] = Math.sqrt(1 - lp[i]*lp[i]);
+ }
+ else
+ hp[i] = 1;
+ }
+ break;
+ }
+ case Shannon:
+ {
+ for (int i = 0; i < distMap.length; i++)
+ {
+ double dist = distMap[i];
+ if (dist < 0.5f)// ||w|| < pi/2
+ {
+ lp[i] = 1;
+ }
+ else // ||w|| >= pi/2
+ {
+ hp[i] = 1;
+ }
+
+ }
+ hp[0] = 0;
+ // lp[0] = 0;
+ break;
+ }
+ case Simoncelli:
+ {
+ for (int i = 0; i < distMap.length; i++)
+ {
+ double dist = distMap[i];
+ if (dist < 0.5f)// ||w|| < pi/2
+ {
+ if (dist <= 0.25f) // ||w|| <= pi/4
+ {
+ lp[i] = 1;
+ }
+ else // ||w|| > pi/4
+ {
+ hp[i] = (double) Math.cos(Math.PI*0.5*Math.log(2*dist)/Math.log(2));
+ lp[i] = (double) Math.cos(Math.PI*0.5*Math.log(4*dist)/Math.log(2));
+ }
+ }
+ else // ||w|| >= pi/2
+ {
+ hp[i] = 1;
+ }
+
+ }
+ hp[0] = 0;
+ // lp[0] = 0;
+ break;
+ }
+ }
+ boolean debug = false;
+ if (debug)
+ {
+ double[] sumSqFilters = new double[hp.length];
+ for (int i = 0; i < sumSqFilters.length; i++)
+ {
+ sumSqFilters[i] = hp[i]*hp[i] + lp[i]*lp[i];
+ }
+ Sequence distSeq = new Sequence();
+ distSeq.setName("dist");
+ IcyBufferedImage realImage = new IcyBufferedImage(width, height, distMap);
+ distSeq.addImage(0, realImage);
+ IcyBufferedImage lpImage = new IcyBufferedImage(width, height, lp);
+ distSeq.addImage(1, lpImage);
+ IcyBufferedImage hpImage = new IcyBufferedImage(width, height, hp);
+ distSeq.addImage(2, hpImage);
+ IcyBufferedImage sumSqImage = new IcyBufferedImage(width, height, sumSqFilters);
+ distSeq.addImage(3, sumSqImage);
+ Icy.getMainInterface().addSequence(distSeq);
+ }
+ return new double[][]{lp, hp};
+ }
+
+ /**
+ * Get prefilter in Fourier domain for removing high frequency coefficients (corners of Fourier domain)
+ * @param waveletType type of radial function for the wavelet
+ * @param width width of the filter
+ * @param height height of the filter
+ * @return 2D Fourier coefficients of the prefilter
+ * */
+ public static double[][] getPreFiltersFFT(IsotropicWaveletType waveletType, int width, int height)
+ {
+ double[] distMap = new double[width*height];
+ int cX = (int) Math.ceil((width + 1)/2);
+ int cY = (int) Math.ceil((height + 1)/2);
+ double normX = (double)(cX*cX);
+ double normY = (double)(cY*cY);
+
+ for (int i = 0; i <= cX; i++)
+ {
+ for (int j = 0; j <= cY; j++)
+ {
+ distMap[i + j*width] = (double) Math.sqrt(i*i/normX + j*j/normY);
+ }
+ for (int j = cY + 1; j < height; j++)
+ {
+ distMap[i + j*width] = (double) Math.sqrt(i*i/normX + (j - 2*cY)*(j - 2*cY)/normY);
+ }
+
+ }
+ for (int i = cX + 1; i < width; i++)
+ {
+ for (int j = 0; j <= cY; j++)
+ {
+ distMap[i + j*width] = (double) Math.sqrt((i - 2*cX)*(i - 2*cX)/normX + j*j/normY);
+ }
+ for (int j = cY + 1; j < height; j++)
+ {
+ distMap[i + j*width] = (double)Math.sqrt((i - 2*cX)*(i - 2*cX)/normX + (j - 2*cY)*(j - 2*cY)/normY);
+ }
+ }
+ double[] hp = new double[distMap.length];
+ double[] lp = new double[distMap.length];
+ switch(waveletType)
+ {
+ case Meyer:
+ {
+ for (int i = 0; i < distMap.length; i++)
+ {
+ double dist = distMap[i]/2;
+ if (dist < 1d/4)
+ {
+ lp[i] = 1;
+ }
+ else if (dist < 1d/2)
+ {
+ lp[i] = Math.cos(0.5*Math.PI*nu(4*dist - 1));
+ }
+ hp[i] = Math.sqrt(1 - lp[i]*lp[i]);
+ }
+ break;
+ }
+ case Papadakis:
+ {
+ for (int i = 0; i < distMap.length; i++)
+ {
+ double dist = distMap[i]/2;
+ if (dist <= 3d/10)
+ {
+ lp[i] = 1;
+ }
+ else if (dist <= 1d/2)
+ {
+ lp[i] = Math.sqrt((1 + Math.cos(5*dist*Math.PI - 3*Math.PI/2))/2);
+ hp[i] = Math.sqrt(1 - lp[i]*lp[i]);
+ }
+ else
+ hp[i] = 1;
+ }
+ break;
+ }
+ case Shannon:
+ {
+ for (int i = 0; i < distMap.length; i++)
+ {
+ double dist = distMap[i];
+ if (dist > 1f)
+ hp[i] = 1;
+ else
+ lp[i] = 1;
+ }
+ break;
+ }
+ case Simoncelli:
+ {
+ for (int i = 0; i < distMap.length; i++)
+ {
+ double dist = distMap[i];
+ if (dist > 1f)
+ hp[i] = 1;
+ else if (dist > 0.5f)// ||w|| > pi/2
+ {
+ hp[i] = (double) Math.cos(Math.PI*0.5*Math.log(dist)/Math.log(2));
+ lp[i] = (double)Math.sqrt(1 - hp[i]*hp[i]);
+ }
+ else
+ {
+ lp[i] = 1;
+ }
+ }
+ break;
+ }
+ }
+ boolean debug = false;
+ if (debug)
+ {
+ Sequence distSeq = new Sequence();
+ distSeq.setName("dist");
+ IcyBufferedImage hpImage = new IcyBufferedImage(width, height, hp);
+ distSeq.addImage(1, hpImage);
+ Icy.getMainInterface().addSequence(distSeq);
+ }
+ return new double[][]{lp, hp};
+ }
+}