IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 8, 2011, 2:56:56 PM (15 years ago)
Author:
eugene
Message:

add unofficial PS1_V4 output format (not yet standard); add clean radial profile to find sky limit; use sky limit to control kron analysis window; fix thread / memory problem for radial photometry; fix object matching for forced psf photometry; fix object matching to get the right matches

Location:
trunk/psphot
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot

  • trunk/psphot/src/psphotRadialApertures.c

    r32348 r32633  
    2626    int num = psphotFileruleCount(config, filerule);
    2727
     28    // skip the chisq image (optionally?)
     29    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
     30    if (!status) chisqNum = -1;
     31
    2832    // loop over the available readouts
    2933    for (int i = 0; i < num; i++) {
     34        if (i == chisqNum) continue; // skip chisq image
     35
    3036        if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe, entry)) {
    3137            psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for %s entry %d", filerule, i);
    … …  
    3541    return true;
    3642}
     43
     44// these values are used by all threads repeatedly (and are not modified)
     45static psVector *aperRadii = NULL;
     46static psVector *aperRadii2 = NULL;
     47static float outerRadius = NAN;
     48static float SN_LIM = NAN;
     49static psImageMaskType maskVal = 0;
    3750
    3851// aperture-like measurements for extended sources
    … …  
    7588        nThreads = 0;
    7689    }
     90
     91    // aperRadii stores the upper bounds of the annuli
     92    // XXX keep the same name here as for the petrosian / elliptical apertures?
     93    aperRadii = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
     94    psAssert (aperRadii, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
     95    psAssert (aperRadii->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
     96
     97    outerRadius = aperRadii->data.F32[aperRadii->n - 1];
     98
     99    // save the R^2 values as well for quicker comparison
     100    aperRadii2 = psVectorAlloc(aperRadii->n, PS_TYPE_F32);
     101    for (int i = 0; i < aperRadii->n; i++) {
     102        aperRadii2->data.F32[i] = PS_SQR(aperRadii->data.F32[i]);
     103    }
     104
     105    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
     106    maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     107    assert (maskVal);
     108
     109    // S/N limit to perform full non-linear fits
     110    SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
    77111
    78112    // source analysis is done in S/N order (brightest first)
    … …  
    117151            psArrayAdd(job->args, 1, cells->data[j]); // sources
    118152            psArrayAdd(job->args, 1, AnalysisRegion);
    119             psArrayAdd(job->args, 1, recipe);
    120 
    121153            PS_ARRAY_ADD_SCALAR(job->args, entry,  PS_TYPE_S32);
    122154            PS_ARRAY_ADD_SCALAR(job->args, nEntry, PS_TYPE_S32);
    123 
    124155            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nradial
    125156
    126 // set this to 0 to run without threading
     157            // set this to 0 to run without threading
    127158# if (1)           
    128159            if (!psThreadJobAddPending(job)) {
    … …  
    138169            }
    139170            psScalar *scalar = NULL;
    140             scalar = job->args->data[6];
     171            scalar = job->args->data[5];
    141172            Nradial += scalar->data.S32;
    142173            psFree(job);
    … …  
    158189            } else {
    159190                psScalar *scalar = NULL;
    160                 scalar = job->args->data[6];
     191                scalar = job->args->data[5];
    161192                Nradial += scalar->data.S32;
    162193            }
    … …  
    166197    psFree (cellGroups);
    167198    psFree(AnalysisRegion);
     199    psFree (aperRadii2);
    168200
    169201    psLogMsg ("psphot", PS_LOG_WARN, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
    … …  
    173205bool psphotRadialApertures_Threaded (psThreadJob *job) {
    174206
    175     bool status;
    176207    int Nradial = 0;
    177208
    … …  
    180211    psArray *sources        = job->args->data[1];
    181212    psRegion *region        = job->args->data[2];
    182     psMetadata *recipe      = job->args->data[3];
    183 
    184     int entry               = PS_SCALAR_VALUE(job->args->data[4],S32); // which psf-matched image are we working on? (0 == unmatched)
    185     int nEntry              = PS_SCALAR_VALUE(job->args->data[5],S32); // total number of psf-matched images + 1 unmatched
    186 
    187     // radMax stores the upper bounds of the annuli
    188     // XXX keep the same name here as for the petrosian / elliptical apertures?
    189     psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
    190     psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
    191     psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
    192 
    193     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
    194     psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
    195     assert (maskVal);
    196 
    197     // S/N limit to perform full non-linear fits
    198     float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
    199 
    200     float outerRadius = radMax->data.F32[radMax->n - 1];
     213    int entry               = PS_SCALAR_VALUE(job->args->data[3],S32); // which psf-matched image are we working on? (0 == unmatched)
     214    int nEntry              = PS_SCALAR_VALUE(job->args->data[4],S32); // total number of psf-matched images + 1 unmatched
     215
     216    // storage for the derived pixel values (these are passed into psphotRadialApertureSource)
     217    psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
     218    psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
     219    psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
    201220
    202221    // choose the sources of interest
    … …  
    204223
    205224        pmSource *source = sources->data[i];
     225
     226        // if we have checked the source validity on the basis of the object set, then
     227        // we either skip these tests below or we skip the source completely
     228        if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_SKIP) continue;
     229        if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_KEEP) goto keepSource;
    206230
    207231        // skip PSF-like and non-astronomical objects
    … …  
    222246        if (source->peak->x > region->x1) continue;
    223247        if (source->peak->y > region->y1) continue;
     248
     249    keepSource:
    224250
    225251        // allocate pmSourceExtendedParameters, if not already defined
    … …  
    240266        }
    241267
    242         // we need to change the view for the radial aperture analysis, but we want to recover exactly
    243         // the original view; the following elements get destroyed by pmSourceRedefinePixels so save them:
    244         psImage *oldMaskObj   = psMemIncrRefCounter(source->maskObj);
    245         psImage *oldModelFlux = psMemIncrRefCounter(source->modelFlux);
    246         psImage *oldPSFimage  = psMemIncrRefCounter(source->psfImage);
    247         psRegion oldRegion    = source->region;
    248 
    249268        Nradial ++;
    250269
    251         // force source image to be a bit larger...
    252         pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
    253 
    254         if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, entry)) {
     270        if (!psphotRadialApertureSource (source, readout, entry, pixRadius2, pixFlux, pixVar)) {
    255271            psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
    256272        } else {
    … …  
    258274        }
    259275
    260         pmSourceRedefinePixelsByRegion (source, readout, oldRegion);
    261         psFree(source->maskObj);   source->maskObj   = oldMaskObj;
    262         psFree(source->modelFlux); source->modelFlux = oldModelFlux;
    263         psFree(source->psfImage);  source->psfImage  = oldPSFimage;
    264        
    265276        // re-subtract the object, leave local sky
    266277        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    267278    }
    268     psScalar *scalar = job->args->data[6];
     279    psScalar *scalar = job->args->data[5];
    269280    scalar->data.S32 = Nradial;
     281
     282    psFree (pixRadius2);
     283    psFree (pixFlux);
     284    psFree (pixVar);
    270285
    271286    return true;
    272287}
    273288
    274 bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
    275 
     289bool psphotRadialApertureSource (pmSource *source, pmReadout *readout, int entry, psVector *pixRadius2, psVector *pixFlux, psVector *pixVar) {
     290                                           
    276291    // if we are a child source, save the results to the parent source radial aperture array
    277292    psArray *radialAperSet = source->radialAper;
    … …  
    285300    radialAperSet->data[entry] = radialAper;
    286301
    287     // storage for the derived pixel values
    288     psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
    289     psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
    290     psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
     302    // find the largest aperture of interest (use only apertures with inner radii <=
     303    // source->skyRadius)
     304    int lastAp = aperRadii->n;
     305    for (int i = 0; i < aperRadii->n; i++) {
     306        if (aperRadii->data.F32[i] < source->skyRadius) continue;
     307        lastAp = i + 1;
     308        break;
     309    }
    291310
    292311    // outer-most radius for initial truncation
    293     float Rmax  = aperRadii->data.F32[aperRadii->n - 1];
     312    float Rmax  = aperRadii->data.F32[lastAp - 1];
    294313    float Rmax2 = PS_SQR(Rmax);
    295314
    296     // store the R^2 values for the apertures
    297     psVector *aperRadii2 = psVectorAlloc(aperRadii->n, PS_TYPE_F32);
    298     for (int i = 0; i < aperRadii->n; i++) {
    299         aperRadii2->data.F32[i] = PS_SQR(aperRadii->data.F32[i]);
    300     }
     315    // in this function, the operatins are relative to the full image (readout->image, etc)
    301316
    302317    float xCM = NAN, yCM = NAN;
    303318    if (pmSourcePositionUseMoments(source)) {
    304         xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
    305         yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
     319        xCM = source->moments->Mx; // index coord of peak in readout
     320        yCM = source->moments->My; // index coord of peak in readout
    306321    } else {
    307         xCM = source->peak->xf - 0.5 - source->pixels->col0; // coord of peak in subimage
    308         yCM = source->peak->yf - 0.5 - source->pixels->row0; // coord of peak in subimage
    309     }
     322        xCM = source->peak->xf; // index coord of peak in readout
     323        yCM = source->peak->yf; // index coord of peak in readout
     324    }
     325
     326    int Nx = readout->image->numCols;
     327    int Ny = readout->image->numRows;
     328
     329    pixRadius2->n = 0;
     330    pixFlux->n = 0;
     331    pixVar->n = 0;
    310332
    311333    // one pass through the pixels to select the valid pixels and calculate R^2
    312     for (int iy = 0; iy < source->pixels->numRows; iy++) {
    313 
    314         float yDiff = iy - yCM;
    315         if (fabs(yDiff) > Rmax) continue;
    316 
    317         float *vPix = source->pixels->data.F32[iy];
    318         float *vWgt = source->variance->data.F32[iy];
    319         psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy];
    320 
    321         for (int ix = 0; ix < source->pixels->numCols; ix++, vPix++, vWgt++) {
    322 
    323             if (vMsk) {
    324                 if (*vMsk & maskVal) {
    325                     vMsk++;
    326                     continue;
    327                 }
    328                 vMsk++;
    329             }
    330             if (isnan(*vPix)) continue;
    331 
    332             float xDiff = ix - xCM;
    333             if (fabs(xDiff) > Rmax) continue;
     334    for (int iy = -Rmax; iy < Rmax + 1; iy++) {
     335
     336        float yDiff = iy + 0.5 + yCM;  // y-coordinate at this offse
     337        int yPix = (int) yDiff;
     338
     339        if (yPix < 0) continue;
     340        if (yPix > Ny - 1) continue;
     341        if (fabs(iy) > Rmax) continue;
     342
     343        float *vPix = readout->image->data.F32[yPix];
     344        float *vWgt = readout->variance->data.F32[yPix];
     345        psImageMaskType  *vMsk = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[yPix];
     346
     347        for (int ix = -Rmax; ix < Rmax + 1; ix++) {
     348
     349            float xDiff = ix + 0.5 + xCM;  // x-coordinate at this offse
     350            int xPix = (int) xDiff;
     351           
     352            if (xPix < 0) continue;
     353            if (xPix > Nx - 1) continue;
     354            if (fabs(ix) > Rmax) continue;
     355           
     356            if (vMsk[xPix] & maskVal) continue;
     357            if (isnan(vPix[xPix])) continue;
    334358
    335359            // radius is just a function of (xDiff, yDiff)
    336             float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
     360            float r2  = PS_SQR(ix) + PS_SQR(iy);
    337361            if (r2 > Rmax2) continue;
    338362
    339363            psVectorAppend(pixRadius2, r2);
    340             psVectorAppend(pixFlux, *vPix);
    341             psVectorAppend(pixVar, *vWgt);
     364            psVectorAppend(pixFlux, vPix[xPix]);
     365            psVectorAppend(pixVar, vWgt[xPix]);
    342366        }
    343367    }
    … …  
    348372    psVector *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
    349373
     374    // init the apertures of interest to 0.0, the rest go to NAN
    350375    psVectorInit (flux,    0.0);
    351376    psVectorInit (fluxStd, 0.0);
    352377    psVectorInit (fluxErr, 0.0);
    353378    psVectorInit (fill,    0.0);
     379    for (int i = lastAp; i < flux->n; i++) {
     380        flux->data.F32[i] = NAN;
     381        fluxStd->data.F32[i] = NAN;
     382        fluxErr->data.F32[i] = NAN;
     383        fill->data.F32[i] = NAN;
     384    }
    354385
    355386    float *rPix2 = pixRadius2->data.F32;
    … …  
    358389        int j = 0;
    359390        float *aRad2 = aperRadii2->data.F32;
    360         for (; (*aRad2 < *rPix2) && (j < aperRadii2->n); j++, aRad2++);
    361         for (; j < aperRadii2->n; j++, aRad2++) {
     391        for (; (*aRad2 < *rPix2) && (j < lastAp); j++, aRad2++);
     392
     393        // XXX I can speed this up by only saving this single aperture
     394        for (; j < lastAp; j++, aRad2++) {
    362395            flux->data.F32[j]    += pixFlux->data.F32[i];
    363396            fluxStd->data.F32[j] += PS_SQR(pixFlux->data.F32[i]);
    … …  
    371404       2) the fractional fill factor (count of valid pixels / effective area of the aperture
    372405       3) the error on the flux within that aperture
    373      */
    374 
    375     for (int i = 0; i < flux->n; i++) {
     406    */
     407
     408    for (int i = 0; i < lastAp; i++) {
    376409        // calculate the total flux for bin 'nOut'
    377410        float Area = M_PI*aperRadii2->data.F32[i];
    … …  
    383416        // XXX report the total flux or the mask-corrected flux?
    384417        // flux->data.F32[i]    = SBmean * Area;
    385         // fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
     418        // fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / Pinx;
    386419
    387420        fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
    … …  
    393426    }
    394427   
     428# if (1)
    395429    radialAper->flux = flux;
    396430    radialAper->fluxStdev = fluxStd;
    397431    radialAper->fluxErr = fluxErr;
    398432    radialAper->fill = fill;
    399 
    400     psFree (aperRadii2);
    401     psFree (pixRadius2);
    402     psFree (pixFlux);
    403     psFree (pixVar);
     433# else
     434    // XXX TEST
     435    psFree(flux);
     436    psFree(fluxStd);
     437    psFree(fluxErr);
     438    psFree(fill);
     439# endif
    404440
    405441    return true;
    406442}
     443
     444/*** below is a test to use a sort to speed this up, not very successfully ***/
    407445
    408446static int nCalls = 0;
    … …  
    495533// *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors
    496534# define COMPARE_VECT(A,B) (radius->data.F32[A] < radius->data.F32[B])
    497 # define SWAP_VECT(TYPE,A,B) { \
    498   float tmp; \
    499   if (A != B) { \
    500     tmp = radius->data.F32[A]; \
    501     radius->data.F32[A] = radius->data.F32[B];  \
    502     radius->data.F32[B] = tmp; \
    503     tmp = pixFlux->data.F32[A]; \
    504     pixFlux->data.F32[A] = pixFlux->data.F32[B]; \
    505     pixFlux->data.F32[B] = tmp; \
    506     tmp = pixVar->data.F32[A]; \
    507     pixVar->data.F32[A] = pixVar->data.F32[B]; \
    508     pixVar->data.F32[B] = tmp; \
    509   } \
    510 }
     535# define SWAP_VECT(TYPE,A,B) {                                  \
     536        float tmp;                                              \
     537        if (A != B) {                                           \
     538            tmp = radius->data.F32[A];                          \
     539            radius->data.F32[A] = radius->data.F32[B];          \
     540            radius->data.F32[B] = tmp;                          \
     541            tmp = pixFlux->data.F32[A];                         \
     542            pixFlux->data.F32[A] = pixFlux->data.F32[B];        \
     543            pixFlux->data.F32[B] = tmp;                         \
     544            tmp = pixVar->data.F32[A];                          \
     545            pixVar->data.F32[A] = pixVar->data.F32[B];          \
     546            pixVar->data.F32[B] = tmp;                          \
     547        }                                                       \
     548    }
    511549
    512550bool psphotRadialAperturesSortFlux (psVector *radius, psVector *pixFlux, psVector *pixVar) {
Note: See TracChangeset for help on using the changeset viewer.