Changeset 6182
- Timestamp:
- Jan 22, 2006, 10:20:47 AM (20 years ago)
- Location:
- branches/eam_rel9_p0/psModules/src/objects
- Files:
-
- 2 edited
-
pmObjects.c (modified) (4 diffs)
-
pmObjects.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_rel9_p0/psModules/src/objects/pmObjects.c
r5992 r6182 6 6 * @author EAM, IfA: significant modifications. 7 7 * 8 * @version $Revision: 1.5.4. 2$ $Name: not supported by cvs2svn $9 * @date $Date: 2006-01- 15 18:21:49$8 * @version $Revision: 1.5.4.3 $ $Name: not supported by cvs2svn $ 9 * @date $Date: 2006-01-22 20:20:47 $ 10 10 * 11 11 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 844 844 psF32 xPeak = source->peak->x; 845 845 psF32 yPeak = source->peak->y; 846 psF32 xOff = source->pixels->col0 - source->peak->x; 847 psF32 yOff = source->pixels->row0 - source->peak->y; 846 848 847 849 // XXX why do I get different results for these two methods of finding Sx? … … 854 856 // XXX EAM : mask == 0 is valid 855 857 858 psF32 **vPix = source->pixels->data.F32; 859 psF32 **vWgt = source->weight->data.F32; 860 psU8 **vMsk = source->mask->data.U8; 861 856 862 for (psS32 row = 0; row < source->pixels->numRows ; row++) { 857 863 for (psS32 col = 0; col < source->pixels->numCols ; col++) { 858 if ((source->mask != NULL) && ( source->mask->data.U8[row][col])) {864 if ((source->mask != NULL) && (vMsk[row][col])) { 859 865 continue; 860 866 } 861 867 862 psF32 xDiff = col + source->pixels->col0 - xPeak; 863 psF32 yDiff = row + source->pixels->row0 - yPeak; 868 // psF32 xDiff = col + source->pixels->col0 - xPeak; 869 // psF32 yDiff = row + source->pixels->row0 - yPeak; 870 871 psF32 xDiff = col + xOff; 872 psF32 yDiff = row + yOff; 864 873 865 874 // XXX EAM : calculate xDiff, yDiff up front; … … 869 878 } 870 879 871 psF32 pDiff = source->pixels->data.F32[row][col] - sky;872 psF32 wDiff = source->weight->data.F32[row][col];880 psF32 pDiff = vPix[row][col] - sky; 881 psF32 wDiff = vWgt[row][col]; 873 882 874 883 // XXX EAM : check for valid S/N in pixel 875 884 // XXX EAM : should this limit be user-defined? 876 if (pDiff / sqrt(source->weight->data.F32[row][col]) < 1) { 885 if (PS_SQR(pDiff) < wDiff) { 886 887 Var += wDiff; 888 Sum += pDiff; 889 890 psF32 xWght = xDiff * pDiff; 891 psF32 yWght = yDiff * pDiff; 892 893 X1 += xWght; 894 Y1 += yWght; 895 XY += xDiff * yWght; 896 897 X2 += xDiff * xWght; 898 Y2 += yDiff * yWght; 899 900 peakPixel = PS_MAX (vPix[row][col], peakPixel); 901 numPixels++; 902 } 903 } 904 905 // if we have less than (1/4) of the possible pixels, force a retry 906 // XXX EAM - the limit is a bit arbitrary. make it user defined? 907 if ((numPixels < 0.75*R2) || (Sum <= 0)) { 908 psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n"); 909 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 910 return (false); 911 } 912 913 psTrace (".psModules.pmSourceMoments", 5, 914 "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", 915 sky, Sum, X1, Y1, X2, Y2, XY, numPixels); 916 917 // 918 // first moment X = X1/Sum + xc 919 // second moment X = sqrt (X2/Sum - (X1/Sum)^2) 920 // Sxy = XY / Sum 921 // 922 x = X1/Sum; 923 y = Y1/Sum; 924 if ((fabs(x) > radius) || (fabs(y) > radius)) { 925 psTrace (".psModules.pmSourceMoments", 3, 926 "large centroid swing; invalid peak %d, %d\n", 927 source->peak->x, source->peak->y); 928 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 929 return (false); 930 } 931 932 source->moments->x = x + xPeak; 933 source->moments->y = y + yPeak; 934 935 // XXX EAM : Sxy needs to have x*y subtracted 936 source->moments->Sxy = XY/Sum - x*y; 937 source->moments->Sum = Sum; 938 source->moments->SN = Sum / sqrt(Var); 939 source->moments->Peak = peakPixel; 940 source->moments->nPixels = numPixels; 941 942 // XXX EAM : these values can be negative, so we need to limit the range 943 source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0)); 944 source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0)); 945 946 psTrace (".psModules.pmSourceMoments", 4, 947 "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", 948 sky, Sum, source->moments->x, source->moments->y, 949 source->moments->Sx, source->moments->Sy, source->moments->Sxy); 950 951 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); 952 return(true); 953 } 954 955 // XXX EAM : I used 956 int pmComparePeakAscend (const void **a, const void **b) { 957 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 958 pmPeak *A = *(pmPeak **)a; 959 pmPeak *B = *(pmPeak **)b; 960 961 psF32 diff; 962 963 diff = A->counts - B->counts; 964 if (diff < FLT_EPSILON) { 965 psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); 966 return (-1); 967 } else if (diff > FLT_EPSILON) { 968 psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); 969 return (+1); 970 } 971 psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); 972 return (0); 973 } 974 975 int pmComparePeakDescend (const void **a, const void **b) { 976 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 977 pmPeak *A = *(pmPeak **)a; 978 pmPeak *B = *(pmPeak **)b; 979 980 psF32 diff; 981 982 diff = A->counts - B->counts; 983 if (diff < FLT_EPSILON) { 984 psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); 985 return (+1); 986 } else if (diff > FLT_EPSILON) { 987 psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); 988 return (-1); 989 } 990 psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); 991 return (0); 992 } 993 994 /****************************************************************************** 995 pmSourcePSFClump(source, metadata): Find the likely PSF clump in the 996 sigma-x, sigma-y plane. return 0,0 clump in case of error. 997 *****************************************************************************/ 998 999 // XXX EAM include a S/N cutoff in selecting the sources? 1000 pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata) { 1001 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1002 1003 # define NPIX 10 1004 # define SCALE 0.1 1005 1006 psArray *peaks = NULL; 1007 pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0}; 1008 pmPSFClump psfClump = emptyClump; 1009 1010 PS_ASSERT_PTR_NON_NULL(sources, emptyClump); 1011 PS_ASSERT_PTR_NON_NULL(metadata, emptyClump); 1012 1013 // find the sigmaX, sigmaY clump 1014 { 1015 psStats *stats = NULL; 1016 psImage *splane = NULL; 1017 int binX, binY; 1018 bool status; 1019 1020 psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX"); 1021 if (!status) 1022 SX_MAX = 10.0; 1023 psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX"); 1024 if (!status) 1025 SY_MAX = 10.0; 1026 1027 // construct a sigma-plane image 1028 // psImageAlloc does zero the data 1029 splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32); 1030 for (int i = 0; i < splane->numRows; i++) 1031 { 1032 memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32)); 1033 } 1034 1035 // place the sources in the sigma-plane image (ignore 0,0 values?) 1036 for (psS32 i = 0 ; i < sources->n ; i++) 1037 { 1038 pmSource *tmpSrc = (pmSource *) sources->data[i]; 1039 if (tmpSrc == NULL) { 1040 continue; 1041 } 1042 if (tmpSrc->moments == NULL) { 1043 continue; 1044 } 1045 1046 // Sx,Sy are limited at 0. a peak at 0,0 is artificial 1047 if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) { 1048 continue; 1049 } 1050 1051 // for the moment, force splane dimensions to be 10x10 image pix 1052 binX = tmpSrc->moments->Sx/SCALE; 1053 if (binX < 0) 1054 continue; 1055 if (binX >= splane->numCols) 1056 continue; 1057 1058 binY = tmpSrc->moments->Sy/SCALE; 1059 if (binY < 0) 1060 continue; 1061 if (binY >= splane->numRows) 1062 continue; 1063 1064 splane->data.F32[binY][binX] += 1.0; 1065 } 1066 1067 // find the peak in this image 1068 stats = psStatsAlloc (PS_STAT_MAX); 1069 stats = psImageStats (stats, splane, NULL, 0); 1070 peaks = pmFindImagePeaks (splane, stats[0].max / 2); 1071 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2); 1072 1073 psFree (splane); 1074 psFree (stats); 1075 1076 } 1077 // XXX EAM : possible errors: 1078 // 1) no peak in splane 1079 // 2) no significant peak in splane 1080 1081 // measure statistics on Sx, Sy if Sx, Sy within range of clump 1082 { 1083 pmPeak *clump; 1084 psF32 minSx, maxSx; 1085 psF32 minSy, maxSy; 1086 psVector *tmpSx = NULL; 1087 psVector *tmpSy = NULL; 1088 psStats *stats = NULL; 1089 1090 // XXX EAM : this lets us takes the single highest peak 1091 psArraySort (peaks, pmComparePeakDescend); 1092 clump = peaks->data[0]; 1093 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts); 1094 1095 // define section window for clump 1096 minSx = clump->x * SCALE - 0.2; 1097 maxSx = clump->x * SCALE + 0.2; 1098 minSy = clump->y * SCALE - 0.2; 1099 maxSy = clump->y * SCALE + 0.2; 1100 1101 tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32); 1102 tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32); 1103 tmpSx->n = 0; 1104 tmpSy->n = 0; 1105 1106 // XXX clip sources based on flux? 1107 // create vectors with Sx, Sy values in window 1108 for (psS32 i = 0 ; i < sources->n ; i++) 1109 { 1110 pmSource *tmpSrc = (pmSource *) sources->data[i]; 1111 1112 if (tmpSrc->moments->Sx < minSx) 1113 continue; 1114 if (tmpSrc->moments->Sx > maxSx) 1115 continue; 1116 if (tmpSrc->moments->Sy < minSy) 1117 continue; 1118 if (tmpSrc->moments->Sy > maxSy) 1119 continue; 1120 tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx; 1121 tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy; 1122 tmpSx->n++; 1123 tmpSy->n++; 1124 if (tmpSx->n == tmpSx->nalloc) { 1125 psVectorRealloc (tmpSx, tmpSx->nalloc + 100); 1126 psVectorRealloc (tmpSy, tmpSy->nalloc + 100); 1127 } 1128 } 1129 1130 // measures stats of Sx, Sy 1131 stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); 1132 1133 stats = psVectorStats (stats, tmpSx, NULL, NULL, 0); 1134 psfClump.X = stats->clippedMean; 1135 psfClump.dX = stats->clippedStdev; 1136 1137 stats = psVectorStats (stats, tmpSy, NULL, NULL, 0); 1138 psfClump.Y = stats->clippedMean; 1139 psfClump.dY = stats->clippedStdev; 1140 1141 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); 1142 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY); 1143 // these values should be pushed on the metadata somewhere 1144 1145 psFree (stats); 1146 psFree (peaks); 1147 psFree (tmpSx); 1148 psFree (tmpSy); 1149 } 1150 1151 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1152 return (psfClump); 1153 } 1154 1155 /****************************************************************************** 1156 pmSourceRoughClass(source, metadata): make a guess at the source 1157 classification. 1158 1159 XXX: push the clump info into the metadata? 1160 1161 XXX: How can this function ever return FALSE? 1162 1163 EAM: I moved S/N calculation to pmSourceMoments, using weight image 1164 *****************************************************************************/ 1165 1166 bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump) { 1167 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1168 1169 psBool rc = true; 1170 1171 int Nsat = 0; 1172 int Next = 0; 1173 int Nstar = 0; 1174 int Npsf = 0; 1175 int Ncr = 0; 1176 int Nsatstar = 0; 1177 // psRegion allArray = psRegionSet (0, 0, 0, 0); 1178 psRegion inner; 1179 1180 // report stats on S/N values for star-like objects 1181 psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); 1182 starsn->n = 0; 1183 1184 // check return status value (do these exist?) 1185 bool status; 1186 psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); 1187 1188 // XXX allow clump size to be scaled relative to sigmas? 1189 // make rough IDs based on clumpX,Y,DX,DY 1190 for (psS32 i = 0 ; i < sources->n ; i++) { 1191 1192 pmSource *tmpSrc = (pmSource *) sources->data[i]; 1193 1194 tmpSrc->peak->class = 0; 1195 1196 psF32 sigX = tmpSrc->moments->Sx; 1197 psF32 sigY = tmpSrc->moments->Sy; 1198 1199 // XXX EAM : can we use the value of SATURATE if mask is NULL? 1200 inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2); 1201 int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED); 1202 1203 // saturated star (size consistent with PSF or larger) 1204 // Nsigma should be user-configured parameter 1205 bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY)); 1206 big = true; 1207 if ((Nsatpix > 1) && big) { 1208 tmpSrc->type = PM_SOURCE_STAR; 1209 tmpSrc->mode = PM_SOURCE_SATSTAR; 1210 Nsatstar ++; 877 1211 continue; 878 1212 } 879 1213 880 Var += wDiff; 881 Sum += pDiff; 882 X1 += xDiff * pDiff; 883 Y1 += yDiff * pDiff; 884 XY += xDiff * yDiff * pDiff; 885 886 X2 += PS_SQR(xDiff) * pDiff; 887 Y2 += PS_SQR(yDiff) * pDiff; 888 889 peakPixel = PS_MAX (source->pixels->data.F32[row][col], peakPixel); 890 numPixels++; 891 } 892 } 893 894 // if we have less than (1/4) of the possible pixels, force a retry 895 // XXX EAM - the limit is a bit arbitrary. make it user defined? 896 if ((numPixels < 0.75*R2) || (Sum <= 0)) { 897 psTrace (".psModules.pmSourceMoments", 3, "no valid pixels for source\n"); 898 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 899 return (false); 900 } 901 902 psTrace (".psModules.pmSourceMoments", 5, 903 "sky: %f Sum: %f X1: %f Y1: %f X2: %f Y2: %f XY: %f Npix: %d\n", 904 sky, Sum, X1, Y1, X2, Y2, XY, numPixels); 905 906 // 907 // first moment X = X1/Sum + xc 908 // second moment X = sqrt (X2/Sum - (X1/Sum)^2) 909 // Sxy = XY / Sum 910 // 911 x = X1/Sum; 912 y = Y1/Sum; 913 if ((fabs(x) > radius) || (fabs(y) > radius)) { 914 psTrace (".psModules.pmSourceMoments", 3, 915 "large centroid swing; invalid peak %d, %d\n", 916 source->peak->x, source->peak->y); 917 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 918 return (false); 919 } 920 921 source->moments->x = x + xPeak; 922 source->moments->y = y + yPeak; 923 924 // XXX EAM : Sxy needs to have x*y subtracted 925 source->moments->Sxy = XY/Sum - x*y; 926 source->moments->Sum = Sum; 927 source->moments->SN = Sum / sqrt(Var); 928 source->moments->Peak = peakPixel; 929 source->moments->nPixels = numPixels; 930 931 // XXX EAM : these values can be negative, so we need to limit the range 932 source->moments->Sx = sqrt(PS_MAX(X2/Sum - PS_SQR(x), 0)); 933 source->moments->Sy = sqrt(PS_MAX(Y2/Sum - PS_SQR(y), 0)); 934 935 psTrace (".psModules.pmSourceMoments", 4, 936 "sky: %f Sum: %f x: %f y: %f Sx: %f Sy: %f Sxy: %f\n", 937 sky, Sum, source->moments->x, source->moments->y, 938 source->moments->Sx, source->moments->Sy, source->moments->Sxy); 939 940 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); 941 return(true); 942 } 943 944 // XXX EAM : I used 945 int pmComparePeakAscend (const void **a, const void **b) 946 { 947 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 948 pmPeak *A = *(pmPeak **)a; 949 pmPeak *B = *(pmPeak **)b; 950 951 psF32 diff; 952 953 diff = A->counts - B->counts; 954 if (diff < FLT_EPSILON) { 955 psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); 956 return (-1); 957 } else if (diff > FLT_EPSILON) { 958 psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); 959 return (+1); 960 } 961 psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); 962 return (0); 963 } 964 965 int pmComparePeakDescend (const void **a, const void **b) 966 { 967 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 968 pmPeak *A = *(pmPeak **)a; 969 pmPeak *B = *(pmPeak **)b; 970 971 psF32 diff; 972 973 diff = A->counts - B->counts; 974 if (diff < FLT_EPSILON) { 975 psTrace(__func__, 3, "---- %s(+1) end ----\n", __func__); 976 return (+1); 977 } else if (diff > FLT_EPSILON) { 978 psTrace(__func__, 3, "---- %s(-1) end ----\n", __func__); 979 return (-1); 980 } 981 psTrace(__func__, 3, "---- %s(0) end ----\n", __func__); 982 return (0); 983 } 984 985 /****************************************************************************** 986 pmSourcePSFClump(source, metadata): Find the likely PSF clump in the 987 sigma-x, sigma-y plane. return 0,0 clump in case of error. 988 *****************************************************************************/ 989 990 // XXX EAM include a S/N cutoff in selecting the sources? 991 pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *metadata) 992 { 993 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 994 995 # define NPIX 10 996 # define SCALE 0.1 997 998 psArray *peaks = NULL; 999 pmPSFClump emptyClump = {0.0, 0.0, 0.0, 0.0}; 1000 pmPSFClump psfClump = emptyClump; 1001 1002 PS_ASSERT_PTR_NON_NULL(sources, emptyClump); 1003 PS_ASSERT_PTR_NON_NULL(metadata, emptyClump); 1004 1005 // find the sigmaX, sigmaY clump 1214 // saturated object (not a star, eg bleed trails, hot pixels) 1215 if (Nsatpix > 1) { 1216 tmpSrc->type = PM_SOURCE_SATURATED; 1217 tmpSrc->mode = PM_SOURCE_DEFAULT; 1218 Nsat ++; 1219 continue; 1220 } 1221 1222 // likely defect (too small to be stellar) (push out to 3 sigma) 1223 // low S/N objects which are small are probably stellar 1224 // only set candidate defects if 1225 if ((sigX < 0.05) || (sigY < 0.05)) { 1226 tmpSrc->type = PM_SOURCE_DEFECT; 1227 tmpSrc->mode = PM_SOURCE_DEFAULT; 1228 Ncr ++; 1229 continue; 1230 } 1231 1232 // likely unsaturated extended source (too large to be stellar) 1233 if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { 1234 tmpSrc->type = PM_SOURCE_EXTENDED; 1235 tmpSrc->mode = PM_SOURCE_DEFAULT; 1236 Next ++; 1237 continue; 1238 } 1239 1240 // the rest are probable stellar objects 1241 starsn->data.F32[starsn->n] = tmpSrc->moments->SN; 1242 starsn->n ++; 1243 Nstar ++; 1244 1245 // PSF star (within 1.5 sigma of clump center, S/N > limit) 1246 psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); 1247 if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) { 1248 tmpSrc->type = PM_SOURCE_STAR; 1249 tmpSrc->mode = PM_SOURCE_PSFSTAR; 1250 Npsf ++; 1251 continue; 1252 } 1253 1254 // random type of star 1255 tmpSrc->type = PM_SOURCE_STAR; 1256 tmpSrc->mode = PM_SOURCE_DEFAULT; 1257 } 1258 1259 { 1260 psStats *stats = NULL; 1261 stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); 1262 stats = psVectorStats (stats, starsn, NULL, NULL, 0); 1263 psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max); 1264 psFree (stats); 1265 psFree (starsn); 1266 } 1267 1268 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); 1269 psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); 1270 psTrace (".pmObjects.pmSourceRoughClass", 2, "Next: %3d\n", Next); 1271 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar); 1272 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); 1273 psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); 1274 1275 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); 1276 return(rc); 1277 } 1278 1279 /** pmSourceDefinePixels() 1280 * 1281 * Define psImage subarrays for the source located at coordinates x,y on the 1282 * image set defined by readout. The pixels defined by this operation consist of 1283 * a square window (of full width 2Radius+1) centered on the pixel which contains 1284 * the given coordinate, in the frame of the readout. The window is defined to 1285 * have limits which are valid within the boundary of the readout image, thus if 1286 * the radius would fall outside the image pixels, the subimage is truncated to 1287 * only consist of valid pixels. If readout->mask or readout->weight are not 1288 * NULL, matching subimages are defined for those images as well. This function 1289 * fails if no valid pixels can be defined (x or y less than Radius, for 1290 * example). This function should be used to define a region of interest around a 1291 * source, including both source and sky pixels. 1292 * 1293 * XXX: must code this. 1294 * 1295 */ 1296 bool pmSourceDefinePixels( 1297 pmSource *mySource, ///< Add comment. 1298 pmReadout *readout, ///< Add comment. 1299 psF32 x, ///< Add comment. 1300 psF32 y, ///< Add comment. 1301 psF32 Radius) ///< Add comment. 1006 1302 { 1007 psStats *stats = NULL; 1008 psImage *splane = NULL; 1009 int binX, binY; 1010 bool status; 1011 1012 psF32 SX_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SX_MAX"); 1013 if (!status) 1014 SX_MAX = 10.0; 1015 psF32 SY_MAX = psMetadataLookupF32 (&status, metadata, "MOMENTS_SY_MAX"); 1016 if (!status) 1017 SY_MAX = 10.0; 1018 1019 // construct a sigma-plane image 1020 // psImageAlloc does zero the data 1021 splane = psImageAlloc (SX_MAX/SCALE, SY_MAX/SCALE, PS_TYPE_F32); 1022 for (int i = 0; i < splane->numRows; i++) 1023 { 1024 memset (splane->data.F32[i], 0, splane->numCols*sizeof(PS_TYPE_F32)); 1025 } 1026 1027 // place the sources in the sigma-plane image (ignore 0,0 values?) 1028 for (psS32 i = 0 ; i < sources->n ; i++) 1029 { 1030 pmSource *tmpSrc = (pmSource *) sources->data[i]; 1031 if (tmpSrc == NULL) { 1303 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1304 psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented. Returning FALSE.\n"); 1305 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1306 return(false); 1307 } 1308 1309 /****************************************************************************** 1310 pmSourceSetPixelsCircle(source, image, radius) 1311 1312 XXX: This was replaced by DefinePixels in SDRS. Remove it. 1313 *****************************************************************************/ 1314 bool pmSourceSetPixelsCircle(pmSource *source, 1315 const psImage *image, 1316 psF32 radius) { 1317 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1318 PS_ASSERT_IMAGE_NON_NULL(image, false); 1319 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); 1320 PS_ASSERT_PTR_NON_NULL(source, false); 1321 PS_ASSERT_PTR_NON_NULL(source->moments, false); 1322 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1323 PS_FLOAT_COMPARE(0.0, radius, false); 1324 1325 // 1326 // We define variables for code readability. 1327 // 1328 // XXX: Since the peak->xy coords are in image, not subImage coords, 1329 // these variables should be renamed for clarity (imageCenterRow, etc). 1330 // 1331 psS32 radiusS32 = (psS32) radius; 1332 psS32 SubImageCenterRow = source->peak->y; 1333 psS32 SubImageCenterCol = source->peak->x; 1334 // XXX EAM : for the circle to stay on the image 1335 // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1) 1336 psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - radiusS32); 1337 psS32 SubImageEndRow = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1); 1338 psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32); 1339 psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1); 1340 1341 // XXX: Must recycle image. 1342 // XXX EAM: this message reflects a programming error we know about. 1343 // i am setting it to a trace message which we can take out 1344 if (source->pixels != NULL) { 1345 psTrace (".psModule.pmObjects.pmSourceSetPixelsCircle", 4, 1346 "WARNING: pmSourceSetPixelsCircle(): image->pixels not NULL. Freeing and reallocating.\n"); 1347 psFree(source->pixels); 1348 } 1349 source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol, 1350 SubImageStartRow, 1351 SubImageEndCol, 1352 SubImageEndRow)); 1353 1354 // XXX: Must recycle image. 1355 if (source->mask != NULL) { 1356 psFree(source->mask); 1357 } 1358 source->mask = psImageAlloc(source->pixels->numCols, 1359 source->pixels->numRows, 1360 PS_TYPE_U8); // XXX EAM : type was F32 1361 1362 // 1363 // Loop through the subimage mask, initialize mask to 0 or 1. 1364 // XXX EAM: valid pixels should have 0, not 1 1365 for (psS32 row = 0 ; row < source->mask->numRows; row++) { 1366 for (psS32 col = 0 ; col < source->mask->numCols; col++) { 1367 1368 if (checkRadius2((psF32) radiusS32, 1369 (psF32) radiusS32, 1370 radius, 1371 (psF32) col, 1372 (psF32) row)) { 1373 source->mask->data.U8[row][col] = 0; 1374 } else { 1375 source->mask->data.U8[row][col] = 1; 1376 } 1377 } 1378 } 1379 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); 1380 return(true); 1381 } 1382 1383 /****************************************************************************** 1384 pmSourceModelGuess(source, model): This function allocates a new 1385 pmModel structure based on the given modelType specified in the argument list. 1386 The corresponding pmModelGuess function is returned, and used to 1387 supply the values of the params array in the pmModel structure. 1388 1389 XXX: Many parameters are based on the src->moments structure, which is in 1390 image, not subImage coords. Therefore, the calls to the model evaluation 1391 functions will be in image, not subImage coords. Remember this. 1392 *****************************************************************************/ 1393 pmModel *pmSourceModelGuess(pmSource *source, 1394 pmModelType modelType) { 1395 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1396 PS_ASSERT_PTR_NON_NULL(source->moments, false); 1397 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1398 1399 pmModel *model = pmModelAlloc(modelType); 1400 1401 pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType); 1402 modelGuessFunc(model, source); 1403 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1404 return(model); 1405 } 1406 1407 /****************************************************************************** 1408 evalModel(source, level, row): a private function which evaluates the 1409 source->modelPSF function at the specified coords. The coords are subImage, not 1410 image coords. 1411 1412 NOTE: The coords are in subImage source->pixel coords, not image coords. 1413 1414 XXX: reverse order of row,col args? 1415 1416 XXX: rename all coords in this file such that their name defines whether 1417 the coords is in subImage or image space. 1418 1419 XXX: This should probably be a public pmModules function. 1420 1421 XXX: Use static vectors for x. 1422 1423 XXX: Figure out if it's (row, col) or (col, row) for the model functions. 1424 1425 XXX: For a while, the first psVectorAlloc() was generating a seg fault during 1426 testing. Try to reproduce that and debug. 1427 *****************************************************************************/ 1428 1429 // XXX EAM : I have made this a public function 1430 // XXX EAM : this now uses a pmModel as the input 1431 // XXX EAM : it was using src->type to find the model, not model->type 1432 psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row) { 1433 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1434 PS_ASSERT_PTR_NON_NULL(image, false); 1435 PS_ASSERT_PTR_NON_NULL(model, false); 1436 PS_ASSERT_PTR_NON_NULL(model->params, false); 1437 1438 // Allocate the x coordinate structure and convert row/col to image space. 1439 // 1440 psVector *x = psVectorAlloc(2, PS_TYPE_F32); 1441 x->data.F32[0] = (psF32) (col + image->col0); 1442 x->data.F32[1] = (psF32) (row + image->row0); 1443 psF32 tmpF; 1444 pmModelFunc modelFunc; 1445 1446 modelFunc = pmModelFunc_GetFunction (model->type); 1447 tmpF = modelFunc (NULL, model->params, x); 1448 psFree(x); 1449 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1450 return(tmpF); 1451 } 1452 1453 /****************************************************************************** 1454 pmSourceContour(src, img, level, mode): For an input subImage, and model, this 1455 routine returns a psArray of coordinates that evaluate to the specified level. 1456 1457 XXX: Probably should remove the "image" argument. 1458 XXX: What type should the output coordinate vectors consist of? col,row? 1459 XXX: Why a pmArray output? 1460 XXX: doex x,y correspond with col,row or row/col? 1461 XXX: What is mode? 1462 XXX: The top, bottom of the contour is not correctly determined. 1463 XXX EAM : this function is using the model for the contour, but it should 1464 be using only the image counts 1465 *****************************************************************************/ 1466 psArray *pmSourceContour(pmSource *source, 1467 const psImage *image, 1468 psF32 level, 1469 pmContourType mode) { 1470 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1471 PS_ASSERT_PTR_NON_NULL(source, false); 1472 PS_ASSERT_PTR_NON_NULL(image, false); 1473 PS_ASSERT_PTR_NON_NULL(source->moments, false); 1474 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1475 PS_ASSERT_PTR_NON_NULL(source->pixels, false); 1476 PS_ASSERT_PTR_NON_NULL(source->modelEXT, false); 1477 // XXX EAM : what is the purpose of modelPSF/modelEXT? 1478 1479 // 1480 // Allocate data for x/y pairs. 1481 // 1482 psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); 1483 psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); 1484 1485 // 1486 // Start at the row with peak pixel, then decrement. 1487 // 1488 psS32 col = source->peak->x; 1489 for (psS32 row = source->peak->y; row>= 0 ; row--) { 1490 // XXX: yVec contain no real information. Do we really need it? 1491 yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); 1492 yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); 1493 1494 // Starting at peak pixel, search leftwards for the column intercept. 1495 psF32 leftIntercept = findValue(source, level, row, col, 0); 1496 if (isnan(leftIntercept)) { 1497 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); 1498 psFree(xVec); 1499 psFree(yVec); 1500 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); 1501 return(NULL); 1502 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); 1503 } 1504 xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; 1505 1506 // Starting at peak pixel, search rightwards for the column intercept. 1507 1508 psF32 rightIntercept = findValue(source, level, row, col, 1); 1509 if (isnan(rightIntercept)) { 1510 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); 1511 psFree(xVec); 1512 psFree(yVec); 1513 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); 1514 return(NULL); 1515 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); 1516 } 1517 psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept); 1518 xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; 1519 1520 // Set starting column for next row 1521 col = (psS32) ((leftIntercept + rightIntercept) / 2.0); 1522 } 1523 // 1524 // Start at the row (+1) with peak pixel, then increment. 1525 // 1526 col = source->peak->x; 1527 for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) { 1528 // XXX: yVec contain no real information. Do we really need it? 1529 yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); 1530 yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); 1531 1532 // Starting at peak pixel, search leftwards for the column intercept. 1533 psF32 leftIntercept = findValue(source, level, row, col, 0); 1534 if (isnan(leftIntercept)) { 1535 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); 1536 psFree(xVec); 1537 psFree(yVec); 1538 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); 1539 return(NULL); 1540 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); 1541 } 1542 xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; 1543 1544 // Starting at peak pixel, search rightwards for the column intercept. 1545 psF32 rightIntercept = findValue(source, level, row, col, 1); 1546 if (isnan(rightIntercept)) { 1547 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); 1548 psFree(xVec); 1549 psFree(yVec); 1550 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); 1551 return(NULL); 1552 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); 1553 } 1554 xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; 1555 1556 // Set starting column for next row 1557 col = (psS32) ((leftIntercept + rightIntercept) / 2.0); 1558 } 1559 1560 // 1561 // Allocate an array for result, store coord vectors there. 1562 // 1563 psArray *tmpArray = psArrayAlloc(2); 1564 tmpArray->data[0] = (psPtr *) yVec; 1565 tmpArray->data[1] = (psPtr *) xVec; 1566 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1567 return(tmpArray); 1568 } 1569 1570 // save a static values so they may be set externally 1571 static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15; 1572 static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1; 1573 1574 bool pmSourceFitModelInit (float nIter, float tol) { 1575 1576 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = nIter; 1577 PM_SOURCE_FIT_MODEL_TOLERANCE = tol; 1578 } 1579 1580 bool pmSourceFitModel (pmSource *source, 1581 pmModel *model, 1582 const bool PSF) { 1583 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1584 PS_ASSERT_PTR_NON_NULL(source, false); 1585 PS_ASSERT_PTR_NON_NULL(source->moments, false); 1586 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1587 PS_ASSERT_PTR_NON_NULL(source->pixels, false); 1588 PS_ASSERT_PTR_NON_NULL(source->mask, false); 1589 PS_ASSERT_PTR_NON_NULL(source->weight, false); 1590 1591 // XXX EAM : is it necessary for the mask & weight to exist? the 1592 // tests below could be conditions (!NULL) 1593 1594 psBool fitStatus = true; 1595 psBool onPic = true; 1596 psBool rc = true; 1597 1598 psVector *params = model->params; 1599 psVector *dparams = model->dparams; 1600 psVector *paramMask = NULL; 1601 1602 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); 1603 1604 int nParams = PSF ? 4 : params->n; 1605 1606 // maximum number of valid pixels 1607 psS32 nPix = source->pixels->numRows * source->pixels->numCols; 1608 1609 // construct the coordinate and value entries 1610 psArray *x = psArrayAlloc(nPix); 1611 psVector *y = psVectorAlloc(nPix, PS_TYPE_F32); 1612 psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32); 1613 1614 nPix = 0; 1615 for (psS32 i = 0; i < source->pixels->numRows; i++) { 1616 for (psS32 j = 0; j < source->pixels->numCols; j++) { 1617 if (source->mask->data.U8[i][j]) { 1618 continue; 1619 } 1620 psVector *coord = psVectorAlloc(2, PS_TYPE_F32); 1621 1622 // Convert i/j to image space: 1623 coord->data.F32[0] = (psF32) (j + source->pixels->col0); 1624 coord->data.F32[1] = (psF32) (i + source->pixels->row0); 1625 x->data[nPix] = (psPtr *) coord; 1626 y->data.F32[nPix] = source->pixels->data.F32[i][j]; 1627 1628 // psMinimizeLMChi2 takes wt = 1/dY^2 1629 if (source->weight->data.F32[i][j] == 0) { 1630 continue; 1631 } 1632 yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; 1633 nPix++; 1634 } 1635 } 1636 if (nPix < nParams + 1) { 1637 psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); 1638 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 1639 model->status = PM_MODEL_BADARGS; 1640 psFree (x); 1641 psFree (y); 1642 psFree (yErr); 1643 return(false); 1644 } 1645 x->n = nPix; 1646 y->n = nPix; 1647 yErr->n = nPix; 1648 1649 psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, 1650 PM_SOURCE_FIT_MODEL_TOLERANCE); 1651 1652 // PSF model only fits first 4 parameters, EXT model fits all 1653 if (PSF) { 1654 paramMask = psVectorAlloc (params->n, PS_TYPE_U8); 1655 for (int i = 0; i < 4; i++) { 1656 paramMask->data.U8[i] = 0; 1657 } 1658 for (int i = 4; i < paramMask->n; i++) { 1659 paramMask->data.U8[i] = 1; 1660 } 1661 } 1662 1663 // Set the parameter range checks 1664 pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); 1665 psVector *beta_lim = NULL; 1666 psVector *params_min = NULL; 1667 psVector *params_max = NULL; 1668 1669 // XXX EAM : in this implementation, I pass in the limits with the covar matrix. 1670 // in the SDRS, I define a new psMinimization which will take these in 1671 psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64); 1672 modelLimits (&beta_lim, ¶ms_min, ¶ms_max); 1673 for (int i = 0; i < params->n; i++) { 1674 covar->data.F64[0][i] = beta_lim->data.F32[i]; 1675 covar->data.F64[1][i] = params_min->data.F32[i]; 1676 covar->data.F64[2][i] = params_max->data.F32[i]; 1677 } 1678 1679 psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); 1680 fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc); 1681 for (int i = 0; i < dparams->n; i++) { 1682 if ((paramMask != NULL) && paramMask->data.U8[i]) 1032 1683 continue; 1033 } 1034 if (tmpSrc->moments == NULL) { 1035 continue; 1036 } 1037 1038 // Sx,Sy are limited at 0. a peak at 0,0 is artificial 1039 if ((fabs(tmpSrc->moments->Sx) < FLT_EPSILON) && (fabs(tmpSrc->moments->Sy) < FLT_EPSILON)) { 1040 continue; 1041 } 1042 1043 // for the moment, force splane dimensions to be 10x10 image pix 1044 binX = tmpSrc->moments->Sx/SCALE; 1045 if (binX < 0) 1046 continue; 1047 if (binX >= splane->numCols) 1048 continue; 1049 1050 binY = tmpSrc->moments->Sy/SCALE; 1051 if (binY < 0) 1052 continue; 1053 if (binY >= splane->numRows) 1054 continue; 1055 1056 splane->data.F32[binY][binX] += 1.0; 1057 } 1058 1059 // find the peak in this image 1060 stats = psStatsAlloc (PS_STAT_MAX); 1061 stats = psImageStats (stats, splane, NULL, 0); 1062 peaks = pmFindImagePeaks (splane, stats[0].max / 2); 1063 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump threshold is %f\n", stats[0].max/2); 1064 1065 psFree (splane); 1066 psFree (stats); 1067 1068 } 1069 // XXX EAM : possible errors: 1070 // 1) no peak in splane 1071 // 2) no significant peak in splane 1072 1073 // measure statistics on Sx, Sy if Sx, Sy within range of clump 1074 { 1075 pmPeak *clump; 1076 psF32 minSx, maxSx; 1077 psF32 minSy, maxSy; 1078 psVector *tmpSx = NULL; 1079 psVector *tmpSy = NULL; 1080 psStats *stats = NULL; 1081 1082 // XXX EAM : this lets us takes the single highest peak 1083 psArraySort (peaks, pmComparePeakDescend); 1084 clump = peaks->data[0]; 1085 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump is at %d, %d (%f)\n", clump->x, clump->y, clump->counts); 1086 1087 // define section window for clump 1088 minSx = clump->x * SCALE - 0.2; 1089 maxSx = clump->x * SCALE + 0.2; 1090 minSy = clump->y * SCALE - 0.2; 1091 maxSy = clump->y * SCALE + 0.2; 1092 1093 tmpSx = psVectorAlloc (sources->n, PS_TYPE_F32); 1094 tmpSy = psVectorAlloc (sources->n, PS_TYPE_F32); 1095 tmpSx->n = 0; 1096 tmpSy->n = 0; 1097 1098 // XXX clip sources based on flux? 1099 // create vectors with Sx, Sy values in window 1100 for (psS32 i = 0 ; i < sources->n ; i++) 1101 { 1102 pmSource *tmpSrc = (pmSource *) sources->data[i]; 1103 1104 if (tmpSrc->moments->Sx < minSx) 1105 continue; 1106 if (tmpSrc->moments->Sx > maxSx) 1107 continue; 1108 if (tmpSrc->moments->Sy < minSy) 1109 continue; 1110 if (tmpSrc->moments->Sy > maxSy) 1111 continue; 1112 tmpSx->data.F32[tmpSx->n] = tmpSrc->moments->Sx; 1113 tmpSy->data.F32[tmpSy->n] = tmpSrc->moments->Sy; 1114 tmpSx->n++; 1115 tmpSy->n++; 1116 if (tmpSx->n == tmpSx->nalloc) { 1117 psVectorRealloc (tmpSx, tmpSx->nalloc + 100); 1118 psVectorRealloc (tmpSy, tmpSy->nalloc + 100); 1119 } 1120 } 1121 1122 // measures stats of Sx, Sy 1123 stats = psStatsAlloc (PS_STAT_CLIPPED_MEAN | PS_STAT_CLIPPED_STDEV); 1124 1125 stats = psVectorStats (stats, tmpSx, NULL, NULL, 0); 1126 psfClump.X = stats->clippedMean; 1127 psfClump.dX = stats->clippedStdev; 1128 1129 stats = psVectorStats (stats, tmpSy, NULL, NULL, 0); 1130 psfClump.Y = stats->clippedMean; 1131 psfClump.dY = stats->clippedStdev; 1132 1133 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump X, Y: %f, %f\n", psfClump.X, psfClump.Y); 1134 psTrace (".pmObjects.pmSourceRoughClass", 2, "clump DX, DY: %f, %f\n", psfClump.dX, psfClump.dY); 1135 // these values should be pushed on the metadata somewhere 1136 1137 psFree (stats); 1138 psFree (peaks); 1139 psFree (tmpSx); 1140 psFree (tmpSy); 1141 } 1142 1143 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1144 return (psfClump); 1145 } 1146 1147 /****************************************************************************** 1148 pmSourceRoughClass(source, metadata): make a guess at the source 1149 classification. 1150 1151 XXX: push the clump info into the metadata? 1152 1153 XXX: How can this function ever return FALSE? 1154 1155 EAM: I moved S/N calculation to pmSourceMoments, using weight image 1156 *****************************************************************************/ 1157 1158 bool pmSourceRoughClass(psArray *sources, psMetadata *metadata, pmPSFClump clump) 1159 { 1160 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1161 1162 psBool rc = true; 1163 1164 int Nsat = 0; 1165 int Next = 0; 1166 int Nstar = 0; 1167 int Npsf = 0; 1168 int Ncr = 0; 1169 int Nsatstar = 0; 1170 // psRegion allArray = psRegionSet (0, 0, 0, 0); 1171 psRegion inner; 1172 1173 // report stats on S/N values for star-like objects 1174 psVector *starsn = psVectorAlloc (sources->n, PS_TYPE_F32); 1175 starsn->n = 0; 1176 1177 // check return status value (do these exist?) 1178 bool status; 1179 psF32 PSF_SN_LIM = psMetadataLookupF32 (&status, metadata, "PSF_SN_LIM"); 1180 1181 // XXX allow clump size to be scaled relative to sigmas? 1182 // make rough IDs based on clumpX,Y,DX,DY 1183 for (psS32 i = 0 ; i < sources->n ; i++) { 1184 1185 pmSource *tmpSrc = (pmSource *) sources->data[i]; 1186 1187 tmpSrc->peak->class = 0; 1188 1189 psF32 sigX = tmpSrc->moments->Sx; 1190 psF32 sigY = tmpSrc->moments->Sy; 1191 1192 // XXX EAM : can we use the value of SATURATE if mask is NULL? 1193 inner = psRegionForSquare (tmpSrc->peak->x - tmpSrc->mask->col0, tmpSrc->peak->y - tmpSrc->mask->row0, 2); 1194 int Nsatpix = psImageCountPixelMask (tmpSrc->mask, inner, PSPHOT_MASK_SATURATED); 1195 1196 // saturated star (size consistent with PSF or larger) 1197 // Nsigma should be user-configured parameter 1198 bool big = (sigX > (clump.X - clump.dX)) && (sigY > (clump.Y - clump.dY)); 1199 big = true; 1200 if ((Nsatpix > 1) && big) { 1201 tmpSrc->type = PM_SOURCE_STAR; 1202 tmpSrc->mode = PM_SOURCE_SATSTAR; 1203 Nsatstar ++; 1204 continue; 1205 } 1206 1207 // saturated object (not a star, eg bleed trails, hot pixels) 1208 if (Nsatpix > 1) { 1209 tmpSrc->type = PM_SOURCE_SATURATED; 1210 tmpSrc->mode = PM_SOURCE_DEFAULT; 1211 Nsat ++; 1212 continue; 1213 } 1214 1215 // likely defect (too small to be stellar) (push out to 3 sigma) 1216 // low S/N objects which are small are probably stellar 1217 // only set candidate defects if 1218 if ((sigX < 0.05) || (sigY < 0.05)) { 1219 tmpSrc->type = PM_SOURCE_DEFECT; 1220 tmpSrc->mode = PM_SOURCE_DEFAULT; 1221 Ncr ++; 1222 continue; 1223 } 1224 1225 // likely unsaturated extended source (too large to be stellar) 1226 if ((sigX > (clump.X + 3*clump.dX)) || (sigY > (clump.Y + 3*clump.dY))) { 1227 tmpSrc->type = PM_SOURCE_EXTENDED; 1228 tmpSrc->mode = PM_SOURCE_DEFAULT; 1229 Next ++; 1230 continue; 1231 } 1232 1233 // the rest are probable stellar objects 1234 starsn->data.F32[starsn->n] = tmpSrc->moments->SN; 1235 starsn->n ++; 1236 Nstar ++; 1237 1238 // PSF star (within 1.5 sigma of clump center, S/N > limit) 1239 psF32 radius = hypot ((sigX-clump.X)/clump.dX, (sigY-clump.Y)/clump.dY); 1240 if ((tmpSrc->moments->SN > PSF_SN_LIM) && (radius < 1.5)) { 1241 tmpSrc->type = PM_SOURCE_STAR; 1242 tmpSrc->mode = PM_SOURCE_PSFSTAR; 1243 Npsf ++; 1244 continue; 1245 } 1246 1247 // random type of star 1248 tmpSrc->type = PM_SOURCE_STAR; 1249 tmpSrc->mode = PM_SOURCE_DEFAULT; 1250 } 1251 1252 { 1253 psStats *stats = NULL; 1254 stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX); 1255 stats = psVectorStats (stats, starsn, NULL, NULL, 0); 1256 psLogMsg ("pmObjects", 3, "SN range: %f - %f\n", stats[0].min, stats[0].max); 1257 psFree (stats); 1258 psFree (starsn); 1259 } 1260 1261 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nstar: %3d\n", Nstar); 1262 psTrace (".pmObjects.pmSourceRoughClass", 2, "Npsf: %3d\n", Npsf); 1263 psTrace (".pmObjects.pmSourceRoughClass", 2, "Next: %3d\n", Next); 1264 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsatstar: %3d\n", Nsatstar); 1265 psTrace (".pmObjects.pmSourceRoughClass", 2, "Nsat: %3d\n", Nsat); 1266 psTrace (".pmObjects.pmSourceRoughClass", 2, "Ncr: %3d\n", Ncr); 1267 1268 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); 1269 return(rc); 1270 } 1271 1272 /** pmSourceDefinePixels() 1273 * 1274 * Define psImage subarrays for the source located at coordinates x,y on the 1275 * image set defined by readout. The pixels defined by this operation consist of 1276 * a square window (of full width 2Radius+1) centered on the pixel which contains 1277 * the given coordinate, in the frame of the readout. The window is defined to 1278 * have limits which are valid within the boundary of the readout image, thus if 1279 * the radius would fall outside the image pixels, the subimage is truncated to 1280 * only consist of valid pixels. If readout->mask or readout->weight are not 1281 * NULL, matching subimages are defined for those images as well. This function 1282 * fails if no valid pixels can be defined (x or y less than Radius, for 1283 * example). This function should be used to define a region of interest around a 1284 * source, including both source and sky pixels. 1285 * 1286 * XXX: must code this. 1287 * 1288 */ 1289 bool pmSourceDefinePixels( 1290 pmSource *mySource, ///< Add comment. 1291 pmReadout *readout, ///< Add comment. 1292 psF32 x, ///< Add comment. 1293 psF32 y, ///< Add comment. 1294 psF32 Radius) ///< Add comment. 1295 { 1296 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1297 psLogMsg(__func__, PS_LOG_WARN, "WARNING: pmSourceDefinePixels() has not been implemented. Returning FALSE.\n"); 1298 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1299 return(false); 1300 } 1301 1302 /****************************************************************************** 1303 pmSourceSetPixelsCircle(source, image, radius) 1304 1305 XXX: This was replaced by DefinePixels in SDRS. Remove it. 1306 *****************************************************************************/ 1307 bool pmSourceSetPixelsCircle(pmSource *source, 1308 const psImage *image, 1309 psF32 radius) 1310 { 1311 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1312 PS_ASSERT_IMAGE_NON_NULL(image, false); 1313 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); 1314 PS_ASSERT_PTR_NON_NULL(source, false); 1315 PS_ASSERT_PTR_NON_NULL(source->moments, false); 1316 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1317 PS_FLOAT_COMPARE(0.0, radius, false); 1318 1319 // 1320 // We define variables for code readability. 1321 // 1322 // XXX: Since the peak->xy coords are in image, not subImage coords, 1323 // these variables should be renamed for clarity (imageCenterRow, etc). 1324 // 1325 psS32 radiusS32 = (psS32) radius; 1326 psS32 SubImageCenterRow = source->peak->y; 1327 psS32 SubImageCenterCol = source->peak->x; 1328 // XXX EAM : for the circle to stay on the image 1329 // XXX EAM : EndRow is *exclusive* of pixel region (ie, last pixel + 1) 1330 psS32 SubImageStartRow = PS_MAX (0, SubImageCenterRow - radiusS32); 1331 psS32 SubImageEndRow = PS_MIN (image->numRows, SubImageCenterRow + radiusS32 + 1); 1332 psS32 SubImageStartCol = PS_MAX (0, SubImageCenterCol - radiusS32); 1333 psS32 SubImageEndCol = PS_MIN (image->numCols, SubImageCenterCol + radiusS32 + 1); 1334 1335 // XXX: Must recycle image. 1336 // XXX EAM: this message reflects a programming error we know about. 1337 // i am setting it to a trace message which we can take out 1338 if (source->pixels != NULL) { 1339 psTrace (".psModule.pmObjects.pmSourceSetPixelsCircle", 4, 1340 "WARNING: pmSourceSetPixelsCircle(): image->pixels not NULL. Freeing and reallocating.\n"); 1341 psFree(source->pixels); 1342 } 1343 source->pixels = psImageSubset((psImage *) image, psRegionSet(SubImageStartCol, 1344 SubImageStartRow, 1345 SubImageEndCol, 1346 SubImageEndRow)); 1347 1348 // XXX: Must recycle image. 1349 if (source->mask != NULL) { 1350 psFree(source->mask); 1351 } 1352 source->mask = psImageAlloc(source->pixels->numCols, 1353 source->pixels->numRows, 1354 PS_TYPE_U8); // XXX EAM : type was F32 1355 1356 // 1357 // Loop through the subimage mask, initialize mask to 0 or 1. 1358 // XXX EAM: valid pixels should have 0, not 1 1359 for (psS32 row = 0 ; row < source->mask->numRows; row++) { 1360 for (psS32 col = 0 ; col < source->mask->numCols; col++) { 1361 1362 if (checkRadius2((psF32) radiusS32, 1363 (psF32) radiusS32, 1364 radius, 1365 (psF32) col, 1366 (psF32) row)) { 1367 source->mask->data.U8[row][col] = 0; 1368 } else { 1369 source->mask->data.U8[row][col] = 1; 1370 } 1371 } 1372 } 1373 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); 1374 return(true); 1375 } 1376 1377 /****************************************************************************** 1378 pmSourceModelGuess(source, model): This function allocates a new 1379 pmModel structure based on the given modelType specified in the argument list. 1380 The corresponding pmModelGuess function is returned, and used to 1381 supply the values of the params array in the pmModel structure. 1382 1383 XXX: Many parameters are based on the src->moments structure, which is in 1384 image, not subImage coords. Therefore, the calls to the model evaluation 1385 functions will be in image, not subImage coords. Remember this. 1386 *****************************************************************************/ 1387 pmModel *pmSourceModelGuess(pmSource *source, 1388 pmModelType modelType) 1389 { 1390 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1391 PS_ASSERT_PTR_NON_NULL(source->moments, false); 1392 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1393 1394 pmModel *model = pmModelAlloc(modelType); 1395 1396 pmModelGuessFunc modelGuessFunc = pmModelGuessFunc_GetFunction(modelType); 1397 modelGuessFunc(model, source); 1398 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1399 return(model); 1400 } 1401 1402 /****************************************************************************** 1403 evalModel(source, level, row): a private function which evaluates the 1404 source->modelPSF function at the specified coords. The coords are subImage, not 1405 image coords. 1406 1407 NOTE: The coords are in subImage source->pixel coords, not image coords. 1408 1409 XXX: reverse order of row,col args? 1410 1411 XXX: rename all coords in this file such that their name defines whether 1412 the coords is in subImage or image space. 1413 1414 XXX: This should probably be a public pmModules function. 1415 1416 XXX: Use static vectors for x. 1417 1418 XXX: Figure out if it's (row, col) or (col, row) for the model functions. 1419 1420 XXX: For a while, the first psVectorAlloc() was generating a seg fault during 1421 testing. Try to reproduce that and debug. 1422 *****************************************************************************/ 1423 1424 // XXX EAM : I have made this a public function 1425 // XXX EAM : this now uses a pmModel as the input 1426 // XXX EAM : it was using src->type to find the model, not model->type 1427 psF32 pmModelEval(pmModel *model, psImage *image, psS32 col, psS32 row) 1428 { 1429 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1430 PS_ASSERT_PTR_NON_NULL(image, false); 1431 PS_ASSERT_PTR_NON_NULL(model, false); 1432 PS_ASSERT_PTR_NON_NULL(model->params, false); 1433 1434 // Allocate the x coordinate structure and convert row/col to image space. 1435 // 1436 psVector *x = psVectorAlloc(2, PS_TYPE_F32); 1437 x->data.F32[0] = (psF32) (col + image->col0); 1438 x->data.F32[1] = (psF32) (row + image->row0); 1439 psF32 tmpF; 1440 pmModelFunc modelFunc; 1441 1442 modelFunc = pmModelFunc_GetFunction (model->type); 1443 tmpF = modelFunc (NULL, model->params, x); 1444 psFree(x); 1445 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1446 return(tmpF); 1447 } 1448 1449 /****************************************************************************** 1450 pmSourceContour(src, img, level, mode): For an input subImage, and model, this 1451 routine returns a psArray of coordinates that evaluate to the specified level. 1452 1453 XXX: Probably should remove the "image" argument. 1454 XXX: What type should the output coordinate vectors consist of? col,row? 1455 XXX: Why a pmArray output? 1456 XXX: doex x,y correspond with col,row or row/col? 1457 XXX: What is mode? 1458 XXX: The top, bottom of the contour is not correctly determined. 1459 XXX EAM : this function is using the model for the contour, but it should 1460 be using only the image counts 1461 *****************************************************************************/ 1462 psArray *pmSourceContour(pmSource *source, 1463 const psImage *image, 1464 psF32 level, 1465 pmContourType mode) 1466 { 1467 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1468 PS_ASSERT_PTR_NON_NULL(source, false); 1469 PS_ASSERT_PTR_NON_NULL(image, false); 1470 PS_ASSERT_PTR_NON_NULL(source->moments, false); 1471 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1472 PS_ASSERT_PTR_NON_NULL(source->pixels, false); 1473 PS_ASSERT_PTR_NON_NULL(source->modelEXT, false); 1474 // XXX EAM : what is the purpose of modelPSF/modelEXT? 1475 1476 // 1477 // Allocate data for x/y pairs. 1478 // 1479 psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); 1480 psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); 1481 1482 // 1483 // Start at the row with peak pixel, then decrement. 1484 // 1485 psS32 col = source->peak->x; 1486 for (psS32 row = source->peak->y; row>= 0 ; row--) { 1487 // XXX: yVec contain no real information. Do we really need it? 1488 yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); 1489 yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); 1490 1491 // Starting at peak pixel, search leftwards for the column intercept. 1492 psF32 leftIntercept = findValue(source, level, row, col, 0); 1493 if (isnan(leftIntercept)) { 1494 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); 1495 psFree(xVec); 1496 psFree(yVec); 1497 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); 1498 return(NULL); 1499 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); 1500 } 1501 xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; 1502 1503 // Starting at peak pixel, search rightwards for the column intercept. 1504 1505 psF32 rightIntercept = findValue(source, level, row, col, 1); 1506 if (isnan(rightIntercept)) { 1507 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); 1508 psFree(xVec); 1509 psFree(yVec); 1510 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); 1511 return(NULL); 1512 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); 1513 } 1514 psTrace(__func__, 4, "The intercepts are (%.2f, %.2f)\n", leftIntercept, rightIntercept); 1515 xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; 1516 1517 // Set starting column for next row 1518 col = (psS32) ((leftIntercept + rightIntercept) / 2.0); 1519 } 1520 // 1521 // Start at the row (+1) with peak pixel, then increment. 1522 // 1523 col = source->peak->x; 1524 for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) { 1525 // XXX: yVec contain no real information. Do we really need it? 1526 yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); 1527 yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); 1528 1529 // Starting at peak pixel, search leftwards for the column intercept. 1530 psF32 leftIntercept = findValue(source, level, row, col, 0); 1531 if (isnan(leftIntercept)) { 1532 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); 1533 psFree(xVec); 1534 psFree(yVec); 1535 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); 1536 return(NULL); 1537 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); 1538 } 1539 xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; 1540 1541 // Starting at peak pixel, search rightwards for the column intercept. 1542 psF32 rightIntercept = findValue(source, level, row, col, 1); 1543 if (isnan(rightIntercept)) { 1544 psError(PS_ERR_UNKNOWN, true, "Could not find contour edge (NAN)"); 1545 psFree(xVec); 1546 psFree(yVec); 1547 psTrace(__func__, 3, "---- %s(NULL) end ----\n", __func__); 1548 return(NULL); 1549 //psLogMsg(__func__, PS_LOG_WARN, "WARNING: Could not find contour edge (NAN)\n"); 1550 } 1551 xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; 1552 1553 // Set starting column for next row 1554 col = (psS32) ((leftIntercept + rightIntercept) / 2.0); 1555 } 1556 1557 // 1558 // Allocate an array for result, store coord vectors there. 1559 // 1560 psArray *tmpArray = psArrayAlloc(2); 1561 tmpArray->data[0] = (psPtr *) yVec; 1562 tmpArray->data[1] = (psPtr *) xVec; 1563 psTrace(__func__, 3, "---- %s() end ----\n", __func__); 1564 return(tmpArray); 1565 } 1566 1567 // XXX EAM : these are better starting values, but should be available from metadata? 1568 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 15 1569 #define PM_SOURCE_FIT_MODEL_TOLERANCE 0.1 1570 bool pmSourceFitModel (pmSource *source, 1571 pmModel *model, 1572 const bool PSF) 1573 { 1574 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1575 PS_ASSERT_PTR_NON_NULL(source, false); 1576 PS_ASSERT_PTR_NON_NULL(source->moments, false); 1577 PS_ASSERT_PTR_NON_NULL(source->peak, false); 1578 PS_ASSERT_PTR_NON_NULL(source->pixels, false); 1579 PS_ASSERT_PTR_NON_NULL(source->mask, false); 1580 PS_ASSERT_PTR_NON_NULL(source->weight, false); 1581 1582 // XXX EAM : is it necessary for the mask & weight to exist? the 1583 // tests below could be conditions (!NULL) 1584 1585 psBool fitStatus = true; 1586 psBool onPic = true; 1587 psBool rc = true; 1588 1589 psVector *params = model->params; 1590 psVector *dparams = model->dparams; 1591 psVector *paramMask = NULL; 1592 1593 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); 1594 1595 int nParams = PSF ? 4 : params->n; 1596 1597 // maximum number of valid pixels 1598 psS32 nPix = source->pixels->numRows * source->pixels->numCols; 1599 1600 // construct the coordinate and value entries 1601 psArray *x = psArrayAlloc(nPix); 1602 psVector *y = psVectorAlloc(nPix, PS_TYPE_F32); 1603 psVector *yErr = psVectorAlloc(nPix, PS_TYPE_F32); 1604 1605 nPix = 0; 1606 for (psS32 i = 0; i < source->pixels->numRows; i++) { 1607 for (psS32 j = 0; j < source->pixels->numCols; j++) { 1608 if (source->mask->data.U8[i][j]) { 1609 continue; 1610 } 1611 psVector *coord = psVectorAlloc(2, PS_TYPE_F32); 1612 1613 // Convert i/j to image space: 1614 coord->data.F32[0] = (psF32) (j + source->pixels->col0); 1615 coord->data.F32[1] = (psF32) (i + source->pixels->row0); 1616 x->data[nPix] = (psPtr *) coord; 1617 y->data.F32[nPix] = source->pixels->data.F32[i][j]; 1618 1619 // psMinimizeLMChi2 takes wt = 1/dY^2 1620 if (source->weight->data.F32[i][j] == 0) { 1621 continue; 1622 } 1623 yErr->data.F32[nPix] = 1.0 / source->weight->data.F32[i][j]; 1624 nPix++; 1625 } 1626 } 1627 if (nPix < nParams + 1) { 1628 psTrace (".pmObjects.pmSourceFitModel", 4, "insufficient valid pixels\n"); 1629 psTrace(__func__, 3, "---- %s(false) end ----\n", __func__); 1630 model->status = PM_MODEL_BADARGS; 1631 psFree (x); 1632 psFree (y); 1633 psFree (yErr); 1634 return(false); 1635 } 1636 x->n = nPix; 1637 y->n = nPix; 1638 yErr->n = nPix; 1639 1640 psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, 1641 PM_SOURCE_FIT_MODEL_TOLERANCE); 1642 1643 // PSF model only fits first 4 parameters, EXT model fits all 1644 if (PSF) { 1645 paramMask = psVectorAlloc (params->n, PS_TYPE_U8); 1646 for (int i = 0; i < 4; i++) { 1647 paramMask->data.U8[i] = 0; 1648 } 1649 for (int i = 4; i < paramMask->n; i++) { 1650 paramMask->data.U8[i] = 1; 1651 } 1652 } 1653 1654 // Set the parameter range checks 1655 pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type); 1656 psVector *beta_lim = NULL; 1657 psVector *params_min = NULL; 1658 psVector *params_max = NULL; 1659 1660 // XXX EAM : in this implementation, I pass in the limits with the covar matrix. 1661 // in the SDRS, I define a new psMinimization which will take these in 1662 psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64); 1663 modelLimits (&beta_lim, ¶ms_min, ¶ms_max); 1664 for (int i = 0; i < params->n; i++) { 1665 covar->data.F64[0][i] = beta_lim->data.F32[i]; 1666 covar->data.F64[1][i] = params_min->data.F32[i]; 1667 covar->data.F64[2][i] = params_max->data.F32[i]; 1668 } 1669 1670 psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n"); 1671 fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc); 1672 for (int i = 0; i < dparams->n; i++) { 1673 if ((paramMask != NULL) && paramMask->data.U8[i]) 1674 continue; 1675 dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); 1676 } 1677 1678 // save the resulting chisq, nDOF, nIter 1679 model->chisq = myMin->value; 1680 model->nIter = myMin->iter; 1681 model->nDOF = y->n - nParams; 1682 1683 // get the Gauss-Newton distance for fixed model parameters 1684 if (paramMask != NULL) { 1685 psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); 1686 psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc); 1687 for (int i = 0; i < dparams->n; i++) { 1688 if (!paramMask->data.U8[i]) 1689 continue; 1690 dparams->data.F32[i] = delta->data.F64[i]; 1691 } 1692 } 1693 1694 // set the model success or failure status 1695 if (!fitStatus) { 1696 model->status = PM_MODEL_NONCONVERGE; 1697 } else { 1698 model->status = PM_MODEL_SUCCESS; 1699 } 1700 1701 // models can go insane: reject these 1702 onPic &= (params->data.F32[2] >= source->pixels->col0); 1703 onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); 1704 onPic &= (params->data.F32[3] >= source->pixels->row0); 1705 onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); 1706 if (!onPic) { 1707 model->status = PM_MODEL_OFFIMAGE; 1708 } 1709 1710 source->mode |= PM_SOURCE_FITTED; 1711 1712 psFree(x); 1713 psFree(y); 1714 psFree(yErr); 1715 psFree(myMin); 1716 psFree(covar); 1717 psFree(paramMask); 1718 1719 rc = (onPic && fitStatus); 1720 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); 1721 return(rc); 1722 } 1723 1724 bool p_pmSourceAddOrSubModel(psImage *image, 1725 psImage *mask, 1726 pmModel *model, 1727 bool center, 1728 bool sky, 1729 bool add 1730 ) 1731 { 1732 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1733 1734 PS_ASSERT_PTR_NON_NULL(model, false); 1735 PS_ASSERT_IMAGE_NON_NULL(image, false); 1736 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); 1737 1738 psVector *x = psVectorAlloc(2, PS_TYPE_F32); 1739 psVector *params = model->params; 1740 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); 1741 psS32 imageCol; 1742 psS32 imageRow; 1743 psF32 skyValue = params->data.F32[0]; 1744 psF32 pixelValue; 1745 1746 for (psS32 i = 0; i < image->numRows; i++) { 1747 for (psS32 j = 0; j < image->numCols; j++) { 1748 if ((mask != NULL) && mask->data.U8[i][j]) 1749 continue; 1750 1751 // XXX: Should you be adding the pixels for the entire subImage, 1752 // or a radius of pixels around it? 1753 1754 // Convert i/j to imace coord space: 1755 // XXX: Make sure you have col/row order correct. 1756 // XXX EAM : 'center' option changes this 1757 // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2] 1758 if (center) { 1759 imageCol = j - 0.5*image->numCols + model->params->data.F32[2]; 1760 imageRow = i - 0.5*image->numRows + model->params->data.F32[3]; 1761 } else { 1762 imageCol = j + image->col0; 1763 imageRow = i + image->row0; 1764 } 1765 1766 x->data.F32[0] = (float) imageCol; 1767 x->data.F32[1] = (float) imageRow; 1768 1769 // set the appropriate pixel value for this coordinate 1770 if (sky) { 1771 pixelValue = modelFunc (NULL, params, x); 1772 } else { 1773 pixelValue = modelFunc (NULL, params, x) - skyValue; 1774 } 1775 1776 1777 // add or subtract the value 1778 if (add 1779 ) { 1780 image->data.F32[i][j] += pixelValue; 1781 } 1782 else { 1783 image->data.F32[i][j] -= pixelValue; 1784 } 1785 } 1786 } 1787 psFree(x); 1788 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); 1789 return(true); 1790 } 1791 1792 1793 1794 /****************************************************************************** 1795 *****************************************************************************/ 1796 bool pmSourceAddModel(psImage *image, 1797 psImage *mask, 1798 pmModel *model, 1799 bool center, 1800 bool sky) 1801 { 1802 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1803 psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true); 1804 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); 1805 return(rc); 1806 } 1807 1808 /****************************************************************************** 1809 *****************************************************************************/ 1810 bool pmSourceSubModel(psImage *image, 1811 psImage *mask, 1812 pmModel *model, 1813 bool center, 1814 bool sky) 1815 { 1816 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1817 psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false); 1818 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); 1819 return(rc); 1820 } 1821 1822 bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) 1823 { 1824 1825 float obsSum = 0; 1826 float fitSum = 0; 1827 float sky = model->params->data.F32[0]; 1828 1829 pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); 1830 fitSum = modelFluxFunc (model->params); 1831 1832 for (int ix = 0; ix < image->numCols; ix++) { 1833 for (int iy = 0; iy < image->numRows; iy++) { 1834 if (mask->data.U8[iy][ix]) 1835 continue; 1836 obsSum += image->data.F32[iy][ix] - sky; 1837 } 1838 } 1839 if (obsSum <= 0) 1840 return false; 1841 if (fitSum <= 0) 1842 return false; 1843 1844 *fitMag = -2.5*log10(fitSum); 1845 *obsMag = -2.5*log10(obsSum); 1846 return (true); 1847 } 1848 1684 dparams->data.F32[i] = sqrt(covar->data.F64[i][i]); 1685 } 1686 1687 // save the resulting chisq, nDOF, nIter 1688 model->chisq = myMin->value; 1689 model->nIter = myMin->iter; 1690 model->nDOF = y->n - nParams; 1691 1692 // get the Gauss-Newton distance for fixed model parameters 1693 if (paramMask != NULL) { 1694 psVector *delta = psVectorAlloc (params->n, PS_TYPE_F64); 1695 psMinimizeGaussNewtonDelta(delta, params, NULL, x, y, yErr, modelFunc); 1696 for (int i = 0; i < dparams->n; i++) { 1697 if (!paramMask->data.U8[i]) 1698 continue; 1699 dparams->data.F32[i] = delta->data.F64[i]; 1700 } 1701 } 1702 1703 // set the model success or failure status 1704 if (!fitStatus) { 1705 model->status = PM_MODEL_NONCONVERGE; 1706 } else { 1707 model->status = PM_MODEL_SUCCESS; 1708 } 1709 1710 // models can go insane: reject these 1711 onPic &= (params->data.F32[2] >= source->pixels->col0); 1712 onPic &= (params->data.F32[2] < source->pixels->col0 + source->pixels->numCols); 1713 onPic &= (params->data.F32[3] >= source->pixels->row0); 1714 onPic &= (params->data.F32[3] < source->pixels->row0 + source->pixels->numRows); 1715 if (!onPic) { 1716 model->status = PM_MODEL_OFFIMAGE; 1717 } 1718 1719 source->mode |= PM_SOURCE_FITTED; 1720 1721 psFree(x); 1722 psFree(y); 1723 psFree(yErr); 1724 psFree(myMin); 1725 psFree(covar); 1726 psFree(paramMask); 1727 1728 rc = (onPic && fitStatus); 1729 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); 1730 return(rc); 1731 } 1732 1733 bool p_pmSourceAddOrSubModel(psImage *image, 1734 psImage *mask, 1735 pmModel *model, 1736 bool center, 1737 bool sky, 1738 bool add 1739 ) { 1740 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1741 1742 PS_ASSERT_PTR_NON_NULL(model, false); 1743 PS_ASSERT_IMAGE_NON_NULL(image, false); 1744 PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); 1745 1746 psVector *x = psVectorAlloc(2, PS_TYPE_F32); 1747 psVector *params = model->params; 1748 pmModelFunc modelFunc = pmModelFunc_GetFunction (model->type); 1749 psS32 imageCol; 1750 psS32 imageRow; 1751 psF32 skyValue = params->data.F32[0]; 1752 psF32 pixelValue; 1753 1754 for (psS32 i = 0; i < image->numRows; i++) { 1755 for (psS32 j = 0; j < image->numCols; j++) { 1756 if ((mask != NULL) && mask->data.U8[i][j]) 1757 continue; 1758 1759 // XXX: Should you be adding the pixels for the entire subImage, 1760 // or a radius of pixels around it? 1761 1762 // Convert i/j to imace coord space: 1763 // XXX: Make sure you have col/row order correct. 1764 // XXX EAM : 'center' option changes this 1765 // XXX EAM : i == numCols/2 -> x = model->params->data.F32[2] 1766 if (center) { 1767 imageCol = j - 0.5*image->numCols + model->params->data.F32[2]; 1768 imageRow = i - 0.5*image->numRows + model->params->data.F32[3]; 1769 } else { 1770 imageCol = j + image->col0; 1771 imageRow = i + image->row0; 1772 } 1773 1774 x->data.F32[0] = (float) imageCol; 1775 x->data.F32[1] = (float) imageRow; 1776 1777 // set the appropriate pixel value for this coordinate 1778 if (sky) { 1779 pixelValue = modelFunc (NULL, params, x); 1780 } else { 1781 pixelValue = modelFunc (NULL, params, x) - skyValue; 1782 } 1783 1784 1785 // add or subtract the value 1786 if (add 1787 ) { 1788 image->data.F32[i][j] += pixelValue; 1789 } 1790 else { 1791 image->data.F32[i][j] -= pixelValue; 1792 } 1793 } 1794 } 1795 psFree(x); 1796 psTrace(__func__, 3, "---- %s(true) end ----\n", __func__); 1797 return(true); 1798 } 1799 1800 1801 1802 /****************************************************************************** 1803 *****************************************************************************/ 1804 bool pmSourceAddModel(psImage *image, 1805 psImage *mask, 1806 pmModel *model, 1807 bool center, 1808 bool sky) { 1809 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1810 psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, true); 1811 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); 1812 return(rc); 1813 } 1814 1815 /****************************************************************************** 1816 *****************************************************************************/ 1817 bool pmSourceSubModel(psImage *image, 1818 psImage *mask, 1819 pmModel *model, 1820 bool center, 1821 bool sky) { 1822 psTrace(__func__, 3, "---- %s() begin ----\n", __func__); 1823 psBool rc = p_pmSourceAddOrSubModel(image, mask, model, center, sky, false); 1824 psTrace(__func__, 3, "---- %s(%d) end ----\n", __func__, rc); 1825 return(rc); 1826 } 1827 1828 bool pmSourcePhotometry (float *fitMag, float *obsMag, pmModel *model, psImage *image, psImage *mask) { 1829 1830 float obsSum = 0; 1831 float fitSum = 0; 1832 float sky = model->params->data.F32[0]; 1833 1834 pmModelFlux modelFluxFunc = pmModelFlux_GetFunction (model->type); 1835 fitSum = modelFluxFunc (model->params); 1836 1837 for (int ix = 0; ix < image->numCols; ix++) { 1838 for (int iy = 0; iy < image->numRows; iy++) { 1839 if (mask->data.U8[iy][ix]) 1840 continue; 1841 obsSum += image->data.F32[iy][ix] - sky; 1842 } 1843 } 1844 if (obsSum <= 0) 1845 return false; 1846 if (fitSum <= 0) 1847 return false; 1848 1849 *fitMag = -2.5*log10(fitSum); 1850 *obsMag = -2.5*log10(obsSum); 1851 return (true); 1852 } 1853 -
branches/eam_rel9_p0/psModules/src/objects/pmObjects.h
r6062 r6182 10 10 * @author GLG, MHPCC 11 11 * 12 * @version $Revision: 1.4.4. 4$ $Name: not supported by cvs2svn $13 * @date $Date: 2006-01-2 0 02:38:28$12 * @version $Revision: 1.4.4.5 $ $Name: not supported by cvs2svn $ 13 * @date $Date: 2006-01-22 20:20:47 $ 14 14 * 15 15 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 475 475 476 476 477 bool pmSourceFitModelInit( 478 float nIter; ///< max number of allowed iterations 479 float tol ///< convergence criterion 480 ); 481 477 482 /** pmSourceFitModel() 478 483 *
Note:
See TracChangeset
for help on using the changeset viewer.
