IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Jun 23, 2009, 7:53:02 AM (17 years ago)
Author:
eugene
Message:

updates to fix moments problem and to do multi-plane photometry

Location:
branches/eam_branches/20090522/psphot
Files:
4 added
6 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/20090522/psphot/src/Makefile.am

    r23808 r24529  
    1616
    1717# Force recompilation of psphotVersion.c, since it gets the version information
    18 psphotVersion.c: psphotVersionDefinitions.h
    19 psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE
    20         -$(RM) psphotVersionDefinitions.h
    21         $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h
    22 FORCE: ;
     18# psphotVersion.c: psphotVersionDefinitions.h
     19# psphotVersionDefinitions.h: psphotVersionDefinitions.h.in FORCE
     20#       -$(RM) psphotVersionDefinitions.h
     21#       $(SED) -e "s|@PSPHOT_VERSION@|\"$(PSPHOT_VERSION)\"|" -e "s|@PSPHOT_BRANCH@|\"$(PSPHOT_BRANCH)\"|" -e "s|@PSPHOT_SOURCE@|\"$(PSPHOT_SOURCE)\"|" psphotVersionDefinitions.h.in > psphotVersionDefinitions.h
     22# FORCE: ;
    2323
    2424libpsphot_la_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    2525libpsphot_la_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    2626
    27 bin_PROGRAMS = psphot psphotTest
     27bin_PROGRAMS = psphot psphotTest psphotMomentsStudy
    2828
    2929psphot_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
    … …  
    3434psphotTest_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
    3535psphotTest_LDADD = libpsphot.la
     36
     37psphotMomentsStudy_CFLAGS = $(PSPHOT_CFLAGS) $(PSMODULE_CFLAGS) $(PSLIB_CFLAGS)
     38psphotMomentsStudy_LDFLAGS = $(PSPHOT_LIBS) $(PSMODULE_LIBS) $(PSLIB_LIBS)
     39psphotMomentsStudy_LDADD = libpsphot.la
    3640
    3741psphot_SOURCES = \
    … …  
    5458psphotTest_SOURCES = \
    5559        psphotTest.c
     60
     61psphotMomentsStudy_SOURCES = \
     62        psphotMomentsStudy.c
    5663
    5764libpsphot_la_SOURCES = \

  • branches/eam_branches/20090522/psphot/src/psphotArguments.c

    r24144 r24529  
    105105    pmConfig *config = pmConfigRead(&argc, argv, PSPHOT_RECIPE);
    106106    if (config == NULL) {
    107         psError(PSPHOT_ERR_CONFIG, false, "Can't read site configuration");
     107      psErrorStackPrint(stderr, "Can't read site configuration");
    108108        exit(PS_EXIT_CONFIG_ERROR);
    109109    }
    … …  
    121121    if ((N = psArgumentGet (argc, argv, "-modeltest"))) {
    122122        if (argc<=N+2) {
    123           psError(PSPHOT_ERR_ARGUMENTS, true,
    124                   "Expected to see 2 more arguments; saw %d", argc - 1);
     123          psErrorStackPrint(stderr, "Expected to see 2 more arguments; saw %d", argc - 1);
    125124          exit(PS_EXIT_CONFIG_ERROR);
    126125        }
    … …  
    137136        if ((N = psArgumentGet (argc, argv, "-model"))) {
    138137            if (argc<=N+1) {
    139               psError(PSPHOT_ERR_ARGUMENTS, true,
    140                       "Expected to see 1 more argument; saw %d", argc - 1);
     138              psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    141139              exit(PS_EXIT_CONFIG_ERROR);
    142140            }
    … …  
    149147        if ((N = psArgumentGet (argc, argv, "-fitmode"))) {
    150148            if (argc<=N+1) {
    151               psError(PSPHOT_ERR_ARGUMENTS, true,
    152                       "Expected to see 1 more argument; saw %d", argc - 1);
     149              psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    153150              exit(PS_EXIT_CONFIG_ERROR);
    154151            }
    … …  
    159156        if ((N = psArgumentGet (argc, argv, "-fitset"))) {
    160157            if (argc<=N+1) {
    161               psError(PSPHOT_ERR_ARGUMENTS, true,
    162                       "Expected to see 1 more argument; saw %d", argc - 1);
     158              psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    163159              exit(PS_EXIT_CONFIG_ERROR);
    164160            }
    … …  
    172168    if ((N = psArgumentGet (argc, argv, "-photcode"))) {
    173169        if (argc<=N+1) {
    174           psError(PSPHOT_ERR_ARGUMENTS, true,
    175                   "Expected to see 1 more argument; saw %d", argc - 1);
     170          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    176171          exit(PS_EXIT_CONFIG_ERROR);
    177172        }
    … …  
    190185    if ((N = psArgumentGet (argc, argv, "-break"))) {
    191186        if (argc<=N+1) {
    192           psError(PSPHOT_ERR_ARGUMENTS, true,
    193                   "Expected to see 1 more argument; saw %d", argc - 1);
     187          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    194188          exit(PS_EXIT_CONFIG_ERROR);
    195189        }
    … …  
    202196    if ((N = psArgumentGet (argc, argv, "-region"))) {
    203197        if (argc<=N+1) {
    204           psError(PSPHOT_ERR_ARGUMENTS, true,
    205                   "Expected to see 1 more argument; saw %d", argc - 1);
     198          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    206199          exit(PS_EXIT_CONFIG_ERROR);
    207200        }
    … …  
    214207    if ((N = psArgumentGet (argc, argv, "-chip"))) {
    215208        if (argc<=N+1) {
    216           psError(PSPHOT_ERR_ARGUMENTS, true,
    217                   "Expected to see 1 more argument; saw %d", argc - 1);
     209          psErrorStackPrint(stderr, "Expected to see 1 more argument; saw %d", argc - 1);
    218210          exit(PS_EXIT_CONFIG_ERROR);
    219211        }

  • branches/eam_branches/20090522/psphot/src/psphotModelTest.c

    r21183 r24529  
    139139
    140140    // get the source moments
    141     status = pmSourceMoments (source, mRADIUS);
     141    status = pmSourceMoments (source, mRADIUS, 0.0, 1.0);
    142142    if (!status) psAbort("psSourceMoments error");
    143143    source->peak->value = source->moments->Peak;

  • branches/eam_branches/20090522/psphot/src/psphotSourceFits.c

    r24188 r24529  
    224224
    225225    // recalculate the source moments using the larger extended-source moments radius
    226     if (!pmSourceMoments (source, EXT_MOMENTS_RAD)) return false;
     226    // at this stage, skip Gaussian windowing, and allow pixels with S/N > 0.5
     227    if (!pmSourceMoments (source, EXT_MOMENTS_RAD, 0.0, 0.5)) return false;
    227228
    228229    psTrace ("psphot", 5, "trying blob...\n");

  • branches/eam_branches/20090522/psphot/src/psphotSourceStats.c

    r23219 r24529  
    8989
    9090            if (!psThreadJobAddPending(job)) {
    91                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     91                psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
    9292                psFree (job);
    9393                return NULL;
    9494            }
    9595            psFree(job);
    96 
    97 # if (0)
    98                 int nfail = 0;
    99                 int nmoments = 0;
    100                 if (!psphotSourceStats_Unthreaded (&nfail, &nmoments, cells->data[j], recipe)) {
    101                     psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
    102                     return NULL;
    103                 }
    104                 Nfail += nfail;
    105                 Nmoments += nmoments;
    106 # endif
    10796        }
    10897
    10998        // wait for the threads to finish and manage results
    11099        if (!psThreadPoolWait (false)) {
    111             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
     100            psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
    112101            return NULL;
    113102        }
    … …  
    149138    psMetadata *recipe              = job->args->data[1];
    150139
    151     float INNER    = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
    152     if (!status) return false;
    153     float MIN_SN   = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
    154     if (!status) return false;
    155     float RADIUS   = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     140    float INNER        = psMetadataLookupF32 (&status, recipe, "SKY_INNER_RADIUS");
     141    if (!status) return false;
     142    float MIN_SN       = psMetadataLookupF32 (&status, recipe, "MOMENTS_SN_MIN");
     143    if (!status) return false;
     144    float RADIUS       = psMetadataLookupF32 (&status, recipe, "PSF_MOMENTS_RADIUS");
     145    if (!status) return false;
     146    float MIN_PIXEL_SN = psMetadataLookupF32 (&status, recipe, "MOMENTS_MIN_PIXEL_SN");
     147    if (!status) return false;
     148    float SIGMA        = psMetadataLookupF32 (&status, recipe, "MOMENTS_GAUSS_SIGMA");
    156149    if (!status) return false;
    157150
    … …  
    202195
    203196        // measure basic source moments
    204         status = pmSourceMoments (source, RADIUS);
     197        status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN);
    205198        if (status) {
    206199            Nmoments ++;
    … …  
    212205        BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
    213206        psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
    214         status = pmSourceMoments (source, BIG_RADIUS);
     207        status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, MIN_PIXEL_SN);
    215208        if (status) {
    216209            source->mode |= PM_SOURCE_MODE_BIG_RADIUS;
    … …  
    292285
    293286        // measure basic source moments
    294         status = pmSourceMoments (source, RADIUS);
     287        status = pmSourceMoments (source, RADIUS, SIGMA, MIN_PIXEL_SN);
    295288        if (status) {
    296289            Nmoments ++;
    … …  
    302295        BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
    303296        psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
    304         status = pmSourceMoments (source, BIG_RADIUS);
     297        status = pmSourceMoments (source, BIG_RADIUS, 3.0*SIGMA, MIN_PIXEL_SN);
    305298        if (status) {
    306299            Nmoments ++;

  • branches/eam_branches/20090522/psphot/src/psphotVisual.c

    r23259 r24529  
    471471    bool status;
    472472    Graphdata graphdata;
    473 
    474     if (!pmVisualIsVisual()) return true;
     473    KapaSection section;
     474
     475    // if (!pmVisualIsVisual()) return true;
    475476
    476477    if (kapa3 == -1) {
    … …  
    485486    KapaClearPlots (kapa3);
    486487    KapaInitGraph (&graphdata);
    487 
    488     // there are N regions: use the first (guaranteed to exist) to get the overall limits
    489     psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, "PSF.CLUMP.REGION.000");
    490 
    491     float psfX  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
    492     float psfY  = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
    493 
    494     // examine sources to set data range
    495     graphdata.xmin = -0.05;
    496     graphdata.ymin = -0.05;
    497     graphdata.xmax = 2.0*psfX;
    498     graphdata.ymax = 2.0*psfY;
    499     KapaSetLimits (kapa3, &graphdata);
    500 
    501     KapaSetFont (kapa3, "helvetica", 14);
    502     KapaBox (kapa3, &graphdata);
    503     KapaSendLabel (kapa3, "&ss&h_x| (pixels)", KAPA_LABEL_XM);
    504     KapaSendLabel (kapa3, "&ss&h_y| (pixels)", KAPA_LABEL_YM);
    505 
     488    KapaSetFont (kapa3, "courier", 14);
     489
     490    float SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM");
     491
     492    // select the max psfX,Y values for the plot limits
     493    float Xmin = 0.0, Xmax = 0.0;
     494    float Ymin = 0.0, Ymax = 0.0;
     495    {
     496        int nRegions = psMetadataLookupS32 (&status, recipe, "PSF.CLUMP.NREGIONS");
     497        for (int n = 0; n < nRegions; n++) {
     498
     499            char regionName[64];
     500            snprintf (regionName, 64, "PSF.CLUMP.REGION.%03d", n);
     501            psMetadata *regionMD = psMetadataLookupPtr (&status, recipe, regionName);
     502
     503            float psfX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.X");
     504            float psfY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.Y");
     505            float psfdX = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DX");
     506            float psfdY = psMetadataLookupF32 (&status, regionMD, "PSF.CLUMP.DY");
     507
     508            float X0 = psfX - 4.0*psfdX;
     509            float X1 = psfX + 4.0*psfdX;
     510            float Y0 = psfY - 4.0*psfdY;
     511            float Y1 = psfY + 4.0*psfdY;
     512
     513            if (isfinite(X0)) { Xmin = PS_MAX(Xmin, X0); }
     514            if (isfinite(X1)) { Xmax = PS_MAX(Xmax, X1); }
     515            if (isfinite(Y0)) { Ymin = PS_MAX(Ymin, Y0); }
     516            if (isfinite(Y1)) { Ymax = PS_MAX(Ymax, Y1); }
     517        }
     518    }
     519
     520    // storage vectors for data to be plotted
    506521    psVector *xBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    507522    psVector *yBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     523    psVector *mBright = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    508524    psVector *xFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    509525    psVector *yFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
     526    psVector *mFaint  = psVectorAllocEmpty (sources->n, PS_TYPE_F32);
    510527
    511528    // construct the vectors
    … …  
    519536        xFaint->data.F32[nF] = source->moments->Mxx;
    520537        yFaint->data.F32[nF] = source->moments->Myy;
     538        mFaint->data.F32[nF] = -2.5*log10(source->moments->Sum);
    521539        nF++;
    522540
    523541        // XXX make this a user-defined cutoff
    524         if (source->moments->SN < 50)
     542        if (source->moments->SN < SN_LIM)
    525543            continue;
    526544
    527545        xBright->data.F32[nB] = source->moments->Mxx;
    528546        yBright->data.F32[nB] = source->moments->Myy;
     547        mBright->data.F32[nB] = -2.5*log10(source->moments->Sum);
    529548        nB++;
    530549    }
    531550    xFaint->n = nF;
    532551    yFaint->n = nF;
     552    mFaint->n = nF;
    533553
    534554    xBright->n = nB;
    535555    yBright->n = nB;
     556    mBright->n = nB;
     557
     558    // three sections: MxxMyy, MagMxx, MagMyy
     559
     560    // first section: MxxMyy
     561    section.dx = 0.75;
     562    section.dy = 0.75;
     563    section.x  = 0.00;
     564    section.y  = 0.00;
     565    section.name = psStringCopy ("MxxMyy");
     566    KapaSetSection (kapa3, &section);
     567    psFree (section.name);
     568
     569    graphdata.color = KapaColorByName ("black");
     570    graphdata.xmin = Xmin;
     571    graphdata.ymin = Ymin;
     572    graphdata.xmax = Xmax;
     573    graphdata.ymax = Ymax;
     574    KapaSetLimits (kapa3, &graphdata);
     575
     576    KapaBox (kapa3, &graphdata);
     577    KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM);
     578    KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM);
    536579
    537580    graphdata.color = KapaColorByName ("black");
    … …  
    551594    KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
    552595
    553     // XXX draw N circles to outline the clumps
     596    // second section: MagMyy
     597    section.dx = 0.75;
     598    section.dy = 0.25;
     599    section.x  = 0.00;
     600    section.y  = 0.80;
     601    section.name = psStringCopy ("MagMyy");
     602    KapaSetSection (kapa3, &section);
     603    psFree (section.name);
     604
     605    graphdata.color = KapaColorByName ("black");
     606    graphdata.xmin = -17.1;
     607    graphdata.xmax =  -7.9;
     608    graphdata.ymin = Ymin;
     609    graphdata.ymax = Ymax;
     610    KapaSetLimits (kapa3, &graphdata);
     611
     612    strcpy (graphdata.labels, "0210");
     613    KapaBox (kapa3, &graphdata);
     614    KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_XP);
     615    KapaSendLabel (kapa3, "M_yy| (pixels)", KAPA_LABEL_YM);
     616
     617    graphdata.color = KapaColorByName ("black");
     618    graphdata.ptype = 0;
     619    graphdata.size = 0.3;
     620    graphdata.style = 2;
     621    KapaPrepPlot (kapa3, nF, &graphdata);
     622    KapaPlotVector (kapa3, nF, mFaint->data.F32, "x");
     623    KapaPlotVector (kapa3, nF, yFaint->data.F32, "y");
     624
     625    graphdata.color = KapaColorByName ("red");
     626    graphdata.ptype = 2;
     627    graphdata.size = 0.5;
     628    graphdata.style = 2;
     629    KapaPrepPlot (kapa3, nB, &graphdata);
     630    KapaPlotVector (kapa3, nB, mBright->data.F32, "x");
     631    KapaPlotVector (kapa3, nB, yBright->data.F32, "y");
     632
     633    // third section: MagMxx
     634    section.dx = 0.25;
     635    section.dy = 0.75;
     636    section.x  = 0.80;
     637    section.y  = 0.00;
     638    section.name = psStringCopy ("MagMxx");
     639    KapaSetSection (kapa3, &section);
     640    psFree (section.name);
     641
     642    graphdata.color = KapaColorByName ("black");
     643    graphdata.xmin = Xmin;
     644    graphdata.xmax = Xmax;
     645    graphdata.ymin =  -7.9;
     646    graphdata.ymax = -17.1;
     647    KapaSetLimits (kapa3, &graphdata);
     648
     649    strcpy (graphdata.labels, "2001");
     650    KapaBox (kapa3, &graphdata);
     651    KapaSendLabel (kapa3, "M_xx| (pixels)", KAPA_LABEL_XM);
     652    KapaSendLabel (kapa3, "inst mag", KAPA_LABEL_YP);
     653
     654    graphdata.color = KapaColorByName ("black");
     655    graphdata.ptype = 0;
     656    graphdata.size = 0.3;
     657    graphdata.style = 2;
     658    KapaPrepPlot (kapa3, nF, &graphdata);
     659    KapaPlotVector (kapa3, nF, xFaint->data.F32, "x");
     660    KapaPlotVector (kapa3, nF, mFaint->data.F32, "y");
     661
     662    graphdata.color = KapaColorByName ("red");
     663    graphdata.ptype = 2;
     664    graphdata.size = 0.5;
     665    graphdata.style = 2;
     666    KapaPrepPlot (kapa3, nB, &graphdata);
     667    KapaPlotVector (kapa3, nB, xBright->data.F32, "x");
     668    KapaPlotVector (kapa3, nB, mBright->data.F32, "y");
     669
     670    // draw N circles to outline the clumps
    554671    {
     672        KapaSelectSection (kapa3, "MxxMyy");
     673
    555674        // draw a circle centered on psfX,Y with size of the psf limit
    556675        psVector *xLimit  = psVectorAlloc (120, PS_TYPE_F32);
    … …  
    562681        graphdata.color = KapaColorByName ("blue");
    563682        graphdata.style = 0;
     683
     684        graphdata.xmin = Xmin;
     685        graphdata.ymin = Ymin;
     686        graphdata.xmax = Xmax;
     687        graphdata.ymax = Ymax;
     688        KapaSetLimits (kapa3, &graphdata);
    564689
    565690        for (int n = 0; n < nRegions; n++) {
    … …  
    626751    psFree (xBright);
    627752    psFree (yBright);
     753    psFree (mBright);
    628754    psFree (xFaint);
    629755    psFree (yFaint);
     756    psFree (mFaint);
    630757
    631758    // pause and wait for user input:
Note: See TracChangeset for help on using the changeset viewer.