IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Apr 4, 2011, 1:16:41 PM (15 years ago)
Author:
eugene
Message:

improvements to the output headers; some code organization; plug leaks; user-specified limits to input seeing ranges

Location:
trunk/ppStack
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppStack

    • Property svn:ignore

      •  

        old new  
        1717autom4te.cache
        1818Doxyfile
         19a.out.dSYM
    • Property svn:mergeinfo deleted

  • trunk/ppStack/src/ppStackConvolve.c

    r31054 r31158  
    44
    55// Update the value of a concept
    6 #define UPDATE_CONCEPT(SOURCE, NAME, VALUE) { \
    7     psMetadataItem *item = psMetadataLookup(SOURCE->concepts, NAME); \
    8     psAssert(item, "Concept should be present"); \
    9     psAssert(item->type == PS_TYPE_F32, "Concept should be F32"); \
    10     item->data.F32 = VALUE; \
    11 }
     6#define UPDATE_CONCEPT(SOURCE, NAME, VALUE) {                           \
     7        psMetadataItem *item = psMetadataLookup(SOURCE->concepts, NAME); \
     8        psAssert(item, "Concept should be present");                    \
     9        psAssert(item->type == PS_TYPE_F32, "Concept should be F32");   \
     10        item->data.F32 = VALUE;                                         \
     11    }
    1212
    1313bool ppStackConvolve(ppStackOptions *options, pmConfig *config)
    … …  
    4747    psVectorInit(renorms, NAN);
    4848
     49    psVector *satValues = psVectorAllocEmpty(num, PS_TYPE_F32); // Renormalisation values for variances
     50
    4951    psList *fpaList = psListAlloc(NULL); // List of input FPAs, for concept averaging
    5052    psList *cellList = psListAlloc(NULL); // List of input cells, for concept averaging
    … …  
    7072            psFree(cellList);
    7173            psFree(target);
     74            psFree(renorms);
     75            psFree(satValues);
    7276            return false;
    7377        }
    … …  
    8791            psFree(cellList);
    8892            psFree(target);
     93            psFree(renorms);
     94            psFree(satValues);
    8995            return false;
    9096        }
    … …  
    106112                psFree(cellList);
    107113                psFree(target);
     114                psFree(renorms);
     115                psFree(satValues);
    108116                return false;
    109117                // Non-fatal errors
    … …  
    120128                    psFree(cellList);
    121129                    psFree(target);
     130                    psFree(renorms);
     131                    psFree(satValues);
    122132                    return false;
    123133                }
    … …  
    166176            psFree(rng);
    167177            psFree(target);
     178            psFree(renorms);
     179            psFree(satValues);
    168180            return false;
    169181        }
    … …  
    177189            psFree(maskHeader);
    178190            psFree(target);
     191            psFree(renorms);
     192            psFree(satValues);
    179193            return false;
    180194        }
    … …  
    186200            psFree(rng);
    187201            psFree(target);
     202            psFree(renorms);
     203            psFree(satValues);
    188204            return false;
    189205        }
    … …  
    221237        if (options->matchZPs) {
    222238            // I think I need to take off the exposure time because we're going to set the new exposure time
     239            // Clarification: the zero point (ZP) in the header should be set such that:
     240            // M_app = m_inst + ZP + 2.5*log(exptime), where exptime in the output is sumExposure
    223241            psMetadataItem *zpItem = psMetadataLookup(inCell->parent->parent->concepts, "FPA.ZP");
     242            float inZP = zpItem->data.F32;
    224243            zpItem->data.F32 += options->norm->data.F32[i] + 2.5*log10(options->sumExposure);
    225244
    … …  
    228247
    229248            expItem = psMetadataLookup(inCell->concepts, "CELL.EXPOSURE");
     249            float inExptime = expItem->data.F32;
    230250            expItem->data.F32 = options->sumExposure;
     251
     252            // flux_out = flux_in * ten(-0.4*norm) -- save the individual saturation values
     253            psMetadataItem *satItem = psMetadataLookup(inCell->concepts, "CELL.SATURATION");
     254            float inSat = satItem->data.F32;
     255            satItem->data.F32 *= pow(10.0, -0.4*options->norm->data.F32[i]);
     256            psVectorAppend (satValues, satItem->data.F32);
     257
     258            psLogMsg("ppStack", PS_LOG_INFO, "image %d mods : zp %f -> %f, exptime %f -> %f, sat %f -> %f",
     259                     i, inZP, zpItem->data.F32, inExptime, expItem->data.F32, inSat, satItem->data.F32);
    231260        }
    232261
    … …  
    237266            psFree(rng);
    238267            psFree(target);
     268            psFree(renorms);
     269            psFree(satValues);
    239270            return false;
    240271        }
    … …  
    257288        psFree(fpaList);
    258289        psFree(cellList);
     290        psFree(renorms);
     291        psFree(satValues);
    259292        return true;
    260293    }
    … …  
    281314        psFree(fpaList);
    282315        psFree(cellList);
    283 
    284         UPDATE_CONCEPT(outFPA,  "FPA.EXPOSURE",  options->sumExposure);
    285         UPDATE_CONCEPT(outCell, "CELL.EXPOSURE", options->sumExposure);
    286         UPDATE_CONCEPT(outCell, "CELL.DARKTIME", NAN);
    287         UPDATE_CONCEPT(outFPA,  "FPA.ZP",        options->zp);
    288 
    289         UPDATE_CONCEPT(unconvFPA,  "FPA.EXPOSURE",  options->sumExposure);
    290         UPDATE_CONCEPT(unconvCell, "CELL.EXPOSURE", options->sumExposure);
    291         UPDATE_CONCEPT(unconvCell, "CELL.DARKTIME", NAN);
    292         UPDATE_CONCEPT(unconvFPA,  "FPA.ZP",        options->zp);
     316       
     317        // The best guess for an output saturation value depends on the recipe.  If we have
     318        // 'safe' on, the we require at least 2 pixels to generate a valid output pixel.  In
     319        // this case, the best value for CELL.SATURATION is the 2nd highest value in the list.
     320        // If not, it should be the higest value in the list
     321        bool mdok = false;
     322        bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels
     323        psVectorSortInPlace(satValues);
     324        float satBest = safe && satValues->n > 1 ? satValues->data.F32[1] : satValues->data.F32[0];
     325
     326        // UPDATE CELL.SATURATION here
     327        UPDATE_CONCEPT(outFPA,  "FPA.EXPOSURE",    options->sumExposure);
     328        UPDATE_CONCEPT(outCell, "CELL.EXPOSURE",   options->sumExposure);
     329        UPDATE_CONCEPT(outCell, "CELL.DARKTIME",   NAN);
     330        UPDATE_CONCEPT(outCell, "CELL.SATURATION", satBest);
     331        UPDATE_CONCEPT(outFPA,  "FPA.ZP",          options->zp);
     332        UPDATE_CONCEPT(outFPA,  "FPA.AIRMASS",     options->airmass);
     333
     334        UPDATE_CONCEPT(unconvFPA,  "FPA.EXPOSURE",    options->sumExposure);
     335        UPDATE_CONCEPT(unconvCell, "CELL.EXPOSURE",   options->sumExposure);
     336        UPDATE_CONCEPT(unconvCell, "CELL.DARKTIME",   NAN);
     337        UPDATE_CONCEPT(unconvCell, "CELL.SATURATION", satBest);
     338        UPDATE_CONCEPT(unconvFPA,  "FPA.ZP",          options->zp);
     339        UPDATE_CONCEPT(unconvFPA,  "FPA.AIRMASS",     options->airmass);
     340
     341        psLogMsg("ppStack", PS_LOG_INFO, "stack adjust metadata values : zp %f, exptime %f, sat %f", options->zp, options->sumExposure, satBest);
    293342
    294343        if (options->stats) {
    … …  
    297346            double time = psTimeToMJD(fpaTime);
    298347            psMetadataAddF64(options->stats, PS_LIST_TAIL, "MJD_OBS", PS_META_DUPLICATE_OK,
    299                                  "Average MJD_OBS of inputs", time);
    300 
    301         }
    302     }
    303 
     348                             "Average MJD_OBS of inputs", time);
     349
     350        }
     351    }
     352
     353    // XXX EAM : this may be overly harsh -- or at least it would be if I (EAM) hadn't changed
     354    // the values of matchChi2 I modified pmSubtraction.c to fit a 2nd order polynomial to the
     355    // star chisq distribution (because of systematic errors in the model being matched).  This
     356    // fit forces the mean value to be 0.0.  Perhaps we can / should exclude images which have
     357    // an excessively high value for the rms?
    304358
    305359    // Reject images out-of-hand on the basis of their match chi^2
    … …  
    316370            psError(PPSTACK_ERR_PROG, false, "Unable to sort vector.");
    317371            psFree(values);
     372            psFree(renorms);
     373            psFree(satValues);
    318374            return false;
    319375        }
    … …  
    356412        psErrorClear();
    357413        psWarning("No good images survived convolution stage.");
     414        psFree(renorms);
     415        psFree(satValues);
    358416        return true;
    359417    }
    … …  
    366424    }
    367425    psFree(renorms);
     426    psFree(satValues);
    368427
    369428    if (options->stats) {
Note: See TracChangeset for help on using the changeset viewer.