IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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

updates from eam_branches/20091201

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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}
Note: See TracChangeset for help on using the changeset viewer.