ij.measure.CurveFitter Maven / Gradle / Ivy
Go to download
Show more of this group Show more artifacts with this name
Show all versions of ij Show documentation
Show all versions of ij Show documentation
ImageJ is an open source Java image processing program inspired by NIH Image for the Macintosh.
package ij.measure;
import ij.*;
import ij.gui.*;
import ij.macro.*;
import ij.util.Tools;
import java.util.Arrays;
import java.util.Hashtable;
/** Curve fitting class based on the Simplex method in the Minimizer class
*
*
* Notes on fitting polynomial functions:
* (i) The range of x values should not be too far from 0, especially for higher-order polynomials.
* For polynomials of 4th order, the average x value should fulfill |xMean| < 2*(xMax-xMin).
* For polynomials of 5th order or higher, the x range should encompass x=0; and for
* 7th and 8th order it is desirable to have x=0 near the center of the x range.
* (ii) In contrast to some fitting algorithms using the normal equations and matrix inversion, the
* simplex algorithm used here can cope with parameters having very different orders of magnitude,
* as long as the coefficients of the polynomial are within a reasonable range (say, 1e-80 to 1e80).
* Thus, it is usually not needed to scale the x values, even for high-order polynomials.
*
* Version history:
*
* 2008-01-21: Modified to do Gaussian fitting by Stefan Woerz (s.woerz at dkfz.de).
* 2012-01-30: Modified for external Minimizer class and UserFunction fit by Michael Schmid.
* - Also fixed incorrect equation string for 'Gamma Variate' & 'Rodbard (NIH Image)',
* - Added 'Inverse Rodbard', 'Exponential (linear regression)', 'Power (linear regression)'
* functions and polynomials of order 5-8. Added 'nicely' sorted list of types.
* - Added absolute error for minimizer to avoid endless minimization if exact fit is possible.
* - Added 'initialParamVariations' (order of magnitude if parameter variations) -
* this is important for safer and better convergence.
* - Linear regression for up to 2 linear parameters, reduces the number of parameters handled
* by the simplex Minimizer and improves convergence. These parameters can be an offset and
* either a linear slope or a factor that the full function is multiplied with.
* 2012-10-07: added GAUSSIAN_NOOFFSET fit type
* 2012-11-20: Bugfix: exception on Gaussian&Rodbard with initial params, bad initial params for Gaussian
* 2013-09-24: Added "Exponential Recovery (no offset)" and "Chapman-Richards" (3-parameter) fit types.
* 2013-10-11: bugfixes, added setStatusAndEsc to show iterations and enable abort by ESC
* 2015-03-26: bugfix, did not use linear regression for RODBARD
*/
public class CurveFitter implements UserFunction{
/** Constants for the built-in fit functions */
public static final int STRAIGHT_LINE=0, POLY2=1, POLY3=2, POLY4=3,
EXPONENTIAL=4, POWER=5, LOG=6, RODBARD=7, GAMMA_VARIATE=8, LOG2=9,
RODBARD2=10, EXP_WITH_OFFSET=11, GAUSSIAN=12, EXP_RECOVERY=13,
INV_RODBARD=14, EXP_REGRESSION=15, POWER_REGRESSION=16,
POLY5=17, POLY6=18, POLY7=19, POLY8=20,
GAUSSIAN_NOOFFSET=21, EXP_RECOVERY_NOOFFSET=22, CHAPMAN=23;
/** Nicer sequence of the built-in function types */
public static final int[] sortedTypes = { STRAIGHT_LINE, POLY2, POLY3, POLY4,
POLY5, POLY6, POLY7, POLY8,
POWER, POWER_REGRESSION,
EXPONENTIAL, EXP_REGRESSION, EXP_WITH_OFFSET,
EXP_RECOVERY, EXP_RECOVERY_NOOFFSET,
LOG, LOG2,
GAUSSIAN, GAUSSIAN_NOOFFSET,
RODBARD, RODBARD2, INV_RODBARD,
GAMMA_VARIATE,CHAPMAN
};
/** Names of the built-in fit functions */
public static final String[] fitList = {"Straight Line","2nd Degree Polynomial",
"3rd Degree Polynomial", "4th Degree Polynomial","Exponential","Power",
"Log","Rodbard", "Gamma Variate", "y = a+b*ln(x-c)","Rodbard (NIH Image)",
"Exponential with Offset","Gaussian", "Exponential Recovery",
"Inverse Rodbard", "Exponential (linear regression)", "Power (linear regression)",
"5th Degree Polynomial","6th Degree Polynomial","7th Degree Polynomial","8th Degree Polynomial",
"Gaussian (no offset)", "Exponential Recovery (no offset)",
"Chapman-Richards"
}; // fList, doFit(), getNumParams() and makeInitialParamsAndVariations() must also be updated
/** Equations of the built-in fit functions */
public static final String[] fList = {
"y = a+bx","y = a+bx+cx^2", //STRAIGHT_LINE,POLY2
"y = a+bx+cx^2+dx^3","y = a+bx+cx^2+dx^3+ex^4",
"y = a*exp(bx)","y = a*x^b", "y = a*ln(bx)", //EXPONENTIAL,POWER,LOG
"y = d+(a-d)/(1+(x/c)^b)", "y = b*(x-a)^c*exp(-(x-a)/d)", //RODBARD,GAMMA_VARIATE
"y = a+b*ln(x-c)", "x = d+(a-d)/(1+(y/c)^b) [y = c*((x-a)/(d-x))^(1/b)]", //LOG2,RODBARD2
"y = a*exp(-bx) + c", "y = a + (b-a)*exp(-(x-c)*(x-c)/(2*d*d))", //EXP_WITH_OFFSET,GAUSSIAN
"y = a*(1-exp(-b*x)) + c", "y = c*((x-a)/(d-x))^(1/b)", //EXP_RECOVERY, INV_RODBARD
"y = a*exp(bx)", "y = a*x^b", //EXP_REGRESSION, POWER_REGRESSION
"y = a+bx+cx^2+dx^3+ex^4+fx^5", "y = a+bx+cx^2+dx^3+ex^4+fx^5+gx^6",
"y = a+bx+cx^2+dx^3+ex^4+fx^5+gx^6+hx^7", "y = a+bx+cx^2+dx^3+ex^4+fx^5+gx^6+hx^7+ix^8",
"y = a*exp(-(x-b)*(x-b)/(2*c*c)))", //GAUSSIAN_NOOFFSET
"y = a*(1-exp(-b*x))", //EXP_RECOVERY_NOOFFSET
"y = a*(1-exp(-b*x))^c" //CHAPMAN
};
/** @deprecated now in the Minimizer class (since ImageJ 1.46f).
* (probably of not much value for anyone anyhow?) */
public static final int IterFactor = 500;
private static final int CUSTOM = 100; // user function defined in macro or plugin
private static final int GAUSSIAN_INTERNAL = 101; // Gaussian with separate offset & multiplier
private static final int RODBARD_INTERNAL = 102; // Rodbard with separate offset & multiplier
private int fitType = -1; // Number of curve type to fit
private double[] xData, yData; // x,y data to fit
private double[] xDataSave, yDataSave; //saved original data after fitting modified data
private int numPoints; // number of data points in actual fit
private double ySign = 0; // remember sign of y data for power-law fit via regression
private double sumY = Double.NaN, sumY2 = Double.NaN; // sum(y), sum(y^2) of the data used for fitting
private int numParams; // number of parameters
private double[] initialParams; // user specified or estimated initial parameters
private double[] initialParamVariations; // estimate of range of parameters
private double[] minimizerInitialParams; // modified initialParams of the minimizer
private double[] minimizerInitialParamVariations; // modified initialParamVariations of the minimizer
private double maxRelError = 1e-10;// maximum relative deviation of sum of residuals^2 from minimum
private long time; // elapsed time in ms
private int customParamCount; // number of parameters of user-supplied function (macro or plugin)
private String customFormula; // equation defined in macro or text
private UserFunction userFunction; // plugin with custom fit function
private Interpreter macro; // holds macro with custom equation
private int macroStartProgramCounter; // equation in macro starts here
private int numRegressionParams;// number of parameters that can be calculated by linear regression
private int offsetParam = -1; // parameter number: function has this parameter as offset
private int factorParam = -1; // index of parameter that is slope or multiplicative factor of the function
private boolean hasSlopeParam; // whether the 'factorParam' is the slope of the function
private double[] finalParams; // parameters obtained by fit
private boolean linearRegressionUsed; // whether linear regression alone was used instead of the minimizer
private boolean restrictPower; // power via linear regression fit: (0,0) requires positive power
private Minimizer minimizer = new Minimizer();
private int minimizerStatus = Minimizer.INITIALIZATION_FAILURE; // status of the minimizer after minimizing
private String errorString; // in case of error before invoking the minimizer
private static String[] sortedFitList; // names like fitList, but in more logical sequence
private static Hashtable namesTable; // converts fitList String into number
/** Construct a new CurveFitter. */
public CurveFitter (double[] xData, double[] yData) {
this.xData = xData;
this.yData = yData;
numPoints = xData.length;
}
/** Perform curve fitting with one of the built-in functions
* doFit(fitType) does the fit quietly
* Use getStatus() and/or getStatusString() to see whether fitting was (probably) successful and
* getParams() to access the result.
*/
public void doFit(int fitType) {
doFit(fitType, false);
}
/** Perform curve fitting with one of the built-in functions
* doFit(fitType, true) pops up a dialog allowing the user to set the initial
* fit parameters and various numbers controlling the Minimizer
* Use getStatus() and/or getStatusString() to see whether fitting was (probably) successful and
* getParams() to access the result.
*/
public void doFit(int fitType, boolean showSettings) {
if (!(fitType>=STRAIGHT_LINE && fitType 0)
minimizerParamsToFullParams(finalParams, false);
}
if (isModifiedFitType(fitType)) //params of actual fit to user params
postProcessModifiedFitType(fitType);
switch (fitType) { //postprocessing for nicer output:
case GAUSSIAN: //Gaussians: width (std deviation) should be >0
finalParams[3] = Math.abs(finalParams[3]); break;
case GAUSSIAN_NOOFFSET:
finalParams[2] = Math.abs(finalParams[2]); break;
}
time = System.currentTimeMillis()-startTime;
}
/** Fit a function defined as a macro String like "y = a + b*x + c*x*x".
* Returns the number of parameters, or 0 in case of a macro syntax error.
*
* For good performance, it is advisable to set also the typical variation range
* of the initial parameters by the
* getMinimizer().setInitialParamVariations(double[]) method (especially if one or
* more of the initialParams are zero).
* Use getStatus() and/or getStatusString() to see whether fitting was (probably) successful and
* getParams() to access the result.
*/
public int doCustomFit(String equation, double[] initialParams, boolean showSettings) {
customFormula = null;
customParamCount = 0;
Program pgm = (new Tokenizer()).tokenize(equation);
if (!pgm.hasWord("y") || !pgm.hasWord("x"))
return 0;
String[] params = {"a","b","c","d","e","f"};
for (int i=0; i= 0;
factorParam = hasSlopeParam ? slopeParam : multiplyParam;
numRegressionParams = 0;
if (offsetParam >=0) numRegressionParams++;
if (factorParam >=0) numRegressionParams++;
}
/** Get number of parameters for current fit formula
* Do not use before 'doFit', because the fit function would be undefined. */
public int getNumParams() {
switch (fitType) {
case STRAIGHT_LINE: return 2;
case POLY2: return 3;
case POLY3: return 4;
case POLY4: return 5;
case POLY5: return 6;
case POLY6: return 7;
case POLY7: return 8;
case POLY8: return 9;
case EXPONENTIAL: case EXP_REGRESSION: return 2;
case POWER: case POWER_REGRESSION: return 2;
case EXP_RECOVERY_NOOFFSET: return 2;
case LOG: return 2;
case LOG2: return 3;
case GAUSSIAN_NOOFFSET: return 3;
case EXP_RECOVERY: return 3;
case CHAPMAN: return 3;
case EXP_WITH_OFFSET: return 3;
case RODBARD: case RODBARD2: case INV_RODBARD: case RODBARD_INTERNAL: return 4;
case GAMMA_VARIATE: return 4;
case GAUSSIAN: case GAUSSIAN_INTERNAL: return 4;
case CUSTOM: return customParamCount;
}
return 0;
}
/** Returns the formula value for parameters 'p' at 'x'.
* Do not use before 'doFit', because the fit function would be undefined. */
public final double f(double x) {
if (finalParams==null)
finalParams = minimizer.getParams();
return f(finalParams, x);
}
/** Returns the formula value for parameters 'p' at 'x'.
* Do not use before 'doFit', because the fit function would be undefined. */
public final double f(double[] p, double x) {
if (fitType!=CUSTOM)
return f(fitType, p, x);
else {
if (macro==null) { // function defined in plugin
return userFunction.userFunction(p, x);
} else { // function defined in macro
macro.setVariable("x", x);
macro.setVariable("a", p[0]);
if (customParamCount>1) macro.setVariable("b", p[1]);
if (customParamCount>2) macro.setVariable("c", p[2]);
if (customParamCount>3) macro.setVariable("d", p[3]);
if (customParamCount>4) macro.setVariable("e", p[4]);
if (customParamCount>5) macro.setVariable("f", p[5]);
macro.run(macroStartProgramCounter);
return macro.getVariable("y");
}
}
}
/** Returns value of built-in 'fitType' formula value for parameters "p" at "x" */
public static double f(int fitType, double[] p, double x) {
switch (fitType) {
case STRAIGHT_LINE:
return p[0] + x*p[1];
case POLY2:
return p[0] + x*(p[1] + x*p[2]);
case POLY3:
return p[0] + x*(p[1] + x*(p[2] + x*p[3]));
case POLY4:
return p[0] + x*(p[1] + x*(p[2] + x*(p[3] + x*p[4])));
case POLY5:
return p[0] + x*(p[1] + x*(p[2] + x*(p[3] + x*(p[4] + x*p[5]))));
case POLY6:
return p[0] + x*(p[1] + x*(p[2] + x*(p[3] + x*(p[4] + x*(p[5] + x*p[6])))));
case POLY7:
return p[0] + x*(p[1] + x*(p[2] + x*(p[3] + x*(p[4] + x*(p[5] + x*(p[6] + x*p[7]))))));
case POLY8:
return p[0] + x*(p[1] + x*(p[2] + x*(p[3] + x*(p[4] + x*(p[5] + x*(p[6] + x*(p[7]+x*p[8])))))));
case EXPONENTIAL:
case EXP_REGRESSION:
return p[0]*Math.exp(p[1]*x);
case EXP_WITH_OFFSET:
return p[0]*Math.exp(-p[1]*x)+p[2];
case EXP_RECOVERY:
return p[0]*(1-Math.exp(-p[1]*x))+p[2];
case EXP_RECOVERY_NOOFFSET:
return p[0]*(1-Math.exp(-p[1]*x));
case CHAPMAN: // a*(1-exp(-b*x))^c
double value = p[0]*(Math.pow((1-(Math.exp(-p[1]*x))), p[2]));
// Log.e("test", "values = " + value);
return value;
case GAUSSIAN:
return p[0]+(p[1]-p[0])*Math.exp(-(x-p[2])*(x-p[2])/(2.0*p[3]*p[3]));
case GAUSSIAN_INTERNAL: // replaces GAUSSIAN for the fitting process
return p[0]+p[1]*Math.exp(-(x-p[2])*(x-p[2])/(2.0*p[3]*p[3]));
case GAUSSIAN_NOOFFSET:
return p[0]*Math.exp(-(x-p[1])*(x-p[1])/(2.0*p[2]*p[2]));
case POWER: // ax^b
case POWER_REGRESSION:
return p[0]*Math.pow(x,p[1]);
case LOG:
if (x == 0.0)
return -1000*p[0];
return p[0]*Math.log(p[1]*x);
case RODBARD: { // d+(a-d)/(1+(x/c)^b)
double ex = Math.pow(x/p[2], p[1]); //(x/c)^b
return p[3]+(p[0]-p[3])/(1.0+ex); }
case RODBARD_INTERNAL: { // d+a/(1+(x/c)^b) , replaces RODBARD of the fitting process
double ex = Math.pow(x/p[2], p[1]); //(x/c)^b
return p[3]+p[0]/(1.0+ex); }
case GAMMA_VARIATE: // b*(x-a)^c*exp(-(x-a)/d)
if (p[0] >= x) return 0.0;
if (p[1] <= 0) return Double.NaN;
if (p[2] <= 0) return Double.NaN;
if (p[3] <= 0) return Double.NaN;
double pw = Math.pow((x - p[0]), p[2]);
double e = Math.exp((-(x - p[0]))/p[3]);
return p[1]*pw*e;
case LOG2:
double tmp = x-p[2];
if (tmp<=0)
return Double.NaN;
return p[0]+p[1]*Math.log(tmp);
case INV_RODBARD: // c*((x-a)/(d-x))^(1/b), the inverse Rodbard function
case RODBARD2: // also used after the 'Rodbard NIH Image' fit
double y;
if (p[3]-x < 2*Double.MIN_VALUE || x=d (singularity) and x
r^2 = 1 - SSE/SSD
where: SSE = sum of the squared errors
SSD = sum of the squared deviations about the mean.
© 2015 - 2025 Weber Informatics LLC | Privacy Policy