use of gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction in project GDSC-SMLM by aherbert.
the class LSQLVMGradientProcedureTest method doubleCreateGaussianData.
/**
* Create random elliptical Gaussian data an returns the data plus an estimate of the parameters.
* Only the chosen parameters are randomised and returned for a maximum of (background, amplitude, angle, xpos,
* ypos, xwidth, ywidth }
*
* @param npeaks
* the npeaks
* @param params
* set on output
* @param randomiseParams
* Set to true to randomise the params
* @return the double[]
*/
private double[] doubleCreateGaussianData(int npeaks, double[] params, boolean randomiseParams) {
int n = blockWidth * blockWidth;
// Generate a 2D Gaussian
ErfGaussian2DFunction func = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(npeaks, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
params[0] = random(Background);
for (int i = 0, j = 1; i < npeaks; i++, j += 6) {
params[j] = random(Signal);
params[j + 2] = random(Xpos);
params[j + 3] = random(Ypos);
params[j + 4] = random(Xwidth);
params[j + 5] = random(Ywidth);
}
double[] y = new double[n];
func.initialise(params);
for (int i = 0; i < y.length; i++) {
// Add random Poisson noise
y[i] = rdg.nextPoisson(func.eval(i));
}
if (randomiseParams) {
params[0] = random(params[0]);
for (int i = 0, j = 1; i < npeaks; i++, j += 6) {
params[j] = random(params[j]);
params[j + 2] = random(params[j + 2]);
params[j + 3] = random(params[j + 3]);
params[j + 4] = random(params[j + 4]);
params[j + 5] = random(params[j + 5]);
}
}
return y;
}
use of gdsc.smlm.function.gaussian.erf.ErfGaussian2DFunction in project GDSC-SMLM by aherbert.
the class FastMLEGradient2ProcedureTest method gradientProcedureComputesSameWithPrecomputed.
@Test
public void gradientProcedureComputesSameWithPrecomputed() {
int iter = 10;
rdg = new RandomDataGenerator(new Well19937c(30051977));
ErfGaussian2DFunction f1 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(1, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
ErfGaussian2DFunction f2 = (ErfGaussian2DFunction) GaussianFunctionFactory.create2D(2, 10, 10, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
double[] a1 = new double[7];
double[] a2 = new double[13];
final double[] x = new double[f1.size()];
final double[] b = new double[f1.size()];
for (int i = 0; i < iter; i++) {
a2[Gaussian2DFunction.BACKGROUND] = rdg.nextUniform(0.1, 0.3);
a2[Gaussian2DFunction.SIGNAL] = rdg.nextUniform(100, 300);
a2[Gaussian2DFunction.X_POSITION] = rdg.nextUniform(3, 5);
a2[Gaussian2DFunction.Y_POSITION] = rdg.nextUniform(3, 5);
a2[Gaussian2DFunction.X_SD] = rdg.nextUniform(1, 1.3);
a2[Gaussian2DFunction.Y_SD] = rdg.nextUniform(1, 1.3);
a2[6 + Gaussian2DFunction.SIGNAL] = rdg.nextUniform(100, 300);
a2[6 + Gaussian2DFunction.X_POSITION] = rdg.nextUniform(5, 7);
a2[6 + Gaussian2DFunction.Y_POSITION] = rdg.nextUniform(5, 7);
a2[6 + Gaussian2DFunction.X_SD] = rdg.nextUniform(1, 1.3);
a2[6 + Gaussian2DFunction.Y_SD] = rdg.nextUniform(1, 1.3);
// Simulate Poisson data
f2.initialise0(a2);
f1.forEach(new ValueProcedure() {
int k = 0;
public void execute(double value) {
x[k++] = (value > 0) ? rdg.nextPoisson(value) : 0;
}
});
// Precompute peak 2 (no background)
a1[Gaussian2DFunction.BACKGROUND] = 0;
for (int j = 1; j < 7; j++) a1[j] = a2[6 + j];
f1.initialise0(a1);
f1.forEach(new ValueProcedure() {
int k = 0;
public void execute(double value) {
b[k++] = value;
}
});
// Reset to peak 1
for (int j = 0; j < 7; j++) a1[j] = a2[j];
// Compute peak 1+2
FastMLEGradient2Procedure p12 = FastMLEGradient2ProcedureFactory.create(x, f2);
p12.computeSecondDerivative(a2);
double[] d11 = Arrays.copyOf(p12.d1, f1.getNumberOfGradients());
double[] d21 = Arrays.copyOf(p12.d2, f1.getNumberOfGradients());
// Compute peak 1+(precomputed 2)
FastMLEGradient2Procedure p1b2 = FastMLEGradient2ProcedureFactory.create(x, PrecomputedGradient2Function.wrapGradient2Function(f1, b));
p1b2.computeSecondDerivative(a1);
double[] d12 = p1b2.d1;
double[] d22 = p1b2.d2;
Assert.assertArrayEquals(" Result: Not same @ " + i, p12.u, p1b2.u, 1e-10);
Assert.assertArrayEquals(" D1: Not same @ " + i, d11, d12, 1e-10);
Assert.assertArrayEquals(" D2: Not same @ " + i, d21, d22, 1e-10);
double[] v1 = p12.computeValue(a2);
double[] v2 = p1b2.computeValue(a1);
Assert.assertArrayEquals(" Value: Not same @ " + i, v1, v2, 1e-10);
double[] d1 = Arrays.copyOf(p12.computeFirstDerivative(a2), f1.getNumberOfGradients());
double[] d2 = p1b2.computeFirstDerivative(a1);
Assert.assertArrayEquals(" 1st derivative: Not same @ " + i, d1, d2, 1e-10);
}
}
Aggregations