Index: trunk/psModules/src/camera/pmReadoutFake.c =================================================================== --- trunk/psModules/src/camera/pmReadoutFake.c (revision 26076) +++ trunk/psModules/src/camera/pmReadoutFake.c (revision 26202) @@ -23,11 +23,14 @@ #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? @@ -47,4 +50,186 @@ } 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; } @@ -86,102 +271,57 @@ psImageInit(readout->image, 0); - 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); + 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); 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 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); - } + } + } 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); return true; + } Index: trunk/psModules/src/camera/pmReadoutFake.h =================================================================== --- trunk/psModules/src/camera/pmReadoutFake.h (revision 26076) +++ trunk/psModules/src/camera/pmReadoutFake.h (revision 26202) @@ -12,4 +12,11 @@ #include #include + +/// Set threading +/// +/// Returns old threading state +bool pmReadoutFakeThreads( + bool new // New threading state + ); /// Generate a fake readout from vectors