IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 9734


Ignore:
Timestamp:
Oct 24, 2006, 2:07:03 PM (20 years ago)
Author:
Paul Price
Message:

Changed definition of psVectorAlloc and psArrayAlloc

Location:
trunk/psphot/src
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotApResid.c

    r9568 r9734  
    44
    55// measure the aperture residual statistics
    6 bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) { 
     6bool psphotApResid (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf) {
    77
    88    int Nfail = 0;
    … …  
    1818    float MAX_AP_OFFSET = psMetadataLookupF32 (&status, recipe, "MAX_AP_OFFSET");
    1919
    20     // set limits on the aperture magnitudes 
     20    // set limits on the aperture magnitudes
    2121    pmSourceMagnitudesInit (recipe);
    2222
    … …  
    2525    psf->growth = pmGrowthCurveAlloc (3.0, 100.0, REF_RADIUS);
    2626    psphotGrowthCurve (readout, psf);
    27    
    28     psVector *mask    = psVectorAlloc (300, PS_TYPE_U8);
    29     psVector *xPos    = psVectorAlloc (300, PS_TYPE_F64);
    30     psVector *yPos    = psVectorAlloc (300, PS_TYPE_F64);
    31     psVector *flux    = psVectorAlloc (300, PS_TYPE_F64);
    32     psVector *r2rflux = psVectorAlloc (300, PS_TYPE_F64);
    33     psVector *apResid = psVectorAlloc (300, PS_TYPE_F64);
    34     psVector *dMag    = psVectorAlloc (300, PS_TYPE_F64);
    35     mask->n = xPos->n = yPos->n = flux->n = r2rflux->n = apResid->n = dMag->n = 0;
     27
     28    psVector *mask    = psVectorAllocEmpty (300, PS_TYPE_U8);
     29    psVector *xPos    = psVectorAllocEmpty (300, PS_TYPE_F64);
     30    psVector *yPos    = psVectorAllocEmpty (300, PS_TYPE_F64);
     31    psVector *flux    = psVectorAllocEmpty (300, PS_TYPE_F64);
     32    psVector *r2rflux = psVectorAllocEmpty (300, PS_TYPE_F64);
     33    psVector *apResid = psVectorAllocEmpty (300, PS_TYPE_F64);
     34    psVector *dMag    = psVectorAllocEmpty (300, PS_TYPE_F64);
    3635    Npsf = 0;
    3736
    3837    // select all good PM_SOURCE_TYPE_STAR entries
    3938    for (int i = 0; i < sources->n; i++) {
    40         source = sources->data[i];
    41         model = source->modelPSF;
    42 
    43         if (source->type != PM_SOURCE_TYPE_STAR) continue;
    44         if (source->mode &  PM_SOURCE_MODE_SATSTAR) continue;
    45         if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
    46         if (source->mode &  PM_SOURCE_MODE_FAIL) continue;
    47         if (source->mode &  PM_SOURCE_MODE_POOR) continue;
    48 
    49         // get growth-corrected, apTrend-uncorrected magnitudes in scaled apertures
    50         // will fail if below S/N threshold or model is missing
    51         if (!pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH)) {
    52             Nskip ++;
    53             continue;
    54         }
    55 
    56         apResid->data.F64[Npsf] = source->apMag - source->psfMag;
    57 
    58         xPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
    59         yPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_YPOS];
    60 
    61         flux->data.F64[Npsf]    = pow(10.0, -0.4*source->psfMag);
    62         r2rflux->data.F64[Npsf] = PS_SQR(model->radiusTMP) / flux->data.F64[Npsf];
    63 
    64         mask->data.U8[Npsf] = 0;
    65 
    66         // XXX sanity clip?
    67         // XXX need to see if all data were tossed?
    68         // XXX need to subtract median?
    69         // XXX need to put this in the config data...
    70         // if (fabs(apResid->data.F64[Npsf]) > 0.2) continue;
    71         if ((MAX_AP_OFFSET > 0) && (fabs(apResid->data.F64[Npsf]) > MAX_AP_OFFSET)) {
    72             Nfail ++;
    73             continue;
    74         }
    75 
    76         dMag->data.F64[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
    77 
    78         psVectorExtend (mask,    100, 1);
    79         psVectorExtend (xPos,    100, 1);
    80         psVectorExtend (yPos,    100, 1);
    81         psVectorExtend (flux,    100, 1);
    82         psVectorExtend (r2rflux, 100, 1);
    83         psVectorExtend (dMag,    100, 1);
    84         psVectorExtend (apResid, 100, 1);
    85         Npsf ++;
    86     }
    87     psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec for %d objects (%d skipped, %d failed, %ld invalid)\n", 
    88               psTimerMark ("psphot"), Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail);
     39        source = sources->data[i];
     40        model = source->modelPSF;
     41
     42        if (source->type != PM_SOURCE_TYPE_STAR) continue;
     43        if (source->mode &  PM_SOURCE_MODE_SATSTAR) continue;
     44        if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
     45        if (source->mode &  PM_SOURCE_MODE_FAIL) continue;
     46        if (source->mode &  PM_SOURCE_MODE_POOR) continue;
     47
     48        // get growth-corrected, apTrend-uncorrected magnitudes in scaled apertures
     49        // will fail if below S/N threshold or model is missing
     50        if (!pmSourceMagnitudes (source, psf, PM_SOURCE_PHOT_GROWTH)) {
     51            Nskip ++;
     52            continue;
     53        }
     54
     55        apResid->data.F64[Npsf] = source->apMag - source->psfMag;
     56
     57        xPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
     58        yPos->data.F64[Npsf]    = model->params->data.F32[PM_PAR_YPOS];
     59
     60        flux->data.F64[Npsf]    = pow(10.0, -0.4*source->psfMag);
     61        r2rflux->data.F64[Npsf] = PS_SQR(model->radiusTMP) / flux->data.F64[Npsf];
     62
     63        mask->data.U8[Npsf] = 0;
     64
     65        // XXX sanity clip?
     66        // XXX need to see if all data were tossed?
     67        // XXX need to subtract median?
     68        // XXX need to put this in the config data...
     69        // if (fabs(apResid->data.F64[Npsf]) > 0.2) continue;
     70        if ((MAX_AP_OFFSET > 0) && (fabs(apResid->data.F64[Npsf]) > MAX_AP_OFFSET)) {
     71            Nfail ++;
     72            continue;
     73        }
     74
     75        dMag->data.F64[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
     76
     77        psVectorExtend (mask,    100, 1);
     78        psVectorExtend (xPos,    100, 1);
     79        psVectorExtend (yPos,    100, 1);
     80        psVectorExtend (flux,    100, 1);
     81        psVectorExtend (r2rflux, 100, 1);
     82        psVectorExtend (dMag,    100, 1);
     83        psVectorExtend (apResid, 100, 1);
     84        Npsf ++;
     85    }
     86    psLogMsg ("psphot.apresid", 4, "measure aperture residuals : %f sec for %d objects (%d skipped, %d failed, %ld invalid)\n",
     87              psTimerMark ("psphot"), Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail);
    8988
    9089    // XXX choose a better value here?
    9190    if (Npsf < 5) {
    92         psError(PSPHOT_ERR_APERTURE, true, "no valid aperture residual sources, giving up");
    93         return false;
    94     }
    95 
    96     // APTREND options : NONE SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL 
     91        psError(PSPHOT_ERR_APERTURE, true, "no valid aperture residual sources, giving up");
     92        return false;
     93    }
     94
     95    // APTREND options : NONE SKYBIAS XY_LIN XY_QUAD SKY_XY_LIN FULL
    9796    // APTREND options are used in the switch block below
    9897    pmPSFApTrendOptions ApTrendOption = DEFAULT_OPTION;
    … …  
    10099    if (status) ApTrendOption = pmPSFApTrendOptionFromName (optionName);
    101100    if (ApTrendOption == PM_PSF_APTREND_ERROR) {
    102         psError(PSPHOT_ERR_APERTURE, true, "invalid aperture residual trend %s", optionName);
    103         return false;
     101        psError(PSPHOT_ERR_APERTURE, true, "invalid aperture residual trend %s", optionName);
     102        return false;
    104103    }
    105104
    … …  
    112111    switch (ApTrendOption) {
    113112      case PM_PSF_APTREND_NONE:
    114         // remove ApTrend fit from pmPSFtry
    115         psf->ApTrend->coeff[0][0][0][0] = 0;
    116         break;
     113        // remove ApTrend fit from pmPSFtry
     114        psf->ApTrend->coeff[0][0][0][0] = 0;
     115        break;
    117116      case PM_PSF_APTREND_CONSTANT:
    118         stats->clipIter = 2;
    119         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    120         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    121         if (psf->ApTrend == NULL) {
    122             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    123             return false;
    124         }
    125         // apply the fit
    126         stats->clipIter = 3;
    127         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    128         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    129         if (psf->ApTrend == NULL) {
    130             psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction");
    131             return false;
    132         }
    133         break;
     117        stats->clipIter = 2;
     118        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     119        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     120        if (psf->ApTrend == NULL) {
     121            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     122            return false;
     123        }
     124        // apply the fit
     125        stats->clipIter = 3;
     126        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     127        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     128        if (psf->ApTrend == NULL) {
     129            psError(PSPHOT_ERR_PHOTOM, false, "Fitting aperture correction");
     130            return false;
     131        }
     132        break;
    134133      case PM_PSF_APTREND_SKYBIAS:
    135         // first clip out objects which are too far from the median
    136         stats->clipIter = 2;
    137         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    138         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    139         if (psf->ApTrend == NULL) {
    140             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    141             return false;
    142         }
    143         // apply the fit
    144         stats->clipIter = 3;
    145         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    146         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    147         if (psf->ApTrend == NULL) {
    148             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    149             return false;
    150         }
    151         break;
     134        // first clip out objects which are too far from the median
     135        stats->clipIter = 2;
     136        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     137        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     138        if (psf->ApTrend == NULL) {
     139            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     140            return false;
     141        }
     142        // apply the fit
     143        stats->clipIter = 3;
     144        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
     145        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     146        if (psf->ApTrend == NULL) {
     147            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
     148            return false;
     149        }
     150        break;
    152151      case PM_PSF_APTREND_SKYSAT:
    153         // first clip out objects which are too far from the median
    154         stats->clipIter = 2;
    155         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    156         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    157         if (psf->ApTrend == NULL) {
    158             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    159             return false;
    160         }
    161         // apply the fit
    162         stats->clipIter = 2;
    163         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    164         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    165         if (psf->ApTrend == NULL) {
    166             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    167             return false;
    168         }
    169         // apply the fit
    170         stats->clipIter = 3;
    171         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT);
    172         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    173         if (psf->ApTrend == NULL) {
    174             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat");
    175             return false;
    176         }
    177         break;
     152        // first clip out objects which are too far from the median
     153        stats->clipIter = 2;
     154        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     155        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     156        if (psf->ApTrend == NULL) {
     157            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     158            return false;
     159        }
     160        // apply the fit
     161        stats->clipIter = 2;
     162        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
     163        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     164        if (psf->ApTrend == NULL) {
     165            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
     166            return false;
     167        }
     168        // apply the fit
     169        stats->clipIter = 3;
     170        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT);
     171        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     172        if (psf->ApTrend == NULL) {
     173            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting skysat");
     174            return false;
     175        }
     176        break;
    178177      case PM_PSF_APTREND_XY_LIN:
    179         // first clip out objects which are too far from the median
    180         stats->clipIter = 2;
    181         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    182         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    183         if (psf->ApTrend == NULL) {
    184             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    185             return false;
    186         }
    187         // apply the fit
    188         stats->clipIter = 3;
    189         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN);
    190         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    191         if (psf->ApTrend == NULL) {
    192             psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN");
    193             return false;
    194         }
    195         break;
     178        // first clip out objects which are too far from the median
     179        stats->clipIter = 2;
     180        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     181        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     182        if (psf->ApTrend == NULL) {
     183            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     184            return false;
     185        }
     186        // apply the fit
     187        stats->clipIter = 3;
     188        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_LIN);
     189        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     190        if (psf->ApTrend == NULL) {
     191            psError(PSPHOT_ERR_PHOTOM, false, "fitting, XY_LIN");
     192            return false;
     193        }
     194        break;
    196195      case PM_PSF_APTREND_XY_QUAD:
    197         // first clip out objects which are too far from the median
    198         stats->clipIter = 2;
    199         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    200         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    201         if (psf->ApTrend == NULL) {
    202             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    203             return false;
    204         }
    205         // apply the fit
    206         stats->clipIter = 3;
    207         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD);
    208         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    209         if (psf->ApTrend == NULL) {
    210             psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD");
    211             return false;
    212         }
    213         break;
     196        // first clip out objects which are too far from the median
     197        stats->clipIter = 2;
     198        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     199        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     200        if (psf->ApTrend == NULL) {
     201            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     202            return false;
     203        }
     204        // apply the fit
     205        stats->clipIter = 3;
     206        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_XY_QUAD);
     207        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     208        if (psf->ApTrend == NULL) {
     209            psError(PSPHOT_ERR_PHOTOM, false, "Fitting XY_QUAD");
     210            return false;
     211        }
     212        break;
    214213      case PM_PSF_APTREND_SKY_XY_LIN:
    215         // first clip out objects which are too far from the median
    216         stats->clipIter = 2;
    217         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    218         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    219         if (psf->ApTrend == NULL) {
    220             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    221             return false;
    222         }
    223         // apply the fit
    224         stats->clipIter = 3;
    225         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN);
    226         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    227         if (psf->ApTrend == NULL) {
    228             psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin");
    229             return false;
    230         }
    231         break;
     214        // first clip out objects which are too far from the median
     215        stats->clipIter = 2;
     216        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     217        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     218        if (psf->ApTrend == NULL) {
     219            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     220            return false;
     221        }
     222        // apply the fit
     223        stats->clipIter = 3;
     224        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKY_XY_LIN);
     225        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     226        if (psf->ApTrend == NULL) {
     227            psError(PSPHOT_ERR_PHOTOM, false, "Fitting sky xy_lin");
     228            return false;
     229        }
     230        break;
    232231      case PM_PSF_APTREND_SKYSAT_XY_LIN:
    233         // first clip out objects which are too far from the median
    234         stats->clipIter = 2;
    235         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    236         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    237         if (psf->ApTrend == NULL) {
    238             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
    239             return false;
    240         }
    241         // apply the fit
    242         stats->clipIter = 3;
    243         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    244         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    245         if (psf->ApTrend == NULL) {
    246             psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
    247             return false;
    248         }
    249         // apply the fit
    250         stats->clipIter = 3;
    251         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN);
    252         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    253         if (psf->ApTrend == NULL) {
    254             psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin");
    255             return false;
    256         }
    257         break;
     232        // first clip out objects which are too far from the median
     233        stats->clipIter = 2;
     234        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     235        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     236        if (psf->ApTrend == NULL) {
     237            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting nothing");
     238            return false;
     239        }
     240        // apply the fit
     241        stats->clipIter = 3;
     242        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
     243        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     244        if (psf->ApTrend == NULL) {
     245            psError(PSPHOT_ERR_PHOTOM, false, "clipping, fitting sky bias");
     246            return false;
     247        }
     248        // apply the fit
     249        stats->clipIter = 3;
     250        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYSAT_XY_LIN);
     251        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     252        if (psf->ApTrend == NULL) {
     253            psError(PSPHOT_ERR_PHOTOM, false, "Fitting skyset xy_lin");
     254            return false;
     255        }
     256        break;
    258257      case PM_PSF_APTREND_ALL:
    259         // first clip out objects which are too far from the median
    260         stats->clipIter = 2;
    261         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
    262         psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    263         if (psf->ApTrend == NULL) {
    264             psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend");
    265             return false;
    266         }
    267         // fit just SkyBias and clip out objects which are too far from the median
    268         stats->clipIter = 2;
    269         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
    270         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    271         if (psf->ApTrend == NULL) {
    272             psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias");
    273             return false;
    274         }
    275         // finally, fit x, y, SkyBias and clip out objects which are too far from the median
    276         stats->clipIter = 3;
    277         pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL);
    278         psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
    279         if (psf->ApTrend == NULL) {
    280             psError(PSPHOT_ERR_PHOTOM, false, "fitting all");
    281             return false;
    282         }
    283         break;
     258        // first clip out objects which are too far from the median
     259        stats->clipIter = 2;
     260        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_CONSTANT);
     261        psf->ApTrend  = psVectorClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     262        if (psf->ApTrend == NULL) {
     263            psError(PSPHOT_ERR_PHOTOM, false, "Failed to measure apTrend");
     264            return false;
     265        }
     266        // fit just SkyBias and clip out objects which are too far from the median
     267        stats->clipIter = 2;
     268        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_SKYBIAS);
     269        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     270        if (psf->ApTrend == NULL) {
     271            psError(PSPHOT_ERR_PHOTOM, false, "fitting skyBias");
     272            return false;
     273        }
     274        // finally, fit x, y, SkyBias and clip out objects which are too far from the median
     275        stats->clipIter = 3;
     276        pmPSFMaskApTrend (psf->ApTrend, PM_PSF_APTREND_ALL);
     277        psf->ApTrend  = psVectorChiClipFitPolynomial4D (psf->ApTrend, stats, mask, PSFTRY_MASK_ALL, apResid, dMag, xPos, yPos, r2rflux, flux);
     278        if (psf->ApTrend == NULL) {
     279            psError(PSPHOT_ERR_PHOTOM, false, "fitting all");
     280            return false;
     281        }
     282        break;
    284283      default:
    285         psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", optionName);
    286         return false;
     284        psError(PSPHOT_ERR_PHOTOM, true, "Unknown APTREND value: %s", optionName);
     285        return false;
    287286    }
    288287
    … …  
    295294    psStats *residStats = psStatsAlloc (PS_STAT_SAMPLE_STDEV);
    296295    for (int i = 0; i < dMag->n; i++) {
    297         if (dMag->data.F64[i] > 0.01) {
    298             mask->data.U8[i] |= 0x02;
    299         }
    300         if (! mask->data.U8[i]) Nkeep ++;
     296        if (dMag->data.F64[i] > 0.01) {
     297            mask->data.U8[i] |= 0x02;
     298        }
     299        if (! mask->data.U8[i]) Nkeep ++;
    301300    }
    302301    residStats  = psVectorStats (residStats, resid, NULL, mask, 0x03);
    … …  
    320319
    321320    psLogMsg ("psphot.apresid", 3, "measure full-frame aperture residual: %f sec\n", psTimerMark ("psphot"));
    322     psLogMsg ("psphot.apresid", 4, "aperture residual: %f +/- %f : %f bias, %f skysat (%d of %d used)\n", 
    323               psf->ApResid, psf->dApResid, psf->skyBias, psf->skySat, Nkeep, Npsf);
    324     psLogMsg ("psphot.apresid", 4, "apresid trends: %f %f %f %f %f\n", 
    325               1e3*psf->ApTrend->coeff[1][0][0][0],
    326               1e6*psf->ApTrend->coeff[2][0][0][0],
    327               1e6*psf->ApTrend->coeff[1][1][0][0],
    328               1e3*psf->ApTrend->coeff[0][1][0][0],
    329               1e6*psf->ApTrend->coeff[0][2][0][0]);
     321    psLogMsg ("psphot.apresid", 4, "aperture residual: %f +/- %f : %f bias, %f skysat (%d of %d used)\n",
     322              psf->ApResid, psf->dApResid, psf->skyBias, psf->skySat, Nkeep, Npsf);
     323    psLogMsg ("psphot.apresid", 4, "apresid trends: %f %f %f %f %f\n",
     324              1e3*psf->ApTrend->coeff[1][0][0][0],
     325              1e6*psf->ApTrend->coeff[2][0][0][0],
     326              1e6*psf->ApTrend->coeff[1][1][0][0],
     327              1e3*psf->ApTrend->coeff[0][1][0][0],
     328              1e6*psf->ApTrend->coeff[0][2][0][0]);
    330329
    331330    psFree (mask);

  • trunk/psphot/src/psphotBasicDeblend.c

    r9270 r9734  
    22
    33// 2006.02.07 : no leaks
    4 bool psphotBasicDeblend (psArray *sources, psMetadata *recipe) { 
     4bool psphotBasicDeblend (psArray *sources, psMetadata *recipe) {
    55
    66    int N;
    … …  
    2525    // we use an index for this so the spatial sorting is kept
    2626    psVector *SN = psVectorAlloc (sources->n, PS_DATA_F32);
    27     SN->n = SN->nalloc;
    2827    for (int i = 0; i < SN->n; i++) {
    29         source = sources->data[i];
    30         SN->data.F32[i] = source->peak->counts;
     28        source = sources->data[i];
     29        SN->data.F32[i] = source->peak->counts;
    3130    }
    3231    psVector *index = psVectorSortIndex (NULL, SN);
    … …  
    3534    // examine sources in decreasing SN order
    3635    for (int i = sources->n - 1; i >= 0; i--) {
    37         N = index->data.U32[i];
    38         source = sources->data[N];
     36        N = index->data.U32[i];
     37        source = sources->data[N];
    3938
    40         if (source->mode & PM_SOURCE_MODE_BLEND) continue;
     39        if (source->mode & PM_SOURCE_MODE_BLEND) continue;
    4140
    42         // temporary array for overlapping objects we find
    43         psArray *overlap = psArrayAlloc (100);
    44         // DROP overlap->n = 0;
     41        // temporary array for overlapping objects we find
     42        psArray *overlap = psArrayAllocEmpty (100);
    4543
    46         // search backwards for overlapping sources
    47         for (int j = N - 1; j >= 0; j--) {
    48             testSource = sources->data[j];
    49             if (testSource->peak->x <  source->pixels->col0) continue;
    50             if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue;
    51             if (testSource->peak->y <  source->pixels->row0) break;
    52             if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) {
    53                 fprintf (stderr, "warning: invalid condition\n");
    54                 continue;
    55             }
    56             psArrayAdd (overlap, 100, testSource);
    57         }
     44        // search backwards for overlapping sources
     45        for (int j = N - 1; j >= 0; j--) {
     46            testSource = sources->data[j];
     47            if (testSource->peak->x <  source->pixels->col0) continue;
     48            if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue;
     49            if (testSource->peak->y <  source->pixels->row0) break;
     50            if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) {
     51                fprintf (stderr, "warning: invalid condition\n");
     52                continue;
     53            }
     54            psArrayAdd (overlap, 100, testSource);
     55        }
    5856
    59         // search forwards for overlapping sources
    60         for (int j = N + 1; j < sources->n; j++) {
    61             testSource = sources->data[j];
    62             if (testSource->peak->x <  source->pixels->col0) continue;
    63             if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue;
    64             if (testSource->peak->y <  source->pixels->row0) {
    65                 fprintf (stderr, "warning: invalid condition\n");
    66                 continue;
    67             }
    68             if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) break;
    69             psArrayAdd (overlap, 100, testSource);
    70         }
     57        // search forwards for overlapping sources
     58        for (int j = N + 1; j < sources->n; j++) {
     59            testSource = sources->data[j];
     60            if (testSource->peak->x <  source->pixels->col0) continue;
     61            if (testSource->peak->x >= source->pixels->col0 + source->pixels->numCols) continue;
     62            if (testSource->peak->y <  source->pixels->row0) {
     63                fprintf (stderr, "warning: invalid condition\n");
     64                continue;
     65            }
     66            if (testSource->peak->y >= source->pixels->row0 + source->pixels->numRows) break;
     67            psArrayAdd (overlap, 100, testSource);
     68        }
    7169
    72         if (overlap->n == 0) {
    73             psFree (overlap);
    74             continue;
    75         }
     70        if (overlap->n == 0) {
     71            psFree (overlap);
     72            continue;
     73        }
    7674
    77         // this source has overlapping neighbors, check for actual blends
    78         // generate source contour (1/4 peak counts)
    79         // set the threshold based on user inputs
    80        
    81         // threshold is fraction of the source peak flux
    82         // image is background subtracted; source->moments->Sky should always be 0.0
    83         threshold = FRACTION * source->moments->Peak;
    84         // threshold is no less than NSIGMA above the local median sigma?
    85         threshold = PS_MAX (threshold, NSIGMA*sqrt(source->moments->dSky));
     75        // this source has overlapping neighbors, check for actual blends
     76        // generate source contour (1/4 peak counts)
     77        // set the threshold based on user inputs
    8678
    87         // generate a basic contour (set of x,y coordinates at-or-below flux level)
    88         psArray *contour = pmSourceContour (source->pixels, source->peak->x, source->peak->y, threshold);
    89         if (contour == NULL) {
    90             psFree (overlap);
    91             continue;
    92         }
     79        // threshold is fraction of the source peak flux
     80        // image is background subtracted; source->moments->Sky should always be 0.0
     81        threshold = FRACTION * source->moments->Peak;
     82        // threshold is no less than NSIGMA above the local median sigma?
     83        threshold = PS_MAX (threshold, NSIGMA*sqrt(source->moments->dSky));
    9384
    94         // the source contour consists of two vectors, xv and yv.  the contour is
    95         // a series of coordinate pairs, (xv[i],yv[i]) & (xv[i+1],yv[i+1]).  both
    96         // coordinate pairs have the same yv[] value, with xv[i] corresponding to
    97         // the left boundary and xv[i+1] corresponding to the right boundary
     85        // generate a basic contour (set of x,y coordinates at-or-below flux level)
     86        psArray *contour = pmSourceContour (source->pixels, source->peak->x, source->peak->y, threshold);
     87        if (contour == NULL) {
     88            psFree (overlap);
     89            continue;
     90        }
    9891
    99         psVector *xv = contour->data[0];
    100         psVector *yv = contour->data[1];
    101         for (int k = 0; k < overlap->n; k++) {
    102             testSource = overlap->data[k];
    103             if (testSource->peak->counts > source->peak->counts) continue;
    104             for (int j = 0; j < xv->n; j+=2) {
    105                 if (fabs(yv->data.F32[j] - testSource->peak->y) > 0.5) continue;
    106                 if (xv->data.F32[j+0] > testSource->peak->x) break;
    107                 if (xv->data.F32[j+1] < testSource->peak->x) break;
     92        // the source contour consists of two vectors, xv and yv.  the contour is
     93        // a series of coordinate pairs, (xv[i],yv[i]) & (xv[i+1],yv[i+1]).  both
     94        // coordinate pairs have the same yv[] value, with xv[i] corresponding to
     95        // the left boundary and xv[i+1] corresponding to the right boundary
    10896
    109                 testSource->mode |= PM_SOURCE_MODE_BLEND;
     97        psVector *xv = contour->data[0];
     98        psVector *yv = contour->data[1];
     99        for (int k = 0; k < overlap->n; k++) {
     100            testSource = overlap->data[k];
     101            if (testSource->peak->counts > source->peak->counts) continue;
     102            for (int j = 0; j < xv->n; j+=2) {
     103                if (fabs(yv->data.F32[j] - testSource->peak->y) > 0.5) continue;
     104                if (xv->data.F32[j+0] > testSource->peak->x) break;
     105                if (xv->data.F32[j+1] < testSource->peak->x) break;
    110106
    111                 // add this to the list of source->blends
    112                 if (source->blends == NULL) {
    113                     source->blends = psArrayAlloc (16);
    114                     // DROP source->blends->n = 0;
    115                 }
    116                 psArrayAdd (source->blends, 16, testSource);
     107                testSource->mode |= PM_SOURCE_MODE_BLEND;
    117108
    118                 Nblend ++;
    119                 j = xv->n;
    120             }
    121         } 
    122         psFree (overlap);
    123         psFree (contour);
     109                // add this to the list of source->blends
     110                if (source->blends == NULL) {
     111                    source->blends = psArrayAllocEmpty (16);
     112                }
     113                psArrayAdd (source->blends, 16, testSource);
     114
     115                Nblend ++;
     116                j = xv->n;
     117            }
     118        }
     119        psFree (overlap);
     120        psFree (contour);
    124121    }
    125122    psLogMsg ("psphot.deblend", 3, "identified %d blended objects (%f sec)\n", Nblend, psTimerMark ("psphot"));

  • trunk/psphot/src/psphotChoosePSF.c

    r9661 r9734  
    5050    pmSourceFitModelInit (15, 0.1, POISSON_ERRORS);
    5151
    52     stars = psArrayAlloc (sources->n);
     52    stars = psArrayAllocEmpty (sources->n);
    5353
    5454    // select the candidate PSF stars (pointers to original sources)
    … …  
    8787    // set up an array to store the results
    8888    psArray *models = psArrayAlloc (list->n);
    89     models->n = list->n;
    9089
    9190    // try each model option listed in config

  • trunk/psphot/src/psphotEnsemblePSF.c

    r9595 r9734  
    2727
    2828    // pre-calculate all model pixels
    29     psArray  *models = psArrayAlloc (sources->n);
    30     psVector *index  = psVectorAlloc (sources->n, PS_TYPE_U32);
     29    psArray  *models = psArrayAllocEmpty (sources->n);
     30    psVector *index  = psVectorAllocEmpty (sources->n, PS_TYPE_U32);
    3131    // DROP models->n = index->n = 0;
    3232
    … …  
    141141    psSparse *sparse = psSparseAlloc (models->n, 100);
    142142    psVector *weight = psVectorAlloc (models->n, PS_TYPE_F32);
    143     weight->n = models->n;
    144143
    145144    for (int i = 0; i < models->n; i++) {

  • trunk/psphot/src/psphotFakeSources.c

    r9270 r9734  
    66
    77    psArray *sources = psArrayAlloc (50);
    8     sources->n = 50;
    98
    109    for (int i = 0; i < sources->n; i++) {
    11         pmSource *source = pmSourceAlloc ();
    12         source->moments = pmMomentsAlloc ();
    13         source->moments->x = 10;
    14         source->moments->y = 10;
    15         source->moments->Sx = 1;
    16         source->moments->Sy = 1;
    17         source->moments->Sxy = 0;
    18         source->moments->Sum = 1000;
    19         source->moments->Peak = 100;
    20         source->moments->Sky = 10;
    21         source->moments->nPixels = 10;
     10        pmSource *source = pmSourceAlloc ();
     11        source->moments = pmMomentsAlloc ();
     12        source->moments->x = 10;
     13        source->moments->y = 10;
     14        source->moments->Sx = 1;
     15        source->moments->Sy = 1;
     16        source->moments->Sxy = 0;
     17        source->moments->Sum = 1000;
     18        source->moments->Peak = 100;
     19        source->moments->Sky = 10;
     20        source->moments->nPixels = 10;
    2221
    23         source->peak = pmPeakAlloc (10, 10, 0, 0);
    24         source->type = PM_SOURCE_TYPE_STAR;
     22        source->peak = pmPeakAlloc (10, 10, 0, 0);
     23        source->type = PM_SOURCE_TYPE_STAR;
    2524
    26         pmModelType modelType = pmModelSetType ("PS_MODEL_QGAUSS");
    27         source->modelPSF = pmSourceModelGuess (source, modelType);
    28         sources->data[i] = source;
     25        pmModelType modelType = pmModelSetType ("PS_MODEL_QGAUSS");
     26        source->modelPSF = pmSourceModelGuess (source, modelType);
     27        sources->data[i] = source;
    2928    }
    3029    return sources;

  • trunk/psphot/src/psphotFitSet.c

    r9270 r9734  
    88    if (f == NULL) return false;
    99
    10     psArray *modelSet = psArrayAlloc (16);
    11     // DROP modelSet->n = 0;
     10    psArray *modelSet = psArrayAllocEmpty (16);
    1211
    1312    while (fscanf (f, "%lf %lf %lf", &x, &y, &Io) == 3) {
    14         pmModel *model = pmModelAlloc (oneModel->type);
    15        
    16         for (psS32 i = 0; i < model->params->n; i++) {
    17             model->params->data.F32[i] = oneModel->params->data.F32[i];
    18             model->dparams->data.F32[i] = oneModel->dparams->data.F32[i];
    19         }
    20         model->params->data.F32[1] = Io;
    21         model->params->data.F32[2] = x;
    22         model->params->data.F32[3] = y;
    23         psArrayAdd (modelSet, 16, model);
     13        pmModel *model = pmModelAlloc (oneModel->type);
     14
     15        for (psS32 i = 0; i < model->params->n; i++) {
     16            model->params->data.F32[i] = oneModel->params->data.F32[i];
     17            model->dparams->data.F32[i] = oneModel->dparams->data.F32[i];
     18        }
     19        model->params->data.F32[1] = Io;
     20        model->params->data.F32[2] = x;
     21        model->params->data.F32[3] = y;
     22        psArrayAdd (modelSet, 16, model);
    2423    }
    2524
    … …  
    3130    // subtract object, leave local sky
    3231    for (int i = 0; i < modelSet->n; i++) {
    33         pmModel *model = modelSet->data[i];
    34         pmModelSub (source->pixels, source->mask, model, false, false);
    35    
    36         fprintf (stderr, "output parameters (obj %d):\n", i);
    37         for (int n = 0; n < model->params->n; n++) {
    38             fprintf (stderr, "%d : %f\n", n, model->params->data.F32[n]);
    39         }
     32        pmModel *model = modelSet->data[i];
     33        pmModelSub (source->pixels, source->mask, model, false, false);
     34
     35        fprintf (stderr, "output parameters (obj %d):\n", i);
     36        for (int n = 0; n < model->params->n; n++) {
     37            fprintf (stderr, "%d : %f\n", n, model->params->data.F32[n]);
     38        }
    4039    }
    4140
    42     // write out 
     41    // write out
    4342    psphotSaveImage (NULL, source->pixels, "resid.fits");
    4443    psphotSaveImage (NULL, source->mask, "mask.fits");

  • trunk/psphot/src/psphotSourceFits.c

    r9595 r9734  
    2121    y = PSF->params->data.F32[3];
    2222
    23     psArray *modelSet = psArrayAlloc (source->blends->n + 1);
     23    psArray *modelSet = psArrayAllocEmpty (source->blends->n + 1);
    2424    psArrayAdd (modelSet, 16, PSF);
    2525
    26     psArray *sourceSet = psArrayAlloc (source->blends->n + 1);
     26    psArray *sourceSet = psArrayAllocEmpty (source->blends->n + 1);
    2727    psArrayAdd (sourceSet, 16, source);
    2828
    … …  
    300300
    301301    modelSet = psArrayAlloc (2);
    302     modelSet->n = 2;
    303302
    304303    DBL = pmModelCopy (PSF);

  • trunk/psphot/src/psphotSourceStats.c

    r9520 r9734  
    22
    33// 2006.02.02 : no leaks
    4 psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, psArray *peaks) 
     4psArray *psphotSourceStats (pmReadout *readout, psMetadata *recipe, psArray *peaks)
    55{
    66    bool     status  = false;
    … …  
    1616    char *breakPt  = psMetadataLookupStr (&status, recipe, "BREAK_POINT");
    1717
    18     sources = psArrayAlloc (peaks->n);
     18    sources = psArrayAllocEmpty (peaks->n);
    1919
    2020    for (int i = 0; i < peaks->n; i++) {
    2121
    22         // create a new source, add peak
    23         pmSource *source = pmSourceAlloc();
    24         source->peak = (pmPeak *)psMemIncrRefCounter(peaks->data[i]);
     22        // create a new source, add peak
     23        pmSource *source = pmSourceAlloc();
     24        source->peak = (pmPeak *)psMemIncrRefCounter(peaks->data[i]);
    2525
    26         // allocate image, weight, mask arrays for each peak (square of radius OUTER)
    27         pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
    28         if (!strcasecmp (breakPt, "PEAKS")) {
    29             psArrayAdd (sources, 100, source);
    30             psFree (source);
    31             continue;
    32         }
     26        // allocate image, weight, mask arrays for each peak (square of radius OUTER)
     27        pmSourceDefinePixels (source, readout, source->peak->x, source->peak->y, OUTER);
     28        if (!strcasecmp (breakPt, "PEAKS")) {
     29            psArrayAdd (sources, 100, source);
     30            psFree (source);
     31            continue;
     32        }
    3333
    34         // XXX skip faint sources?
     34        // XXX skip faint sources?
    3535
    36         // measure a local sky value
    37         // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
    38         status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER);
    39         if (!status) {
    40           psFree (source);
    41           continue;
    42         }
     36        // measure a local sky value
     37        // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
     38        status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER);
     39        if (!status) {
     40          psFree (source);
     41          continue;
     42        }
    4343
    44         // measure the local sky variance (needed if noise is not sqrt(signal))
    45         // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
    46         status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER);
    47         if (!status) {
    48           psFree (source);
    49           continue;
    50         }
     44        // measure the local sky variance (needed if noise is not sqrt(signal))
     45        // XXX EAM : this should use ROBUST not SAMPLE median, but it is broken
     46        status = pmSourceLocalSkyVariance (source, PS_STAT_SAMPLE_MEDIAN, INNER);
     47        if (!status) {
     48          psFree (source);
     49          continue;
     50        }
    5151
    52         // measure basic source moments
    53         status = pmSourceMoments (source, RADIUS);
    54         if (status) {
    55             // add to the source array
    56             psArrayAdd (sources, 100, source);
    57             psFree (source);
    58             continue;
    59         }
     52        // measure basic source moments
     53        status = pmSourceMoments (source, RADIUS);
     54        if (status) {
     55            // add to the source array
     56            psArrayAdd (sources, 100, source);
     57            psFree (source);
     58            continue;
     59        }
    6060
    61         // if no valid pixels, or massive swing, likely saturated source,
    62         // try a much larger box
    63         BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
    64         psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
    65         status = pmSourceMoments (source, BIG_RADIUS);
    66         if (status) {
    67             // add to the source array
    68             psArrayAdd (sources, 100, source);
    69             psFree (source);
    70             continue;
    71         }
     61        // if no valid pixels, or massive swing, likely saturated source,
     62        // try a much larger box
     63        BIG_RADIUS = PS_MIN (INNER, 3*RADIUS);
     64        psTrace ("psphot", 4, "retrying moments for %d, %d\n", source->peak->x, source->peak->y);
     65        status = pmSourceMoments (source, BIG_RADIUS);
     66        if (status) {
     67            // add to the source array
     68            psArrayAdd (sources, 100, source);
     69            psFree (source);
     70            continue;
     71        }
    7272
    73         psFree (source);
    74         continue;
     73        psFree (source);
     74        continue;
    7575    }
    76      
     76
    7777    psLogMsg ("psphot", 3, "%ld moments: %f sec\n", sources->n, psTimerMark ("psphot"));
    7878

  • trunk/psphot/src/psphotTest.c

    r9270 r9734  
    1515int main (int argc, char **argv) {
    1616
    17     psRegion region = {0,0,0,0};        // a region representing the entire array
     17    psRegion region = {0,0,0,0};        // a region representing the entire array
    1818    psphotTestArguments (&argc, argv);
    1919
    … …  
    5656    mdi->data.S32 = 0;
    5757    mdi->type = PS_DATA_S32;
    58    
     58
    5959    // create a test image
    6060    // psImage *tmpimage = psImageAlloc (10, 10, PS_DATA_F32);
    6161
    6262    // create a test table
    63     table = psArrayAlloc (10);
    64 table->n = 0;
     63    table = psArrayAllocEmpty (10);
    6564
    6665    for (int i = 0; i < 10; i++) {
    67         row = psMetadataAlloc ();
    68         psMetadataAdd (row, PS_LIST_TAIL, "ROW",   PS_DATA_S32,    "", i);
    69         psMetadataAdd (row, PS_LIST_TAIL, "FROW",  PS_TYPE_F32,    "", 0.1*i);
    70         psMetadataAdd (row, PS_LIST_TAIL, "DUMMY", PS_DATA_STRING, "", "test line");
    71    
    72         table->data[i] = row;
     66        row = psMetadataAlloc ();
     67        psMetadataAdd (row, PS_LIST_TAIL, "ROW",   PS_DATA_S32,    "", i);
     68        psMetadataAdd (row, PS_LIST_TAIL, "FROW",  PS_TYPE_F32,    "", 0.1*i);
     69        psMetadataAdd (row, PS_LIST_TAIL, "DUMMY", PS_DATA_STRING, "", "test line");
     70
     71        table->data[i] = row;
    7372    }
    7473    table->n = 10;
    75    
     74
    7675    psMetadata *theader = psMetadataAlloc ();
    7776    psMetadataAdd (theader, PS_LIST_HEAD, "EXTNAME", PS_DATA_STRING, "extension name", "SMPFILE");
Note: See TracChangeset for help on using the changeset viewer.