Changeset 3511
- Timestamp:
- Mar 25, 2005, 1:22:30 PM (21 years ago)
- File:
-
- 1 edited
-
trunk/psModules/src/pmObjects.c (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/pmObjects.c
r3510 r3511 5 5 * @author GLG, MHPCC 6 6 * 7 * @version $Revision: 1. 7$ $Name: not supported by cvs2svn $8 * @date $Date: 2005-03-25 23: 18:35$7 * @version $Revision: 1.8 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2005-03-25 23:22:30 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 762 762 tmpSrc->peak->class = 0; 763 763 764 psF32 sigX = 0.0; 765 psF32 sigY = 0.0; 766 // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE. 767 psF32 clumpX = 0.0; 768 psF32 clumpDX = 0.0; 769 psF32 clumpY = 0.0; 770 psF32 clumpDY = 0.0; 771 764 772 if (tmpSrc->moments->Peak > SATURATE) { 765 773 tmpSrc->peak->class|= PS_SOURCE_SATURATED; … … 779 787 tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR; 780 788 } 781 // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE.782 psF32 clumpX = 0.0;783 psF32 clumpDX = 0.0;784 psF32 clumpY = 0.0;785 psF32 clumpYX = 0.0;786 789 // XXX: The SDRS is not real clear on how to calculate sigX, sigY. 787 psF32 sigX = 0.0 788 psF32 sigY = 0.0 789 if ((fabs(sigX - clumpX) < clumpDX) && 790 (fabs(sigY - clumpY) < clumpDY)) { 791 tmpSrc->peak->class|= PS_SOURCE_PSFSTAR; 792 } 793 794 if ((sigX < (clumpX - clumpDX)) && 795 (sigY < (clumpY - clumpDY))) 796 tmpSrc->peak->class|= PS_SOURCE_DEFECT; 790 if ((fabs(sigX - clumpX) < clumpDX) && 791 (fabs(sigY - clumpY) < clumpDY)) { 792 tmpSrc->peak->class|= PS_SOURCE_PSFSTAR; 793 } 794 795 if ((sigX < (clumpX - clumpDX)) && 796 (sigY < (clumpY - clumpDY))) 797 tmpSrc->peak->class|= PS_SOURCE_DEFECT; 797 798 } 798 799 799 800 if ((sigX > (clumpX + clumpDX)) && 800 (sigY > (clumpY + clumpDY))) 801 (sigY > (clumpY + clumpDY))) { 801 802 tmpSrc->peak->class|= PS_SOURCE_GALAXY; 802 } 803 804 if (tmpSrc->peak->class == 0) { 805 tmpSrc->peak->class|= PS_SOURCE_OTHER; 806 } 807 } 808 } 809 810 return(rc); 811 } 812 813 814 815 /****************************************************************************** 816 pmSourceSetPixelCircle(source, image, radius) 817 818 XXX: Why boolean output? 819 820 XXX: Why are we checking source->moments for NULL? Should the circle be 821 centered on the centroid or the peak? 822 823 XXX: The circle will have a diameter of (1+radius). This is different from 824 the pmSourceSetLocal() function. 825 *****************************************************************************/ 826 bool pmSourceSetPixelCircle(psSource *source, 827 const psImage *image, 828 psF32 radius) 829 { 830 PS_IMAGE_CHECK_NULL(image, false); 831 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); 832 PS_PTR_CHECK_NULL(source, false); 833 // PS_PTR_CHECK_NULL(source->moments, false); 834 PS_PTR_CHECK_NULL(source->peak, false); 835 PS_FLOAT_COMPARE(0.0, radius, false); 836 837 // 838 // We define variables for code readability. 839 // 840 psS32 radiusS32 = (psS32) radius; 841 psS32 SubImageCenterRow = source->peak->y; 842 psS32 SubImageCenterCol = source->peak->x; 843 psS32 SubImageStartRow = SubImageCenterRow - radiusS32; 844 psS32 SubImageEndRow = SubImageCenterRow + radiusS32; 845 psS32 SubImageStartCol = SubImageCenterCol - radiusS32; 846 psS32 SubImageEndCol = SubImageCenterCol + radiusS32; 847 848 if (SubImageStartRow < 0) { 849 psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", 850 SubImageStartRow); 851 return(false); 852 } 853 if (SubImageEndRow+1 >= image->numRows) { 854 psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", 855 SubImageEndRow); 856 return(false); 857 } 858 if (SubImageStartCol < 0) { 859 psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n", 860 SubImageStartCol); 861 return(false); 862 } 863 if (SubImageEndCol+1 >= image->numCols) { 864 psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", 865 SubImageEndCol); 866 return(false); 867 } 868 869 // XXX: Must recycle image. 870 if (source->pixels != NULL) { 871 psLogMsg(__func__, PS_LOG_WARN, 872 "WARNING: pmSourceSetPixelCircle(): image->pixels not NULL. Freeing and reallocating.\n"); 873 psFree(source->pixels); 874 } 875 source->pixels = psImageSubset((psImage *) image, 876 SubImageStartCol, 877 SubImageStartRow, 878 SubImageEndCol+1, 879 SubImageEndRow+1); 880 881 // XXX: Must recycle image. 882 if (source->mask != NULL) { 883 psFree(source->mask); 884 } 885 source->mask = psImageAlloc(1 + 2 * radiusS32, 886 1 + 2 * radiusS32, 887 PS_TYPE_F32); 888 889 // 890 // Loop through the subimage mask, initialize mask to 0 or 1. 891 // 892 for (psS32 row = 0 ; row < source->mask->numRows; row++) { 893 for (psS32 col = 0 ; col < source->mask->numCols; col++) { 894 895 if (CheckRadius2((psF32) radiusS32, 896 (psF32) radiusS32, 897 radius, 898 (psF32) col, 899 (psF32) row)) { 900 source->mask->data.U8[row][col] = 1; 901 } else { 902 source->mask->data.U8[row][col] = 0; 903 } 904 } 905 } 906 return(true); 907 } 908 909 910 /****************************************************************************** 911 pmSourceModelGuess(source, image, model): This function allocates a new 912 psModel structure and store it in the psSource data structure specified in the 913 argument list. The model type is specified in the argument list. The params 914 array in that psModel structure are allocated, and then set to the appropriate 915 values. This function returns true if everything was successful. 916 917 XXX: Many of the initial parameters are set to 0.0 since I don't know what 918 the appropiate initial guesses are. 919 *****************************************************************************/ 920 bool pmSourceModelGuess(psSource *source, 921 const psImage *image, 922 psModelType model) 923 { 924 PS_PTR_CHECK_NULL(source, false); 925 PS_PTR_CHECK_NULL(source->moments, false); 926 PS_PTR_CHECK_NULL(source->peak, false); 927 PS_IMAGE_CHECK_NULL(image, false); 928 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); 929 if (source->models != NULL) { 930 psLogMsg(__func__, PS_LOG_WARN, "WARNING: source->models was non-NULL; calling psFree(source->models).\n"); 931 psFree(source->models); 932 } 933 source->models = pmModelAlloc(PS_MODEL_UNDEFINED); 934 935 switch (model) { 936 case PS_MODEL_GAUSS: 937 source->models->type = PS_MODEL_GAUSS; 938 source->models->Nparams = 7; 939 source->models->params = (psF32 *) psAlloc(7 * sizeof(psF32)); 940 source->models->dparams = (psF32 *) psAlloc(7 * sizeof(psF32)); 941 for (psS32 i = 0 ; i < 7 ; i++) { 942 source->models->params[i] = 0.0; 943 source->models->dparams[i] = 0.0; 944 } 945 source->models->params[0] = source->moments->Sky; 946 source->models->params[1] = source->peak->counts - source->moments->Sky; 947 source->models->params[2] = source->moments->x; 948 source->models->params[3] = source->moments->y; 949 source->models->params[4] = sqrt(2.0) / source->moments->Sx; 950 source->models->params[5] = sqrt(2.0) / source->moments->Sy; 951 source->models->params[6] = source->moments->Sxy; 952 source->models->chisq = 0.0; 953 return(true); 954 case PS_MODEL_PGAUSS: 955 source->models->type = PS_MODEL_PGAUSS; 956 source->models->Nparams = 7; 957 source->models->params = (psF32 *) psAlloc(7 * sizeof(psF32)); 958 source->models->dparams = (psF32 *) psAlloc(7 * sizeof(psF32)); 959 for (psS32 i = 0 ; i < 7 ; i++) { 960 source->models->params[i] = 0.0; 961 source->models->dparams[i] = 0.0; 962 } 963 source->models->params[0] = source->moments->Sky; 964 source->models->params[1] = source->peak->counts - source->moments->Sky; 965 source->models->params[2] = source->moments->x; 966 source->models->params[3] = source->moments->y; 967 source->models->params[4] = sqrt(2.0) / source->moments->Sx; 968 source->models->params[5] = sqrt(2.0) / source->moments->Sy; 969 source->models->params[6] = source->moments->Sxy; 970 source->models->chisq = 0.0; 971 return(true); 972 case PS_MODEL_TWIST_GAUSS: 973 source->models->type = PS_MODEL_TWIST_GAUSS; 974 source->models->Nparams = 11; 975 source->models->params = (psF32 *) psAlloc(11 * sizeof(psF32)); 976 source->models->dparams = (psF32 *) psAlloc(11 * sizeof(psF32)); 977 for (psS32 i = 0 ; i < 11 ; i++) { 978 source->models->params[i] = 0.0; 979 source->models->dparams[i] = 0.0; 980 } 981 982 source->models->params[0] = source->moments->Sky; 983 source->models->params[1] = source->peak->counts - source->moments->Sky; 984 source->models->params[2] = source->moments->x; 985 source->models->params[3] = source->moments->y; 986 // XXX: What are these? 987 // source->models->params[4] = SxInner; 988 // source->models->params[5] = SyInner; 989 // source->models->params[6] = SxyInner; 990 // source->models->params[7] = SxOuter; 991 // source->models->params[8] = SyOuter; 992 // source->models->params[9] = SxyOuter; 993 // source->models->params[10] = N; 994 995 source->models->chisq = 0.0; 996 return(true); 997 case PS_MODEL_WAUSS: 998 999 source->models->params[0] = source->moments->Sky; 1000 source->models->params[1] = source->peak->counts - source->moments->Sky; 1001 source->models->params[2] = source->moments->x; 1002 source->models->params[3] = source->moments->y; 1003 source->models->params[4] = sqrt(2.0) / source->moments->Sx; 1004 source->models->params[5] = sqrt(2.0) / source->moments->Sy; 1005 source->models->params[6] = source->moments->Sxy; 1006 // XXX: What are these? 1007 // source->models->params[7] = B2; 1008 // source->models->params[8] = B3; 1009 1010 source->models->type = PS_MODEL_WAUSS; 1011 source->models->Nparams = 9; 1012 source->models->params = (psF32 *) psAlloc(9 * sizeof(psF32)); 1013 source->models->dparams = (psF32 *) psAlloc(9 * sizeof(psF32)); 1014 for (psS32 i = 0 ; i < 9 ; i++) { 1015 source->models->params[i] = 0.0; 1016 source->models->dparams[i] = 0.0; 1017 } 1018 source->models->chisq = 0.0; 1019 return(true); 1020 case PS_MODEL_SERSIC: 1021 source->models->type = PS_MODEL_SERSIC; 1022 source->models->Nparams = 8; 1023 source->models->params = (psF32 *) psAlloc(8 * sizeof(psF32)); 1024 source->models->dparams = (psF32 *) psAlloc(8 * sizeof(psF32)); 1025 for (psS32 i = 0 ; i < 8 ; i++) { 1026 source->models->params[i] = 0.0; 1027 source->models->dparams[i] = 0.0; 1028 } 1029 1030 source->models->params[0] = source->moments->Sky; 1031 source->models->params[1] = source->peak->counts - source->moments->Sky; 1032 source->models->params[2] = source->moments->x; 1033 source->models->params[3] = source->moments->y; 1034 source->models->params[4] = sqrt(2.0) / source->moments->Sx; 1035 source->models->params[5] = sqrt(2.0) / source->moments->Sy; 1036 source->models->params[6] = source->moments->Sxy; 1037 // XXX: What are these? 1038 //source->models->params[7] = Nexp; 1039 1040 source->models->chisq = 0.0; 1041 return(true); 1042 case PS_MODEL_SERSIC_CORE: 1043 source->models->type = PS_MODEL_SERSIC_CORE; 1044 source->models->Nparams = 12; 1045 source->models->params = (psF32 *) psAlloc(12 * sizeof(psF32)); 1046 source->models->dparams = (psF32 *) psAlloc(12 * sizeof(psF32)); 1047 for (psS32 i = 0 ; i < 12 ; i++) { 1048 source->models->params[i] = 0.0; 1049 source->models->dparams[i] = 0.0; 1050 } 1051 1052 source->models->params[0] = source->moments->Sky; 1053 source->models->params[1] = source->peak->counts - source->moments->Sky; 1054 source->models->params[2] = source->moments->x; 1055 source->models->params[3] = source->moments->y; 1056 // XXX: What are these? 1057 //source->models->params[4] SxInner; 1058 //source->models->params[5] SyInner; 1059 //source->models->params[6] SxyInner; 1060 //source->models->params[7] Zd; 1061 //source->models->params[8] SxOuter; 1062 //source->models->params[9] SyOuter; 1063 //source->models->params[10] = SxyOuter; 1064 //source->models->params[11] = Nexp; 1065 1066 source->models->chisq = 0.0; 1067 return(true); 1068 default: 1069 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1070 return(false); 1071 } 1072 } 1073 1074 /****************************************************************************** 1075 evalModel(source, level, row): a private function which evaluates the 1076 source->model function at the specified coords. 1077 1078 NOTE: The coords are in source->pixel coords, not image coords. 1079 1080 XXX: reverse order of row,col args? 1081 1082 XXX: This should probably be a public function. 1083 1084 XXX: We don't care about the model derivatives here. So, for the individual 1085 model evaluation functions, we should probably separate the derivative 1086 calculation. 1087 1088 XXX: Get rid of deriv. 1089 1090 XXX: Use static vectors for x. 1091 1092 XXX: Figure out if it's (row, col) or (col, row) for the model functions. 1093 *****************************************************************************/ 1094 psF32 evalModel(psSource *src, 1095 psU32 row, 1096 psU32 col) 1097 { 1098 // XXX: You won't need this once you separate the function evaluation. 1099 // from the derivative calculation. 1100 psVector *deriv = psVectorAlloc(src->models->Nparams, PS_TYPE_F32); 1101 // XXX: The following step will not be necessary if the models->params 1102 // member is a psVector. 1103 psVector *params = psVectorAlloc(src->models->Nparams, PS_TYPE_F32); 1104 for (psS32 i = 0 ; i < src->models->Nparams ; i++) { 1105 params->data.F32[i] = src->models->params[i]; 1106 } 1107 psVector *x = psVectorAlloc(2, PS_TYPE_F32); 1108 x->data.F32[0] = (psF32) (col + src->pixels->col0); 1109 x->data.F32[1] = (psF32) (row + src->pixels->row0); 1110 psF32 tmpF; 1111 1112 1113 switch (src->models->type) { 1114 case PS_MODEL_GAUSS: 1115 tmpF = pmMinLM_Gauss2D(deriv, params, x); 1116 break; 1117 case PS_MODEL_PGAUSS: 1118 tmpF = pmMinLM_PsuedoGauss2D(deriv, params, x); 1119 break; 1120 case PS_MODEL_TWIST_GAUSS: 1121 tmpF = pmMinLM_TwistGauss2D(deriv, params, x); 1122 break; 1123 case PS_MODEL_WAUSS: 1124 tmpF = pmMinLM_Wauss2D(deriv, params, x); 1125 break; 1126 case PS_MODEL_SERSIC: 1127 tmpF = pmMinLM_Sersic(deriv, params, x); 1128 break; 1129 case PS_MODEL_SERSIC_CORE: 1130 tmpF = pmMinLM_SersicCore(deriv, params, x); 1131 break; 1132 default: 1133 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1134 return(NAN); 1135 } 1136 1137 psFree(deriv); 1138 psFree(params); 1139 psFree(x); 1140 return(tmpF); 1141 } 1142 1143 /****************************************************************************** 1144 findValue(source, level, row, col, dir): a private function which determines 1145 the column coordinate of the model function which has the value "level". If 1146 dir equals 0, then you loop leftwards from the peak pixel, otherwise, 1147 rightwards. 1148 1149 XXX: reverse order of row,col args? 1150 1151 XXX: The result is returned in subImage coords, not image coords. 1152 *****************************************************************************/ 1153 psF32 findValue(psSource *source, 1154 psF32 level, 1155 psU32 row, 1156 psU32 col, 1157 psU32 dir) 1158 { 1159 // Ensure that the starting column is allowable. 1160 if (!((0 <= col) && (col < source->pixels->numCols))) { 1161 psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range"); 1162 return(NAN); 1163 } 1164 1165 psF32 oldValue = evalModel(source, row, col); 1166 if (oldValue == level) { 1167 return(((psF32) col)); 1168 } 1169 1170 // 1171 // We define variables incr and lastColumn so that we can use the same loop 1172 // whether we are stepping leftwards, or rightwards. 1173 // 1174 psS32 incr; 1175 psS32 lastColumn; 1176 if (dir == 0) { 1177 incr = -1; 1178 lastColumn = -1; 1179 } else { 1180 incr = 1; 1181 lastColumn = source->pixels->numCols; 1182 } 1183 col+=incr; 1184 1185 while (col != lastColumn) { 1186 psF32 newValue = evalModel(source, row, col); 1187 if (oldValue == level) { 1188 return(col); 1189 } 1190 1191 if ((newValue <= level) && (level <= oldValue)) { 1192 // This is simple linear interpolation. 1193 return( ((psF32) col) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) ); 1194 } 1195 1196 if ((oldValue <= level) && (level <= newValue)) { 1197 // This is simple linear interpolation. 1198 return( ((psF32) col) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) ); 1199 } 1200 1201 col+=incr; 1202 } 1203 1204 return(NAN); 1205 } 1206 /****************************************************************************** 1207 XXX: Probably should remove the "image" argument. 1208 XXX: What type should the output coordinate vectors consist of? 1209 XXX: Why a pmArray output? 1210 XXX: doex x,y correspond with col,row or row/col? 1211 *****************************************************************************/ 1212 psArray *pmSourceContour(psSource *source, 1213 const psImage *image, 1214 psF32 level, 1215 pmContourType mode) 1216 { 1217 PS_PTR_CHECK_NULL(source, false); 1218 PS_PTR_CHECK_NULL(source->moments, false); 1219 PS_PTR_CHECK_NULL(source->peak, false); 1220 PS_PTR_CHECK_NULL(source->pixels, false); 1221 PS_PTR_CHECK_NULL(source->models, false); 1222 1223 // 1224 // Allocate data for x/y pairs. 1225 // 1226 psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); 1227 psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); 1228 1229 // 1230 // Start at the row with peak pixel, then decrement. 1231 // 1232 psS32 col = source->peak->x; 1233 for (psS32 row = source->peak->y; row>= 0 ; row--) { 1234 // XXX: yVec contain no real information. Do we really need it? 1235 yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); 1236 yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); 1237 1238 // Starting at peak pixel, search leftwards for the column intercept. 1239 psF32 leftIntercept = findValue(source, level, row, col, 0); 1240 xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; 1241 1242 // Starting at peak pixel, search rightwards for the column intercept. 1243 psF32 rightIntercept = findValue(source, level, row, col, 1); 1244 xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; 1245 1246 // Set starting column for next row 1247 col = (psS32) ((leftIntercept + rightIntercept) / 2.0); 1248 } 1249 // 1250 // Start at the row (+1) with peak pixel, then increment. 1251 // 1252 col = source->peak->x; 1253 for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) { 1254 // XXX: yVec contain no real information. Do we really need it? 1255 yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); 1256 yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); 1257 1258 // Starting at peak pixel, search leftwards for the column intercept. 1259 psF32 leftIntercept = findValue(source, level, row, col, 0); 1260 xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; 1261 1262 // Starting at peak pixel, search rightwards for the column intercept. 1263 psF32 rightIntercept = findValue(source, level, row, col, 1); 1264 xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; 1265 1266 // Set starting column for next row 1267 col = (psS32) ((leftIntercept + rightIntercept) / 2.0); 1268 } 1269 1270 // 1271 // Allocate an array for result, store coord vectors there. 1272 // 1273 psArray *tmpArray = psArrayAlloc(2); 1274 tmpArray->data[0] = (psPtr *) yVec; 1275 tmpArray->data[1] = (psPtr *) xVec; 1276 return(tmpArray); 1277 } 1278 1279 psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x); 1280 psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); 1281 psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv, psVector *params, psArray *x); 1282 psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); 1283 psVector *p_pmMinLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x); 1284 psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x); 1285 1286 //XXX: What should these values be? 1287 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 100 1288 #define PM_SOURCE_FIT_MODEL_TOLERANCE 1.0 1289 /****************************************************************************** 1290 pmSourceFitModel(source, image): must create the appropiate arguments to the 1291 LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions. 1292 1293 XXX: should there be a mask value? 1294 XXX: Probably should remove the "image" argument. 1295 *****************************************************************************/ 1296 bool pmSourceFitModel(psSource *source, 1297 const psImage *image) 1298 { 1299 PS_PTR_CHECK_NULL(source, false); 1300 PS_PTR_CHECK_NULL(source->moments, false); 1301 PS_PTR_CHECK_NULL(source->peak, false); 1302 PS_PTR_CHECK_NULL(source->pixels, false); 1303 PS_PTR_CHECK_NULL(source->models, false); 1304 PS_IMAGE_CHECK_NULL(image, false); 1305 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); 1306 1307 psBool rc; 1308 psS32 count = 0; 1309 for (psS32 i = 0 ; i < source->pixels->numRows ; i++) { 1310 for (psS32 j = 0 ; j < source->pixels->numCols ; j++) { 1311 if (source->mask->data.U8[i][j] == 0) { 1312 count++; 1313 } 1314 } 1315 } 1316 psArray *x = psArrayAlloc(count); 1317 psVector *y = psVectorAlloc(count, PS_TYPE_F32); 1318 for (psS32 i = 0 ; i < source->pixels->numRows ; i++) { 1319 for (psS32 j = 0 ; j < source->pixels->numCols ; j++) { 1320 if (source->mask->data.U8[i][j] == 0) { 1321 psVector *coord = psVectorAlloc(2, PS_TYPE_F32); 1322 // XXX: Should we use the subimage offsets here, or does it not matter? 1323 coord->data.F32[0] = (psF32) (i); 1324 coord->data.F32[1] = (psF32) (j); 1325 x->data[count] = (psPtr *) coord; 1326 y->data.F32[count] = source->pixels->data.F32[i][j]; 1327 } 1328 } 1329 } 1330 1331 psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, 1332 PM_SOURCE_FIT_MODEL_TOLERANCE); 1333 1334 psVector *params = psVectorAlloc(source->models->Nparams, PS_TYPE_F32); 1335 1336 switch (source->models->type) { 1337 case PS_MODEL_GAUSS: 1338 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1339 NULL, (psMinimizeLMChi2Func) p_pmMinLM_Gauss2D_Vec); 1340 break; 1341 case PS_MODEL_PGAUSS: 1342 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1343 NULL, (psMinimizeLMChi2Func) p_pmMinLM_PsuedoGauss2D_Vec); 1344 break; 1345 case PS_MODEL_TWIST_GAUSS: 1346 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1347 NULL, (psMinimizeLMChi2Func) p_pmMinLM_Wauss2D_Vec); 1348 break; 1349 case PS_MODEL_WAUSS: 1350 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1351 NULL, (psMinimizeLMChi2Func) p_pmMinLM_TwistGauss2D_Vec); 1352 break; 1353 case PS_MODEL_SERSIC: 1354 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1355 NULL, (psMinimizeLMChi2Func) p_pmMinLM_Sersic_Vec); 1356 break; 1357 case PS_MODEL_SERSIC_CORE: 1358 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1359 NULL, (psMinimizeLMChi2Func) p_pmMinLM_SersicCore_Vec); 1360 break; 1361 default: 1362 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1363 rc = false; 1364 } 1365 1366 psFree(x); 1367 psFree(y); 1368 psFree(myMin); 1369 psFree(params); 1370 return(rc); 1371 } 1372 1373 bool p_pmSourceAddOrSubModel(psImage *image, 1374 psSource *src, 1375 bool center, 1376 psS32 flag) 1377 { 1378 PS_PTR_CHECK_NULL(src, false); 1379 PS_PTR_CHECK_NULL(src->moments, false); 1380 PS_PTR_CHECK_NULL(src->peak, false); 1381 PS_PTR_CHECK_NULL(src->pixels, false); 1382 PS_PTR_CHECK_NULL(src->models, false); 1383 PS_IMAGE_CHECK_NULL(image, false); 1384 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); 1385 1386 // XXX: make static, document. 1387 // We use a length of 100 here since it is larger enough for all models. 1388 #define MAX_PARAMS 1389 1390 psVector *deriv = psVectorAlloc(MAX_PARAMS, PS_TYPE_F32); 1391 psVector *params = psVectorAlloc(MAX_PARAMS, PS_TYPE_F32); 1392 psVector *x = psVectorAlloc(2, PS_TYPE_F32); 1393 if (src->models->Nparams > MAX_PARAMS) { 1394 psError(PS_ERR_UNKNOWN, true, "Internal error: increase MAX_PARAMS to %d", src->models->Nparams); 1395 } 1396 for (psS32 i = 0 ; i < src->models->Nparams ; i++) { 1397 params->data.F32[i] = src->models->params[i]; 1398 } 1399 1400 for (psS32 i = 0 ; i < src->pixels->numRows ; i++) { 1401 for (psS32 j = 0 ; j < src->pixels->numCols ; j++) { 1402 psF32 pixelValue; 1403 // XXX: Should you use offsets here? 1404 // XXX: Make sure you have col/row order correct. 1405 // XXX: Should you be adding the pixels for the entire subImage, 1406 // or a radius of pixels around it? 1407 1408 x->data.F32[0] = (float) j; 1409 x->data.F32[1] = (float) i; 1410 switch (src->models->type) { 1411 case PS_MODEL_GAUSS: 1412 pixelValue = pmMinLM_Gauss2D(deriv, params, x); 1413 break; 1414 case PS_MODEL_PGAUSS: 1415 pixelValue = pmMinLM_PsuedoGauss2D(deriv, params, x); 1416 break; 1417 case PS_MODEL_TWIST_GAUSS: 1418 pixelValue = pmMinLM_TwistGauss2D(deriv, params, x); 1419 break; 1420 case PS_MODEL_WAUSS: 1421 pixelValue = pmMinLM_Wauss2D(deriv, params, x); 1422 break; 1423 case PS_MODEL_SERSIC: 1424 pixelValue = pmMinLM_Sersic(deriv, params, x); 1425 break; 1426 case PS_MODEL_SERSIC_CORE: 1427 pixelValue = pmMinLM_SersicCore(deriv, params, x); 1428 break; 1429 default: 1430 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1431 psFree(x); 1432 psFree(deriv); 1433 psFree(params); 1434 return(false); 1435 } 1436 if (flag == 1) { 1437 pixelValue = -pixelValue; 1438 } 1439 1440 // XXX: Must figure out how to calculate the image coordinates and 1441 // how to use the boolean "center" flag. 1442 psS32 imageRow = i + src->pixels->row0; 1443 psS32 imageCol = j + src->pixles->col0; 1444 1445 image->data.F32[imageRow][imageCol]+= pixelValue; 1446 } 1447 } 1448 psFree(x); 1449 psFree(deriv); 1450 psFree(params); 1451 return(true); 1452 } 1453 1454 1455 1456 /****************************************************************************** 1457 *****************************************************************************/ 1458 bool pmSourceAddModel(psImage *image, 1459 psSource *src, 1460 bool center) 1461 { 1462 return(p_pmSourceAddOrSubModel(image, src, center, 0)); 1463 } 1464 1465 /****************************************************************************** 1466 *****************************************************************************/ 1467 bool pmSourceSubModel(psImage *image, 1468 psSource *src, 1469 bool center) 1470 { 1471 return(p_pmSourceAddOrSubModel(image, src, center, 1)); 1472 } 1473 1474 1475 // XXX: Put this is psConstants.h 1476 #define PS_VECTOR_CHECK_SIZE(VEC1, N, RVAL) \ 1477 if (VEC1->n != N) { \ 1478 psError(PS_ERR_BAD_PARAMETER_SIZE, true, \ 1479 "psVector %s has size %d, should be %d.", \ 1480 #VEC1, VEC1->n, N); \ 1481 return(RVAL); \ 1482 } 1483 1484 1485 /****************************************************************************** 1486 pmMinLM_Gauss2D(*deriv, *params, *x): the argument "x" contains a single "x, 1487 y" coordinate pair. This function computes the gaussian, specified by the 1488 parameters in "params" at that x,y point and returns the value. The 1489 derivatives are also caculated and returned in the "deriv" argument. 1490 1491 params->data.F32[0] = So; 1492 params->data.F32[1] = Zo; 1493 params->data.F32[2] = Xo; 1494 params->data.F32[3] = Yo; 1495 params->data.F32[4] = sqrt(2.0) / SigmaX; 1496 params->data.F32[5] = sqrt(2.0) / SigmaY; 1497 params->data.F32[6] = Sxy; 1498 1499 XXX: Consider getting rid of the parameter checks since this might consume 1500 a significant fraction of this function CPU time. 1501 *****************************************************************************/ 1502 psF32 pmMinLM_Gauss2D(psVector *deriv, 803 } 804 805 if (tmpSrc->peak->class == 0) { 806 tmpSrc->peak->class|= PS_SOURCE_OTHER; 807 } 808 } 809 810 return(rc); 811 } 812 813 814 815 /****************************************************************************** 816 pmSourceSetPixelCircle(source, image, radius) 817 818 XXX: Why boolean output? 819 820 XXX: Why are we checking source->moments for NULL? Should the circle be 821 centered on the centroid or the peak? 822 823 XXX: The circle will have a diameter of (1+radius). This is different from 824 the pmSourceSetLocal() function. 825 *****************************************************************************/ 826 bool pmSourceSetPixelCircle(psSource *source, 827 const psImage *image, 828 psF32 radius) 829 { 830 PS_IMAGE_CHECK_NULL(image, false); 831 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); 832 PS_PTR_CHECK_NULL(source, false); 833 // PS_PTR_CHECK_NULL(source->moments, false); 834 PS_PTR_CHECK_NULL(source->peak, false); 835 PS_FLOAT_COMPARE(0.0, radius, false); 836 837 // 838 // We define variables for code readability. 839 // 840 psS32 radiusS32 = (psS32) radius; 841 psS32 SubImageCenterRow = source->peak->y; 842 psS32 SubImageCenterCol = source->peak->x; 843 psS32 SubImageStartRow = SubImageCenterRow - radiusS32; 844 psS32 SubImageEndRow = SubImageCenterRow + radiusS32; 845 psS32 SubImageStartCol = SubImageCenterCol - radiusS32; 846 psS32 SubImageEndCol = SubImageCenterCol + radiusS32; 847 848 if (SubImageStartRow < 0) { 849 psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n", 850 SubImageStartRow); 851 return(false); 852 } 853 if (SubImageEndRow+1 >= image->numRows) { 854 psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n", 855 SubImageEndRow); 856 return(false); 857 } 858 if (SubImageStartCol < 0) { 859 psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n", 860 SubImageStartCol); 861 return(false); 862 } 863 if (SubImageEndCol+1 >= image->numCols) { 864 psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n", 865 SubImageEndCol); 866 return(false); 867 } 868 869 // XXX: Must recycle image. 870 if (source->pixels != NULL) { 871 psLogMsg(__func__, PS_LOG_WARN, 872 "WARNING: pmSourceSetPixelCircle(): image->pixels not NULL. Freeing and reallocating.\n"); 873 psFree(source->pixels); 874 } 875 source->pixels = psImageSubset((psImage *) image, 876 SubImageStartCol, 877 SubImageStartRow, 878 SubImageEndCol+1, 879 SubImageEndRow+1); 880 881 // XXX: Must recycle image. 882 if (source->mask != NULL) { 883 psFree(source->mask); 884 } 885 source->mask = psImageAlloc(1 + 2 * radiusS32, 886 1 + 2 * radiusS32, 887 PS_TYPE_F32); 888 889 // 890 // Loop through the subimage mask, initialize mask to 0 or 1. 891 // 892 for (psS32 row = 0 ; row < source->mask->numRows; row++) { 893 for (psS32 col = 0 ; col < source->mask->numCols; col++) { 894 895 if (CheckRadius2((psF32) radiusS32, 896 (psF32) radiusS32, 897 radius, 898 (psF32) col, 899 (psF32) row)) { 900 source->mask->data.U8[row][col] = 1; 901 } else { 902 source->mask->data.U8[row][col] = 0; 903 } 904 } 905 } 906 return(true); 907 } 908 909 910 /****************************************************************************** 911 pmSourceModelGuess(source, image, model): This function allocates a new 912 psModel structure and store it in the psSource data structure specified in the 913 argument list. The model type is specified in the argument list. The params 914 array in that psModel structure are allocated, and then set to the appropriate 915 values. This function returns true if everything was successful. 916 917 XXX: Many of the initial parameters are set to 0.0 since I don't know what 918 the appropiate initial guesses are. 919 *****************************************************************************/ 920 bool pmSourceModelGuess(psSource *source, 921 const psImage *image, 922 psModelType model) 923 { 924 PS_PTR_CHECK_NULL(source, false); 925 PS_PTR_CHECK_NULL(source->moments, false); 926 PS_PTR_CHECK_NULL(source->peak, false); 927 PS_IMAGE_CHECK_NULL(image, false); 928 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); 929 if (source->models != NULL) { 930 psLogMsg(__func__, PS_LOG_WARN, "WARNING: source->models was non-NULL; calling psFree(source->models).\n"); 931 psFree(source->models); 932 } 933 source->models = pmModelAlloc(PS_MODEL_UNDEFINED); 934 935 switch (model) { 936 case PS_MODEL_GAUSS: 937 source->models->type = PS_MODEL_GAUSS; 938 source->models->Nparams = 7; 939 source->models->params = (psF32 *) psAlloc(7 * sizeof(psF32)); 940 source->models->dparams = (psF32 *) psAlloc(7 * sizeof(psF32)); 941 for (psS32 i = 0 ; i < 7 ; i++) { 942 source->models->params[i] = 0.0; 943 source->models->dparams[i] = 0.0; 944 } 945 source->models->params[0] = source->moments->Sky; 946 source->models->params[1] = source->peak->counts - source->moments->Sky; 947 source->models->params[2] = source->moments->x; 948 source->models->params[3] = source->moments->y; 949 source->models->params[4] = sqrt(2.0) / source->moments->Sx; 950 source->models->params[5] = sqrt(2.0) / source->moments->Sy; 951 source->models->params[6] = source->moments->Sxy; 952 source->models->chisq = 0.0; 953 return(true); 954 case PS_MODEL_PGAUSS: 955 source->models->type = PS_MODEL_PGAUSS; 956 source->models->Nparams = 7; 957 source->models->params = (psF32 *) psAlloc(7 * sizeof(psF32)); 958 source->models->dparams = (psF32 *) psAlloc(7 * sizeof(psF32)); 959 for (psS32 i = 0 ; i < 7 ; i++) { 960 source->models->params[i] = 0.0; 961 source->models->dparams[i] = 0.0; 962 } 963 source->models->params[0] = source->moments->Sky; 964 source->models->params[1] = source->peak->counts - source->moments->Sky; 965 source->models->params[2] = source->moments->x; 966 source->models->params[3] = source->moments->y; 967 source->models->params[4] = sqrt(2.0) / source->moments->Sx; 968 source->models->params[5] = sqrt(2.0) / source->moments->Sy; 969 source->models->params[6] = source->moments->Sxy; 970 source->models->chisq = 0.0; 971 return(true); 972 case PS_MODEL_TWIST_GAUSS: 973 source->models->type = PS_MODEL_TWIST_GAUSS; 974 source->models->Nparams = 11; 975 source->models->params = (psF32 *) psAlloc(11 * sizeof(psF32)); 976 source->models->dparams = (psF32 *) psAlloc(11 * sizeof(psF32)); 977 for (psS32 i = 0 ; i < 11 ; i++) { 978 source->models->params[i] = 0.0; 979 source->models->dparams[i] = 0.0; 980 } 981 982 source->models->params[0] = source->moments->Sky; 983 source->models->params[1] = source->peak->counts - source->moments->Sky; 984 source->models->params[2] = source->moments->x; 985 source->models->params[3] = source->moments->y; 986 // XXX: What are these? 987 // source->models->params[4] = SxInner; 988 // source->models->params[5] = SyInner; 989 // source->models->params[6] = SxyInner; 990 // source->models->params[7] = SxOuter; 991 // source->models->params[8] = SyOuter; 992 // source->models->params[9] = SxyOuter; 993 // source->models->params[10] = N; 994 995 source->models->chisq = 0.0; 996 return(true); 997 case PS_MODEL_WAUSS: 998 999 source->models->params[0] = source->moments->Sky; 1000 source->models->params[1] = source->peak->counts - source->moments->Sky; 1001 source->models->params[2] = source->moments->x; 1002 source->models->params[3] = source->moments->y; 1003 source->models->params[4] = sqrt(2.0) / source->moments->Sx; 1004 source->models->params[5] = sqrt(2.0) / source->moments->Sy; 1005 source->models->params[6] = source->moments->Sxy; 1006 // XXX: What are these? 1007 // source->models->params[7] = B2; 1008 // source->models->params[8] = B3; 1009 1010 source->models->type = PS_MODEL_WAUSS; 1011 source->models->Nparams = 9; 1012 source->models->params = (psF32 *) psAlloc(9 * sizeof(psF32)); 1013 source->models->dparams = (psF32 *) psAlloc(9 * sizeof(psF32)); 1014 for (psS32 i = 0 ; i < 9 ; i++) { 1015 source->models->params[i] = 0.0; 1016 source->models->dparams[i] = 0.0; 1017 } 1018 source->models->chisq = 0.0; 1019 return(true); 1020 case PS_MODEL_SERSIC: 1021 source->models->type = PS_MODEL_SERSIC; 1022 source->models->Nparams = 8; 1023 source->models->params = (psF32 *) psAlloc(8 * sizeof(psF32)); 1024 source->models->dparams = (psF32 *) psAlloc(8 * sizeof(psF32)); 1025 for (psS32 i = 0 ; i < 8 ; i++) { 1026 source->models->params[i] = 0.0; 1027 source->models->dparams[i] = 0.0; 1028 } 1029 1030 source->models->params[0] = source->moments->Sky; 1031 source->models->params[1] = source->peak->counts - source->moments->Sky; 1032 source->models->params[2] = source->moments->x; 1033 source->models->params[3] = source->moments->y; 1034 source->models->params[4] = sqrt(2.0) / source->moments->Sx; 1035 source->models->params[5] = sqrt(2.0) / source->moments->Sy; 1036 source->models->params[6] = source->moments->Sxy; 1037 // XXX: What are these? 1038 //source->models->params[7] = Nexp; 1039 1040 source->models->chisq = 0.0; 1041 return(true); 1042 case PS_MODEL_SERSIC_CORE: 1043 source->models->type = PS_MODEL_SERSIC_CORE; 1044 source->models->Nparams = 12; 1045 source->models->params = (psF32 *) psAlloc(12 * sizeof(psF32)); 1046 source->models->dparams = (psF32 *) psAlloc(12 * sizeof(psF32)); 1047 for (psS32 i = 0 ; i < 12 ; i++) { 1048 source->models->params[i] = 0.0; 1049 source->models->dparams[i] = 0.0; 1050 } 1051 1052 source->models->params[0] = source->moments->Sky; 1053 source->models->params[1] = source->peak->counts - source->moments->Sky; 1054 source->models->params[2] = source->moments->x; 1055 source->models->params[3] = source->moments->y; 1056 // XXX: What are these? 1057 //source->models->params[4] SxInner; 1058 //source->models->params[5] SyInner; 1059 //source->models->params[6] SxyInner; 1060 //source->models->params[7] Zd; 1061 //source->models->params[8] SxOuter; 1062 //source->models->params[9] SyOuter; 1063 //source->models->params[10] = SxyOuter; 1064 //source->models->params[11] = Nexp; 1065 1066 source->models->chisq = 0.0; 1067 return(true); 1068 default: 1069 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1070 return(false); 1071 } 1072 } 1073 1074 /****************************************************************************** 1075 evalModel(source, level, row): a private function which evaluates the 1076 source->model function at the specified coords. 1077 1078 NOTE: The coords are in source->pixel coords, not image coords. 1079 1080 XXX: reverse order of row,col args? 1081 1082 XXX: This should probably be a public function. 1083 1084 XXX: We don't care about the model derivatives here. So, for the individual 1085 model evaluation functions, we should probably separate the derivative 1086 calculation. 1087 1088 XXX: Get rid of deriv. 1089 1090 XXX: Use static vectors for x. 1091 1092 XXX: Figure out if it's (row, col) or (col, row) for the model functions. 1093 *****************************************************************************/ 1094 psF32 evalModel(psSource *src, 1095 psU32 row, 1096 psU32 col) 1097 { 1098 // XXX: You won't need this once you separate the function evaluation. 1099 // from the derivative calculation. 1100 psVector *deriv = psVectorAlloc(src->models->Nparams, PS_TYPE_F32); 1101 // XXX: The following step will not be necessary if the models->params 1102 // member is a psVector. 1103 psVector *params = psVectorAlloc(src->models->Nparams, PS_TYPE_F32); 1104 for (psS32 i = 0 ; i < src->models->Nparams ; i++) { 1105 params->data.F32[i] = src->models->params[i]; 1106 } 1107 psVector *x = psVectorAlloc(2, PS_TYPE_F32); 1108 x->data.F32[0] = (psF32) (col + src->pixels->col0); 1109 x->data.F32[1] = (psF32) (row + src->pixels->row0); 1110 psF32 tmpF; 1111 1112 1113 switch (src->models->type) { 1114 case PS_MODEL_GAUSS: 1115 tmpF = pmMinLM_Gauss2D(deriv, params, x); 1116 break; 1117 case PS_MODEL_PGAUSS: 1118 tmpF = pmMinLM_PsuedoGauss2D(deriv, params, x); 1119 break; 1120 case PS_MODEL_TWIST_GAUSS: 1121 tmpF = pmMinLM_TwistGauss2D(deriv, params, x); 1122 break; 1123 case PS_MODEL_WAUSS: 1124 tmpF = pmMinLM_Wauss2D(deriv, params, x); 1125 break; 1126 case PS_MODEL_SERSIC: 1127 tmpF = pmMinLM_Sersic(deriv, params, x); 1128 break; 1129 case PS_MODEL_SERSIC_CORE: 1130 tmpF = pmMinLM_SersicCore(deriv, params, x); 1131 break; 1132 default: 1133 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1134 return(NAN); 1135 } 1136 1137 psFree(deriv); 1138 psFree(params); 1139 psFree(x); 1140 return(tmpF); 1141 } 1142 1143 /****************************************************************************** 1144 findValue(source, level, row, col, dir): a private function which determines 1145 the column coordinate of the model function which has the value "level". If 1146 dir equals 0, then you loop leftwards from the peak pixel, otherwise, 1147 rightwards. 1148 1149 XXX: reverse order of row,col args? 1150 1151 XXX: The result is returned in subImage coords, not image coords. 1152 *****************************************************************************/ 1153 psF32 findValue(psSource *source, 1154 psF32 level, 1155 psU32 row, 1156 psU32 col, 1157 psU32 dir) 1158 { 1159 // Ensure that the starting column is allowable. 1160 if (!((0 <= col) && (col < source->pixels->numCols))) { 1161 psError(PS_ERR_UNKNOWN, true, "Starting column outside subImage range"); 1162 return(NAN); 1163 } 1164 1165 psF32 oldValue = evalModel(source, row, col); 1166 if (oldValue == level) { 1167 return(((psF32) col)); 1168 } 1169 1170 // 1171 // We define variables incr and lastColumn so that we can use the same loop 1172 // whether we are stepping leftwards, or rightwards. 1173 // 1174 psS32 incr; 1175 psS32 lastColumn; 1176 if (dir == 0) { 1177 incr = -1; 1178 lastColumn = -1; 1179 } else { 1180 incr = 1; 1181 lastColumn = source->pixels->numCols; 1182 } 1183 col+=incr; 1184 1185 while (col != lastColumn) { 1186 psF32 newValue = evalModel(source, row, col); 1187 if (oldValue == level) { 1188 return(col); 1189 } 1190 1191 if ((newValue <= level) && (level <= oldValue)) { 1192 // This is simple linear interpolation. 1193 return( ((psF32) col) + ((psF32) incr) * ((level - newValue) / (oldValue - newValue)) ); 1194 } 1195 1196 if ((oldValue <= level) && (level <= newValue)) { 1197 // This is simple linear interpolation. 1198 return( ((psF32) col) + ((psF32) incr) * ((level - oldValue) / (newValue - oldValue)) ); 1199 } 1200 1201 col+=incr; 1202 } 1203 1204 return(NAN); 1205 } 1206 /****************************************************************************** 1207 XXX: Probably should remove the "image" argument. 1208 XXX: What type should the output coordinate vectors consist of? 1209 XXX: Why a pmArray output? 1210 XXX: doex x,y correspond with col,row or row/col? 1211 *****************************************************************************/ 1212 psArray *pmSourceContour(psSource *source, 1213 const psImage *image, 1214 psF32 level, 1215 pmContourType mode) 1216 { 1217 PS_PTR_CHECK_NULL(source, false); 1218 PS_PTR_CHECK_NULL(source->moments, false); 1219 PS_PTR_CHECK_NULL(source->peak, false); 1220 PS_PTR_CHECK_NULL(source->pixels, false); 1221 PS_PTR_CHECK_NULL(source->models, false); 1222 1223 // 1224 // Allocate data for x/y pairs. 1225 // 1226 psVector *xVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); 1227 psVector *yVec = psVectorAlloc(2 * source->pixels->numRows, PS_TYPE_F32); 1228 1229 // 1230 // Start at the row with peak pixel, then decrement. 1231 // 1232 psS32 col = source->peak->x; 1233 for (psS32 row = source->peak->y; row>= 0 ; row--) { 1234 // XXX: yVec contain no real information. Do we really need it? 1235 yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); 1236 yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); 1237 1238 // Starting at peak pixel, search leftwards for the column intercept. 1239 psF32 leftIntercept = findValue(source, level, row, col, 0); 1240 xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; 1241 1242 // Starting at peak pixel, search rightwards for the column intercept. 1243 psF32 rightIntercept = findValue(source, level, row, col, 1); 1244 xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; 1245 1246 // Set starting column for next row 1247 col = (psS32) ((leftIntercept + rightIntercept) / 2.0); 1248 } 1249 // 1250 // Start at the row (+1) with peak pixel, then increment. 1251 // 1252 col = source->peak->x; 1253 for (psS32 row = 1 + source->peak->y; row < source->pixels->numRows ; row++) { 1254 // XXX: yVec contain no real information. Do we really need it? 1255 yVec->data.F32[row] = (psF32) (source->pixels->row0 + row); 1256 yVec->data.F32[row+yVec->n] = (psF32) (source->pixels->row0 + row); 1257 1258 // Starting at peak pixel, search leftwards for the column intercept. 1259 psF32 leftIntercept = findValue(source, level, row, col, 0); 1260 xVec->data.F32[row] = ((psF32) source->pixels->col0) + leftIntercept; 1261 1262 // Starting at peak pixel, search rightwards for the column intercept. 1263 psF32 rightIntercept = findValue(source, level, row, col, 1); 1264 xVec->data.F32[row+xVec->n] = ((psF32) source->pixels->col0) + rightIntercept; 1265 1266 // Set starting column for next row 1267 col = (psS32) ((leftIntercept + rightIntercept) / 2.0); 1268 } 1269 1270 // 1271 // Allocate an array for result, store coord vectors there. 1272 // 1273 psArray *tmpArray = psArrayAlloc(2); 1274 tmpArray->data[0] = (psPtr *) yVec; 1275 tmpArray->data[1] = (psPtr *) xVec; 1276 return(tmpArray); 1277 } 1278 1279 psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x); 1280 psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); 1281 psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv, psVector *params, psArray *x); 1282 psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); 1283 psVector *p_pmMinLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x); 1284 psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x); 1285 1286 //XXX: What should these values be? 1287 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 100 1288 #define PM_SOURCE_FIT_MODEL_TOLERANCE 1.0 1289 /****************************************************************************** 1290 pmSourceFitModel(source, image): must create the appropiate arguments to the 1291 LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions. 1292 1293 XXX: should there be a mask value? 1294 XXX: Probably should remove the "image" argument. 1295 *****************************************************************************/ 1296 bool pmSourceFitModel(psSource *source, 1297 const psImage *image) 1298 { 1299 PS_PTR_CHECK_NULL(source, false); 1300 PS_PTR_CHECK_NULL(source->moments, false); 1301 PS_PTR_CHECK_NULL(source->peak, false); 1302 PS_PTR_CHECK_NULL(source->pixels, false); 1303 PS_PTR_CHECK_NULL(source->models, false); 1304 PS_IMAGE_CHECK_NULL(image, false); 1305 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); 1306 1307 psBool rc; 1308 psS32 count = 0; 1309 for (psS32 i = 0 ; i < source->pixels->numRows ; i++) { 1310 for (psS32 j = 0 ; j < source->pixels->numCols ; j++) { 1311 if (source->mask->data.U8[i][j] == 0) { 1312 count++; 1313 } 1314 } 1315 } 1316 psArray *x = psArrayAlloc(count); 1317 psVector *y = psVectorAlloc(count, PS_TYPE_F32); 1318 for (psS32 i = 0 ; i < source->pixels->numRows ; i++) { 1319 for (psS32 j = 0 ; j < source->pixels->numCols ; j++) { 1320 if (source->mask->data.U8[i][j] == 0) { 1321 psVector *coord = psVectorAlloc(2, PS_TYPE_F32); 1322 // XXX: Should we use the subimage offsets here, or does it not matter? 1323 coord->data.F32[0] = (psF32) (i); 1324 coord->data.F32[1] = (psF32) (j); 1325 x->data[count] = (psPtr *) coord; 1326 y->data.F32[count] = source->pixels->data.F32[i][j]; 1327 } 1328 } 1329 } 1330 1331 psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, 1332 PM_SOURCE_FIT_MODEL_TOLERANCE); 1333 1334 psVector *params = psVectorAlloc(source->models->Nparams, PS_TYPE_F32); 1335 1336 switch (source->models->type) { 1337 case PS_MODEL_GAUSS: 1338 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1339 NULL, (psMinimizeLMChi2Func) p_pmMinLM_Gauss2D_Vec); 1340 break; 1341 case PS_MODEL_PGAUSS: 1342 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1343 NULL, (psMinimizeLMChi2Func) p_pmMinLM_PsuedoGauss2D_Vec); 1344 break; 1345 case PS_MODEL_TWIST_GAUSS: 1346 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1347 NULL, (psMinimizeLMChi2Func) p_pmMinLM_Wauss2D_Vec); 1348 break; 1349 case PS_MODEL_WAUSS: 1350 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1351 NULL, (psMinimizeLMChi2Func) p_pmMinLM_TwistGauss2D_Vec); 1352 break; 1353 case PS_MODEL_SERSIC: 1354 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1355 NULL, (psMinimizeLMChi2Func) p_pmMinLM_Sersic_Vec); 1356 break; 1357 case PS_MODEL_SERSIC_CORE: 1358 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1359 NULL, (psMinimizeLMChi2Func) p_pmMinLM_SersicCore_Vec); 1360 break; 1361 default: 1362 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1363 rc = false; 1364 } 1365 1366 psFree(x); 1367 psFree(y); 1368 psFree(myMin); 1369 psFree(params); 1370 return(rc); 1371 } 1372 1373 bool p_pmSourceAddOrSubModel(psImage *image, 1374 psSource *src, 1375 bool center, 1376 psS32 flag) 1377 { 1378 PS_PTR_CHECK_NULL(src, false); 1379 PS_PTR_CHECK_NULL(src->moments, false); 1380 PS_PTR_CHECK_NULL(src->peak, false); 1381 PS_PTR_CHECK_NULL(src->pixels, false); 1382 PS_PTR_CHECK_NULL(src->models, false); 1383 PS_IMAGE_CHECK_NULL(image, false); 1384 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); 1385 1386 // XXX: make static, document. 1387 // We use a length of 100 here since it is larger enough for all models. 1388 #define MAX_PARAMS 1389 1390 psVector *deriv = psVectorAlloc(MAX_PARAMS, PS_TYPE_F32); 1391 psVector *params = psVectorAlloc(MAX_PARAMS, PS_TYPE_F32); 1392 psVector *x = psVectorAlloc(2, PS_TYPE_F32); 1393 if (src->models->Nparams > MAX_PARAMS) { 1394 psError(PS_ERR_UNKNOWN, true, "Internal error: increase MAX_PARAMS to %d", src->models->Nparams); 1395 } 1396 for (psS32 i = 0 ; i < src->models->Nparams ; i++) { 1397 params->data.F32[i] = src->models->params[i]; 1398 } 1399 1400 for (psS32 i = 0 ; i < src->pixels->numRows ; i++) { 1401 for (psS32 j = 0 ; j < src->pixels->numCols ; j++) { 1402 psF32 pixelValue; 1403 // XXX: Should you use offsets here? 1404 // XXX: Make sure you have col/row order correct. 1405 // XXX: Should you be adding the pixels for the entire subImage, 1406 // or a radius of pixels around it? 1407 1408 x->data.F32[0] = (float) j; 1409 x->data.F32[1] = (float) i; 1410 switch (src->models->type) { 1411 case PS_MODEL_GAUSS: 1412 pixelValue = pmMinLM_Gauss2D(deriv, params, x); 1413 break; 1414 case PS_MODEL_PGAUSS: 1415 pixelValue = pmMinLM_PsuedoGauss2D(deriv, params, x); 1416 break; 1417 case PS_MODEL_TWIST_GAUSS: 1418 pixelValue = pmMinLM_TwistGauss2D(deriv, params, x); 1419 break; 1420 case PS_MODEL_WAUSS: 1421 pixelValue = pmMinLM_Wauss2D(deriv, params, x); 1422 break; 1423 case PS_MODEL_SERSIC: 1424 pixelValue = pmMinLM_Sersic(deriv, params, x); 1425 break; 1426 case PS_MODEL_SERSIC_CORE: 1427 pixelValue = pmMinLM_SersicCore(deriv, params, x); 1428 break; 1429 default: 1430 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1431 psFree(x); 1432 psFree(deriv); 1433 psFree(params); 1434 return(false); 1435 } 1436 if (flag == 1) { 1437 pixelValue = -pixelValue; 1438 } 1439 1440 // XXX: Must figure out how to calculate the image coordinates and 1441 // how to use the boolean "center" flag. 1442 psS32 imageRow = i + src->pixels->row0; 1443 psS32 imageCol = j + src->pixles->col0; 1444 1445 image->data.F32[imageRow][imageCol]+= pixelValue; 1446 } 1447 } 1448 psFree(x); 1449 psFree(deriv); 1450 psFree(params); 1451 return(true); 1452 } 1453 1454 1455 1456 /****************************************************************************** 1457 *****************************************************************************/ 1458 bool pmSourceAddModel(psImage *image, 1459 psSource *src, 1460 bool center) 1461 { 1462 return(p_pmSourceAddOrSubModel(image, src, center, 0)); 1463 } 1464 1465 /****************************************************************************** 1466 *****************************************************************************/ 1467 bool pmSourceSubModel(psImage *image, 1468 psSource *src, 1469 bool center) 1470 { 1471 return(p_pmSourceAddOrSubModel(image, src, center, 1)); 1472 } 1473 1474 1475 // XXX: Put this is psConstants.h 1476 #define PS_VECTOR_CHECK_SIZE(VEC1, N, RVAL) \ 1477 if (VEC1->n != N) { \ 1478 psError(PS_ERR_BAD_PARAMETER_SIZE, true, \ 1479 "psVector %s has size %d, should be %d.", \ 1480 #VEC1, VEC1->n, N); \ 1481 return(RVAL); \ 1482 } 1483 1484 1485 /****************************************************************************** 1486 pmMinLM_Gauss2D(*deriv, *params, *x): the argument "x" contains a single "x, 1487 y" coordinate pair. This function computes the gaussian, specified by the 1488 parameters in "params" at that x,y point and returns the value. The 1489 derivatives are also caculated and returned in the "deriv" argument. 1490 1491 params->data.F32[0] = So; 1492 params->data.F32[1] = Zo; 1493 params->data.F32[2] = Xo; 1494 params->data.F32[3] = Yo; 1495 params->data.F32[4] = sqrt(2.0) / SigmaX; 1496 params->data.F32[5] = sqrt(2.0) / SigmaY; 1497 params->data.F32[6] = Sxy; 1498 1499 XXX: Consider getting rid of the parameter checks since this might consume 1500 a significant fraction of this function CPU time. 1501 *****************************************************************************/ 1502 psF32 pmMinLM_Gauss2D(psVector *deriv, 1503 psVector *params, 1504 psVector *x) 1505 { 1506 PS_VECTOR_CHECK_NULL(deriv, NAN); 1507 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1508 PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); 1509 PS_VECTOR_CHECK_NULL(params, NAN); 1510 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1511 PS_VECTOR_CHECK_SIZE(params, 7, NAN); 1512 PS_VECTOR_CHECK_NULL(x, NAN); 1513 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1514 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1515 1516 psF32 X = x->data.F32[0] - params->data.F32[2]; 1517 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1518 psF32 px = params->data.F32[4]*X; 1519 psF32 py = params->data.F32[5]*Y; 1520 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1521 psF32 r = exp(-z); 1522 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1523 1524 psF32 q = params->data.F32[1]*r; 1525 deriv->data.F32[0] = +1.0; 1526 deriv->data.F32[1] = +r; 1527 deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y); 1528 deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X); 1529 deriv->data.F32[4] = -2.0*q*px*X; 1530 deriv->data.F32[5] = -2.0*q*py*Y; 1531 deriv->data.F32[6] = -q*X*Y; 1532 1533 return(f); 1534 } 1535 1536 /****************************************************************************** 1537 p_pmMinLM_Gauss2D_Vec(*deriv, *params, *x): this function wraps the above 1538 function in a form that is usable in the LM minimization routines. 1539 *****************************************************************************/ 1540 psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv, 1541 psVector *params, 1542 psArray *x) 1543 { 1544 PS_IMAGE_CHECK_NULL(deriv, NULL); 1545 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1546 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1547 PS_VECTOR_CHECK_NULL(params, NULL); 1548 PS_VECTOR_CHECK_EMPTY(params, NULL); 1549 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1550 PS_PTR_CHECK_NULL(x, NULL); 1551 if (deriv->numRows != x->n) { 1552 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1553 } 1554 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1555 // XXX: use static memory here. 1556 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1557 1558 for (psS32 i = 0 ; i < x->n ; i++) { 1559 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1560 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1561 } 1562 1563 tmpVec->data.F32[i] = pmMinLM_Gauss2D(tmpRow, 1564 params, 1565 (psVector *) x->data[i]); 1566 1567 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1568 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1569 } 1570 } 1571 1572 psFree(tmpRow); 1573 return(tmpVec); 1574 } 1575 1576 1577 /****************************************************************************** 1578 params->data.F32[0] = So; 1579 params->data.F32[1] = Zo; 1580 params->data.F32[2] = Xo; 1581 params->data.F32[3] = Yo; 1582 params->data.F32[4] = sqrt(2) / SigmaX; 1583 params->data.F32[5] = sqrt(2) / SigmaY; 1584 params->data.F32[6] = Sxy; 1585 *****************************************************************************/ 1586 psF32 pmMinLM_PsuedoGauss2D(psVector *deriv, 1503 1587 psVector *params, 1504 1588 psVector *x) 1505 { 1506 PS_VECTOR_CHECK_NULL(deriv, NAN); 1507 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1508 PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); 1509 PS_VECTOR_CHECK_NULL(params, NAN); 1510 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1511 PS_VECTOR_CHECK_SIZE(params, 7, NAN); 1512 PS_VECTOR_CHECK_NULL(x, NAN); 1513 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1514 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1515 1516 psF32 X = x->data.F32[0] - params->data.F32[2]; 1517 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1518 psF32 px = params->data.F32[4]*X; 1519 psF32 py = params->data.F32[5]*Y; 1520 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1521 psF32 r = exp(-z); 1522 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1523 1524 psF32 q = params->data.F32[1]*r; 1525 deriv->data.F32[0] = +1.0; 1526 deriv->data.F32[1] = +r; 1527 deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y); 1528 deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X); 1529 deriv->data.F32[4] = -2.0*q*px*X; 1530 deriv->data.F32[5] = -2.0*q*py*Y; 1531 deriv->data.F32[6] = -q*X*Y; 1532 1533 return(f); 1534 } 1535 1536 /****************************************************************************** 1537 p_pmMinLM_Gauss2D_Vec(*deriv, *params, *x): this function wraps the above 1538 function in a form that is usable in the LM minimization routines. 1539 *****************************************************************************/ 1540 psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv, 1589 { 1590 PS_VECTOR_CHECK_NULL(deriv, NAN); 1591 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1592 PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); 1593 PS_VECTOR_CHECK_NULL(params, NAN); 1594 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1595 PS_VECTOR_CHECK_SIZE(params, 7, NAN); 1596 PS_VECTOR_CHECK_NULL(x, NAN); 1597 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1598 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1599 1600 psF32 X = x->data.F32[0] - params->data.F32[2]; 1601 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1602 psF32 px = params->data.F32[4]*X; 1603 psF32 py = params->data.F32[5]*Y; 1604 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1605 psF32 t = 1 + z + 0.5*z*z; 1606 psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */ 1607 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1608 1609 // 1610 // note difference from a pure gaussian: q = params->data.F32[1]*r 1611 // 1612 1613 psF32 q = params->data.F32[1]*r*r*t; 1614 deriv->data.F32[0] = +1.0; 1615 deriv->data.F32[1] = +r; 1616 deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); 1617 deriv->data.F32[3] = q * 1618 (2.0*py*params->data.F32[5] + params->data.F32[6]*X); 1619 deriv->data.F32[4] = -2.0*q*px*X; 1620 deriv->data.F32[5] = -2.0*q*py*Y; 1621 deriv->data.F32[6] = -q*X*Y; 1622 1623 return(f); 1624 } 1625 1626 /****************************************************************************** 1627 p_pmMinLM_PsuedoGauss2D_Vec(*deriv, *params, *x): this function wraps the 1628 above function in a form that is usable in the LM minimization routines. 1629 *****************************************************************************/ 1630 psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv, 1541 1631 psVector *params, 1542 1632 psArray *x) 1543 { 1544 PS_IMAGE_CHECK_NULL(deriv, NULL); 1545 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1546 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1547 PS_VECTOR_CHECK_NULL(params, NULL); 1548 PS_VECTOR_CHECK_EMPTY(params, NULL); 1549 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1550 PS_PTR_CHECK_NULL(x, NULL); 1551 if (deriv->numRows != x->n) { 1552 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1553 } 1554 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1555 // XXX: use static memory here. 1556 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1557 1558 for (psS32 i = 0 ; i < x->n ; i++) { 1559 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1560 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1561 } 1562 1563 tmpVec->data.F32[i] = pmMinLM_Gauss2D(tmpRow, 1564 params, 1565 (psVector *) x->data[i]); 1566 1567 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1568 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1569 } 1570 } 1571 1572 psFree(tmpRow); 1573 return(tmpVec); 1574 } 1575 1576 1577 /****************************************************************************** 1578 params->data.F32[0] = So; 1579 params->data.F32[1] = Zo; 1580 params->data.F32[2] = Xo; 1581 params->data.F32[3] = Yo; 1582 params->data.F32[4] = sqrt(2) / SigmaX; 1583 params->data.F32[5] = sqrt(2) / SigmaY; 1584 params->data.F32[6] = Sxy; 1585 *****************************************************************************/ 1586 psF32 pmMinLM_PsuedoGauss2D(psVector *deriv, 1587 psVector *params, 1588 psVector *x) 1589 { 1590 PS_VECTOR_CHECK_NULL(deriv, NAN); 1591 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1592 PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); 1593 PS_VECTOR_CHECK_NULL(params, NAN); 1594 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1595 PS_VECTOR_CHECK_SIZE(params, 7, NAN); 1596 PS_VECTOR_CHECK_NULL(x, NAN); 1597 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1598 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1599 1600 psF32 X = x->data.F32[0] - params->data.F32[2]; 1601 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1602 psF32 px = params->data.F32[4]*X; 1603 psF32 py = params->data.F32[5]*Y; 1604 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1605 psF32 t = 1 + z + 0.5*z*z; 1606 psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */ 1607 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1608 1609 // 1610 // note difference from a pure gaussian: q = params->data.F32[1]*r 1611 // 1612 1613 psF32 q = params->data.F32[1]*r*r*t; 1614 deriv->data.F32[0] = +1.0; 1615 deriv->data.F32[1] = +r; 1616 deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); 1617 deriv->data.F32[3] = q * 1618 (2.0*py*params->data.F32[5] + params->data.F32[6]*X); 1619 deriv->data.F32[4] = -2.0*q*px*X; 1620 deriv->data.F32[5] = -2.0*q*py*Y; 1621 deriv->data.F32[6] = -q*X*Y; 1622 1623 return(f); 1624 } 1625 1626 /****************************************************************************** 1627 p_pmMinLM_PsuedoGauss2D_Vec(*deriv, *params, *x): this function wraps the 1628 above function in a form that is usable in the LM minimization routines. 1629 *****************************************************************************/ 1630 psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv, 1631 psVector *params, 1632 psArray *x) 1633 { 1634 PS_IMAGE_CHECK_NULL(deriv, NULL); 1635 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1636 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1637 PS_VECTOR_CHECK_NULL(params, NULL); 1638 PS_VECTOR_CHECK_EMPTY(params, NULL); 1639 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1640 PS_PTR_CHECK_NULL(x, NULL); 1641 if (deriv->numRows != x->n) { 1642 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1643 } 1644 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1645 // XXX: use static memory here. 1646 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1647 1648 for (psS32 i = 0 ; i < x->n ; i++) { 1649 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1650 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1651 } 1652 1653 tmpVec->data.F32[i] = pmMinLM_PsuedoGauss2D(tmpRow, 1654 params, 1655 (psVector *) x->data[i]); 1656 1657 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1658 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1659 } 1660 } 1661 1662 psFree(tmpRow); 1663 return(tmpVec); 1664 } 1665 1666 1667 1668 1669 /****************************************************************************** 1670 params->data.F32[0] = So; 1671 params->data.F32[1] = Zo; 1672 params->data.F32[2] = Xo; 1673 params->data.F32[3] = Yo; 1674 params->data.F32[4] = Sx; 1675 params->data.F32[5] = Sy; 1676 params->data.F32[6] = Sxy; 1677 params->data.F32[7] = B2; 1678 params->data.F32[8] = B3; 1679 *****************************************************************************/ 1680 psF32 pmMinLM_Wauss2D(psVector *deriv, 1681 psVector *params, 1682 psVector *x) 1683 { 1684 PS_VECTOR_CHECK_NULL(deriv, NAN); 1685 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1686 PS_VECTOR_CHECK_SIZE(deriv, 9, NAN); 1687 PS_VECTOR_CHECK_NULL(params, NAN); 1688 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1689 PS_VECTOR_CHECK_SIZE(params, 9, NAN); 1690 PS_VECTOR_CHECK_NULL(x, NAN); 1691 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1692 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1693 1694 psF32 X = x->data.F32[0] - params->data.F32[2]; 1695 psF32 Y = x->data.F32[1] - params->data.F32[2]; 1696 psF32 px = params->data.F32[4]*X; 1697 psF32 py = params->data.F32[5]*Y; 1698 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1699 psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0); 1700 psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */ 1701 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1702 1703 // 1704 // note difference from gaussian: q = params->data.F32[1]*r 1705 // 1706 psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0)); 1707 deriv->data.F32[0] = +1.0; 1708 deriv->data.F32[1] = +r; 1709 deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); 1710 deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); 1711 deriv->data.F32[4] = -2.0*q*px*X; 1712 deriv->data.F32[5] = -2.0*q*py*Y; 1713 deriv->data.F32[6] = -q*X*Y; 1714 deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t; 1715 deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0; 1716 // 1717 // The values of 100 dampen the swing of params->data.F32[7,8] */ 1718 // 1719 1720 return(f); 1721 } 1722 1723 /****************************************************************************** 1724 p_pmMinLM_Wauss2D_Vec(*deriv, *params, *x): this function wraps the above 1725 function in a form that is usable in the LM minimization routines. 1726 *****************************************************************************/ 1727 psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv, 1728 psVector *params, 1729 psArray *x) 1730 { 1731 PS_IMAGE_CHECK_NULL(deriv, NULL); 1732 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1733 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1734 PS_VECTOR_CHECK_NULL(params, NULL); 1735 PS_VECTOR_CHECK_EMPTY(params, NULL); 1736 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1737 PS_PTR_CHECK_NULL(x, NULL); 1738 if (deriv->numRows != x->n) { 1739 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1740 } 1741 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1742 // XXX: use static memory here. 1743 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1744 1745 for (psS32 i = 0 ; i < x->n ; i++) { 1746 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1747 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1748 } 1749 1750 tmpVec->data.F32[i] = pmMinLM_Wauss2D(tmpRow, 1751 params, 1752 (psVector *) x->data[i]); 1753 1754 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1755 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1756 } 1757 } 1758 1759 psFree(tmpRow); 1760 return(tmpVec); 1761 } 1762 1763 1764 1765 1766 1767 1768 // XXX: What should these be? 1769 #define FFACTOR 1.0 1770 #define FSCALE 1.0 1771 /****************************************************************************** 1772 params->data.F32[0] = So; 1773 params->data.F32[1] = Zo; 1774 params->data.F32[2] = Xo; 1775 params->data.F32[3] = Yo; 1776 params->data.F32[4] = SxInner; 1777 params->data.F32[5] = SyInner; 1778 params->data.F32[6] = SxyInner; 1779 params->data.F32[7] = SxOuter; 1780 params->data.F32[8] = SyOuter; 1781 params->data.F32[9] = SxyOuter; 1782 params->data.F32[10] = N; 1783 *****************************************************************************/ 1784 psF32 pmMinLM_TwistGauss2D(psVector *deriv, 1785 psVector *params, 1786 psVector *x) 1787 { 1788 PS_VECTOR_CHECK_NULL(deriv, NAN); 1789 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1790 PS_VECTOR_CHECK_SIZE(deriv, 11, NAN); 1791 PS_VECTOR_CHECK_NULL(params, NAN); 1792 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1793 PS_VECTOR_CHECK_SIZE(params, 11, NAN); 1794 PS_VECTOR_CHECK_NULL(x, NAN); 1795 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1796 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1797 1798 psF32 X = x->data.F32[0] - params->data.F32[2]; 1799 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1800 psF32 px1 = params->data.F32[4]*X; 1801 psF32 py1 = params->data.F32[5]*Y; 1802 psF32 px2 = params->data.F32[7]*X; 1803 psF32 py2 = params->data.F32[8]*Y; 1804 psF32 z1 = 0.5*PS_SQR(px1) + 0.5*PS_SQR(py1) + params->data.F32[4]*X*Y; 1805 psF32 z2 = 0.5*PS_SQR(px2) + 0.5*PS_SQR(py2) + params->data.F32[9]*X*Y; 1806 psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10])); 1807 1808 1809 psF32 f = params->data.F32[5]*r + params->data.F32[6]; 1810 psF32 q1 = params->data.F32[5]*PS_SQR(r); 1811 psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0)); 1812 deriv->data.F32[0] = +1.0; 1813 deriv->data.F32[1] = +r; 1814 deriv->data.F32[2] = q1*(2.0*px1*params->data.F32[4] + params->data.F32[6]*Y) + q2*(2*px2*params->data.F32[7] + params->data.F32[9]*Y); 1815 deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X); 1816 1817 // 1818 // These fudge factors impede the growth of params->data.F32[4] beyond 1819 // params->data.F32[7]. 1820 // 1821 psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]); 1822 psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; 1823 deriv->data.F32[4] = -2.0*q1*px1*X*f2; 1824 1825 // 1826 // These fudge factors impede the growth of params->data.F32[5] beyond 1827 // params->data.F32[8]. 1828 // 1829 f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]); 1830 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; 1831 deriv->data.F32[5] = -2.0*q1*py1*Y*f2; 1832 deriv->data.F32[6] = -q1*X*Y; 1833 deriv->data.F32[7] = -2.0*q2*px2*X; 1834 deriv->data.F32[8] = -2.0*q2*py2*Y; 1835 deriv->data.F32[9] = -q2*X*Y; 1836 deriv->data.F32[10] = -q1*log(z2); 1837 1838 return(f); 1839 } 1840 1841 /****************************************************************************** 1842 p_pmMinLM_TwistGauss2D_Vec(*deriv, *params, *x): this function wraps the above 1843 function in a form that is usable in the LM minimization routines. 1844 *****************************************************************************/ 1845 psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv, 1846 psVector *params, 1847 psArray *x) 1848 { 1849 PS_IMAGE_CHECK_NULL(deriv, NULL); 1850 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1851 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1852 PS_VECTOR_CHECK_NULL(params, NULL); 1853 PS_VECTOR_CHECK_EMPTY(params, NULL); 1854 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1855 PS_PTR_CHECK_NULL(x, NULL); 1856 if (deriv->numRows != x->n) { 1857 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1858 } 1859 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1860 // XXX: use static memory here. 1861 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1862 1863 for (psS32 i = 0 ; i < x->n ; i++) { 1864 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1865 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1866 } 1867 1868 tmpVec->data.F32[i] = pmMinLM_TwistGauss2D(tmpRow, 1869 params, 1870 (psVector *) x->data[i]); 1871 1872 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1873 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1874 } 1875 } 1876 1877 psFree(tmpRow); 1878 return(tmpVec); 1879 } 1880 1881 1882 1883 /****************************************************************************** 1884 float Sersic() 1885 params->data.F32[0] = So; 1886 params->data.F32[1] = Zo; 1887 params->data.F32[2] = Xo; 1888 params->data.F32[3] = Yo; 1889 params->data.F32[4] = Sx; 1890 params->data.F32[5] = Sy; 1891 params->data.F32[6] = Sxy; 1892 params->data.F32[7] = Nexp; 1893 *****************************************************************************/ 1894 psF32 pmMinLM_Sersic(psVector *deriv, 1633 { 1634 PS_IMAGE_CHECK_NULL(deriv, NULL); 1635 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1636 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1637 PS_VECTOR_CHECK_NULL(params, NULL); 1638 PS_VECTOR_CHECK_EMPTY(params, NULL); 1639 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1640 PS_PTR_CHECK_NULL(x, NULL); 1641 if (deriv->numRows != x->n) { 1642 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1643 } 1644 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1645 // XXX: use static memory here. 1646 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1647 1648 for (psS32 i = 0 ; i < x->n ; i++) { 1649 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1650 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1651 } 1652 1653 tmpVec->data.F32[i] = pmMinLM_PsuedoGauss2D(tmpRow, 1654 params, 1655 (psVector *) x->data[i]); 1656 1657 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1658 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1659 } 1660 } 1661 1662 psFree(tmpRow); 1663 return(tmpVec); 1664 } 1665 1666 1667 1668 1669 /****************************************************************************** 1670 params->data.F32[0] = So; 1671 params->data.F32[1] = Zo; 1672 params->data.F32[2] = Xo; 1673 params->data.F32[3] = Yo; 1674 params->data.F32[4] = Sx; 1675 params->data.F32[5] = Sy; 1676 params->data.F32[6] = Sxy; 1677 params->data.F32[7] = B2; 1678 params->data.F32[8] = B3; 1679 *****************************************************************************/ 1680 psF32 pmMinLM_Wauss2D(psVector *deriv, 1681 psVector *params, 1682 psVector *x) 1683 { 1684 PS_VECTOR_CHECK_NULL(deriv, NAN); 1685 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1686 PS_VECTOR_CHECK_SIZE(deriv, 9, NAN); 1687 PS_VECTOR_CHECK_NULL(params, NAN); 1688 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1689 PS_VECTOR_CHECK_SIZE(params, 9, NAN); 1690 PS_VECTOR_CHECK_NULL(x, NAN); 1691 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1692 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1693 1694 psF32 X = x->data.F32[0] - params->data.F32[2]; 1695 psF32 Y = x->data.F32[1] - params->data.F32[2]; 1696 psF32 px = params->data.F32[4]*X; 1697 psF32 py = params->data.F32[5]*Y; 1698 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1699 psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0); 1700 psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */ 1701 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1702 1703 // 1704 // note difference from gaussian: q = params->data.F32[1]*r 1705 // 1706 psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0)); 1707 deriv->data.F32[0] = +1.0; 1708 deriv->data.F32[1] = +r; 1709 deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); 1710 deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); 1711 deriv->data.F32[4] = -2.0*q*px*X; 1712 deriv->data.F32[5] = -2.0*q*py*Y; 1713 deriv->data.F32[6] = -q*X*Y; 1714 deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t; 1715 deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0; 1716 // 1717 // The values of 100 dampen the swing of params->data.F32[7,8] */ 1718 // 1719 1720 return(f); 1721 } 1722 1723 /****************************************************************************** 1724 p_pmMinLM_Wauss2D_Vec(*deriv, *params, *x): this function wraps the above 1725 function in a form that is usable in the LM minimization routines. 1726 *****************************************************************************/ 1727 psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv, 1728 psVector *params, 1729 psArray *x) 1730 { 1731 PS_IMAGE_CHECK_NULL(deriv, NULL); 1732 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1733 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1734 PS_VECTOR_CHECK_NULL(params, NULL); 1735 PS_VECTOR_CHECK_EMPTY(params, NULL); 1736 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1737 PS_PTR_CHECK_NULL(x, NULL); 1738 if (deriv->numRows != x->n) { 1739 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1740 } 1741 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1742 // XXX: use static memory here. 1743 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1744 1745 for (psS32 i = 0 ; i < x->n ; i++) { 1746 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1747 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1748 } 1749 1750 tmpVec->data.F32[i] = pmMinLM_Wauss2D(tmpRow, 1751 params, 1752 (psVector *) x->data[i]); 1753 1754 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1755 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1756 } 1757 } 1758 1759 psFree(tmpRow); 1760 return(tmpVec); 1761 } 1762 1763 1764 1765 1766 1767 1768 // XXX: What should these be? 1769 #define FFACTOR 1.0 1770 #define FSCALE 1.0 1771 /****************************************************************************** 1772 params->data.F32[0] = So; 1773 params->data.F32[1] = Zo; 1774 params->data.F32[2] = Xo; 1775 params->data.F32[3] = Yo; 1776 params->data.F32[4] = SxInner; 1777 params->data.F32[5] = SyInner; 1778 params->data.F32[6] = SxyInner; 1779 params->data.F32[7] = SxOuter; 1780 params->data.F32[8] = SyOuter; 1781 params->data.F32[9] = SxyOuter; 1782 params->data.F32[10] = N; 1783 *****************************************************************************/ 1784 psF32 pmMinLM_TwistGauss2D(psVector *deriv, 1895 1785 psVector *params, 1896 1786 psVector *x) 1897 { 1898 PS_VECTOR_CHECK_NULL(deriv, NAN); 1899 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1900 PS_VECTOR_CHECK_SIZE(deriv, 8, NAN); 1901 PS_VECTOR_CHECK_NULL(params, NAN); 1902 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1903 PS_VECTOR_CHECK_SIZE(params, 8, NAN); 1904 PS_VECTOR_CHECK_NULL(x, NAN); 1905 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1906 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1907 1908 psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); 1909 return(0.0); 1910 } 1911 /****************************************************************************** 1912 p_pmMinLM_Sersic_Vec(*deriv, *params, *x): this function wraps the above 1913 function in a form that is usable in the LM minimization routines. 1914 *****************************************************************************/ 1915 psVector *p_pmMinLM_Sersic_Vec(psImage *deriv, 1787 { 1788 PS_VECTOR_CHECK_NULL(deriv, NAN); 1789 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1790 PS_VECTOR_CHECK_SIZE(deriv, 11, NAN); 1791 PS_VECTOR_CHECK_NULL(params, NAN); 1792 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1793 PS_VECTOR_CHECK_SIZE(params, 11, NAN); 1794 PS_VECTOR_CHECK_NULL(x, NAN); 1795 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1796 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1797 1798 psF32 X = x->data.F32[0] - params->data.F32[2]; 1799 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1800 psF32 px1 = params->data.F32[4]*X; 1801 psF32 py1 = params->data.F32[5]*Y; 1802 psF32 px2 = params->data.F32[7]*X; 1803 psF32 py2 = params->data.F32[8]*Y; 1804 psF32 z1 = 0.5*PS_SQR(px1) + 0.5*PS_SQR(py1) + params->data.F32[4]*X*Y; 1805 psF32 z2 = 0.5*PS_SQR(px2) + 0.5*PS_SQR(py2) + params->data.F32[9]*X*Y; 1806 psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10])); 1807 1808 1809 psF32 f = params->data.F32[5]*r + params->data.F32[6]; 1810 psF32 q1 = params->data.F32[5]*PS_SQR(r); 1811 psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0)); 1812 deriv->data.F32[0] = +1.0; 1813 deriv->data.F32[1] = +r; 1814 deriv->data.F32[2] = q1*(2.0*px1*params->data.F32[4] + params->data.F32[6]*Y) + q2*(2*px2*params->data.F32[7] + params->data.F32[9]*Y); 1815 deriv->data.F32[3] = q1*(2.0*py1*params->data.F32[5] + params->data.F32[6]*X) + q2*(2*py2*params->data.F32[8] + params->data.F32[9]*X); 1816 1817 // 1818 // These fudge factors impede the growth of params->data.F32[4] beyond 1819 // params->data.F32[7]. 1820 // 1821 psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]); 1822 psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; 1823 deriv->data.F32[4] = -2.0*q1*px1*X*f2; 1824 1825 // 1826 // These fudge factors impede the growth of params->data.F32[5] beyond 1827 // params->data.F32[8]. 1828 // 1829 f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]); 1830 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; 1831 deriv->data.F32[5] = -2.0*q1*py1*Y*f2; 1832 deriv->data.F32[6] = -q1*X*Y; 1833 deriv->data.F32[7] = -2.0*q2*px2*X; 1834 deriv->data.F32[8] = -2.0*q2*py2*Y; 1835 deriv->data.F32[9] = -q2*X*Y; 1836 deriv->data.F32[10] = -q1*log(z2); 1837 1838 return(f); 1839 } 1840 1841 /****************************************************************************** 1842 p_pmMinLM_TwistGauss2D_Vec(*deriv, *params, *x): this function wraps the above 1843 function in a form that is usable in the LM minimization routines. 1844 *****************************************************************************/ 1845 psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv, 1916 1846 psVector *params, 1917 1847 psArray *x) 1918 { 1919 PS_IMAGE_CHECK_NULL(deriv, NULL); 1920 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1921 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1922 PS_VECTOR_CHECK_NULL(params, NULL); 1923 PS_VECTOR_CHECK_EMPTY(params, NULL); 1924 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1925 PS_PTR_CHECK_NULL(x, NULL); 1926 if (deriv->numRows != x->n) { 1927 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1928 } 1929 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1930 // XXX: use static memory here. 1931 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1932 1933 for (psS32 i = 0 ; i < x->n ; i++) { 1934 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1935 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1936 } 1937 1938 tmpVec->data.F32[i] = pmMinLM_Sersic(tmpRow, 1939 params, 1940 (psVector *) x->data[i]); 1941 1942 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1943 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1944 } 1945 } 1946 1947 psFree(tmpRow); 1948 return(tmpVec); 1949 } 1950 1951 /****************************************************************************** 1952 float SersicBulge() 1953 params->data.F32[0] So; 1954 params->data.F32[1] Zo; 1955 params->data.F32[2] Xo; 1956 params->data.F32[3] Yo; 1957 params->data.F32[4] SxInner; 1958 params->data.F32[5] SyInner; 1959 params->data.F32[6] SxyInner; 1960 params->data.F32[7] Zd; 1961 params->data.F32[8] SxOuter; 1962 params->data.F32[9] SyOuter; 1963 params->data.F32[10] = SxyOuter; 1964 params->data.F32[11] = Nexp; 1965 *****************************************************************************/ 1966 psF32 pmMinLM_SersicCore(psVector *deriv, 1848 { 1849 PS_IMAGE_CHECK_NULL(deriv, NULL); 1850 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1851 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1852 PS_VECTOR_CHECK_NULL(params, NULL); 1853 PS_VECTOR_CHECK_EMPTY(params, NULL); 1854 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1855 PS_PTR_CHECK_NULL(x, NULL); 1856 if (deriv->numRows != x->n) { 1857 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1858 } 1859 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1860 // XXX: use static memory here. 1861 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1862 1863 for (psS32 i = 0 ; i < x->n ; i++) { 1864 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1865 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1866 } 1867 1868 tmpVec->data.F32[i] = pmMinLM_TwistGauss2D(tmpRow, 1869 params, 1870 (psVector *) x->data[i]); 1871 1872 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1873 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1874 } 1875 } 1876 1877 psFree(tmpRow); 1878 return(tmpVec); 1879 } 1880 1881 1882 1883 /****************************************************************************** 1884 float Sersic() 1885 params->data.F32[0] = So; 1886 params->data.F32[1] = Zo; 1887 params->data.F32[2] = Xo; 1888 params->data.F32[3] = Yo; 1889 params->data.F32[4] = Sx; 1890 params->data.F32[5] = Sy; 1891 params->data.F32[6] = Sxy; 1892 params->data.F32[7] = Nexp; 1893 *****************************************************************************/ 1894 psF32 pmMinLM_Sersic(psVector *deriv, 1895 psVector *params, 1896 psVector *x) 1897 { 1898 PS_VECTOR_CHECK_NULL(deriv, NAN); 1899 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1900 PS_VECTOR_CHECK_SIZE(deriv, 8, NAN); 1901 PS_VECTOR_CHECK_NULL(params, NAN); 1902 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1903 PS_VECTOR_CHECK_SIZE(params, 8, NAN); 1904 PS_VECTOR_CHECK_NULL(x, NAN); 1905 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1906 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1907 1908 psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); 1909 return(0.0); 1910 } 1911 /****************************************************************************** 1912 p_pmMinLM_Sersic_Vec(*deriv, *params, *x): this function wraps the above 1913 function in a form that is usable in the LM minimization routines. 1914 *****************************************************************************/ 1915 psVector *p_pmMinLM_Sersic_Vec(psImage *deriv, 1967 1916 psVector *params, 1968 psVector *x) 1969 { 1970 PS_VECTOR_CHECK_NULL(deriv, NAN); 1971 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1972 PS_VECTOR_CHECK_SIZE(deriv, 12, NAN); 1973 PS_VECTOR_CHECK_NULL(params, NAN); 1974 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1975 PS_VECTOR_CHECK_SIZE(params, 12, NAN); 1976 PS_VECTOR_CHECK_NULL(x, NAN); 1977 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1978 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1979 1980 psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); 1981 return(0.0); 1982 } 1983 /****************************************************************************** 1984 p_pmMinLM_SersicCore_Vec(*deriv, *params, *x): this function wraps the above 1985 function in a form that is usable in the LM minimization routines. 1986 *****************************************************************************/ 1987 psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, 1988 psVector *params, 1989 psArray *x) 1990 { 1991 PS_IMAGE_CHECK_NULL(deriv, NULL); 1992 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1993 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1994 PS_VECTOR_CHECK_NULL(params, NULL); 1995 PS_VECTOR_CHECK_EMPTY(params, NULL); 1996 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1997 PS_PTR_CHECK_NULL(x, NULL); 1998 if (deriv->numRows != x->n) { 1999 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 2000 } 2001 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 2002 // XXX: use static memory here. 2003 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 2004 2005 for (psS32 i = 0 ; i < x->n ; i++) { 2006 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 2007 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 2008 } 2009 2010 tmpVec->data.F32[i] = pmMinLM_SersicCore(tmpRow, 2011 params, 2012 (psVector *) x->data[i]); 2013 2014 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 2015 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 2016 } 2017 } 2018 2019 psFree(tmpRow); 2020 return(tmpVec); 2021 } 2022 2023 2024 2025 /****************************************************************************** 2026 *****************************************************************************/ 2027 psF32 pmMinLM_PsuedoSersic(psVector *deriv, 2028 psVector *params, 2029 psVector *x) 2030 { 2031 return(0.0); 2032 } 1917 psArray *x) 1918 { 1919 PS_IMAGE_CHECK_NULL(deriv, NULL); 1920 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1921 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1922 PS_VECTOR_CHECK_NULL(params, NULL); 1923 PS_VECTOR_CHECK_EMPTY(params, NULL); 1924 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1925 PS_PTR_CHECK_NULL(x, NULL); 1926 if (deriv->numRows != x->n) { 1927 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1928 } 1929 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1930 // XXX: use static memory here. 1931 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1932 1933 for (psS32 i = 0 ; i < x->n ; i++) { 1934 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1935 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1936 } 1937 1938 tmpVec->data.F32[i] = pmMinLM_Sersic(tmpRow, 1939 params, 1940 (psVector *) x->data[i]); 1941 1942 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1943 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1944 } 1945 } 1946 1947 psFree(tmpRow); 1948 return(tmpVec); 1949 } 1950 1951 /****************************************************************************** 1952 float SersicBulge() 1953 params->data.F32[0] So; 1954 params->data.F32[1] Zo; 1955 params->data.F32[2] Xo; 1956 params->data.F32[3] Yo; 1957 params->data.F32[4] SxInner; 1958 params->data.F32[5] SyInner; 1959 params->data.F32[6] SxyInner; 1960 params->data.F32[7] Zd; 1961 params->data.F32[8] SxOuter; 1962 params->data.F32[9] SyOuter; 1963 params->data.F32[10] = SxyOuter; 1964 params->data.F32[11] = Nexp; 1965 *****************************************************************************/ 1966 psF32 pmMinLM_SersicCore(psVector *deriv, 1967 psVector *params, 1968 psVector *x) 1969 { 1970 PS_VECTOR_CHECK_NULL(deriv, NAN); 1971 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1972 PS_VECTOR_CHECK_SIZE(deriv, 12, NAN); 1973 PS_VECTOR_CHECK_NULL(params, NAN); 1974 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1975 PS_VECTOR_CHECK_SIZE(params, 12, NAN); 1976 PS_VECTOR_CHECK_NULL(x, NAN); 1977 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1978 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1979 1980 psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); 1981 return(0.0); 1982 } 1983 /****************************************************************************** 1984 p_pmMinLM_SersicCore_Vec(*deriv, *params, *x): this function wraps the above 1985 function in a form that is usable in the LM minimization routines. 1986 *****************************************************************************/ 1987 psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, 1988 psVector *params, 1989 psArray *x) 1990 { 1991 PS_IMAGE_CHECK_NULL(deriv, NULL); 1992 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1993 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1994 PS_VECTOR_CHECK_NULL(params, NULL); 1995 PS_VECTOR_CHECK_EMPTY(params, NULL); 1996 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1997 PS_PTR_CHECK_NULL(x, NULL); 1998 if (deriv->numRows != x->n) { 1999 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 2000 } 2001 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 2002 // XXX: use static memory here. 2003 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 2004 2005 for (psS32 i = 0 ; i < x->n ; i++) { 2006 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 2007 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 2008 } 2009 2010 tmpVec->data.F32[i] = pmMinLM_SersicCore(tmpRow, 2011 params, 2012 (psVector *) x->data[i]); 2013 2014 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 2015 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 2016 } 2017 } 2018 2019 psFree(tmpRow); 2020 return(tmpVec); 2021 } 2022 2023 2024 2025 /****************************************************************************** 2026 *****************************************************************************/ 2027 psF32 pmMinLM_PsuedoSersic(psVector *deriv, 2028 psVector *params, 2029 psVector *x) 2030 { 2031 return(0.0); 2032 }
Note:
See TracChangeset
for help on using the changeset viewer.
