use of org.apache.commons.math3.optim.nonlinear.scalar.ObjectiveFunction in project vcell by virtualcell.
the class FitBleachSpotOp method fitToGaussian.
static GaussianFitResults fitToGaussian(double init_center_i, double init_center_j, double init_radius2, FloatImage image) {
//
// do some optimization on the image (fitting to a Gaussian)
// set initial guesses from ROI operation.
//
ISize imageSize = image.getISize();
final int num_i = imageSize.getX();
final int num_j = imageSize.getY();
final float[] floatPixels = image.getFloatPixels();
//
// initial guess based on previous fit of ROI
// do gaussian fit in index space for center and standard deviation (later to translate it back to world coordinates)
//
final int window_size = (int) Math.sqrt(init_radius2) * 4;
// final int window_min_i = 0; // (int) Math.max(0, Math.floor(init_center_i - window_size/2));
// final int window_max_i = num_i-1; // (int) Math.min(num_i-1, Math.ceil(init_center_i + window_size/2));
// final int window_min_j = 0; // (int) Math.max(0, Math.floor(init_center_j - window_size/2));
// final int window_max_j = num_j-1; // (int) Math.min(num_j-1, Math.ceil(init_center_j + window_size/2));
final int window_min_i = (int) Math.max(0, Math.floor(init_center_i - window_size / 2));
final int window_max_i = (int) Math.min(num_i - 1, Math.ceil(init_center_i + window_size / 2));
final int window_min_j = (int) Math.max(0, Math.floor(init_center_j - window_size / 2));
final int window_max_j = (int) Math.min(num_j - 1, Math.ceil(init_center_j + window_size / 2));
final int PARAM_INDEX_CENTER_I = 0;
final int PARAM_INDEX_CENTER_J = 1;
final int PARAM_INDEX_K = 2;
final int PARAM_INDEX_HIGH = 3;
final int PARAM_INDEX_RADIUS_SQUARED = 4;
final int NUM_PARAMETERS = 5;
double[] initParameters = new double[NUM_PARAMETERS];
initParameters[PARAM_INDEX_CENTER_I] = init_center_i;
initParameters[PARAM_INDEX_CENTER_J] = init_center_j;
initParameters[PARAM_INDEX_HIGH] = 1.0;
initParameters[PARAM_INDEX_K] = 10;
initParameters[PARAM_INDEX_RADIUS_SQUARED] = init_radius2;
PowellOptimizer optimizer = new PowellOptimizer(1e-4, 1e-1);
MultivariateFunction func = new MultivariateFunction() {
@Override
public double value(double[] point) {
double center_i = point[PARAM_INDEX_CENTER_I];
double center_j = point[PARAM_INDEX_CENTER_J];
double high = point[PARAM_INDEX_HIGH];
double K = point[PARAM_INDEX_K];
double radius2 = point[PARAM_INDEX_RADIUS_SQUARED];
double error2 = 0;
for (int j = window_min_j; j <= window_max_j; j++) {
// double y = j - center_j;
double y = j;
for (int i = window_min_i; i <= window_max_i; i++) {
// double x = i - center_i;
double x = i;
double modelValue = high - FastMath.exp(-K * FastMath.exp(-2 * (x * x + y * y) / radius2));
double imageValue = floatPixels[j * num_i + i];
double error = modelValue - imageValue;
error2 += error * error;
}
}
System.out.println(new GaussianFitResults(center_i, center_j, radius2, K, high, error2));
return error2;
}
};
PointValuePair pvp = optimizer.optimize(new ObjectiveFunction(func), new InitialGuess(initParameters), new MaxEval(100000), GoalType.MINIMIZE);
double[] fittedParamValues = pvp.getPoint();
double fitted_center_i = fittedParamValues[PARAM_INDEX_CENTER_I];
double fitted_center_j = fittedParamValues[PARAM_INDEX_CENTER_J];
double fitted_radius2 = fittedParamValues[PARAM_INDEX_RADIUS_SQUARED];
double fitted_K = fittedParamValues[PARAM_INDEX_K];
double fitted_high = fittedParamValues[PARAM_INDEX_HIGH];
double objectiveFunctionValue = pvp.getValue();
return new GaussianFitResults(fitted_center_i, fitted_center_j, fitted_radius2, fitted_K, fitted_high, objectiveFunctionValue);
}
Aggregations