Changeset 27838 for branches/tap_branches/psLib

branches/tap_branches

Property svn:mergeinfo changed

/branches/pap_delete (added)	merged: 27537-27541
/trunk (added)	merged: 25902,25904,25906-25910,25913-25916,25918-25923,25925-25937,25939-25947,25949,25951-25953,25957-25958,25961-25963,25965-25966,25971-25974,25976-26005,26008,26012-26051,26053-26054,26067-26068,26070-26118,26138-26162,26169-26194,26201-26218,26242,26256-26266,26271-26294,26296-26297,26327-26329,26331,26333-26338,26349-26371,26376-26400,26403-26412,26437-26461,26465-26466,26477-26478,26481-26485,26489,26492-26500,26510-26521,26526-26529,26538,26544-26546,26548-26550,26554-26561,26563-26571,26584,26600-26602,26615-26622,26624,26629,26636-26641,26644,26651-26656,26658,26662-26665,26670,26672,26676-26680,26682-26683,26692,26694,26698-26699,26709-26710,26736,26744-26745,26749,26766,26796-26798,26815,26826,26829,26831,26848,26851,26856,26884,26886-27086,27088-27260,27262,27273-27280,27282-27290,27293-27294,27297-27304,27307-27311,27313-27335,27341-27349,27354-27356,27358-27367,27369-27371,27387-27389,27391-27396,27399-27400,27402-27404,27407-27408,27417-27420,27426-27427,27431-27436,27438-27439,27441-27486,27491,27503-27506,27510-27516,27518-27522,27524-27525,27527,27529-27534,27536,27542-27547,27549-27558,27560-27620,27631-27634,27638-27646,27652,27657-27692,27694-27705,27710-27712,27716-27725,27728-27731,27734-27767,27773-27783,27785-27786,27788-27790,27793-27813,27816-27822,27832-27834
/branches/czw_branch/cleanup	removed
/branches/eam_branches/20090522	removed
/branches/eam_branches/20090715	removed
/branches/eam_branches/20090820	removed
/branches/pap	removed
/branches/pap_mops	removed

branches/tap_branches/psLib

Property svn:mergeinfo set to (toggle deleted branches)

/branches/czw_branch/cleanup/psLib	merged	eligible
/branches/pap/psLib	merged	eligible
/trunk/psLib	merged	eligible
/branches/eam_branches/20090522/psLib	24238-24573
/branches/eam_branches/20090715/psLib	24799-25750
/branches/eam_branches/20090820/psLib	25139-25874
/branches/pap_mops/psLib	25137-25255

branches/tap_branches/psLib/share/Makefile.am

r10093	r27838
7	7	tab5.2b.dat \
8	8	tab5.2c.dat \
9		finals2000A.dat \
10		finals_all.dat
	9	finals2000A.dat
11	10

branches/tap_branches/psLib/src/fits/psFits.c

-              r25478
+              r27838
             int thisErrno = errno;      // Error number
             char errorBuf[MAX_STRING_LENGTH], *errorMsg;
 #if ((_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) && ! _GNU_SOURCE)
+#if (((_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) && ! _GNU_SOURCE) ||  __APPLE__)
             strerror_r(thisErrno, errorBuf, MAX_STRING_LENGTH);
             errorMsg = errorBuf;
 …
                 int thisErrno = errno;
                 char errorBuf[64], *errorMsg;
 # if ((_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) && ! _GNU_SOURCE)
+#if (((_POSIX_C_SOURCE >= 200112L || _XOPEN_SOURCE >= 600) && ! _GNU_SOURCE) ||  __APPLE__)
                 strerror_r (errno, errorBuf, 64);
                 errorMsg = errorBuf;

branches/tap_branches/psLib/src/fits/psFitsHeader.c

-              r25087
+              r27838
             // to preserve NAXISn etc for reference, so we don't do this.
+            if (keywordInList(name, noWriteFitsKeys) ||
+                (keyStarts && keywordStartsWith(name, noWriteFitsKeyStarts))) {
+                psTrace("psLib.fits", 3, "Not writing FITS keyword %s", name);
+                continue;
+            }
             // Options for compression
             if (compressing) {
 …
             } else if (keywordInList(name, noWriteCompressedKeys) ||
                        (keyStarts && keywordStartsWith(name, noWriteCompressedKeyStarts))) {
-                psTrace("psLib.fits", 3, "Not writing FITS keyword %s", name);
-                continue;
+            }
-            if (keywordInList(name, noWriteFitsKeys) ||
-                (keyStarts && keywordStartsWith(name, noWriteFitsKeyStarts))) {
                 psTrace("psLib.fits", 3, "Not writing FITS keyword %s", name);
                 continue;

branches/tap_branches/psLib/src/fits/psFitsImage.c

-              r25002
+              r27838
     psImage *inImage;                   // Image to read in
     if (floatType == PS_FITS_FLOAT_NONE) {
+        inImage = psImageRecycle(outImage, numCols, numRows, info->psDatatype);
+        outImage = psMemIncrRefCounter(inImage);
+        if (!outImage || outImage->type.type == info->psDatatype) {
+            outImage = psImageRecycle(outImage, numCols, numRows, info->psDatatype);
+            inImage = psMemIncrRefCounter(outImage);
+        } else {
+            outImage = psImageRecycle(outImage, numCols, numRows, outImage->type.type);
+            inImage = psImageAlloc(numCols, numRows, info->psDatatype);
+        }
     } else {
         inImage = psImageAlloc(numCols, numRows, info->psDatatype);
 …
         return NULL;
+    }
-    psFree(info);
     if (floatType != PS_FITS_FLOAT_NONE) {
         outImage = psFitsFloatImageFromDisk(outImage, inImage, floatType);
+    }
+    } else if (outImage->type.type != info->psDatatype) {
+        outImage = psImageCopy(outImage, inImage, outImage->type.type);
+    }
+    psFree(info);
     psFree(inImage);

branches/tap_branches/psLib/src/fits/psFitsScale.c

-              r25439
+              r27838
 #include "psImage.h"
 #include "psFits.h"
+#include "psStats.h"
+#include "psImageStats.h"
 #include "psImageBackground.h"
-#include "psStats.h"
 #include "psImageStructManip.h"
 …
 #define MEAN_STAT PS_STAT_ROBUST_MEDIAN // Statistic to use for mean
 #define STDEV_STAT PS_STAT_ROBUST_STDEV // Statistic to use for stdev
+#define DESPERATE_MEAN_STAT PS_STAT_SAMPLE_MEDIAN // Statistic to use for mean when deperate
+#define DESPERATE_STDEV_STAT PS_STAT_SAMPLE_QUARTILE // Statistic to use for stdev when desperate
 …
     psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
     psStats *stats = psStatsAlloc(MEAN_STAT | STDEV_STAT); // Statistics object
+    double mean, stdev;                                    // Mean and standard deviation
     if (!psImageBackground(stats, NULL, image, mask, maskVal, rng)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to perform statistics on image");
+        psFree(rng);
+        psFree(stats);
+        return false;
+        // It could be because the image is entirely masked, in which case we don't want to error
+        bool good = false;              // Any good pixels?
+// Find good pixels in an image, by image type
+#define GOOD_PIXELS_CASE(TYPE) \
+      case PS_TYPE_##TYPE: \
+        for (int y = 0; y < image->numRows && !good; y++) { \
+            for (int x = 0; x < image->numCols && !good; x++) { \
+                if (mask && (mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskVal)) { \
+                    continue; \
+                } \
+                if (!isfinite(image->data.TYPE[y][x])) { \
+                    continue; \
+                } \
+                good = true; \
+            } \
+        } \
+        break;
+        switch (image->type.type) {
+            GOOD_PIXELS_CASE(F32);
+            GOOD_PIXELS_CASE(F64);
+          default:
+            psAbort("Unsupported case: %x", image->type.type);
+        }
+        if (!good) {
+            psLogMsg("psLib.fits", PS_LOG_DETAIL, "Image has no good pixels, setting BSCALE = 1, BZERO = 0");
+            psErrorClear();
+            *bscale = 1.0;
+            *bzero = 0.0;
+            psFree(rng);
+            psFree(stats);
+            return true;
+        }
+        // There are some good pixels in there somewhere; psImageBackground just didn't find them
+        psLogMsg("psLib.fits", PS_LOG_DETAIL,
+                 "Couldn't measure background statistics for image quantisation; retrying.");
+        psErrorClear();
+        // Retry using all the available pixels
+        stats->nSubsample = image->numCols * image->numRows + 1;
+        if (!psImageStats(stats, image, mask, maskVal)) {
+            psLogMsg("psLib.fits", PS_LOG_DETAIL,
+                     "Couldn't measure background statistics for image quantisation (attempt 2); retrying.");
+            psErrorClear();
+            // Retry with desperate statistic
+            stats->options = DESPERATE_MEAN_STAT | DESPERATE_STDEV_STAT;
+            if (!psImageStats(stats, image, mask, maskVal)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to measure background statistics for image");
+                psFree(rng);
+                psFree(stats);
+                return false;
+            } else {
+                // Desperate retry
+                mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT);
+                stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT);
+            }
+        } else {
+            // Retry with all available pixels
+            mean = psStatsGetValue(stats, MEAN_STAT);
+            stdev = psStatsGetValue(stats, STDEV_STAT);
+        }
+    } else {
+        // First attempt
+        mean = psStatsGetValue(stats, MEAN_STAT);
+        stdev = psStatsGetValue(stats, STDEV_STAT);
+    }
     psFree(rng);
-    double mean = psStatsGetValue(stats, MEAN_STAT); // Mean
-    double stdev = psStatsGetValue(stats, STDEV_STAT); // Standard deviation
     psFree(stats);
     if (!isfinite(mean) || !isfinite(stdev)) {
 …
                 } else { \
                     value = (value - zero) * scale; \
                     if (options->fuzz) { \
+                    if (options->fuzz && (value - (int)value != 0.0)) { \
                        /* Add random factor [-0.5,0.5): adds a variance of 1/12, */ \
                        /* but preserves the expectation value */ \

branches/tap_branches/psLib/src/fits/psFitsTable.c

-              r24512
+              r27838
         switch (typecode) {
+           // TBYTE and TSHORT fall though to read into psS32
           case TBYTE:
           case TSHORT:
-          case TLONGLONG:
             READ_TABLE_ROW_CASE(TLONG, long, S32, S32);
+            READ_TABLE_ROW_CASE(TLONGLONG, long, S64, S64);
             READ_TABLE_ROW_CASE(TFLOAT, float, F32, F32);
             READ_TABLE_ROW_CASE(TDOUBLE, double, F64, F64);

branches/tap_branches/psLib/src/imageops/psImageBackground.c

-              r24481
+              r27838
 #include "psRandom.h"
 #include "psError.h"
+static int nFailures = 0;
+void psImageBackgroundInit() {
+    nFailures = 0;
+    return;
+}
 // XXX allow the user to choose the stats method?
 …
     long ny = image->numRows;
+    const int Npixels = nx*ny;                // Total number of pixels
+    psImage *bad = psImageAlloc(nx, ny, PS_TYPE_U8); // Image with bad pixels
+    psImageInit(bad, 0);
+    int Npixels = 0;                    // Total number of pixels
+    for (int y = 0; y < ny; y++) {
+        for (int x = 0; x < nx; x++) {
+            if (!isfinite(image->data.F32[y][x]) ||
+                (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] & maskValue)) {
+                bad->data.U8[y][x] = 0xFF;
+            }
+            Npixels++;
+        }
+    }
     const int Nsubset = (stats->nSubsample == 0) ? Npixels : PS_MIN(stats->nSubsample, Npixels); // Number of pixels in subset
 …
     long n = 0;                         // Number of actual pixels in subset
     if (Nsubset >= Npixels) {
         // if we have an image smaller than Nsubset, just loop over the image pixels
         for (int iy = 0; iy < ny; iy++) {
             for (int ix = 0; ix < nx; ix++) {
                 if (!isfinite(image->data.F32[iy][ix]) || (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue)) {
                     continue;
+                }
+        // if we have an image smaller than Nsubset, just loop over the (good) image pixels
+        for (int iy = 0; iy < ny; iy++) {
+            for (int ix = 0; ix < nx; ix++) {
+                if (bad->data.U8[iy][ix]) {
+                    continue;
+                }
                 float value = image->data.F32[iy][ix];
                 min = PS_MIN(value, min);
                 max = PS_MAX(value, max);
                 values->data.F32[n] = value;
                 n++;
+            }
+        }
+                float value = image->data.F32[iy][ix];
+                min = PS_MIN(value, min);
+                max = PS_MAX(value, max);
+                values->data.F32[n] = value;
+                n++;
+            }
+        }
     } else {
+        for (long i = 0; i < Nsubset; i++) {
+            double frnd = psRandomUniform(rng);
+            int pixel = Npixels * frnd;
+            int ix = pixel % nx;
+            int iy = pixel / nx;
+        // Subsample all pixels
+        // This is not optimal, since there may be a large masked fraction that leaves us with few good pixels.
+        // In that case, you should have set Nsubset.......
+        Npixels = nx * ny;
+        for (long i = 0; i < Nsubset; i++) {
+            double frnd = psRandomUniform(rng);
+            int pixel = Npixels * frnd;
+            int ix = pixel % nx;
+            int iy = pixel / nx;
             if (!isfinite(image->data.F32[iy][ix]) || (mask && mask->data.PS_TYPE_IMAGE_MASK_DATA[iy][ix] & maskValue)) {
                 continue;
+            }
+            if (bad->data.U8[iy][ix]) {
+                continue;
+            }
             float value = image->data.F32[iy][ix];
             min = PS_MIN(value, min);
             max = PS_MAX(value, max);
             values->data.F32[n] = value;
             n++;
+        }
+            float value = image->data.F32[iy][ix];
+            min = PS_MIN(value, min);
+            max = PS_MAX(value, max);
+            values->data.F32[n] = value;
+            n++;
+        }
+    }
+    psFree(bad);
     if (n < 0.01*Nsubset) {
+        psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL,
+                 "Unable to measure image background: too few data points (%ld)", n);
+        if ((nFailures < 3) || (nFailures % 100 == 0)) {
+            psLogMsg("psLib.psImageBackground", PS_LOG_DETAIL, "Unable to measure image background: too few data points (%ld) (%d failures)", n, nFailures);
+        }
+        nFailures ++;
         psFree(values);
         return false;
 …
         if (!psVectorSort(values, values)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to sort values.\n");
+            if ((nFailures < 3) || (nFailures % 100 == 0)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to sort values.(%d failures)\n", nFailures);
+            }
+            nFailures ++;
             psFree(values);
             return false;
 …
                 fclose (f);
+            }
+            psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for image background "
+                    "(%dx%d, (row0,col0) = (%d,%d)",
+                    image->numRows, image->numCols, image->row0, image->col0);
+            if ((nFailures < 3) || (nFailures % 100 == 0)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to measure statistics for image background "
+                        "(%dx%d, (row0,col0) = (%d,%d) (%d failures)\n",
+                        image->numRows, image->numCols, image->row0, image->col0, nFailures);
+            }
+            nFailures ++;
             psFree(values);
             return false;

branches/tap_branches/psLib/src/imageops/psImageBackground.h

r21183	r27838
22	22	#include <psRandom.h>
23	23
	24	void psImageBackgroundInit();
	25
24	26	// Get the background for an image
25	27	bool psImageBackground(psStats *stats, // desired measurement and options

branches/tap_branches/psLib/src/imageops/psImageConvolve.c

-              r25383
+              r27838
     return kernel;
+}
+bool psKernelTruncate(psKernel *kernel, float frac)
+{
+    PS_ASSERT_KERNEL_NON_NULL(kernel, false);
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(frac, 0.0, false);
+    PS_ASSERT_FLOAT_LESS_THAN(frac, 1.0, false);
+    if (frac == 0.0) {
+        // Nothing to do
+        return true;
+    }
+    int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Bounds
+    int maxSize = PS_MAX(PS_MAX(PS_MAX(-xMin, xMax), -yMin), yMax); // Maximum size
+    // Determine the threshold
+    // Summing absolute values because large negative deviations have power as well
+    double sumKernel = 0.0;             // Sum of the kernel
+    for (int y = yMin; y <= yMax; y++) {
+        for (int x = xMin; x <= xMax; x++) {
+            sumKernel += fabsf(kernel->kernel[y][x]);
+        }
+    }
+    float threshold = sumKernel * (1.0 - frac); // Threshold for truncation
+    // Find truncation size
+    int truncate = 0;                     // Truncation radius
+    for (int radius = 1; truncate == 0 && radius < maxSize; radius++) {
+        int uMin = PS_MAX(-radius, xMin);
+        int uMax = PS_MIN(radius, xMax);
+        int vMin = PS_MAX(-radius, yMin);
+        int vMax = PS_MIN(radius, yMax);
+        int r2 = PS_SQR(radius);
+        double sum = 0.0;
+        for (int v = vMin; v <= vMax; v++) {
+            int v2 = PS_SQR(v);
+            for (int u = uMin; u <= uMax; u++) {
+                int u2 = PS_SQR(u);
+                if (u2 + v2 <= r2) {
+                    sum += fabsf(kernel->kernel[v][u]);
+                }
+            }
+        }
+        if (sum > threshold) {
+            // This is the truncation radius
+            truncate = radius;
+        }
+    }
+    if (truncate == maxSize) {
+        // No truncation possible
+        return true;
+    }
+    // Truncate the kernel
+    {
+        int uMin = PS_MAX(-truncate, xMin);
+        int uMax = PS_MIN(truncate, xMax);
+        int vMin = PS_MAX(-truncate, yMin);
+        int vMax = PS_MIN(truncate, yMax);
+        int r2 = PS_SQR(truncate);
+        for (int v = vMin; v <= vMax; v++) {
+            int v2 = PS_SQR(v);
+            for (int u = uMin; u <= uMax; u++) {
+                int u2 = PS_SQR(u);
+                if (u2 + v2 > r2) {
+                    kernel->kernel[v][u] = 0.0;
+                }
+            }
+        }
+    }
+    kernel->xMin = PS_MAX(-truncate, kernel->xMin);
+    kernel->xMax = PS_MIN(truncate, kernel->xMax);
+    kernel->yMin = PS_MAX(-truncate, kernel->yMin);
+    kernel->yMax = PS_MIN(truncate, kernel->yMax);
+    return true;
+    }
 psImage *psImageConvolveDirect(psImage *out, const psImage *in, const psKernel *kernel)

branches/tap_branches/psLib/src/imageops/psImageConvolve.h

-              r25383
+              r27838
 );
+/// Truncate a kernel
+///
+/// Truncates the outer parts of the kernel where the contribution is below the nominated fraction of the
+/// total kernel.
+bool psKernelTruncate(
+    psKernel *in,                       ///< Kernel to be truncated
+    float frac                          ///< Fraction for truncation threshold
+    );
 /// Convolve an image with a kernel, using a direct convolution
 ///

branches/tap_branches/psLib/src/imageops/psImageCovariance.c

-              r24832
+              r27838
 #include "psTrace.h"
 #include "psBinaryOp.h"
+#include "psScalar.h"
+#include "psThread.h"
 #include "psImageCovariance.h"
+static bool threaded = false;           // Run threaded?
 psKernel *psImageCovarianceNone(void)
 …
+}
+psKernel *psImageCovarianceCalculate(const psKernel *kernel, const psKernel *covariance)
+{
+    PS_ASSERT_KERNEL_NON_NULL(kernel, NULL);
+    // See http://en.wikipedia.org/wiki/Error_propagation
+    //
+    // If
+    //     f_k = sum_i A_ik x_i
+    // is a set of functions, then the covariance matrix for f is given by:
+    //     M^f_ij = sum_k sum_l A_ik M^x_kl A_lj
+    // where M^x is the covariance matrix for x.
+    // Note that the errors in f are correlated (covariance) even if the errors in x are not.
+    psKernel *covar;                    // Covariance matrix to use
+    if (covariance) {
+        covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const
+    } else {
+        covar = psImageCovarianceNone();
+    }
+    // Check for non-finite elements
+    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
+        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+            if (!isfinite(kernel->kernel[y][x])) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                        "Non-finite covariance matrix element in kernel at %d,%d", x, y);
+                psFree(covar);
+                return NULL;
+            }
+        }
+    }
+    for (int y = covar->yMin; y <= covar->yMax; y++) {
+        for (int x = covar->xMin; x <= covar->xMax; x++) {
+            if (!isfinite(covar->kernel[y][x])) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                        "Non-finite covariance matrix element in covariance matrix at %d,%d", x, y);
+                psFree(covar);
+                return NULL;
+            }
+        }
+    }
+    // The above (double) sum for the covariance matrix means that, for each point in the output covariance
+    // matrix, we need to work out all combinations of getting to the central point via a kernel, input
+    // covariance matrix and another kernel.  This means that the resultant covariance matrix has twice the
+    // size of the kernel plus the size of the input covariance matrix.
+    int xMin = kernel->xMin - kernel->xMax + covar->xMin, xMax = kernel->xMax - kernel->xMin + covar->xMax;
+    int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax;
+    psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output
+/// Calculation of covariance matrix element when convolving
+static float imageCovarianceCalculate(const psKernel *covar, // Original covariance matrix
+                                      const psKernel *kernel, // Convolution kernel
+                                      int x, int y            // Coordinates in output covariance matrix
+                                      )
+{
+    psAssert(covar, "Require covariance matrix");
+    psAssert(kernel, "Require kernel");
     // Need to go:
 …
     // from the source coordinate), we take the smallest possible (because everything else is zero outside).
+    double total = 0.0;                 // Total covariance
+    // Range for v
+    int vMin = PS_MAX(kernel->yMin + covar->yMin, y + kernel->yMin);
+    int vMax = PS_MIN(kernel->yMax + covar->yMax, y + kernel->yMax);
+    // Range for u
+    int uMin = PS_MAX(kernel->xMin + covar->xMin, x + kernel->xMin);
+    int uMax = PS_MIN(kernel->xMax + covar->xMax, x + kernel->xMax);
+    double sum = 0.0;           // Sum for value of covariance matrix at (x,y)
+    for (int v = vMin; v <= vMax; v++) {
+        // Range for q
+        int qMin = PS_MAX(v + covar->yMin, kernel->yMin);
+        int qMax = PS_MIN(v + covar->yMax, kernel->yMax);
+        for (int u = uMin; u <= uMax; u++) {
+            // Range for p
+            int pMin = PS_MAX(u + covar->xMin, kernel->xMin);
+            int pMax = PS_MIN(u + covar->xMax, kernel->xMax);
+            double xyuvValue = kernel->kernel[v-y][u-x]; // Value for (x,y) --> (u,v)
+            double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0)
+            for (int q = qMin; q <= qMax; q++) {
+                for (int p = pMin; p <= pMax; p++) {
+                    uvpqValue += (double)covar->kernel[q-v][p-u] * (double)kernel->kernel[q][p];
+                }
+            }
+            sum += xyuvValue * uvpqValue;
+        }
+    }
+    return sum;
+}
+/// Thread entry point for calculation of covariance matrix element when convolving
+static bool imageCovarianceCalculateThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert(job->args, "No job arguments");
+    psAssert(job->args->n == 5, "Wrong number of job arguments: %ld", job->args->n);
+    psKernel *out = job->args->data[0]; // Output covariance matrix
+    const psKernel *covar = job->args->data[1]; // Input covariance matrix
+    const psKernel *kernel = job->args->data[2]; // Convolution kernel
+    int x = PS_SCALAR_VALUE(job->args->data[3], S32); // x coordinate in output covariance matrix
+    int y = PS_SCALAR_VALUE(job->args->data[4], S32); // y coordinate in output covariance matrix
+    out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
+    return true;
+}
+psKernel *psImageCovarianceCalculate(const psKernel *kernel, const psKernel *covariance)
+{
+    PS_ASSERT_KERNEL_NON_NULL(kernel, NULL);
+    // See http://en.wikipedia.org/wiki/Error_propagation
+    //
+    // If
+    //     f_k = sum_i A_ik x_i
+    // is a set of functions, then the covariance matrix for f is given by:
+    //     M^f_ij = sum_k sum_l A_ik M^x_kl A_lj
+    // where M^x is the covariance matrix for x.
+    // Note that the errors in f are correlated (covariance) even if the errors in x are not.
+    psKernel *covar;                    // Covariance matrix to use
+    if (covariance) {
+        covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const
+    } else {
+        covar = psImageCovarianceNone();
+    }
+    // Check for non-finite elements
+    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
+        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+            if (!isfinite(kernel->kernel[y][x])) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                        "Non-finite covariance matrix element in kernel at %d,%d", x, y);
+                psFree(covar);
+                return NULL;
+            }
+        }
+    }
+    for (int y = covar->yMin; y <= covar->yMax; y++) {
+        for (int x = covar->xMin; x <= covar->xMax; x++) {
+            if (!isfinite(covar->kernel[y][x])) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                        "Non-finite covariance matrix element in covariance matrix at %d,%d", x, y);
+                psFree(covar);
+                return NULL;
+            }
+        }
+    }
+    // The above (double) sum for the covariance matrix means that, for each point in the output covariance
+    // matrix, we need to work out all combinations of getting to the central point via a kernel, input
+    // covariance matrix and another kernel.  This means that the resultant covariance matrix has twice the
+    // size of the kernel plus the size of the input covariance matrix.
+    int xMin = kernel->xMin - kernel->xMax + covar->xMin, xMax = kernel->xMax - kernel->xMin + covar->xMax;
+    int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax;
+    psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output
     for (int y = yMin; y <= yMax; y++) {
-        // Range for v
-        int vMin = PS_MAX(kernel->yMin + covar->yMin, y + kernel->yMin);
-        int vMax = PS_MIN(kernel->yMax + covar->yMax, y + kernel->yMax);
         for (int x = xMin; x <= xMax; x++) {
+            // Range for u
+            int uMin = PS_MAX(kernel->xMin + covar->xMin, x + kernel->xMin);
+            int uMax = PS_MIN(kernel->xMax + covar->xMax, x + kernel->xMax);
+            double sum = 0.0;           // Sum for value of covariance matrix at (x,y)
+            for (int v = vMin; v <= vMax; v++) {
+                // Range for q
+                int qMin = PS_MAX(v + covar->yMin, kernel->yMin);
+                int qMax = PS_MIN(v + covar->yMax, kernel->yMax);
+                for (int u = uMin; u <= uMax; u++) {
+                    // Range for p
+                    int pMin = PS_MAX(u + covar->xMin, kernel->xMin);
+                    int pMax = PS_MIN(u + covar->xMax, kernel->xMax);
+                    double xyuvValue = kernel->kernel[v-y][u-x]; // Value for (x,y) --> (u,v)
+                    double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0)
+                    for (int q = qMin; q <= qMax; q++) {
+                        for (int p = pMin; p <= pMax; p++) {
+                            uvpqValue += (double)covar->kernel[q-v][p-u] * (double)kernel->kernel[q][p];
+                        }
+                    }
+                    sum += xyuvValue * uvpqValue;
+            if (threaded) {
+                psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE");
+                psArrayAdd(job->args, 1, out);
+                psArrayAdd(job->args, 1, covar);
+                psArrayAdd(job->args, 1, (psKernel*)kernel); // Casting away const
+                PS_ARRAY_ADD_SCALAR(job->args, x, PS_TYPE_S32);
+                PS_ARRAY_ADD_SCALAR(job->args, y, PS_TYPE_S32);
+                if (!psThreadJobAddPending(job)) {
+                    psFree(covar);
+                    return NULL;
+                }
+            }
+            out->kernel[y][x] = sum;
+            total += sum;
+        }
+    }
+    psTrace("psLib.imageops", 3, "Total covariance: %lf ; Central variance: %f\n", total, out->kernel[0][0]);
+                psFree(job);
+            } else {
+                out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
+            }
+        }
+    }
     psFree(covar);
+    if (threaded && !psThreadPoolWait(true)) {
+        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+        return false;
+    }
     return out;
+}
+float psImageCovarianceCalculateFactor(const psKernel *kernel, const psKernel *covariance)
+{
+    psKernel *covar;                    // Covariance matrix to use
+    if (covariance) {
+        covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const
+    } else {
+        covar = psImageCovarianceNone();
+    }
+    // Check for non-finite elements
+    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
+        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
+            if (!isfinite(kernel->kernel[y][x])) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                        "Non-finite covariance matrix element in kernel at %d,%d", x, y);
+                psFree(covar);
+                return NAN;
+            }
+        }
+    }
+    for (int y = covar->yMin; y <= covar->yMax; y++) {
+        for (int x = covar->xMin; x <= covar->xMax; x++) {
+            if (!isfinite(covar->kernel[y][x])) {
+                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                        "Non-finite covariance matrix element in covariance matrix at %d,%d", x, y);
+                psFree(covar);
+                return NAN;
+            }
+        }
+    }
+    float factor = imageCovarianceCalculate(covar, kernel, 0, 0); // Covariance factor
+    psFree(covar);
+    return factor;
+}
+// Calculation of covariance matrix element when binning
+static float imageCovarianceBin(const psKernel *covar, // Original covariance matrix
+                                int bin,               // Binning factor
+                                float binVal,          // Convolution kernel value for binning
+                                int x, int y           // Coordinates in output covariance matrix
+    )
+{
+    psAssert(covar, "Require covariance matrix");
+    psAssert(bin > 0 && binVal > 0, "Require binning: %d %f", bin, binVal);
+    int binMin = -(bin - 1) / 2, binMax = bin / 2; // Range of "kernel"
+    // Range for v
+    int vMin = PS_MAX(binMin + covar->yMin, bin * y + binMin);
+    int vMax = PS_MIN(binMax + covar->yMax, bin * y + binMax);
+    // Range for u
+    int uMin = PS_MAX(binMin + covar->xMin, bin * x + binMin);
+    int uMax = PS_MIN(binMax + covar->xMax, bin * x + binMax);
+    double sum = 0.0;           // Sum for value of covariance matrix at (x,y)
+    for (int v = vMin; v <= vMax; v++) {
+        // Range for q
+        int qMin = PS_MAX(v + covar->yMin, binMin);
+        int qMax = PS_MIN(v + covar->yMax, binMax);
+        for (int u = uMin; u <= uMax; u++) {
+            // Range for p
+            int pMin = PS_MAX(u + covar->xMin, binMin);
+            int pMax = PS_MIN(u + covar->xMax, binMax);
+            double xyuvValue = binVal; // Value for (x,y) --> (u,v)
+            double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0)
+            for (int q = qMin; q <= qMax; q++) {
+                for (int p = pMin; p <= pMax; p++) {
+                    uvpqValue += (double)covar->kernel[q-v][p-u] * (double)binVal;
+                }
+            }
+            sum += xyuvValue * uvpqValue;
+        }
+    }
+    return sum;
+}
+/// Thread entry point for calculation of covariance matrix element when binning
+static bool imageCovarianceBinThread(psThreadJob *job)
+{
+    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
+    psAssert(job->args, "No job arguments");
+    psAssert(job->args->n == 6, "Wrong number of job arguments: %ld", job->args->n);
+    psKernel *out = job->args->data[0]; // Output covariance matrix
+    const psKernel *covar = job->args->data[1]; // Input covariance matrix
+    int bin = PS_SCALAR_VALUE(job->args->data[2], S32); // Binning factor
+    float binVal = PS_SCALAR_VALUE(job->args->data[3], F32); // Convolution kernel value for binning
+    int x = PS_SCALAR_VALUE(job->args->data[4], S32); // x coordinate in output covariance matrix
+    int y = PS_SCALAR_VALUE(job->args->data[5], S32); // y coordinate in output covariance matrix
+    out->kernel[y][x] = imageCovarianceBin(covar, bin, binVal, x, y);
+    return true;
+}
 psKernel *psImageCovarianceBin(int bin, const psKernel *covariance, bool average)
 …
     psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output
-    double total = 0.0;                 // Total covariance
     for (int y = yMin; y <= yMax; y++) {
-        // Range for v
-        int vMin = PS_MAX(binMin + covar->yMin, bin * y + binMin);
-        int vMax = PS_MIN(binMax + covar->yMax, bin * y + binMax);
         for (int x = xMin; x <= xMax; x++) {
+            // Range for u
+            int uMin = PS_MAX(binMin + covar->xMin, bin * x + binMin);
+            int uMax = PS_MIN(binMax + covar->xMax, bin * x + binMax);
+            double sum = 0.0;           // Sum for value of covariance matrix at (x,y)
+            for (int v = vMin; v <= vMax; v++) {
+                // Range for q
+                int qMin = PS_MAX(v + covar->yMin, binMin);
+                int qMax = PS_MIN(v + covar->yMax, binMax);
+                for (int u = uMin; u <= uMax; u++) {
+                    // Range for p
+                    int pMin = PS_MAX(u + covar->xMin, binMin);
+                    int pMax = PS_MIN(u + covar->xMax, binMax);
+                    double xyuvValue = binVal; // Value for (x,y) --> (u,v)
+                    double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0)
+                    for (int q = qMin; q <= qMax; q++) {
+                        for (int p = pMin; p <= pMax; p++) {
+                            uvpqValue += (double)covar->kernel[q-v][p-u] * (double)binVal;
+                        }
+                    }
+                    sum += xyuvValue * uvpqValue;
+            if (threaded) {
+                psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_COVARIANCE_BIN");
+                psArrayAdd(job->args, 1, out);
+                psArrayAdd(job->args, 1, covar);
+                PS_ARRAY_ADD_SCALAR(job->args, bin, PS_TYPE_S32);
+                PS_ARRAY_ADD_SCALAR(job->args, binVal, PS_TYPE_F32);
+                PS_ARRAY_ADD_SCALAR(job->args, x, PS_TYPE_S32);
+                PS_ARRAY_ADD_SCALAR(job->args, y, PS_TYPE_S32);
+                if (!psThreadJobAddPending(job)) {
+                    psFree(covar);
+                    return NULL;
+                }
+            }
+            out->kernel[y][x] = sum;
+            total += sum;
+        }
+    }
+    psTrace("psLib.imageops", 3, "Total covariance: %lf ; Central variance: %f\n", total, out->kernel[0][0]);
+                psFree(job);
+            } else {
+                out->kernel[y][x] = imageCovarianceBin(covar, bin, binVal, x, y);
+            }
+        }
+    }
     psFree(covar);
+    if (threaded && !psThreadPoolWait(true)) {
+        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
+        return false;
+    }
     return out;
+}
 …
+}
+psKernel *psImageCovarianceSum(const psArray *array)
+float psImageCovarianceFactorForAperture(const psKernel *covar, float radius)
+{
+    if (!covar) return 1.0;
+    float Sum = 0.0;
+    for (int y = covar->yMin; y <= covar->yMax; y++) {
+        if (y < -radius) continue;
+        if (y > +radius) continue;
+        for (int x = covar->xMin; x <= covar->xMax; x++) {
+            if (x < -radius) continue;
+            if (x > +radius) continue;
+            if (hypot(x, y) > radius) continue;
+            psAssert (isfinite(covar->kernel[y][x]), "invalid NAN in covariance matrix");
+            Sum += covar->kernel[y][x];
+        }
+    }
+    return Sum;
+}
+float psImageCovarianceSum(const psKernel *covar)
+{
+    PS_ASSERT_KERNEL_NON_NULL(covar, NAN);
+    int xMin = covar->xMin, xMax = covar->xMax, yMin = covar->yMin, yMax = covar->yMax; // Range for covariance
+    double sum = 0.0;                                                                   // Sum of covariance
+    for (int y = yMin; y <= yMax; y++) {
+        for (int x = xMin; x <= xMax; x++) {
+            sum += covar->kernel[y][x];
+        }
+    }
+    return sum;
+}
+psKernel *psImageCovarianceAverage(const psArray *array)
+{
     PS_ASSERT_ARRAY_NON_NULL(array, NULL);
     PS_ASSERT_ARRAY_NON_EMPTY(array, NULL);
+    int xMin = INT_MAX, xMax = INT_MIN, yMin = INT_MAX, yMax = INT_MIN; // Range for covariance
+    int num = 0;                        // Number of good matrices to sum
+    for (int i = 0; i < array->n; i++) {
+        psKernel *covar = array->data[i]; // Covariance matrix
+        if (!covar) {
+            continue;
+        }
+        xMin = PS_MIN(xMin, covar->xMin);
+        xMax = PS_MAX(xMax, covar->xMax);
+        yMin = PS_MIN(yMin, covar->yMin);
+        yMax = PS_MAX(yMax, covar->yMax);
+        num++;
+    }
+    if (num == 0) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No covariance matrices supplied for summation");
+        return NULL;
+    }
+    psKernel *sum = psKernelAlloc(xMin, xMax, yMin, yMax); // Summed covariance
+    for (int i = 0; i < array->n; i++) {
+        psKernel *covar = array->data[i]; // Covariance matrix
+        if (!covar) {
+            continue;
+        }
+        for (int y = covar->yMin; y <= covar->yMax; y++) {
+            for (int x = covar->xMin; x <= covar->xMax; x++) {
+                if (!isfinite(covar->kernel[y][x])) {
+                    psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                            "Non-finite covariance matrix element at %d,%d for input %d",
+                            x, y, i);
+                    psFree(sum);
+                    return NULL;
+                }
+                sum->kernel[y][x] += covar->kernel[y][x];
+            }
+        }
+    }
+    return sum;
+}
+psKernel *psImageCovarianceAverage(const psArray *array)
+{
+    PS_ASSERT_ARRAY_NON_NULL(array, NULL);
+    PS_ASSERT_ARRAY_NON_EMPTY(array, NULL);
+    int num = 0;                        // Number of good matrices to average
+    for (int i = 0; i < array->n; i++) {
+        psKernel *covar = array->data[i]; // Covariance matrix
+        if (covar) {
+            num++;
+        }
+    }
+    if (num == 0) {
+        psError(PS_ERR_BAD_PARAMETER_SIZE, true, "No covariance matrices supplied for averaging.");
+        return NULL;
+    }
+    psKernel *sum = psImageCovarianceSum(array); // Sum of covariances
+    psBinaryOp(sum->image, sum->image, "/", psScalarAlloc(num, PS_TYPE_F32));
+    return sum;
+    psVector *weights = psVectorAlloc(array->n, PS_TYPE_F32); // Weights to apply
+    psVectorInit(weights, 1.0);
+    psKernel *out = psImageCovarianceAverageWeighted(array, weights);
+    psFree(weights);
+    return out;
+}
 …
     psKernel *sum = psKernelAlloc(xMin, xMax, yMin, yMax); // Summed covariance
+    psImageInit(sum->image, 0.0);
     for (int i = 0; i < array->n; i++) {
         psKernel *covar = array->data[i]; // Covariance matrix
 …
     return true;
+}
+bool psImageCovarianceSetThreads(bool set)
+{
+    bool old = threaded;                // Old value
+    if (set && !threaded) {
+        {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 5);
+            task->function = &imageCovarianceCalculateThread;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
+        {
+            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_BIN", 6);
+            task->function = &imageCovarianceBinThread;
+            psThreadTaskAdd(task);
+            psFree(task);
+        }
+    } else if (!set && threaded) {
+        psThreadTaskRemove("PSLIB_IMAGE_COVARIANCE_CALCULATE");
+        psThreadTaskRemove("PSLIB_IMAGE_COVARIANCE_BIN");
+    }
+    threaded = set;
+    return old;
+}
+bool psImageCovarianceGetThreads(void)
+{
+    return threaded;
+}

branches/tap_branches/psLib/src/imageops/psImageCovariance.h

-              r24832
+              r27838
     );
+/// Sum multiple covariance pseudo-matrices
+psKernel *psImageCovarianceSum(
+    const psArray *array                ///< Array of covariance pseudo-matrices
+/// Return the pixel-to-pixel covariance factor following calculation
+///
+/// This doesn't require calculation of the entire covariance matrix, so is much faster.
+float psImageCovarianceCalculateFactor(
+    const psKernel *kernel,             ///< Convolution kernel
+    const psKernel *covariance          ///< Current covariance pseudo-matrix
+    );
+/// Return the covariance factor for an aperture of a given radius
+float psImageCovarianceFactorForAperture(const psKernel *covar, float radius);
+/// Return the sum of the covariance pseudo-matrix
+float psImageCovarianceSum(
+    const psKernel *covariance          ///< Covariance pseudo-matrix
     );
 …
     );
+/// Control threading for image covariance functions
+///
+/// Returns old threading status
+bool psImageCovarianceSetThreads(bool threaded ///< Run image covariance functions threaded?
+    );
+/// Return whether image covariance functions are threaded
+bool psImageCovarianceGetThreads(void);
 /// @}

branches/tap_branches/psLib/src/imageops/psImageInterpolate.c

-              r23231
+              r27838
+{
     // Casting away const
     psFree((psImage*)interp->image);
     psFree((psImage*)interp->mask);
     psFree((psImage*)interp->variance);
     psFree((psImage*)interp->kernel);
     psFree((psImage*)interp->kernel2);
     psFree((psVector*)interp->sumKernel2);
+    psFree(interp->image);
+    psFree(interp->mask);
+    psFree(interp->variance);
+    psFree(interp->kernel);
+    psFree(interp->kernel2);
+    psFree(interp->sumKernel2);
+}

branches/tap_branches/psLib/src/jpeg/psImageJpeg.c

-              r24220
+              r27838
             if (isfinite(row[i])) {
                 pixel = PS_JPEG_SCALEVALUE(row[i],zero,scale);
                 outPix[0] = Rpix[pixel];
                 outPix[1] = Gpix[pixel];
                 outPix[2] = Bpix[pixel];
+                outPix[0] = Rpix[pixel];
+                outPix[1] = Gpix[pixel];
+                outPix[2] = Bpix[pixel];
             } else {
+                // XXX NAN value should be set per-color map
+                outPix[0] = 0xff;
+                outPix[1] = 0x00;
+                outPix[2] = 0xff;
+            }
+        }
+        jpeg_write_scanlines(&cinfo, jpegLineList, 1);
+                // XXX NAN value should be set per-color map
+                outPix[0] = 0xff;
+                outPix[1] = 0x00;
+                outPix[2] = 0xff;
+            }
+        }
+        if (jpeg_write_scanlines(&cinfo, jpegLineList, 1) == 0) {
+            psError(PS_ERR_IO, true, "Unable to write line %d to JPEG", j);
+            psFree(jpegLine);
+            return false;
+        }
+    }
     jpeg_finish_compress(&cinfo);
+    fclose(f);
+    if (fclose(f) == EOF) {
+        psError(PS_ERR_IO, true, "Failed to close %s", filename);
+        psFree(jpegLine);
+        return false;
+    }
     jpeg_destroy_compress(&cinfo);

branches/tap_branches/psLib/src/math/psBinaryOp.c

r17050	r27838
379	379	}
380	380
381		psMathType* psBinaryOp(psPtr out, ~~const psPtr in1, const char *op, const~~ psPtr in2)
	381	psMathType* psBinaryOp(psPtr out, psPtr in1, const char *op, psPtr in2)
382	382	{
383	383

branches/tap_branches/psLib/src/math/psBinaryOp.h

-              r11248
+              r27838
 psMathType* psBinaryOp(
     psPtr out,                         ///< Output type, either psImage or psVector.
     const psPtr in1,                   ///< First input, either psImage or psVector.
+    psPtr in1,                   ///< First input, either psImage or psVector.
     const char *op,                    ///< Operator.
     const psPtr in2                    ///< Second input, either psImage or psVector.
+    psPtr in2                    ///< Second input, either psImage or psVector.
 );

branches/tap_branches/psLib/src/math/psHistogram.c

-              r21183
+              r27838
     psTrace("psLib.math", 3, "---- %s() begin  ----\n", __func__);
     psTrace("psLib.math", 5, "(lower, upper, n) is (%f, %f, %d)\n", lower, upper, n);
     PS_ASSERT_INT_POSITIVE(n, NULL);
     PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(upper, lower, NULL);
+    psAssert(n > 0, "Number of bins must be positive");
+    psAssert(upper >= lower, "Bounds must be sensical");
     // Allocate memory for the new histogram structure.  If there are N bins, then there are N+1 bounds to
 …
 static void histogramFree(psHistogram* myHist)
+{
     psFree((void *)myHist->bounds);
+    psFree(myHist->bounds);
     psFree(myHist->nums);
+}
 …
                     double binSize = (out->bounds->data.F32[out->nums->n] - out->bounds->data.F32[0]) / (float) out->nums->n; // Histogram bin size
                     binNum = (inF32->data.F32[i] - out->bounds->data.F32[0]) / binSize;
                     binNum = PS_MAX (binNum, 0);
                     binNum = PS_MIN (binNum, numBins - 1);
                     // value is in bin 'i' if bound[i] <= value < bound[i]
                     // we may slightly overshoot or undershoot.  creep up or down on the true bin
                     if (inF32->data.F32[i] < out->bounds->data.F32[binNum]) {
                         psTrace("psLib.math", 6, "missed target bin, adjusting: %f vs %f to %f\n", inF32->data.F32[i], out->bounds->data.F32[binNum], out->bounds->data.F32[binNum+1]);
                         while ((inF32->data.F32[i] < out->bounds->data.F32[binNum]) && (binNum > 0)) {
                             binNum --;
+                        }
+                    }
                     if (inF32->data.F32[i] >= out->bounds->data.F32[binNum+1]) {
                         psTrace("psLib.math", 6, "missed target bin, adjusting: %f vs %f to %f\n", inF32->data.F32[i], out->bounds->data.F32[binNum], out->bounds->data.F32[binNum+1]);
                         while ((inF32->data.F32[i] >= out->bounds->data.F32[binNum+1]) && (binNum < numBins - 1)) {
                             binNum ++;
+                        }
+                    }
+                    binNum = PS_MAX (binNum, 0);
+                    binNum = PS_MIN (binNum, numBins - 1);
+                    // value is in bin 'i' if bound[i] <= value < bound[i]
+                    // we may slightly overshoot or undershoot.  creep up or down on the true bin
+                    if (inF32->data.F32[i] < out->bounds->data.F32[binNum]) {
+                        psTrace("psLib.math", 6, "missed target bin, adjusting: %f vs %f to %f\n", inF32->data.F32[i], out->bounds->data.F32[binNum], out->bounds->data.F32[binNum+1]);
+                        while ((inF32->data.F32[i] < out->bounds->data.F32[binNum]) && (binNum > 0)) {
+                            binNum --;
+                        }
+                    }
+                    if (inF32->data.F32[i] >= out->bounds->data.F32[binNum+1]) {
+                        psTrace("psLib.math", 6, "missed target bin, adjusting: %f vs %f to %f\n", inF32->data.F32[i], out->bounds->data.F32[binNum], out->bounds->data.F32[binNum+1]);
+                        while ((inF32->data.F32[i] >= out->bounds->data.F32[binNum+1]) && (binNum < numBins - 1)) {
+                            binNum ++;
+                        }
+                    }
                     if (errorsF32) {
 …
                     // correct bin number requires a bit more work.
                     tmpScalar.data.F32 = inF32->data.F32[i];
                     psVectorBinaryDisectResult result;
+                    psVectorBinaryDisectResult result;
                     binNum = psVectorBinaryDisect(&result, out->bounds, &tmpScalar);
                     if (result != PS_BINARY_DISECT_PASS) {
                         continue;
+                    }
                     if (errorsF32 != NULL) {
                         if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errors->data.F32[i])) {
                             psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram "
                                      "bins with the errors vector.\n");
+                        }
                     } else {
                         out->nums->data.F32[binNum] += 1.0;
+                    }
+                    if (result != PS_BINARY_DISECT_PASS) {
+                        continue;
+                    }
+                    if (errorsF32 != NULL) {
+                        if (!UpdateHistogramBins(binNum, out, inF32->data.F32[i], errors->data.F32[i])) {
+                            psLogMsg(__func__, PS_LOG_WARN, "WARNING: Failed to update the histogram "
+                                     "bins with the errors vector.\n");
+                        }
+                    } else {
+                        out->nums->data.F32[binNum] += 1.0;
+                    }
+                }
+            }

branches/tap_branches/psLib/src/math/psMatrix.c

-              r24122
+              r27838
 LHS_NAME.data  = RHS_NAME;
+/*****************************************************************************/
+/* FILE STATIC FUNCTIONS                                                     */
+/*****************************************************************************/
+static void  psVectorToGslVector(gsl_vector *outGslVector, const psVector *inVector);
+static void gslVectorToPsVector(psVector *outVector, gsl_vector *inGslVector);
+static void  psImageToGslMatrix(gsl_matrix *outGslMatrix, const psImage *inImage);
+static void gslMatrixToPsImage(psImage *outImage, gsl_matrix *inGslMatrix);
+////////////////////////////////////////////////////////////////////////////////
+// Conversion functions
+////////////////////////////////////////////////////////////////////////////////
+// gsl_vector holds *doubles*, so we can directly copy F64, but need to convert F32
 /** Static function to copy psF32 or psF64 vector data to a GSL vector */
+static void  psVectorToGslVector(gsl_vector *outGslVector,
+                                 const psVector *inVector)
+{
+    psU32 i = 0;
+    psU32 n = 0;
+    n = inVector->n;
+    for(i=0; i<n; i++) {
+        if(inVector->type.type == PS_TYPE_F32) {
+            outGslVector->data[i] = (psF64)inVector->data.F32[i];
+static void vectorPStoGSL(gsl_vector *out, const psVector *in)
+{
+    psAssert(out->size == in->n, "Sizes don't match!");
+    long n = in->n;                     // Size of input
+    switch (in->type.type) {
+      case PS_TYPE_F32:
+        for (long i = 0; i < n; i++) {
+            out->data[i] = in->data.F32[i];
+        }
+        break;
+      case PS_TYPE_F64:
+        memcpy(out->data, in->data.F64, n * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+        break;
+      default:
+        psAbort("Unsupported vector type: %x\n", in->type.type);
+    }
+    return;
+}
+/** Static function to copy GSL vector data to a psF32 or psF64 vector */
+static void vectorGSLtoPS(psVector *out, const gsl_vector *in)
+{
+    psAssert(in->size == out->n, "Sizes don't match!");
+    long n = out->n;                    // Size of output
+    switch (out->type.type) {
+      case PS_TYPE_F32:
+        for (long i = 0; i < n; i++) {
+            out->data.F32[i] = in->data[i];
+        }
+        break;
+      case PS_TYPE_F64:
+        memcpy(out->data.F64, in->data, n * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+        break;
+      default:
+        psAbort("Unsupported vector type: %x\n", out->type.type);
+    }
+    return;
+}
+/** Static function to copy psF32 or psF64 image data to a GSL matrix */
+static void matrixPStoGSL(gsl_matrix *out, const psImage *in)
+{
+    psAssert(out->size1 == in->numRows && out->size2 == in->numCols, "Sizes don't match!");
+    int numCols = in->numCols, numRows = in->numRows; // Size of matrix
+    switch (in->type.type) {
+      case PS_TYPE_F32:
+        for (int y = 0, i = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++, i++) {
+                out->data[i] = in->data.F32[y][x];
+            }
+        }
+        break;
+      case PS_TYPE_F64:
+        if (in->parent|| out->tda != out->size1) {
+            for (int y = 0, i = 0; y < numRows; y++, i += numCols) {
+                memcpy(&out->data[i], in->data.F64[y], numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+            }
         } else {
+            outGslVector->data[i] = inVector->data.F64[i];
+        }
+    }
+}
+/** Static function to copy GSL vector data to a psF32 or psF64 vector */
+static void gslVectorToPsVector(psVector *outVector,
+                                gsl_vector *inGslVector)
+{
+    psU32 i = 0;
+    psU32 n = 0;
+    n = outVector->n;
+    for(i=0; i<n; i++) {
+        if(outVector->type.type == PS_TYPE_F32) {
+            outVector->data.F32[i] = (psF32)inGslVector->data[i];
+            memcpy(out->data, in->p_rawDataBuffer, numCols * numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+        }
+        break;
+      default:
+        psAbort("Unsupported vector type: %x\n", in->type.type);
+    }
+    return;
+}
+/** Static function to copy GSL matrix data to a psF32 or psF64 image */
+static void matrixGSLtoPS(psImage *out, const gsl_matrix *in)
+{
+    psAssert(in->size1 == out->numRows && in->size2 == out->numCols, "Sizes don't match!");
+    int numCols = out->numCols, numRows = out->numRows; // Size of matrix
+    switch (out->type.type) {
+      case PS_TYPE_F32:
+        for (int y = 0, i = 0; y < numRows; y++) {
+            for (int x = 0; x < numCols; x++, i++) {
+                out->data.F32[y][x] = in->data[i];
+            }
+        }
+        break;
+      case PS_TYPE_F64:
+        if (out->parent || in->tda != in->size1) {
+            for (int y = 0, i = 0; y < numRows; y++, i += numCols) {
+                memcpy(out->data.F64[y], &in->data[i], numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+            }
         } else {
+            outVector->data.F64[i] = inGslVector->data[i];
+        }
+    }
+}
+/** Static function to copy psF32 or psF64 image data to a GSL matrix */
+static void  psImageToGslMatrix(gsl_matrix *outGslMatrix,
+                                const psImage *inImage)
+{
+    psU32 i = 0;
+    psU32 j = 0;
+    psU32 numRows = 0;
+    psU32 numCols = 0;
+    numRows = inImage->numRows;
+    numCols = inImage->numCols;
+    if(inImage->type.type == PS_TYPE_F32) {
+        for(i=0; i<numRows; i++) {
+            for(j=0; j<numCols; j++) {
+                outGslMatrix->data[i*numCols+j] = inImage->data.F32[i][j];
+            }
+        }
+    } else {
+        for(i=0; i<numRows; i++) {
+            for(j=0; j<numCols; j++) {
+                outGslMatrix->data[i*numCols+j] = inImage->data.F64[i][j];
+            }
+        }
+    }
+}
+/** Static function to copy GSL matrix data to a psF32 or psF64 image */
+static void gslMatrixToPsImage(psImage *outImage,
+                               gsl_matrix *inGslMatrix)
+{
+    psU32 i = 0;
+    psU32 j = 0;
+    psU32 numRows = 0;
+    psU32 numCols = 0;
+    numRows = outImage->numRows;
+    numCols = outImage->numCols;
+    if(outImage->type.type == PS_TYPE_F32) {
+        for(i=0; i<numRows; i++) {
+            for(j=0; j<numCols; j++) {
+                outImage->data.F32[i][j] = inGslMatrix->data[i*numCols+j];
+            }
+        }
+    } else {
+        for(i=0; i<numRows; i++) {
+            for(j=0; j<numCols; j++) {
+                outImage->data.F64[i][j] = inGslMatrix->data[i*numCols+j];
+            }
+        }
+    }
+            memcpy(out->p_rawDataBuffer, in->data, numCols * numRows * PSELEMTYPE_SIZEOF(PS_TYPE_F64));
+        }
+        break;
+      default:
+        psAbort("Unsupported vector type: %x\n", out->type.type);
+    }
+    return;
+}
 …
     // Copy psImage data into GSL matrix data
     psImageToGslMatrix(lu, in);
+    matrixPStoGSL(lu, in);
     // Calculate LU decomposition
 …
     // Copy GSL matrix data to psImage data
     gslMatrixToPsImage(out, lu);
+    matrixGSLtoPS(out, lu);
     // Free GSL data
 …
     // Initialize GSL data
     lu = gsl_matrix_alloc(numRows, numCols);
     psImageToGslMatrix(lu, LU);
+    matrixPStoGSL(lu, LU);
     b = gsl_vector_alloc(RHS->n);
     psVectorToGslVector(b, RHS);
+    vectorPStoGSL(b, RHS);
     x = gsl_vector_alloc(RHS->n);
 …
     // Copy GSL vector data to psVector data
     gslVectorToPsVector(out, x);
+    vectorGSLtoPS(out, x);
     // Free GSL data
 …
     // Initialize GSL data
     lu = gsl_matrix_alloc(numRows, numCols);
     psImageToGslMatrix(lu, LU);
+    matrixPStoGSL(lu, LU);
     permGSL.size = perm->n;
 …
     // Copy GSL vector data to psVector data
     gslMatrixToPsImage(out, inverse);
+    matrixGSLtoPS(out, inverse);
     // Free GSL data
 …
     switch (a->type.type) {
       case PS_TYPE_F32: {
           psF32 **values = a->data.F32; /* Dereference */
           int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
           for (int i = 0; i < numRows; i++) {
               for (int j = 0; j < numCols; j++) {
                   if (!isfinite(values[i][j])) {
                       // psError(PS_ERR_BAD_PARAMETER_VALUE, 3,
                       // "Input matrix contains non-finite elements: matrix[%d][%d] is %.2f\n",
                       // i, j, values[i][j]);
                       return false;
+                  }
+              }
+          }
           break;
+          psF32 **values = a->data.F32; /* Dereference */
+          int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
+          for (int i = 0; i < numRows; i++) {
+              for (int j = 0; j < numCols; j++) {
+                  if (!isfinite(values[i][j])) {
+                      // psError(PS_ERR_BAD_PARAMETER_VALUE, 3,
+                      // "Input matrix contains non-finite elements: matrix[%d][%d] is %.2f\n",
+                      // i, j, values[i][j]);
+                      return false;
+                  }
+              }
+          }
+          break;
+      }
       case PS_TYPE_F64: {
           psF64 **values = a->data.F64; /* Dereference */
           int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
           for (int i = 0; i < numRows; i++) {
               for (int j = 0; j < numCols; j++) {
                   if (!isfinite(values[i][j])) {
                       // psError(PS_ERR_BAD_PARAMETER_VALUE, 3,
                       // "Input matrix contains non-finite elements: matrix[%d][%d] is %.2f\n",
                       // i, j, values[i][j]);
                       return false;
+                  }
+              }
+          }
           break;
+          psF64 **values = a->data.F64; /* Dereference */
+          int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
+          for (int i = 0; i < numRows; i++) {
+              for (int j = 0; j < numCols; j++) {
+                  if (!isfinite(values[i][j])) {
+                      // psError(PS_ERR_BAD_PARAMETER_VALUE, 3,
+                      // "Input matrix contains non-finite elements: matrix[%d][%d] is %.2f\n",
+                      // i, j, values[i][j]);
+                      return false;
+                  }
+              }
+          }
+          break;
+      }
         // MATRIX_CHECK_NONFINITE_CASE(F32, a);
         // MATRIX_CHECK_NONFINITE_CASE(F64, a);
+        // MATRIX_CHECK_NONFINITE_CASE(F32, a);
+        // MATRIX_CHECK_NONFINITE_CASE(F64, a);
       default:
         psAbort("Should never get here.");
+        psAbort("Should never get here.");
+    }
 …
     switch (a->type.type) {
       case PS_TYPE_F32: {
           psF32 **values = a->data.F32; /* Dereference */
           int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
           for (int i = 0; i < numRows; i++) {
               for (int j = 0; j < numCols; j++) {
                   if (!isfinite(values[i][j])) {
                       return false;
+                  }
+              }
+          }
           break;
+          psF32 **values = a->data.F32; /* Dereference */
+          int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
+          for (int i = 0; i < numRows; i++) {
+              for (int j = 0; j < numCols; j++) {
+                  if (!isfinite(values[i][j])) {
+                      return false;
+                  }
+              }
+          }
+          break;
+      }
       case PS_TYPE_F64: {
           psF64 **values = a->data.F64; /* Dereference */
           int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
           for (int i = 0; i < numRows; i++) {
               for (int j = 0; j < numCols; j++) {
                   if (!isfinite(values[i][j])) {
                       return false;
+                  }
+              }
+          }
           break;
+          psF64 **values = a->data.F64; /* Dereference */
+          int numCols = a->numCols, numRows = a->numRows; /* Size of matrix */
+          for (int i = 0; i < numRows; i++) {
+              for (int j = 0; j < numCols; j++) {
+                  if (!isfinite(values[i][j])) {
+                      return false;
+                  }
+              }
+          }
+          break;
+      }
       default:
         psAbort("Should never get here.");
+    }
+        psAbort("Should never get here.");
+    }
     // Following the algorithm laid out by Press et al., we loop along the matrix diagonal, but
     // we do not operate on the diagonal elements in order.  Instead, we are looking for the
 …
     if (a->type.type == PS_TYPE_F32) {
         psF32 **A = a->data.F32;
         psF32  *B = b->data.F32;
         int *colIndex = colIndexV->data.S32;
         int *rowIndex = rowIndexV->data.S32;
         int *pivot    = pivotV->data.S32;
         psF32 growth = 1.0;
         for (int diag = 0; diag < nSquare; diag++) {
             psF32 maxval = 0.0;
             int maxrow = 0;
             int maxcol = 0;
             // search for the next pivot
             for (int row = 0; row < nSquare; row++) {
                 if (!isfinite(A[row][diag])) goto escape;
                 // if we have already operated on this row (pivot[row] is true), skip it
                 if (pivot[row]) continue;
                 // if we have not yet operated on this row (pivot[row] is false), look for pivot for this row
                 for (int col = 0; col < nSquare; col++) {
                     if (pivot[col]) continue;
                     if (fabs (A[row][col]) < maxval) continue;
                     maxval = fabs (A[row][col]);
                     maxrow = row;
                     maxcol = col;
+                }
+            }
             // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix
             if (pivot[maxcol]) goto escape;
             pivot[maxcol] = 1;
             // if the selected pivot is off the diagonal, do a row swap
             if (maxrow != maxcol) {
                 for (int col = 0; col < nSquare; col++) PS_SWAP (A[maxrow][col], A[maxcol][col]);
                 PS_SWAP (B[maxrow], B[maxcol]);
+            }
             rowIndex[diag] = maxrow;
             colIndex[diag] = maxcol;
             if (A[maxcol][maxcol] == 0.0) goto escape;
             // Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow.
             // Here we are going to raise an error if the dynamic range is too large.
             /* rescale by pivot reciprocal */
             psF32 tmpval = 1.0 / A[maxcol][maxcol];
             A[maxcol][maxcol] = 1.0;
             for (int col = 0; col < nSquare; col++) A[maxcol][col] *= tmpval;
             B[maxcol] *= tmpval;
             // check for ill-conditioned matrix: measure the pivot growth and trigger on over/under flow
             growth *= tmpval;
             psTrace ("psLib.math", 4, "growth : %e\n", growth);
             if (fabs(growth) > MAX_RANGE) goto escape;
             /* adjust the elements above the pivot */
             for (int row = 0; row < nSquare; row++) {
                 if (row == maxcol) continue;
                 tmpval = A[row][maxcol];
                 A[row][maxcol] = 0.0;
                 for (int col = 0; col < nSquare; col++) A[row][col] -= A[maxcol][col]*tmpval;
                 B[row] -= B[maxcol]*tmpval;
+            }
+        }
         // swap back the inverse matrix based on the row swaps above
         for (int col = nSquare - 1; col >= 0; col--) {
             if (rowIndex[col] != colIndex[col]) {
                 for (int row = 0; row < nSquare; row++) PS_SWAP (A[row][rowIndex[col]], A[row][colIndex[col]]);
+            }
+        }
+        psF32 **A = a->data.F32;
+        psF32  *B = b->data.F32;
+        int *colIndex = colIndexV->data.S32;
+        int *rowIndex = rowIndexV->data.S32;
+        int *pivot    = pivotV->data.S32;
+        psF32 growth = 1.0;
+        for (int diag = 0; diag < nSquare; diag++) {
+            psF32 maxval = 0.0;
+            int maxrow = 0;
+            int maxcol = 0;
+            // search for the next pivot
+            for (int row = 0; row < nSquare; row++) {
+                if (!isfinite(A[row][diag])) goto escape;
+                // if we have already operated on this row (pivot[row] is true), skip it
+                if (pivot[row]) continue;
+                // if we have not yet operated on this row (pivot[row] is false), look for pivot for this row
+                for (int col = 0; col < nSquare; col++) {
+                    if (pivot[col]) continue;
+                    if (fabs (A[row][col]) < maxval) continue;
+                    maxval = fabs (A[row][col]);
+                    maxrow = row;
+                    maxcol = col;
+                }
+            }
+            // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix
+            if (pivot[maxcol]) goto escape;
+            pivot[maxcol] = 1;
+            // if the selected pivot is off the diagonal, do a row swap
+            if (maxrow != maxcol) {
+                for (int col = 0; col < nSquare; col++) PS_SWAP (A[maxrow][col], A[maxcol][col]);
+                PS_SWAP (B[maxrow], B[maxcol]);
+            }
+            rowIndex[diag] = maxrow;
+            colIndex[diag] = maxcol;
+            if (A[maxcol][maxcol] == 0.0) goto escape;
+            // Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow.
+            // Here we are going to raise an error if the dynamic range is too large.
+            /* rescale by pivot reciprocal */
+            psF32 tmpval = 1.0 / A[maxcol][maxcol];
+            A[maxcol][maxcol] = 1.0;
+            for (int col = 0; col < nSquare; col++) A[maxcol][col] *= tmpval;
+            B[maxcol] *= tmpval;
+            // check for ill-conditioned matrix: measure the pivot growth and trigger on over/under flow
+            growth *= tmpval;
+            psTrace ("psLib.math", 4, "growth : %e\n", growth);
+            if (fabs(growth) > MAX_RANGE) goto escape;
+            /* adjust the elements above the pivot */
+            for (int row = 0; row < nSquare; row++) {
+                if (row == maxcol) continue;
+                tmpval = A[row][maxcol];
+                A[row][maxcol] = 0.0;
+                for (int col = 0; col < nSquare; col++) A[row][col] -= A[maxcol][col]*tmpval;
+                B[row] -= B[maxcol]*tmpval;
+            }
+        }
+        // swap back the inverse matrix based on the row swaps above
+        for (int col = nSquare - 1; col >= 0; col--) {
+            if (rowIndex[col] != colIndex[col]) {
+                for (int row = 0; row < nSquare; row++) PS_SWAP (A[row][rowIndex[col]], A[row][colIndex[col]]);
+            }
+        }
     } else {
         psF64 **A = a->data.F64;
         psF64  *B = b->data.F64;
         int *colIndex = colIndexV->data.S32;
         int *rowIndex = rowIndexV->data.S32;
         int *pivot    = pivotV->data.S32;
         psF64 growth = 1.0;
         for (int diag = 0; diag < nSquare; diag++) {
             psF64 maxval = 0.0;
             int maxrow = 0;
             int maxcol = 0;
             // search for the next pivot
             for (int row = 0; row < nSquare; row++) {
                 if (!isfinite(A[row][diag])) goto escape;
                 // if we have already operated on this row (pivot[row] is true), skip it
                 if (pivot[row]) continue;
                 // if we have not yet operated on this row (pivot[row] is false), look for pivot for this row
                 for (int col = 0; col < nSquare; col++) {
                     if (pivot[col]) continue;
                     if (fabs (A[row][col]) < maxval) continue;
                     maxval = fabs (A[row][col]);
                     maxrow = row;
                     maxcol = col;
+                }
+            }
             // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix
             if (pivot[maxcol]) goto escape;
             pivot[maxcol] = 1;
             // if the selected pivot is off the diagonal, do a row swap
             if (maxrow != maxcol) {
                 for (int col = 0; col < nSquare; col++) PS_SWAP (A[maxrow][col], A[maxcol][col]);
                 PS_SWAP (B[maxrow], B[maxcol]);
+            }
             rowIndex[diag] = maxrow;
             colIndex[diag] = maxcol;
             if (A[maxcol][maxcol] == 0.0) goto escape;
             // Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow.
             // Here we are going to raise an error if the dynamic range is too large.
             /* rescale by pivot reciprocal */
             psF64 tmpval = 1.0 / A[maxcol][maxcol];
             A[maxcol][maxcol] = 1.0;
             for (int col = 0; col < nSquare; col++) A[maxcol][col] *= tmpval;
             B[maxcol] *= tmpval;
             // check for ill-conditioned matrix: measure the pivot growth and trigger on over/under flow
             growth *= tmpval;
             psTrace ("psLib.math", 4, "growth : %e\n", growth);
             if (fabs(growth) > MAX_RANGE) goto escape;
             /* adjust the elements above the pivot */
             for (int row = 0; row < nSquare; row++) {
                 if (row == maxcol) continue;
                 tmpval = A[row][maxcol];
                 A[row][maxcol] = 0.0;
                 for (int col = 0; col < nSquare; col++) A[row][col] -= A[maxcol][col]*tmpval;
                 B[row] -= B[maxcol]*tmpval;
+            }
+        }
         // swap back the inverse matrix based on the row swaps above
         for (int col = nSquare - 1; col >= 0; col--) {
             if (rowIndex[col] != colIndex[col]) {
                 for (int row = 0; row < nSquare; row++) PS_SWAP (A[row][rowIndex[col]], A[row][colIndex[col]]);
+            }
+        }
+        psF64 **A = a->data.F64;
+        psF64  *B = b->data.F64;
+        int *colIndex = colIndexV->data.S32;
+        int *rowIndex = rowIndexV->data.S32;
+        int *pivot    = pivotV->data.S32;
+        psF64 growth = 1.0;
+        for (int diag = 0; diag < nSquare; diag++) {
+            psF64 maxval = 0.0;
+            int maxrow = 0;
+            int maxcol = 0;
+            // search for the next pivot
+            for (int row = 0; row < nSquare; row++) {
+                if (!isfinite(A[row][diag])) goto escape;
+                // if we have already operated on this row (pivot[row] is true), skip it
+                if (pivot[row]) continue;
+                // if we have not yet operated on this row (pivot[row] is false), look for pivot for this row
+                for (int col = 0; col < nSquare; col++) {
+                    if (pivot[col]) continue;
+                    if (fabs (A[row][col]) < maxval) continue;
+                    maxval = fabs (A[row][col]);
+                    maxrow = row;
+                    maxcol = col;
+                }
+            }
+            // if pivot[maxcol] is set, we have already done this row: this implies a singular matrix
+            if (pivot[maxcol]) goto escape;
+            pivot[maxcol] = 1;
+            // if the selected pivot is off the diagonal, do a row swap
+            if (maxrow != maxcol) {
+                for (int col = 0; col < nSquare; col++) PS_SWAP (A[maxrow][col], A[maxcol][col]);
+                PS_SWAP (B[maxrow], B[maxcol]);
+            }
+            rowIndex[diag] = maxrow;
+            colIndex[diag] = maxcol;
+            if (A[maxcol][maxcol] == 0.0) goto escape;
+            // Kahan replaces the 0.0 pivot with epsilon*(largest element in column) + underFlow.
+            // Here we are going to raise an error if the dynamic range is too large.
+            /* rescale by pivot reciprocal */
+            psF64 tmpval = 1.0 / A[maxcol][maxcol];
+            A[maxcol][maxcol] = 1.0;
+            for (int col = 0; col < nSquare; col++) A[maxcol][col] *= tmpval;
+            B[maxcol] *= tmpval;
+            // check for ill-conditioned matrix: measure the pivot growth and trigger on over/under flow
+            growth *= tmpval;
+            psTrace ("psLib.math", 4, "growth : %e\n", growth);
+            if (fabs(growth) > MAX_RANGE) goto escape;
+            /* adjust the elements above the pivot */
+            for (int row = 0; row < nSquare; row++) {
+                if (row == maxcol) continue;
+                tmpval = A[row][maxcol];
+                A[row][maxcol] = 0.0;
+                for (int col = 0; col < nSquare; col++) A[row][col] -= A[maxcol][col]*tmpval;
+                B[row] -= B[maxcol]*tmpval;
+            }
+        }
+        // swap back the inverse matrix based on the row swaps above
+        for (int col = nSquare - 1; col >= 0; col--) {
+            if (rowIndex[col] != colIndex[col]) {
+                for (int row = 0; row < nSquare; row++) PS_SWAP (A[row][rowIndex[col]], A[row][colIndex[col]]);
+            }
+        }
+    }
 …
     lu = gsl_matrix_alloc(numRows, numCols);
     inv = gsl_matrix_alloc(numRows, numCols);
     psImageToGslMatrix(lu, in);
+    matrixPStoGSL(lu, in);
     // Invert data and calculate determinant
 …
     if (determinant) {
       // XXX this is getting the wrong value: is it the wrong calculation?
       // it disagrees with the results of
+      // it disagrees with the results of
       // det = (psF32)gsl_linalg_LU_det(lu, signum);
       // used in psMatrixDeterminatn
 …
     // Copy GSL matrix data to psImage data
     gslMatrixToPsImage(out, inv);
+    matrixGSLtoPS(out, inv);
     // Free GSL structs
 …
     perm = gsl_permutation_alloc(numRows);
     lu = gsl_matrix_alloc(numRows, numCols);
     psImageToGslMatrix(lu, in);
+    matrixPStoGSL(lu, in);
     // Calculate determinant
     gsl_linalg_LU_decomp(lu, perm, &signum);
     det = (psF32)gsl_linalg_LU_det(lu, signum);
+    det = (psF32)gsl_linalg_LU_lndet(lu);
     // Free GSL structs
 …
     inGSL = gsl_matrix_alloc(numRows, numCols);
     psImageToGslMatrix(inGSL, in);
+    matrixPStoGSL(inGSL, in);
     outGSL = gsl_matrix_alloc(numRows, numCols);
 …
     // Copy GSL matrix data to psImage data
     gslMatrixToPsImage(out, outGSL);
+    matrixGSLtoPS(out, outGSL);
     // Free GSL structs
 …
+}
+psVector *psMatrixSolveSVD(psVector *out, const psImage *matrix, const psVector *vector, float thresh)
+{
+    #define psMatrixSolveSVD_EXIT {psFree(out); return NULL; }
+    PS_ASSERT_GENERAL_IMAGE_NON_NULL(matrix, psMatrixSolveSVD_EXIT);
+    PS_CHECK_DIMEN_AND_TYPE(matrix, PS_DIMEN_IMAGE, psMatrixSolveSVD_EXIT);
+    PS_ASSERT_GENERAL_VECTOR_NON_NULL(vector, psMatrixSolveSVD_EXIT);
+    PS_CHECK_DIMEN_AND_TYPE(vector, PS_DIMEN_VECTOR, psMatrixSolveSVD_EXIT);
+    int numCols = matrix->numCols, numRows = matrix->numRows; // Size of matrix
+    // Decompose matrix: A = U S V^T
+    gsl_matrix *A = gsl_matrix_alloc(numRows, numCols); // Input matrix in GSL-speak; becomes matrix U
+    gsl_matrix *V = gsl_matrix_alloc(numCols, numCols); // Untransposed matrix V
+    gsl_vector *S = gsl_vector_alloc(numCols);          // Singular values
+    gsl_vector *work = gsl_vector_alloc(numCols);       // Work space for GSL
+    matrixPStoGSL(A, matrix);
+    int gslStatus = 0;                  // Status of GSL
+    if ((gslStatus = gsl_linalg_SV_decomp(A, V, S, work))) {
+        const char *err = gsl_strerror(gslStatus);
+        psError(PS_ERR_UNKNOWN, true, "Unable to decompose matrix: %s", err);
+        gsl_matrix_free(A);
+        gsl_matrix_free(V);
+        gsl_vector_free(S);
+        gsl_vector_free(work);
+        return NULL;
+    }
+    gsl_vector_free(work);
+    if (isfinite(thresh) && thresh > 0.0) {
+        // Trim the singular values
+        double total = 0.0;             // Total of singular values
+        for (int i = 0; i < numCols; i++) {
+            total += gsl_vector_get(S, i);
+        }
+        thresh *= total;
+        for (int i = 0; i < numCols; i++) {
+            double value = gsl_vector_get(S, i); // Singular value
+            if (value < thresh) {
+                psTrace("psLib.math", 5, "Trimming singular value %d: %lg", i, value);
+                gsl_vector_set(S, i, 0.0);
+#if 0
+                for (int j = 0; j < numCols; j++) {
+                    // Being thorough; probably unnecessary
+                    gsl_matrix_set(V, j, i, 0.0);
+                    gsl_matrix_set(A, j, i, 0.0);
+                }
+#endif
+            } else {
+                psTrace("psLib.math", 5, "Singular value %d: %lg", i, value);
+            }
+        }
+    }
+    // Solve system (or minimise least-squares if overconstrained): Ax = b
+    gsl_vector *b = gsl_vector_alloc(numCols); // Vector b
+    gsl_vector *x = gsl_vector_alloc(numCols); // Solution
+    vectorPStoGSL(b, vector);
+    if ((gslStatus = gsl_linalg_SV_solve(A, V, S, b, x))) {
+        const char *err = gsl_strerror(gslStatus);
+        psError(PS_ERR_UNKNOWN, true, "Unable to solve matrix equation: %s", err);
+        gsl_matrix_free(A);
+        gsl_matrix_free(V);
+        gsl_vector_free(S);
+        gsl_vector_free(b);
+        gsl_vector_free(x);
+        return NULL;
+    }
+    gsl_matrix_free(A);
+    gsl_matrix_free(V);
+    gsl_vector_free(S);
+    gsl_vector_free(b);
+    out = psVectorRecycle(out, numCols, PS_TYPE_F64);
+    vectorGSLtoPS(out, x);
+    gsl_vector_free(x);
+    return out;
+}
 // This code supplied by Andy Becker (becker@astro.washington.edu)
 psImage *psMatrixSVD(psImage* evec, psVector* eval, const psImage* in)
+psImage *psMatrixSVD_old(psImage* evec, psVector* eval, const psImage* in)
+{
     #define psMatrixSVD_EXIT {psFree(evec); psFree(eval); return NULL;}
 …
     // Copy psImage data into GSL matrix data
     psImageToGslMatrix(A, in);
+    matrixPStoGSL(A, in);
     // Calculate SVD decomposition
 …
     // Copy GSL matrix data to psImage data
     gslMatrixToPsImage(evec, V);
     gslVectorToPsVector(eval, S);
+    matrixGSLtoPS(evec, V);
+    vectorGSLtoPS(eval, S);
     // Take the square root of eval
 …
     return evec;
+}
+// this is basically a wrapper for the gsl function: gsl_linalg_SV_decomp() SVD decomposes
+// matrix A based on the following equation: A = U w V^T .  This function (as usual for SVD
+// implementations) returns V not V^T.  U and V are returned to images; w is returned to a
+// vector representing the diagonal of w.  The input image A is not modified.  U, V, and w may
+// be supplied as NULL or may be allocated; their lengths are set here to match the
+// dimensionality of A.  XXX there is no error handling for the gsl functions (anywhere in
+// psMatrix.c)
+bool psMatrixSVD(psImage **U, psVector **w, psImage **V, const psImage *A)
+{
+    // Error checks  Missing one for eval
+    PS_ASSERT_PTR_NON_NULL(U, false);
+    PS_ASSERT_PTR_NON_NULL(w, false);
+    PS_ASSERT_PTR_NON_NULL(V, false);
+    PS_ASSERT_PTR_NON_NULL(A, false);
+    // A is provided with size Nx,Ny = numCols,numRows
+    // U has size Nx,Ny
+    // V has size Nx,Nx
+    // w has size Nx
+    // Initialize data
+    int numRows = A->numRows;
+    int numCols = A->numCols;
+    *U = psImageRecycle(*U,  numCols, numRows, A->type.type);
+    *V = psImageRecycle(*V,  numCols, numCols, A->type.type);
+    *w = psVectorRecycle(*w, numCols, A->type.type);
+    gsl_matrix *Agsl = gsl_matrix_alloc(numRows, numCols);
+    gsl_matrix *Vgsl = gsl_matrix_alloc(numCols, numCols);
+    gsl_vector *Sgsl = gsl_vector_alloc(numCols);
+    gsl_vector *work = gsl_vector_alloc(numCols);
+    // Copy psImage data into GSL matrix data
+    matrixPStoGSL(Agsl, A);
+    // Calculate SVD decomposition
+    gsl_linalg_SV_decomp(Agsl, Vgsl, Sgsl, work);
+    // Copy GSL matrix data to psImage data
+    matrixGSLtoPS(*V, Vgsl);
+    matrixGSLtoPS(*U, Agsl);  // gsl_linalg_SV_decomp replaces A with U
+    vectorGSLtoPS(*w, Sgsl);
+    // Free GSL data
+    gsl_matrix_free(Agsl);
+    gsl_matrix_free(Vgsl);
+    gsl_vector_free(Sgsl);
+    gsl_vector_free(work);
+    return true;
+}

branches/tap_branches/psLib/src/math/psMatrix.h

-              r24084
+              r27838
  */
 psImage *psMatrixLUInvert(
     psImage *out,                  ///< place result here if not NULL
     const psImage* LU,             ///< LU-decomposed matrix.
     const psVector* perm           ///< Permutation vector resulting from psMatrixLUD function.
+    psImage *out,                  ///< place result here if not NULL
+    const psImage* LU,             ///< LU-decomposed matrix.
+    const psVector* perm           ///< Permutation vector resulting from psMatrixLUD function.
 );
 …
 );
+/// Single value decomposition, provided by Andy Becker
+psImage *psMatrixSVD(psImage* evec, psVector* eval, const psImage* in);
+/// Solve a matrix equation using Singular Value Decomposition
+///
+/// Solves Ax = b for x
+psVector *psMatrixSolveSVD(
+    psVector *solution,                 ///< Solution to output, or NULL
+    const psImage *matrix,              ///< Matrix to be solved
+    const psVector *vector,             ///< Vector of values
+    float thresh                        ///< Threshold relative to maximum for trimming singular values
+    );
+/// Single value decomposition (original by Andy Becker, updated by EAM)
+bool psMatrixSVD(psImage **U, psVector **w, psImage **V, const psImage *A);
 /// @}

branches/tap_branches/psLib/src/math/psMinimizeLMM.c

-              r24089
+              r27838
 // XXX check that the GJ solver works:
 # if (TESTGJ)
+# if (TESTGJ)
     psImage *out = psImageAlloc (alpha->numRows, alpha->numCols, PS_TYPE_F32);
     for (int oy = 0; oy < out->numRows; oy++) {
         for (int ox = 0; ox < out->numCols; ox++) {
             float value = 0;
             for (int i = 0; i < alpha->numCols; i++) {
                 value += alpha->data.F32[i][ox]*Alpha->data.F32[oy][i];
+            }
             out->data.F32[oy][ox] = value;
+        }
+        for (int ox = 0; ox < out->numCols; ox++) {
+            float value = 0;
+            for (int i = 0; i < alpha->numCols; i++) {
+                value += alpha->data.F32[i][ox]*Alpha->data.F32[oy][i];
+            }
+            out->data.F32[oy][ox] = value;
+        }
+    }
     psVector *vect = psVectorAlloc (beta->n, PS_TYPE_F32);
     for (int oy = 0; oy < vect->n; oy++) {
         float value = 0;
         for (int i = 0; i < alpha->numCols; i++) {
             value += alpha->data.F32[oy][i]*Beta->data.F32[i];
+        }
         vect->data.F32[oy] = value;
+    }
+        float value = 0;
+        for (int i = 0; i < alpha->numCols; i++) {
+            value += alpha->data.F32[oy][i]*Beta->data.F32[i];
+        }
+        vect->data.F32[oy] = value;
+    }
     psFree (out);
     psFree (vect);
 …
     if (isnan(chisq)) {
         psTrace ("psLib.math", 5, "psMinLM_SetABX() returned a NAN chisq.\n");
         psVectorInit (delta, NAN);
+        psVectorInit (delta, NAN);
         retValue = false;
+    }
 …
     if (!status) {
         psTrace ("psLib.math", 5, "psMinLM_GuessABP() returned FALSE.\n");
         psVectorInit (delta, NAN);
+        psVectorInit (delta, NAN);
         retValue = false;
+    }
 …
     PS_ASSERT_VECTOR_NON_NULL(dy, NAN);
     PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, false);
+    PS_ASSERT_VECTOR_TYPE(params, PS_TYPE_F32, NAN);
     if (paramMask) {
         PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, false);
+        PS_ASSERT_VECTOR_TYPE(paramMask, PS_TYPE_VECTOR_MASK, NAN);
+    }
 …
         chisq += PS_SQR(delta) * dy->data.F32[i];
         assert (!isnan(dy->data.F32[i]));
         assert (!isnan(delta));
         assert (!isnan(chisq));
+        if (isnan(dy->data.F32[i])) return NAN;
+        if (isnan(delta)) return NAN;
+        if (isnan(chisq)) return NAN;
         // we track alpha,beta and params,deriv separately

branches/tap_branches/psLib/src/math/psStats.c

-              r25884
+              r27838
     // Iterate to get the best bin size; an iteration limit is enforced at the bottom of the loop.
     for (int iterate = 1; iterate > 0; iterate++) {
+        psTrace(TRACE, 6,
+                "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n",
+                iterate);
+        psTrace(TRACE, 6, "-------------------- Iterating on Bin size.  Iteration number %d --------------------\n", iterate);
+        if (iterate >= PS_ROBUST_MAX_ITERATIONS) {
+          // This occurs when a large number of the values are identical --- a bin size cannot be found
+          // that will spread out the distribution.  Therefore, set what we can, and fall over
+          // gracefully.
+          COUNT_WARNING(10, 100, "Maximum number of iterations (%d) exceeded.", PS_ROBUST_MAX_ITERATIONS);
+          goto escape;
+        }
         // Get the minimum and maximum values
 …
         psTrace(TRACE, 6, "Initial robust bin size is %.2f\n", binSize);
+        // ADD step 0: Construct the histogram with the specified bin size.  NOTE: we can not specify the bin
+        // size precisely since the argument to psHistogramAlloc() is the number of bins, not the binSize.  If
+        // we get here, we know that binSize != 0.0.
+        long numBins = (max - min) / binSize; // Number of bins
+        // ADD step 0: Construct the histogram with the specified bin size.  NOTE: we can
+        // not specify the bin size precisely since the argument to psHistogramAlloc() is
+        // the number of bins, not the binSize.  If we get here, we know that binSize !=
+        // 0.0.  We can also have a floating-point round-off error such that the last bin
+        // of the histogram does not correspond exactly with the value of 'max'.  Let's be
+        // a bit generous and extend the histogram by two bins in either direction
+        long numBins = 4 + (max - min) / binSize; // Number of bins
         psTrace(TRACE, 6, "Numbins is %ld\n", numBins);
         psTrace(TRACE, 6, "Creating a robust histogram from data range (%.2f - %.2f)\n", min, max);
         // Generate the histogram
         histogram = psHistogramAlloc(min, max, numBins);
+        histogram = psHistogramAlloc(min - 2.0*binSize, max + 2.0*binSize, numBins);
         // XXXXX we need to consider this step if errors -> variance
         if (!psVectorHistogram(histogram, myVector, errors, mask, maskVal)) {
 …
             psFree(statsMinMax);
             psFree(mask);
             return false;
+        }
 …
             PS_VECTOR_PRINT_F32(histogram->bounds);
             PS_VECTOR_PRINT_F32(histogram->nums);
+        }
+        // perversity check: if most of the values land in a single bin, then we probably
+        // have a perverse case (eg, small number of points at extremely large / small
+        // values; nearly bi-modal distribution).  if so, keep only points within 5? 10?
+        // bins of that excess bin:
+        int nMaxBin = 0;
+        int iMaxBin = 0;
+        for (long i = 1; i < histogram->nums->n; i++) {
+            if (histogram->nums->data.F32[i] > nMaxBin) {
+                nMaxBin = histogram->nums->data.F32[i];
+                iMaxBin = i;
+            }
+        }
+        if (nMaxBin > numValid / 2) {
+            float minKeep = histogram->bounds->data.F32[iMaxBin] - 10*binSize;
+            float maxKeep = histogram->bounds->data.F32[iMaxBin + 1] + 10*binSize;
+            int nInvalid = 0;
+            for (long i = 0; i < myVector->n; i++) {
+                // skip the already-masked values
+                if (mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & maskVal) continue;
+                bool invalid = false;
+                invalid |= (myVector->data.F32[i] <= minKeep);
+                invalid |= (myVector->data.F32[i] >= maxKeep);
+                invalid |= (!isfinite(myVector->data.F32[i]));
+                if (!invalid) continue;
+                mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = maskVal;
+                nInvalid ++;
+            }
+            if (nInvalid) {
+              psTrace(TRACE, 6, "data is concentrated in a single bin, masking %d extreme outliers and retrying\n", nInvalid);
+              psFree(histogram);
+              psFree(cumulative);
+              histogram = NULL;
+              cumulative = NULL;
+              continue;
+            }
+            // if we did not mask anything, give up.
+        }
 …
     stats->robustN50 = N50;
     psTrace(TRACE, 6, "The robustN50 is %ld.\n", N50);
+    psTrace(TRACE, 6, "The robust median and stdev are %f, %f\n", stats->robustMedian, stats->robustStdev);
     // Clean up
 …
             COUNT_WARNING(10, 100, "Failed to calculate the min/max of the input vector.\n");
             psFree(statsMinMax);
+            psFree(histogram);
             goto escape;
+        }
 …
             if (!status) {
                 psErrorClear();
+                psErrorClear();
                 COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
                 psFree(poly);
 …
             if (!status) {
                 psErrorClear();
+                psErrorClear();
                 COUNT_WARNING(10, 100, "Failed to fit a gaussian to the robust histogram.\n");
                 psFree(poly);

branches/tap_branches/psLib/src/math/psUnaryOp.c

r17050	r27838
211	211	}
212	212
213		psMathType* psUnaryOp(psPtr out, ~~const~~ psPtr in, const char *op)
	213	psMathType* psUnaryOp(psPtr out, psPtr in, const char *op)
214	214	{
215	215	#define psUnaryOp_EXIT { \

branches/tap_branches/psLib/src/math/psUnaryOp.h

r11248	r27838
59	59	psMathType* psUnaryOp(
60	60	psPtr out, ///< Output type, either psImage or psVector.
61		~~const psPtr in,~~ ///< Input, either psImage or psVector.
	61	psPtr in, ///< Input, either psImage or psVector.
62	62	const char *op ///< Operator.
63	63	);

branches/tap_branches/psLib/src/mathtypes/psImage.h

r19056	r27838
66	66	#define P_PSIMAGE_SET_ROW0(img,r0) {(int)&img->row0 = r0;}
67	67	#define P_PSIMAGE_SET_TYPE(img,t) {(psMathType)&img->type = t;}
	68	#define P_PSIMAGE_GET_TYPE(img) ((img)->type->type)
68	69
69	70	/** Create an image of the specified size and type.

branches/tap_branches/psLib/src/mathtypes/psVector.c

r24886	r27838
729	729	char line[1024];
730	730
731		sprintf (line, "vector: %s\n", name);
	731	sprintf (line, "# vector: %s\n", name);
732	732	if (write(fd, line, strlen(line))) {;} //ignore return value
733	733

branches/tap_branches/psLib/src/pslib_strict.h

r23149	r27838
102	102	#include "psType.h"
103	103	#include "psArray.h"
104		#include "psBit~~Set~~.h"
	104	#include "psBits.h"
105	105	#include "psHash.h"
106	106	#include "psList.h"

branches/tap_branches/psLib/src/sys/psErrorCodes.c.in

r11675	r27838
67	67	static void freeErrorDescription(psErrorDescription* err)
68	68	{
69		psFree(~~(void *)~~err->description);
	69	psFree(err->description);
70	70	}
71	71

branches/tap_branches/psLib/src/sys/psMemory.h

-              r23305
+              r27838
 /** Free memory.  This operates much like free().
+ *
+ *
  *  @see psAlloc, psRealloc
+ *  note: we cast ptr to (void *) in case we are supplied a const pointer.
  */
 #ifdef DOXYGEN
 …
 #ifndef SWIG
 #define psFree(ptr) \
         psMemDecrRefCounter(ptr)
+    ptr = psMemDecrRefCounter((void *)ptr);
 #endif // ifndef SWIG
 #endif // ifdef DOXYGEN

branches/tap_branches/psLib/src/sys/psTrace.c

r20546	r27838
119	119
120	120	psMemSetPersistent((psPtr)comp->name,false);
121		psFree(~~(void *)~~comp->name);
	121	psFree(comp->name);
122	122	}
123	123

branches/tap_branches/psLib/src/sys/psType.c

-              r11617
+              r27838
 #include "psType.h"
 #include "psBitSet.h"
+#include "psBits.h"
 #include "psFits.h"
 #include "psPixels.h"
 …
+        }
         break;
     case PS_DATA_BITSET:
         if (psMemCheckBitSet(ptr)) {
+    case PS_DATA_BITS:
+        if (psMemCheckBits(ptr)) {
             return true;
+        }

branches/tap_branches/psLib/src/sys/psType.h

r25256	r27838
107	107	PS_DATA_STRING = 0x10000, ///< psString (char *)
108	108	PS_DATA_ARRAY, ///< psArray
109		PS_DATA_BITS~~ET, ///< psBitSet~~
	109	PS_DATA_BITS, ///< psBits
110	110	PS_DATA_CUBE, ///< psCube
111	111	PS_DATA_FITS, ///< psFits

branches/tap_branches/psLib/src/types/Makefile.am

-              r23148
+              r27838
 libpslibtypes_la_SOURCES = \
         psArray.c \
         psBitSet.c \
+        psBits.c \
         psHash.c \
         psList.c \
 …
 pkginclude_HEADERS = \
         psArray.h \
         psBitSet.h \
+        psBits.h \
         psHash.h \
         psList.h \

branches/tap_branches/psLib/src/types/psArray.c

r15714	r27838
166	166	// drop an item from the array and free it
167	167	bool psArrayRemoveData(psArray* array,
168		~~const~~ psPtr data)
	168	psPtr data)
169	169	{
170	170	PS_ASSERT_ARRAY_NON_NULL(array, false);

branches/tap_branches/psLib/src/types/psArray.h

r19056	r27838
186	186	bool psArrayRemoveData(
187	187	psArray* array, ///< array to operate on
188		~~const psPtr data~~ ///< the data pointer to remove from psArray
	188	psPtr data ///< the data pointer to remove from psArray
189	189	);
190	190

branches/tap_branches/psLib/src/types/psList.c

-              r18955
+              r27838
+    }
+    // XXX remove this as an error
     if (location < 0 || location >= (int)list->n) {
         psError(PS_ERR_BAD_PARAMETER_VALUE, true,
 …
     psListIterator* iterator = list->iterators->data[0];
+    // XXX remove this as an eror
     if (! psListIteratorSet(iterator,location)) {
         psError(PS_ERR_BAD_PARAMETER_VALUE, true,

branches/tap_branches/psLib/src/types/psLookupTable.c

-              r17447
+              r27838
 #include "psString.h"
 #include "psError.h"
+#include "psString.h"
+#include "psSlurp.h"
 #include "psLookupTable.h"
 …
         char *end = NULL; \
         ps##TYPE value = FUNC(strValue, &end, 0); \
+        if (*end != '\0' && !isspace(*end)) { \
+        if (*end != '\0' && *end != '\n' && !isspace(*end)) { \
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Characters left over after parsing %s: %s", \
+                strValue, end); \
             *status = PS_PARSE_ERROR_VALUE; \
         } \
 …
         char *end = NULL; \
         ps##TYPE value = FUNC(strValue, &end); \
+        if (*end != '\0' && !isspace(*end)) { \
+        if (*end != '\0' && *end != '\n' && !isspace(*end)) { \
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Characters left over after parsing %s: %s", \
+                strValue, end); \
             *status = PS_PARSE_ERROR_VALUE; \
         } \
 …
+}
+psArray *psVectorsReadFromFile(const char *filename,
+                               const char *format)
+psArray *psVectorsReadFromFile(const char *filename, const char *format)
+{
     PS_ASSERT_STRING_NON_EMPTY(filename, NULL);
     PS_ASSERT_STRING_NON_EMPTY(format, NULL);
+    psArray*          outputArray = NULL;
+    psVector*         colVector   = NULL;
+    char*             strValue    = NULL;
+    char*             strNum      = NULL;
+    char*             line        = NULL;
+    char*             linePtr     = NULL;
+    int               numCols     = 0;
+    int               numRows     = 0;
+    FILE*             fp          = NULL;
+    const char*       tempFormat  = NULL;
+    psParseErrorType  parseStatus = PS_PARSE_SUCCESS;
+    // Create temp pointer which can then be used several times
+    tempFormat = format;
+    // Create output array and set array elements to zero
+    outputArray = psArrayAllocEmpty(INITIAL_NUM);
+    // Parse the format string to determine how many vectors
+    // and whether the format string is valid
+    while ((strValue = getToken((char**)&tempFormat, " \t", &parseStatus))) {
+        // Check for %d format sub string
+    psArray *outputArray = psArrayAllocEmpty(INITIAL_NUM); // Array of vectors to return
+    psParseErrorType parseStatus = PS_PARSE_SUCCESS; // Status of parsing
+    // Parse the format string to determine how many vectors and whether the format string is valid
+    const char *tempFormat = format;    // Pointer into format
+    psString strValue;                  // Format of interest
+    int numCols = 0;                    // Number of columns found in format
+    while ((strValue = getToken((char**)&tempFormat, " \t", &parseStatus)) &&
+           parseStatus == PS_PARSE_SUCCESS) {
+        if (strstr(strValue,"\%*") != 0) {
+            // Don't increase number of columns
+            continue;
+        }
+        psElemType type;                // Type specified
         if (strcmp(strValue,"\%d") == 0 ) {
+            numCols++;
+            colVector = psVectorAlloc(1,PS_TYPE_S32);
+            outputArray = psArrayAdd(outputArray, ARRAY_STRIDE, colVector);
+            psFree(colVector);
+            type = PS_TYPE_S32;
         } else if (strcmp(strValue,"\%ld") == 0) {
+            numCols++;
+            colVector = psVectorAlloc(1,PS_TYPE_S64);
+            outputArray = psArrayAdd(outputArray, ARRAY_STRIDE, colVector);
+            psFree(colVector);
+            type = PS_TYPE_S64;
         } else if (strcmp(strValue,"\%f") == 0) {
+            numCols++;
+            colVector = psVectorAlloc(1,PS_TYPE_F32);
+            outputArray = psArrayAdd(outputArray, ARRAY_STRIDE, colVector);
+            psFree(colVector);
+            type = PS_TYPE_F32;
         } else if (strcmp(strValue,"\%lf") == 0) {
+            numCols++;
+            colVector = psVectorAlloc(1,PS_TYPE_F64);
+            outputArray = psArrayAdd(outputArray, ARRAY_STRIDE, colVector);
+            psFree(colVector);
+        } else if (strstr(strValue,"\%*") != 0) {
+            // Don't increase number of columns
+            type = PS_TYPE_F64;
         } else {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true,
+                    "Invalid format specifier");
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Invalid format specifier: %s", strValue);
             psFree(strValue);
+            numCols = 0;
+            break;
+        }
+            psFree(outputArray);
+            return NULL;
+        }
+        psVector *colVector = psVectorAllocEmpty(1, type); // Vector for type
+        outputArray = psArrayAdd(outputArray, ARRAY_STRIDE, colVector);
+        psFree(colVector);
+        numCols++;
         psFree(strValue);
+    }
+    if (parseStatus != PS_PARSE_SUCCESS) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Failed to parse format at column %d: %s",
+                numCols, strValue);
+        psFree(strValue);
+        psFree(outputArray);
+        return NULL;
+    }
+    if (numCols == 0) {
+        // Format string parse error detected
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Format string was not parsed sucessfully");
+        psFree(outputArray);
+        return NULL;
+    }
     // If the format string was parsed successfully and return numCols the
     // prepare to open file and read values
+    if (numCols > 0) {
+        // Open specified file
+        if ((fp=fopen(filename, "r")) == NULL) {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, _("Failed to open file %s."),
+                    filename);
+            psFree(outputArray);
+            return NULL;
+        } else {
+            // Initialize array index
+            int arrayIndex = 0;
+            // Create reusable line for continuous read
+            line = (char*)psAlloc(MAX_STRING_LENGTH*sizeof(char));
+            // Loop through file to get numRows, numCols, and column data types
+            while ((fgets(line, MAX_STRING_LENGTH, fp) != NULL) &&
+                    (parseStatus == PS_PARSE_SUCCESS))   {
+                // Copy pointer to line for parsing
+                linePtr = line;
+                // If line is not a comment or blank, then extract data
+                if (!ignoreLine(linePtr)) {
+                    numRows++;
+                    // Copy format pointer for parsing
+                    tempFormat = format;
+                    arrayIndex = 0;
+                    parseStatus = PS_PARSE_SUCCESS;
+                    // Loop through format and line strings to get values in text table file
+                    while ((strValue=getToken((char**)&tempFormat," \t",&parseStatus))
+                            && (strNum=getToken((char**)&linePtr," \t",&parseStatus)) ) {
+                        // Set column vector
+                        colVector = outputArray->data[arrayIndex];
+                        // Set column entries based on format string defining the type
+                        if (strcmp(strValue,"\%d") == 0 ) {
+                            colVector = psVectorRecycle(colVector, numRows,
+                                                        colVector->type.type);
+                            parseValue(colVector,numRows-1,strNum,&parseStatus);
+                            arrayIndex++;
+                        } else if (strcmp(strValue,"\%ld") == 0) {
+                            colVector = psVectorRecycle(colVector, numRows,
+                                                        colVector->type.type);
+                            parseValue(colVector,numRows-1,strNum,&parseStatus);
+                            arrayIndex++;
+                        } else if (strcmp(strValue,"\%f") == 0) {
+                            colVector = psVectorRecycle(colVector, numRows,
+                                                        colVector->type.type);
+                            parseValue(colVector,numRows-1,strNum,&parseStatus);
+                            arrayIndex++;
+                        } else if (strcmp(strValue,"\%lf") == 0) {
+                            colVector = psVectorRecycle(colVector, numRows,
+                                                        colVector->type.type);
+                            parseValue(colVector,numRows-1,strNum,&parseStatus);
+                            arrayIndex++;
+                        } else if (strstr(strValue,"\%*") != 0) {
+                            // Don't increase number of columns
+                        }
+                        psFree(strValue);
+                        psFree(strNum);
+                        // If the file line was not parsed successful report
+                        // error and return NULL
+                        if (parseStatus != PS_PARSE_SUCCESS) {
+                            psError(PS_ERR_UNKNOWN, true,
+                                    "Parsing text file failed.");
+                            fclose(fp);
+                            psFree(outputArray);
+                            psFree(line);
+                            return NULL;
+                        }
+                    }
+                    if (strValue != NULL && strNum == NULL) {
+                        psError(PS_ERR_UNKNOWN, true,
+                                "Parsing text file failed - missing table value(s).");
+                        fclose(fp);
+                        psFree(outputArray);
+                        psFree(line);
+                        psFree(strValue);
+                        return NULL;
+                    }
+                }  // ignore line
+    psString file = psSlurpFilename(filename); // Contents of file
+    if (!file) {
+        psError(psErrorCodeLast(), false, "Unable to read file of vectors");
+        psFree(outputArray);
+        return NULL;
+    }
+    psArray *lines = psStringSplitArray(file, "\n", false); // Lines of file
+    psFree(file);
+    long numRows = 0;                                  // Number of rows
+    for (long i = 0; i < lines->n; i++) {
+        psString line = lines->data[i]; // Line of interest
+        if (ignoreLine(line)) {
+            continue;
+        }
+        numRows++;
+        char *linePtr = line;           // Pointer into line
+        // Copy format pointer for parsing
+        const char *tempFormat = format; // Pointer into format
+        long arrayIndex = 0;            // Index in array
+        parseStatus = PS_PARSE_SUCCESS;
+        // Loop through format and line strings to get values in text table file
+        char *strNum;                   // Number within line
+        while ((strValue=getToken((char**)&tempFormat," \t",&parseStatus)) &&
+               (strNum=getToken((char**)&linePtr," \t",&parseStatus)) &&
+               parseStatus == PS_PARSE_SUCCESS) {
+            if (strstr(strValue,"\%*") != 0) {
+                continue;
+            }
+            //Return NULL for an empty table
+            if (numRows == 0) {
+                psError(PS_ERR_UNKNOWN, true,
+                        "Parsing text file failed - input table is empty.");
+                fclose(fp);
+            // Set column vector
+            psVector *colVector = outputArray->data[arrayIndex]; // Column vector of interest
+            outputArray->data[arrayIndex] = colVector = psVectorRecycle(colVector, numRows,
+                                                                        colVector->type.type);
+            parseValue(colVector, numRows - 1, strNum, &parseStatus);
+            arrayIndex++;
+            if (parseStatus != PS_PARSE_SUCCESS) {
+                psError(PS_ERR_UNKNOWN, false, "Parsing text file failed: %s as %s", strNum, strValue);
                 psFree(outputArray);
+                psFree(line);
+                psFree(lines);
+                psFree(strNum);
+                psFree(strValue);
                 return NULL;
+            }
+            // Read on the lines in the file - close file pointer
+            fclose(fp);
+            psFree(line);
+        }
+    } else {
+        // Format string parse error detected
+        psError(PS_ERR_UNKNOWN, true,
+                "Format string was not parsed sucessfully");
+            psFree(strValue);
+            psFree(strNum);
+        }
+        if (strValue != NULL && strNum == NULL) {
+            psError(PS_ERR_UNKNOWN, true,
+                    "Parsing text file failed - missing table value(s).");
+            psFree(outputArray);
+            psFree(lines);
+            psFree(strValue);
+            return NULL;
+        }
+    }
+    if (parseStatus != PS_PARSE_SUCCESS) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Failed to parse format at column %d: %s",
+                numCols, strValue);
+        psFree(strValue);
         psFree(outputArray);
         return NULL;
+    }
+    psFree(lines);
     // Return populated array

branches/tap_branches/psLib/src/types/psMetadata.c

-              r25383
+              r27838
+    }
     case     PS_DATA_ARRAY:                     // psArray
     case     PS_DATA_BITSET:                    // psBitSet
+    case     PS_DATA_BITS:                      // psBits
     case     PS_DATA_CUBE:                      // psCube
     case     PS_DATA_FITS:                      // psFits
 …
 // may need to extend this to change the keyname in the copy
+bool psMetadataItemSupplement(psMetadata *out,
+bool psMetadataItemSupplement(bool *status,
+                              psMetadata *out,
                               const psMetadata *in,
                               const char *key)
 …
     psMetadataItem *item = psMetadataLookup(in, key);
     if (!item) {
+        psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Could not find '%s' in metadata.\n", key);
+        if (status) {
+            *status = false;
+        } else {
+            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Could not find '%s' in metadata.\n", key);
+        }
         return false;
+    }
 …
+}
+psMetadataItem* psMetadataLookup(const psMetadata *md,
+// Get entry by index
+psMetadataItem *psMetadataGetIndex(psMetadata *md, long location)
+{
+    PS_ASSERT_METADATA_NON_NULL(md, false);
+    return psListGet(md->list, location);
+}
+psMetadataItem *psMetadataLookup(const psMetadata *md,
                                  const char *key)
+{
 …
     PS_ASSERT_METADATA_NON_NULL(md,NULL);
+    // XXX remove this as an error
     entry = (psMetadataItem*) psListGet(md->list, location);
     if (entry == NULL) {

branches/tap_branches/psLib/src/types/psMetadata.h

r25383	r27838
545	545	*/
546	546	bool psMetadataItemSupplement(
	547	bool *status, ///< if supplied, returns true/false if key is found (suppresses the error)
547	548	psMetadata *out, ///< output Metadata container for copying.
548	549	const psMetadata *in, ///< Metadata collection from which to copy.

branches/tap_branches/psLib/src/types/psMetadataConfig.c

-              r23859
+              r27838
         return false;
+    }
+    fprintf(file, "%s", fileString);
+    if (fprintf(file, "%s", fileString) != strlen(fileString)) {
+        psError(PS_ERR_IO, true, "Failed to write contents of configuration file %s", filename);
+        psFree(fileString);
+        fclose(file);
+        return false;
+    }
     psFree(fileString);
     if (fclose(file) == EOF) {

branches/tap_branches/psLib/test/types/Makefile.am

r18145	r27838
29	29	tap_psPixels_all \
30	30	tap_psHash_all \
31		tap_psBit~~Set~~_all \
	31	tap_psBits_all \
32	32	tap_psList_all \
33	33	tap_psLookupTable_all \

Context Navigation

Legend:

Download in other formats: