Index: trunk/psphot/src/psphotImageQuality.c =================================================================== --- trunk/psphot/src/psphotImageQuality.c (revision 19907) +++ trunk/psphot/src/psphotImageQuality.c (revision 19951) @@ -1,10 +1,17 @@ # include "psphotInternal.h" +# define USE_SAMPLE 1 /* sample vs robust stats */ // selecting the 'good' stars (likely to be psf stars), measure the M_cn, M_sn terms for n = 2,3,4 bool psphotImageQuality (psMetadata *recipe, psArray *sources) { + psVector *FWHM_MAJOR = psVectorAllocEmpty (100, PS_TYPE_F32); + psVector *FWHM_MINOR = psVectorAllocEmpty (100, PS_TYPE_F32); + psVector *M2 = psVectorAllocEmpty (100, PS_TYPE_F32); psVector *M3 = psVectorAllocEmpty (100, PS_TYPE_F32); psVector *M4 = psVectorAllocEmpty (100, PS_TYPE_F32); + + psVector *M2c = psVectorAllocEmpty (100, PS_TYPE_F32); + psVector *M2s = psVectorAllocEmpty (100, PS_TYPE_F32); for (int i = 0; i < sources->n; i++) { @@ -36,8 +43,19 @@ // 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 + psEllipseMoments emoments; + emoments.x2 = moments->Mxx; + emoments.xy = moments->Mxy; + emoments.y2 = moments->Myy; + psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0); + + psVectorAppend (FWHM_MAJOR, 2.355*axes.major); + psVectorAppend (FWHM_MINOR, 2.355*axes.minor); + // the moments->Mnnn terms contain the straight sums: for example: moments->Mxxx contains \sum x^3 / r float M_c2 = moments->Mxx - moments->Myy; float M_s2 = 2*moments->Mxy; psVectorAppend (M2, hypot(M_c2, M_s2)); + psVectorAppend (M2c, M_c2); + psVectorAppend (M2s, M_s2); float M_c3 = moments->Mxxx - 3.0*moments->Mxyy; @@ -50,11 +68,87 @@ } - psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psMetadataAddS32 (recipe, PS_LIST_TAIL, "IQ_NSTAR", PS_META_REPLACE, "Number of stars used for IQ measurements", M2->n); + +# if (USE_SAMPLE) + psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV | PS_STAT_SAMPLE_QUARTILE); + + psVectorStats (stats, FWHM_MAJOR, NULL, NULL, 0); + float fwhm_major = stats->sampleMean; + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW1", PS_META_REPLACE, "FWHM of Major Axis from moments", stats->sampleMean); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW1_E", PS_META_REPLACE, "FWHM scatter (Major) from moments", stats->sampleStdev); + + psVectorStats (stats, FWHM_MINOR, NULL, NULL, 0); + float fwhm_minor = stats->sampleMean; + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW2", PS_META_REPLACE, "FWHM of Minor Axis from moments", stats->sampleMean); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW2_E", PS_META_REPLACE, "FWHM scatter (Minor) from moments", stats->sampleStdev); + + psVectorStats (stats, M2, NULL, NULL, 0); + float vM2 = stats->sampleMean; + float dM2 = stats->sampleStdev; + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2", PS_META_REPLACE, "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2_ER", PS_META_REPLACE, "Stdev of M_2", dM2); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2_LQ", PS_META_REPLACE, "Lower Quartile of M_2", stats->sampleLQ); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2_UQ", PS_META_REPLACE, "Upper Quartile of M_2", stats->sampleUQ); + + psVectorStats (stats, M2c, NULL, NULL, 0); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C", PS_META_REPLACE, "M_2c = sum f r^2 cos(2phi) / sum f", stats->sampleMean); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C_E", PS_META_REPLACE, "Stdev of M_2c", stats->sampleStdev); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C_L", PS_META_REPLACE, "Lower Quartile of M_2c", stats->sampleLQ); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C_U", PS_META_REPLACE, "Upper Quartile of M_2c", stats->sampleUQ); + + psVectorStats (stats, M2s, NULL, NULL, 0); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S", PS_META_REPLACE, "M_2s = sum f r^2 cos(2phi) / sum f", stats->sampleMean); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S_E", PS_META_REPLACE, "Stdev of M_2s", stats->sampleStdev); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S_L", PS_META_REPLACE, "Lower Quartile of M_2s", stats->sampleLQ); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S_U", PS_META_REPLACE, "Upper Quartile of M_2s", stats->sampleUQ); + + psVectorStats (stats, M3, NULL, NULL, 0); + float vM3 = stats->sampleMean; + float dM3 = stats->sampleStdev; + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3", PS_META_REPLACE, "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3_ER", PS_META_REPLACE, "Stdev of M_3", dM3); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3_LQ", PS_META_REPLACE, "Lower Quartile of M_3", stats->sampleLQ); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3_UQ", PS_META_REPLACE, "Upper Quartile of M_3", stats->sampleUQ); + + psVectorStats (stats, M4, NULL, NULL, 0); + float vM4 = stats->sampleMean; + float dM4 = stats->sampleStdev; + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4", PS_META_REPLACE, "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4_ER", PS_META_REPLACE, "Stdev of M_4", dM4); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4_LQ", PS_META_REPLACE, "Lower Quartile of M_4", stats->sampleLQ); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4_UQ", PS_META_REPLACE, "Upper Quartile of M_4", stats->sampleUQ); + +# else + psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV | PS_STAT_ROBUST_QUARTILE); + + psVectorStats (stats, FWHM_MAJOR, NULL, NULL, 0); + float fwhm_major = stats->robustMedian; + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW1", PS_META_REPLACE, "FWHM of Major Axis from moments", stats->robustMedian); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW1_E", PS_META_REPLACE, "FWHM scatter (Major) from moments", stats->robustStdev); + + psVectorStats (stats, FWHM_MINOR, NULL, NULL, 0); + float fwhm_minor = stats->robustMedian; + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW2", PS_META_REPLACE, "FWHM of Minor Axis from moments", stats->robustMedian); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_FW2_E", PS_META_REPLACE, "FWHM scatter (Minor) from moments", stats->robustStdev); psVectorStats (stats, M2, NULL, NULL, 0); float vM2 = stats->robustMedian; float dM2 = stats->robustStdev; - psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2", PS_META_REPLACE, "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_D_M2", PS_META_REPLACE, "Stdev of M_2", dM2); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2", PS_META_REPLACE, "M_2 = sqrt (M_c2^2 + M_s2^2)", vM2); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2_ER", PS_META_REPLACE, "Stdev of M_2", dM2); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2_LQ", PS_META_REPLACE, "Lower Quartile of M_2", stats->robustLQ); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2_UQ", PS_META_REPLACE, "Upper Quartile of M_2", stats->robustUQ); + + psVectorStats (stats, M2c, NULL, NULL, 0); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C", PS_META_REPLACE, "M_2c = sum f r^2 cos(2phi) / sum f", stats->robustMedian); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C_E", PS_META_REPLACE, "Stdev of M_2c", stats->robustStdev); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C_L", PS_META_REPLACE, "Lower Quartile of M_2c", stats->robustLQ); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2C_U", PS_META_REPLACE, "Upper Quartile of M_2c", stats->robustUQ); + + psVectorStats (stats, M2s, NULL, NULL, 0); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S", PS_META_REPLACE, "M_2s = sum f r^2 cos(2phi) / sum f", stats->robustMedian); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S_E", PS_META_REPLACE, "Stdev of M_2s", stats->robustStdev); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S_L", PS_META_REPLACE, "Lower Quartile of M_2s", stats->robustLQ); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M2S_U", PS_META_REPLACE, "Upper Quartile of M_2s", stats->robustUQ); psVectorStats (stats, M3, NULL, NULL, 0); @@ -62,5 +156,7 @@ float dM3 = stats->robustStdev; psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3", PS_META_REPLACE, "M_3 = sqrt (M_c3^2 + M_s3^2)", vM3); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_D_M3", PS_META_REPLACE, "Stdev of M_3", dM3); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3_ER", PS_META_REPLACE, "Stdev of M_3", dM3); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3_LQ", PS_META_REPLACE, "Lower Quartile of M_3", stats->robustLQ); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M3_UQ", PS_META_REPLACE, "Upper Quartile of M_3", stats->robustUQ); psVectorStats (stats, M4, NULL, NULL, 0); @@ -68,12 +164,20 @@ float dM4 = stats->robustStdev; psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4", PS_META_REPLACE, "M_4 = sqrt (M_c4^2 + M_s4^2)", vM4); - psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_D_M4", PS_META_REPLACE, "Stdev of M_4", dM4); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4_ER", PS_META_REPLACE, "Stdev of M_4", dM4); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4_LQ", PS_META_REPLACE, "Lower Quartile of M_4", stats->robustLQ); + psMetadataAddF32 (recipe, PS_LIST_TAIL, "IQ_M4_UQ", PS_META_REPLACE, "Upper Quartile of M_4", stats->robustUQ); +# endif + 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); + + psFree (FWHM_MAJOR); + psFree (FWHM_MINOR); psFree (M2); + psFree (M2c); + psFree (M2s); psFree (M3); psFree (M4); psFree (stats); - 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); return true; }