package features;
import Utils.Utils;
import ij.ImagePlus;
import ij.process.ByteProcessor;
import ij.process.FloatProcessor;
import ij.process.ImageProcessor;
import ij.process.ShortProcessor;
import java.util.concurrent.ForkJoinPool;
import java.util.concurrent.RecursiveAction;
import math3d.Eigensystem2x2Float;
import math3d.Eigensystem3x3Float;
import math3d.JacobiFloat;
public class ForkEigenValuesAtPoint2D extends RecursiveAction {
float[] evalues = new float[3];
protected static ForkEigenValuesAtPoint2D.FloatArray2D data2D;
protected boolean fixUp = false;
protected boolean orderOnAbsoluteSize = true;
protected boolean normalize = false;
float sepX = 1.0F;
float sepY = 1.0F;
int width = 0;
int height = 0;
ImagePlus original;
int threshold;
public float[] sliceFinal;
float minResult = Float.MAX_VALUE;
float maxResult = Float.MIN_VALUE;
int mLength;
int mStartA;
int mStartB;
protected double sigma;
protected static int sThreshold = 5;
public ForkEigenValuesAtPoint2D(
ImagePlus original, ForkEigenValuesAtPoint2D.FloatArray2D data, int start, int length, double sigma, int threshold, float[] dst
) {
this.original = original;
this.width = original.getWidth();
this.height = original.getHeight();
data2D = data;
this.mLength = length;
this.mStartA = this.mStartB = start;
if (this.mStartA == 0) {
this.mStartA = 1;
}
if (this.mStartB + this.mLength >= this.width - 1) {
--this.mLength;
}
this.sigma = sigma;
this.threshold = threshold;
this.sliceFinal = dst;
}
ForkEigenValuesAtPoint2D() {
}
private void computeDirectly() {
long count = 0L;
long total = (long)(this.height * this.width);
ImageProcessor fp = this.original.getProcessor();
for(int y = 1; y < this.height - 1; ++y) {
for(int x = this.mStartA; x < this.mStartB + this.mLength; ++x) {
if (fp.getPixelValue(x, y) > (float)this.threshold) {
boolean real = this.hessianEigenvaluesAtPoint2D(x, y, true, this.evalues, this.normalize, false, this.sepX, this.sepY);
int index = y * this.width + x;
float value = 0.0F;
if (real) {
value = measureFromEvalues2D(this.evalues, 1);
}
this.sliceFinal[index] = value;
}
++count;
}
}
}
@Override
protected void compute() {
if (this.mLength < sThreshold) {
this.computeDirectly();
} else {
int split = this.mLength / 2;
int remainder = this.mLength % 2;
int secondHalf = split + remainder;
invokeAll(
new ForkEigenValuesAtPoint2D(this.original, data2D, this.mStartB, split, this.sigma, this.threshold, this.sliceFinal),
new ForkEigenValuesAtPoint2D(this.original, data2D, this.mStartB + split, secondHalf, this.sigma, this.threshold, this.sliceFinal)
);
}
}
public float[] computeEigenvalues(ImagePlus original, double sigma, int _threshold) {
//System.out.println("_threshold= " + _threshold + "\ttreshold = " + this.threshold);
data2D = this.ImageToFloatArray(original.getProcessor());
int[] histogram2 = Utils.getFloatHistogram2(original);
this.threshold = Utils.findHistogramMax(histogram2) + 2;
/*
if (!Utils.isReleaseVersion) {
System.out.println("treshold = " + this.threshold);
}
*/
if (_threshold <= 0) {
this.threshold = 3;
} else {
this.threshold = _threshold;
}
//System.out.println("_threshold= " + _threshold + "\ttreshold = " + this.threshold);
this.sliceFinal = new float[data2D.data.length];
ForkEigenValuesAtPoint2D fe = new ForkEigenValuesAtPoint2D(original, data2D, 0, original.getWidth(), sigma, this.threshold, this.sliceFinal);
ForkJoinPool pool = new ForkJoinPool();
long startTime = System.currentTimeMillis();
pool.invoke(fe);
long endTime = System.currentTimeMillis();
return this.sliceFinal;
}
public boolean hessianEigenvaluesAtPoint2D(
int x, int y, boolean orderOnAbsoluteSize, float[] result, boolean normalize, boolean fixUp, float sepX, float sepY
) {
if (fixUp) {
if (x == 0) {
x = 1;
}
if (x == data2D.width - 1) {
x = data2D.width - 2;
}
if (y == 0) {
y = 1;
}
if (y == data2D.height - 1) {
y = data2D.height - 2;
}
}
float[][] hessianMatrix = this.computeHessianMatrix2DFloat(data2D, x, y, this.sigma, sepX, sepY);
float[] eigenValues = this.findEigenValuesForMatrix(hessianMatrix);
if (eigenValues == null) {
return false;
} else {
float e0 = eigenValues[0];
float e1 = eigenValues[1];
float e0c = orderOnAbsoluteSize ? Math.abs(e0) : e0;
float e1c = orderOnAbsoluteSize ? Math.abs(e1) : e1;
if (e0c <= e1c) {
result[0] = e0;
result[1] = e1;
} else {
result[0] = e1;
result[1] = e0;
}
if (normalize) {
float divideBy = Math.abs(result[1]);
result[0] /= divideBy;
result[1] /= divideBy;
}
return true;
}
}
public static float measureFromEvalues2D(float[] evalues, int vesselness) {
float measure = 0.0F;
return evalues[1] >= 0.0F ? 0.0F : Math.abs(evalues[1]);
}
public long[] hessianEigenvaluesAtPoint2D2(
int x, int y, boolean orderOnAbsoluteSize, float[] result, boolean normalize, boolean fixUp, float sepX, float sepY
) {
long start = System.currentTimeMillis();
if (fixUp) {
if (x == 0) {
x = 1;
}
if (x == data2D.width - 1) {
x = data2D.width - 2;
}
if (y == 0) {
y = 1;
}
if (y == data2D.height - 1) {
y = data2D.height - 2;
}
}
long hessianMatrixStart = System.currentTimeMillis();
float[][] hessianMatrix = this.computeHessianMatrix2DFloat(data2D, x, y, this.sigma, sepX, sepY);
long hessianMatrixEnd = System.currentTimeMillis();
float[] eigenValues = this.findEigenValuesForMatrix(hessianMatrix);
long eigenValuesEnd = System.currentTimeMillis();
if (eigenValues == null) {
return null;
} else {
float e0 = eigenValues[0];
float e1 = eigenValues[1];
float e0c = orderOnAbsoluteSize ? Math.abs(e0) : e0;
float e1c = orderOnAbsoluteSize ? Math.abs(e1) : e1;
if (e0c <= e1c) {
result[0] = e0;
result[1] = e1;
} else {
result[0] = e1;
result[1] = e0;
}
if (normalize) {
float divideBy = Math.abs(result[1]);
result[0] /= divideBy;
result[1] /= divideBy;
}
long end1 = System.currentTimeMillis();
long z1 = hessianMatrixEnd - hessianMatrixStart;
long z2 = eigenValuesEnd - hessianMatrixEnd;
long z3 = end1 - eigenValuesEnd;
long end = System.currentTimeMillis();
long z4 = end - start;
return new long[]{z1, z2, z3, z4};
}
}
public float[][] computeHessianMatrix2DFloat(ForkEigenValuesAtPoint2D.FloatArray2D laPlace, int x, int y, double sigma, float sepX, float sepY) {
float[][] hessianMatrix = new float[2][2];
float temp = 2.0F * laPlace.get(x, y);
hessianMatrix[0][0] = laPlace.get(x + 1, y) - temp + laPlace.get(x - 1, y);
hessianMatrix[1][1] = laPlace.get(x, y + 1) - temp + laPlace.get(x, y - 1);
hessianMatrix[0][1] = hessianMatrix[1][0] = (
(laPlace.get(x + 1, y + 1) - laPlace.get(x - 1, y + 1)) / 2.0F - (laPlace.get(x + 1, y - 1) - laPlace.get(x - 1, y - 1)) / 2.0F
)
/ 2.0F;
for(int i = 0; i < 2; ++i) {
for(int j = 0; j < 2; ++j) {
hessianMatrix[i][j] = (float)((double)hessianMatrix[i][j] * sigma * sigma);
}
}
return hessianMatrix;
}
public float[] findEigenValuesForMatrix(float[][] matrix) {
if (matrix.length == 3 && matrix[0].length == 3) {
Eigensystem3x3Float e = new Eigensystem3x3Float(matrix);
boolean result = e.findEvalues();
return result ? e.getEvaluesCopy() : null;
} else if (matrix.length == 2 && matrix[0].length == 2) {
Eigensystem2x2Float e = new Eigensystem2x2Float(matrix);
boolean result = e.findEvalues();
return result ? e.getEvaluesCopy() : null;
} else {
JacobiFloat jc = new JacobiFloat(matrix, 50);
return jc.getEigenValues();
}
}
public ForkEigenValuesAtPoint2D.FloatArray2D ImageToFloatArray(ImageProcessor ip) {
Object pixelArray = ip.getPixels();
int size = ip.getWidth() * ip.getHeight();
ForkEigenValuesAtPoint2D.FloatArray2D image = null;
if (ip instanceof ByteProcessor) {
image = new ForkEigenValuesAtPoint2D.FloatArray2D(ip.getWidth(), ip.getHeight());
byte[] pixels = (byte[])pixelArray;
for(int i = 0; i < size; i++) {
image.data[i] = (float)(pixels[i] & 0xff);
}
} else if (ip instanceof ShortProcessor) {
image = new ForkEigenValuesAtPoint2D.FloatArray2D(ip.getWidth(), ip.getHeight());
short[] pixels = (short[])pixelArray;
for(int i = 0; i < size; i++) {
image.data[i] = (float)(pixels[i] & 0xffff);
}
} else if (ip instanceof FloatProcessor) {
image = new ForkEigenValuesAtPoint2D.FloatArray2D(ip.getWidth(), ip.getHeight());
System.arraycopy((float[])pixelArray, 0, image.data, 0, size);
} else {
if (!Utils.isReleaseVersion) {
System.err.println("RGB images not supported");
}
}
return image;
}
public static abstract class FloatArray {
public float[] data = null;
public abstract ForkEigenValuesAtPoint2D.FloatArray clone();
}
public static class FloatArray2D extends ForkEigenValuesAtPoint2D.FloatArray {
public float[] data = null;
public int width = 0;
public int height = 0;
public FloatArray2D(int width, int height) {
this.data = new float[width * height];
this.width = width;
this.height = height;
}
public FloatArray2D(float[] data, int width, int height) {
this.data = data;
this.width = width;
this.height = height;
}
public ForkEigenValuesAtPoint2D.FloatArray2D clone() {
ForkEigenValuesAtPoint2D.FloatArray2D clone = new ForkEigenValuesAtPoint2D.FloatArray2D(this.width, this.height);
System.arraycopy(this.data, 0, clone.data, 0, this.data.length);
return clone;
}
public int getPos(int x, int y) {
return x + this.width * y;
}
public float get(int x, int y) {
return this.data[this.getPos(x, y)];
}
public float getMirror(int x, int y) {
if (x >= this.width) {
x = this.width - (x - this.width + 2);
}
if (y >= this.height) {
y = this.height - (y - this.height + 2);
}
if (x < 0) {
int tmp = 0;
for(int dir = 1; x < 0; ++x) {
tmp += dir;
if (tmp == this.width - 1 || tmp == 0) {
dir *= -1;
}
}
x = tmp;
}
if (y < 0) {
int tmp = 0;
for(int dir = 1; y < 0; ++y) {
tmp += dir;
if (tmp == this.height - 1 || tmp == 0) {
dir *= -1;
}
}
y = tmp;
}
return this.data[this.getPos(x, y)];
}
public float getZero(int x, int y) {
if (x >= this.width) {
return 0.0F;
} else if (y >= this.height) {
return 0.0F;
} else if (x < 0) {
return 0.0F;
} else {
return y < 0 ? 0.0F : this.data[this.getPos(x, y)];
}
}
public void set(float value, int x, int y) {
this.data[this.getPos(x, y)] = value;
}
}
}