Index: trunk/psModules/src/camera/pmReadoutFake.c =================================================================== --- trunk/psModules/src/camera/pmReadoutFake.c (revision 26202) +++ trunk/psModules/src/camera/pmReadoutFake.c (revision 26216) @@ -23,14 +23,11 @@ #include "pmSourceUtils.h" #include "pmModelUtils.h" -#include "pmSourceGroups.h" #include "pmReadoutFake.h" +#define MODEL_TYPE "PS_MODEL_RGAUSS" // Type of model to use #define MAX_AXIS_RATIO 20.0 // Maximum axis ratio for PSF model #define MODEL_MASK (PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE | \ PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_LIMITS) // Mask to apply to models - - -static bool threaded = false; // Running threaded? @@ -50,186 +47,4 @@ } return pmPSF_AxesToModel(params, axes); -} - -/// Generate fake sources on a readout -static bool readoutFake(pmReadout *readout, // Readout of interest - const pmSourceGroups *groups, // Source groups - const psVector *x, // x coordinates - const psVector *y, // y coordinates - const psVector *mag, // Magnitudes - const psVector *xOffset, // Offsets in x - const psVector *yOffset, // Offsets in y - const pmPSF *psf, // PSF - float minFlux, // Minimum flux - float radius, // Minimum radius - bool circularise, // Circularise PSF? - bool normalisePeak, // Normalise sources for peak? - int groupIndex, // Group index - int cellIndex // Cell index - ) -{ - psArray *cells = groups->groups->data[groupIndex]; // Cells in group - psVector *cellSources = cells->data[cellIndex]; // Sources in cell - - for (int i = 0; i < cellSources->n; i++) { - int index = cellSources->data.S32[i]; // Index for source of interest - float flux = powf(10.0, -0.4 * mag->data.F32[index]); // Flux of source - float xSrc = x->data.F32[index], ySrc = y->data.F32[index]; // Coordinates of source - - if (normalisePeak) { - // Normalise flux - pmModel *normModel = pmModelFromPSFforXY(psf, xSrc, ySrc, 1.0); // Model for normalisation - if (!normModel || (normModel->flags & MODEL_MASK)) { - psFree(normModel); - continue; - } - // check that all params are valid: - bool validParams = true; - for (int j = 0; validParams && (j < normModel->params->n); j++) { - switch (j) { - case PM_PAR_SKY: - case PM_PAR_I0: - case PM_PAR_XPOS: - case PM_PAR_YPOS: - continue; - default: - if (!isfinite(normModel->params->data.F32[j])) { - validParams = false; - } - } - } - if (!validParams) { - psFree(normModel); - continue; - } - if (circularise && !circulariseModel(normModel)) { - psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model."); - psFree(normModel); - return false; - } - - flux /= normModel->modelFlux(normModel->params); - psFree(normModel); - } - - pmModel *fakeModel = pmModelFromPSFforXY(psf, xSrc, ySrc, flux); - if (!fakeModel || (fakeModel->flags & MODEL_MASK)) { - psFree(fakeModel); - continue; - } - // check that all params are valid: - bool validParams = true; - for (int j = 0; validParams && (j < fakeModel->params->n); j++) { - switch (j) { - case PM_PAR_SKY: - case PM_PAR_I0: - case PM_PAR_XPOS: - case PM_PAR_YPOS: - continue; - default: - if (!isfinite(fakeModel->params->data.F32[j])) { - validParams = false; - } - } - } - if (!validParams) { - psFree(fakeModel); - continue; - } - if (circularise && !circulariseModel(fakeModel)) { - psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model."); - psFree(fakeModel); - return false; - } - - psTrace("psModules.camera", 10, "Adding source at %f,%f with flux %f\n", - fakeModel->params->data.F32[PM_PAR_XPOS], fakeModel->params->data.F32[PM_PAR_YPOS], - fakeModel->params->data.F32[PM_PAR_I0]); - - pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate - fakeSource->peak = pmPeakAlloc(xSrc, ySrc, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE); - float fakeRadius = 1.0; // Radius of fake source - if (isfinite(minFlux)) { - fakeRadius = PS_MAX(fakeRadius, fakeModel->modelRadius(fakeModel->params, minFlux)); - } - if (radius > 0) { - fakeRadius = PS_MAX(fakeRadius, radius); - } - - if (xOffset) { - if (!pmSourceDefinePixels(fakeSource, readout, xSrc + xOffset->data.S32[index], - ySrc + yOffset->data.S32[index], fakeRadius)) { - psErrorClear(); - continue; - } - if (!pmModelAddWithOffset(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0, - xOffset->data.S32[index], yOffset->data.S32[index])) { - psErrorClear(); - continue; - } - } else { - if (!pmSourceDefinePixels(fakeSource, readout, xSrc, ySrc, fakeRadius)) { - psErrorClear(); - continue; - } - if (!pmModelAdd(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0)) { - psErrorClear(); - continue; - } - } - psFree(fakeSource); - psFree(fakeModel); - } - - return true; -} - -/// Thread job for readoutFake() -static bool readoutFakeThread(psThreadJob *job) -{ - PS_ASSERT_THREAD_JOB_NON_NULL(job, false); - - psArray *args = job->args; // Arguments - - pmReadout *readout = args->data[0]; // Readout of interest - const pmSourceGroups *groups = args->data[1]; // Source groups - const psVector *x = args->data[2]; // x coordinates - const psVector *y = args->data[3]; // y coordinates - const psVector *mag = args->data[4]; // Magnitudes - const psVector *xOffset = args->data[5]; // Offsets in x - const psVector *yOffset = args->data[6]; // Offsets in y - const pmPSF *psf = args->data[7]; // PSF - float minFlux = PS_SCALAR_VALUE(args->data[8], F32); // Minimum flux - float radius = PS_SCALAR_VALUE(args->data[9], F32); // Minimum radius - bool circularise = PS_SCALAR_VALUE(args->data[10], U8); // Circularise PSF? - bool normalisePeak = PS_SCALAR_VALUE(args->data[11], U8); // Normalise for peak? - int groupIndex = PS_SCALAR_VALUE(args->data[12], S32); // Group index - int cellIndex = PS_SCALAR_VALUE(args->data[13], S32); // Cell index - - return readoutFake(readout, groups, x, y, mag, xOffset, yOffset, psf, minFlux, radius, circularise, - normalisePeak, groupIndex, cellIndex); -} - - -bool pmReadoutFakeThreads(bool new) -{ - bool old = threaded; // Old status, to return - - if (!old && new) { - threaded = true; - - { - psThreadTask *task = psThreadTaskAlloc("PSMODULES_READOUT_FAKE", 14); - task->function = &readoutFakeThread; - psThreadTaskAdd(task); - psFree(task); - } - - } else if (old && !new) { - threaded = false; - psThreadTaskRemove("PSMODULES_READOUT_FAKE"); - } - - return old; } @@ -271,57 +86,102 @@ psImageInit(readout->image, 0); - int numThreads = threaded ? psThreadPoolSize() : 0; // Number of threads - pmSourceGroups *groups = pmSourceGroupsFromVectors(readout, x, y, numThreads); // Groups of sources - if (!groups) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate source groups"); - return false; - } - - if (threaded) { - for (int i = 0; i < groups->groups->n; i++) { - psArray *cells = groups->groups->data[i]; // Cell with sources - for (int j = 0; j < cells->n; j++) { - psThreadJob *job = psThreadJobAlloc("PSMODULES_READOUT_FAKE"); - psArray *args = job->args; - psArrayAdd(args, 1, readout); - psArrayAdd(args, 1, groups); - // Casting away const to add to array - psArrayAdd(args, 1, (psVector*)x); - psArrayAdd(args, 1, (psVector*)y); - psArrayAdd(args, 1, (psVector*)mag); - psArrayAdd(args, 1, (psVector*)xOffset); - psArrayAdd(args, 1, (psVector*)yOffset); - psArrayAdd(args, 1, (pmPSF*)psf); - PS_ARRAY_ADD_SCALAR(args, minFlux, PS_TYPE_F32); - PS_ARRAY_ADD_SCALAR(args, radius, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(args, circularise, PS_TYPE_U8); - PS_ARRAY_ADD_SCALAR(args, normalisePeak, PS_TYPE_U8); - PS_ARRAY_ADD_SCALAR(args, i, PS_TYPE_S32); - PS_ARRAY_ADD_SCALAR(args, j, PS_TYPE_S32); - - if (!psThreadJobAddPending(job)) { - psFree(job); - psFree(groups); - return false; - } - psFree(job); - } - if (!psThreadPoolWait(true)) { - psError(PS_ERR_UNKNOWN, false, "Error waiting for threads."); - psFree(groups); + for (long i = 0; i < numSources; i++) { + float flux = powf(10.0, -0.4 * mag->data.F32[i]); // Flux of source + float xSrc = x->data.F32[i], ySrc = y->data.F32[i]; // Coordinates of source + + if (normalisePeak) { + // Normalise flux + pmModel *normModel = pmModelFromPSFforXY(psf, xSrc, ySrc, 1.0); // Model for normalisation + if (!normModel || (normModel->flags & MODEL_MASK)) { + psFree(normModel); + continue; + } + // check that all params are valid: + bool validParams = true; + for (int n = 0; validParams && (n < normModel->params->n); n++) { + if (n == PM_PAR_SKY) continue; + if (n == PM_PAR_I0) continue; + if (n == PM_PAR_XPOS) continue; + if (n == PM_PAR_YPOS) continue; + if (!isfinite(normModel->params->data.F32[n])) validParams = false; + } + if (!validParams) { + psFree(normModel); + continue; + } + if (circularise && !circulariseModel(normModel)) { + psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model."); + psFree(normModel); return false; } - } - } else if (!readoutFake(readout, groups, x, y, mag, xOffset, yOffset, psf, minFlux, radius, circularise, - normalisePeak, 0, 0)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate fake sources on readout"); - psFree(groups); - return false; - } - - psFree(groups); + + flux /= normModel->modelFlux(normModel->params); + psFree(normModel); + } + + pmModel *fakeModel = pmModelFromPSFforXY(psf, xSrc, ySrc, flux); + if (!fakeModel || (fakeModel->flags & MODEL_MASK)) { + psFree(fakeModel); + continue; + } + // check that all params are valid: + bool validParams = true; + for (int n = 0; validParams && (n < fakeModel->params->n); n++) { + if (n == PM_PAR_SKY) continue; + if (n == PM_PAR_I0) continue; + if (n == PM_PAR_XPOS) continue; + if (n == PM_PAR_YPOS) continue; + if (!isfinite(fakeModel->params->data.F32[n])) validParams = false; + } + if (!validParams) { + psFree(fakeModel); + continue; + } + if (circularise && !circulariseModel(fakeModel)) { + psError(PS_ERR_UNKNOWN, false, "Unable to circularise PSF model."); + psFree(fakeModel); + return false; + } + + psTrace("psModules.camera", 10, "Adding source at %f,%f with flux %f\n", + fakeModel->params->data.F32[PM_PAR_XPOS], fakeModel->params->data.F32[PM_PAR_YPOS], + fakeModel->params->data.F32[PM_PAR_I0]); + + pmSource *fakeSource = pmSourceAlloc(); // Fake source to generate + fakeSource->peak = pmPeakAlloc(xSrc, ySrc, fakeModel->params->data.F32[PM_PAR_I0], PM_PEAK_LONE); + float fakeRadius = 1.0; // Radius of fake source + if (isfinite(minFlux)) { + fakeRadius = PS_MAX(fakeRadius, fakeModel->modelRadius(fakeModel->params, minFlux)); + } + if (radius > 0) { + fakeRadius = PS_MAX(fakeRadius, radius); + } + + if (xOffset) { + if (!pmSourceDefinePixels(fakeSource, readout, xSrc + xOffset->data.S32[i], + ySrc + yOffset->data.S32[i], fakeRadius)) { + psErrorClear(); + continue; + } + if (!pmModelAddWithOffset(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0, + xOffset->data.S32[i], yOffset->data.S32[i])) { + psErrorClear(); + continue; + } + } else { + if (!pmSourceDefinePixels(fakeSource, readout, xSrc, ySrc, fakeRadius)) { + psErrorClear(); + continue; + } + if (!pmModelAdd(fakeSource->pixels, NULL, fakeModel, PM_MODEL_OP_FULL, 0)) { + psErrorClear(); + continue; + } + } + psFree(fakeSource); + psFree(fakeModel); + } return true; - } Index: trunk/psModules/src/camera/pmReadoutFake.h =================================================================== --- trunk/psModules/src/camera/pmReadoutFake.h (revision 26202) +++ trunk/psModules/src/camera/pmReadoutFake.h (revision 26216) @@ -12,11 +12,4 @@ #include #include - -/// Set threading -/// -/// Returns old threading state -bool pmReadoutFakeThreads( - bool new // New threading state - ); /// Generate a fake readout from vectors