Changeset 35767

Timestamp:

Jul 3, 2013, 2:30:22 PM (13 years ago)

Author:

eugene

Message:

add new psMinimizeLMChi2 function with new convergence criterion; add threaded version of psImageSmoothNoMask

Location:

trunk/psLib

Files:

• src/imageops/psImageConvolve.c (modified) (4 diffs)
• src/imageops/psImageConvolve.h (modified) (4 diffs)
• src/math/psMinimizeLMM.c (modified) (6 diffs)
• src/math/psMinimizeLMM.h (modified) (3 diffs)
• test/imageops/Makefile.am (modified) (1 diff)
• test/imageops/tap_psImageSmoothMask_Threaded.c (copied) (copied from branches/eam_branches/ipp-20130509/psLib/test/imageops/tap_psImageSmoothMask_Threaded.c )
• test/math/Makefile.am (modified) (1 diff)
• test/math/tap_psMinimizeLMM_Alt.c (copied) (copied from branches/eam_branches/ipp-20130509/psLib/test/math/tap_psMinimizeLMM_Alt.c )
• test/math/tap_psMinimizeLMM_Trail.c (copied) (copied from branches/eam_branches/ipp-20130509/psLib/test/math/tap_psMinimizeLMM_Trail.c )
• test/pstap/src/pstap.h (modified) (1 diff)

trunk/psLib/src/imageops/psImageConvolve.c

@@ -1231 +1231 @@
 }
 
+/*********************** smooth with mask, threaded version ***********************/
+
 bool psImageSmoothMask_ScanRows_F32(psImage *calculation, psImage *calcMask, const psImage *image,
                                     const psImage *mask, psImageMaskType maskVal, psVector *gaussNorm,
@@ -1484 +1486 @@
 
 
+/*********************** smooth no-mask ***********************/
+
+bool psImageSmoothNoMask_ScanRows_F32(psImage *calculation, const psImage *image,
+                                      psVector *gaussNorm, float minGauss, 
+                                      int size, int rowStart, int rowStop)
+{
+    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
+
+    int numCols = image->numCols;
+    int xLast = numCols - 1; // Last index
+
+    for (int j = rowStart; j < rowStop; j++) {
+        for (int i = 0; i < numCols; i++) {
+            int xMin = PS_MAX(i - size, 0);
+            int xMax = PS_MIN(i + size, xLast);
+            const psF32 *imageData = &image->data.F32[j][xMin];
+            int uMin = - PS_MIN(i, size); /* Minimum kernel index */
+            const psF32 *gaussData = &gauss[uMin];
+            double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */
+            for (int x = xMin; x <= xMax; x++, imageData++, gaussData++) {
+                sumIG += *imageData * *gaussData;
+                sumG += *gaussData;
+            }
+            if (sumG > minGauss) {
+                /* BW */
+                calculation->data.F32[i][j] = sumIG / sumG;
+            } 
+        }
+    }
+    return true;
+}
+
+bool psImageSmoothNoMask_ScanCols_F32(psImage *output, psImage *calculation, 
+                                      psVector *gaussNorm, float minGauss,
+                                      int size, int colStart, int colStop)
+{
+    const psF32 *gauss = &gaussNorm->data.F32[size]; // Gaussian convolution kernel
+
+    int numRows = output->numRows;
+    int yLast = numRows - 1; // Last index
+
+    for (int i = colStart; i < colStop; i++) {
+        for (int j = 0; j < numRows; j++) {
+            int yMin = PS_MAX(j - size, 0);
+            int yMax = PS_MIN(j + size, yLast);
+            const psF32 *imageData = &calculation->data.F32[i][yMin]; /* BW */
+            int vMin = - PS_MIN(j, size); /* Minimum kernel index */
+            const psF32 *gaussData = &gauss[vMin];
+            double sumIG = 0.0, sumG = 0.0; /* Sums for convolution */
+            for (int y = yMin; y <= yMax; y++, imageData++, gaussData++) {
+                sumIG += *imageData * *gaussData;
+                sumG += *gaussData;
+            }
+            output->data.F32[j][i] = (sumG > minGauss) ? sumIG / sumG : NAN;
+        }
+    }
+    return true;
+}
+
+bool psImageSmoothNoMask_ScanRows_F32_Threaded(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert (job->args, "programming error: job args not alloced?");
+    psAssert (job->args->n == 7, "programming error: job Nargs mismatch");
+
+    psImage *calculation  = job->args->data[0]; // calculation image
+    const psImage *image  = job->args->data[1]; // input image
+
+    psVector *gaussNorm   = job->args->data[2]; // gauss kernel
+    float minGauss        = PS_SCALAR_VALUE(job->args->data[3],F32);
+    int size              = PS_SCALAR_VALUE(job->args->data[4],S32);
+    int rowStart          = PS_SCALAR_VALUE(job->args->data[5],S32);
+    int rowStop           = PS_SCALAR_VALUE(job->args->data[6],S32);
+    return psImageSmoothNoMask_ScanRows_F32(calculation, image,
+                                          gaussNorm, minGauss, size,
+                                          rowStart, rowStop);
+}
+
+bool psImageSmoothNoMask_ScanCols_F32_Threaded(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert (job->args, "programming error: job args not alloced?");
+    psAssert (job->args->n == 7, "programming error: job Nargs mismatch");
+
+    psImage *output       = job->args->data[0]; // input image
+    psImage *calculation  = job->args->data[1]; // calculation image
+
+    psVector *gaussNorm   = job->args->data[2]; // gauss kernel
+    float minGauss        = PS_SCALAR_VALUE(job->args->data[3],F32);
+    int size              = PS_SCALAR_VALUE(job->args->data[4],S32);
+    int colStart          = PS_SCALAR_VALUE(job->args->data[5],S32);
+    int colStop           = PS_SCALAR_VALUE(job->args->data[6],S32);
+    return psImageSmoothNoMask_ScanCols_F32(output, calculation,
+                                          gaussNorm, minGauss, size,
+                                          colStart, colStop);
+}
+
+psImage *psImageSmoothNoMask_Threaded(psImage *output, const psImage *image, float sigma, float numSigma, float minGauss)
+{
+    PS_ASSERT_IMAGE_NON_NULL(image, NULL);
+
+    int numCols = image->numCols, numRows = image->numRows; // Size of image
+
+    int nThreads = psThreadPoolSize();
+    int scanCols = nThreads ? (numCols / 4 / nThreads) : numCols;
+    int scanRows = nThreads ? (numRows / 4 / nThreads) : numRows;
+    int size = sigma * numSigma + 0.5;  // Half-size of kernel
+
+    // Generate normalized gaussian
+    psVector *gaussNorm = NULL;
+    IMAGE_SMOOTH_GAUSS(gaussNorm, size, sigma, F32);
+
+    switch (image->type.type) {
+        // Smooth an image with masked pixels
+        // The calculation and calcMask images are *deliberately* backwards (row,col instead of col,row or
+        // [col][row] instead of [row][col]) to optimise the iteration over rows; such instances are marked "BW"
+
+      case PS_TYPE_F32: {
+          psImage *calculation = psImageAlloc(numRows, numCols, PS_TYPE_F32); /* Calculation image; BW */
+
+          if (threaded) {
+              /** Smooth in X direction **/
+              for (int rowStart = 0; rowStart < numRows; rowStart+=scanRows) {
+                  int rowStop = PS_MIN (rowStart + scanRows, numRows);
+
+                  // allocate a job, construct the arguments for this job
+                  psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANROWS");
+                  psArrayAdd(job->args, 0, calculation);
+                  psArrayAdd(job->args, 0, (psImage *) image); // cast away const
+                  psArrayAdd(job->args, 0, gaussNorm);
+                  PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32);
+                  PS_ARRAY_ADD_SCALAR(job->args, size,     PS_TYPE_S32);
+                  PS_ARRAY_ADD_SCALAR(job->args, rowStart, PS_TYPE_S32);
+                  PS_ARRAY_ADD_SCALAR(job->args, rowStop,  PS_TYPE_S32);
+                  // -> psImageSmoothNoMask_ScanRows_F32 (calculation, image, gauss, minGauss, size, rowStart, rowStop);
+
+                  // if threading is not active, we simply run the job and return
+                  if (!psThreadJobAddPending(job)) {
+                      psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+                      psFree(calculation);
+                      psFree(gaussNorm);
+                      return false;
+                  }
+              }
+              // wait here for the threaded jobs to finish (NOP if threading is not active)
+              if (!psThreadPoolWait(true, true)) {
+                  psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+                  psFree(calculation);
+                  psFree(gaussNorm);
+                  return false;
+              }
+          } else if (!psImageSmoothNoMask_ScanRows_F32(calculation, image, gaussNorm, minGauss, size, 0, numRows)) {
+              psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+              psFree(calculation);
+              psFree(gaussNorm);
+              return false;
+          }
+
+          output = psImageRecycle(output, numCols, numRows, PS_TYPE_F32);
+
+          /** Smooth in Y direction  **/
+          if (threaded) {
+              for (int colStart = 0; colStart < numCols; colStart+=scanCols) {
+                  int colStop = PS_MIN (colStart + scanCols, numCols);
+
+                  // allocate a job, construct the arguments for this job
+                  psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANCOLS");
+                  psArrayAdd(job->args, 0, output);
+                  psArrayAdd(job->args, 0, calculation);
+                  psArrayAdd(job->args, 0, gaussNorm);
+                  PS_ARRAY_ADD_SCALAR(job->args, minGauss, PS_TYPE_F32);
+                  PS_ARRAY_ADD_SCALAR(job->args, size,     PS_TYPE_S32);
+                  PS_ARRAY_ADD_SCALAR(job->args, colStart, PS_TYPE_S32);
+                  PS_ARRAY_ADD_SCALAR(job->args, colStop,  PS_TYPE_S32);
+                  // -> psImageSmoothMask_ScanCols_F32 (output, calculation, gauss, minGauss, size, colStart, colStop);
+
+                  // if threading is not active, we simply run the job and return
+                  if (!psThreadJobAddPending(job)) {
+                      psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+                      psFree(calculation);
+                      psFree(gaussNorm);
+                      return false;
+                  }
+              }
+
+              // wait here for the threaded jobs to finish (NOP if threading is not active)
+              if (!psThreadPoolWait(true, true)) {
+                  psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+                  psFree(calculation);
+                  psFree(gaussNorm);
+                  return false;
+              }
+          } else if (!psImageSmoothNoMask_ScanCols_F32(output, calculation, gaussNorm, minGauss, size, 0, numCols)) {
+              psError(PS_ERR_UNKNOWN, false, "Unable to smooth image");
+              psFree(calculation);
+              psFree(gaussNorm);
+              return false;
+          }
+
+          psFree(calculation);
+          psFree(gaussNorm);
+
+          return output;
+      }
+      default:
+        psFree(gaussNorm);
+        char *typeStr;
+        PS_TYPE_NAME(typeStr,image->type.type);
+        psError(PS_ERR_BAD_PARAMETER_TYPE, true,
+                _("Specified psImage type, %s, is not supported."),
+                typeStr);
+        return false;
+    }
+
+    psAbort("Should never reach here.");
+    return false;
+}
+
+/****************************************************************************************/
+
 bool psImageSmoothMaskF32(psImage *image, psImage *mask, psImageMaskType maskVal,
                           double  sigma, double  Nsigma)
@@ -1900 +2122 @@
         }
         {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANROWS", 7);
+            task->function = &psImageSmoothNoMask_ScanRows_F32_Threaded;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
+        {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTH_NOMASK_SCANCOLS", 7);
+            task->function = &psImageSmoothNoMask_ScanCols_F32_Threaded;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
+        {
             psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_SMOOTHMASK_PIXELS", 9);
             task->function = &psImageSmoothMaskPixelsThread;
@@ -1909 +2143 @@
         psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_SCANROWS");
         psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_SCANCOLS");
+        psThreadTaskRemove("PSLIB_IMAGE_SMOOTH_NOMASK_SCANROWS");
+        psThreadTaskRemove("PSLIB_IMAGE_SMOOTH_NOMASK_SCANCOLS");
         psThreadTaskRemove("PSLIB_IMAGE_SMOOTHMASK_PIXELS");
     }

trunk/psLib/src/imageops/psImageConvolve.h

@@ -242 +242 @@
     );
 
-/// Smooth an image by parts using 1D Gaussian independently in x and y, allowing for masked pixels
+/// Smooth an imageby parts using 1D Gaussian independently in x and y, allowing for
+/// MASKED PIXELS
 ///
 /// Applies a circularly symmetric Gaussian smoothing first in x and then in y
@@ -256 +257 @@
     );
 
-/// Smooth particular pixels on an image, allowing for masked pixels
+/// Smooth particular pixels on an image, allowing for MASKED PIXELS
 ///
 /// Applies a circularly symmetric Gaussian smoothing first in x and then in y
@@ -271 +272 @@
     );
 
+/// Smooth an image by parts using 1D Gaussian independently in x and y, allowing for
+/// MASKED PIXELS : THREADED VERSION
+///
+/// Applies a circularly symmetric Gaussian smoothing first in x and then in y
+/// directions with just a vector.  This process is 2N faster than 2D convolutions (in general).
 psImage *psImageSmoothMask_Threaded(psImage *output,
                                     const psImage *image,
@@ -279 +285 @@
                                     float minGauss);
 
+/// Smooth an image by parts using 1D Gaussian independently in x and y, allowing for
+/// MASKED PIXELS : THREADED VERSION
+///
+/// Applies a circularly symmetric Gaussian smoothing first in x and then in y
+/// directions with just a vector.  This process is 2N faster than 2D convolutions (in general).
+psImage *psImageSmoothNoMask_Threaded(psImage *output,
+                                    const psImage *image,
+                                    float sigma,
+                                    float numSigma,
+                                    float minGauss);
+
+/// Smooth an image by parts IN-SITU using 1D Gaussian independently in x and y, allowing
+/// for MASKED PIXELS
 bool psImageSmoothMaskF32(
     psImage *image,                    ///< the image to be smoothed

trunk/psLib/src/math/psMinimizeLMM.c

@@ -325 +325 @@
         chisq += PS_SQR(delta) * dy->data.F32[i];
 
+        // XXX remove this later:
+        // psVector *tmp = x->data[i];
+        // fprintf (stderr, "%f %f  %f %f  %f\n", tmp->data.F32[0], tmp->data.F32[1], y->data.F32[i], dy->data.F32[i], ymodel);
+
         if (isnan(dy->data.F32[i])) goto escape;
         if (isnan(delta)) goto escape;
@@ -360 +364 @@
 }
 
+/******************************************************************************
+psMinimizeLMChi2():  wrapper to call either _Old or _Alt
+  *****************************************************************************/
+bool psMinimizeLMChi2(
+    psMinimization *min,
+    psImage *covar,
+    psVector *params,
+    psMinConstraint *constraint,
+    const psArray *x,
+    const psVector *y,
+    const psVector *yWt,
+    psMinimizeLMChi2Func func)
+{
+    bool status = psMinimizeLMChi2_Alt(
+        min,
+        covar,
+        params,
+        constraint,
+        x,
+        y,
+        yWt,
+        func);
+    return status;
+}
 
 /******************************************************************************
@@ -370 +398 @@
 XXX Make an F64 version?
   *****************************************************************************/
-bool psMinimizeLMChi2(
+bool psMinimizeLMChi2_Old(
     psMinimization *min,
     psImage *covar,
@@ -507 +535 @@
         psF32 rho = (min->value - Chisq) / dLinear;
 
-        psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %f, nDOF: %d\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda, nDOF);
-
-        psTrace("psLib.math.dLinear", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %f\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda);
+        psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %g, nDOF: %d\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda, nDOF);
+
+        psTrace("psLib.math.dLinear", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %g\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda);
 
         // dump some useful info if trace is defined
@@ -580 +608 @@
 }
 
+/******************************************************************************
+psMinimizeLMChi2(): This routine takes a function-pointer (func) which calculates an arbitrary
+function and it's derivatives and minimizes the chi-squared match between that function at the
+specified points and the specified value at those points.
+
+The original version of this function used a convergence criterion based on the change in
+chisq.  this has problems since it depends on the choice of points used to measure the fit.
+(consider a gaussian on a background : it 100 pixels are used -- and some or most contribute to
+the chisq -- and the delta-chisq is 10%, then the same change in model fit will yield a
+delta-chisq of 1% if 1000 pixels are used (all but the 100 measuring the background)).
+
+This implementation uses changes to the parameters and stops if they are no longer significant.
+
+This requires F32 input data; all internal calls use F32.
+  *****************************************************************************/
+bool psMinimizeLMChi2_Alt(
+    psMinimization *min,
+    psImage *covar,
+    psVector *params,
+    psMinConstraint *constraint,
+    const psArray *x,
+    const psVector *y,
+    const psVector *yWt,
+    psMinimizeLMChi2Func func)
+{
+    psTrace("psLib.math", 3, "---- begin ----\n");
+    PS_ASSERT_PTR_NON_NULL(min, false);
+    PS_ASSERT_VECTOR_NON_NULL(params, false);
+    PS_ASSERT_VECTOR_NON_EMPTY(params, false);
+    PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false);
+    psVector *paramMask = NULL;
+    if (constraint != NULL) {
+        paramMask = constraint->paramMask;
+        if (paramMask != NULL) {
+            PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, false);
+            PS_ASSERT_VECTORS_SIZE_EQUAL(params, paramMask, false);
+        }
+    }
+    PS_ASSERT_PTR_NON_NULL(x, false);
+    for (psS32 i = 0 ; i < x->n ; i++) {
+        psVector *coord = (psVector *) (x->data[i]);
+        PS_ASSERT_VECTOR_NON_NULL(coord, false);
+        PS_ASSERT_VECTOR_TYPE(coord, PS_TYPE_F32, false);
+    }
+    PS_ASSERT_VECTOR_NON_NULL(y, false);
+    PS_ASSERT_VECTOR_NON_EMPTY(y, false);
+    PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, false);
+    PS_ASSERT_VECTORS_SIZE_EQUAL(x, y, false);
+    if (yWt != NULL) {
+        PS_ASSERT_VECTOR_TYPE(yWt, PS_TYPE_F32, false);
+        PS_ASSERT_VECTORS_SIZE_EQUAL(y, yWt, false);
+    }
+    PS_ASSERT_PTR_NON_NULL(func, false);
+
+    psMinimizeLMLimitFunc checkLimits = NULL;
+    if (constraint) {
+        checkLimits = constraint->checkLimits;
+    }
+
+    // this function has test and current values for several things
+    // the current best value is in lower case
+    // the next guess value is in upper case
+
+    // allocate internal arrays (current vs Guess)
+    psImage *Alpha = NULL;
+    psVector *Beta = NULL;
+
+    // Alpha & Beta only contain elements to represent the unmasked parameters
+    if (!psMinLM_AllocAB (&Alpha, &Beta, params, paramMask)) {
+        psAbort ("programming error: no unmasked parameters to be fit\n");
+    }
+
+    int nFitParams = Beta->n;
+    psImage *alpha   = psImageAlloc(nFitParams, nFitParams, PS_TYPE_F32);
+    psVector *beta   = psVectorAlloc(nFitParams, PS_TYPE_F32);
+    psVector *Params = psVectorAlloc(params->n, PS_TYPE_F32);
+
+    psVector *dy     = NULL;
+    psF32 Chisq = 0.0;
+    psF32 lambda = 0.001;
+    psF32 dLinear = 0.0;
+    psF32 nu = 2.0;
+
+    // the user provides the error or NULL.  we need to convert
+    // to appropriate weights
+    if (yWt != NULL) {
+        dy = (psVector *) yWt;
+    } else {
+        dy = psVectorAlloc(y->n, PS_TYPE_F32);
+        psVectorInit(dy, 1.0);
+    }
+
+    // number of degrees of freedom for this fit
+    int nDOF = dy->n - nFitParams;
+
+    // calculate initial alpha and beta, set chisq (min->value)
+    min->value = psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
+    if (isnan(min->value)) {
+        min->iter = min->maxIter;
+        psFree(alpha);
+        psFree(Alpha);
+        psFree(beta);
+        psFree(Beta);
+        psFree(Params);
+        return(false);
+    }
+    // dump some useful info if trace is defined
+    if (psTraceGetLevel("psLib.math") >= 6) {
+        p_psImagePrint(psTraceGetDestination(), alpha, "alpha guess (0)");
+        p_psVectorPrint(psTraceGetDestination(), beta, "beta guess (0)");
+    }
+    if (psTraceGetLevel("psLib.math") >= 5) {
+        p_psVectorPrint(psTraceGetDestination(), params, "params guess (0)");
+    }
+
+    bool done = (min->iter >= min->maxIter);
+    while (!done) {
+        psTrace("psLib.math", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
+
+        if (min->chisqConvergence) {
+          psTrace("psLib.math", 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol);
+        } else {
+          psTrace("psLib.math", 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->rParSigma, min->minTol*nFitParams, min->maxTol*nFitParams);
+        }
+
+        // set a new guess for Alpha, Beta, Params
+        if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) {
+            min->iter ++;
+            if (min->iter >=  min->maxIter) break;
+            lambda *= 10.0;
+            // ALT? lambda *= 2.0;
+            continue;
+        }
+
+        // dump some useful info if trace is defined
+        if (psTraceGetLevel("psLib.math") >= 6) {
+            p_psImagePrint(psTraceGetDestination(), Alpha, "Alpha guess (1)");
+            p_psVectorPrint(psTraceGetDestination(), Beta, "Beta guess (1)");
+            p_psVectorPrint(psTraceGetDestination(), beta, "beta current (1)");
+        }
+        if (psTraceGetLevel("psLib.math") >= 5) {
+            p_psVectorPrint(psTraceGetDestination(), Params, "params guess (1)");
+        }
+
+        // calculate the parameter change (rParDelta) and error radius (rParSigma)
+        //    rParDelta : radius of parameter change;
+        //    rParSigma : radius of parameter error 
+        
+        // note that (before SetABX) Alpha[i][i] is the covariance matrix and
+        // Beta is the actual parameter change for this pass
+
+        // note that Alpha & Beta only represent unmasked parameters, while params and Params have all 
+
+        // dParSigma = Alpha[i][i] : error (squared) on parameter i
+        // dParDelta = Params->data.F32[i] - params->data.F32[i]     : change on parameter i
+        float rParSigma = 0.0;
+        for (int j = 0, J = 0; j < Params->n; j++) {
+            if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) {
+                continue;
+            }
+            rParSigma += PS_SQR(Params->data.F32[j] - params->data.F32[j]) / Alpha->data.F32[J][J];
+            J++;
+        }
+        rParSigma = sqrt(rParSigma);
+        psTrace("psLib.math", 5, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter);
+        // fprintf (stderr, "rParSigma: %f, Niter: %d\n", rParSigma, min->iter);
+
+        // calculate Chisq for new guess, update Alpha & Beta
+        Chisq = psMinLM_SetABX(Alpha, Beta, Params, paramMask, x, y, dy, func);
+        if (isnan(Chisq)) {
+            min->iter ++;
+            if (min->iter >= min->maxIter) break;
+            lambda *= 10.0;
+            // ALT lambda *= 2.0;
+            continue;
+        }
+
+        // convergence criterion:
+        // compare the delta (min->value - Chisq) with the
+        // expected delta from the linear model (dLinear)
+        // accept new guess if it is an improvement (rho > 0), or else increase lambda
+        psF32 rho = (min->value - Chisq) / dLinear;
+
+        psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %g, nDOF: %d\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda, nDOF);
+
+        psTrace("psLib.math.dLinear", 5, "last chisq: %f, new chisq %f, delta: %f, dLinear: %f, rho: %f, lambda: %g\n", min->value, Chisq, min->lastDelta, dLinear, rho, lambda);
+
+        // dump some useful info if trace is defined
+        if (psTraceGetLevel("psLib.math") >= 6) {
+            p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (2)");
+            p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (2)");
+        }
+
+        // change in chisq/nDOF since last minimum
+        min->lastDelta = (min->value - Chisq) / nDOF;
+
+        // rho is positive if the new chisq is smaller; allow for some insignificant change (slight negative rho)
+
+        // XXX the old version of lambda changes:
+        // XXX : Madsen gives suggestion for better use of rho
+        // rho is positive if the new chisq is smaller
+        if (rho >= -1e-6) {
+            min->value = Chisq;
+            alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F32);
+            beta   = psVectorCopy(beta, Beta, PS_TYPE_F32);
+            params = psVectorCopy(params, Params, PS_TYPE_F32);
+        } 
+        switch (min->gainFactorMode) {
+          case 0:
+            if (rho >= -1e-6) {
+              lambda *= 0.25;
+            } else {
+              lambda *= 10.0;
+            }
+            break;
+
+          case 1:
+            // adjust the gain ratio (lambda) based on rho
+            if (rho < 0.25) {
+              lambda *= 2.0;
+            } 
+            if (rho > 0.75) {
+              lambda *= 0.333;
+            }
+            break;
+
+          case 2:
+            if (rho > 0.0) {
+              lambda *= PS_MAX(0.33, (1.0 - pow(2.0*rho - 1.0, 3.0)));
+              nu = 2.0;
+            } else {
+              lambda *= nu;
+              nu *= 2.0;
+            }
+            break;
+        }
+        min->iter++;
+
+        // ending conditions:
+        // 1) hard limit : too many iterations
+        done = (min->iter >= min->maxIter);
+        
+        // 2) require deltaChi > 1e-6 (ie, chisq is decreasing, but accept an insignificant change)
+        if (min->lastDelta < -1e-6) {
+            continue;
+        }
+
+        // save this value in case we stop iterating
+        min->rParSigma = rParSigma;
+
+        // 2) require chisqDOF < maxChisqDOF (if maxChisqDOF is not NAN)
+        // keep iterating regardless of rParSigma in this case
+        float chisqDOF = Chisq / nDOF;
+        if (isfinite(min->maxChisqDOF) && (chisqDOF > min->maxChisqDOF)) {
+            continue;
+        }
+
+        // delta-chisq or rParSigma ?
+        if (min->chisqConvergence) {
+          done |= (min->lastDelta < min->minTol);
+        } else {
+          done |= (rParSigma < min->minTol*nFitParams);
+        }
+    }
+    psTrace("psLib.math", 5, "chisq: %f, last delta: %f, rParSigma: %f, Niter: %d\n", min->value, min->lastDelta, min->rParSigma, min->iter);
+
+    // construct & return the covariance matrix (if requested)
+    if (covar != NULL) {
+        if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {
+            psTrace ("psLib.math", 5, "failure to calculate covariance matrix\n");
+        }
+        // set covar values which are not masked
+        psImageInit (covar, 0.0);
+        for (int j = 0, J = 0; j < params->n; j++) {
+            if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[j])) {
+                covar->data.F32[j][j] = 1.0;
+                continue;
+            }
+            for (int k = 0, K = 0; k < params->n; k++) {
+                if (paramMask && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[k])) continue;
+                covar->data.F32[j][k] = Alpha->data.F32[J][K];
+                K++;
+            }
+            J++;
+        }
+    }
+
+    // free the internal temporary data
+    psFree(alpha);
+    psFree(Alpha);
+    psFree(beta);
+    psFree(Beta);
+    psFree(Params);
+    if (yWt == NULL) {
+        psFree(dy);
+    }
+
+    // if the last improvement was at least as good as maxTol, accept the fit:
+    if (min->chisqConvergence) {
+      if (min->lastDelta <= min->maxTol) {
+        psTrace("psLib.math", 6, "---- end (true) ----\n");
+        return(true);
+      }
+    } else {
+      if (min->rParSigma <= min->maxTol*nFitParams) {
+        psTrace("psLib.math", 6, "---- end (true) ----\n");
+        return(true);
+      }
+    }
+    psTrace("psLib.math", 6, "---- end (false) ----\n");
+    return(false);
+}
+
 bool psMinLM_AllocAB (psImage **Alpha, psVector **Beta, const psVector *params, const psVector *paramMask) {
 
@@ -625 +966 @@
     min->iter = 0;
     min->lastDelta = NAN;
+    min->rParSigma = NAN;
     min->maxChisqDOF = NAN;
 
+    // we default to the old algorithm for convergence
+    min->chisqConvergence = true;
+    min->gainFactorMode = 0;
+    
     return(min);
 }

trunk/psLib/src/math/psMinimizeLMM.h

@@ -30 +30 @@
 #include "psConstants.h"
 
+# define PS_MINIMIZE_LMM_GAIN_FACTOR_MODE 0
+# define PS_MINIMIZE_LMM_CHISQ_CONVERGENCE 1
+
 #define PS_DETERMINE_BRACKET_STEP_SIZE 0.10
 #define PS_MAX_LMM_ITERATIONS 100
@@ -82 +85 @@
     int iter;                          ///< Number of iterations to date
     float lastDelta;                   ///< The last difference for the fit
+    float rParSigma;                   ///< last fractional change in the parameters
     float maxChisqDOF;                 ///< for Chisq minimization, require that we reach here before checking tolerance
+    int gainFactorMode;
+    bool chisqConvergence;
 }
 psMinimization;
@@ -124 +130 @@
  */
 bool psMinimizeLMChi2(
+    psMinimization *min,               ///< Minimization specification
+    psImage *covar,                    ///< Covariance matrix
+    psVector *params,                  ///< "Best Guess" for the parameters that minimize func
+    psMinConstraint *constraint, ///< Constraints on the parameters
+    const psArray *x,                  ///< Measurement ordinates of multiple vectors
+    const psVector *y,                 ///< Measurement coordinates
+    const psVector *yWt,               ///< Errors in the measurement coordinates
+    psMinimizeLMChi2Func func          ///< Specified function
+);
+
+/** Minimizes a specified function based on the Levenberg-Marquardt method.
+ *
+ *  @return bool:   True if successful.
+ */
+bool psMinimizeLMChi2_Old(
+    psMinimization *min,               ///< Minimization specification
+    psImage *covar,                    ///< Covariance matrix
+    psVector *params,                  ///< "Best Guess" for the parameters that minimize func
+    psMinConstraint *constraint, ///< Constraints on the parameters
+    const psArray *x,                  ///< Measurement ordinates of multiple vectors
+    const psVector *y,                 ///< Measurement coordinates
+    const psVector *yWt,               ///< Errors in the measurement coordinates
+    psMinimizeLMChi2Func func          ///< Specified function
+);
+
+/** Minimizes a specified function based on the Levenberg-Marquardt method
+    Uses alternative convergence criterion.
+ *
+ *  @return bool:   True if successful.
+ */
+bool psMinimizeLMChi2_Alt(
     psMinimization *min,               ///< Minimization specification
     psImage *covar,                    ///< Covariance matrix

trunk/psLib/test/imageops/Makefile.am

@@ -17 +17 @@
         tap_psImagePixelManip \
         tap_psImageSmooth \
+        tap_psImageSmoothMask_Threaded \
         tap_psImageSmooth_PreAlloc \
         tap_psImageStructManip \

trunk/psLib/test/math/Makefile.am

@@ -48 +48 @@
         tap_psMatrixVectorArithmetic04 \
         tap_psMinimizePowell \
+        tap_psMinimizeLMM \
+        tap_psMinimizeLMM_Alt \
+        tap_psMinimizeLMM_Trail \
         tap_psSpline1D \
         tap_psPolynomialMD \

trunk/psLib/test/pstap/src/pstap.h

@@ -144 +144 @@
     } \
 }
+
+# define ok_float      is_float     
+# define ok_float_tol  is_float_tol 
+# define ok_double     is_double    
+# define ok_double_tol is_double_tol
+# define ok_str        is_str         
+# define ok_strn       is_strn      
+# define ok_int        is_int         
+# define ok_long       is_long      
+# define ok_bool       is_bool