IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26794


Ignore:
Timestamp:
Feb 5, 2010, 3:14:03 PM (16 years ago)
Author:
Paul Price
Message:

Threading psImageCovarianceCalculate and psImageCovarianceBin.

Location:
branches/eam_branches/20091201
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20091201/ppSmooth/src/ppSmoothReadout.c

    r26008 r26794  
    5151    psImageSmoothMask_Threaded(input->variance, input->variance, input->mask, maskVal, sigma * M_SQRT1_2, nSigma, minGauss);
    5252    psLogMsg("ppSmooth", PS_LOG_MINUTIA, "smooth variance: %f sec\n", psTimerMark("ppSmooth"));
     53    psImageConvolveSetThreads(oldThreads);
    5354
    5455    // determine covariance matrix for this smoothing, replace existing kernel
     56    oldThreads = psImageCovarianceSetThreads(true);
    5557    psKernel *kernel = psImageSmoothKernel(sigma, nSigma); // Kernel used for smoothing
    5658    psKernel *covar = psImageCovarianceCalculate(kernel, input->covariance); // Covariance matrix
    … …  
    5860    input->covariance = covar;
    5961    psFree(kernel);
    60    
    61     // Calculate correction factor for the covariance produced by the (potentially multiple) smoothing
    6262    float factor = 1.0 / psImageCovarianceFactor(covar);
     63    psImageCovarianceSetThreads(oldThreads);
    6364
    6465    // record the effective area and significance scaling factor
    … …  
    6667    psMetadataAddF32(recipe, PS_LIST_TAIL, "EFFECTIVE_AREA", PS_META_REPLACE, "Effective Area", effArea);
    6768    psMetadataAddF32(recipe, PS_LIST_TAIL, "SIGNIFICANCE_SCALE_FACTOR", PS_META_REPLACE, "Signicance scale factor", factor);
    68 
    69     psImageConvolveSetThreads(oldThreads);
    7069
    7170    return true;

  • branches/eam_branches/20091201/ppStack/src/ppStackMatch.c

    r26769 r26794  
    278278
    279279            psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
     280            bool oldThreads = psImageCovarianceSetThreads(true);              // Old thread setting
    280281            psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix
     282            psImageCovarianceSetThreads(oldThreads);
    281283            psFree(readout->covariance);
    282284            readout->covariance = covar;

  • branches/eam_branches/20091201/psLib/src/imageops/psImageCovariance.c

    r26704 r26794  
    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 float psImageCovarianceCalculateFactor(const psKernel *kernel, const psKernel *covariance)
    27 {
    28     psKernel *covar;                    // Covariance matrix to use
    29     if (covariance) {
    30         covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const
    31     } else {
    32         covar = psImageCovarianceNone();
    33     }
    34 
    35     // Check for non-finite elements
    36     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
    37         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
    38             if (!isfinite(kernel->kernel[y][x])) {
    39                 psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    40                         "Non-finite covariance matrix element in kernel at %d,%d", x, y);
    41                 psFree(covar);
    42                 return NAN;
    43             }
    44         }
    45     }
    46     for (int y = covar->yMin; y <= covar->yMax; y++) {
    47         for (int x = covar->xMin; x <= covar->xMax; x++) {
    48             if (!isfinite(covar->kernel[y][x])) {
    49                 psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    50                         "Non-finite covariance matrix element in covariance matrix at %d,%d", x, y);
    51                 psFree(covar);
    52                 return NAN;
    53             }
    54         }
    55     }
    56 
    57     int x = 0, y = 0;                   // Coordinates on the output covariance matrix
    58 
    59     // Range for v
    60     int vMin = PS_MAX(kernel->yMin + covar->yMin, y + kernel->yMin);
    61     int vMax = PS_MIN(kernel->yMax + covar->yMax, y + kernel->yMax);
    62     // Range for u
    63     int uMin = PS_MAX(kernel->xMin + covar->xMin, x + kernel->xMin);
    64     int uMax = PS_MIN(kernel->xMax + covar->xMax, x + kernel->xMax);
    65 
    66     double sum = 0.0;           // Sum for value of covariance matrix at (x,y)
    67     for (int v = vMin; v <= vMax; v++) {
    68         // Range for q
    69         int qMin = PS_MAX(v + covar->yMin, kernel->yMin);
    70         int qMax = PS_MIN(v + covar->yMax, kernel->yMax);
    71         for (int u = uMin; u <= uMax; u++) {
    72             // Range for p
    73             int pMin = PS_MAX(u + covar->xMin, kernel->xMin);
    74             int pMax = PS_MIN(u + covar->xMax, kernel->xMax);
    75 
    76             double xyuvValue = kernel->kernel[v-y][u-x]; // Value for (x,y) --> (u,v)
    77 
    78             double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0)
    79             for (int q = qMin; q <= qMax; q++) {
    80                 for (int p = pMin; p <= pMax; p++) {
    81                     uvpqValue += (double)covar->kernel[q-v][p-u] * (double)kernel->kernel[q][p];
    82                 }
    83             }
    84             sum += xyuvValue * uvpqValue;
    85         }
    86     }
    87 
    88     psFree(covar);
    89     return sum;
    90 }
    91 
    92 psKernel *psImageCovarianceCalculate(const psKernel *kernel, const psKernel *covariance)
    93 {
    94     PS_ASSERT_KERNEL_NON_NULL(kernel, NULL);
    95 
    96     // See http://en.wikipedia.org/wiki/Error_propagation
    97     //
    98     // If
    99     //     f_k = sum_i A_ik x_i
    100     // is a set of functions, then the covariance matrix for f is given by:
    101     //     M^f_ij = sum_k sum_l A_ik M^x_kl A_lj
    102     // where M^x is the covariance matrix for x.
    103     // Note that the errors in f are correlated (covariance) even if the errors in x are not.
    104 
    105     psKernel *covar;                    // Covariance matrix to use
    106     if (covariance) {
    107         covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const
    108     } else {
    109         covar = psImageCovarianceNone();
    110     }
    111 
    112     // Check for non-finite elements
    113     for (int y = kernel->yMin; y <= kernel->yMax; y++) {
    114         for (int x = kernel->xMin; x <= kernel->xMax; x++) {
    115             if (!isfinite(kernel->kernel[y][x])) {
    116                 psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    117                         "Non-finite covariance matrix element in kernel at %d,%d", x, y);
    118                 psFree(covar);
    119                 return NULL;
    120             }
    121         }
    122     }
    123     for (int y = covar->yMin; y <= covar->yMax; y++) {
    124         for (int x = covar->xMin; x <= covar->xMax; x++) {
    125             if (!isfinite(covar->kernel[y][x])) {
    126                 psError(PS_ERR_BAD_PARAMETER_VALUE, true,
    127                         "Non-finite covariance matrix element in covariance matrix at %d,%d", x, y);
    128                 psFree(covar);
    129                 return NULL;
    130             }
    131         }
    132     }
    133 
    134     // The above (double) sum for the covariance matrix means that, for each point in the output covariance
    135     // matrix, we need to work out all combinations of getting to the central point via a kernel, input
    136     // covariance matrix and another kernel.  This means that the resultant covariance matrix has twice the
    137     // size of the kernel plus the size of the input covariance matrix.
    138     int xMin = kernel->xMin - kernel->xMax + covar->xMin, xMax = kernel->xMax - kernel->xMin + covar->xMax;
    139     int yMin = kernel->yMin - kernel->yMax + covar->yMin, yMax = kernel->yMax - kernel->yMin + covar->yMax;
    140     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");
    14139
    14240    // Need to go:
    … …  
    15452    // from the source coordinate), we take the smallest possible (because everything else is zero outside).
    15553
    156     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
    157156    for (int y = yMin; y <= yMax; y++) {
    158         // Range for v
    159         int vMin = PS_MAX(kernel->yMin + covar->yMin, y + kernel->yMin);
    160         int vMax = PS_MIN(kernel->yMax + covar->yMax, y + kernel->yMax);
    161157        for (int x = xMin; x <= xMax; x++) {
    162             // Range for u
    163             int uMin = PS_MAX(kernel->xMin + covar->xMin, x + kernel->xMin);
    164             int uMax = PS_MIN(kernel->xMax + covar->xMax, x + kernel->xMax);
    165 
    166             double sum = 0.0;           // Sum for value of covariance matrix at (x,y)
    167             for (int v = vMin; v <= vMax; v++) {
    168                 // Range for q
    169                 int qMin = PS_MAX(v + covar->yMin, kernel->yMin);
    170                 int qMax = PS_MIN(v + covar->yMax, kernel->yMax);
    171                 for (int u = uMin; u <= uMax; u++) {
    172                     // Range for p
    173                     int pMin = PS_MAX(u + covar->xMin, kernel->xMin);
    174                     int pMax = PS_MIN(u + covar->xMax, kernel->xMax);
    175 
    176                     double xyuvValue = kernel->kernel[v-y][u-x]; // Value for (x,y) --> (u,v)
    177 
    178                     double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0)
    179                     for (int q = qMin; q <= qMax; q++) {
    180                         for (int p = pMin; p <= pMax; p++) {
    181                             uvpqValue += (double)covar->kernel[q-v][p-u] * (double)kernel->kernel[q][p];
    182                         }
    183                     }
    184                     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;
    185168                }
    186             }
    187             out->kernel[y][x] = sum;
    188             total += sum;
    189         }
    190     }
    191     psTrace("psLib.imageops", 3, "Total covariance: %lf ; Central variance: %f\n", total, out->kernel[0][0]);
    192 
     169                psFree(job);
     170            } else {
     171                out->kernel[y][x] = imageCovarianceCalculate(covar, kernel, x, y);
     172            }
     173        }
     174    }
    193175    psFree(covar);
     176
    194177    return out;
    195178}
     179
     180float psImageCovarianceCalculateFactor(const psKernel *kernel, const psKernel *covariance)
     181{
     182    psKernel *covar;                    // Covariance matrix to use
     183    if (covariance) {
     184        covar = psMemIncrRefCounter((psKernel*)covariance); // Casting away const
     185    } else {
     186        covar = psImageCovarianceNone();
     187    }
     188
     189    // Check for non-finite elements
     190    for (int y = kernel->yMin; y <= kernel->yMax; y++) {
     191        for (int x = kernel->xMin; x <= kernel->xMax; x++) {
     192            if (!isfinite(kernel->kernel[y][x])) {
     193                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     194                        "Non-finite covariance matrix element in kernel at %d,%d", x, y);
     195                psFree(covar);
     196                return NAN;
     197            }
     198        }
     199    }
     200    for (int y = covar->yMin; y <= covar->yMax; y++) {
     201        for (int x = covar->xMin; x <= covar->xMax; x++) {
     202            if (!isfinite(covar->kernel[y][x])) {
     203                psError(PS_ERR_BAD_PARAMETER_VALUE, true,
     204                        "Non-finite covariance matrix element in covariance matrix at %d,%d", x, y);
     205                psFree(covar);
     206                return NAN;
     207            }
     208        }
     209    }
     210
     211    float factor = imageCovarianceCalculate(covar, kernel, 0, 0); // Covariance factor
     212    psFree(covar);
     213    return factor;
     214}
     215
     216// Calculation of covariance matrix element when binning
     217static float imageCovarianceBin(const psKernel *covar, // Original covariance matrix
     218                                int bin,               // Binning factor
     219                                float binVal,          // Convolution kernel value for binning
     220                                int x, int y           // Coordinates in output covariance matrix
     221    )
     222{
     223    psAssert(covar, "Require covariance matrix");
     224    psAssert(bin > 0 && binVal > 0, "Require binning: %d %f", bin, binVal);
     225
     226    int binMin = -(bin - 1) / 2, binMax = bin / 2; // Range of "kernel"
     227
     228    // Range for v
     229    int vMin = PS_MAX(binMin + covar->yMin, bin * y + binMin);
     230    int vMax = PS_MIN(binMax + covar->yMax, bin * y + binMax);
     231    // Range for u
     232    int uMin = PS_MAX(binMin + covar->xMin, bin * x + binMin);
     233    int uMax = PS_MIN(binMax + covar->xMax, bin * x + binMax);
     234
     235    double sum = 0.0;           // Sum for value of covariance matrix at (x,y)
     236    for (int v = vMin; v <= vMax; v++) {
     237        // Range for q
     238        int qMin = PS_MAX(v + covar->yMin, binMin);
     239        int qMax = PS_MIN(v + covar->yMax, binMax);
     240        for (int u = uMin; u <= uMax; u++) {
     241            // Range for p
     242            int pMin = PS_MAX(u + covar->xMin, binMin);
     243            int pMax = PS_MIN(u + covar->xMax, binMax);
     244
     245            double xyuvValue = binVal; // Value for (x,y) --> (u,v)
     246
     247            double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0)
     248            for (int q = qMin; q <= qMax; q++) {
     249                for (int p = pMin; p <= pMax; p++) {
     250                    uvpqValue += (double)covar->kernel[q-v][p-u] * (double)binVal;
     251                }
     252            }
     253            sum += xyuvValue * uvpqValue;
     254        }
     255    }
     256
     257    return sum;
     258}
     259
     260/// Thread entry point for calculation of covariance matrix element when binning
     261static bool imageCovarianceBinThread(psThreadJob *job)
     262{
     263    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     264    psAssert(job->args, "No job arguments");
     265    psAssert(job->args->n == 6, "Wrong number of job arguments: %ld", job->args->n);
     266
     267    psKernel *out = job->args->data[0]; // Output covariance matrix
     268    const psKernel *covar = job->args->data[1]; // Input covariance matrix
     269    int bin = PS_SCALAR_VALUE(job->args->data[2], S32); // Binning factor
     270    float binVal = PS_SCALAR_VALUE(job->args->data[3], F32); // Convolution kernel value for binning
     271    int x = PS_SCALAR_VALUE(job->args->data[4], S32); // x coordinate in output covariance matrix
     272    int y = PS_SCALAR_VALUE(job->args->data[5], S32); // y coordinate in output covariance matrix
     273
     274    out->kernel[y][x] = imageCovarianceBin(covar, bin, binVal, x, y);
     275
     276    return true;
     277}
     278
    196279
    197280psKernel *psImageCovarianceBin(int bin, const psKernel *covariance, bool average)
    … …  
    230313    psKernel *out = psKernelAlloc(xMin, xMax, yMin, yMax); // Covariance matrix for output
    231314
    232     double total = 0.0;                 // Total covariance
    233315    for (int y = yMin; y <= yMax; y++) {
    234         // Range for v
    235         int vMin = PS_MAX(binMin + covar->yMin, bin * y + binMin);
    236         int vMax = PS_MIN(binMax + covar->yMax, bin * y + binMax);
    237316        for (int x = xMin; x <= xMax; x++) {
    238             // Range for u
    239             int uMin = PS_MAX(binMin + covar->xMin, bin * x + binMin);
    240             int uMax = PS_MIN(binMax + covar->xMax, bin * x + binMax);
    241 
    242             double sum = 0.0;           // Sum for value of covariance matrix at (x,y)
    243             for (int v = vMin; v <= vMax; v++) {
    244                 // Range for q
    245                 int qMin = PS_MAX(v + covar->yMin, binMin);
    246                 int qMax = PS_MIN(v + covar->yMax, binMax);
    247                 for (int u = uMin; u <= uMax; u++) {
    248                     // Range for p
    249                     int pMin = PS_MAX(u + covar->xMin, binMin);
    250                     int pMax = PS_MIN(u + covar->xMax, binMax);
    251 
    252                     double xyuvValue = binVal; // Value for (x,y) --> (u,v)
    253 
    254                     double uvpqValue = 0.0; // Value for (u,v) --> (p,q) --> (0,0)
    255                     for (int q = qMin; q <= qMax; q++) {
    256                         for (int p = pMin; p <= pMax; p++) {
    257                             uvpqValue += (double)covar->kernel[q-v][p-u] * (double)binVal;
    258                         }
    259                     }
    260                     sum += xyuvValue * uvpqValue;
     317            if (threaded) {
     318                psThreadJob *job = psThreadJobAlloc("PSLIB_IMAGE_COVARIANCE_BIN");
     319                psArrayAdd(job->args, 1, out);
     320                psArrayAdd(job->args, 1, covar);
     321                PS_ARRAY_ADD_SCALAR(job->args, bin, PS_TYPE_S32);
     322                PS_ARRAY_ADD_SCALAR(job->args, binVal, PS_TYPE_F32);
     323                PS_ARRAY_ADD_SCALAR(job->args, x, PS_TYPE_S32);
     324                PS_ARRAY_ADD_SCALAR(job->args, y, PS_TYPE_S32);
     325                if (!psThreadJobAddPending(job)) {
     326                    psFree(covar);
     327                    return NULL;
    261328                }
    262             }
    263             out->kernel[y][x] = sum;
    264             total += sum;
    265         }
    266     }
    267     psTrace("psLib.imageops", 3, "Total covariance: %lf ; Central variance: %f\n", total, out->kernel[0][0]);
    268 
     329                psFree(job);
     330            } else {
     331                out->kernel[y][x] = imageCovarianceBin(covar, bin, binVal, x, y);
     332            }
     333        }
     334    }
    269335    psFree(covar);
    270336
    … …  
    493559    return true;
    494560}
     561
     562
     563bool psImageCovarianceSetThreads(bool set)
     564{
     565    bool old = threaded;                // Old value
     566    if (set && !threaded) {
     567        {
     568            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_CALCULATE", 5);
     569            task->function = &imageCovarianceCalculateThread;
     570            psThreadTaskAdd(task);
     571            psFree(task);
     572        }
     573        {
     574            psThreadTask *task = psThreadTaskAlloc("PSLIB_IMAGE_COVARIANCE_BIN", 6);
     575            task->function = &imageCovarianceBinThread;
     576            psThreadTaskAdd(task);
     577            psFree(task);
     578        }
     579    } else if (!set && threaded) {
     580        psThreadTaskRemove("PSLIB_IMAGE_COVARIANCE_CALCULATE");
     581        psThreadTaskRemove("PSLIB_IMAGE_COVARIANCE_BIN");
     582    }
     583    threaded = set;
     584    return old;
     585}
     586
     587bool psImageCovarianceGetThreads(void)
     588{
     589    return threaded;
     590}

  • branches/eam_branches/20091201/psLib/src/imageops/psImageCovariance.h

    r26704 r26794  
    9090    );
    9191
     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);
    92100
    93101/// @}

  • branches/eam_branches/20091201/psModules/src/imcombine/pmSubtraction.c

    r26769 r26794  
    13251325
    13261326    // Calculate covariances
    1327     // This can take a while, so we only do it for a single instance
    1328     // XXX psImageCovarianceCalculate could be multithreaded
     1327    // This can be fairly involved, so we only do it for a single instance
     1328    // Enable threads for covariance calculation, since we're not threading on top of it.
     1329    oldThreads = psImageCovarianceSetThreads(true);
    13291330    if (kernels->mode == PM_SUBTRACTION_MODE_1 || kernels->mode == PM_SUBTRACTION_MODE_DUAL) {
    13301331        psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel
    … …  
    13401341        psFree(kernel);
    13411342    }
     1343    psImageCovarianceSetThreads(oldThreads);
    13421344
    13431345    // Copy anything that wasn't convolved

  • branches/eam_branches/20091201/pswarp/src/pswarpTransformReadout.c

    r25043 r26794  
    107107        psImageInit(output->mask, maskBad);
    108108    }
     109
     110    // Ensure threading is off for the covariance calculation, since we are threading on a different level.
     111    psImageCovarianceSetThreads(false);
    109112
    110113    // create jobs and supply them to the threads
Note: See TracChangeset for help on using the changeset viewer.