Changeset 28998

Timestamp:

Aug 20, 2010, 11:47:48 AM (16 years ago)

Author:

eugene

Message:

added pre-ffted kernel convolutions, changed minimization to accept maxTol and minTol

Location:

trunk/psLib

Files:

• src/fft/psImageFFT.c (modified) (1 diff)
• src/fft/psImageFFT.h (modified) (2 diffs)
• src/jpeg/psImageJpeg.c (modified) (2 diffs)
• src/math/psMinimizeLMM.c (modified) (9 diffs)
• src/math/psMinimizeLMM.h (modified) (3 diffs)
• src/math/psMinimizePowell.c (modified) (4 diffs)
• src/math/psStats.c (modified) (1 diff)
• test/math/Makefile.am (modified) (1 diff)
• test/math/data/polyMD.dat (copied) (copied from branches/eam_branches/ipp-20100621/psLib/test/math/data/polyMD.dat )
• test/math/data/polyMD.pro (copied) (copied from branches/eam_branches/ipp-20100621/psLib/test/math/data/polyMD.pro )
• test/math/tap_psPolynomialMD_sampleDark.c (copied) (copied from branches/eam_branches/ipp-20100621/psLib/test/math/tap_psPolynomialMD_sampleDark.c )

trunk/psLib/src/fft/psImageFFT.c

@@ -532 +532 @@
     return out;
 }
+
+void psKernelFFTFree(psKernelFFT *kernel) {
+
+    if (!kernel->fft) return;
+    
+    bool threaded = psThreadPoolSize(); // Are we running threaded?
+
+    FFTW_LOCK;
+    fftwf_free(kernel->fft);
+    FFTW_UNLOCK;
+
+    return;
+}
+
+psKernelFFT *psKernelFFTAlloc(const psKernel *input) {
+
+    psKernelFFT *kernel = psAlloc(sizeof(psKernelFFT));
+    psMemSetDeallocator(kernel, (psFreeFunc)psKernelFFTFree);
+
+    kernel->fft = NULL;
+
+    kernel->xMin = input->xMin;
+    kernel->xMax = input->xMax;
+    kernel->yMin = input->yMin;
+    kernel->yMax = input->yMax;
+
+    kernel->paddedCols = 0;
+    kernel->paddedRows = 0;
+
+    return kernel;
+}
+
+// generate the FFTed kernel image appropriate to be used for convolution of the given input image
+psKernelFFT *psImageConvolveKernelInit(const psImage *in, const psKernel *kernel)
+{
+    PS_ASSERT_IMAGE_NON_NULL(in, NULL);
+    PS_ASSERT_IMAGE_TYPE(in, PS_TYPE_F32, NULL);
+    PS_ASSERT_KERNEL_NON_NULL(kernel, NULL);
+
+    bool threaded = psThreadPoolSize(); // Are we running threaded?
+
+    int numCols = in->numCols, numRows = in->numRows; // Size of image
+    int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Kernel sizes
+
+    // Need to pad the kernel (and later the input image) to protect from wrap-around effects
+    if (xMax - xMin > numCols || yMax - yMin > numRows) {
+        // Cannot pad the image if the kernel is larger.
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true,
+                _("Kernel cannot extend further than input image size (%dx%d vs %dx%d)."),
+                xMax, yMax, numCols, numRows);
+        return NULL;
+    }
+
+    int paddedCols = numCols + PS_MAX(-xMin, xMax); // Number of columns in padded image
+    int paddedRows = numRows + PS_MAX(-yMin, yMax); // Number of rows in padded image
+
+#if CONVOLVE_FFT_BINARY_SIZE
+    // Make the size an integer power of two
+    {
+        int twoCols, twoRows;           // Size that is a factor of two
+        for (twoCols = 1; twoCols <= paddedCols && twoCols < INT_MAX - 1; twoCols <<= 1); // No action
+        for (twoRows = 1; twoRows <= paddedRows && twoRows < INT_MAX - 1; twoRows <<= 1); // No action
+        if (paddedCols > twoCols || paddedRows > twoRows) {
+            psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Unable to scale size (%dx%d) up to factor of two",
+                    paddedCols, paddedRows);
+            return NULL;
+        }
+        paddedCols = twoCols;
+        paddedRows = twoRows;
+    }
+#endif
+
+    int numPadded = paddedCols * paddedRows; // Number of pixels in padded image
+
+    // Create data array containing the padded kernel
+    FFTW_LOCK;
+    psF32 *data = fftwf_malloc(numPadded * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Data for FFTW
+    FFTW_UNLOCK;
+
+    // Generate the padded kernel image.  We could generate the padded kernel image using
+    // memcpy, but by going pixel by pixel we can apply the normalisation that corrects for the
+    // FFT renormalisation.  By applying it to the kernel here, we save applying it to the
+    // output image(s).
+
+    // copy kernel into data array
+    psF32 *dataPtr = data;              // Pointer into FFTW data
+    float norm = 1.0 / (float)(paddedRows * paddedCols); // Normalisation to correct for FFT
+
+    int xNegMin = PS_MIN(-1, xMin), xNegMax = PS_MIN(-1, xMax); // Min and max for x when negative
+    int xPosMin = PS_MAX(0, xMin), xPosMax = PS_MAX(0, xMax); // Min and max for x when positive
+    int yNegMin = PS_MIN(-1, yMin), yNegMax = PS_MIN(-1, yMax); // Min and max for x when negative
+    int yPosMin = PS_MAX(0, yMin), yPosMax = PS_MAX(0, yMax); // Min and max for x when positive
+    int blankCols = xNegMin + paddedCols - xPosMax - 1; // Number of columns between kernel extrema
+    int blankRows = (yNegMin + paddedRows - yPosMax - 1) * paddedCols; // Rows between kernel extrema
+    size_t blankColBytes = blankCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes in blankCols
+    for (int y = yPosMin; y <= yPosMax; y++) {
+        // y is positive
+        for (int x = xPosMin; x <= xPosMax; x++, dataPtr++) {
+            // x is positive
+            *dataPtr = kernel->kernel[y][x] * norm;
+        }
+        // Columns between kernel extrema
+        memset(dataPtr, 0, blankColBytes);
+        dataPtr += blankCols;
+        for (int x = xNegMin; x <= xNegMax; x++, dataPtr++) {
+            // x is negative
+            *dataPtr = kernel->kernel[y][x] * norm;
+        }
+    }
+    // Rows between kernel extrema
+    memset(dataPtr, 0, blankRows * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+    dataPtr += blankRows;
+    for (int y = yNegMin; y <= yNegMax; y++) {
+        // y is negative
+        for (int x = xPosMin; x <= xPosMax; x++, dataPtr++) {
+            // x is positive
+            *dataPtr = kernel->kernel[y][x] * norm;
+        }
+        // Columns between kernel extrema
+        memset(dataPtr, 0, blankColBytes);
+        dataPtr += blankCols;
+        for (int x = xNegMin; x <= xNegMax; x++, dataPtr++) {
+            // x is negative
+            *dataPtr = kernel->kernel[y][x] * norm;
+        }
+    }
+
+#if 0
+    {
+        // Use this for inspecting the result of copying the kernel
+        psImage *test = psImageAlloc(paddedCols, paddedRows, PS_TYPE_F32);
+        psFree(test->p_rawDataBuffer);
+        test->p_rawDataBuffer = data;
+        test->data.V[0] = test->p_rawDataBuffer;
+        for (int y = 1; y < paddedRows; y++) {
+            test->data.V[y] = (psPtr)((int8_t *)test->data.V[y - 1] + paddedCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+        }
+        // View image here or write to disk
+        test->p_rawDataBuffer = NULL;
+        psFree(test);
+    }
+#endif
+
+    // Do the forward FFT
+    // Note that the FFT images have different size from the input
+    FFTW_LOCK;
+    fftwf_complex *fft = fftwf_malloc((paddedCols/2 + 1) * paddedRows * sizeof(fftwf_complex)); // FFT
+    FFTW_UNLOCK;
+
+    FFTW_LOCK;
+    fftwf_plan plan = fftwf_plan_dft_r2c_2d(paddedRows, paddedCols, data, fft, FFTW_PLAN_RIGOR);
+    FFTW_UNLOCK;
+
+    fftwf_execute(plan);
+
+    FFTW_LOCK;
+    fftwf_destroy_plan(plan);
+    fftwf_free(data);
+    FFTW_UNLOCK;
+
+    psKernelFFT *output = psKernelFFTAlloc(kernel);
+    output->fft = fft;
+    output->paddedCols = paddedCols;
+    output->paddedRows = paddedRows;
+
+    return output;
+}
+
+psImage *psImageConvolveKernel(psImage *out, const psImage *in, const psImage *mask, psImageMaskType maskVal, const psKernelFFT *kernel)
+{
+    PS_ASSERT_IMAGE_NON_NULL(in, NULL);
+    PS_ASSERT_IMAGE_TYPE(in, PS_TYPE_F32, NULL);
+    PS_ASSERT_PTR_NON_NULL(kernel, NULL);
+    PS_ASSERT_PTR_NON_NULL(kernel->fft, NULL);
+    if (mask) {
+        PS_ASSERT_IMAGE_NON_NULL(mask, NULL);
+        PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_IMAGE_MASK, NULL);
+        PS_ASSERT_IMAGES_SIZE_EQUAL(mask, in, NULL);
+    }
+
+    bool threaded = psThreadPoolSize(); // Are we running threaded?
+
+    int numCols = in->numCols, numRows = in->numRows; // Size of image
+    int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Kernel sizes
+
+    // Need to pad the input image to protect from wrap-around effects
+    if (xMax - xMin > numCols || yMax - yMin > numRows) {
+        // Cannot pad the image if the kernel is larger.
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true, _("Kernel cannot extend further than input image size (%dx%d vs %dx%d)."), xMax, yMax, numCols, numRows);
+        return NULL;
+    }
+
+    int paddedCols = numCols + PS_MAX(-xMin, xMax); // Number of columns in padded image
+    int paddedRows = numRows + PS_MAX(-yMin, yMax); // Number of rows in padded image
+
+#if CONVOLVE_FFT_BINARY_SIZE
+    // Make the size an integer power of two
+    {
+        int twoCols, twoRows;           // Size that is a factor of two
+        for (twoCols = 1; twoCols <= paddedCols && twoCols < INT_MAX - 1; twoCols <<= 1); // No action
+        for (twoRows = 1; twoRows <= paddedRows && twoRows < INT_MAX - 1; twoRows <<= 1); // No action
+        if (paddedCols > twoCols || paddedRows > twoRows) {
+            psError(PS_ERR_BAD_PARAMETER_SIZE, true, "Unable to scale size (%dx%d) up to factor of two",
+                    paddedCols, paddedRows);
+            return NULL;
+        }
+        paddedCols = twoCols;
+        paddedRows = twoRows;
+    }
+#endif
+
+    if (paddedCols != kernel->paddedCols) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true, _("Image is inconsistent with allocated kernel"));
+        return NULL;
+    }   
+    if (paddedRows != kernel->paddedRows) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true, _("Image is inconsistent with allocated kernel"));
+        return NULL;
+    }   
+
+    int numPadded = paddedCols * paddedRows; // Number of pixels in padded image
+
+    // Create data array containing the padded image
+    FFTW_LOCK;
+    psF32 *data = fftwf_malloc(numPadded * PSELEMTYPE_SIZEOF(PS_TYPE_F32)); // Data for FFTW
+    FFTW_UNLOCK;
+    psF32 *dataPtr = data;              // Pointer into FFTW data
+    psF32 **imageData = in->data.F32;   // Pointer into image data
+
+    // Copy image into data array
+    size_t goodBytes = numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes per image row
+    size_t padBytes = (paddedCols - numCols) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to pad
+    for (int y = 0; y < numRows; y++, dataPtr += paddedCols, imageData++) {
+        memcpy(dataPtr, *imageData, goodBytes);
+        memset(dataPtr + numCols, 0, padBytes);
+    }
+    memset(dataPtr, 0, (paddedRows - numRows) * paddedCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+
+#if 0
+    {
+        // Use this for inspecting the result of copying the image
+        psImage *test = psImageAlloc(paddedCols, paddedRows, PS_TYPE_F32);
+        psFree(test->p_rawDataBuffer);
+        test->p_rawDataBuffer = data;
+        test->data.V[0] = test->p_rawDataBuffer;
+        for (int y = 1; y < paddedRows; y++) {
+            test->data.V[y] = (psPtr)((int8_t *)test->data.V[y - 1] +
+                                      paddedCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32));
+        }
+        // View image here
+        test->p_rawDataBuffer = NULL;
+        psFree(test);
+    }
+#endif
+
+    // Mask bad pixels (which may be NANs), lest they infect everything
+    if (mask && maskVal) {
+        for (int y = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++) {
+                if (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal) {
+                    data[x + paddedCols * y] = 0;
+                }
+            }
+        }
+    }
+
+    // Do the forward FFT
+    // Note that the FFT images have different size from the input
+    FFTW_LOCK;
+    fftwf_complex *fft = fftwf_malloc((paddedCols/2 + 1) * paddedRows * sizeof(fftwf_complex)); // FFT
+    FFTW_UNLOCK;
+
+    int fftCols = paddedCols/2 + 1, fftRows = paddedRows; // Size of FFT images
+    int fftPixels = fftCols * fftRows;  // Number of pixels in FFT image
+
+    FFTW_LOCK;
+    fftwf_plan forward = fftwf_plan_dft_r2c_2d(paddedRows, paddedCols, data, fft, FFTW_PLAN_RIGOR);
+    FFTW_UNLOCK;
+
+    fftwf_execute(forward);
+
+    FFTW_LOCK;
+    fftwf_destroy_plan(forward);
+    FFTW_UNLOCK;
+
+    fftwf_complex *kernelFFT = kernel->fft;
+
+    // Multiply the two transforms
+    for (int i = 0; i < fftPixels; i++) {
+        // (a + bi) * (c + di) = (ac - bd) + (bc + ad)i
+#if !defined(FFTW_NO_Complex) && defined(_Complex_I) && defined(complex) && defined(I)
+        // C99 complex support
+        fft[i] *= kernelFFT[i];
+#else
+        // FFTW's backup complex support
+        float imageReal = fft[i][0], imageImag = fft[i][1];
+        float kernelReal = kernelFFT[i][0], kernelImag = kernelFFT[i][1];
+        fft[i][0] = imageReal * kernelReal - imageImag * kernelImag;
+        fft[i][1] = imageImag * kernelReal + imageReal * kernelImag;
+#endif
+    }
+
+    // Do the backward FFT
+    FFTW_LOCK;
+    fftwf_plan backward = fftwf_plan_dft_c2r_2d(paddedRows, paddedCols, fft, data, FFTW_PLAN_RIGOR);
+    FFTW_UNLOCK;
+
+    fftwf_execute(backward);
+
+    FFTW_LOCK;
+    fftwf_destroy_plan(backward);
+    fftwf_free(fft);
+    FFTW_UNLOCK;
+
+    // Copy into the target, without the padding
+    out = psImageRecycle(out, numCols, numRows, PS_TYPE_F32);
+    psF32 **outData = out->data.F32;    // Pointer into output
+    dataPtr = data;                     // Reset to start
+    for (int y = 0; y < numRows; y++, outData++, dataPtr += paddedCols) {
+        memcpy(*outData, dataPtr, goodBytes);
+    }
+
+    FFTW_LOCK;
+    fftwf_free(data);
+    FFTW_UNLOCK;
+
+    return out;
+}

trunk/psLib/src/fft/psImageFFT.h

@@ -15 +15 @@
 #include "psImage.h"
 #include "psImageConvolve.h"
+
+typedef struct {
+    void *fft;
+    int xMin;                          ///< Most negative x index
+    int yMin;                          ///< Most negative y index
+    int xMax;                          ///< Most positive x index
+    int yMax;                          ///< Most positive y index
+    int paddedCols;
+    int paddedRows;
+} psKernelFFT;
 
 /// @addtogroup MathOps Mathematical Operations
@@ -71 +81 @@
     const psImage *mask,                ///< Corresponding mask
     psImageMaskType maskVal,            ///< Value to mask
-    const psKernel *kernel              ///< kernel to colvolve with
+    const psKernel *kernel              ///< kernel to convolve with
+    );
+
+/// Allocate an the psKernelFFT structure
+psKernelFFT *psKernelFFTAlloc(
+    const psKernel *input               ///< kernel to convolve with
 );
+
+/// Generate an FFT'ed kernel suitable for convolution with the given image
+psKernelFFT *psImageConvolveKernelInit(
+    const psImage *in,                  ///< representative image to convolve
+    const psKernel *kernel              ///< kernel to convolve with
+    );
+
+/// Convolve the given image with the pre-FFT'ed kernel
+psImage *psImageConvolveKernel(
+    psImage *out,                       ///< Output image, or NULL
+    const psImage *in,                  ///< Image to convolve
+    const psImage *mask,                ///< Corresponding mask
+    psImageMaskType maskVal,            ///< Value to mask
+    const psKernelFFT *kernel           ///< kernel to convolve with
+    );
 
 /// @}

trunk/psLib/src/jpeg/psImageJpeg.c

@@ -7 +7 @@
 #include <string.h>
 
+#include <kapa.h>
 #include "psMemory.h"
 #include "psImage.h"
@@ -134 +135 @@
     return psImageJpegColormapSet (map, "greyscale");
 }
+
+// XXX need to fix library references for this (psLib does not depend on libkapa)
+# if (0)
+// XXX Add colormap bar with scale (min -> max)
+// XXX Add option to plot the source overlay (pass in bDrawBuffer populated with points?)
+// XXX need to update bDraw APIs to pass in/out structure and avoid the local static 
+bool psImageJpegNew(const psImageJpegColormap *map, const psImage *image, const char *filename,
+                 float min, float max)
+{
+    PS_ASSERT_PTR_NON_NULL(map, false);
+    PS_ASSERT_IMAGE_NON_NULL(image, false);
+    PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false);
+    PS_ASSERT_VECTOR_NON_NULL(map->red, false);
+    PS_ASSERT_VECTOR_NON_NULL(map->green, false);
+    PS_ASSERT_VECTOR_NON_NULL(map->blue, false);
+    PS_ASSERT_PTR_NON_NULL(filename, false);
+    PS_ASSERT_INT_POSITIVE(strlen(filename), false);
+    PS_ASSERT_FLOAT_REAL(min, false);
+    PS_ASSERT_FLOAT_REAL(max, false);
+
+    float zero, scale;
+    struct jpeg_compress_struct cinfo;
+    struct jpeg_error_mgr jerr;
+
+    long pixel;
+    JSAMPLE *jpegLine;   // Points to data for current line
+    JSAMPROW jpegLineList[1];  // pointer to JSAMPLE row[s]
+    JSAMPLE *jpegImage;
+    JSAMPLE *outPix;
+
+    /* JPEG init calls */
+    cinfo.err = jpeg_std_error (&jerr);
+    jpeg_create_compress (&cinfo);
+
+    /* open file, prep for jpeg */
+    FILE *f = fopen(filename, "w");
+    if (!f) {
+        psError(PS_ERR_IO, true, "failed to open %s for output\n", filename);
+        return false;
+    }
+    jpeg_stdio_dest(&cinfo, f);
+
+    /* set up color jpeg buffers */
+    int quality = 75;
+    cinfo.image_width = image->numCols; // image width and height, in pixels
+    cinfo.image_height = image->numRows;
+    cinfo.input_components = 3;
+    cinfo.in_color_space = JCS_RGB;
+    jpeg_set_defaults (&cinfo);
+    jpeg_set_quality (&cinfo, quality, true); // limit to baseline-JPEG values
+    jpeg_start_compress (&cinfo, true);
+
+    psU8 *Rpix = map->red->data.U8;
+    psU8 *Gpix = map->green->data.U8;
+    psU8 *Bpix = map->blue->data.U8;
+
+    if (max == min) {
+        zero = min - 0.1;
+        scale = 256.0/0.2;
+    } else {
+        zero = min;
+        scale = 256.0/(max - min);
+    }
+
+    int dx = image->numCols;
+    int dy = image->numRows;
+
+    // output image buffer and line buffer
+    jpegLine = psAlloc (3*dx*sizeof(JSAMPLE));
+    jpegImage = psAlloc (3*dx*dy*sizeof(JSAMPLE));
+
+    // first copy the image data into the output buffer 
+    for (int j = 0; j < dy; j++) {
+        psF32 *row = image->data.F32[j];
+
+        outPix = jpegLine;
+        for (int i = 0; i < dx; i++, outPix += 3) {
+            if (isfinite(row[i])) {
+                pixel = PS_JPEG_SCALEVALUE(row[i],zero,scale);
+                outPix[0] = Rpix[pixel];
+                outPix[1] = Gpix[pixel];
+                outPix[2] = Bpix[pixel];
+            } else {
+                // XXX NAN value should be set per-color map
+                outPix[0] = 0x00;
+                outPix[1] = 0xff;
+                outPix[2] = 0x00;
+            }
+        }
+        memcpy (&jpegImage[j*3*dx], jpegLine, 3*dx);
+    }
+
+    bDrawBuffer *bdbuf = bDrawBufferCreate(dx, dy);
+    bDrawSetBuffer(bdbuf);
+    bDrawColor red = KapaColorByName("red");
+    bDrawSetStyle (red, 1, 0);
+    bDrawCircle(40.0, 20.0, 3.0);
+
+    {
+        int Npalette;
+        png_color *palette = KapaPNGPalette (&Npalette);
+        bDrawColor white = KapaColorByName ("white");
+        for (int j = 0; j < dy; j++) {
+            for (int i = 0; i < dx; i++) {
+                bDrawColor color = bdbuf[0].pixels[j][i];
+                if (color == white) continue;
+                jpegImage[j*3*dx + 3*i + 0] = palette[color].red;
+                jpegImage[j*3*dx + 3*i + 1] = palette[color].green;
+                jpegImage[j*3*dx + 3*i + 2] = palette[color].blue;
+            }
+        }
+    }
+    bDrawBufferFree (bdbuf);
+
+    // write out the image buffer
+    for (int j = 0; j < image->numRows; j++) {
+        jpegLineList[0] = &jpegImage[j*3*dx];
+        if (jpeg_write_scanlines(&cinfo, jpegLineList, 1) == 0) {
+            psError(PS_ERR_IO, true, "Unable to write line %d to JPEG", j);
+            psFree(jpegLine);
+            psFree(jpegImage);
+            fclose(f);
+            return false;
+        }
+    }
+
+    jpeg_finish_compress(&cinfo);
+    if (fclose(f) == EOF) {
+        psError(PS_ERR_IO, true, "Failed to close %s", filename);
+        psFree(jpegLine);
+        psFree(jpegImage);
+        return false;
+    }
+    jpeg_destroy_compress(&cinfo);
+
+    psFree(jpegLine);
+    psFree(jpegImage);
+    return true;
+}
+# endif
 
 bool psImageJpeg(const psImageJpegColormap *map, const psImage *image, const char *filename,

trunk/psLib/src/math/psMinimizeLMM.c

@@ -441 +441 @@
     }
 
+    // number of degrees of freedom for this fit
+    int nDOF = dy->n - params->n;
+
     // calculate initial alpha and beta, set chisq (min->value)
     min->value = psMinLM_SetABX(alpha, beta, params, paramMask, x, y, dy, func);
@@ -456 +459 @@
     }
 
-    // iterate until the tolerance is reached, or give up
-    while ((min->iter < min->maxIter) && ((min->lastDelta > min->tol) || !isfinite(min->lastDelta))) {
+    // iterate until: (a) nIter = min->iter or (b) (chisq / ndof) < maxChisq and deltaChisq < minTol (but don't stop unless Chisq is finite)
+    bool done = (min->iter >= min->maxIter);
+    while (!done) {
         psTrace("psLib.math", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
-        psTrace("psLib.math", 5, "Last delta is %f.  Min->tol is %f.\n", min->lastDelta, min->tol);
+        psTrace("psLib.math", 5, "Last delta is %f.  stop if < %f, accept if < %f\n", min->lastDelta, min->minTol, min->maxTol);
 
         // 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;
             continue;
@@ -482 +487 @@
         if (isnan(Chisq)) {
             min->iter ++;
+            if (min->iter >= min->maxIter) break;
             lambda *= 10.0;
             continue;
@@ -492 +498 @@
         psF32 rho = (min->value - Chisq) / dLinear;
 
-        psTrace("psLib.math", 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.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: %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);
 
         // dump some useful info if trace is defined
@@ -506 +510 @@
         /* rho is positive if the new chisq is smaller; allow for some insignificant change (slight negative rho) */
         if (rho >= -1e-6) {
-            min->lastDelta = (min->value - Chisq) / (dy->n - params->n);
+            min->lastDelta = (min->value - Chisq) / nDOF;
             min->value = Chisq;
             alpha  = psImageCopy(alpha, Alpha, PS_TYPE_F32);
@@ -516 +520 @@
         }
         min->iter++;
+
+        done = (min->iter >= min->maxIter);
+        
+        // check for convergence:
+        float chisqDOF = Chisq / nDOF;
+        if (!isfinite(min->maxChisqDOF) || ((chisqDOF < min->maxChisqDOF) && isfinite(min->lastDelta))) {
+            done |= (min->lastDelta < min->minTol);
+        }
     }
     psTrace("psLib.math", 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
@@ -549 +561 @@
         psFree(dy);
     }
-    if (min->iter == min->maxIter) {
-        psTrace("psLib.math", 3, "---- end (false) ----\n");
-        return(false);
-    }
-    psTrace("psLib.math", 3, "---- end (true) ----\n");
-    return(true);
+
+    // if the last improvement was at least as good as maxTol, accept the fit:
+    if (min->lastDelta <= min->maxTol) {
+        psTrace("psLib.math", 6, "---- end (true) ----\n");
+        return(true);
+    }
+    psTrace("psLib.math", 6, "---- end (false) ----\n");
+    return(false);
 }
 
@@ -588 +602 @@
 }
 
-psMinimization *psMinimizationAlloc(int maxIter,
-                                    float tol)
+psMinimization *psMinimizationAlloc(int maxIter, float minTol, float maxTol)
 {
     PS_ASSERT_INT_NONNEGATIVE(maxIter, NULL);
@@ -595 +608 @@
     psMinimization *min = psAlloc(sizeof(psMinimization));
     psMemSetDeallocator(min, (psFreeFunc)minimizationFree);
+
     P_PSMINIMIZATION_SET_MAXITER(min,maxIter);
-    P_PSMINIMIZATION_SET_TOL(min,tol);
+    P_PSMINIMIZATION_SET_MIN_TOL(min,minTol);
+    P_PSMINIMIZATION_SET_MAX_TOL(min,maxTol);
+
     min->value = 0.0;
     min->iter = 0;
     min->lastDelta = NAN;
+    min->maxChisqDOF = NAN;
 
     return(min);

trunk/psLib/src/math/psMinimizeLMM.h

@@ -33 +33 @@
 #define PS_MAX_LMM_ITERATIONS 100
 #define PS_MAX_MINIMIZE_ITERATIONS 100
-#define P_PSMINIMIZATION_SET_MAXITER(m,val) *(int*)&m->maxIter = val
-        #define P_PSMINIMIZATION_SET_TOL(m,val) *(float*)&m->tol = val
+#define P_PSMINIMIZATION_SET_MAXITER(m,val) { *(int*)&m->maxIter  = val; }
+#define P_PSMINIMIZATION_SET_MIN_TOL(m,val) { *(float*)&m->minTol = val; }
+#define P_PSMINIMIZATION_SET_MAX_TOL(m,val) { *(float*)&m->maxTol = val; }
 
                 typedef enum {
@@ -75 +76 @@
 typedef struct
 {
-    const int maxIter;                 ///< Convergence limit
-    const float tol;                   ///< Error Tolerance
+    const int maxIter;                  ///< Convergence limit
+    const float minTol;                 ///< Convergence Tolerance (stop if we reach this value)
+    const float maxTol;                 ///< Max Tolerance (accept fit if last improvement was this good)
     float value;                       ///< Value of function at minimum
     int iter;                          ///< Number of iterations to date
     float lastDelta;                   ///< The last difference for the fit
+    float maxChisqDOF;                 ///< for Chisq minimization, require that we reach here before checking tolerance
 }
 psMinimization;
@@ -102 +105 @@
 psMinimization *psMinimizationAlloc(
     int maxIter,                       ///< Number of minimization iterations to perform.
-    float tol                          ///< Requested error tolerance
+    float minTol,                      ///< stop if tolerance is less than this
+    float maxTol                       ///< accept fit if tolerance is less than this
 ) PS_ATTR_MALLOC;
 

trunk/psLib/src/math/psMinimizePowell.c

@@ -457 +457 @@
 
         mul = bracket->data.F32[1];
-        if ((fabs(a-b) < min->tol) && (fabs(b-c) < min->tol)) {
+        if ((fabs(a-b) < min->minTol) && (fabs(b-c) < min->minTol)) {
             PS_VECTOR_ADD_MULTIPLE(params, paramMask, line, params, mul);
             min->value = func(params, coords);
@@ -524 +524 @@
     psTrace("psLib.math", 4, "---- psMinimizePowell() begin ----\n");
     psTrace("psLib.math", 6, "min->maxIter is %d\n", min->maxIter);
-    psTrace("psLib.math", 6, "min->tol is %f\n", min->tol);
+    psTrace("psLib.math", 6, "min->minTol is %f\n", min->minTol);
 
     if (paramMask == NULL) {
@@ -573 +573 @@
         for (i=0;i<numDims;i++) {
             if (myParamMask->data.PS_TYPE_VECTOR_MASK_DATA[i] == 0) {
+
                 P_PSMINIMIZATION_SET_MAXITER((&dummyMin),PS_MINIMIZE_POWELL_LINEMIN_MAX_ITERATIONS);
-                *(float*)&dummyMin.tol = PS_MINIMIZE_POWELL_LINEMIN_ERROR_TOLERANCE;
+                P_PSMINIMIZATION_SET_MIN_TOL((&dummyMin),PS_MINIMIZE_POWELL_LINEMIN_ERROR_TOLERANCE);
+
                 mul = LineMin(&dummyMin, Q, ((psVector *) v->data[i]),
                               myParamMask, coords, func);
@@ -677 +679 @@
 
         // 8: Go to step 3 until the change is less than some tolerance.
-        if (fabs(baseFuncVal - currFuncVal) <= min->tol) {
+        if (fabs(baseFuncVal - currFuncVal) <= min->minTol) {
             psFree(v);
             psFree(pQP);

trunk/psLib/src/math/psStats.c

@@ -1211 +1211 @@
 
         // Fit a Gaussian to the data.
-        psMinimization *minimizer = psMinimizationAlloc(100, 0.01); // The minimizer information
+        psMinimization *minimizer = psMinimizationAlloc(100, 0.01, 1.0); // The minimizer information
         psVector *params = psVectorAlloc(2, PS_TYPE_F32); // Parameters for the Gaussian
         // Initial guess for the mean (index 0) and var (index 1).

trunk/psLib/test/math/Makefile.am

@@ -49 +49 @@
         tap_psMinimizePowell \
         tap_psSpline1D \
-        tap_psPolynomialMD
+        tap_psPolynomialMD \
+        tap_psPolynomialMD_sampleDark
 
 #       tap_psRandom