IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 37626 for trunk/psphot/src


Ignore:
Timestamp:
Nov 18, 2014, 10:35:52 AM (12 years ago)
Author:
bills
Message:

Use footprints to select background pixels instead of a threshold.
Cache sersic model before subtracting it

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotGalaxyParams.c

    r37466 r37626  
    11# include "psphotInternal.h"
     2
     3#define USE_FOOTPRINTS 1
    24
    35// measure some properties on the exended objects
    … …  
    9597    }
    9698
     99    float fitNsigmaConv = psMetadataLookupF32 (&status, recipe, "EXT_FIT_NSIGMA_CONV"); // number of sigma for the convolutio
     100    if (!status || !isfinite(fitNsigmaConv) || fitNsigmaConv <= 0) {
     101        fitNsigmaConv = 5.0;
     102    }
     103    int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
     104
     105    psImage *footprintImage = NULL;
     106#ifdef USE_FOOTPRINTS
     107    footprintImage = psImageAlloc(readout->image->numCols, readout->image->numRows, PS_TYPE_S32);
     108    psImageInit(footprintImage, 0);
     109    pmSetFootprintArrayIDsForImage(footprintImage, detections->footprints, true);
     110#endif
     111
     112
    97113    // source analysis is done in S/N order (brightest first)
    98114    sources = psArraySort (sources, pmSourceSortByFlux);
    … …  
    125141
    126142            PS_ARRAY_ADD_SCALAR(job->args, skynoise, PS_TYPE_F32);
     143            PS_ARRAY_ADD_SCALAR(job->args, fitNsigmaConv, PS_TYPE_F32);
     144            PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_F32);
     145            psArrayAdd(job->args, 1, footprintImage);
    127146
    128147            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
    129148            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npetro
    130149            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nannuli
     150
    131151
    132152// set this to 0 to run without threading
    … …  
    168188            } else {
    169189                psScalar *scalar = NULL;
    170                 scalar = job->args->data[5];
     190                scalar = job->args->data[8];
    171191                Next += scalar->data.S32;
    172                 scalar = job->args->data[6];
     192                scalar = job->args->data[9];
    173193                Npetro += scalar->data.S32;
    174                 scalar = job->args->data[7];
     194                scalar = job->args->data[10];
    175195                Nannuli += scalar->data.S32;
    176196            }
    … …  
    180200    psFree (cellGroups);
    181201    psFree(AnalysisRegion);
     202    psFree(footprintImage);
    182203
    183204    psLogMsg ("psphot", PS_LOG_WARN, "galaxy parameter analysis: %f sec for %d objects\n", psTimerMark ("psphot.galaxy"), Next);
    … …  
    211232    int Npetro = 0;
    212233    int Nannuli = 0;
     234    int Ngood = 0;
    213235
    214236    // arguments: readout, sources, models, region, psfSize, maskVal, markVal
    … …  
    218240    psMetadata *recipe      = job->args->data[3];
    219241
     242#ifndef USE_FOOTPRINTS
    220243    float skynoise          = PS_SCALAR_VALUE(job->args->data[4],F32);
     244#endif
     245    float fitNsigmaConv     = PS_SCALAR_VALUE(job->args->data[5],F32);
     246    float psfSize           = PS_SCALAR_VALUE(job->args->data[6],F32);
     247    psImage *footprintImage = job->args->data[7];
    221248
    222249    bool status;
    … …  
    227254    assert (maskVal);
    228255
     256    // don't try and use bad models
     257    pmModelStatus badModel = PM_MODEL_STATUS_NONE;
     258    badModel |= PM_MODEL_STATUS_BADARGS;
     259    badModel |= PM_MODEL_STATUS_OFFIMAGE;
     260    badModel |= PM_MODEL_STATUS_NAN_CHISQ;
     261    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GUESS;
     262    badModel |= PM_MODEL_SERSIC_PCM_FAIL_GRID;
     263    badModel |= PM_MODEL_PCM_FAIL_GUESS;
     264
    229265    pmModelType sersicModelType = pmModelClassGetType("PS_MODEL_SERSIC");
    230266
    231267    psVector *bgPixels = NULL;
    232268    psVector *bgPixelsFlip = NULL;
     269#define dontMAKE_MARK_IMAGE 1
    233270#ifdef MAKE_MARK_IMAGE
    234271    psImage *markImage = NULL;
    … …  
    239276    #define MARK_BG 2
    240277    psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS);
     278
     279#define dontDUMP_IMAGES 1
     280#ifdef DUMP_IMAGES
     281    char *filebase = "/data/ippc18.0/home/bills/debug/psphot/images";
     282    char filename[80];
     283    psMetadata *dummyHeader = psMetadataAlloc();
     284#endif
    241285
    242286    // choose the sources of interest
    … …  
    260304        if (!source->extpars) SKIP (Nskip6);
    261305        if (!source->modelFits) SKIP (Nskip7);
     306
     307        // Flux cut to reduce sources to those that match SDSS
     308        // if (source->moments->KronFlux < 2e5) continue;
    262309
    263310        pmModel *sersicModel = NULL;
    … …  
    270317        }
    271318        if (!sersicModel) SKIP (Nskip8);
     319        if (sersicModel->flags & badModel) SKIP (Nskip8);
    272320        Next ++;
    273321
    … …  
    283331
    284332        if (!isExtended || extModel != sersicModel) {
     333            // XXX: I'm not doing this right. For now just look at objects for which
    285334            source->modelEXT = psMemIncrRefCounter (sersicModel);
    286335            source->type = PM_SOURCE_TYPE_EXTENDED;
    … …  
    288337            // do we need this?
    289338            source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
    290         }
    291        
    292         // force source image to be a bit larger...
    293         float radius = source->peak->xf - source->pixels->col0;
    294         radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
    295         radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
    296         radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
    297         pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
    298 
    299         // Do the work here
     339            if (source->modelEXT->isPCM) {
     340                pmPCMCacheModel (source, maskVal, psfSize, fitNsigmaConv);
     341            } else {
     342                pmSourceCacheModel (source, maskVal);
     343            }
     344        }
    300345
    301346        psF32 *PAR = sersicModel->params->data.F32;
    … …  
    306351        psEllipseAxes axes = pmPSF_ModelToAxes (PAR, sersicModel->class->useReff);
    307352        source->extpars->ghalfLightRadius = axes.major;
     353        // and we will integrate out to 2 * half light radius
    308354        psF32 R = 2 * axes.major;
    309355        psF32 R2 = PS_SQR(R);
     356
     357        float radius = source->peak->xf - source->pixels->col0;
     358        radius = PS_MAX (radius, source->peak->yf - source->pixels->row0);
     359        radius = PS_MAX (radius, source->pixels->numRows - source->peak->yf + source->pixels->row0);
     360        radius = PS_MAX (radius, source->pixels->numCols - source->peak->xf + source->pixels->col0);
     361        // Make image at least 10 * half light radius
     362        radius = PS_MAX (radius, 5*axes.major);
     363        pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 2*radius);
     364
     365
     366#ifdef DUMP_IMAGES
     367        sprintf(filename, "%s/%04d-pix.fits", filebase, source->id);
     368
     369        psphotSaveImage(dummyHeader, source->pixels, filename);
     370#endif
     371
    310372
    311373        // Create the markImage and accumlate the Flux
    … …  
    336398#endif
    337399
     400#ifdef USE_FOOTPRINTS
     401        psS32 **vFootprint = footprintImage->data.S32;
     402#else
    338403        // XXX use recipe value instead of hard coded 1.5
    339         psF32 threshold = skynoise * 1.5;
     404        psF32 threshold = skynoise * 3;
     405#endif
    340406
    341407        psF32 sumI = 0;
     408        psS32 row0 = source->pixels->row0;
     409        psS32 col0 = source->pixels->col0;
    342410        for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    343411            psF32 yDiff = row - yCM;
     412            // y coordinate of mirror pixel
     413            int yFlip = yCM - yDiff;
     414            if (yFlip < 0) continue;
     415            if (yFlip >= source->pixels->numRows) continue;
     416            psS32 yImage = row + row0;
     417            for (psS32 col = 0; col < source->pixels->numCols ; col++) {
     418                // check mask and value for this pixel
     419                if (vMsk && (vMsk[row][col] & maskVal)) continue;
     420                psF32 pixel = vPix[row][col];
     421                if (isnan(pixel)) continue;
     422                // valid pixel
     423
     424                psS32 xImage = col + col0;
     425
     426                psF32 xDiff = col - xCM;
     427
     428                psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff);
     429                if (r2 <= R2) {
     430                    sumI += pixel;
     431                }
     432
     433#ifdef USE_FOOTPRINTS
     434                if (vFootprint[yImage][xImage]) {
     435#else
     436                if (fabs(pixel) >= threshold) {
     437#endif
     438                    SET_MARK(MARK_SOURCE);
     439                    continue;
     440                }
     441                // we have a background pixel.
     442                // Check whether it's mirror is a background pixel as well.
     443
     444                // x coordinate of mirror pixel
     445                int xFlip = xCM - xDiff;
     446                if (xFlip < 0) continue;
     447                if (xFlip >= source->pixels->numCols) continue;
     448                // check mask and value for mirror pixel
     449                if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
     450                psF32 mirrorPixel = vPix[yFlip][xFlip];
     451                if (isnan(mirrorPixel)) continue;
     452
     453#ifdef USE_FOOTPRINTS
     454                if (vFootprint[yFlip + row0][xFlip + col0] == 0) {
     455#else
     456                if (fabs(mirrorPixel) < threshold) {
     457#endif
     458                    // This is a good background pixel with a good mirror Save the values
     459                    SET_MARK(MARK_BG);
     460                    psVectorAppend(bgPixels, pixel);
     461                    psVectorAppend(bgPixelsFlip, mirrorPixel);
     462                }
     463            }
     464        }
     465#ifdef DUMP_IMAGES
     466        sprintf(filename, "%s/%04d-mrk.fits", filebase, source->id);
     467
     468        psphotSaveImage(dummyHeader, markImage, filename);
     469#endif
     470
     471        // subtract model to create the residual image
     472        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     473        if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
     474            // Hmm something has gone wrong
     475            fprintf(stderr, "failed to subtract sersic model for source %d\n", source->id);
     476            goto restoreModel;
     477        }
     478#ifdef DUMP_IMAGES
     479        sprintf(filename, "%s/%04d-res.fits", filebase, source->id);
     480
     481        psphotSaveImage(dummyHeader, source->pixels, filename);
     482#endif
     483
     484        // don't know if these pointers could change in the background subtraction so grab them again
     485        vPix = source->pixels->data.F32;
     486        vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
     487        // Accumulate Rt and Ra
     488        psF32 sumRt = 0;
     489        psF32 sumRa = 0;
     490        int numPixels = 0;
     491        for (psS32 row = 0; row < source->pixels->numRows ; row++) {
     492            psF32 yDiff = row - yCM;
     493            if (fabs(yDiff) > R) continue;
    344494            // y coordinate of mirror pixel
    345495            int yFlip = yCM - yDiff;
    … …  
    351501                psF32 pixel = vPix[row][col];
    352502                if (isnan(pixel)) continue;
    353                 // valid pixel
    354503
    355504                psF32 xDiff = col - xCM;
    356 
    357505                psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff);
    358                 if (r2 <= R2) {
    359                     sumI += pixel;
    360                 }
    361 
    362 
    363                 if (fabs(pixel) >= threshold) {
    364                     SET_MARK(MARK_SOURCE);
    365                     continue;
    366                 }
    367                 // we have a background pixel.
    368                 // Check whether it's mirror is a background pixel as well.
     506                if (r2 > R2) continue;
    369507
    370508                // x coordinate of mirror pixel
    … …  
    376514                psF32 mirrorPixel = vPix[yFlip][xFlip];
    377515                if (isnan(mirrorPixel)) continue;
    378                 if (fabs(mirrorPixel) < threshold) {
    379                     // This is a good background pixel with a good mirror Save the values
    380                     SET_MARK(MARK_BG);
    381                     psVectorAppend(bgPixels, pixel);
    382                     psVectorAppend(bgPixelsFlip, mirrorPixel);
    383                 }
    384             }
    385         }
    386 
    387         // subtract model to create the residual image
    388         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    389         if (!(source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED)) {
    390             // HMM something has gone wrong
    391             goto restoreModel;
    392         }
    393 
    394         // don't know if these could change in the background subtraction
    395         vPix = source->pixels->data.F32;
    396         vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA;
    397         // Accumulate Rt and Ra
    398         psF32 sumRt = 0;
    399         psF32 sumRa = 0;
    400         int numPixels = 0;
    401         for (psS32 row = 0; row < source->pixels->numRows ; row++) {
    402             psF32 yDiff = row - yCM;
    403             // y coordinate of mirror pixel
    404             int yFlip = yCM - yDiff;
    405             if (yFlip < 0) continue;
    406             if (yFlip >= source->pixels->numRows) continue;
    407             for (psS32 col = 0; col < source->pixels->numCols ; col++) {
    408                 // check mask and value for this pixel
    409                 if (vMsk && (vMsk[row][col] & maskVal)) continue;
    410                 psF32 pixel = vPix[row][col];
    411                 if (isnan(pixel)) continue;
    412 
    413                 psF32 xDiff = col - xCM;
    414                 psF32 r2 = PS_SQR(xDiff) + PS_SQR(yDiff);
    415                 if (r2 > R2) continue;
    416 
    417                 // x coordinate of mirror pixel
    418                 int xFlip = xCM - xDiff;
    419                 if (xFlip < 0) continue;
    420                 if (xFlip >= source->pixels->numCols) continue;
    421                 // check mask and value for mirror pixel
    422                 if (vMsk && (vMsk[yFlip][xFlip] & maskVal)) continue;
    423                 psF32 mirrorPixel = vPix[yFlip][xFlip];
    424                 if (isnan(mirrorPixel)) continue;
    425516
    426517                sumRt += fabs(pixel + mirrorPixel);
    … …  
    429520            }
    430521        }
     522
     523        if (bgPixels->n < numPixels) {
     524            fprintf(stderr, "could not find enough background pixels for source %d found %ld need %d\n", source->id, bgPixels->n, numPixels);
     525            goto restoreModel;
     526        }
     527        Ngood++;
     528
    431529        // and calculate Bt and Ba
    432530        psF32 sumBt = 0;
    … …  
    439537            double rand = psRandomUniform(rng);
    440538            long index = floor(rand * length);
    441             if (index >= length) continue;
     539            if (index >= length) {
     540                fprintf(stderr, "random index %ld is larger than vector length %ld\n",
     541                    index, length);
     542                continue;
     543            }
    442544            sumBt += fabs(vB[index] + vBFlip[index]);
    443545            sumBa += fabs(vB[index] - vBFlip[index]);
    … …  
    449551        source->extpars->gS2 = source->extpars->gRT + source->extpars->gRA;
    450552
     553#ifdef PRINT_STUFF
     554        if (Ngood % 40 == 1) {
     555            printf("   id  sumI rHL rad   sumI     gRT   gRA       sumRt      sumBt         numPixels numBGPixels\n");
     556        }
     557        printf("%5d %5.2e %3.1f %10.1f %6.3f %6.3f %5.2e %5.2e %10d %6ld\n",
     558            source->id, sumI, source->extpars->ghalfLightRadius, radius, source->extpars->gRT, source->extpars->gRA, sumRt, sumBt, numPixels, bgPixels->n);
     559#endif
    451560
    452561        // Now put things back the way that they were.
    … …  
    468577            pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    469578        }
    470 
    471579    }
    472580
    473581#ifdef MAKE_MARK_IMAGE
    474582    psFree(markImage);
     583#endif
     584#ifdef DUMP_IMAGES
     585    psFree(dummyHeader);
     586    // psphotSaveImage raises an error
     587    psErrorClear();
    475588#endif
    476589    psFree(bgPixels)
    … …  
    481594
    482595    // change the value of a scalar on the array (wrap this and put it in psArray.h)
    483     scalar = job->args->data[5];
     596    scalar = job->args->data[8];
    484597    scalar->data.S32 = Next;
    485598
    486     scalar = job->args->data[6];
     599    scalar = job->args->data[9];
    487600    scalar->data.S32 = Npetro;
    488601
    489     scalar = job->args->data[7];
     602    scalar = job->args->data[10];
    490603    scalar->data.S32 = Nannuli;
    491604
Note: See TracChangeset for help on using the changeset viewer.