IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 20, 2009, 1:21:38 PM (16 years ago)
Author:
bills
Message:

back out changes from -r26202

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r26202 r26216  
    2323#include "pmSourceUtils.h"
    2424#include "pmModelUtils.h"
    25 #include "pmSourceGroups.h"
    2625
    2726#include "pmReadoutFake.h"
    2827
     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 
    34 static bool threaded = false;           // Running threaded?
    3532
    3633
    … …  
    5047    }
    5148    return pmPSF_AxesToModel(params, axes);
    52 }
    53 
    54 /// Generate fake sources on a readout
    55 static 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()
    188 static 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 
    214 bool 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;
    23449}
    23550
    … …  
    27186    psImageInit(readout->image, 0);
    27287
    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);
     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);
    312115                return false;
    313116            }
    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);
     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    }
    323184
    324185    return true;
    325 
    326186}
    327187
Note: See TracChangeset for help on using the changeset viewer.