IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26892


Ignore:
Timestamp:
Feb 10, 2010, 7:27:29 PM (16 years ago)
Author:
eugene
Message:

updates from eam_branches/20091201

Location:
trunk/psLib/src
Files:
20 edited

Legend:

Unmodified
Added
Removed

  • trunk/psLib/src/imageops/psImageConvolve.c

    r25383 r26892  
    246246    return kernel;
    247247}
     248
     249bool psKernelTruncate(psKernel *kernel, float frac)
     250{
     251    PS_ASSERT_KERNEL_NON_NULL(kernel, false);
     252    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(frac, 0.0, false);
     253    PS_ASSERT_FLOAT_LESS_THAN(frac, 1.0, false);
     254
     255    if (frac == 0.0) {
     256        // Nothing to do
     257        return true;
     258    }
     259
     260    int xMin = kernel->xMin, xMax = kernel->xMax, yMin = kernel->yMin, yMax = kernel->yMax; // Bounds
     261    int maxSize = PS_MAX(PS_MAX(PS_MAX(-xMin, xMax), -yMin), yMax); // Maximum size
     262
     263    // Determine the threshold
     264    // Summing absolute values because large negative deviations have power as well
     265    double sumKernel = 0.0;             // Sum of the kernel
     266    for (int y = yMin; y <= yMax; y++) {
     267        for (int x = xMin; x <= xMax; x++) {
     268            sumKernel += fabsf(kernel->kernel[y][x]);
     269        }
     270    }
     271
     272    float threshold = sumKernel * (1.0 - frac); // Threshold for truncation
     273
     274    // Find truncation size
     275    int truncate = 0;                     // Truncation radius
     276    for (int radius = 1; truncate == 0 && radius < maxSize; radius++) {
     277        int uMin = PS_MAX(-radius, xMin);
     278        int uMax = PS_MIN(radius, xMax);
     279        int vMin = PS_MAX(-radius, yMin);
     280        int vMax = PS_MIN(radius, yMax);
     281        int r2 = PS_SQR(radius);
     282        double sum = 0.0;
     283        for (int v = vMin; v <= vMax; v++) {
     284            int v2 = PS_SQR(v);
     285            for (int u = uMin; u <= uMax; u++) {
     286                int u2 = PS_SQR(u);
     287                if (u2 + v2 <= r2) {
     288                    sum += fabsf(kernel->kernel[v][u]);
     289                }
     290            }
     291        }
     292        if (sum > threshold) {
     293            // This is the truncation radius
     294            truncate = radius;
     295        }
     296    }
     297    if (truncate == maxSize) {
     298        // No truncation possible
     299        return true;
     300    }
     301
     302    // Truncate the kernel
     303    {
     304        int uMin = PS_MAX(-truncate, xMin);
     305        int uMax = PS_MIN(truncate, xMax);
     306        int vMin = PS_MAX(-truncate, yMin);
     307        int vMax = PS_MIN(truncate, yMax);
     308        int r2 = PS_SQR(truncate);
     309        for (int v = vMin; v <= vMax; v++) {
     310            int v2 = PS_SQR(v);
     311            for (int u = uMin; u <= uMax; u++) {
     312                int u2 = PS_SQR(u);
     313                if (u2 + v2 > r2) {
     314                    kernel->kernel[v][u] = 0.0;
     315                }
     316            }
     317        }
     318    }
     319    kernel->xMin = PS_MAX(-truncate, kernel->xMin);
     320    kernel->xMax = PS_MIN(truncate, kernel->xMax);
     321    kernel->yMin = PS_MAX(-truncate, kernel->yMin);
     322    kernel->yMax = PS_MIN(truncate, kernel->yMax);
     323
     324    return true;
     325    }
     326
     327
    248328
    249329psImage *psImageConvolveDirect(psImage *out, const psImage *in, const psKernel *kernel)

  • trunk/psLib/src/imageops/psImageConvolve.h

    r25383 r26892  
    138138);
    139139
     140/// Truncate a kernel
     141///
     142/// Truncates the outer parts of the kernel where the contribution is below the nominated fraction of the
     143/// total kernel.
     144bool psKernelTruncate(
     145    psKernel *in,                       ///< Kernel to be truncated
     146    float frac                          ///< Fraction for truncation threshold
     147    );
     148
     149
    140150/// Convolve an image with a kernel, using a direct convolution
    141151///

  • trunk/psLib/src/imageops/psImageCovariance.c

    r25994 r26892  
    1414#include "psTrace.h"
    1515#include "psBinaryOp.h"
     16#include "psScalar.h"
     17#include "psThread.h"
    1618
    1719#include "psImageCovariance.h"
     20
     21static bool threaded = false;           // Run threaded?
     22
    1823
    1924psKernel *psImageCovarianceNone(void)
    … …  
    2429}
    2530
    26 
    27 psKernel *psImageCovarianceCalculate(const psKernel *kernel, const psKernel *covariance)
    28 {
    29     PS_ASSERT_KERNEL_NON_NULL(kernel, NULL);
    30 
    31     // See http://en.wikipedia.org/wiki/Error_propagation
    32     //
    33     // If
    34     //     f_k = sum_i A_ik x_i
    35     // is a set of functions, then the covariance matrix for f is given by:
    36     //     M^f_ij = sum_k sum_l A_ik M^x_kl A_lj
    37     // where M^x is the covariance matrix for x.
    38     // Note that the errors in f are correlated (covariance) even if the errors in x are not.
    39 
    40     psKernel *covar;                    // Covariance matrix to use
    41     if (covariance) {
    42         covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const
    43     } else {
    44         covar = psImageCovarianceNone();
    45     }
    46 
    47     // Check for non-finite elements
    48     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
    49         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
    50             if (!isfinite(kernel->kernel[y][x])) {
    51                 psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    52                         "Non-finite covariance matrix element in kernel at %d,%d", x, y);
    53                 psFree(covar);
    54                 return NULL;
    55             }
    56         }
    57     }
    58     for (int y = covar->yMin; y <= covar->yMax; y++) {
    59         for (int x = covar->xMin; x <= covar->xMax; x++) {
    60             if (!isfinite(covar->kernel[y][x])) {
    61                 psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    62                         "Non-finite covariance matrix element in covariance matrix at %d,%d", x, y);
    63                 psFree(covar);
    64                 return NULL;
    65             }
    66         }
    67     }
    68 
    69     // The above (double) sum for the covariance matrix means that, for each point in the output covariance
    70     // matrix, we need to work out all combinations of getting to the central point via a kernel, input
    71     // covariance matrix and another kernel.  This means that the resultant covariance matrix has twice the
    72     // size of the kernel plus the size of the input covariance matrix.
    73     int xMin = kernel->xMin - kernel->xMax + covar->xMin, xMax = kernel->xMax - kernel->xMin + covar->xMax;
    74     int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax;
    75     psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output
     31/// Calculation of covariance matrix element when convolving
     32static float imageCovarianceCalculate(const psKernel *covar, // Original covariance matrix
     33                                      const psKernel *kernel, // Convolution kernel
     34                                      int x, int y            // Coordinates in output covariance matrix
     35                                      )
     36{
     37    psAssert(covar, "Require covariance matrix");
     38    psAssert(kernel, "Require kernel");
    7639
    7740    // Need to go:
    … …  
    8952    // from the source coordinate), we take the smallest possible (because everything else is zero outside).
    9053
    91     double total = 0.0;                 // Total covariance
     54    // Range for v
     55    int vMin = PS_MAX(kernel->yMin + covar->yMin, y + kernel->yMin);
     56    int vMax = PS_MIN(kernel->yMax + covar->yMax, y + kernel->yMax);
     57    // Range for u
     58    int uMin = PS_MAX(kernel->xMin + covar->xMin, x + kernel->xMin);
     59    int uMax = PS_MIN(kernel->xMax + covar->xMax, x + kernel->xMax);
     60
     61    double sum = 0.0;           // Sum for value of covariance matrix at (x,y)
     62    for (int v = vMin; v <= vMax; v++) {
     63        // Range for q
     64        int qMin = PS_MAX(v + covar->yMin, kernel->yMin);
     65        int qMax = PS_MIN(v + covar->yMax, kernel->yMax);
     66        for (int u = uMin; u <= uMax; u++) {
     67            // Range for p
     68            int pMin = PS_MAX(u + covar->xMin, kernel->xMin);
     69            int pMax = PS_MIN(u + covar->xMax, kernel->xMax);
     70
     71            double xyuvValue = kernel->kernel[v-y][u-x]; // Value for (x,y) --> (u,v)
     72
     73            double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0)
     74            for (int q = qMin; q <= qMax; q++) {
     75                for (int p = pMin; p <= pMax; p++) {
     76                    uvpqValue += (double)covar->kernel[q-v][p-u] * (double)kernel->kernel[q][p];
     77                }
     78            }
     79            sum += xyuvValue * uvpqValue;
     80        }
     81    }
     82
     83    return sum;
     84}
     85
     86/// Thread entry point for calculation of covariance matrix element when convolving
     87static bool imageCovarianceCalculateThread(psThreadJob *job)
     88{
     89    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     90    psAssert(job->args, "No job arguments");
     91    psAssert(job->args->n == 5, "Wrong number of job arguments: %ld", job->args->n);
     92
     93    psKernel *out = job->args->data[0]; // Output covariance matrix
     94    const psKernel *covar = job->args->data[1]; // Input covariance matrix
     95    const psKernel *kernel = job->args->data[2]; // Convolution kernel
     96    int x = PS_SCALAR_VALUE(job->args->data[3], S32); // x coordinate in output covariance matrix
     97    int y = PS_SCALAR_VALUE(job->args->data[4], S32); // y coordinate in output covariance matrix
     98
     99    out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
     100
     101    return true;
     102}
     103
     104
     105
     106psKernel *psImageCovarianceCalculate(const psKernel *kernel, const psKernel *covariance)
     107{
     108    PS_ASSERT_KERNEL_NON_NULL(kernel, NULL);
     109
     110    // See http://en.wikipedia.org/wiki/Error_propagation
     111    //
     112    // If
     113    //     f_k = sum_i A_ik x_i
     114    // is a set of functions, then the covariance matrix for f is given by:
     115    //     M^f_ij = sum_k sum_l A_ik M^x_kl A_lj
     116    // where M^x is the covariance matrix for x.
     117    // Note that the errors in f are correlated (covariance) even if the errors in x are not.
     118
     119    psKernel *covar;                    // Covariance matrix to use
     120    if (covariance) {
     121        covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const
     122    } else {
     123        covar = psImageCovarianceNone();
     124    }
     125
     126    // Check for non-finite elements
     127    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
     128        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
     129            if (!isfinite(kernel->kernel[y][x])) {
     130                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     131                        "Non-finite covariance matrix element in kernel at %d,%d", x, y);
     132                psFree(covar);
     133                return NULL;
     134            }
     135        }
     136    }
     137    for (int y = covar->yMin; y <= covar->yMax; y++) {
     138        for (int x = covar->xMin; x <= covar->xMax; x++) {
     139            if (!isfinite(covar->kernel[y][x])) {
     140                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     141                        "Non-finite covariance matrix element in covariance matrix at %d,%d", x, y);
     142                psFree(covar);
     143                return NULL;
     144            }
     145        }
     146    }
     147
     148    // The above (double) sum for the covariance matrix means that, for each point in the output covariance
     149    // matrix, we need to work out all combinations of getting to the central point via a kernel, input
     150    // covariance matrix and another kernel.  This means that the resultant covariance matrix has twice the
     151    // size of the kernel plus the size of the input covariance matrix.
     152    int xMin = kernel->xMin - kernel->xMax + covar->xMin, xMax = kernel->xMax - kernel->xMin + covar->xMax;
     153    int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax;
     154    psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output
     155
    92156    for (int y = yMin; y <= yMax; y++) {
    93         // Range for v
    94         int vMin = PS_MAX(kernel->yMin + covar->yMin, y + kernel->yMin);
    95         int vMax = PS_MIN(kernel->yMax + covar->yMax, y + kernel->yMax);
    96157        for (int x = xMin; x <= xMax; x++) {
    97             // Range for u
    98             int uMin = PS_MAX(kernel->xMin + covar->xMin, x + kernel->xMin);
    99             int uMax = PS_MIN(kernel->xMax + covar->xMax, x + kernel->xMax);
    100 
    101             double sum = 0.0;           // Sum for value of covariance matrix at (x,y)
    102             for (int v = vMin; v <= vMax; v++) {
    103                 // Range for q
    104                 int qMin = PS_MAX(v + covar->yMin, kernel->yMin);
    105                 int qMax = PS_MIN(v + covar->yMax, kernel->yMax);
    106                 for (int u = uMin; u <= uMax; u++) {
    107                     // Range for p
    108                     int pMin = PS_MAX(u + covar->xMin, kernel->xMin);
    109                     int pMax = PS_MIN(u + covar->xMax, kernel->xMax);
    110 
    111                     double xyuvValue = kernel->kernel[v-y][u-x]; // Value for (x,y) --> (u,v)
    112 
    113                     double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0)
    114                     for (int q = qMin; q <= qMax; q++) {
    115                         for (int p = pMin; p <= pMax; p++) {
    116                             uvpqValue += (double)covar->kernel[q-v][p-u] * (double)kernel->kernel[q][p];
    117                         }
    118                     }
    119                     sum += xyuvValue * uvpqValue;
     158            if (threaded) {
     159                psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE");
     160                psArrayAdd(job->args, 1, out);
     161                psArrayAdd(job->args, 1, covar);
     162                psArrayAdd(job->args, 1, (psKernel*)kernel); // Casting away const
     163                PS_ARRAY_ADD_SCALAR(job->args, x, PS_TYPE_S32);
     164                PS_ARRAY_ADD_SCALAR(job->args, y, PS_TYPE_S32);
     165                if (!psThreadJobAddPending(job)) {
     166                    psFree(covar);
     167                    return NULL;
    120168                }
    121             }
    122             out->kernel[y][x] = sum;
    123             total += sum;
    124         }
    125     }
    126     psTrace("psLib.imageops", 3, "Total covariance: %lf ; Central variance: %f\n", total, out->kernel[0][0]);
    127 
     169                psFree(job);
     170            } else {
     171                out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
     172            }
     173        }
     174    }
    128175    psFree(covar);
     176
     177    if (threaded && !psThreadPoolWait(true)) {
     178        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
     179        return false;
     180    }
     181
    129182    return out;
    130183}
     184
     185float psImageCovarianceCalculateFactor(const psKernel *kernel, const psKernel *covariance)
     186{
     187    psKernel *covar;                    // Covariance matrix to use
     188    if (covariance) {
     189        covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const
     190    } else {
     191        covar = psImageCovarianceNone();
     192    }
     193
     194    // Check for non-finite elements
     195    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
     196        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
     197            if (!isfinite(kernel->kernel[y][x])) {
     198                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     199                        "Non-finite covariance matrix element in kernel at %d,%d", x, y);
     200                psFree(covar);
     201                return NAN;
     202            }
     203        }
     204    }
     205    for (int y = covar->yMin; y <= covar->yMax; y++) {
     206        for (int x = covar->xMin; x <= covar->xMax; x++) {
     207            if (!isfinite(covar->kernel[y][x])) {
     208                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     209                        "Non-finite covariance matrix element in covariance matrix at %d,%d", x, y);
     210                psFree(covar);
     211                return NAN;
     212            }
     213        }
     214    }
     215
     216    float factor = imageCovarianceCalculate(covar, kernel, 0, 0); // Covariance factor
     217    psFree(covar);
     218    return factor;
     219}
     220
     221// Calculation of covariance matrix element when binning
     222static float imageCovarianceBin(const psKernel *covar, // Original covariance matrix
     223                                int bin,               // Binning factor
     224                                float binVal,          // Convolution kernel value for binning
     225                                int x, int y           // Coordinates in output covariance matrix
     226    )
     227{
     228    psAssert(covar, "Require covariance matrix");
     229    psAssert(bin > 0 && binVal > 0, "Require binning: %d %f", bin, binVal);
     230
     231    int binMin = -(bin - 1) / 2, binMax = bin / 2; // Range of "kernel"
     232
     233    // Range for v
     234    int vMin = PS_MAX(binMin + covar->yMin, bin * y + binMin);
     235    int vMax = PS_MIN(binMax + covar->yMax, bin * y + binMax);
     236    // Range for u
     237    int uMin = PS_MAX(binMin + covar->xMin, bin * x + binMin);
     238    int uMax = PS_MIN(binMax + covar->xMax, bin * x + binMax);
     239
     240    double sum = 0.0;           // Sum for value of covariance matrix at (x,y)
     241    for (int v = vMin; v <= vMax; v++) {
     242        // Range for q
     243        int qMin = PS_MAX(v + covar->yMin, binMin);
     244        int qMax = PS_MIN(v + covar->yMax, binMax);
     245        for (int u = uMin; u <= uMax; u++) {
     246            // Range for p
     247            int pMin = PS_MAX(u + covar->xMin, binMin);
     248            int pMax = PS_MIN(u + covar->xMax, binMax);
     249
     250            double xyuvValue = binVal; // Value for (x,y) --> (u,v)
     251
     252            double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0)
     253            for (int q = qMin; q <= qMax; q++) {
     254                for (int p = pMin; p <= pMax; p++) {
     255                    uvpqValue += (double)covar->kernel[q-v][p-u] * (double)binVal;
     256                }
     257            }
     258            sum += xyuvValue * uvpqValue;
     259        }
     260    }
     261
     262    return sum;
     263}
     264
     265/// Thread entry point for calculation of covariance matrix element when binning
     266static bool imageCovarianceBinThread(psThreadJob *job)
     267{
     268    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     269    psAssert(job->args, "No job arguments");
     270    psAssert(job->args->n == 6, "Wrong number of job arguments: %ld", job->args->n);
     271
     272    psKernel *out = job->args->data[0]; // Output covariance matrix
     273    const psKernel *covar = job->args->data[1]; // Input covariance matrix
     274    int bin = PS_SCALAR_VALUE(job->args->data[2], S32); // Binning factor
     275    float binVal = PS_SCALAR_VALUE(job->args->data[3], F32); // Convolution kernel value for binning
     276    int x = PS_SCALAR_VALUE(job->args->data[4], S32); // x coordinate in output covariance matrix
     277    int y = PS_SCALAR_VALUE(job->args->data[5], S32); // y coordinate in output covariance matrix
     278
     279    out->kernel[y][x] = imageCovarianceBin(covar, bin, binVal, x, y);
     280
     281    return true;
     282}
     283
    131284
    132285psKernel *psImageCovarianceBin(int bin, const psKernel *covariance, bool average)
    … …  
    165318    psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output
    166319
    167     double total = 0.0;                 // Total covariance
    168320    for (int y = yMin; y <= yMax; y++) {
    169         // Range for v
    170         int vMin = PS_MAX(binMin + covar->yMin, bin * y + binMin);
    171         int vMax = PS_MIN(binMax + covar->yMax, bin * y + binMax);
    172321        for (int x = xMin; x <= xMax; x++) {
    173             // Range for u
    174             int uMin = PS_MAX(binMin + covar->xMin, bin * x + binMin);
    175             int uMax = PS_MIN(binMax + covar->xMax, bin * x + binMax);
    176 
    177             double sum = 0.0;           // Sum for value of covariance matrix at (x,y)
    178             for (int v = vMin; v <= vMax; v++) {
    179                 // Range for q
    180                 int qMin = PS_MAX(v + covar->yMin, binMin);
    181                 int qMax = PS_MIN(v + covar->yMax, binMax);
    182                 for (int u = uMin; u <= uMax; u++) {
    183                     // Range for p
    184                     int pMin = PS_MAX(u + covar->xMin, binMin);
    185                     int pMax = PS_MIN(u + covar->xMax, binMax);
    186 
    187                     double xyuvValue = binVal; // Value for (x,y) --> (u,v)
    188 
    189                     double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0)
    190                     for (int q = qMin; q <= qMax; q++) {
    191                         for (int p = pMin; p <= pMax; p++) {
    192                             uvpqValue += (double)covar->kernel[q-v][p-u] * (double)binVal;
    193                         }
    194                     }
    195                     sum += xyuvValue * uvpqValue;
     322            if (threaded) {
     323                psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_COVARIANCE_BIN");
     324                psArrayAdd(job->args, 1, out);
     325                psArrayAdd(job->args, 1, covar);
     326                PS_ARRAY_ADD_SCALAR(job->args, bin, PS_TYPE_S32);
     327                PS_ARRAY_ADD_SCALAR(job->args, binVal, PS_TYPE_F32);
     328                PS_ARRAY_ADD_SCALAR(job->args, x, PS_TYPE_S32);
     329                PS_ARRAY_ADD_SCALAR(job->args, y, PS_TYPE_S32);
     330                if (!psThreadJobAddPending(job)) {
     331                    psFree(covar);
     332                    return NULL;
    196333                }
    197             }
    198             out->kernel[y][x] = sum;
    199             total += sum;
    200         }
    201     }
    202     psTrace("psLib.imageops", 3, "Total covariance: %lf ; Central variance: %f\n", total, out->kernel[0][0]);
    203 
     334                psFree(job);
     335            } else {
     336                out->kernel[y][x] = imageCovarianceBin(covar, bin, binVal, x, y);
     337            }
     338        }
     339    }
    204340    psFree(covar);
    205341
     342    if (threaded && !psThreadPoolWait(true)) {
     343        psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
     344        return false;
     345    }
     346
    206347    return out;
    207348}
    … …  
    219360
    220361    for (int y = covar->yMin; y <= covar->yMax; y++) {
    221         if (y < -radius) continue;
    222         if (y > +radius) continue;
    223         for (int x = covar->xMin; x <= covar->xMax; x++) {
    224             if (x < -radius) continue;
    225             if (x > +radius) continue;
    226 
    227             if (hypot(x, y) > radius) continue;
    228 
    229             psAssert (isfinite(covar->kernel[y][x]), "invalid NAN in covariance matrix");
    230             Sum += covar->kernel[y][x];
    231         }
     362        if (y < -radius) continue;
     363        if (y > +radius) continue;
     364        for (int x = covar->xMin; x <= covar->xMax; x++) {
     365            if (x < -radius) continue;
     366            if (x > +radius) continue;
     367
     368            if (hypot(x, y) > radius) continue;
     369
     370            psAssert (isfinite(covar->kernel[y][x]), "invalid NAN in covariance matrix");
     371            Sum += covar->kernel[y][x];
     372        }
    232373    }
    233374
    … …  
    428569    return true;
    429570}
     571
     572
     573bool psImageCovarianceSetThreads(bool set)
     574{
     575    bool old = threaded;                // Old value
     576    if (set && !threaded) {
     577        {
     578            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 5);
     579            task->function = &imageCovarianceCalculateThread;
     580            psThreadTaskAdd(task);
     581            psFree(task);
     582        }
     583        {
     584            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_BIN", 6);
     585            task->function = &imageCovarianceBinThread;
     586            psThreadTaskAdd(task);
     587            psFree(task);
     588        }
     589    } else if (!set && threaded) {
     590        psThreadTaskRemove("PSLIB_IMAGE_COVARIANCE_CALCULATE");
     591        psThreadTaskRemove("PSLIB_IMAGE_COVARIANCE_BIN");
     592    }
     593    threaded = set;
     594    return old;
     595}
     596
     597bool psImageCovarianceGetThreads(void)
     598{
     599    return threaded;
     600}

  • trunk/psLib/src/imageops/psImageCovariance.h

    r25994 r26892  
    4747    );
    4848
     49/// Return the pixel-to-pixel covariance factor following calculation
     50///
     51/// This doesn't require calculation of the entire covariance matrix, so is much faster.
     52float psImageCovarianceCalculateFactor(
     53    const psKernel *kernel,             ///< Convolution kernel
     54    const psKernel *covariance          ///< Current covariance pseudo-matrix
     55    );
     56
     57
    4958/// Return the covariance factor for an aperture of a given radius
    5059float psImageCovarianceFactorForAperture(const psKernel *covar, float radius);
    … …  
    8190    );
    8291
     92/// Control threading for image covariance functions
     93///
     94/// Returns old threading status
     95bool psImageCovarianceSetThreads(bool threaded ///< Run image covariance functions threaded?
     96    );
     97
     98/// Return whether image covariance functions are threaded
     99bool psImageCovarianceGetThreads(void);
    83100
    84101/// @}

  • trunk/psLib/src/imageops/psImageInterpolate.c

    r23231 r26892  
    197197{
    198198    // Casting away const
    199     psFree((psImage*)interp->image);
    200     psFree((psImage*)interp->mask);
    201     psFree((psImage*)interp->variance);
    202     psFree((psImage*)interp->kernel);
    203     psFree((psImage*)interp->kernel2);
    204     psFree((psVector*)interp->sumKernel2);
     199    psFree(interp->image);
     200    psFree(interp->mask);
     201    psFree(interp->variance);
     202    psFree(interp->kernel);
     203    psFree(interp->kernel2);
     204    psFree(interp->sumKernel2);
    205205}
    206206

  • trunk/psLib/src/math/psBinaryOp.c

    r17050 r26892  
    379379}
    380380
    381 psMathType* psBinaryOp(psPtr out, const psPtr in1, const char *op, const psPtr in2)
     381psMathType* psBinaryOp(psPtr out, psPtr in1, const char *op, psPtr in2)
    382382{
    383383

  • trunk/psLib/src/math/psBinaryOp.h

    r11248 r26892  
    5555psMathType* psBinaryOp(
    5656    psPtr out,                         ///< Output type, either psImage or psVector.
    57     const psPtr in1,                   ///< First input, either psImage or psVector.
     57    psPtr in1,                   ///< First input, either psImage or psVector.
    5858    const char *op,                    ///< Operator.
    59     const psPtr in2                    ///< Second input, either psImage or psVector.
     59    psPtr in2                    ///< Second input, either psImage or psVector.
    6060);
    6161

  • trunk/psLib/src/math/psHistogram.c

    r21183 r26892  
    127127static void histogramFree(psHistogram* myHist)
    128128{
    129     psFree((void *)myHist->bounds);
     129    psFree(myHist->bounds);
    130130    psFree(myHist->nums);
    131131}

  • trunk/psLib/src/math/psMatrix.c

    r26001 r26892  
    10331033}
    10341034
    1035 psVector *psMatrixSolveSVD(psVector *out, const psImage *matrix, const psVector *vector)
     1035psVector *psMatrixSolveSVD(psVector *out, const psImage *matrix, const psVector *vector, float thresh)
    10361036{
    10371037    #define psMatrixSolveSVD_EXIT {psFree(out); return NULL; }
    … …  
    10631063    gsl_vector_free(work);
    10641064
     1065    if (isfinite(thresh) && thresh > 0.0) {
     1066        // Trim the singular values
     1067        double total = 0.0;             // Total of singular values
     1068        for (int i = 0; i < numCols; i++) {
     1069            total += gsl_vector_get(S, i);
     1070        }
     1071        thresh *= total;
     1072        for (int i = 0; i < numCols; i++) {
     1073            double value = gsl_vector_get(S, i); // Singular value
     1074            if (value < thresh) {
     1075                psTrace("psLib.math", 5, "Trimming singular value %d: %lg", i, value);
     1076                gsl_vector_set(S, i, 0.0);
     1077#if 0
     1078                for (int j = 0; j < numCols; j++) {
     1079                    // Being thorough; probably unnecessary
     1080                    gsl_matrix_set(V, j, i, 0.0);
     1081                    gsl_matrix_set(A, j, i, 0.0);
     1082                }
     1083#endif
     1084            } else {
     1085                psTrace("psLib.math", 5, "Singular value %d: %lg", i, value);
     1086            }
     1087        }
     1088    }
     1089
    10651090    // Solve system (or minimise least-squares if overconstrained): Ax = b
    10661091    gsl_vector *b = gsl_vector_alloc(numCols); // Vector b
    … …  
    10931118}
    10941119
    1095 
    1096 
    10971120// This code supplied by Andy Becker (becker@astro.washington.edu)
    1098 psImage *psMatrixSVD(psImage* evec, psVector* eval, const psImage* in)
     1121psImage *psMatrixSVD_old(psImage* evec, psVector* eval, const psImage* in)
    10991122{
    11001123    #define psMatrixSVD_EXIT {psFree(evec); psFree(eval); return NULL;}
    … …  
    11451168    return evec;
    11461169}
     1170
     1171// this is basically a wrapper for the gsl function: gsl_linalg_SV_decomp() SVD decomposes
     1172// matrix A based on the following equation: A = U w V^T .  This function (as usual for SVD
     1173// implementations) returns V not V^T.  U and V are returned to images; w is returned to a
     1174// vector representing the diagonal of w.  The input image A is not modified.  U, V, and w may
     1175// be supplied as NULL or may be allocated; their lengths are set here to match the
     1176// dimensionality of A.  XXX there is no error handling for the gsl functions (anywhere in
     1177// psMatrix.c)
     1178bool psMatrixSVD(psImage **U, psVector **w, psImage **V, const psImage *A)
     1179{
     1180    // Error checks  Missing one for eval
     1181    PS_ASSERT_PTR_NON_NULL(U, false);
     1182    PS_ASSERT_PTR_NON_NULL(w, false);
     1183    PS_ASSERT_PTR_NON_NULL(V, false);
     1184    PS_ASSERT_PTR_NON_NULL(A, false);
     1185
     1186    // A is provided with size Nx,Ny = numCols,numRows
     1187    // U has size Nx,Ny
     1188    // V has size Nx,Nx
     1189    // w has size Nx
     1190
     1191    // Initialize data
     1192    int numRows = A->numRows;
     1193    int numCols = A->numCols;
     1194
     1195    *U = psImageRecycle(*U,  numCols, numRows, A->type.type);
     1196    *V = psImageRecycle(*V,  numCols, numCols, A->type.type);
     1197    *w = psVectorRecycle(*w, numCols, A->type.type);
     1198
     1199    gsl_matrix *Agsl = gsl_matrix_alloc(numRows, numCols);
     1200    gsl_matrix *Vgsl = gsl_matrix_alloc(numCols, numCols);
     1201    gsl_vector *Sgsl = gsl_vector_alloc(numCols);
     1202    gsl_vector *work = gsl_vector_alloc(numCols);
     1203
     1204    // Copy psImage data into GSL matrix data
     1205    matrixPStoGSL(Agsl, A);
     1206
     1207    // Calculate SVD decomposition
     1208    gsl_linalg_SV_decomp(Agsl, Vgsl, Sgsl, work);
     1209
     1210    // Copy GSL matrix data to psImage data
     1211    matrixGSLtoPS(*V, Vgsl);
     1212    matrixGSLtoPS(*U, Agsl);  // gsl_linalg_SV_decomp replaces A with U
     1213    vectorGSLtoPS(*w, Sgsl);
     1214
     1215    // Free GSL data
     1216    gsl_matrix_free(Agsl);
     1217    gsl_matrix_free(Vgsl);
     1218    gsl_vector_free(Sgsl);
     1219    gsl_vector_free(work);
     1220
     1221    return true;
     1222}
     1223

  • trunk/psLib/src/math/psMatrix.h

    r26001 r26892  
    192192    psVector *solution,                 ///< Solution to output, or NULL
    193193    const psImage *matrix,              ///< Matrix to be solved
    194     const psVector *vector              ///< Vector of values
     194    const psVector *vector,             ///< Vector of values
     195    float thresh                        ///< Threshold relative to maximum for trimming singular values
    195196    );
    196197
    197 
    198 /// Single value decomposition, provided by Andy Becker
    199 psImage *psMatrixSVD(psImage* evec, psVector* eval, const psImage* in);
     198/// Single value decomposition (original by Andy Becker, updated by EAM)
     199bool psMatrixSVD(psImage **U, psVector **w, psImage **V, const psImage *A);
    200200
    201201/// @}

  • trunk/psLib/src/math/psUnaryOp.c

    r17050 r26892  
    211211}
    212212
    213 psMathType* psUnaryOp(psPtr out, const psPtr in, const char *op)
     213psMathType* psUnaryOp(psPtr out, psPtr in, const char *op)
    214214{
    215215    #define psUnaryOp_EXIT { \

  • trunk/psLib/src/math/psUnaryOp.h

    r11248 r26892  
    5959psMathType* psUnaryOp(
    6060    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.
    6262    const char *op                     ///< Operator.
    6363);

  • trunk/psLib/src/mathtypes/psVector.c

    r24886 r26892  
    729729    char line[1024];
    730730
    731     sprintf (line, "vector: %s\n", name);
     731    sprintf (line, "# vector: %s\n", name);
    732732    if (write(fd, line, strlen(line))) {;} //ignore return value
    733733

  • trunk/psLib/src/sys/psErrorCodes.c.in

    r11675 r26892  
    6767static void freeErrorDescription(psErrorDescription* err)
    6868{
    69     psFree((void *)err->description);
     69    psFree(err->description);
    7070}
    7171

  • trunk/psLib/src/sys/psMemory.h

    r23305 r26892  
    326326
    327327/** Free memory.  This operates much like free().
    328  *
     328 * 
    329329 *  @see psAlloc, psRealloc
     330 *  note: we cast ptr to (void *) in case we are supplied a const pointer.
    330331 */
    331332#ifdef DOXYGEN
    … …  
    336337#ifndef SWIG
    337338#define psFree(ptr) \
    338         psMemDecrRefCounter(ptr)
     339    ptr = psMemDecrRefCounter((void *)ptr);
    339340#endif // ifndef SWIG
    340341#endif // ifdef DOXYGEN

  • trunk/psLib/src/sys/psTrace.c

    r20546 r26892  
    119119
    120120    psMemSetPersistent((psPtr)comp->name,false);
    121     psFree((void *)comp->name);
     121    psFree(comp->name);
    122122}
    123123

  • trunk/psLib/src/types/psArray.c

    r15714 r26892  
    166166// drop an item from the array and free it
    167167bool psArrayRemoveData(psArray* array,
    168                        const psPtr data)
     168                       psPtr data)
    169169{
    170170    PS_ASSERT_ARRAY_NON_NULL(array, false);

  • trunk/psLib/src/types/psArray.h

    r19056 r26892  
    186186bool psArrayRemoveData(
    187187    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
    189189);
    190190

  • trunk/psLib/src/types/psMetadata.c

    r25383 r26892  
    622622
    623623// may need to extend this to change the keyname in the copy
    624 bool psMetadataItemSupplement(psMetadata *out,
     624bool psMetadataItemSupplement(bool *status,
     625                              psMetadata *out,
    625626                              const psMetadata *in,
    626627                              const char *key)
    … …  
    632633    psMetadataItem *item = psMetadataLookup(in, key);
    633634    if (!item) {
    634         psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Could not find '%s' in metadata.\n", key);
     635        if (status) {
     636            *status = false;
     637        } else {
     638            psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Could not find '%s' in metadata.\n", key);
     639        }
    635640        return false;
    636641    }
    … …  
    968973}
    969974
    970 psMetadataItem* psMetadataLookup(const psMetadata *md,
     975// Get entry by index
     976psMetadataItem *psMetadataGetIndex(psMetadata *md, long location)
     977{
     978    PS_ASSERT_METADATA_NON_NULL(md, false);
     979    return psListGet(md->list, location);
     980}
     981
     982psMetadataItem *psMetadataLookup(const psMetadata *md,
    971983                                 const char *key)
    972984{

  • trunk/psLib/src/types/psMetadata.h

    r25383 r26892  
    545545 */
    546546bool psMetadataItemSupplement(
     547    bool *status,                       ///< if supplied, returns true/false if key is found (suppresses the error)
    547548    psMetadata *out,                   ///< output Metadata container for copying.
    548549    const psMetadata *in,              ///< Metadata collection from which to copy.
Note: See TracChangeset for help on using the changeset viewer.