org.scijava.ops.image.create.DefaultCreateKernel2ndDerivBiGauss Maven / Gradle / Ivy
Show all versions of scijava-ops-image Show documentation
/*
* #%L
* Image processing operations for SciJava Ops.
* %%
* Copyright (C) 2014 - 2024 SciJava developers.
* %%
* Redistribution and use in source and binary forms, with or without
* modification, are permitted provided that the following conditions are met:
*
* 1. Redistributions of source code must retain the above copyright notice,
* this list of conditions and the following disclaimer.
* 2. Redistributions in binary form must reproduce the above copyright notice,
* this list of conditions and the following disclaimer in the documentation
* and/or other materials provided with the distribution.
*
* THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
* AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
* IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
* ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR CONTRIBUTORS BE
* LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
* CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
* SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
* INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
* CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
* ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
* POSSIBILITY OF SUCH DAMAGE.
* #L%
*/
package org.scijava.ops.image.create;
import java.util.function.BiFunction;
import org.scijava.ops.image.logic.Default;
import net.imglib2.Cursor;
import net.imglib2.Dimensions;
import net.imglib2.FinalInterval;
import net.imglib2.RandomAccessibleInterval;
import net.imglib2.img.Img;
import net.imglib2.type.Type;
import net.imglib2.type.numeric.ComplexType;
import net.imglib2.view.Views;
/**
* Creates 2nd derivative of an isotropic BiGauss kernel with the pair of sigmas
* specification.
*
* The BiGauss kernel is a composition of two standard Gauss kernels. If we were
* to assume 1D kernel centered at zero (0), an inner kernel of the BiGauss,
* with its shape given with sigmas[0], would span from -sigmas[0] till
* +sigmas[0]; outer kernel, with its shape given with sigmas[1], surrounds the
* inner, e.g. for the positive side, from sigmas[0] till sigmas[0]+2*sigmas[1]
* and its center having at sigmas[0]-sigmas[1]. That is, the inner Gauss exist
* up to its inflection points from which the filter takes shape of the outer
* Gauss. Both kernels are, however, appropriately scaled and shifted to obtain
* a smooth BiGauss kernel.
*
* Note that the kernel is always isotropic. The second parameter gives
* dimensionality of the created kernel.
*
* All values are in units of pixels.
*
* Literature:C. Xiao, M. Staring, Y. Wang, D.P. Shamonin, B.C. Stoel.
* Multiscale Bi-Gaussian Filter for Adjacent Curvilinear Structures Detection
* with Application to Vasculature Images. IEEE TMI, vol. 22, no. 1, 2013.
*
* @author Vladimír Ulman
*/
public final class DefaultCreateKernel2ndDerivBiGauss {
private DefaultCreateKernel2ndDerivBiGauss() {
// Prevent instantiation of static utility class
}
public static , C extends ComplexType>
RandomAccessibleInterval createKernel(final double[] sigmas,
final Integer dimensionality, final C typeVar,
final BiFunction> createImgFunc)
{
// both sigmas must be available
if (sigmas.length < 2) throw new IllegalArgumentException(
"Two sigmas (for inner and outer Gauss)" + " must be supplied.");
// both sigmas must be reasonable
if (sigmas[0] <= 0 || sigmas[1] <= 0) throw new IllegalArgumentException(
"Input sigmas must be both positive.");
// dimension as well...
if (dimensionality <= 0) throw new IllegalArgumentException(
"Input dimensionality must both positive.");
// the size and center of the output image
final long[] dims = new long[dimensionality];
final long[] centre = new long[dimensionality];
// time-saver... (must hold now: dimensionality > 0)
// NB: size of the image is 2px wider than for 0th order BiGauss to have
// some space for smooth approach-to-zero at the kernel image borders
dims[0] = Math.max(3, 2 * (int) (sigmas[0] + 3 * sigmas[1] + 0.5) + 1) + 2;
centre[0] = (int) (dims[0] / 2);
// fill the size and center arrays
for (int d = 1; d < dims.length; d++) {
dims[d] = dims[0];
centre[d] = centre[0];
}
// prepare some scaling constants
/*
* //orig full math version: final double k = (sigmas[1]/sigmas[0]) *
* (sigmas[1]/sigmas[0]); //eq. (6) final double[] C = {
* 1.0/(2.50663*sigmas[0]*sigmas[0]*sigmas[0]),
* 1.0/(2.50663*sigmas[1]*sigmas[1]*sigmas[1]) }; //2.50663 = sqrt(2*PI)
*/
// less math version:
// note that originally there was C[0] for inner Gauss, k*C[1] for outer
// Gauss
// we get rid of k by using new C[0] and C[1]:
final double[] C = { 1.0 / (2.50663 * sigmas[0] * sigmas[0] * sigmas[0]),
1.0 / (2.50663 * sigmas[1] * sigmas[0] * sigmas[0]) };
// prepare squared input sigmas
final double sigmasSq[] = { sigmas[0] * sigmas[0], sigmas[1] * sigmas[1] };
// prepare the output image
final RandomAccessibleInterval out =
(RandomAccessibleInterval) createImgFunc.apply(new FinalInterval(dims),
(T) typeVar);
// fill the output image
final Cursor cursor = Views.iterable(out).cursor();
while (cursor.hasNext()) {
cursor.fwd();
// obtain the current coordinate (use dims to store it)
cursor.localize(dims);
// calculate distance from the image centre
double dist = 0.; // TODO: can JVM reuse this var or is it allocated again
// and again (and
// multipling in the memory)?
for (int d = 0; d < dims.length; d++) {
final double dx = dims[d] - centre[d];
dist += dx * dx;
}
// dist = Math.sqrt(dist); -- gonna work with squared distance
// which of the two Gaussians should we use?
double val = 0.;
if (dist < sigmasSq[0]) {
// the inner one
val = dist / sigmasSq[0] - 1.0;
val *= C[0] * Math.exp(-0.5 * dist / sigmasSq[0]);
}
else {
// the outer one, get new distance first:
dist = Math.sqrt(dist) - (sigmas[0] - sigmas[1]);
dist *= dist;
val = dist / sigmasSq[1] - 1.0;
val *= C[1] * Math.exp(-0.5 * dist / sigmasSq[1]);
}
// compose the real value finally
cursor.get().setReal(val);
}
return out;
}
}