IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 20263


Ignore:
Timestamp:
Oct 18, 2008, 2:26:45 PM (18 years ago)
Author:
beaumont
Message:

Updates to psastro to include plotting calls

Location:
branches/cnb_branch_20081011/psastro/src
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • branches/cnb_branch_20081011/psastro/src/Makefile.am

    r20067 r20263  
    5353        psastroErrorCodes.c         \
    5454        psastroVersion.c            \
     55        psastroVisual.c             \
    5556        psastroDefineFiles.c        \
    5657        psastroAnalysis.c           \

  • branches/cnb_branch_20081011/psastro/src/psastro.h

    r20036 r20263  
    7575psString          psastroVersionLong(void);
    7676
     77// psastroVisual functions
     78bool psastroSetVisual (bool mode);
     79bool psastroVisualClose();
     80bool psastroVisualPlotLuminosityFunction (psVector *lnMag, psVector *Mag, pmLumFunc *lumFunc, pmLumFunc *rawFunc);
     81bool psastroVisualPlotRawStars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe);
     82bool psastroVisualPlotRefStars (psArray *refstars, psMetadata *recipe);
     83bool psastroVisualPlotRemoveClumps (psArray *input, psImage *count, int scale, float limit);
     84bool psastroVisualPlotOneChipFit (psArray *rawstars, psArray *refstars, psArray *match, psMetadata *recipe);
     85bool psastroVisualPlotAstromGuessCheck (psVector *cornerPo, psVector *cornerQo, psVector *cornerPn, psVector *cornerQn, psVector *cornerPd, psVector *cornerQd);
     86
    7787// demo plots
    7888bool              psastroPlotRawstars (psArray *rawstars, pmFPA *fpa, pmChip *chip, psMetadata *recipe);

  • branches/cnb_branch_20081011/psastro/src/psastroArguments.c

    r20067 r20263  
    3737        psArgumentRemove (N, &argc, argv);
    3838    }
    39    
     39
    4040    // apply the chip correction based on the reference astrometry?
    4141    if ((N = psArgumentGet (argc, argv, "-fixchips"))) {
    … …  
    5151    status = pmConfigFileSetsMD (config->arguments, &argc, argv, "ASTROM.MODEL", "-astrommodel", "-astrommodellist");
    5252    if (status) {
    53         // if supplied, assume -fixchips is desired
     53        // if supplied, assume -fixchips is desired
    5454        psMetadataAddBool (config->arguments, PS_LIST_TAIL, "PSASTRO.FIX.CHIPS", PS_META_REPLACE, "", true);
    5555    }
    … …  
    8181    }
    8282
     83    // run in visual mode?
     84    if ((N = psArgumentGet (argc, argv, "-visual"))) {
     85        psArgumentRemove (N, &argc, argv);
     86        psastroSetVisual (true);
     87    }
     88
    8389    // dump the configuration to a file?
    8490    if ((N = psArgumentGet (argc, argv, "-dumpconfig"))) {

  • branches/cnb_branch_20081011/psastro/src/psastroAstromGuess.c

    r19977 r20263  
    2929    psMetadata *recipe  = psMetadataLookupPtr (NULL, config->recipes, PSASTRO_RECIPE);
    3030    if (!recipe) {
    31         psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!");
    32         return false;
     31        psError(PSASTRO_ERR_CONFIG, true, "Can't find PSASTRO recipe!");
     32        return false;
    3333    }
    3434
    … …  
    3636    bool useModel = psMetadataLookupBool (&status, config->arguments, "PSASTRO.USE.MODEL");
    3737    if (!status) {
    38         useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL");
     38        useModel = psMetadataLookupBool (&status, recipe, "PSASTRO.USE.MODEL");
    3939    }
    4040
    … …  
    4242    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    4343    if (!input) {
    44         psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
    45         return false;
     44        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     45        return false;
    4646    }
    4747
    … …  
    4949    double pixelScale = psMetadataLookupF32 (&status, recipe, "PSASTRO.PIXEL.SCALE");
    5050    if (!status) {
    51         psError(PS_ERR_IO, true, "Failed to lookup pixel scale");
    52         return false;
    53     } 
     51        psError(PS_ERR_IO, true, "Failed to lookup pixel scale");
     52        return false;
     53    }
    5454
    5555    psVector *cornerL = psVectorAllocEmpty (100, PS_TYPE_F32);
    … …  
    6767    bool bilevelAstrometry = false;
    6868    if (!useModel) {
    69         psastroAstromGuessSetFPA (fpa, &bilevelAstrometry);
     69        psastroAstromGuessSetFPA (fpa, &bilevelAstrometry);
    7070    }
    7171
    … …  
    7575        if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
    7676
    77         if (!useModel) {
    78             if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue;
    79         }
     77        if (!useModel) {
     78            if (!psastroAstromGuessSetChip (fpa, chip, view, pixelScale, bilevelAstrometry)) continue;
     79        }
    8080
    8181        if (newFPA) {
    8282            newFPA = false;
    83             while (fpa->toSky->R <        0) fpa->toSky->R += 2.0*M_PI;
    84             while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI;
     83            while (fpa->toSky->R <        0) fpa->toSky->R += 2.0*M_PI;
     84            while (fpa->toSky->R > 2.0*M_PI) fpa->toSky->R -= 2.0*M_PI;
    8585            RAminSky = fpa->toSky->R - M_PI;
    8686            RAmaxSky = fpa->toSky->R + M_PI;
    8787        }
    8888
    89         // report and save the current best guess for the chip 0,0 pixel coordinates
    90         {
    91             psPlane ptCH, ptFP, ptTP;
    92             psSphere ptSky;
    93 
    94             ptCH.x = 0;
    95             ptCH.y = 0;
    96             psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
    97             psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
    98             psDeproject (&ptSky, &ptTP, fpa->toSky);
    99             psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
    100 
    101             psVectorAppend (cornerL, ptFP.x);
    102             psVectorAppend (cornerM, ptFP.y);
    103             psVectorAppend (cornerP, ptTP.x);
    104             psVectorAppend (cornerQ, ptTP.y);
    105             psVectorAppend (cornerR, ptSky.r);
    106             psVectorAppend (cornerD, ptSky.d);
    107         }
     89        // report and save the current best guess for the chip 0,0 pixel coordinates
     90        {
     91            psPlane ptCH, ptFP, ptTP;
     92            psSphere ptSky;
     93
     94            ptCH.x = 0;
     95            ptCH.y = 0;
     96            psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
     97            psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
     98            psDeproject (&ptSky, &ptTP, fpa->toSky);
     99            psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
     100
     101            psVectorAppend (cornerL, ptFP.x);
     102            psVectorAppend (cornerM, ptFP.y);
     103            psVectorAppend (cornerP, ptTP.x);
     104            psVectorAppend (cornerQ, ptTP.y);
     105            psVectorAppend (cornerR, ptSky.r);
     106            psVectorAppend (cornerD, ptSky.d);
     107        }
    108108
    109109        // apply the new WCS guess data to all of the data in the readouts
    … …  
    119119                if (rawstars == NULL) { continue; }
    120120
    121                 *nStars += rawstars->n;
     121                *nStars += rawstars->n;
    122122                for (int i = 0; i < rawstars->n; i++) {
    123123                    pmAstromObj *raw = rawstars->data[i];
    … …  
    138138                }
    139139
    140                 // dump or plot the resulting projected positions
    141                 if (psTraceGetLevel("psastro.dump") > 0) {
    142                     psastroDumpRawstars (rawstars, fpa, chip);
    143                 }
    144 
    145                 if (psTraceGetLevel("psastro.plot") > 0) {
    146                     psastroPlotRawstars (rawstars, fpa, chip, recipe);
    147                 }
     140                // dump or plot the resulting projected positions
     141                if (psTraceGetLevel("psastro.dump") > 0) {
     142                    psastroDumpRawstars (rawstars, fpa, chip);
     143                }
     144
     145                psastroVisualPlotRawStars(rawstars, fpa, chip, recipe);
     146
     147                if (psTraceGetLevel("psastro.plot") > 0) {
     148                    psastroPlotRawstars (rawstars, fpa, chip, recipe);
     149                }
    148150            }
    149151        }
    … …  
    155157    psMetadataAddS32 (recipe, PS_LIST_TAIL, "NTOTSTAR",  PS_META_REPLACE, "", *nStars);
    156158    if (*nStars == 0) {
    157         psLogMsg ("psastro", 2, "no sources available for astrometry\n");
    158         psFree (view);
    159         return true;
    160     }
    161 
    162     psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n", 
    163               DEG_RAD*RAmin, DEG_RAD*DECmin,
    164               DEG_RAD*RAmax, DEG_RAD*DECmax);
     159        psLogMsg ("psastro", 2, "no sources available for astrometry\n");
     160        psFree (view);
     161        return true;
     162    }
     163
     164    psLogMsg ("psastro", 2, "loaded raw data from %f,%f to %f,%f\n",
     165              DEG_RAD*RAmin, DEG_RAD*DECmin,
     166              DEG_RAD*RAmax, DEG_RAD*DECmax);
    165167
    166168    psMetadataAddF32 (recipe, PS_LIST_TAIL, "RA_MIN",  PS_META_REPLACE, "", RAmin);
    … …  
    201203    pmHDU *hdu = pmFPAviewThisHDU (view, fpa);
    202204    if (bilevelAstrometry) {
    203         if (!pmAstromReadBilevelChip (chip, hdu->header)) {
    204             psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
    205             return false;
    206         }
     205        if (!pmAstromReadBilevelChip (chip, hdu->header)) {
     206            psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
     207            return false;
     208        }
    207209    } else {
    208         if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) {
    209             psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
    210             return false;
    211         }
     210        if (!pmAstromReadWCS (fpa, chip, hdu->header, pixelScale)) {
     211            psWarning("Could not get WCS information from header for chip %d, skipping", view->chip);
     212            return false;
     213        }
    212214    }
    213215    return true;
    … …  
    223225    // load mosaic-level astrometry?
    224226    if (phu) {
    225         char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
    226         if (ctype) {
    227             *bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
    228         }
     227        char *ctype = psMetadataLookupStr (NULL, phu->header, "CTYPE1");
     228        if (ctype) {
     229            *bilevelAstrometry = !strcmp (&ctype[4], "-DIS");
     230        }
    229231    }
    230232    if (*bilevelAstrometry) {
    231         pmAstromReadBilevelMosaic (fpa, phu->header);
    232     } 
     233        pmAstromReadBilevelMosaic (fpa, phu->header);
     234    }
    233235    psFree (view);
    234236    return true;
    … …  
    243245    pmFPAfile *input = psMetadataLookupPtr (NULL, config->files, "PSASTRO.INPUT");
    244246    if (!input) {
    245         psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
    246         return false;
     247        psError(PSASTRO_ERR_CONFIG, true, "Can't find input data");
     248        return false;
    247249    }
    248250
    … …  
    273275        if (!chip->process || !chip->file_exists || !chip->data_exists) { continue; }
    274276
    275         psPlane ptCH, ptFP, ptTP;
    276         psSphere ptSky;
    277 
    278         ptCH.x = 0;
    279         ptCH.y = 0;
    280         psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
    281         psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
    282         psDeproject (&ptSky, &ptTP, fpa->toSky);
    283         psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
    284 
    285         // new corner locations based on the calibrated astrometry
    286         psVectorAppend (cornerLn, ptFP.x);
    287         psVectorAppend (cornerMn, ptFP.y);
    288         psVectorAppend (cornerPn, ptTP.x);
    289         psVectorAppend (cornerQn, ptTP.y);
    290         psVectorAppend (cornerRn, ptSky.r);
    291         psVectorAppend (cornerDn, ptSky.d);
     277        psPlane ptCH, ptFP, ptTP;
     278        psSphere ptSky;
     279
     280        ptCH.x = 0;
     281        ptCH.y = 0;
     282        psPlaneTransformApply (&ptFP, chip->toFPA, &ptCH);
     283        psPlaneTransformApply (&ptTP, fpa->toTPA, &ptFP);
     284        psDeproject (&ptSky, &ptTP, fpa->toSky);
     285        psLogMsg ("psastro", 3, "0,0 pix for chip %3d = %f,%f\n", view->chip, DEG_RAD*ptSky.r, DEG_RAD*ptSky.d);
     286
     287        // new corner locations based on the calibrated astrometry
     288        psVectorAppend (cornerLn, ptFP.x);
     289        psVectorAppend (cornerMn, ptFP.y);
     290        psVectorAppend (cornerPn, ptTP.x);
     291        psVectorAppend (cornerQn, ptTP.y);
     292        psVectorAppend (cornerRn, ptSky.r);
     293        psVectorAppend (cornerDn, ptSky.d);
    292294    }
    293295
    … …  
    298300
    299301    for (int i = 0; i < cornerRo->n; i++) {
    300        
    301         psPlane ptTP;
    302         psSphere ptSky;
    303 
    304         ptSky.r = cornerRo->data.F32[i];
    305         ptSky.d = cornerDo->data.F32[i];
    306 
    307         psProject (&ptTP, &ptSky, fpa->toSky);
    308         psVectorAppend (cornerPs, ptTP.x);
    309         psVectorAppend (cornerQs, ptTP.y);
     302
     303        psPlane ptTP;
     304        psSphere ptSky;
     305
     306        ptSky.r = cornerRo->data.F32[i];
     307        ptSky.d = cornerDo->data.F32[i];
     308
     309        psProject (&ptTP, &ptSky, fpa->toSky);
     310        psVectorAppend (cornerPs, ptTP.x);
     311        psVectorAppend (cornerQs, ptTP.y);
    310312    }
    311313
    … …  
    313315    map->x->coeffMask[1][1] = PS_POLY_MASK_SET;
    314316    map->y->coeffMask[1][1] = PS_POLY_MASK_SET;
    315    
     317
    316318    psVectorFitPolynomial2D (map->x, NULL, 0, cornerPn, NULL, cornerPs, cornerQs);
    317319    psVectorFitPolynomial2D (map->y, NULL, 0, cornerQn, NULL, cornerPs, cornerQs);
    318    
     320
    319321    // apply the linear fit...
    320322    psVector *cornerPf = psPolynomial2DEvalVector (map->x, cornerPs, cornerQs);
    … …  
    325327    psVector *cornerQd = (psVector *) psBinaryOp (NULL, cornerQn, "-", cornerQf);
    326328
     329    psastroVisualPlotAstromGuessCheck (cornerPo, cornerQo, cornerPn, cornerQn, cornerPd, cornerQd);
     330
    327331    psStats *statsP = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    328332    psStats *statsQ = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV);
    … …  
    333337    float angle = atan2 (map->y->coeff[1][0], map->x->coeff[1][0]);
    334338    float scale = hypot (map->y->coeff[1][0], map->x->coeff[1][0]);
    335    
     339
    336340    psLogMsg ("psastro", 3, "boresite offset  : %f,%f\n", map->x->coeff[0][0], map->y->coeff[0][0]);
    337341    psLogMsg ("psastro", 3, "boresite angle   : %f, scale: %f", angle*PS_DEG_RAD, scale);
    … …  
    341345    psMetadata *header = psMetadataLookupMetadata (&status, input->fpa->analysis, "PSASTRO.HEADER");
    342346    if (!header) {
    343         header = psMetadataAlloc();
    344         psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", header);
    345         psFree (header);  // drop this reference
     347        header = psMetadataAlloc();
     348        psMetadataAddMetadata (input->fpa->analysis, PS_LIST_TAIL, "PSASTRO.HEADER",  PS_META_REPLACE, "psastro header stats", header);
     349        psFree (header);  // drop this reference
    346350    }
    347351
    … …  
    354358
    355359    if (0) {
    356         FILE *f = fopen ("corners.dat", "w");
    357         for (int i = 0; i < cornerRo->n; i++) {
    358             fprintf (f, "%10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f  |  %10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f\n",
    359                      cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i],
    360                      cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]);
    361         }
    362         fclose (f);
     360        FILE *f = fopen ("corners.dat", "w");
     361        for (int i = 0; i < cornerRo->n; i++) {
     362            fprintf (f, "%10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f  |  %10.6f %10.6f  %9.2f %9.2f  %9.2f %9.2f\n",
     363                     cornerRn->data.F32[i], cornerDn->data.F32[i], cornerPn->data.F32[i], cornerQn->data.F32[i], cornerLn->data.F32[i], cornerMn->data.F32[i],
     364                     cornerRo->data.F32[i], cornerDo->data.F32[i], cornerPo->data.F32[i], cornerQo->data.F32[i], cornerLo->data.F32[i], cornerMo->data.F32[i]);
     365        }
     366        fclose (f);
    363367    }
    364368
    … …  
    381385    psFree (map);
    382386    psFree (view);
    383    
     387
    384388
    385389    return true;

  • branches/cnb_branch_20081011/psastro/src/psastroCleanup.c

    r14641 r20263  
    1111    pmConceptsDone ();
    1212    pmConfigDone ();
     13    psastroVisualClose ();
    1314    fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, stdout, false), "psastro");
    1415    // fprintf (stderr, "found %d leaks at %s\n", psMemCheckLeaks (0, NULL, NULL, false), "psastro");

  • branches/cnb_branch_20081011/psastro/src/psastroLoadRefstars.c

    r20067 r20263  
    138138        psastroDumpRefstars (refstars, "refstars.dat");
    139139    }
     140
     141    psastroVisualPlotRefStars (refstars, recipe);
    140142
    141143    if (psTraceGetLevel("psastro.plot") > 0) {

  • branches/cnb_branch_20081011/psastro/src/psastroLuminosityFunction.c

    r20037 r20263  
    129129    lumFunc->sPeak = sPeak;
    130130
    131 #if 0
    132     psastroLuminosityFunctionPlot(lnMag, Mag, lumFunc, rawFunc);
    133 #endif
     131    psastroVisualPlotLuminosityFunction(lnMag, Mag, lumFunc, rawFunc);
    134132
    135133    psFree (lnMag);

  • branches/cnb_branch_20081011/psastro/src/psastroOneChipFit.c

    r16070 r20263  
    1111
    1212    // default value for match/fit : radius is in pixels
    13     REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32); 
     13    REQUIRED_RECIPE_VALUE (double RADIUS, "PSASTRO.MATCH.RADIUS", F32);
    1414
    1515    // run the match/fit sequence NITER times
    16     REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32); 
     16    REQUIRED_RECIPE_VALUE (int nIter, "PSASTRO.MATCH.FIT.NITER", S32);
    1717
    1818    // correct radius to FP units (physical pixel scale in microns per pixel)
    19     REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32); 
     19    REQUIRED_RECIPE_VALUE (double pixelScale, "PSASTRO.PIXEL.SCALE", F32);
    2020    RADIUS *= pixelScale;
    2121
    … …  
    3232
    3333    for (int iter = 0; iter < nIter; iter++) {
    34        
    35         char name[128];
    3634
    37         sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter);
    38         float radius = psMetadataLookupF32 (&status, recipe, name);
    39         radius *= pixelScale;
    40         if (!status || (radius == 0.0)) {
    41             radius = RADIUS;
    42         }
     35        char name[128];
     36
     37        sprintf (name, "PSASTRO.MATCH.RADIUS.N%d", iter);
     38        float radius = psMetadataLookupF32 (&status, recipe, name);
     39        radius *= pixelScale;
     40        if (!status || (radius == 0.0)) {
     41            radius = RADIUS;
     42        }
    4343
    4444
    45         // use small radius to match stars
    46         match = pmAstromRadiusMatchFP (rawstars, refstars, radius);
    47         if (match == NULL) {
    48             psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
    49             return false;
    50         }
     45        // use small radius to match stars
     46        match = pmAstromRadiusMatchFP (rawstars, refstars, radius);
     47        if (match == NULL) {
     48            psLogMsg ("psastro", 3, "failed to find radius-matched sources\n");
     49            return false;
     50        }
    5151
    52         // modify the order to correspond to the actual number of matched stars:
    53         if ((match->n < 11) && (order >= 3)) order = 2;
    54         if ((match->n <  7) && (order >= 2)) order = 1;
    55         if ((match->n <  4) && (order >= 1)) order = 0;
    56         if (order < 1) {
    57             psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n);
    58             psFree (match);
    59             return false;
    60         }
     52        // modify the order to correspond to the actual number of matched stars:
     53        if ((match->n < 11) && (order >= 3)) order = 2;
     54        if ((match->n <  7) && (order >= 2)) order = 1;
     55        if ((match->n <  4) && (order >= 1)) order = 0;
     56        if (order < 1) {
     57            psLogMsg ("psastro", 3, "insufficient stars or invalid order: %ld stars", match->n);
     58            psFree (match);
     59            return false;
     60        }
    6161
    62         // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
    63         psFree (chip->toFPA);
    64         chip->toFPA = psPlaneTransformAlloc (order, order);
    65         for (int i = 0; i <= chip->toFPA->x->nX; i++) {
    66             for (int j = 0; j <= chip->toFPA->x->nY; j++) {
    67                 if (i + j > order) {
    68                     chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET;
    69                     chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
    70                 }
    71             }
    72         }
     62        // create output toFPA; set masks appropriate to the Elixir DVO astrometry format
     63        psFree (chip->toFPA);
     64        chip->toFPA = psPlaneTransformAlloc (order, order);
     65        for (int i = 0; i <= chip->toFPA->x->nX; i++) {
     66            for (int j = 0; j <= chip->toFPA->x->nY; j++) {
     67                if (i + j > order) {
     68                    chip->toFPA->x->coeffMask[i][j] = PS_POLY_MASK_SET;
     69                    chip->toFPA->y->coeffMask[i][j] = PS_POLY_MASK_SET;
     70                }
     71            }
     72        }
    7373
    74         // XXX allow statitic to be set by the user
    75         // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
    76         fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    77         fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA");
    78         fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER");
     74        // XXX allow statitic to be set by the user
     75        // fitStats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV);
     76        fitStats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     77        fitStats->clipSigma = psMetadataLookupF32 (&status, recipe, "PSASTRO.CHIP.NSIGMA");
     78        fitStats->clipIter = psMetadataLookupS32 (&status, recipe, "PSASTRO.CHIP.NITER");
    7979
    80         // improved fit for astrometric terms
    81         results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
    82         if (!results) {
    83             psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
    84             psFree (match);
    85             psFree (fitStats);
    86             return false;
    87         }
    88    
    89         // determine fromFPA transformation and apply new transformation to raw & ref stars
    90         psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
    91    
    92         // toSky converts from FPA & TPA units (microns) to sky units (radians)
    93         float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
    94         // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
    95         float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
    96         int astNstar = results->yStats->clippedNvalues;
    97         psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
     80        // improved fit for astrometric terms
     81        results = pmAstromMatchFit (chip->toFPA, rawstars, refstars, match, fitStats);
     82        if (!results) {
     83            psLogMsg ("psastro", 3, "failed to perform the matched fit\n");
     84            psFree (match);
     85            psFree (fitStats);
     86            return false;
     87        }
    9888
    99         if (iter < nIter - 1) {
    100             psFree (fitStats);
    101             psFree (results);
    102             psFree (match);
    103         }
     89        // determine fromFPA transformation and apply new transformation to raw & ref stars
     90        psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
     91
     92        // toSky converts from FPA & TPA units (microns) to sky units (radians)
     93        float plateScale = 0.5*(fpa->toSky->Xs + fpa->toSky->Ys)*3600.0*PM_DEG_RAD;
     94        // float astError = 0.5*(results->xStats->clippedStdev + results->yStats->clippedStdev) * plateScale;
     95        float astError = 0.5*(results->xStats->robustStdev + results->yStats->robustStdev) * plateScale;
     96        int astNstar = results->yStats->clippedNvalues;
     97        psLogMsg ("psastro", PS_LOG_INFO, "pass %d, error: %f arcsec, Nstars: %d", iter, astError, astNstar);
     98
     99        if (iter < nIter - 1) {
     100            psFree (fitStats);
     101            psFree (results);
     102            psFree (match);
     103        }
    104104    }
    105105
    … …  
    122122    if (astError > maxError) {
    123123        psLogMsg("psastro", PS_LOG_INFO, "residual error is too large, failed to find a solution: %f > %f", astError, maxError);
    124         validSolution = false;
     124        validSolution = false;
    125125    }
    126126    if (astNstar < minNstar) {
    127127        psLogMsg("psastro", PS_LOG_INFO, "solution uses too few stars: %d < %d", astNstar, minNstar);
    128         validSolution = false;
     128        validSolution = false;
    129129    }
    130130
    … …  
    133133    psMetadataAddF32 (updates, PS_LIST_TAIL, "CERROR",   PS_META_REPLACE, "astrometry error (arcsec)", astError);
    134134    if (validSolution) {
    135         psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
    136         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
     135        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", astError/sqrt(astNstar));
     136        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", astNstar);
    137137    } else {
    138         psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
    139         psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
     138        psMetadataAddF32 (updates, PS_LIST_TAIL, "CPRECISE", PS_META_REPLACE, "astrometry precision (arcsec)", 0.0);
     139        psMetadataAddS32 (updates, PS_LIST_TAIL, "NASTRO",   PS_META_REPLACE, "number of astrometry stars", 0);
    140140    }
    141141    psMetadataAddF32 (updates, PS_LIST_TAIL, "EQUINOX",  PS_META_REPLACE, "equinox of ref catalog", 2000.0); // XXX this is bogus: should be defined based on equinox of refstars
    … …  
    143143    // XXX drop from here : determine fromFPA transformation and apply new transformation to raw & ref stars
    144144    // psastroUpdateChipToFPA (fpa, chip, rawstars, refstars);
    145    
     145
    146146    // XXX check if we correctly applied the new transformation:
    147147    if (psTraceGetLevel("psastro.dump") > 0) {
    148         psastroDumpRawstars (rawstars, fpa, chip);
    149         psastroDumpMatchedStars ("match.dat", rawstars, refstars, match);
    150         psastroDumpStars (refstars, "refstars.cal.dat");
     148        psastroDumpRawstars (rawstars, fpa, chip);
     149        psastroDumpMatchedStars ("match.dat", rawstars, refstars, match);
     150        psastroDumpStars (refstars, "refstars.cal.dat");
    151151    }
    152152
     153    psastroVisualPlotOneChipFit (rawstars, refstars, match, recipe);
     154
    153155    if (psTraceGetLevel("psastro.plot") > 0) {
    154         psastroPlotOneChipFit (rawstars, refstars, match, recipe);
     156        psastroPlotOneChipFit (rawstars, refstars, match, recipe);
    155157    }
    156158

  • branches/cnb_branch_20081011/psastro/src/psastroRemoveClumps.c

    r17109 r20263  
    1111    float Ymax = obj->FP->y;
    1212    for (int i = 0; i < input->n; i++) {
    13         obj = (pmAstromObj *)input->data[i];
    14         if (!isfinite(obj->FP->x)) continue;
    15         if (!isfinite(obj->FP->y)) continue;
    16         Xmin = PS_MIN (Xmin, obj->FP->x);
    17         Xmax = PS_MAX (Xmax, obj->FP->x);
    18         Ymin = PS_MIN (Ymin, obj->FP->y);
    19         Ymax = PS_MAX (Ymax, obj->FP->y);
     13        obj = (pmAstromObj *)input->data[i];
     14        if (!isfinite(obj->FP->x)) continue;
     15        if (!isfinite(obj->FP->y)) continue;
     16        Xmin = PS_MIN (Xmin, obj->FP->x);
     17        Xmax = PS_MAX (Xmax, obj->FP->x);
     18        Ymin = PS_MIN (Ymin, obj->FP->y);
     19        Ymax = PS_MAX (Ymax, obj->FP->y);
    2020    }
    2121
    … …  
    2727    // accumulate 2D histogram in image
    2828    for (int i = 0; i < input->n; i++) {
    29         obj = (pmAstromObj *)input->data[i];
    30         if (!isfinite(obj->FP->x)) continue;
    31         if (!isfinite(obj->FP->y)) continue;
    32         int Xi = PS_MIN (PS_MAX((obj->FP->x - Xmin) / scale + 5, 0), count->numCols);
    33         int Yi = PS_MIN (PS_MAX((obj->FP->y - Ymin) / scale + 5, 0), count->numRows);
    34         count->data.U32[Yi][Xi] ++;
     29        obj = (pmAstromObj *)input->data[i];
     30        if (!isfinite(obj->FP->x)) continue;
     31        if (!isfinite(obj->FP->y)) continue;
     32        int Xi = PS_MIN (PS_MAX((obj->FP->x - Xmin) / scale + 5, 0), count->numCols);
     33        int Yi = PS_MIN (PS_MAX((obj->FP->y - Ymin) / scale + 5, 0), count->numRows);
     34        count->data.U32[Yi][Xi] ++;
    3535    }
    3636
    … …  
    3838    psStats *stats = psStatsAlloc (PS_STAT_MAX | PS_STAT_MAX | PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
    3939    if (!psImageStats(stats, count, NULL, 0)) {
    40         psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
    41         psFree(stats);
    42         psFree(count);
    43         return NULL;
     40        psError(PS_ERR_UNKNOWN, false, "Unable to get image statistics.\n");
     41        psFree(stats);
     42        psFree(count);
     43        return NULL;
    4444    }
    4545
    4646    if (stats->max < 1) {
    47         psError(PS_ERR_UNKNOWN, false, "no valid sources in image\n");
    48         psFree(stats);
    49         psFree(count);
    50         return NULL;
     47        psError(PS_ERR_UNKNOWN, false, "no valid sources in image\n");
     48        psFree(stats);
     49        psFree(count);
     50        return NULL;
    5151    }
    5252
    … …  
    5555    psTrace ("psastro", 4, "skipping stars in cells with more than %f stars\n", limit);
    5656
     57    psastroVisualPlotRemoveClumps (input, count, scale, limit);
     58
    5759    // find and exclude objects in bad pixels
    5860    psArray *output = psArrayAllocEmpty (input->n);
    5961    for (int i = 0; i < input->n; i++) {
    60         obj = (pmAstromObj *)input->data[i];
    61         if (!isfinite(obj->FP->x)) continue;
    62         if (!isfinite(obj->FP->y)) continue;
    63         int Xi = PS_MIN (PS_MAX((obj->FP->x - Xmin) / scale + 5, 0), count->numCols);
    64         int Yi = PS_MIN (PS_MAX((obj->FP->y - Ymin) / scale + 5, 0), count->numRows);
    65         if (count->data.U32[Yi][Xi] > limit) continue;
    66         psArrayAdd (output, 16, obj);
     62        obj = (pmAstromObj *)input->data[i];
     63        if (!isfinite(obj->FP->x)) continue;
     64        if (!isfinite(obj->FP->y)) continue;
     65        int Xi = PS_MIN (PS_MAX((obj->FP->x - Xmin) / scale + 5, 0), count->numCols);
     66        int Yi = PS_MIN (PS_MAX((obj->FP->y - Ymin) / scale + 5, 0), count->numRows);
     67        if (count->data.U32[Yi][Xi] > limit) continue;
     68        psArrayAdd (output, 16, obj);
    6769    }
    6870
    … …  
    7779for (int iy = 0; iy < count->numRows; iy++) {
    7880    for (int ix = 0; ix < count->numCols; ix++) {
    79         if (count->data.U32[iy][ix] > limit) {
    80             psPlane *pixel = psPlaneAlloc();
    81             pixel->x = ix;
    82             pixel->y = iy;
    83             psArrayAdd (badpix, 16, pixel);
    84             psFree (pixel);
    85         }
     81        if (count->data.U32[iy][ix] > limit) {
     82            psPlane *pixel = psPlaneAlloc();
     83            pixel->x = ix;
     84            pixel->y = iy;
     85            psArrayAdd (badpix, 16, pixel);
     86            psFree (pixel);
     87        }
    8688    }
    8789}
Note: See TracChangeset for help on using the changeset viewer.