Curve Fitting using Levenberg-Marquardt or similar

[Andreas Pohlmann]

Hi,

Does anyone know how I can implement a curve fit (based on Levenberg-Marquardt, or similar) in ImageJ?

*** THE PROBLEM ***
I am trying to convert a Matlab program that I wrote into an IJ plugin, so it would be usable by anyone (not just people with a Matlab license). I have to do some curve fitting of a function with several exponentials and 5 parameters (2 compartment model) - As I understand, not something the in IJ available simplex algorithm would be suitable for.

*** WHAT I FOUND SO FAR ***
I couldn't find any suitable plugin so I looked for a java implementation of LM. Wikipedia pointed me to open source Java implementations offered by
J.P. Lewis (http://www.idiom.com/~zilla/Computer/Javanumeric/index.html)
and J. Holopainen (http://users.utu.fi/jaolho).
I had a look at those but couldn't figure out how to use the code. I don't have a huge experience with IJ plugins (have written 2, usually write macros) and I am not a java expert (have average knowledge of C++).

*** HELP? ***
I was hoping that there is someone out there that has implemented a LM fit in IJ or java and could help me.
I would be also interested in any suitable alternative to LM - the fit is quite robust as all parameters control rather independent aspects of the curve (x0, y0, up slope, down slope, peak height).

If I don't find a solution I would have to ditch the project. That would be a pity as I spent 2 yrs developing the Matlab program and I would really like people to be able to use it.

Thanks a lot for your time and help!!

Andreas

__________________________________________________________________________
Erweitern Sie FreeMail zu einem noch leistungsstärkeren E-Mail-Postfach!
Mehr Infos unter http://freemail.web.de/home/landingpad/?mc=021131

[ctrueden]

Hi Andreas,

I have successfully used the LMA library you linked
(http://users.utu.fi/jaolho) to perform curve fitting in Java, but not
within ImageJ. You have to create a subclass of LMAFunction that
defines the function you want to fit (in my case, an exponential of
the form: y(t) = a*e^(-b*t) + c). Then you construct an LMA object and
call the fit() method.

Below is the relevant code from my application. It won't compile on
its own, but you could use or adapt the ExpFunction subclass I wrote,
if it is of the proper form for you. Note that I have only gotten it
working with a single exponential component -- attempts to fit
multiple exponential components thus far have resulted in failure, but
I have undoubtably made an error somewhere (any LMA experts on the
list?).

Hope it helps,
Curtis

----------

import jaolho.data.lma.LMA;
import jaolho.data.lma.LMAFunction;
import jaolho.data.lma.implementations.JAMAMatrix;

...

    double[][] fitResults = null;
    int numExp = 1; // number of exponentials to fit
    if (adjustPeaks && doRefit) {
      // perform exponential curve fitting: y(x) = a * e^(-b*t) + c
      progress.setNote("Fitting curves");
      fitResults = new double[channels][];
      ExpFunction func = new ExpFunction(numExp);
      float[] params = new float[3 * numExp];
      if (numExp == 1) {
        params[0] = maxVal;
        params[1] = 1;
        params[2] = 0;
      }
      else if (numExp == 2) {
        params[0] = maxVal / 2;
        params[1] = 0.8f;
        params[2] = 0;
        params[0] = maxVal / 2;
        params[1] = 2;
        params[2] = 0;
      }
//      for (int i=0; i<numExp; i++) {
//        // initial guess for (a, b, c)
//        int e = 3 * i;
//        params[e] = (numExp - i) * maxVal / (numExp + 1);
//        params[e + 1] = 1;
//        params[e + 2] = 0;
//      }
      int num = timeBins - maxPeak;
      float[] xVals = new float[num];
      for (int i=0; i<num; i++) xVals[i] = i;
      float[] yVals = new float[num];
      float[] weights = new float[num];
      Arrays.fill(weights, 1); // no weighting
      log("Computing fit parameters: y(t) = a * e^(-t/" + TAU + ") + c");
      for (int c=0, cc=0; c<channels; c++) {
        if (!cVisible[c]) {
          fitResults[c] = null;
          continue;
        }
        log("\tChannel #" + (c + 1) + ":");
        System.arraycopy(samps, timeBins * cc + maxPeak, yVals, 0, num);
        LMA lma = null;
        lma = new LMA(func, params, new float[][] {xVals, yVals},
          weights, new JAMAMatrix(params.length, params.length));
        lma.fit();
        log("\t\titerations=" + lma.iterationCount);
        log("\t\tchi2=" + lma.chi2);
        for (int i=0; i<numExp; i++) {
          int e = 3 * i;
          log("\t\ta" + i + "=" + lma.parameters[e]);
          log("\t\t" + TAU + i + "=" + (1 / lma.parameters[e + 1]));
          log("\t\tc" + i + "=" + lma.parameters[e + 2]);
        }
        fitResults[c] = lma.parameters;

...

  /**
   * A summed exponential function of the form:
   * y(t) = a1*e^(-b1*t) + ... + an*e^(-bn*t) + c.
   */
  public class ExpFunction extends LMAFunction {
    /** Number of exponentials to fit. */
    private int numExp = 1;

    /** Constructs a function with the given number of summed exponentials. */
    public ExpFunction(int num) { numExp = num; }

    public double getY(double x, double[] a) {
      double sum = 0;
      for (int i=0; i<numExp; i++) {
        int e = 3 * i;
        sum += a[e] * Math.exp(-a[e + 1] * x) + a[e + 2];
      }
      return sum;
    }

    public double getPartialDerivate(double x, double[] a, int parameterIndex) {
      int e = parameterIndex / 3;
      int off = parameterIndex % 3;
      switch (off) {
        case 0: return Math.exp(-a[e + 1] * x);
        case 1: return -a[e] * x * Math.exp(-a[e + 1] * x);
        case 2: return 1;
        default:
          throw new RuntimeException("No such parameter index: " +
            parameterIndex);
      }
    }
  }

On 11/9/06, Andreas Pohlmann <[hidden email]> wrote:

> Hi,
>
> Does anyone know how I can implement a curve fit (based on Levenberg-Marquardt, or similar) in ImageJ?
>
> *** THE PROBLEM ***
> I am trying to convert a Matlab program that I wrote into an IJ plugin, so it would be usable by anyone (not just people with a Matlab license). I have to do some curve fitting of a function with several exponentials and 5 parameters (2 compartment model) - As I understand, not something the in IJ available simplex algorithm would be suitable for.
>
> *** WHAT I FOUND SO FAR ***
> I couldn't find any suitable plugin so I looked for a java implementation of LM. Wikipedia pointed me to open source Java implementations offered by
> J.P. Lewis (http://www.idiom.com/~zilla/Computer/Javanumeric/index.html)
> and J. Holopainen (http://users.utu.fi/jaolho).
> I had a look at those but couldn't figure out how to use the code. I don't have a huge experience with IJ plugins (have written 2, usually write macros) and I am not a java expert (have average knowledge of C++).
>
> *** HELP? ***
> I was hoping that there is someone out there that has implemented a LM fit in IJ or java and could help me.
> I would be also interested in any suitable alternative to LM - the fit is quite robust as all parameters control rather independent aspects of the curve (x0, y0, up slope, down slope, peak height).
>
> If I don't find a solution I would have to ditch the project. That would be a pity as I spent 2 yrs developing the Matlab program and I would really like people to be able to use it.
>
> Thanks a lot for your time and help!!
>
> Andreas
>
> __________________________________________________________________________
> Erweitern Sie FreeMail zu einem noch leistungsstärkeren E-Mail-Postfach!
> Mehr Infos unter http://freemail.web.de/home/landingpad/?mc=021131
>

[Robert Dougherty]

Andreas,

I've had good luck with the optimization package from the Forest Service:
http://www1.fpl.fs.fed.us/optimization.html.  I turn my least-squares
problems into single function minimization problems (sum the squares of the
errors) and call Uncmin_f77.optif0_f77.  For some reason this has worked out
better than my previous efforts in FORTRAN using the LM code SNLS1E.  Watch
out for the 1-based indexing.

Bob

[Simon Roussel]

Hi,
QuickvolII is an ImageJ plugin.
If you look at http://www.quickvol.com/resources/QuickVolII_manual.pdf,
one can read: "Levenberg-Marquardt implementation is currently being
developed to address these problems".
.....
Cheers,
Simon