package Algorithms;
import Utils.*;
import java.util.Arrays;
import java.util.concurrent.ExecutionException;
public class Tubeness
{
public static final int IN_PLACE_THRESHOLD = 5;
public static void computeTubenessImage(
ISliceRunner sliceRunner,
int maxWorkers,
byte[] output,
float[] image,
int width,
int height,
int brightestChannel,
double[] sigma,
int nSigmas
) throws ExecutionException
{
final int area = width * height;
float[] gaussianOutput = FloatBufferPool.acquireAsIs(area);
float[] maxEigenOutput = FloatBufferPool.acquireZeroed(area);
Params params = new Params();
double maxResult = 0.0;
for (int s = 0; s < nSigmas; s++) {
computeGaussianFastMirror(
gaussianOutput,
image,
width,
height,
sigma[s]
);
double highestValue = computeEigenvalues(
params,
sliceRunner,
maxWorkers,
maxEigenOutput,
gaussianOutput,
width,
height,
sigma[s]
);
maxResult = Math.max(maxResult, highestValue);
}
// Help out the GC by nulling the second references to our recycled buffers
params.nullify();
FloatBufferPool.release(gaussianOutput);
float factor = maxResult > 0.0 ? (float)(256.0 / maxResult) : 1.0f;
//System.out.println("maxResult: " + maxResult + ", factor: " + factor);
for (int y = 1; y < height-1; y++) {
for (int x = 1; x < width-1; x++)
output[x+width*y] = (byte)(Math.min(maxEigenOutput[x+width*y] * factor, 255.0f));
}
FloatBufferPool.release(maxEigenOutput);
if (width >= 2) {
for (int y = 1; y < height-1; y++) {
output[width*y] = output[1+width*y];
output[width*(y+1)-1] = output[width*(y+1)-2];
}
}
if (height >= 2) {
for (int x = 1; x < width-1; x++)
output[x] = output[x+width];
for (int x = 1; x < width-1; x++)
output[x+width*(height-1)] = output[x+width*(height-2)];
}
if (width >= 2 && height >= 2) {
output[0] = output[width+1];
output[width-1] = output[2*width-2];
output[width*(height-1)] = output[1+width*(height-2)];
output[width*height-1] = output[width*(height-1)-2];
}
}
private static void computeGaussianFastMirror(
float[] output,
float[] image,
int width,
int height,
double sigma
) {
int kernelSize = determineGaussianKernelSize(sigma);
float[] kernel = FloatBufferPool.acquireAsIs(kernelSize);
populateGaussianKernel1D(kernel, sigma);
final int ksHalf = kernelSize / 2;
final int edgeSizeX = Math.min(ksHalf + 1, width / 2);
final int edgeSizeY = Math.min(ksHalf + 1, height / 2);
for (int y = 0; y < height; y++) {
for (int x = 0; x < edgeSizeX; x++) {
float avg = 0.0f;
for (int f = -ksHalf; f <= ksHalf; f++) {
int xf = Math.abs(x+f) % (2*width-2);
int xfm = Math.min(xf, 2*width-2 - xf);
avg += kernel[f + ksHalf] * image[xfm + width * y];
}
output[x + width * y] = avg;
}
for (int x = edgeSizeX; x < width - edgeSizeX; x++) {
float avg = 0.0f;
for (int f = -ksHalf; f <= ksHalf; f++)
avg += kernel[f + ksHalf] * image[x+f + width * y];
output[x + width * y] = avg;
}
for (int x = width - edgeSizeX; x < width; x++) {
float avg = 0.0f;
for (int f = -ksHalf; f <= ksHalf; f++) {
int xf = Math.abs(x+f) % (2*width-2);
int xfm = Math.min(xf, 2*width-2 - xf);
avg += kernel[f + ksHalf] * image[xfm + width * y];
}
output[x + width * y] = avg;
}
}
float[] tempColumn = FloatBufferPool.acquireAsIs(height);
for (int x = 0; x < width; x++) {
for (int y = 0; y < edgeSizeY; y++) {
float avg = 0.0f;
for (int f = -ksHalf; f <= ksHalf; f++) {
int yf = Math.abs(y+f) % (2*height-2);
int yfm = Math.min(yf, 2*height-2 - yf);
avg += kernel[f + ksHalf] * output[x + width * yfm];
}
tempColumn[y] = avg;
}
for (int y = edgeSizeY; y < height - edgeSizeY; y++) {
float avg = 0.0f;
for (int f = -ksHalf; f <= ksHalf; f++)
avg += kernel[f + ksHalf] * output[x + width * (y+f)];
tempColumn[y] = avg;
}
for (int y = height - edgeSizeY; y < height; y++) {
float avg = 0.0f;
for (int f = -ksHalf; f <= ksHalf; f++) {
int yf = Math.abs(y+f) % (2*height-2);
int yfm = Math.min(yf, 2*height-2 - yf);
avg += kernel[f + ksHalf] * output[x + width * yfm];
}
tempColumn[y] = avg;
}
for (int y = 0; y < height; y++)
output[x + width * y] = tempColumn[y];
}
FloatBufferPool.release(tempColumn);
FloatBufferPool.release(kernel);
}
private static int determineGaussianKernelSize(double sigma)
{
if (sigma <= 0.0F)
return 3;
double twoSigmaSq = 2.0 * sigma * sigma;
int size = Math.max(3, 2 * (int)(3.0 * sigma + 0.5) + 1);
return size;
}
private static void populateGaussianKernel1D(float[] kernel, double sigma)
{
if (sigma <= 0.0F) {
kernel[0] = 0.0f;
kernel[1] = 1.0f;
kernel[2] = 0.0f;
return;
}
double twoSigmaSq = 2.0 * sigma * sigma;
int size = Math.max(3, 2 * (int)(3.0 * sigma + 0.5) + 1);
for (int x = size / 2; x >= 0; --x) {
float val = (float)Math.exp((double)(-x * x) / twoSigmaSq);
kernel[size / 2 - x] = val;
kernel[size / 2 + x] = val;
}
float sum = 0.0F;
for (int i = 0; i < size; ++i)
sum += kernel[i];
for (int i = 0; i < size; ++i)
kernel[i] /= sum;
}
static double computeEigenvalues(
Params params,
ISliceRunner runner,
int maxWorkers,
float[] output,
float[] input,
int width,
int height,
double sigma
) throws ExecutionException
{
params.setup(input, width, height, sigma, 3, output);
runner.runSlices(
params,
maxWorkers,
width,
IN_PLACE_THRESHOLD - 1
);
return params.finalMaximum;
}
static float findSecondHessianEigenvalueAtPoint2D(float[] data, int width, int x, int y, double sigma)
{
double s2 = sigma * sigma;
int pos = x + width*y;
float dblCenter = 2.0F * data[pos];
float corners = (
((data[pos+1+width] - data[pos-1+width]) / 2.0f) -
((data[pos+1-width] - data[pos-1-width]) / 2.0f)
) / 2.0f;
// mA = matrix[0][0]
// mB = matrix[0][1]
// mC = matrix[1][0]
// mD = matrix[1][1]
double mA = (data[pos+1] - dblCenter + data[pos-1]) * s2;
double mB = corners * s2;
double mC = mB;
double mD = (data[pos+width] - dblCenter + data[pos-width]) * s2;
double a = 1.0;
double b = -(mA + mD);
double c = mA * mD - mB * mB;
double discriminant = b * b - 4.0 * a * c;
if (discriminant < 0.0)
return 0.0f;
double dR = Math.sqrt(discriminant);
float e0 = (float)((-b + dR) / (2.0 * a));
float e1 = (float)((-b - dR) / (2.0 * a));
// return the second value, accounting for inverse
return Math.abs(e0) <= Math.abs(e1) ? e1 : e0;
}
static class Params implements ISliceCompute
{
public float[] data;
public int width;
public int height;
public double sigma;
public float threshold;
public float[] output;
public DoubleVector maximums = new DoubleVector();
public double finalMaximum;
public void setup(float[] data, int width, int height, double sigma, float threshold, float[] output)
{
this.data = data;
this.width = width;
this.height = height;
this.sigma = sigma;
this.threshold = threshold;
this.output = output;
this.finalMaximum = 0.0;
}
public void nullify()
{
this.data = null;
this.output = null;
this.maximums = null;
}
@Override
public void initSlices(int nSlices)
{
this.maximums.resize(nSlices);
Arrays.fill(this.maximums.buf, 0, nSlices, 0.0);
}
@Override
public Object computeSlice(int sliceIdx, int start, int length)
{
int end = start + length;
if (start <= 0)
start = 1;
if (end >= width)
end = width - 1;
for (int y = 1; y < height - 1; ++y) {
for (int x = start; x < end; ++x) {
int index = x + width*y;
if (data[index] > threshold) {
float ev = findSecondHessianEigenvalueAtPoint2D(data, width, x, y, sigma);
if (ev < 0.0f) {
output[index] = Math.max(output[index], -ev);
maximums.buf[sliceIdx] = Math.max(maximums.buf[sliceIdx], -ev);
}
}
}
}
return null;
}
@Override
public void finishSlice(ISliceCompute.Result res)
{
finalMaximum = Math.max(finalMaximum, maximums.buf[res.idx]);
}
}
}