IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 19951


Ignore:
Timestamp:
Oct 7, 2008, 10:07:40 AM (18 years ago)
Author:
eugene
Message:

fix names of IQ stats, add M2s, M2c, FWHM major, minor

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/psphot/src/psphotImageQuality.c

    r19907 r19951  
    11# include "psphotInternal.h"
     2# define USE_SAMPLE 1 /* sample vs robust stats */
    23
    34// selecting the 'good' stars (likely to be psf stars), measure the M_cn, M_sn terms for n = 2,3,4
    45bool psphotImageQuality (psMetadata *recipe, psArray *sources) {
    56
     7    psVector *FWHM_MAJOR = psVectorAllocEmpty (100, PS_TYPE_F32);
     8    psVector *FWHM_MINOR = psVectorAllocEmpty (100, PS_TYPE_F32);
     9
    610    psVector *M2 = psVectorAllocEmpty (100, PS_TYPE_F32);
    711    psVector *M3 = psVectorAllocEmpty (100, PS_TYPE_F32);
    812    psVector *M4 = psVectorAllocEmpty (100, PS_TYPE_F32);
     13
     14    psVector *M2c = psVectorAllocEmpty (100, PS_TYPE_F32);
     15    psVector *M2s = psVectorAllocEmpty (100, PS_TYPE_F32);
    916
    1017    for (int i = 0; i < sources->n; i++) {
    … …  
    3643        // 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
    3744
     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
    3854        // the moments->Mnnn terms contain the straight sums: for example: moments->Mxxx contains \sum x^3 / r
    3955        float M_c2 = moments->Mxx - moments->Myy;
    4056        float M_s2 = 2*moments->Mxy;
    4157        psVectorAppend (M2, hypot(M_c2, M_s2));
     58        psVectorAppend (M2c, M_c2);
     59        psVectorAppend (M2s, M_s2);
    4260       
    4361        float M_c3 = moments->Mxxx - 3.0*moments->Mxyy;
    … …  
    5068    }   
    5169
    52     psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     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)
     73    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_QUARTILE);
     74
     75    psVectorStats (stats, FWHM_MAJOR, NULL, NULL, 0);
     76    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);
     81    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);
     86    float vM2 = stats->sampleMean;
     87    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);
     106    float vM3 = stats->sampleMean;
     107    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);
     114    float vM4 = stats->sampleMean;
     115    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
     122    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE);
     123
     124    psVectorStats (stats, FWHM_MAJOR, NULL, NULL, 0);
     125    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);
     130    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);
    53133
    54134    psVectorStats (stats, M2, NULL, NULL, 0);
    55135    float vM2 = stats->robustMedian;
    56136    float dM2 = stats->robustStdev;
    57     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2",   PS_META_REPLACE, "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2);
    58     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_D_M2", PS_META_REPLACE, "Stdev of M_2", dM2);
     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);
    59153
    60154    psVectorStats (stats, M3, NULL, NULL, 0);
    … …  
    62156    float dM3 = stats->robustStdev;
    63157    psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3",   PS_META_REPLACE, "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3);
    64     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_D_M3", PS_META_REPLACE, "Stdev of M_3", dM3);
     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);
    65161
    66162    psVectorStats (stats, M4, NULL, NULL, 0);
    … …  
    68164    float dM4 = stats->robustStdev;
    69165    psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4",   PS_META_REPLACE, "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4);
    70     psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_D_M4", PS_META_REPLACE, "Stdev of M_4", dM4);
     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
    71170
     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);
    72175    psFree (M2);
     176    psFree (M2c);
     177    psFree (M2s);
    73178    psFree (M3);
    74179    psFree (M4);
    75180    psFree (stats);
    76181
    77     psLogMsg ("psphot", PS_LOG_INFO, "Image Quality Stats : M_2 : %f +/- %f, M_3 : %f +/- %f, M_4 : %f +/- %f\n", vM2, dM2, vM3, dM3, vM4, dM4);
    78182    return true;
    79183}
Note: See TracChangeset for help on using the changeset viewer.