IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 26202


Ignore:
Timestamp:
Nov 19, 2009, 3:42:18 PM (16 years ago)
Author:
Paul Price
Message:

Threading for pmReadoutFake is working and active in psphot efficiency measurement and ppStack PSF matching.

Location:
trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack/src/ppStackMatch.c

    r26080 r26202  
    344344                                                       footprint); // Filtered list of sources
    345345
     346            bool oldThreads = pmReadoutFakeThreads(true); // Old threading state
    346347            if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows,
    347348                                          stampSources, SOURCE_MASK, NULL, NULL, options->psf,
    … …  
    353354                return false;
    354355            }
     356            pmReadoutFakeThreads(oldThreads);
    355357
    356358            fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK);

  • trunk/psModules/src/camera/pmReadoutFake.c

    r26076 r26202  
    2323#include "pmSourceUtils.h"
    2424#include "pmModelUtils.h"
     25#include "pmSourceGroups.h"
    2526
    2627#include "pmReadoutFake.h"
    2728
    28 #define MODEL_TYPE "PS_MODEL_RGAUSS"    // Type of model to use
    2929#define MAX_AXIS_RATIO 20.0             // Maximum axis ratio for PSF model
    3030#define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \
    3131                    PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models
     32
     33
     34static bool threaded = false;           // Running threaded?
    3235
    3336
    … …  
    4750    }
    4851    return pmPSF_AxesToModel(params, axes);
     52}
     53
     54/// Generate fake sources on a readout
     55static bool readoutFake(pmReadout *readout, // Readout of interest
     56                        const pmSourceGroups *groups, // Source groups
     57                        const psVector *x,        // x coordinates
     58                        const psVector *y,        // y coordinates
     59                        const psVector *mag,      // Magnitudes
     60                        const psVector *xOffset,  // Offsets in x
     61                        const psVector *yOffset,  // Offsets in y
     62                        const pmPSF *psf,         // PSF
     63                        float minFlux,            // Minimum flux
     64                        float radius,             // Minimum radius
     65                        bool circularise,         // Circularise PSF?
     66                        bool normalisePeak,       // Normalise sources for peak?
     67                        int groupIndex,           // Group index
     68                        int cellIndex             // Cell index
     69                        )
     70{
     71    psArray *cells = groups->groups->data[groupIndex]; // Cells in group
     72    psVector *cellSources = cells->data[cellIndex];    // Sources in cell
     73
     74    for (int i = 0; i < cellSources->n; i++) {
     75        int index = cellSources->data.S32[i];                       // Index for source of interest
     76        float flux = powf(10.0, -0.4 * mag->data.F32[index]);       // Flux of source
     77        float xSrc = x->data.F32[index], ySrc = y->data.F32[index]; // Coordinates of source
     78
     79        if (normalisePeak) {
     80            // Normalise flux
     81            pmModel *normModel = pmModelFromPSFforXY(psf, xSrc, ySrc, 1.0); // Model for normalisation
     82            if (!normModel || (normModel->flags & MODEL_MASK)) {
     83                psFree(normModel);
     84                continue;
     85            }
     86            // check that all params are valid:
     87            bool validParams = true;
     88            for (int j = 0; validParams && (j < normModel->params->n); j++) {
     89                switch (j) {
     90                  case PM_PAR_SKY:
     91                  case PM_PAR_I0:
     92                  case PM_PAR_XPOS:
     93                  case PM_PAR_YPOS:
     94                    continue;
     95                  default:
     96                    if (!isfinite(normModel->params->data.F32[j])) {
     97                        validParams = false;
     98                    }
     99                }
     100            }
     101            if (!validParams) {
     102                psFree(normModel);
     103                continue;
     104            }
     105            if (circularise && !circulariseModel(normModel)) {
     106                psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
     107                psFree(normModel);
     108                return false;
     109            }
     110
     111            flux /= normModel->modelFlux(normModel->params);
     112            psFree(normModel);
     113        }
     114
     115        pmModel *fakeModel = pmModelFromPSFforXY(psf, xSrc, ySrc, flux);
     116        if (!fakeModel || (fakeModel->flags & MODEL_MASK)) {
     117            psFree(fakeModel);
     118            continue;
     119        }
     120        // check that all params are valid:
     121        bool validParams = true;
     122        for (int j = 0; validParams && (j < fakeModel->params->n); j++) {
     123            switch (j) {
     124              case PM_PAR_SKY:
     125              case PM_PAR_I0:
     126              case PM_PAR_XPOS:
     127              case PM_PAR_YPOS:
     128                continue;
     129              default:
     130                if (!isfinite(fakeModel->params->data.F32[j])) {
     131                    validParams = false;
     132                }
     133            }
     134        }
     135        if (!validParams) {
     136            psFree(fakeModel);
     137            continue;
     138        }
     139        if (circularise && !circulariseModel(fakeModel)) {
     140            psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
     141            psFree(fakeModel);
     142            return false;
     143        }
     144
     145        psTrace("psModules.camera", 10, "Adding source at %f,%f with flux %f\n",
     146                fakeModel->params->data.F32[PM_PAR_XPOS], fakeModel->params->data.F32[PM_PAR_YPOS],
     147                fakeModel->params->data.F32[PM_PAR_I0]);
     148
     149        pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate
     150        fakeSource->peak = pmPeakAlloc(xSrc, ySrc, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE);
     151        float fakeRadius = 1.0;         // Radius of fake source
     152        if (isfinite(minFlux)) {
     153            fakeRadius = PS_MAX(fakeRadius, fakeModel->modelRadius(fakeModel->params, minFlux));
     154        }
     155        if (radius > 0) {
     156            fakeRadius = PS_MAX(fakeRadius, radius);
     157        }
     158
     159        if (xOffset) {
     160            if (!pmSourceDefinePixels(fakeSource, readout, xSrc + xOffset->data.S32[index],
     161                                      ySrc + yOffset->data.S32[index], fakeRadius)) {
     162                psErrorClear();
     163                continue;
     164            }
     165            if (!pmModelAddWithOffset(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0,
     166                                      xOffset->data.S32[index], yOffset->data.S32[index])) {
     167                psErrorClear();
     168                continue;
     169            }
     170        } else {
     171            if (!pmSourceDefinePixels(fakeSource, readout, xSrc, ySrc, fakeRadius)) {
     172                psErrorClear();
     173                continue;
     174            }
     175            if (!pmModelAdd(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0)) {
     176                psErrorClear();
     177                continue;
     178            }
     179        }
     180        psFree(fakeSource);
     181        psFree(fakeModel);
     182    }
     183
     184    return true;
     185}
     186
     187/// Thread job for readoutFake()
     188static bool readoutFakeThread(psThreadJob *job)
     189{
     190    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     191
     192    psArray *args = job->args;          // Arguments
     193
     194    pmReadout *readout = args->data[0];     // Readout of interest
     195    const pmSourceGroups *groups = args->data[1]; // Source groups
     196    const psVector *x = args->data[2];        // x coordinates
     197    const psVector *y = args->data[3];        // y coordinates
     198    const psVector *mag = args->data[4];      // Magnitudes
     199    const psVector *xOffset = args->data[5];  // Offsets in x
     200    const psVector *yOffset = args->data[6];  // Offsets in y
     201    const pmPSF *psf = args->data[7];         // PSF
     202    float minFlux = PS_SCALAR_VALUE(args->data[8], F32); // Minimum flux
     203    float radius = PS_SCALAR_VALUE(args->data[9], F32);  // Minimum radius
     204    bool circularise = PS_SCALAR_VALUE(args->data[10], U8); // Circularise PSF?
     205    bool normalisePeak = PS_SCALAR_VALUE(args->data[11], U8); // Normalise for peak?
     206    int groupIndex = PS_SCALAR_VALUE(args->data[12], S32); // Group index
     207    int cellIndex = PS_SCALAR_VALUE(args->data[13], S32);  // Cell index
     208
     209    return readoutFake(readout, groups, x, y, mag, xOffset, yOffset, psf, minFlux, radius, circularise,
     210                       normalisePeak, groupIndex, cellIndex);
     211}
     212
     213
     214bool pmReadoutFakeThreads(bool new)
     215{
     216    bool old = threaded;                // Old status, to return
     217
     218    if (!old && new) {
     219        threaded = true;
     220
     221        {
     222            psThreadTask *task = psThreadTaskAlloc("PSMODULES_READOUT_FAKE", 14);
     223            task->function = &readoutFakeThread;
     224            psThreadTaskAdd(task);
     225            psFree(task);
     226        }
     227
     228    } else if (old && !new) {
     229        threaded = false;
     230        psThreadTaskRemove("PSMODULES_READOUT_FAKE");
     231    }
     232
     233    return old;
    49234}
    50235
    … …  
    86271    psImageInit(readout->image, 0);
    87272
    88     for (long i = 0; i < numSources; i++) {
    89         float flux = powf(10.0, -0.4 * mag->data.F32[i]); // Flux of source
    90         float xSrc = x->data.F32[i], ySrc = y->data.F32[i]; // Coordinates of source
    91 
    92         if (normalisePeak) {
    93             // Normalise flux
    94             pmModel *normModel = pmModelFromPSFforXY(psf, xSrc, ySrc, 1.0); // Model for normalisation
    95             if (!normModel || (normModel->flags & MODEL_MASK)) {
    96                 psFree(normModel);
    97                 continue;
    98             }
    99             // check that all params are valid:
    100             bool validParams = true;
    101             for (int n = 0; validParams && (n < normModel->params->n); n++) {
    102                 if (n == PM_PAR_SKY) continue;
    103                 if (n == PM_PAR_I0) continue;
    104                 if (n == PM_PAR_XPOS) continue;
    105                 if (n == PM_PAR_YPOS) continue;
    106                 if (!isfinite(normModel->params->data.F32[n])) validParams = false;
    107             }
    108             if (!validParams) {
    109                 psFree(normModel);
    110                 continue;
    111             }           
    112             if (circularise && !circulariseModel(normModel)) {
    113                 psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
    114                 psFree(normModel);
     273    int numThreads = threaded ? psThreadPoolSize() : 0; // Number of threads
     274    pmSourceGroups *groups = pmSourceGroupsFromVectors(readout, x, y, numThreads); // Groups of sources
     275    if (!groups) {
     276        psError(PS_ERR_UNKNOWN, false, "Unable to generate source groups");
     277        return false;
     278    }
     279
     280    if (threaded) {
     281        for (int i = 0; i < groups->groups->n; i++) {
     282            psArray *cells = groups->groups->data[i]; // Cell with sources
     283            for (int j = 0; j < cells->n; j++) {
     284                psThreadJob *job = psThreadJobAlloc("PSMODULES_READOUT_FAKE");
     285                psArray *args = job->args;
     286                psArrayAdd(args, 1, readout);
     287                psArrayAdd(args, 1, groups);
     288                // Casting away const to add to array
     289                psArrayAdd(args, 1, (psVector*)x);
     290                psArrayAdd(args, 1, (psVector*)y);
     291                psArrayAdd(args, 1, (psVector*)mag);
     292                psArrayAdd(args, 1, (psVector*)xOffset);
     293                psArrayAdd(args, 1, (psVector*)yOffset);
     294                psArrayAdd(args, 1, (pmPSF*)psf);
     295                PS_ARRAY_ADD_SCALAR(args, minFlux, PS_TYPE_F32);
     296                PS_ARRAY_ADD_SCALAR(args, radius, PS_TYPE_S32);
     297                PS_ARRAY_ADD_SCALAR(args, circularise, PS_TYPE_U8);
     298                PS_ARRAY_ADD_SCALAR(args, normalisePeak, PS_TYPE_U8);
     299                PS_ARRAY_ADD_SCALAR(args, i, PS_TYPE_S32);
     300                PS_ARRAY_ADD_SCALAR(args, j, PS_TYPE_S32);
     301
     302                if (!psThreadJobAddPending(job)) {
     303                    psFree(job);
     304                    psFree(groups);
     305                    return false;
     306                }
     307                psFree(job);
     308            }
     309            if (!psThreadPoolWait(true)) {
     310                psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
     311                psFree(groups);
    115312                return false;
    116313            }
    117 
    118             flux /= normModel->modelFlux(normModel->params);
    119             psFree(normModel);
    120         }
    121 
    122         pmModel *fakeModel = pmModelFromPSFforXY(psf, xSrc, ySrc, flux);
    123         if (!fakeModel || (fakeModel->flags & MODEL_MASK)) {
    124             psFree(fakeModel);
    125             continue;
    126         }
    127         // check that all params are valid:
    128         bool validParams = true;
    129         for (int n = 0; validParams && (n < fakeModel->params->n); n++) {
    130             if (n == PM_PAR_SKY) continue;
    131             if (n == PM_PAR_I0) continue;
    132             if (n == PM_PAR_XPOS) continue;
    133             if (n == PM_PAR_YPOS) continue;
    134             if (!isfinite(fakeModel->params->data.F32[n])) validParams = false;
    135         }
    136         if (!validParams) {
    137             psFree(fakeModel);
    138             continue;
    139         }               
    140         if (circularise && !circulariseModel(fakeModel)) {
    141             psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
    142             psFree(fakeModel);
    143             return false;
    144         }
    145 
    146         psTrace("psModules.camera", 10, "Adding source at %f,%f with flux %f\n",
    147                 fakeModel->params->data.F32[PM_PAR_XPOS], fakeModel->params->data.F32[PM_PAR_YPOS],
    148                 fakeModel->params->data.F32[PM_PAR_I0]);
    149 
    150         pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate
    151         fakeSource->peak = pmPeakAlloc(xSrc, ySrc, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE);
    152         float fakeRadius = 1.0;         // Radius of fake source
    153         if (isfinite(minFlux)) {
    154             fakeRadius = PS_MAX(fakeRadius, fakeModel->modelRadius(fakeModel->params, minFlux));
    155         }
    156         if (radius > 0) {
    157             fakeRadius = PS_MAX(fakeRadius, radius);
    158         }
    159 
    160         if (xOffset) {
    161             if (!pmSourceDefinePixels(fakeSource, readout, xSrc + xOffset->data.S32[i],
    162                                       ySrc + yOffset->data.S32[i], fakeRadius)) {
    163                 psErrorClear();
    164                 continue;
    165             }
    166             if (!pmModelAddWithOffset(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0,
    167                                       xOffset->data.S32[i], yOffset->data.S32[i])) {
    168                 psErrorClear();
    169                 continue;
    170             }
    171         } else {
    172             if (!pmSourceDefinePixels(fakeSource, readout, xSrc, ySrc, fakeRadius)) {
    173                 psErrorClear();
    174                 continue;
    175             }
    176             if (!pmModelAdd(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0)) {
    177                 psErrorClear();
    178                 continue;
    179             }
    180         }
    181         psFree(fakeSource);
    182         psFree(fakeModel);
    183     }
     314        }
     315    } else if (!readoutFake(readout, groups, x, y, mag, xOffset, yOffset, psf, minFlux, radius, circularise,
     316                            normalisePeak, 0, 0)) {
     317        psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources on readout");
     318        psFree(groups);
     319        return false;
     320    }
     321
     322    psFree(groups);
    184323
    185324    return true;
     325
    186326}
    187327

  • trunk/psModules/src/camera/pmReadoutFake.h

    r25383 r26202  
    1212#include <pmPSF.h>
    1313#include <pmSourceMasks.h>
     14
     15/// Set threading
     16///
     17/// Returns old threading state
     18bool pmReadoutFakeThreads(
     19    bool new                            // New threading state
     20    );
    1421
    1522/// Generate a fake readout from vectors

  • trunk/psphot/src/psphotEfficiency.c

    r25890 r26202  
    128128    psFree(rng);
    129129
     130    bool oldThreads = pmReadoutFakeThreads(true); // Old threading status
     131
    130132    pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout
    131133    if (!pmReadoutFakeFromVectors(fakeRO, numCols, numRows, xAll, yAll, magAll,
    … …  
    141143    psFree(xAll);
    142144    psFree(yAll);
     145
     146    pmReadoutFakeThreads(oldThreads);
    143147
    144148    psBinaryOp(ro->image, ro->image, "+", fakeRO->image);
Note: See TracChangeset for help on using the changeset viewer.