IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 20170


Ignore:
Timestamp:
Oct 15, 2008, 9:59:10 AM (18 years ago)
Author:
Paul Price
Message:

Having psphotReadoutCleanup return false when there's an error still in the stack requires that possible sources of errors be checked.

Location:
trunk/psphot/src
Files:
2 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotImageQuality.c

    r19951 r20170  
    1 # include "psphotInternal.h"
    2 # define USE_SAMPLE 1 /* sample vs robust stats */
     1#include "psphotInternal.h"
     2#define USE_SAMPLE 1 /* sample vs robust stats */
     3
     4// XXX Lots of code duplication here
    35
    46// selecting the 'good' stars (likely to be psf stars), measure the M_cn, M_sn terms for n = 2,3,4
    5 bool psphotImageQuality (psMetadata *recipe, psArray *sources) {
    6 
    7     psVector *FWHM_MAJOR = psVectorAllocEmpty (100, PS_TYPE_F32);
    8     psVector *FWHM_MINOR = psVectorAllocEmpty (100, PS_TYPE_F32);
    9 
    10     psVector *M2 = psVectorAllocEmpty (100, PS_TYPE_F32);
    11     psVector *M3 = psVectorAllocEmpty (100, PS_TYPE_F32);
    12     psVector *M4 = psVectorAllocEmpty (100, PS_TYPE_F32);
    13 
    14     psVector *M2c = psVectorAllocEmpty (100, PS_TYPE_F32);
    15     psVector *M2s = psVectorAllocEmpty (100, PS_TYPE_F32);
    16 
     7bool psphotImageQuality(psMetadata *recipe, psArray *sources)
     8{
     9    psVector *FWHM_MAJOR = psVectorAllocEmpty(100, PS_TYPE_F32);
     10    psVector *FWHM_MINOR = psVectorAllocEmpty(100, PS_TYPE_F32);
     11
     12    psVector *M2 = psVectorAllocEmpty(100, PS_TYPE_F32);
     13    psVector *M3 = psVectorAllocEmpty(100, PS_TYPE_F32);
     14    psVector *M4 = psVectorAllocEmpty(100, PS_TYPE_F32);
     15
     16    psVector *M2c = psVectorAllocEmpty(100, PS_TYPE_F32);
     17    psVector *M2s = psVectorAllocEmpty(100, PS_TYPE_F32);
     18
     19    int num = 0;                        // Number of good sources
    1720    for (int i = 0; i < sources->n; i++) {
    1821
    19         pmSource *source = sources->data[i];
    20         if (!source) continue;
    21 
    22         // select by S/N?
    23         // select by PSFSTAR mode?
    24         // ??
    25         if (source->type != PM_SOURCE_TYPE_STAR) continue;
    26         if (!(source->mode & PM_SOURCE_MODE_PSFSTAR)) continue;
    27 
    28         pmMoments *moments = source->moments;
    29         if (!moments) continue;
    30 
    31         // if (source->moments->SN < XXX) continue;
    32 
    33         // M_cn = \sum (f r^2 cos(n t)) / \sum (f)
    34         // M_sn = \sum (f r^2 sin(n t)) / \sum (f)
    35 
    36         // r^2 cos(2t) = r^2 cos^2 t - r^2 sin^2 t = x^2 - y^2
    37         // r^2 sin(2t) = r^2 2 cos t sin t = 2 x y
    38 
    39         // r^2 cos(3t) = r^2 cos^3 t - r^2 3 cos t sin^2 t = (x^3 - 3 x y^2) / r
    40         // r^2 sin(3t) = r^2 3 cos^2 t sin t - sin^3 t = (3 x^2 y - y^3) / r
    41 
    42         // r^2 cos(4t) = r^2 cos^4 t - r^2 6 cos^2 t sin^2 t + r^2 sin^4 t = (x^4 - 6 x^2 y^2 + y^4) / r^2
    43         // r^2 sin(4t) = r^2 4 cos^3 t sin t - 4 sin^3 t cos t = (4 x^3 y - 4 y^3 x) / r^2
    44 
    45         psEllipseMoments emoments;
    46         emoments.x2 = moments->Mxx;
    47         emoments.xy = moments->Mxy;
    48         emoments.y2 = moments->Myy;
    49         psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0);
    50 
    51         psVectorAppend (FWHM_MAJOR, 2.355*axes.major);
    52         psVectorAppend (FWHM_MINOR, 2.355*axes.minor);
    53 
    54         // the moments->Mnnn terms contain the straight sums: for example: moments->Mxxx contains \sum x^3 / r
    55         float M_c2 = moments->Mxx - moments->Myy;
    56         float M_s2 = 2*moments->Mxy;
    57         psVectorAppend (M2, hypot(M_c2, M_s2));
    58         psVectorAppend (M2c, M_c2);
    59         psVectorAppend (M2s, M_s2);
    60        
    61         float M_c3 = moments->Mxxx - 3.0*moments->Mxyy;
    62         float M_s3 = 3.0*moments->Mxxy - moments->Myyy;
    63         psVectorAppend (M3, hypot (M_c3, M_s3));
    64        
    65         float M_c4 = moments->Mxxxx - 6.0*moments->Mxxyy + moments->Myyyy;
    66         float M_s4 = 4.0*moments->Mxxxy - 4.0*moments->Mxyyy;
    67         psVectorAppend (M4, hypot (M_c4, M_s4));
    68     }   
    69 
    70     psMetadataAddS32 (recipe, PS_LIST_TAIL, "IQ_NSTAR", PS_META_REPLACE, "Number of stars used for IQ measurements", M2->n);
    71 
    72 # if (USE_SAMPLE)
     22        pmSource *source = sources->data[i];
     23        if (!source) {
     24            continue;
     25        }
     26
     27        // select by S/N?
     28        // select by PSFSTAR mode?
     29        // ??
     30        if (source->type != PM_SOURCE_TYPE_STAR || !(source->mode & PM_SOURCE_MODE_PSFSTAR)) {
     31            continue;
     32        }
     33
     34        pmMoments *moments = source->moments;
     35        if (!moments) {
     36            continue;
     37        }
     38
     39        // if (source->moments->SN < XXX) continue;
     40
     41        // M_cn = \sum (f r^2 cos(n t)) / \sum (f)
     42        // M_sn = \sum (f r^2 sin(n t)) / \sum (f)
     43
     44        // r^2 cos(2t) = r^2 cos^2 t - r^2 sin^2 t = x^2 - y^2
     45        // r^2 sin(2t) = r^2 2 cos t sin t = 2 x y
     46
     47        // r^2 cos(3t) = r^2 cos^3 t - r^2 3 cos t sin^2 t = (x^3 - 3 x y^2) / r
     48        // r^2 sin(3t) = r^2 3 cos^2 t sin t - sin^3 t = (3 x^2 y - y^3) / r
     49
     50        // r^2 cos(4t) = r^2 cos^4 t - r^2 6 cos^2 t sin^2 t + r^2 sin^4 t = (x^4 - 6 x^2 y^2 + y^4) / r^2
     51        // r^2 sin(4t) = r^2 4 cos^3 t sin t - 4 sin^3 t cos t = (4 x^3 y - 4 y^3 x) / r^2
     52
     53        num++;
     54
     55        psEllipseMoments emoments;
     56        emoments.x2 = moments->Mxx;
     57        emoments.xy = moments->Mxy;
     58        emoments.y2 = moments->Myy;
     59        psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0);
     60
     61        psVectorAppend (FWHM_MAJOR, 2.355*axes.major);
     62        psVectorAppend (FWHM_MINOR, 2.355*axes.minor);
     63
     64        // the moments->Mnnn terms contain the straight sums: for example: moments->Mxxx contains \sum x^3 / r
     65        float M_c2 = moments->Mxx - moments->Myy;
     66        float M_s2 = 2*moments->Mxy;
     67        psVectorAppend (M2, hypot(M_c2, M_s2));
     68        psVectorAppend (M2c, M_c2);
     69        psVectorAppend (M2s, M_s2);
     70
     71        float M_c3 = moments->Mxxx - 3.0*moments->Mxyy;
     72        float M_s3 = 3.0*moments->Mxxy - moments->Myyy;
     73        psVectorAppend (M3, hypot (M_c3, M_s3));
     74
     75        float M_c4 = moments->Mxxxx - 6.0*moments->Mxxyy + moments->Myyyy;
     76        float M_s4 = 4.0*moments->Mxxxy - 4.0*moments->Mxyyy;
     77        psVectorAppend (M4, hypot (M_c4, M_s4));
     78    }
     79
     80    if (num == 0) {
     81        psError(PSPHOT_ERR_PHOTOM, false,
     82                "Unable to find sources from which to measure image quality");
     83        return false;
     84    }
     85
     86    psMetadataAddS32(recipe, PS_LIST_TAIL, "IQ_NSTAR", PS_META_REPLACE,
     87                     "Number of stars used for IQ measurements", M2->n);
     88
     89#if (USE_SAMPLE)
    7390    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_QUARTILE);
    7491
    75     psVectorStats (stats, FWHM_MAJOR, NULL, NULL, 0);
     92    if (!psVectorStats(stats, FWHM_MAJOR, NULL, NULL, 0)) {
     93        psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
     94        goto FAIL;
     95    }
     96
    7697    float fwhm_major = stats->sampleMean;
    77     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW1",   PS_META_REPLACE, "FWHM of Major Axis from moments", stats->sampleMean);
    78     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW1_E", PS_META_REPLACE, "FWHM scatter (Major) from moments", stats->sampleStdev);
    79 
    80     psVectorStats (stats, FWHM_MINOR, NULL, NULL, 0);
     98    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW1",   PS_META_REPLACE,
     99                     "FWHM of Major Axis from moments", stats->sampleMean);
     100    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW1_E", PS_META_REPLACE,
     101                     "FWHM scatter (Major) from moments", stats->sampleStdev);
     102
     103    if (!psVectorStats(stats, FWHM_MINOR, NULL, NULL, 0)) {
     104        psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
     105        goto FAIL;
     106    }
    81107    float fwhm_minor = stats->sampleMean;
    82     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW2",   PS_META_REPLACE, "FWHM of Minor Axis from moments", stats->sampleMean);
    83     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW2_E", PS_META_REPLACE, "FWHM scatter (Minor) from moments", stats->sampleStdev);
    84 
    85     psVectorStats (stats, M2, NULL, NULL, 0);
     108    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW2",   PS_META_REPLACE,
     109                     "FWHM of Minor Axis from moments", stats->sampleMean);
     110    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW2_E", PS_META_REPLACE,
     111                     "FWHM scatter (Minor) from moments", stats->sampleStdev);
     112
     113    if (!psVectorStats(stats, M2, NULL, NULL, 0)) {
     114        psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
     115        goto FAIL;
     116    }
    86117    float vM2 = stats->sampleMean;
    87118    float dM2 = stats->sampleStdev;
    88     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2",    PS_META_REPLACE, "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2);
    89     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2_ER", PS_META_REPLACE, "Stdev of M_2", dM2);
    90     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2_LQ", PS_META_REPLACE, "Lower Quartile of M_2", stats->sampleLQ);
    91     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2_UQ", PS_META_REPLACE, "Upper Quartile of M_2", stats->sampleUQ);
    92 
    93     psVectorStats (stats, M2c, NULL, NULL, 0);
    94     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C",   PS_META_REPLACE, "M_2c = sum f r^2 cos(2phi) / sum f", stats->sampleMean);
    95     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C_E", PS_META_REPLACE, "Stdev of M_2c", stats->sampleStdev);
    96     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C_L", PS_META_REPLACE, "Lower Quartile of M_2c", stats->sampleLQ);
    97     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C_U", PS_META_REPLACE, "Upper Quartile of M_2c", stats->sampleUQ);
    98 
    99     psVectorStats (stats, M2s, NULL, NULL, 0);
    100     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S",   PS_META_REPLACE, "M_2s = sum f r^2 cos(2phi) / sum f", stats->sampleMean);
    101     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S_E", PS_META_REPLACE, "Stdev of M_2s", stats->sampleStdev);
    102     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S_L", PS_META_REPLACE, "Lower Quartile of M_2s", stats->sampleLQ);
    103     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S_U", PS_META_REPLACE, "Upper Quartile of M_2s", stats->sampleUQ);
    104 
    105     psVectorStats (stats, M3, NULL, NULL, 0);
     119    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2",    PS_META_REPLACE,
     120                     "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2);
     121    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_ER", PS_META_REPLACE,
     122                     "Stdev of M_2", dM2);
     123    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_LQ", PS_META_REPLACE,
     124                     "Lower Quartile of M_2", stats->sampleLQ);
     125    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_UQ", PS_META_REPLACE,
     126                     "Upper Quartile of M_2", stats->sampleUQ);
     127
     128    if (!psVectorStats(stats, M2c, NULL, NULL, 0)) {
     129        psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
     130        goto FAIL;
     131    }
     132    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C",   PS_META_REPLACE,
     133                     "M_2c = sum f r^2 cos(2phi) / sum f", stats->sampleMean);
     134    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_E", PS_META_REPLACE,
     135                     "Stdev of M_2c", stats->sampleStdev);
     136    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_L", PS_META_REPLACE,
     137                     "Lower Quartile of M_2c", stats->sampleLQ);
     138    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_U", PS_META_REPLACE,
     139                     "Upper Quartile of M_2c", stats->sampleUQ);
     140
     141    if (!psVectorStats(stats, M2s, NULL, NULL, 0)) {
     142        psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
     143        goto FAIL;
     144    }
     145    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S",   PS_META_REPLACE,
     146                     "M_2s = sum f r^2 cos(2phi) / sum f", stats->sampleMean);
     147    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_E", PS_META_REPLACE,
     148                     "Stdev of M_2s", stats->sampleStdev);
     149    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_L", PS_META_REPLACE,
     150                     "Lower Quartile of M_2s", stats->sampleLQ);
     151    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_U", PS_META_REPLACE,
     152                     "Upper Quartile of M_2s", stats->sampleUQ);
     153
     154    if (!psVectorStats(stats, M3, NULL, NULL, 0)) {
     155        psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
     156        goto FAIL;
     157    }
    106158    float vM3 = stats->sampleMean;
    107159    float dM3 = stats->sampleStdev;
    108     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3",   PS_META_REPLACE, "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3);
    109     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3_ER", PS_META_REPLACE, "Stdev of M_3", dM3);
    110     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3_LQ", PS_META_REPLACE, "Lower Quartile of M_3", stats->sampleLQ);
    111     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3_UQ", PS_META_REPLACE, "Upper Quartile of M_3", stats->sampleUQ);
    112 
    113     psVectorStats (stats, M4, NULL, NULL, 0);
     160    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3",   PS_META_REPLACE,
     161                     "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3);
     162    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_ER", PS_META_REPLACE,
     163                     "Stdev of M_3", dM3);
     164    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_LQ", PS_META_REPLACE,
     165                     "Lower Quartile of M_3", stats->sampleLQ);
     166    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_UQ", PS_META_REPLACE,
     167                     "Upper Quartile of M_3", stats->sampleUQ);
     168
     169    if (!psVectorStats(stats, M4, NULL, NULL, 0)) {
     170        psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
     171        goto FAIL;
     172    }
    114173    float vM4 = stats->sampleMean;
    115174    float dM4 = stats->sampleStdev;
    116     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4",   PS_META_REPLACE, "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4);
    117     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4_ER", PS_META_REPLACE, "Stdev of M_4", dM4);
    118     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4_LQ", PS_META_REPLACE, "Lower Quartile of M_4", stats->sampleLQ);
    119     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4_UQ", PS_META_REPLACE, "Upper Quartile of M_4", stats->sampleUQ);
    120 
    121 # else
     175    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4",   PS_META_REPLACE,
     176                     "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4);
     177    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_ER", PS_META_REPLACE,
     178                     "Stdev of M_4", dM4);
     179    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_LQ", PS_META_REPLACE,
     180                     "Lower Quartile of M_4", stats->sampleLQ);
     181    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_UQ", PS_META_REPLACE,
     182                     "Upper Quartile of M_4", stats->sampleUQ);
     183
     184#else
    122185    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
    123186
    124     psVectorStats (stats, FWHM_MAJOR, NULL, NULL, 0);
     187    if (!psVectorStats(stats, FWHM_MAJOR, NULL, NULL, 0)) {
     188        psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
     189        goto FAIL;
     190    }
    125191    float fwhm_major = stats->robustMedian;
    126     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW1",   PS_META_REPLACE, "FWHM of Major Axis from moments", stats->robustMedian);
    127     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW1_E", PS_META_REPLACE, "FWHM scatter (Major) from moments", stats->robustStdev);
    128 
    129     psVectorStats (stats, FWHM_MINOR, NULL, NULL, 0);
     192    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW1",   PS_META_REPLACE,
     193                     "FWHM of Major Axis from moments", stats->robustMedian);
     194    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW1_E", PS_META_REPLACE,
     195                     "FWHM scatter (Major) from moments", stats->robustStdev);
     196
     197    if (!psVectorStats(stats, FWHM_MINOR, NULL, NULL, 0)) {
     198        psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
     199        goto FAIL;
     200    }
    130201    float fwhm_minor = stats->robustMedian;
    131     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW2",   PS_META_REPLACE, "FWHM of Minor Axis from moments", stats->robustMedian);
    132     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW2_E", PS_META_REPLACE, "FWHM scatter (Minor) from moments", stats->robustStdev);
    133 
    134     psVectorStats (stats, M2, NULL, NULL, 0);
     202    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW2",   PS_META_REPLACE,
     203                     "FWHM of Minor Axis from moments", stats->robustMedian);
     204    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_FW2_E", PS_META_REPLACE,
     205                     "FWHM scatter (Minor) from moments", stats->robustStdev);
     206
     207    if (!psVectorStats(stats, M2, NULL, NULL, 0)) {
     208        psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
     209        goto FAIL;
     210    }
    135211    float vM2 = stats->robustMedian;
    136212    float dM2 = stats->robustStdev;
    137     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2",    PS_META_REPLACE, "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2);
    138     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2_ER", PS_META_REPLACE, "Stdev of M_2", dM2);
    139     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2_LQ", PS_META_REPLACE, "Lower Quartile of M_2", stats->robustLQ);
    140     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2_UQ", PS_META_REPLACE, "Upper Quartile of M_2", stats->robustUQ);
    141 
    142     psVectorStats (stats, M2c, NULL, NULL, 0);
    143     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C",   PS_META_REPLACE, "M_2c = sum f r^2 cos(2phi) / sum f", stats->robustMedian);
    144     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C_E", PS_META_REPLACE, "Stdev of M_2c", stats->robustStdev);
    145     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C_L", PS_META_REPLACE, "Lower Quartile of M_2c", stats->robustLQ);
    146     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C_U", PS_META_REPLACE, "Upper Quartile of M_2c", stats->robustUQ);
    147 
    148     psVectorStats (stats, M2s, NULL, NULL, 0);
    149     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S",   PS_META_REPLACE, "M_2s = sum f r^2 cos(2phi) / sum f", stats->robustMedian);
    150     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S_E", PS_META_REPLACE, "Stdev of M_2s", stats->robustStdev);
    151     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S_L", PS_META_REPLACE, "Lower Quartile of M_2s", stats->robustLQ);
    152     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S_U", PS_META_REPLACE, "Upper Quartile of M_2s", stats->robustUQ);
    153 
    154     psVectorStats (stats, M3, NULL, NULL, 0);
     213    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2",    PS_META_REPLACE,
     214                     "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2);
     215    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_ER", PS_META_REPLACE,
     216                     "Stdev of M_2", dM2);
     217    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_LQ", PS_META_REPLACE,
     218                     "Lower Quartile of M_2", stats->robustLQ);
     219    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2_UQ", PS_META_REPLACE,
     220                     "Upper Quartile of M_2", stats->robustUQ);
     221
     222    if (!psVectorStats(stats, M2c, NULL, NULL, 0)) {
     223        psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
     224        goto FAIL;
     225    }
     226    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C",   PS_META_REPLACE,
     227                     "M_2c = sum f r^2 cos(2phi) / sum f", stats->robustMedian);
     228    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_E", PS_META_REPLACE,
     229                     "Stdev of M_2c", stats->robustStdev);
     230    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_L", PS_META_REPLACE,
     231                     "Lower Quartile of M_2c", stats->robustLQ);
     232    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2C_U", PS_META_REPLACE,
     233                     "Upper Quartile of M_2c", stats->robustUQ);
     234
     235    if (!psVectorStats(stats, M2s, NULL, NULL, 0)) {
     236        psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
     237        goto FAIL;
     238    }
     239    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S",   PS_META_REPLACE,
     240                     "M_2s = sum f r^2 cos(2phi) / sum f", stats->robustMedian);
     241    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_E", PS_META_REPLACE,
     242                     "Stdev of M_2s", stats->robustStdev);
     243    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_L", PS_META_REPLACE,
     244                     "Lower Quartile of M_2s", stats->robustLQ);
     245    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M2S_U", PS_META_REPLACE,
     246                     "Upper Quartile of M_2s", stats->robustUQ);
     247
     248    if (!psVectorStats(stats, M3, NULL, NULL, 0)) {
     249        psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
     250        goto FAIL;
     251    }
    155252    float vM3 = stats->robustMedian;
    156253    float dM3 = stats->robustStdev;
    157     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3",   PS_META_REPLACE, "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3);
    158     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3_ER", PS_META_REPLACE, "Stdev of M_3", dM3);
    159     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3_LQ", PS_META_REPLACE, "Lower Quartile of M_3", stats->robustLQ);
    160     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3_UQ", PS_META_REPLACE, "Upper Quartile of M_3", stats->robustUQ);
    161 
    162     psVectorStats (stats, M4, NULL, NULL, 0);
     254    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3",   PS_META_REPLACE,
     255                     "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3);
     256    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_ER", PS_META_REPLACE,
     257                     "Stdev of M_3", dM3);
     258    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_LQ", PS_META_REPLACE,
     259                     "Lower Quartile of M_3", stats->robustLQ);
     260    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M3_UQ", PS_META_REPLACE,
     261                     "Upper Quartile of M_3", stats->robustUQ);
     262
     263    if (!psVectorStats(stats, M4, NULL, NULL, 0)) {
     264        psError(PSPHOT_ERR_UNKNOWN, false, "Unable to perform statistics to measure image quality");
     265        goto FAIL;
     266    }
    163267    float vM4 = stats->robustMedian;
    164268    float dM4 = stats->robustStdev;
    165     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4",   PS_META_REPLACE, "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4);
    166     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4_ER", PS_META_REPLACE, "Stdev of M_4", dM4);
    167     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4_LQ", PS_META_REPLACE, "Lower Quartile of M_4", stats->robustLQ);
    168     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4_UQ", PS_META_REPLACE, "Upper Quartile of M_4", stats->robustUQ);
    169 # endif
    170 
    171     psLogMsg ("psphot", PS_LOG_INFO, "Image Quality Stats from %ld psf stars : FWHM (major, minor) : %f, %f  M_2 : %f +/- %f, M_3 : %f +/- %f, M_4 : %f +/- %f\n", M2->n, fwhm_major, fwhm_minor, vM2, dM2, vM3, dM3, vM4, dM4);
    172 
    173     psFree (FWHM_MAJOR);
    174     psFree (FWHM_MINOR);
    175     psFree (M2);
    176     psFree (M2c);
    177     psFree (M2s);
    178     psFree (M3);
    179     psFree (M4);
    180     psFree (stats);
     269    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4",   PS_META_REPLACE,
     270                     "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4);
     271    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_ER", PS_META_REPLACE,
     272                     "Stdev of M_4", dM4);
     273    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_LQ", PS_META_REPLACE,
     274                     "Lower Quartile of M_4", stats->robustLQ);
     275    psMetadataAddF32(recipe, PS_LIST_TAIL, "IQ_M4_UQ", PS_META_REPLACE,
     276                     "Upper Quartile of M_4", stats->robustUQ);
     277#endif
     278
     279    psLogMsg ("psphot", PS_LOG_INFO,
     280              "Image Quality Stats from %ld psf stars : FWHM (major, minor) :"
     281              "%f, %f  M_2 : %f +/- %f, M_3 : %f +/- %f, M_4 : %f +/- %f\n",
     282              M2->n, fwhm_major, fwhm_minor, vM2, dM2, vM3, dM3, vM4, dM4);
     283
     284    psFree(FWHM_MAJOR);
     285    psFree(FWHM_MINOR);
     286    psFree(M2);
     287    psFree(M2c);
     288    psFree(M2s);
     289    psFree(M3);
     290    psFree(M4);
     291    psFree(stats);
    181292
    182293    return true;
     294
     295 FAIL:
     296    psFree(FWHM_MAJOR);
     297    psFree(FWHM_MINOR);
     298    psFree(M2);
     299    psFree(M2c);
     300    psFree(M2s);
     301    psFree(M3);
     302    psFree(M4);
     303    psFree(stats);
     304
     305    return false;
    183306}

  • trunk/psphot/src/psphotReadout.c

    r19911 r20170  
    6868        return psphotReadoutCleanup (config, readout, recipe, detections, psf, NULL);
    6969    }
    70  
     70
    7171    // construct sources and measure basic stats
    7272    psArray *sources = psphotSourceStats (readout, recipe, detections);
    … …  
    9595    }
    9696
    97     psphotImageQuality (recipe, sources);
     97    if (!psphotImageQuality (recipe, sources)) {
     98        psLogMsg("psphot", 3, "failed to measure image quality");
     99        return psphotReadoutCleanup(config, readout, recipe, detections, psf, sources);
     100    }
    98101
    99102    // if we were not supplied a PSF, choose one here
Note: See TracChangeset for help on using the changeset viewer.