IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 10, 2010, 4:24:46 PM (16 years ago)
Author:
eugene
Message:

updates from trunk

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20091201/psModules/src/camera/pmReadoutFake.c

    r26260 r26879  
    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?
     35
     36
    3237
    3338
    … …  
    4752    }
    4853    return pmPSF_AxesToModel(params, axes);
     54}
     55
     56/// Generate fake sources on a readout
     57static bool readoutFake(pmReadout *readout, // Readout of interest
     58                        const pmSourceGroups *groups, // Source groups
     59                        const psVector *x,        // x coordinates
     60                        const psVector *y,        // y coordinates
     61                        const psVector *mag,      // Magnitudes
     62                        const psVector *xOffset,  // Offsets in x
     63                        const psVector *yOffset,  // Offsets in y
     64                        const pmPSF *psf,         // PSF
     65                        float minFlux,            // Minimum flux
     66                        float radius,             // Minimum radius
     67                        bool circularise,         // Circularise PSF?
     68                        bool normalisePeak,       // Normalise sources for peak?
     69                        int groupIndex,           // Group index
     70                        int cellIndex             // Cell index
     71                        )
     72{
     73    psArray *cells = groups->groups->data[groupIndex]; // Cells in group
     74    psVector *cellSources = cells->data[cellIndex];    // Sources in cell
     75
     76    for (int i = 0; i < cellSources->n; i++) {
     77        int index = cellSources->data.S32[i];                       // Index for source of interest
     78        float flux = powf(10.0, -0.4 * mag->data.F32[index]);       // Flux of source
     79        float xSrc = x->data.F32[index], ySrc = y->data.F32[index]; // Coordinates of source
     80
     81        if (normalisePeak) {
     82            // Normalise flux
     83            pmModel *normModel = pmModelFromPSFforXY(psf, xSrc, ySrc, 1.0); // Model for normalisation
     84            if (!normModel || (normModel->flags & MODEL_MASK)) {
     85                psFree(normModel);
     86                continue;
     87            }
     88            // check that all params are valid:
     89            bool validParams = true;
     90            for (int j = 0; validParams && (j < normModel->params->n); j++) {
     91                switch (j) {
     92                  case PM_PAR_SKY:
     93                  case PM_PAR_I0:
     94                  case PM_PAR_XPOS:
     95                  case PM_PAR_YPOS:
     96                    continue;
     97                  default:
     98                    if (!isfinite(normModel->params->data.F32[j])) {
     99                        validParams = false;
     100                    }
     101                }
     102            }
     103            if (!validParams) {
     104                psFree(normModel);
     105                continue;
     106            }
     107            if (circularise && !circulariseModel(normModel)) {
     108                psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
     109                psFree(normModel);
     110                return false;
     111            }
     112
     113            flux /= normModel->modelFlux(normModel->params);
     114            psFree(normModel);
     115        }
     116
     117        pmModel *fakeModel = pmModelFromPSFforXY(psf, xSrc, ySrc, flux);
     118        if (!fakeModel || (fakeModel->flags & MODEL_MASK)) {
     119            psFree(fakeModel);
     120            continue;
     121        }
     122        // check that all params are valid:
     123        bool validParams = true;
     124        for (int j = 0; validParams && (j < fakeModel->params->n); j++) {
     125            switch (j) {
     126              case PM_PAR_SKY:
     127              case PM_PAR_I0:
     128              case PM_PAR_XPOS:
     129              case PM_PAR_YPOS:
     130                continue;
     131              default:
     132                if (!isfinite(fakeModel->params->data.F32[j])) {
     133                    validParams = false;
     134                }
     135            }
     136        }
     137        if (!validParams) {
     138            psFree(fakeModel);
     139            continue;
     140        }
     141        if (circularise && !circulariseModel(fakeModel)) {
     142            psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model.");
     143            psFree(fakeModel);
     144            return false;
     145        }
     146
     147        psTrace("psModules.camera", 10, "Adding source at %f,%f with flux %f\n",
     148                fakeModel->params->data.F32[PM_PAR_XPOS], fakeModel->params->data.F32[PM_PAR_YPOS],
     149                fakeModel->params->data.F32[PM_PAR_I0]);
     150
     151        pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate
     152        fakeSource->peak = pmPeakAlloc(xSrc, ySrc, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE);
     153        float fakeRadius = 1.0;         // Radius of fake source
     154        if (isfinite(minFlux)) {
     155            fakeRadius = PS_MAX(fakeRadius, fakeModel->modelRadius(fakeModel->params, minFlux));
     156        }
     157        if (radius > 0) {
     158            fakeRadius = PS_MAX(fakeRadius, radius);
     159        }
     160
     161        if (xOffset) {
     162            if (!pmSourceDefinePixels(fakeSource, readout, xSrc + xOffset->data.S32[index],
     163                                      ySrc + yOffset->data.S32[index], fakeRadius)) {
     164                psErrorClear();
     165                continue;
     166            }
     167            if (!pmModelAddWithOffset(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0,
     168                                      xOffset->data.S32[index], yOffset->data.S32[index])) {
     169                psErrorClear();
     170                continue;
     171            }
     172        } else {
     173            if (!pmSourceDefinePixels(fakeSource, readout, xSrc, ySrc, fakeRadius)) {
     174                psErrorClear();
     175                continue;
     176            }
     177            if (!pmModelAdd(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0)) {
     178                psErrorClear();
     179                continue;
     180            }
     181        }
     182        psFree(fakeSource);
     183        psFree(fakeModel);
     184    }
     185
     186    return true;
     187}
     188
     189/// Thread job for readoutFake()
     190static bool readoutFakeThread(psThreadJob *job)
     191{
     192    PS_ASSERT_THREAD_JOB_NON_NULL(job, false);
     193
     194    psArray *args = job->args;          // Arguments
     195
     196    pmReadout *readout = args->data[0];     // Readout of interest
     197    const pmSourceGroups *groups = args->data[1]; // Source groups
     198    const psVector *x = args->data[2];        // x coordinates
     199    const psVector *y = args->data[3];        // y coordinates
     200    const psVector *mag = args->data[4];      // Magnitudes
     201    const psVector *xOffset = args->data[5];  // Offsets in x
     202    const psVector *yOffset = args->data[6];  // Offsets in y
     203    const pmPSF *psf = args->data[7];         // PSF
     204    float minFlux = PS_SCALAR_VALUE(args->data[8], F32); // Minimum flux
     205    float radius = PS_SCALAR_VALUE(args->data[9], F32);  // Minimum radius
     206    bool circularise = PS_SCALAR_VALUE(args->data[10], U8); // Circularise PSF?
     207    bool normalisePeak = PS_SCALAR_VALUE(args->data[11], U8); // Normalise for peak?
     208    int groupIndex = PS_SCALAR_VALUE(args->data[12], S32); // Group index
     209    int cellIndex = PS_SCALAR_VALUE(args->data[13], S32);  // Cell index
     210
     211    return readoutFake(readout, groups, x, y, mag, xOffset, yOffset, psf, minFlux, radius, circularise,
     212                       normalisePeak, groupIndex, cellIndex);
     213}
     214
     215
     216bool pmReadoutFakeThreads(bool new)
     217{
     218    bool old = threaded;                // Old status, to return
     219
     220    if (!old && new) {
     221        threaded = true;
     222
     223        {
     224            psThreadTask *task = psThreadTaskAlloc("PSMODULES_READOUT_FAKE", 14);
     225            task->function = &readoutFakeThread;
     226            psThreadTaskAdd(task);
     227            psFree(task);
     228        }
     229
     230    } else if (old && !new) {
     231        threaded = false;
     232        psThreadTaskRemove("PSMODULES_READOUT_FAKE");
     233    }
     234
     235    return old;
    49236}
    50237
    … …  
    86273    psImageInit(readout->image, 0);
    87274
    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);
     275    int numThreads = threaded ? psThreadPoolSize() : 0; // Number of threads
     276    pmSourceGroups *groups = pmSourceGroupsFromVectors(readout, x, y, numThreads); // Groups of sources
     277    if (!groups) {
     278        psError(PS_ERR_UNKNOWN, false, "Unable to generate source groups");
     279        return false;
     280    }
     281
     282    if (threaded) {
     283        for (int i = 0; i < groups->groups->n; i++) {
     284            psArray *cells = groups->groups->data[i]; // Cell with sources
     285            for (int j = 0; j < cells->n; j++) {
     286                psThreadJob *job = psThreadJobAlloc("PSMODULES_READOUT_FAKE");
     287                psArray *args = job->args;
     288                psArrayAdd(args, 1, readout);
     289                psArrayAdd(args, 1, groups);
     290                // Casting away const to add to array
     291                psArrayAdd(args, 1, (psVector*)x);
     292                psArrayAdd(args, 1, (psVector*)y);
     293                psArrayAdd(args, 1, (psVector*)mag);
     294                psArrayAdd(args, 1, (psVector*)xOffset);
     295                psArrayAdd(args, 1, (psVector*)yOffset);
     296                psArrayAdd(args, 1, (pmPSF*)psf);
     297                PS_ARRAY_ADD_SCALAR(args, minFlux, PS_TYPE_F32);
     298                PS_ARRAY_ADD_SCALAR(args, radius, PS_TYPE_S32);
     299                PS_ARRAY_ADD_SCALAR(args, circularise, PS_TYPE_U8);
     300                PS_ARRAY_ADD_SCALAR(args, normalisePeak, PS_TYPE_U8);
     301                PS_ARRAY_ADD_SCALAR(args, i, PS_TYPE_S32);
     302                PS_ARRAY_ADD_SCALAR(args, j, PS_TYPE_S32);
     303
     304                if (!psThreadJobAddPending(job)) {
     305                    psFree(job);
     306                    psFree(groups);
     307                    return false;
     308                }
     309                psFree(job);
     310            }
     311            if (!psThreadPoolWait(true)) {
     312                psError(PS_ERR_UNKNOWN, false, "Error waiting for threads.");
     313                psFree(groups);
    115314                return false;
    116315            }
    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);
     316        }
     317    } else if (!readoutFake(readout, groups, x, y, mag, xOffset, yOffset, psf, minFlux, radius, circularise,
     318                            normalisePeak, 0, 0)) {
     319        psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources on readout");
     320        psFree(groups);
     321        return false;
     322    }
     323
     324    psFree(groups);
     325
     326// Set a concept value
     327#define CONCEPT_SET_S32(CONCEPTS, NAME, OLD, NEW) { \
     328        psMetadataItem *item = psMetadataLookup(CONCEPTS, NAME); \
     329        psAssert(item->type == PS_TYPE_S32, "Incorrect type: %x", item->type); \
     330        if (item->data.S32 == OLD) { \
     331            item->data.S32 = NEW; \
     332        } \
     333    }
     334
     335    if (readout->parent) {
     336        CONCEPT_SET_S32(readout->parent->concepts, "CELL.XPARITY", 0, 1);
     337        CONCEPT_SET_S32(readout->parent->concepts, "CELL.YPARITY", 0, 1);
     338        CONCEPT_SET_S32(readout->parent->concepts, "CELL.XBIN", 0, 1);
     339        CONCEPT_SET_S32(readout->parent->concepts, "CELL.YBIN", 0, 1);
    183340    }
    184341
    185342    return true;
     343
    186344}
    187345
Note: See TracChangeset for help on using the changeset viewer.