IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 38302 for trunk


Ignore:
Timestamp:
May 22, 2015, 8:40:53 AM (11 years ago)
Author:
eugene
Message:

re-cache default model if we used a different model in this analysis

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotGalaxyParams.c

    r37881 r38302  
    3434    for (int i = 0; i < num; i++) {
    3535        if (!psphotGalaxyParamsReadout (config, view, filerule, i, recipe)) {
    36             psError (PSPHOT_ERR_CONFIG, false, "failed on galaxy parameter measurements for %s entry %d", filerule, i);
     36            psError (PSPHOT_ERR_CONFIG, false, "failed on galaxy parameter measurements for %s entry %d", filerule, i);
    3737            return false;
    3838        }
    … …  
    7474    psAssert (readout, "missing readout?");
    7575
     76    // XXX EAM TEST
     77    // psphotSaveImage(NULL, readout->image, "dump.v0.fits");
     78
    7679    pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
    7780    psAssert (detections, "missing detections?");
    … …  
    8891    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
    8992    if (!status) {
    90         nThreads = 0;
     93        nThreads = 0;
    9194    }
    9295
    … …  
    100103    float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolutio
    101104    if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) {
    102         fitNsigmaConv = 5.0;
     105        fitNsigmaConv = 5.0;
    103106    }
    104107    int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
    … …  
    111114#endif
    112115
    113 
    114116    // source analysis is done in S/N order (brightest first)
    115117    sources = psArraySort (sources, pmSourceSortByFlux);
    … …  
    129131    for (int i = 0; i < cellGroups->n; i++) {
    130132
    131         psArray *cells = cellGroups->data[i];
    132 
    133         for (int j = 0; j < cells->n; j++) {
    134 
    135             // allocate a job -- if threads are not defined, this just runs the job
    136             psThreadJob *job = psThreadJobAlloc ("PSPHOT_GALAXY_PARAMS");
    137 
    138             psArrayAdd(job->args, 1, readout);
    139             psArrayAdd(job->args, 1, cells->data[j]); // sources
    140             psArrayAdd(job->args, 1, AnalysisRegion);
    141             psArrayAdd(job->args, 1, recipe);
    142 
    143             PS_ARRAY_ADD_SCALAR(job->args, skynoise, PS_TYPE_F32);
    144             PS_ARRAY_ADD_SCALAR(job->args, fitNsigmaConv, PS_TYPE_F32);
    145             PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_F32);
    146             psArrayAdd(job->args, 1, footprintImage);
    147 
    148             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
    149             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npetro
    150             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nannuli
     133        psArray *cells = cellGroups->data[i];
     134
     135        for (int j = 0; j < cells->n; j++) {
     136
     137            // allocate a job -- if threads are not defined, this just runs the job
     138            psThreadJob *job = psThreadJobAlloc ("PSPHOT_GALAXY_PARAMS");
     139
     140            psArrayAdd(job->args, 1, readout);
     141            psArrayAdd(job->args, 1, cells->data[j]); // sources
     142            psArrayAdd(job->args, 1, AnalysisRegion);
     143            psArrayAdd(job->args, 1, recipe);
     144
     145            PS_ARRAY_ADD_SCALAR(job->args, skynoise, PS_TYPE_F32);
     146            PS_ARRAY_ADD_SCALAR(job->args, fitNsigmaConv, PS_TYPE_F32);
     147            PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_F32);
     148            psArrayAdd(job->args, 1, footprintImage);
     149
     150            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
     151            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npetro
     152            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nannuli
    151153
    152154
    153155// set this to 0 to run without threading
    154156# if (1)           
    155             if (!psThreadJobAddPending(job)) {
    156                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     157            if (!psThreadJobAddPending(job)) {
     158                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    157159                psFree(AnalysisRegion);
    158                 return false;
    159             }
     160                return false;
     161            }
    160162# else
    161163            if (!psphotGalaxyParams_Threaded(job)) {
    … …  
    175177        }
    176178
    177         // wait for the threads to finish and manage results
    178         if (!psThreadPoolWait (false, true)) {
    179             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     179        // wait for the threads to finish and manage results
     180        if (!psThreadPoolWait (false, true)) {
     181            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    180182            psFree(AnalysisRegion);
    181             return false;
    182         }
    183 
    184         // we have only supplied one type of job, so we can assume the types here
    185         psThreadJob *job = NULL;
    186         while ((job = psThreadJobGetDone()) != NULL) {
    187             if (job->args->n < 1) {
    188                 fprintf (stderr, "error with job\n");
    189             } else {
     183            return false;
     184        }
     185
     186        // we have only supplied one type of job, so we can assume the types here
     187        psThreadJob *job = NULL;
     188        while ((job = psThreadJobGetDone()) != NULL) {
     189            if (job->args->n < 1) {
     190                fprintf (stderr, "error with job\n");
     191            } else {
    190192                psScalar *scalar = NULL;
    191193                scalar = job->args->data[8];
    … …  
    195197                scalar = job->args->data[10];
    196198                Nannuli += scalar->data.S32;
    197             }
    198             psFree(job);
     199            }
     200            psFree(job);
    199201        }
    200202    }
    … …  
    202204    psFree(AnalysisRegion);
    203205    psFree(footprintImage);
     206
     207    // XXX EAM TEST
     208    // psphotSaveImage(NULL, readout->image, "dump.v1.fits");
    204209
    205210    psLogMsg ("psphot", PS_LOG_WARN, "galaxy parameter analysis: %f sec for %d objects\n", psTimerMark ("psphot.galaxy"), Next);
    … …  
    281286#endif
    282287
    283     #define MARK_SOURCE 0
    284     #define MARK_MASK 1
    285     #define MARK_BG 2
     288#define MARK_SOURCE 0
     289#define MARK_MASK 1
     290#define MARK_BG 2
    286291    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
    287292
    … …  
    301306
    302307        bool testObject = FALSE;
    303         // testObject = testObject || ((fabs(source->peak->xf - 5788) < 10) && (fabs(source->peak->yf - 5711) < 10));
    304         // testObject = testObject || ((fabs(source->peak->xf -   96) < 10) && (fabs(source->peak->yf - 5550) < 10));
     308        // testObject = testObject || ((fabs(source->peak->xf - 2468) < 10) && (fabs(source->peak->yf - 3955) < 10));
     309        // testObject = testObject || ((fabs(source->peak->xf -  432) < 10) && (fabs(source->peak->yf -  781) < 10));
    305310        // testObject = testObject || ((fabs(source->peak->xf -  300) < 10) && (fabs(source->peak->yf - 5988) < 10));
    306311        if (testObject) {
    307           fprintf (stderr, "test object: %f, %f\n", source->peak->xf, source->peak->yf);
    308           testObject = TRUE;
     312            fprintf (stderr, "test object: %f, %f\n", source->peak->xf, source->peak->yf);
     313            testObject = TRUE;
    309314        }
    310315
    … …  
    319324        if (source->peak->y > region->y1) SKIP (Nskip5);
    320325
    321         // Should we skip objects that do not have an extended model?
    322         if (!source->extpars) SKIP (Nskip6);
    323         if (!source->modelFits) SKIP (Nskip7);
    324 
    325         // Flux cut to reduce sources to those that match SDSS
    326         // if (source->moments->KronFlux < 2e5) continue;
    327 
    328         pmModel *sersicModel = NULL;
    329         for (int j = 0; j < source->modelFits->n; j++) {
    330             pmModel *model = source->modelFits->data[j];
    331             if (model->type == sersicModelType) {
    332                 sersicModel = model;
    333                 break;
    334             }
    335         }
    336         if (!sersicModel) SKIP (Nskip8);
    337         if (sersicModel->flags & badModel) SKIP (Nskip8);
    338 
    339         // check out the model params and skip those that are too large
    340         psF32 *PAR = sersicModel->params->data.F32;
    341         psF32 Xo = PAR[PM_PAR_XPOS];
    342         psF32 Yo = PAR[PM_PAR_YPOS];
    343         psF32 Io = PAR[PM_PAR_I0];
    344         //
    345         // The half light radius to the major axis of the model
    346         psEllipseAxes axes = pmPSF_ModelToAxes (PAR, sersicModel->class->useReff);
    347         if (axes.major >= 100) SKIP (Nskip8);
    348 
    349         source->extpars->ghalfLightRadius = axes.major;
     326        // Should we skip objects that do not have an extended model?
     327        if (!source->extpars) SKIP (Nskip6);
     328        if (!source->modelFits) SKIP (Nskip7);
     329
     330        // Flux cut to reduce sources to those that match SDSS
     331        // if (source->moments->KronFlux < 2e5) continue;
     332
     333        pmModel *sersicModel = NULL;
     334        for (int j = 0; j < source->modelFits->n; j++) {
     335            pmModel *model = source->modelFits->data[j];
     336            if (model->type == sersicModelType) {
     337                sersicModel = model;
     338                break;
     339            }
     340        }
     341        if (!sersicModel) SKIP (Nskip8);
     342        if (sersicModel->flags & badModel) SKIP (Nskip8);
     343
     344        // check out the model params and skip those that are too large
     345        psF32 *PAR = sersicModel->params->data.F32;
     346        psF32 Xo = PAR[PM_PAR_XPOS];
     347        psF32 Yo = PAR[PM_PAR_YPOS];
     348        psF32 Io = PAR[PM_PAR_I0];
     349        //
     350        // The half light radius to the major axis of the model
     351        psEllipseAxes axes = pmPSF_ModelToAxes (PAR, sersicModel->class->useReff);
     352        if (axes.major >= 100) SKIP (Nskip8);
     353
     354        source->extpars->ghalfLightRadius = axes.major;
    350355        Next ++;
    351356
    352         // if sersicModel is not modelEXT make it so
    353         pmSourceMode saveMode = source->mode;
    354         pmSourceType saveType = source->type;
    355         bool isExtended = saveType == PM_SOURCE_TYPE_EXTENDED;
    356         pmModel *extModel = source->modelEXT;
    357 
    358         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    359             pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    360         }
    361 
    362         if (!isExtended || extModel != sersicModel) {
    363             // XXX: I'm not doing this right. For now just look at objects for which
    364             source->modelEXT = psMemIncrRefCounter (sersicModel);
    365             source->type = PM_SOURCE_TYPE_EXTENDED;
    366             source->mode |= PM_SOURCE_MODE_EXTMODEL;
    367             // do we need this?
    368             source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
    369         }
    370 
    371 
    372         // and we will integrate out to 2 * half light radius
    373         psF32 R = 2 * axes.major;
    374         psF32 R2 = PS_SQR(R);
     357        // if sersicModel is not modelEXT make it so
     358        pmSourceMode saveMode = source->mode;
     359        pmSourceType saveType = source->type;
     360        bool isExtended = saveType == PM_SOURCE_TYPE_EXTENDED;
     361        pmModel *extModel = source->modelEXT;
     362
     363        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     364            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     365        }
     366
     367        if (!isExtended || extModel != sersicModel) {
     368            // XXX: I'm not doing this right. For now just look at objects for which
     369            source->modelEXT = psMemIncrRefCounter (sersicModel);
     370            source->type = PM_SOURCE_TYPE_EXTENDED;
     371            source->mode |= PM_SOURCE_MODE_EXTMODEL;
     372            // do we need this?
     373            source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
     374        }
     375
     376
     377        // and we will integrate out to 2 * half light radius
     378        psF32 R = 2 * axes.major;
     379        psF32 R2 = PS_SQR(R);
    375380
    376381        float radius = source->peak->xf - source->pixels->col0;
    … …  
    378383        radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
    379384        radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
    380         // Make image at least 10 * half light radius
    381         radius = PS_MAX (radius, 5*axes.major);
     385        // Make image at least 10 * half light radius
     386        radius = PS_MAX (radius, 5*axes.major);
    382387        pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 2*radius);
    383388
    384         // RedefinePixels frees modelFlux
    385         if (source->modelEXT->isPCM) {
    386             pmPCMCacheModel (source, maskVal, psfSize, fitNsigmaConv);
    387         } else {
    388             pmSourceCacheModel (source, maskVal);
    389         }
     389        // RedefinePixels frees modelFlux
     390        if (source->modelEXT->isPCM) {
     391            pmPCMCacheModel (source, maskVal, psfSize, fitNsigmaConv);
     392        } else {
     393            pmSourceCacheModel (source, maskVal);
     394        }
    390395
    391396#ifdef DUMP_IMAGES
    392         sprintf(filename, "%s/%04d-pix.fits", filebase, source->id);
    393 
    394         psphotSaveImage(dummyHeader, source->pixels, filename);
    395 #endif
    396 
    397 
    398         // Create the markImage and accumlate the Flux
    399 
    400         // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
    401         // xCM, yCM from pixel coords to pixel index here.
    402         psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
    403         psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
    404 
    405         psF32 **vPix = source->pixels->data.F32;
    406         psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    407 
    408         bgPixels = psVectorRecycle(bgPixels, source->pixels->numCols*source->pixels->numRows, PS_TYPE_F32);
    409         bgPixels->n = 0;
    410         bgPixelsFlip = psVectorRecycle(bgPixelsFlip, source->pixels->numCols*source->pixels->numRows, PS_TYPE_F32);
    411         bgPixelsFlip->n = 0;
     397        sprintf(filename, "%s/%04d-pix.fits", filebase, source->id);
     398        psphotSaveImage(dummyHeader, source->pixels, filename);
     399#endif
     400
     401        // Create the markImage and accumlate the Flux
     402
     403        // center of mass in subimage.  Note: the calculation below uses pixel index, so we correct
     404        // xCM, yCM from pixel coords to pixel index here.
     405        psF32 xCM = Xo - 0.5 - source->pixels->col0; // coord of peak in subimage
     406        psF32 yCM = Yo - 0.5 - source->pixels->row0; // coord of peak in subimage
     407
     408        psF32 **vPix = source->pixels->data.F32;
     409        psImageMaskType **vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     410
     411        bgPixels = psVectorRecycle(bgPixels, source->pixels->numCols*source->pixels->numRows, PS_TYPE_F32);
     412        bgPixels->n = 0;
     413        bgPixelsFlip = psVectorRecycle(bgPixelsFlip, source->pixels->numCols*source->pixels->numRows, PS_TYPE_F32);
     414        bgPixelsFlip->n = 0;
    412415
    413416#ifdef MAKE_MARK_IMAGE
    414         // The markImage may no longer be needed
    415         markImage = psImageRecycle(markImage, source->pixels->numCols, source->pixels->numRows, PS_TYPE_U8);
    416 
    417         // assume all pixels are masked to start
    418         psImageInit(markImage, MARK_MASK);
    419         psU8 **vMark = markImage->data.U8;
    420         #define SET_MARK(_val) vMark[row][col] = _val
     417        // The markImage may no longer be needed
     418        markImage = psImageRecycle(markImage, source->pixels->numCols, source->pixels->numRows, PS_TYPE_U8);
     419
     420        // assume all pixels are masked to start
     421        psImageInit(markImage, MARK_MASK);
     422        psU8 **vMark = markImage->data.U8;
     423#define SET_MARK(_val) vMark[row][col] = _val
    421424#else
    422         #define SET_MARK(_val)
     425#define SET_MARK(_val)
    423426#endif
    424427
    425428#if (USE_FOOTPRINTS)
    426         psS32 **vFootprint = footprintImage->data.S32;
     429        psS32 **vFootprint = footprintImage->data.S32;
    427430#else
    428         // XXX use recipe value instead of hard coded 1.5
    429         psF32 threshold = skynoise * 3;
    430 #endif
    431 
    432         psF32 sumI = 0;
    433         psS32 row0 = source->pixels->row0;
    434         psS32 col0 = source->pixels->col0;
    435         for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    436             psF32 yDiff = row - yCM;
    437             // y coordinate of mirror pixel
    438             int yFlip = yCM - yDiff;
    439             if (yFlip < 0) continue;
    440             if (yFlip >= source->pixels->numRows) continue;
    441             psS32 yImage = row + row0;
    442             for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    443                 // check mask and value for this pixel
    444                 if (vMsk && (vMsk[row][col] & maskVal)) continue;
    445                 psF32 pixel = vPix[row][col];
    446                 if (isnan(pixel)) continue;
    447                 // valid pixel
    448 
    449                 psS32 xImage = col + col0;
    450 
    451                 psF32 xDiff = col - xCM;
    452 
    453                 psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff);
    454                 if (r2 <= R2) {
    455                     sumI += pixel;
    456                 }
     431        // XXX use recipe value instead of hard coded 1.5
     432        psF32 threshold = skynoise * 3;
     433#endif
     434
     435        psF32 sumI = 0;
     436        psS32 row0 = source->pixels->row0;
     437        psS32 col0 = source->pixels->col0;
     438        for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     439            psF32 yDiff = row - yCM;
     440            // y coordinate of mirror pixel
     441            int yFlip = yCM - yDiff;
     442            if (yFlip < 0) continue;
     443            if (yFlip >= source->pixels->numRows) continue;
     444            psS32 yImage = row + row0;
     445            for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     446                // check mask and value for this pixel
     447                if (vMsk && (vMsk[row][col] & maskVal)) continue;
     448                psF32 pixel = vPix[row][col];
     449                if (isnan(pixel)) continue;
     450                // valid pixel
     451
     452                psS32 xImage = col + col0;
     453
     454                psF32 xDiff = col - xCM;
     455
     456                psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff);
     457                if (r2 <= R2) {
     458                    sumI += pixel;
     459                }
    457460
    458461#if (USE_FOOTPRINTS)
    459462                // this is a pixel in the object
    460                 if (vFootprint[yImage][xImage] == source->peak->footprint->id) {
    461                     SET_MARK(MARK_SOURCE);
    462                     continue;
     463                if (vFootprint[yImage][xImage] == source->peak->footprint->id) {
     464                    SET_MARK(MARK_SOURCE);
     465                    continue;
    463466                }
    464467                // this is a pixel in another object
    465                 if (vFootprint[yImage][xImage] > 0.0) continue;
     468                if (vFootprint[yImage][xImage] > 0.0) continue;
    466469#else
    467470                // this is a pixel in an object (this or the other)
    468                 if (fabs(pixel) >= threshold) {
    469                   continue
    470                 }
    471 #endif
    472                 // we have a background pixel.
    473                 // Check whether it's mirror is a background pixel as well.
    474 
    475                 // x coordinate of mirror pixel
    476                 int xFlip = xCM - xDiff;
    477                 if (xFlip < 0) continue;
    478                 if (xFlip >= source->pixels->numCols) continue;
    479                 // check mask and value for mirror pixel
    480                 if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
    481                 psF32 mirrorPixel = vPix[yFlip][xFlip];
    482                 if (isnan(mirrorPixel)) continue;
    483 
     471                if (fabs(pixel) >= threshold) {
     472                    continue;
     473                }
     474#endif
     475                // we have a background pixel.
     476                // Check whether it's mirror is a background pixel as well.
     477       
     478                // x coordinate of mirror pixel
     479                int xFlip = xCM - xDiff;
     480                if (xFlip < 0) continue;
     481                if (xFlip >= source->pixels->numCols) continue;
     482                // check mask and value for mirror pixel
     483                if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
     484                psF32 mirrorPixel = vPix[yFlip][xFlip];
     485                if (isnan(mirrorPixel)) continue;
     486
     487       
    484488#if (USE_FOOTPRINTS)
    485                 if (vFootprint[yFlip + row0][xFlip + col0] == 0) {
     489                bool goodMirror = (vFootprint[yFlip + row0][xFlip + col0] == 0);
    486490#else
    487                 if (fabs(mirrorPixel) < threshold) {
    488 #endif
    489                     // This is a good background pixel with a good mirror Save the values
    490                     SET_MARK(MARK_BG);
    491                     psVectorAppend(bgPixels, pixel);
    492                     psVectorAppend(bgPixelsFlip, mirrorPixel);
    493                 }
    494             }
    495         }
     491                bool goodMirror = (fabs(mirrorPixel) < threshold);
     492#endif
     493                if (goodMirror) {
     494                    // This is a good background pixel with a good mirror Save the values
     495                    SET_MARK(MARK_BG);
     496                    psVectorAppend(bgPixels, pixel);
     497                    psVectorAppend(bgPixelsFlip, mirrorPixel);
     498                }
     499            }
     500        }
    496501#ifdef DUMP_IMAGES
    497         sprintf(filename, "%s/%04d-mrk.fits", filebase, source->id);
    498 
    499         psphotSaveImage(dummyHeader, markImage, filename);
    500 #endif
    501 
    502         // subtract model to create the residual image
     502        sprintf(filename, "%s/%04d-mrk.fits", filebase, source->id);
     503        psphotSaveImage(dummyHeader, markImage, filename);
     504#endif
     505
     506        // subtract model to create the residual image
    503507        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    504         if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
    505             // Hmm something has gone wrong
    506             fprintf(stderr, "failed to subtract sersic model for source %d\n", source->id);
    507             goto restoreModel;
    508         }
     508        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
     509            // Hmm something has gone wrong
     510            fprintf(stderr, "failed to subtract sersic model for source %d\n", source->id);
     511            goto restoreModel;
     512        }
    509513#ifdef DUMP_IMAGES
    510         sprintf(filename, "%s/%04d-res.fits", filebase, source->id);
    511 
    512         psphotSaveImage(dummyHeader, source->pixels, filename);
    513 #endif
    514 
    515         // don't know if these pointers could change in the background subtraction so grab them again
    516         vPix = source->pixels->data.F32;
    517         vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    518         psF32 **vVar = source->variance->data.F32;
    519         psF32 **vModel = source->modelFlux->data.F32;
    520         // Accumulate Rt and Ra
    521         psF32 sumRt = 0;
    522         psF32 sumRa = 0;
    523         int numPixels = 0;
    524         psF32 sumB = 0;
    525         psF32 sumModel = 0;
    526         int numB = 0;
    527         for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    528             psF32 yDiff = row - yCM;
    529             if (fabs(yDiff) > R) continue;
    530             // y coordinate of mirror pixel
    531             int yFlip = yCM - yDiff;
    532             if (yFlip < 0) continue;
    533             if (yFlip >= source->pixels->numRows) continue;
    534             for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    535                 // check mask and value for this pixel
    536                 if (vMsk && (vMsk[row][col] & maskVal)) continue;
    537                 psF32 pixel = vPix[row][col];
    538                 if (isnan(pixel)) continue;
    539 
    540                 psF32 xDiff = col - xCM;
    541                 psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff);
    542                 if (r2 > R2) continue;
    543 
    544                 // Measurement of bumpiness parameter excludes the pixels
    545                 // within 2 pixels of the center
    546                 if (r2 > 4) {
    547                     psF32 variance = vVar[row][col];
    548                     if (isfinite(variance)) {
    549                         psF32 modelPixel = vModel[row][col];
    550                         if (!isfinite(modelPixel)) {
    551                             // XXX: Can this happen in a good model?
    552                             fprintf(stderr, "non-finite modelPixel %4d (%4d, %4d)\n", source->id, col, row);
    553                         }
    554                         numB++;
    555                         sumB += PS_SQR(pixel) - variance;
    556                         sumModel += Io * modelPixel;
    557                     }
    558                 }
    559 
    560                 // x coordinate of mirror pixel
    561                 int xFlip = xCM - xDiff;
    562                 if (xFlip < 0) continue;
    563                 if (xFlip >= source->pixels->numCols) continue;
    564                 // check mask and value for mirror pixel
    565                 if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
    566                 psF32 mirrorPixel = vPix[yFlip][xFlip];
    567                 if (isnan(mirrorPixel)) continue;
    568 
    569                 sumRt += fabs(pixel + mirrorPixel);
    570                 sumRa += fabs(pixel - mirrorPixel);
    571                 ++numPixels;
    572             }
    573         }
    574 
    575         if (numB > 0) {
    576             source->extpars->gbumpy = sqrt(sumB/numB)  / (sumModel / numB);
    577         }
    578 
    579         if (bgPixels->n < numPixels) {
    580             fprintf(stderr, "could not find enough background pixels for source %d found %ld need %d\n", source->id, bgPixels->n, numPixels);
    581             goto restoreModel;
    582         }
    583         Ngood++;
    584 
    585         // and calculate Bt and Ba
    586         psF32 sumBt = 0;
    587         psF32 sumBa = 0;
    588         psF32 *vB = bgPixels->data.F32;
    589         psF32 *vBFlip = bgPixelsFlip->data.F32;
    590         long    length = bgPixels->n;
    591         // Randomly select numPixels background pixels
    592         for (int n = 0; n < numPixels; n++) {
    593             double rand = psRandomUniform(rng);
    594             long index = floor(rand * length);
    595             if (index >= length) {
    596                 fprintf(stderr, "random index %ld is larger than vector length %ld\n",
    597                     index, length);
    598                 continue;
    599             }
    600             sumBt += fabs(vB[index] + vBFlip[index]);
    601             sumBa += fabs(vB[index] - vBFlip[index]);
    602         }
    603 
    604         // finally save the results
    605         source->extpars->gRT = 0.5 * (sumRt - sumBt)/sumI;
    606         source->extpars->gRA = 0.5 * (sumRa - sumBa)/sumI;
    607         source->extpars->gS2 = source->extpars->gRT + source->extpars->gRA;
     514        sprintf(filename, "%s/%04d-res.fits", filebase, source->id);
     515        psphotSaveImage(dummyHeader, source->pixels, filename);
     516#endif
     517
     518        // don't know if these pointers could change in the background subtraction so grab them again
     519        vPix = source->pixels->data.F32;
     520        vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     521        psF32 **vVar = source->variance->data.F32;
     522        psF32 **vModel = source->modelFlux->data.F32;
     523        // Accumulate Rt and Ra
     524        psF32 sumRt = 0;
     525        psF32 sumRa = 0;
     526        int numPixels = 0;
     527        psF32 sumB = 0;
     528        psF32 sumModel = 0;
     529        int numB = 0;
     530        for (psS32 row = 0; row < source->pixels->numRows; row++) {
     531            psF32 yDiff = row - yCM;
     532            if (fabs(yDiff) > R) continue;
     533            // y coordinate of mirror pixel
     534            int yFlip = yCM - yDiff;
     535            if (yFlip < 0) continue;
     536            if (yFlip >= source->pixels->numRows) continue;
     537            for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     538                // check mask and value for this pixel
     539                if (vMsk && (vMsk[row][col] & maskVal)) continue;
     540                psF32 pixel = vPix[row][col];
     541                if (isnan(pixel)) continue;
     542
     543                psF32 xDiff = col - xCM;
     544                psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff);
     545                if (r2 > R2) continue;
     546
     547                // Measurement of bumpiness parameter excludes the pixels
     548                // within 2 pixels of the center
     549                if (r2 > 4) {
     550                    psF32 variance = vVar[row][col];
     551                    if (isfinite(variance)) {
     552                        psF32 modelPixel = vModel[row][col];
     553                        if (!isfinite(modelPixel)) {
     554                            // XXX: Can this happen in a good model?
     555                            fprintf(stderr, "non-finite modelPixel %4d (%4d, %4d)\n", source->id, col, row);
     556                        }
     557                        numB++;
     558                        sumB += PS_SQR(pixel) - variance;
     559                        sumModel += Io * modelPixel;
     560                    }
     561                }
     562
     563                // x coordinate of mirror pixel
     564                int xFlip = xCM - xDiff;
     565                if (xFlip < 0) continue;
     566                if (xFlip >= source->pixels->numCols) continue;
     567                // check mask and value for mirror pixel
     568                if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
     569                psF32 mirrorPixel = vPix[yFlip][xFlip];
     570                if (isnan(mirrorPixel)) continue;
     571
     572                sumRt += fabs(pixel + mirrorPixel);
     573                sumRa += fabs(pixel - mirrorPixel);
     574                ++numPixels;
     575            }
     576        }
     577
     578        if (numB > 0) {
     579            source->extpars->gbumpy = sqrt(sumB/numB)  / (sumModel / numB);
     580        }
     581
     582        if (bgPixels->n < numPixels) {
     583            fprintf(stderr, "could not find enough background pixels for source %d found %ld need %d\n", source->id, bgPixels->n, numPixels);
     584            goto restoreModel;
     585        }
     586        Ngood++;
     587
     588        // and calculate Bt and Ba
     589        psF32 sumBt = 0;
     590        psF32 sumBa = 0;
     591        psF32 *vB = bgPixels->data.F32;
     592        psF32 *vBFlip = bgPixelsFlip->data.F32;
     593        long    length = bgPixels->n;
     594        // Randomly select numPixels background pixels
     595        for (int n = 0; n < numPixels; n++) {
     596            double rand = psRandomUniform(rng);
     597            long index = floor(rand * length);
     598            if (index >= length) {
     599                fprintf(stderr, "random index %ld is larger than vector length %ld\n",
     600                        index, length);
     601                continue;
     602            }
     603            sumBt += fabs(vB[index] + vBFlip[index]);
     604            sumBa += fabs(vB[index] - vBFlip[index]);
     605        }
     606
     607        // finally save the results
     608        source->extpars->gRT = 0.5 * (sumRt - sumBt)/sumI;
     609        source->extpars->gRA = 0.5 * (sumRa - sumBa)/sumI;
     610        source->extpars->gS2 = source->extpars->gRT + source->extpars->gRA;
    608611
    609612#if (PRINT_STUFF)
    610         if (Ngood % 40 == 1) {
    611           fprintf (stderr, "   id  sumI rHL rad   sumI     gRT   gRA       sumRt      sumBt         numPixels numBGPixels\n");
    612         }
    613         fprintf (stderr, "%5d %5.1f %5.1f : %5.2e %4.1f %10.1f %6.3f %6.3f %5.2e %5.2e %10d %6ld\n",
     613        if (Ngood % 40 == 1) {
     614            fprintf (stderr, "   id  sumI rHL rad   sumI     gRT   gRA       sumRt      sumBt         numPixels numBGPixels\n");
     615        }
     616        fprintf (stderr, "%5d %5.1f %5.1f : %5.2e %4.1f %10.1f %6.3f %6.3f %5.2e %5.2e %10d %6ld\n",
    614617                 source->id, source->peak->xf, source->peak->yf,
    615618                 sumI, source->extpars->ghalfLightRadius,
    … …  
    619622#endif
    620623
    621         // Now put things back the way that they were.
    622 restoreModel:
     624        // Now put things back the way that they were.
     625    restoreModel:
    623626
    624627        if (testObject) {
    625           fprintf (stderr, "done with %d\n", source->id);
     628            fprintf (stderr, "done with %d\n", source->id);
    626629        }
    627630
    628631# define USE_BEST_FIT TRUE
    629632        // XXX TEST: Keep the sersic model regardless
    630         if (USE_BEST_FIT && ((extModel != sersicModel) || !isExtended)) {
    631             // add the sersic model back in
    632             if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
    633                 pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
    634             }
    635 
    636             // replace the extended model
    637             if (extModel != sersicModel) {
    638                 psMemDecrRefCounter (sersicModel);
    639                 source->modelEXT = extModel;
    640             }
    641             source->mode = saveMode;
    642             source->type = saveType;
    643             pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    644         }
     633        if (USE_BEST_FIT && ((extModel != sersicModel) || !isExtended)) {
     634            // add the sersic model back in
     635            if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
     636                pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
     637            }
     638
     639            // replace the extended model
     640            if (extModel != sersicModel) {
     641                psMemDecrRefCounter (sersicModel);
     642                source->modelEXT = extModel;
     643            }
     644            source->mode = saveMode;
     645            source->type = saveType;
     646
     647            // RedefinePixels frees modelFlux
     648            if (isExtended && source->modelEXT->isPCM) {
     649                pmPCMCacheModel (source, maskVal, psfSize, fitNsigmaConv);
     650            } else {
     651                pmSourceCacheModel (source, maskVal);
     652            }
     653
     654            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     655        }
    645656    }
    646657
    … …  
    650661#ifdef DUMP_IMAGES
    651662    psFree(dummyHeader);
    652     // psphotSaveImage raises an error
     663// psphotSaveImage raises an error
    653664    psErrorClear();
    654665#endif
    655666    psFree(bgPixels)
    656     psFree(bgPixelsFlip)
    657     psFree(rng);
     667        psFree(bgPixelsFlip)
     668        psFree(rng);
    658669
    659670    psScalar *scalar = NULL;
    660671
    661     // change the value of a scalar on the array (wrap this and put it in psArray.h)
     672// change the value of a scalar on the array (wrap this and put it in psArray.h)
    662673    scalar = job->args->data[8];
    663674    scalar->data.S32 = Next;
Note: See TracChangeset for help on using the changeset viewer.