IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 25433


Ignore:
Timestamp:
Sep 17, 2009, 2:26:32 PM (17 years ago)
Author:
eugene
Message:

change radius for extended sources to use footprint; clean up some of the visualizations; plug some leaks

Location:
branches/eam_branches/20090715/psphot/src
Files:
22 edited

Legend:

Unmodified
Added
Removed

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

    r25409 r25433  
    147147        psphotPetrosianStats.c         \
    148148        psphotPetrosianAnalysis.c      \
    149         pmPetrosian.c
     149        pmPetrosian.c                  \
    150150        psphotEfficiency.c
    151151

  • branches/eam_branches/20090715/psphot/src/pmPetrosian.h

    r25275 r25433  
    3535bool psphotPetrosianFreeVectors(pmPetrosian *petrosian);
    3636
    37 bool psphotPetrosianProfile (pmSource *source, float skynoise);
    38 bool psphotRadialProfilesByAngles (pmPetrosian *petro, pmSource *source, int Nsec, float Rmax);
    39 float psphotRadiusFromProfile (psVector *radius, psVector *flux, float fluxMin, float fluxMax);
    40 bool psphotRadiiFromProfiles (pmPetrosian *petrosian, float fluxMin, float fluxMax);
     37bool psphotPetrosianProfile (pmReadout *readout, pmSource *source, float skynoise);
     38bool psphotRadialProfilesByAngles (pmSource *source, pmPetrosian *petro, int Nsec, float Rmax);
     39float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax);
     40bool psphotRadiiFromProfiles (pmSource *source, pmPetrosian *petrosian, float fluxMin, float fluxMax);
    4141bool psphotEllipticalProfile (pmSource *source, pmPetrosian *petrosian);
    42 bool psphotEllipticalContour (pmPetrosian *petrosian);
     42bool psphotEllipticalContour (pmSource *source, pmPetrosian *petrosian);
    4343bool psphotPetrosianRadialBins (pmSource *source, pmPetrosian *petrosian, float radiusMax, float skynoise);
    44 bool psphotPetrosianStats (pmPetrosian *petrosian);
     44bool psphotPetrosianStats (pmSource *source, pmPetrosian *petrosian);
    4545
    4646bool psphotPetrosianSortPair (psVector *index, psVector *extra);

  • branches/eam_branches/20090715/psphot/src/psphot.h

    r25409 r25433  
    110110bool            psphotInitRadiusEXT (psMetadata *recipe, pmModelType type);
    111111bool            psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal);
     112float           psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal);
    112113
    113114// output functions

  • branches/eam_branches/20090715/psphot/src/psphotBlendFit.c

    r24029 r25433  
    374374    *nfail = Nfail;
    375375
     376    // moments are modified by the fit; re-display
     377    psphotVisualPlotMoments (recipe, sources);
     378    psphotVisualShowResidualImage (readout);
     379
    376380    return true;
    377381}

  • branches/eam_branches/20090715/psphot/src/psphotEfficiency.c

    r25409 r25433  
    134134        psFree(fakeRO);
    135135        psFree(magAll);
     136        psFree(xAll);
     137        psFree(yAll);
    136138        return false;
    137139    }
    138140    psFree(magAll);
     141    psFree(xAll);
     142    psFree(yAll);
    139143
    140144    psBinaryOp(ro->image, ro->image, "+", fakeRO->image);

  • branches/eam_branches/20090715/psphot/src/psphotEllipticalContour.c

    r25178 r25433  
    11# include "psphotInternal.h"
    2 
    3 // XXX consistency : theta in radians here and in calling functions
    42
    53// model parameters
    … …  
    75psF32 psphotEllipticalContourFunc (psVector *deriv, const psVector *params, const psVector *coord);
    86
    9 bool psphotEllipticalContour (pmPetrosian *petrosian) {
     7bool psphotEllipticalContour (pmSource *source, pmPetrosian *petrosian) {
    108
    119    // use LMM to fit theta vs radius to an ellipse
    … …  
    1917    // arrays to hold the data to be fitted
    2018    // we fit x and y vs theta in separate passes.
    21     psArray *x = psArrayAlloc(2*radius->n);
    22     psVector *y = psVectorAlloc(2*radius->n, PS_TYPE_F32);
    23     psVector *yErr = psVectorAlloc(2*radius->n, PS_TYPE_F32);
     19    psArray *x = psArrayAllocEmpty(2*radius->n);
     20    psVector *y = psVectorAllocEmpty(2*radius->n, PS_TYPE_F32);
     21    psVector *yErr = psVectorAllocEmpty(2*radius->n, PS_TYPE_F32);
    2422
    2523    int n = 0;
    2624    for (int i = 0; i < radius->n; i++) {
     25        if (!isfinite(radius->data.F32[i])) continue;
    2726
    2827        psVector *coord = NULL;
    … …  
    5049        Rmax = MAX (Rmax, radius->data.F32[i]);
    5150    }   
    52     assert (x->n == n);
    53     assert (y->n == n);
     51    x->n = n;
     52    y->n = n;
     53    yErr->n = n;
     54
     55    if (n < 4) {
     56        psFree (x);
     57        psFree (y);
     58        psFree (yErr);
     59        return false;
     60    }
    5461
    5562    psVector *params = psVectorAlloc (3, PS_TYPE_F32);
    … …  
    6168
    6269    // XXX for now, no parameter masks, skip checkLimits
     70    // XXX might help to add a limit to the angle or re-parameterize the ellipse in terms of Rxx, Ryy, Rxy
    6371    // constraint->checkLimits = psastroModelBoresiteLimits;
    6472
    … …  
    8492    }
    8593
    86     fprintf (stderr, "# fitted values:\n");
    87     fprintf (stderr, "Po:  %f\n", params->data.F32[PAR_PHI]*PS_DEG_RAD);
    88     fprintf (stderr, "Ep:  %f\n", params->data.F32[PAR_EPSILON]);
    89     fprintf (stderr, "Rm:  %f\n", params->data.F32[PAR_RMIN]);
    90 
     94    psTrace ("psphot", 4, "# fitted values:\n");
     95    psTrace ("psphot", 4, "Phi:   %f\n", petrosian->axes.theta*PS_DEG_RAD);
     96    psTrace ("psphot", 4, "Rmaj:  %f\n", petrosian->axes.major);
     97    psTrace ("psphot", 4, "Rmin:  %f\n", petrosian->axes.minor);
     98   
    9199    // show the results
    92     psphotPetrosianVisualEllipticalContour (petrosian);
     100    // psphotPetrosianVisualEllipticalContour (petrosian);
    93101
    94102    psFree (x);
    95103    psFree (y);
    96104    psFree (yErr);
     105    psFree (params);
     106    psFree (covar);
     107    psFree (myMin);
     108    psFree (constraint);
    97109
    98110    return true;

  • branches/eam_branches/20090715/psphot/src/psphotEllipticalProfile.c

    r25363 r25433  
    3232    float Syy = shape.sy;
    3333
     34    // XXX drop these two vectors?
    3435    psVector *radiusRaw = psVectorAllocEmpty(100, PS_TYPE_F32);
    3536    psVector *fluxRaw = psVectorAllocEmpty(100, PS_TYPE_F32);
    … …  
    6566
    6667    // psphotPetrosianVisualProfileRadii (radius, flux, radiusRaw, fluxRaw, 0.0);
    67     psphotPetrosianVisualProfileByAngle (radius, flux);
     68    // psphotPetrosianVisualProfileByAngle (radius, flux);
     69
     70    psFree (radiusRaw);
     71    psFree (fluxRaw);
    6872    return true;
    6973}

  • branches/eam_branches/20090715/psphot/src/psphotFitSourcesLinear.c

    r25409 r25433  
    239239
    240240    psphotVisualShowResidualImage (readout);
    241     psphotVisualShowFlags (sources);
     241    // psphotVisualShowFlags (sources);
    242242
    243243    return true;

  • branches/eam_branches/20090715/psphot/src/psphotImageQuality.c

    r23989 r25433  
    279279#endif
    280280
    281     psLogMsg ("psphot", PS_LOG_INFO, "Image Quality Stats from %ld psf stars : FWHM (major, minor) : %f, %f\n",
     281    psLogMsg ("psphot", PS_LOG_INFO, "Image Quality Stats from %ld psf stars : FWHM (major, minor) [pixels]: %f, %f\n",
    282282              M2->n, fwhm_major, fwhm_minor);
    283283
    284     psLogMsg ("psphot", PS_LOG_INFO, "M_2 : %f +/- %f, M_3 : %f +/- %f, M_4 : %f +/- %f\n",
     284    psLogMsg ("psphot", PS_LOG_INFO, "M_2 : %f +/- %f, M_3 : %f +/- %f, M_4 : %f +/- %f  [pixels^n]\n",
    285285              vM2, dM2, vM3, dM3, vM4, dM4);
    286286

  • branches/eam_branches/20090715/psphot/src/psphotPSFConvModel.c

    r21183 r25433  
    6868    psVector *dparams = modelConv->dparams;
    6969
     70    // XXX this needs to be changed to use the new psphotSetRadiusEXT function
    7071    psphotCheckRadiusEXT (readout, source, modelConv, markVal);
    7172

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianAnalysis.c

    r25275 r25433  
    1212    // XXX temporary user-supplied systematic sky noise measurement (derive from background model)
    1313    float skynoise = psMetadataLookupF32 (&status, recipe, "SKY.NOISE");
     14
     15# if (0)
     16    // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
     17    // XXX use this to set skynoise
     18    pmReadout *backModel = psphotSelectBackground (config, view);
     19    pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);
     20# endif
    1421
    1522    // S/N limit to perform full non-linear fits
    … …  
    3340        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
    3441        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
     42        if (source->mode & PM_SOURCE_MODE_DEFECT) continue;
     43        if (source->mode & PM_SOURCE_MODE_SATSTAR) continue;
     44        if (!(source->mode & PM_SOURCE_MODE_EXT_LIMIT)) continue;
    3545
    3646        // limit selection to some SN limit
    3747        assert (source->peak); // how can a source not have a peak?
    3848        if (source->peak->SN < SN_LIM) continue;
    39         if (source->extNsigma < 10.0) continue; // XXX this should not be hardwired
    4049
    4150        // limit selection by analysis region
    … …  
    5059        }
    5160
    52         psphotPetrosianProfile (source, skynoise);
     61        psphotPetrosianProfile (readout, source, skynoise);
     62
     63        pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     64        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    5365    }
     66
     67    psphotVisualShowResidualImage (readout);
    5468
    5569    // pause and wait for user input:

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianProfile.c

    r25275 r25433  
    77// structure to something the pmRadialProfile
    88
    9 bool psphotPetrosianProfile (pmSource *source, float skynoise) {
     9bool psphotPetrosianProfile (pmReadout *readout, pmSource *source, float skynoise) {
    1010
    11   // container to hold results from the radial profile analysis
    12   pmPetrosian *petrosian = pmPetrosianAlloc();
     11    // container to hold results from the radial profile analysis
     12    pmPetrosian *petrosian = pmPetrosianAlloc();
    1313
    14   int Nsec = 24;
    15   float Rmax = 200;
    16   float fluxMin = 0.0;
    17   float fluxMax = source->peak->flux;
     14    int Nsec = 24;
     15    float Rmax = 200;
     16    float fluxMin = 0.0;
     17    float fluxMax = source->peak->flux;
    1818
    19   // generate a series of radial profiles at Nsec evenly spaced angles.  the profile flux
    20   // is measured by interpolation for small radii; for large radii, the pixels in a box
    21   // are averaged to increase the S/N (XXX not yet done)
    22   if (!psphotRadialProfilesByAngles (petrosian, source, Nsec, Rmax)) {
    23     psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian");
    24     return false;
    25   }
     19    // generate a series of radial profiles at Nsec evenly spaced angles.  the profile flux
     20    // is measured by interpolation for small radii; for large radii, the pixels in a box
     21    // are averaged to increase the S/N (XXX not yet done)
     22    if (!psphotRadialProfilesByAngles (source, petrosian, Nsec, Rmax)) {
     23        psError (PS_ERR_UNKNOWN, false, "failed to measure radial profile for petrosian");
     24        psFree (petrosian);
     25        return false;
     26    }
    2627
    27   // use the radial profiles to determine the radius of a given isophote.  this isophote
    28   // is used to determine the elliptical shape of the object, so it has a relatively high
    29   // value (nominally 50% of the peak)
    30   if (!psphotRadiiFromProfiles (petrosian, fluxMin, fluxMax)) {
    31     psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles");
    32     return false;
    33   }
     28    // use the radial profiles to determine the radius of a given isophote.  this isophote
     29    // is used to determine the elliptical shape of the object, so it has a relatively high
     30    // value (nominally 50% of the peak)
     31    if (!psphotRadiiFromProfiles (source, petrosian, fluxMin, fluxMax)) {
     32        psError (PS_ERR_UNKNOWN, false, "failed to measure isophotal radii from profiles");
     33        psFree (petrosian);
     34        return false;
     35    }
    3436
    35   // convert the isophotal radius vs angle measurements to an elliptical contour
    36   if (!psphotEllipticalContour (petrosian)) {
    37     psError (PS_ERR_UNKNOWN, false, "failed to measure elliptical contour");
    38     return false;
    39   }
     37    // convert the isophotal radius vs angle measurements to an elliptical contour
     38    if (!psphotEllipticalContour (source, petrosian)) {
     39        psLogMsg ("psphot", 3, "failed to measure elliptical contour");
     40        psFree (petrosian);
     41        return false;
     42    }
    4043 
    41   // generate a single, normalized radial profile following the elliptical contours.
    42   // the radius is normalized by the axis ratio so that on the major axis, 1 pixel = 1 pixel
    43   if (!psphotEllipticalProfile (source, petrosian)) {
    44     psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
    45     return false;
    46   }
     44    // generate a single, normalized radial profile following the elliptical contours.
     45    // the radius is normalized by the axis ratio so that on the major axis, 1 pixel = 1 pixel
     46    if (!psphotEllipticalProfile (source, petrosian)) {
     47        psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
     48        psFree (petrosian);
     49        return false;
     50    }
    4751 
    48   // integrate the radial profile for radial bins defined for the petrosian measurement:
    49   // SB_i (r_i) where \alpha r_i < r < \beta r_i
    50   if (!psphotPetrosianRadialBins (source, petrosian, Rmax, skynoise)) {
    51     psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
    52     return false;
    53   }
     52    // integrate the radial profile for radial bins defined for the petrosian measurement:
     53    // SB_i (r_i) where \alpha r_i < r < \beta r_i
     54    if (!psphotPetrosianRadialBins (source, petrosian, Rmax, skynoise)) {
     55        psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
     56        psFree (petrosian);
     57        return false;
     58    }
    5459 
    55   // use the SB_i from above to calculate the petrosian radius and the flux within that radius
    56   if (!psphotPetrosianStats (petrosian)) {
    57     psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
    58     return false;
    59   }
     60    // use the SB_i from above to calculate the petrosian radius and the flux within that radius
     61    if (!psphotPetrosianStats (source, petrosian)) {
     62        psError (PS_ERR_UNKNOWN, false, "failed to generate elliptical profile");
     63        psFree (petrosian);
     64        return false;
     65    }
    6066 
    61   // XXX this will only work in the psphot context, not the psphotPetrosianStudy...
    62   // XXX add the petrosian to the pmSource structure...
    63   psphotVisualShowPetrosian (source, petrosian);
     67    // XXX this will only work in the psphot context, not the psphotPetrosianStudy...
     68    // XXX add the petrosian to the pmSource structure...
     69    // psphotVisualShowResidualImage (readout);
     70    psphotVisualShowPetrosian (source, petrosian);
    6471
    65   psphotPetrosianFreeVectors(petrosian);
     72    psphotPetrosianFreeVectors(petrosian);
    6673
    67   fprintf (stdout, "\n petrosian radius: %f\n flux: %f\n axis ratio: %f\n angle: %f\n",
    68            petrosian->petrosianRadius, petrosian->petrosianFlux, petrosian->axes.minor/petrosian->axes.major, PS_DEG_RAD*petrosian->axes.theta);
     74    psTrace ("psphot", 3, "source at %f,%f: petrosian radius: %f, flux: %f, axis ratio: %f, angle: %f",
     75             source->peak->xf, source->peak->yf, petrosian->petrosianRadius, petrosian->petrosianFlux, petrosian->axes.minor/petrosian->axes.major, PS_DEG_RAD*petrosian->axes.theta);
    6976
    70   return true;
     77    psFree (petrosian);
     78    return true;
    7179}

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianRadialBins.c

    r25275 r25433  
    1818   
    1919    float skyModelErrorSQ = PS_SQR(skynoise);
    20 
    21 # if (0)
    22     // if backModel or backStdev are missing, the values of sky and/or skyErr will be set to NAN
    23     pmReadout *backModel = psphotSelectBackground (config, view);
    24     pmReadout *backStdev = psphotSelectBackgroundStdev (config, view);
    25 # endif
    2620
    2721    psVector *radius = petrosian->radiusElliptical;
    … …  
    10498        binArea->data.F32[i] = M_PI * (rMax2 - rMin2);
    10599
    106         if (0) {
    107             fprintf (stderr, "%3d  %5.1f %5.1f : %5.1f : %5.1f %5.1f\n",
     100        psTrace ("psphot", 6, "%3d  %5.1f %5.1f : %5.1f : %5.1f %5.1f\n",
    108101                 i, radAlp->data.F32[i], radMin->data.F32[i], binRad->data.F32[i],
    109102                 radMax->data.F32[i], radBet->data.F32[i]);
    110         }
    111103    }
    112104
    … …  
    131123        if (radius->data.F32[i] > Rmax) {
    132124            // calculate the value for the nOut bin
    133             psVectorStats (stats, values, NULL, NULL, 0);
     125            float value, dvalue;
     126            if (values->n > 0) {
     127                psVectorStats (stats, values, NULL, NULL, 0);
     128                value = stats->robustMedian;
     129                dvalue = stats->robustStdev;
     130            } else {
     131                value = NAN;
     132                dvalue = NAN;
     133            }
    134134            // binSB->data.F32[nOut] = stats->sampleMedian;
    135             binSB->data.F32[nOut] = stats->robustMedian;
    136             binSBstdev->data.F32[nOut] = sqrt(PS_SQR(stats->robustStdev) / values->n + skyModelErrorSQ);
     135            binSB->data.F32[nOut] = value;
     136            binSBstdev->data.F32[nOut] = sqrt(PS_SQR(dvalue) / values->n + skyModelErrorSQ);
    137137            // binSB->data.F32[nOut] = stats->fittedMean;
    138138            // binSBstdev->data.F32[nOut] = sqrt(PS_SQR(stats->fittedStdev) / values->n + skyModelErrorSQ);
    … …  
    142142            // residual flux, but the sky from the sky model)
    143143
    144             if (1) {
    145                 fprintf (stderr, "%3d  %5.1f %5.1f : %5.1f  %5.2f\n",
    146                          nOut, radAlp->data.F32[nOut], radBet->data.F32[nOut], binSB->data.F32[nOut], binSBstdev->data.F32[nOut]);
    147             }
     144            psTrace ("psphot", 5, "%3d  %5.1f %5.1f : %5.1f  %5.2f\n",
     145                     nOut, radAlp->data.F32[nOut], radBet->data.F32[nOut], binSB->data.F32[nOut], binSBstdev->data.F32[nOut]);
    148146
    149147            nOut ++;
    … …  
    170168    petrosian->binSBstdev = binSBstdev;
    171169
    172     psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, source->peak->flux, 0.0);
     170    // psphotPetrosianVisualProfileRadii (radius, flux, binRad, binSB, source->peak->flux, 0.0);
    173171
    174172    psFree(radMin);

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianStats.c

    r25275 r25433  
    88float InterpolateValues (float X0, float Y0, float X1, float Y1, float X);
    99
    10 bool psphotPetrosianStats (pmPetrosian *petrosian) {
     10bool psphotPetrosianStats (pmSource *source, pmPetrosian *petrosian) {
    1111
    1212    float petRadius, petFlux;
    … …  
    7777        psVectorAppend(refRadius, binRad->data.F32[i]);
    7878
    79         if (1) {
    80             fprintf (stderr, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f\n",
    81                      i, refRadius->data.F32[nOut],
    82                      binSB->data.F32[i], binSBstdev->data.F32[i],
    83                      meanSB->data.F32[nOut], meanSBerr,
    84                      petRatio->data.F32[nOut], petRatioErr->data.F32[nOut],
    85                      fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], areaInner);
    86         }
     79        psTrace ("psphot", 4, "%3d : %5.2f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.3f %5.3f : %5.1f %5.1f\n",
     80                 i, refRadius->data.F32[nOut],
     81                 binSB->data.F32[i], binSBstdev->data.F32[i],
     82                 meanSB->data.F32[nOut], meanSBerr,
     83                 petRatio->data.F32[nOut], petRatioErr->data.F32[nOut],
     84                 fluxSum->data.F32[nOut], sqrt(fluxSumErr2->data.F32[nOut]), areaSum->data.F32[nOut], areaInner);
    8785   
    8886        // anytime we transition below the PETROSIAN_RATIO, calculate the radius and flux
    … …  
    149147
    150148    psFree(fluxSum);
     149    psFree(fluxSumErr2);
     150    psFree(refRadius);
    151151    psFree(petRatio);
     152    psFree(petRatioErr);
    152153    psFree(meanSB);
    153154    psFree(areaSum);

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianStudy.c

    r25275 r25433  
    136136  }
    137137
    138   psphotPetrosianProfile (source, skynoise);
     138  psphotPetrosianProfile (readout, source, skynoise);
    139139
    140140  psFree (source);

  • branches/eam_branches/20090715/psphot/src/psphotPetrosianVisual.c

    r25275 r25433  
    1919
    2020static int kapa = -1;
     21static int kapa2 = -1;
    2122
    2223// if no valid data is supplied (NULL or n <- 0), leave limits as they were
    … …  
    5354
    5455    // return true;
    55     if (!pmVisualIsVisual()) return true;
    56 
    57     if (kapa == -1) {
    58         kapa = KapaOpenNamedSocket ("kapa", "psphot:plots");
    59         if (kapa == -1) {
     56    // if (!pmVisualIsVisual()) return true;
     57
     58    if (kapa2 == -1) {
     59        kapa2 = KapaOpenNamedSocket ("kapa", "psphot:plots");
     60        if (kapa2 == -1) {
    6061            fprintf (stderr, "failure to open kapa; visual mode disabled\n");
    6162            pmVisualSetVisual(false);
    … …  
    6465    }
    6566
    66     KapaClearPlots (kapa);
     67    KapaClearPlots (kapa2);
    6768    KapaInitGraph (&graphdata);
    68     KapaSetFont (kapa, "courier", 14);
     69    KapaSetFont (kapa2, "courier", 14);
    6970
    7071    pmVisualLimitsFromVectors (&graphdata, radius, flux);
    71     KapaSetLimits (kapa, &graphdata);
    72 
    73     KapaBox (kapa, &graphdata);
    74     KapaSendLabel (kapa, "radius", KAPA_LABEL_XM);
    75     KapaSendLabel (kapa, "flux", KAPA_LABEL_YM);
    76 
    77     graphdata.color = KapaColorByName ("black");
    78     graphdata.style = 2;
    79     graphdata.ptype = 0;
    80     graphdata.size = 1.0;
    81     KapaPrepPlot (kapa, radius->n, &graphdata);
    82     KapaPlotVector (kapa, radius->n, radius->data.F32, "x");
    83     KapaPlotVector (kapa, radius->n, flux->data.F32, "y");
     72    KapaSetLimits (kapa2, &graphdata);
     73
     74    KapaBox (kapa2, &graphdata);
     75    KapaSendLabel (kapa2, "radius", KAPA_LABEL_XM);
     76    KapaSendLabel (kapa2, "flux", KAPA_LABEL_YM);
     77
     78    graphdata.color = KapaColorByName ("black");
     79    graphdata.style = 2;
     80    graphdata.ptype = 0;
     81    graphdata.size = 1.0;
     82    KapaPrepPlot (kapa2, radius->n, &graphdata);
     83    KapaPlotVector (kapa2, radius->n, radius->data.F32, "x");
     84    KapaPlotVector (kapa2, radius->n, flux->data.F32, "y");
    8485
    8586    // pause and wait for user input:

  • branches/eam_branches/20090715/psphot/src/psphotRadialProfileByAngle.c

    r25361 r25433  
    11# include "psphotInternal.h"
    22
    3 // Given a source at (x,y), generate a collection of radial profiles at even angular
    4 // separations
    5 
     3// Given a source at (x,y), generate a collection of radial profiles at even angular separations
     4
     5// These functions are used to calculate the stats in a rectangle at arbitrary orientation.
     6// XXX Move these elsewhere (psLib?)
    67float psphotMeanSectorValue (psImage *image, float x, float y, float dL, float dW, float theta);
    78psVector *psphotBoxValues (psImage *image, float x0, float y0, float dL, float dW, float theta);
    … …  
    910psVector *psphotLineValuesBresen (psImage *image, int X1, int Y1, int X2, int Y2, int dW, int swapcoords);
    1011
    11 bool psphotRadialProfilesByAngles (pmPetrosian *petrosian, pmSource *source, int Nsec, float Rmax) {
     12bool psphotRadialProfilesByAngles (pmSource *source, pmPetrosian *petrosian, int Nsec, float Rmax) {
    1213
    1314    // we want to have an even number of sectors so we can do 180 deg symmetrizing
    … …  
    124125
    125126    psVector *values = psphotBoxValues (image, x, y, dL, dW, theta);
     127    if (!values) goto escape;
     128    if (!values->n) goto escape;
    126129   
    127130    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
    … …  
    134137   
    135138    return value;
     139
     140escape:
     141    psFree(values);
     142    return NAN;   
    136143}
    137144

  • branches/eam_branches/20090715/psphot/src/psphotRadiiFromProfiles.c

    r25275 r25433  
    33// Given the Petrosian data (radii, fluxes) determine the radius for each profile at the desisred isophote
    44
    5 bool psphotRadiiFromProfiles (pmPetrosian *petrosian, float fluxMin, float fluxMax) {
     5bool psphotRadiiFromProfiles (pmSource *source, pmPetrosian *petrosian, float fluxMin, float fluxMax) {
    66
    77  petrosian->isophotalRadii = psVectorAlloc(petrosian->theta->n, PS_TYPE_F32);
    … …  
    1010      psVector *radii = petrosian->radii->data[i];
    1111      psVector *fluxes = petrosian->fluxes->data[i];
    12       float radius = psphotRadiusFromProfile (radii, fluxes, fluxMin, fluxMax);
     12      float radius = psphotRadiusFromProfile (source, radii, fluxes, fluxMin, fluxMax);
    1313
    1414      // psphotPetrosianVisualProfileByAngle (radii, fluxes, radius);
    … …  
    2020}
    2121
    22 float psphotRadiusFromProfile (psVector *radius, psVector *flux, float fluxMin, float fluxMax) {
     22float psphotRadiusFromProfile (pmSource *source, psVector *radius, psVector *flux, float fluxMin, float fluxMax) {
    2323
    2424    // 'flux' is a noisy sample of the galaxy radial profile at points 'radius'
    … …  
    4949            psVectorAppend (values, radius->data.F32[i]);
    5050        }
    51         psVectorStats (stats, values, NULL, NULL, 0);
     51        if (values->n > 1) {
     52            psVectorStats (stats, values, NULL, NULL, 0);
    5253
    53         // if we have a valid range, rebin with bin size 1/2 of median radius
    54         if (isfinite(stats->sampleMedian)) {
    55             Rbin = MAX(1, 0.5*stats->sampleMedian);
     54            // if we have a valid range, rebin with bin size 1/2 of median radius
     55            if (isfinite(stats->sampleMedian)) {
     56                Rbin = MAX(1, 0.5*stats->sampleMedian);
     57            }
    5658        }
     59        psFree (values);
     60        psFree (stats);
    5761    }
    5862    Rbin = 3;
    … …  
    9094                // calculate the value for the nOut bin
    9195                // XXX need to fix this as well psStats (stats, values);
    92                 psVectorStats (stats, values, NULL, NULL, 0);
    93                 fluxBinned->data.F32[nOut] = stats->sampleMedian;
     96                float value;
     97                if (values->n > 0) {
     98                    psVectorStats (stats, values, NULL, NULL, 0);
     99                    value = stats->sampleMedian;
     100                } else {
     101                    value = NAN;
     102                }
     103                fluxBinned->data.F32[nOut] = value;
    94104                nOut ++;
    95105                radiusBinned->data.F32[nOut] = (nOut + 0.5)*Rbin;
    … …  
    118128            // XXX is there a macro in psLib that does this interpolation?
    119129            if (i == 0) {
    120                 psWarning ("impossible condition f[0] < Fo");
    121                 continue;
    122             } else {
    123                 Ro = radiusBinned->data.F32[i-1] + (radiusBinned->data.F32[i] - radiusBinned->data.F32[i-1]) * (Fo - fluxBinned->data.F32[i-1]) / (fluxBinned->data.F32[i] - fluxBinned->data.F32[i-1]);
    124             }
     130                psLogMsg ("psphot", 3, "bogus radial profile for ..., skipping");
     131                psFree (fluxBinned);
     132                psFree (radiusBinned);
     133                return NAN;
     134            }
     135            Ro = radiusBinned->data.F32[i-1] + (radiusBinned->data.F32[i] - radiusBinned->data.F32[i-1]) * (Fo - fluxBinned->data.F32[i-1]) / (fluxBinned->data.F32[i] - fluxBinned->data.F32[i-1]);
    125136            above = FALSE;
    126137        }

  • branches/eam_branches/20090715/psphot/src/psphotRadiusChecks.c

    r21183 r25433  
    8686
    8787// call this function whenever you (re)-define the EXT model
     88float psphotSetRadiusEXT (pmReadout *readout, pmSource *source, psImageMaskType markVal) {
     89
     90    psAssert (source, "source not defined??");
     91    psAssert (source->peak, "peak not defined??");
     92
     93    pmPeak *peak = source->peak;
     94
     95    // set the radius based on the footprint:
     96    if (!peak->footprint) goto escape;
     97    pmFootprint *footprint = peak->footprint;
     98    if (!footprint->spans) goto escape;
     99    if (footprint->spans->n < 1) goto escape;
     100
     101    // find the max radius
     102    float radius = 0.0;
     103    for (int j = 0; j < footprint->spans->n; j++) {
     104        pmSpan *span = footprint->spans->data[j];
     105
     106        float dY  = span->y  - peak->yf;
     107        float dX0 = span->x0 - peak->xf;
     108        float dX1 = span->x1 - peak->xf;
     109
     110        radius = PS_MAX (radius, hypot(dY, dX0));
     111        radius = PS_MAX (radius, hypot(dY, dX1));
     112    }
     113
     114    radius += EXT_FIT_PADDING;
     115    if (isnan(radius)) psAbort("error in radius");
     116
     117    // redefine the pixels if needed
     118    pmSourceRedefinePixels (source, readout, peak->xf, peak->yf, radius);
     119
     120    // set the mask to flag the excluded pixels
     121    psImageKeepCircle (source->maskObj, peak->xf, peak->yf, radius, "OR", markVal);
     122    return radius;
     123
     124escape:
     125    return NAN;
     126    // bool result = psphotCheckRadiusEXT_old (readout, source, model, markVal);
     127    // return result;
     128}
     129
     130// call this function whenever you (re)-define the EXT model
    88131bool psphotCheckRadiusEXT (pmReadout *readout, pmSource *source, pmModel *model, psImageMaskType markVal) {
     132
     133    psAbort ("do not use this function");
    89134
    90135    psF32 *PAR = model->params->data.F32;

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

    r24592 r25433  
    223223    if (source->peak->SN < EXT_MIN_SN) return false;
    224224
     225    float radius = psphotSetRadiusEXT (readout, source, markVal);
     226
     227    // XXX note that this changes the source moments that are published...
    225228    // recalculate the source moments using the larger extended-source moments radius
    226229    // at this stage, skip Gaussian windowing, and do not clip pixels by S/N
    227     if (!pmSourceMoments (source, EXT_MOMENTS_RAD, 0.0, 0.0)) return false;
     230    // XXX use the deteremined value of PSF_MOMENTS_RADIUS and SIGMA and scale up?
     231    // XXX increase the moments radius iteratively until the measured sigma / input < f?
     232    // XXX use the footprint somehow to judge both radius and aperture?
     233    if (!pmSourceMoments (source, radius, 0.0, 0.0)) return false;
    228234
    229235    psTrace ("psphot", 5, "trying blob...\n");
    … …  
    277283
    278284    // both models failed; reject them both
     285    // XXX -- change type flags to psf in this case and keep original moments?
    279286    psFree (EXT);
    280287    psFree (DBL);
    … …  
    287294    // save new model
    288295    source->modelEXT = EXT;
     296    source->modelEXT->radiusFit = radius;
    289297    source->type = PM_SOURCE_TYPE_EXTENDED;
    290298    source->mode |= PM_SOURCE_MODE_EXTMODEL;
    … …  
    294302    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
    295303    source->tmpFlags |=  PM_SOURCE_TMPF_SUBTRACTED;
     304
     305# if (PS_TRACE_ON)   
    296306    psTrace ("psphot", 5, "blob as EXT: %f %f\n", EXT->params->data.F32[PM_PAR_XPOS], EXT->params->data.F32[PM_PAR_YPOS]);
     307    if (psTraceGetLevel("psphot") >= 6) {
     308        psLogMsg ("psphot", 1, "source 2:\n");
     309        for (int i = 0; i < source->modelEXT->params->n; i++) {
     310            psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelEXT->params->data.F32[i], source->modelEXT->dparams->data.F32[i]);
     311        }
     312    }
     313# endif
     314
    297315    return true;
    298316
    … …  
    306324    source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
    307325    source->mode     |= PM_SOURCE_MODE_PAIR;
     326    source->modelPSF->radiusFit = radius;
    308327
    309328    // copy most data from the primary source (modelEXT, blends stay NULL)
    310     // XXX use pmSourceCopy?
    311329    pmSource *newSrc = pmSourceCopy (source);
    312330    newSrc->modelPSF = psMemIncrRefCounter (DBL->data[1]);
     331    newSrc->modelPSF->radiusFit = radius;
    313332
    314333    // build cached models and subtract
    … …  
    317336    pmSourceCacheModel (newSrc, maskVal);
    318337    pmSourceSub (newSrc, PM_MODEL_OP_FULL, maskVal);
     338
     339# if (PS_TRACE_ON)   
    319340    psTrace ("psphot", 5, "blob as DBL: %f %f\n", ONE->params->data.F32[PM_PAR_XPOS], ONE->params->data.F32[PM_PAR_YPOS]);
     341    if (psTraceGetLevel("psphot") >= 6) {
     342        psLogMsg ("psphot", 1, "source 1:\n");
     343        for (int i = 0; i < newSrc->modelPSF->params->n; i++) {
     344            psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, newSrc->modelPSF->params->data.F32[i], newSrc->modelPSF->dparams->data.F32[i]);
     345        }
     346        psLogMsg ("psphot", 1, "source 2:\n");
     347        for (int i = 0; i < source->modelPSF->params->n; i++) {
     348            psLogMsg ("psphot", 1, "PAR %d : %f +/- %f\n", i, source->modelPSF->params->data.F32[i], source->modelPSF->dparams->data.F32[i]);
     349        }
     350        psphotVisualShowResidualImage (readout);
     351    }
     352# endif
    320353
    321354    psArrayAdd (newSources, 100, newSrc);
    … …  
    356389    // save the PSF model from the Ensemble fit
    357390    PSF = source->modelPSF;
    358     psphotCheckRadiusPSFBlend (readout, source, PSF, markVal, 8.0);
     391    // psphotCheckRadiusPSFBlend (readout, source, PSF, markVal, 8.0);
    359392    if (isnan(PSF->params->data.F32[1])) psAbort("nan in dbl fit");
    360393
    … …  
    390423
    391424    // if (isnan(EXT->params->data.F32[1])) psAbort("nan in ext fit");
    392 
    393     psphotCheckRadiusEXT (readout, source, EXT, markVal);
     425    // psphotCheckRadiusEXT (readout, source, EXT, markVal);
    394426
    395427    if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
    … …  
    401433    return (EXT);
    402434}
    403 

  • branches/eam_branches/20090715/psphot/src/psphotSourceSize.c

    r25390 r25433  
    194194bool psphotSourceSizePSF (psphotSourceSizeOptions *options, psArray *sources, pmPSF *psf) {
    195195
    196     // select the psf stars
    197     psArray *psfstars = psArrayAllocEmpty (100);
    198 
     196    // select stats from the psf stars
    199197    psVector *Ap = psVectorAllocEmpty (100, PS_TYPE_F32);
    200198    psVector *ApErr = psVectorAllocEmpty (100, PS_TYPE_F32);
    … …  
    205203        pmSource *source = sources->data[i];
    206204        if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
    207         psArrayAdd (psfstars, 100, source);
    208205
    209206        // XXX can we test if psfMag is set and calculate only if needed?
    … …  
    229226    options->ApSysErr = psVectorSystematicError(dAp, ApErr, 0.05);
    230227    fprintf (stderr, "psf-sum: %f +/- %f\n", options->ApResid, options->ApSysErr);
     228
     229    psFree (Ap);
     230    psFree (ApErr);
     231    psFree (stats);
     232    psFree (dAp);
    231233
    232234    return true;

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

    r25388 r25433  
    166166bool psphotVisualShowImage (pmReadout *readout) {
    167167
    168     if (!pmVisualIsVisual()) return true;
     168    // if (!pmVisualIsVisual()) return true;
    169169
    170170    int kapa = psphotKapaChannel (1);
    … …  
    20002000bool psphotVisualShowResidualImage (pmReadout *readout) {
    20012001
    2002     if (!pmVisualIsVisual()) return true;
     2002    // if (!pmVisualIsVisual()) return true;
    20032003
    20042004    int myKapa = psphotKapaChannel (1);
    … …  
    21002100    if (source == NULL) return true;
    21012101
    2102     if (!pmVisualIsVisual()) return true;
    2103 
    2104     int kapa = psphotKapaChannel (3);
     2102    // if (!pmVisualIsVisual()) return true;
     2103
     2104    int kapa = psphotKapaChannel (1);
    21052105    if (kapa == -1) return false;
    21062106
Note: See TracChangeset for help on using the changeset viewer.