Changeset 3510
- Timestamp:
- Mar 25, 2005, 1:18:35 PM (21 years ago)
- Location:
- trunk/psModules/src
- Files:
-
- 2 edited
-
pmObjects.c (modified) (4 diffs)
-
pmObjects.h (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/pmObjects.c
r3498 r3510 5 5 * @author GLG, MHPCC 6 6 * 7 * @version $Revision: 1. 6$ $Name: not supported by cvs2svn $8 * @date $Date: 2005-03-2 4 22:36:52$7 * @version $Revision: 1.7 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2005-03-25 23:18:35 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 740 740 float SNlim, 741 741 const psRegion *valid) 742 XXX: Waiting for sample code from IfA.743 744 XXX: Most code this.742 XXX: The sigX and sigY stuff in the SDRS is unclear. 743 744 XXX: How can this function ever return FALSE? 745 745 *****************************************************************************/ 746 746 #define SATURATE 0.0 … … 760 760 psSource *tmpSrc = (psSource *) source->data[i]; 761 761 PS_PTR_CHECK_NULL(tmpSrc->moments, false); 762 tmpSrc->peak->class = 0; 762 763 763 764 if (tmpSrc->moments->Peak > SATURATE) { 764 tmpSrc->peak->class 765 tmpSrc->peak->class|= PS_SOURCE_SATURATED; 765 766 } else { 766 767 // XXX: gleen these from the metadata: keywords GAIN and READ_NOISE. … … 773 774 PS_SQRT_F32(S + (A * B) + ((A * readNoise * readNoise) / PS_SQRT_F32(gain))); 774 775 if (SN < FAINT_SN_LIM) { 775 tmpSrc->peak->class 776 tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR; 776 777 } 777 778 if (SN < PSF_SN_LIM) { 778 tmpSrc->peak->class 779 tmpSrc->peak->class|= PS_SOURCE_FAINTSTAR; 779 780 } 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 // 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; 780 797 } 781 } 782 783 return(rc); 798 799 if ((sigX > (clumpX + clumpDX)) && 800 (sigY > (clumpY + clumpDY))) 801 tmpSrc->peak->class|= PS_SOURCE_GALAXY; 802 } 803 804 if (tmpSrc->peak->class == 0) { 805 tmpSrc->peak->class|= PS_SOURCE_OTHER; 806 } 784 807 } 785 786 787 788 /******************************************************************************789 pmSourceSetPixelCircle(source, image, radius)790 791 XXX: Why boolean output?792 793 XXX: Why are we checking source->moments for NULL? Should the circle be794 centered on the centroid or the peak?795 796 XXX: The circle will have a diameter of (1+radius). This is different from797 the pmSourceSetLocal() function.798 *****************************************************************************/799 bool pmSourceSetPixelCircle(psSource *source,800 const psImage *image,801 psF32 radius)802 {803 PS_IMAGE_CHECK_NULL(image, false);804 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false);805 PS_PTR_CHECK_NULL(source, false);806 // PS_PTR_CHECK_NULL(source->moments, false);807 PS_PTR_CHECK_NULL(source->peak, false);808 PS_FLOAT_COMPARE(0.0, radius, false);809 810 //811 // We define variables for code readability.812 //813 psS32 radiusS32 = (psS32) radius;814 psS32 SubImageCenterRow = source->peak->y;815 psS32 SubImageCenterCol = source->peak->x;816 psS32 SubImageStartRow = SubImageCenterRow - radiusS32;817 psS32 SubImageEndRow = SubImageCenterRow + radiusS32;818 psS32 SubImageStartCol = SubImageCenterCol - radiusS32;819 psS32 SubImageEndCol = SubImageCenterCol + radiusS32;820 821 if (SubImageStartRow < 0) {822 psError(PS_ERR_UNKNOWN, true, "Sub image startRow is outside image boundaries (%d).\n",823 SubImageStartRow);824 return(false);825 }826 if (SubImageEndRow+1 >= image->numRows) {827 psError(PS_ERR_UNKNOWN, true, "Sub image endRow is outside image boundaries (%d).\n",828 SubImageEndRow);829 return(false);830 }831 if (SubImageStartCol < 0) {832 psError(PS_ERR_UNKNOWN, true, "Sub image startCol is outside image boundaries (%d).\n",833 SubImageStartCol);834 return(false);835 }836 if (SubImageEndCol+1 >= image->numCols) {837 psError(PS_ERR_UNKNOWN, true, "Sub image endCol is outside image boundaries (%d).\n",838 SubImageEndCol);839 return(false);840 }841 842 // XXX: Must recycle image.843 if (source->pixels != NULL) {844 psLogMsg(__func__, PS_LOG_WARN,845 "WARNING: pmSourceSetPixelCircle(): image->pixels not NULL. Freeing and reallocating.\n");846 psFree(source->pixels);847 }848 source->pixels = psImageSubset((psImage *) image,849 SubImageStartCol,850 SubImageStartRow,851 SubImageEndCol+1,852 SubImageEndRow+1);853 854 // XXX: Must recycle image.855 if (source->mask != NULL) {856 psFree(source->mask);857 }858 source->mask = psImageAlloc(1 + 2 * radiusS32, 1 + 2 * radiusS32, PS_TYPE_F32);859 860 //861 // Loop through the subimage mask, initialize mask to 0 or 1.862 //863 for (psS32 row = 0 ; row < source->mask->numRows; row++) {864 for (psS32 col = 0 ; col < source->mask->numCols; col++) {865 866 if (CheckRadius2((psF32) radiusS32,867 (psF32) radiusS32,868 radius,869 (psF32) col,870 (psF32) row)) {871 source->mask->data.U8[row][col] = 1;872 } else {873 source->mask->data.U8[row][col] = 1;874 }875 }876 }877 878 /*879 for (psS32 row = SubImageCenterRow - radiusS32; row <= SubImageCenterRow + radiusS32; row++) {880 for (psS32 col = SubImageCenterCol - radiusS32; col <= SubImageCenterCol + radiusS32; col++) {881 if (CheckRadius(source->peak, radius, (psF32) col, (psF32) row)) {882 source->mask->data.U8[row-SubImageCenterRow][col-SubImageCenterCol] = 1;883 }884 }885 }886 */887 888 return(true);889 808 } 890 809 891 892 /****************************************************************************** 893 pmSourceModelGuess(source, image, model): This function allocates a new 894 psModel structure and store it in the psSource data structure specified in the 895 argument list. The model type is specified in the argument list. The params 896 array in that psModel structure are allocated, and then set to the appropriate 897 values. This function returns true if everything was successful. 898 899 XXX: Many of the initial parameters are set to 0.0 since I don't know what 900 the appropiate initial guesses are. 901 *****************************************************************************/ 902 bool pmSourceModelGuess(psSource *source, 903 const psImage *image, 904 psModelType model) 905 { 906 PS_PTR_CHECK_NULL(source, false); 907 PS_PTR_CHECK_NULL(source->moments, false); 908 PS_PTR_CHECK_NULL(source->peak, false); 909 PS_IMAGE_CHECK_NULL(image, false); 910 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); 911 if (source->models != NULL) { 912 psLogMsg(__func__, PS_LOG_WARN, "WARNING: source->models was non-NULL; calling psFree(source->models).\n"); 913 psFree(source->models); 914 } 915 source->models = pmModelAlloc(PS_MODEL_UNDEFINED); 916 917 switch (model) { 918 case PS_MODEL_GAUSS: 919 source->models->type = PS_MODEL_GAUSS; 920 source->models->Nparams = 7; 921 source->models->params = (psF32 *) psAlloc(7 * sizeof(psF32)); 922 source->models->dparams = (psF32 *) psAlloc(7 * sizeof(psF32)); 923 for (psS32 i = 0 ; i < 7 ; i++) { 924 source->models->params[i] = 0.0; 925 source->models->dparams[i] = 0.0; 926 } 927 source->models->params[0] = source->moments->Sky; 928 source->models->params[1] = source->peak->counts - source->moments->Sky; 929 source->models->params[2] = source->moments->x; 930 source->models->params[3] = source->moments->y; 931 source->models->params[4] = sqrt(2.0) / source->moments->Sx; 932 source->models->params[5] = sqrt(2.0) / source->moments->Sy; 933 source->models->params[6] = source->moments->Sxy; 934 source->models->chisq = 0.0; 935 return(true); 936 case PS_MODEL_PGAUSS: 937 source->models->type = PS_MODEL_PGAUSS; 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_TWIST_GAUSS: 955 source->models->type = PS_MODEL_TWIST_GAUSS; 956 source->models->Nparams = 11; 957 source->models->params = (psF32 *) psAlloc(11 * sizeof(psF32)); 958 source->models->dparams = (psF32 *) psAlloc(11 * sizeof(psF32)); 959 for (psS32 i = 0 ; i < 11 ; i++) { 960 source->models->params[i] = 0.0; 961 source->models->dparams[i] = 0.0; 962 } 963 964 source->models->params[0] = source->moments->Sky; 965 source->models->params[1] = source->peak->counts - source->moments->Sky; 966 source->models->params[2] = source->moments->x; 967 source->models->params[3] = source->moments->y; 968 // XXX: What are these? 969 // source->models->params[4] = SxInner; 970 // source->models->params[5] = SyInner; 971 // source->models->params[6] = SxyInner; 972 // source->models->params[7] = SxOuter; 973 // source->models->params[8] = SyOuter; 974 // source->models->params[9] = SxyOuter; 975 // source->models->params[10] = N; 976 977 source->models->chisq = 0.0; 978 return(true); 979 case PS_MODEL_WAUSS: 980 981 source->models->params[0] = source->moments->Sky; 982 source->models->params[1] = source->peak->counts - source->moments->Sky; 983 source->models->params[2] = source->moments->x; 984 source->models->params[3] = source->moments->y; 985 source->models->params[4] = sqrt(2.0) / source->moments->Sx; 986 source->models->params[5] = sqrt(2.0) / source->moments->Sy; 987 source->models->params[6] = source->moments->Sxy; 988 // XXX: What are these? 989 // source->models->params[7] = B2; 990 // source->models->params[8] = B3; 991 992 source->models->type = PS_MODEL_WAUSS; 993 source->models->Nparams = 9; 994 source->models->params = (psF32 *) psAlloc(9 * sizeof(psF32)); 995 source->models->dparams = (psF32 *) psAlloc(9 * sizeof(psF32)); 996 for (psS32 i = 0 ; i < 9 ; i++) { 997 source->models->params[i] = 0.0; 998 source->models->dparams[i] = 0.0; 999 } 1000 source->models->chisq = 0.0; 1001 return(true); 1002 case PS_MODEL_SERSIC: 1003 source->models->type = PS_MODEL_SERSIC; 1004 source->models->Nparams = 8; 1005 source->models->params = (psF32 *) psAlloc(8 * sizeof(psF32)); 1006 source->models->dparams = (psF32 *) psAlloc(8 * sizeof(psF32)); 1007 for (psS32 i = 0 ; i < 8 ; i++) { 1008 source->models->params[i] = 0.0; 1009 source->models->dparams[i] = 0.0; 1010 } 1011 1012 source->models->params[0] = source->moments->Sky; 1013 source->models->params[1] = source->peak->counts - source->moments->Sky; 1014 source->models->params[2] = source->moments->x; 1015 source->models->params[3] = source->moments->y; 1016 source->models->params[4] = sqrt(2.0) / source->moments->Sx; 1017 source->models->params[5] = sqrt(2.0) / source->moments->Sy; 1018 source->models->params[6] = source->moments->Sxy; 1019 // XXX: What are these? 1020 //source->models->params[7] = Nexp; 1021 1022 source->models->chisq = 0.0; 1023 return(true); 1024 case PS_MODEL_SERSIC_CORE: 1025 source->models->type = PS_MODEL_SERSIC_CORE; 1026 source->models->Nparams = 12; 1027 source->models->params = (psF32 *) psAlloc(12 * sizeof(psF32)); 1028 source->models->dparams = (psF32 *) psAlloc(12 * sizeof(psF32)); 1029 for (psS32 i = 0 ; i < 12 ; i++) { 1030 source->models->params[i] = 0.0; 1031 source->models->dparams[i] = 0.0; 1032 } 1033 1034 source->models->params[0] = source->moments->Sky; 1035 source->models->params[1] = source->peak->counts - source->moments->Sky; 1036 source->models->params[2] = source->moments->x; 1037 source->models->params[3] = source->moments->y; 1038 // XXX: What are these? 1039 //source->models->params[4] SxInner; 1040 //source->models->params[5] SyInner; 1041 //source->models->params[6] SxyInner; 1042 //source->models->params[7] Zd; 1043 //source->models->params[8] SxOuter; 1044 //source->models->params[9] SyOuter; 1045 //source->models->params[10] = SxyOuter; 1046 //source->models->params[11] = Nexp; 1047 1048 source->models->chisq = 0.0; 1049 return(true); 1050 default: 1051 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1052 return(false); 1053 } 1054 } 1055 1056 1057 /****************************************************************************** 1058 1059 *****************************************************************************/ 1060 psArray *pmSourceContour(psSource *source, 1061 const psImage *image, 1062 psF32 level, 1063 psS32 mode) 1064 { 1065 return(NULL); 1066 } 1067 1068 psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv, psVector *params, psArray *x); 1069 psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); 1070 psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv, psVector *params, psArray *x); 1071 psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv, psVector *params, psArray *x); 1072 psVector *p_pmMinLM_Sersic_Vec(psImage *deriv, psVector *params, psArray *x); 1073 psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, psVector *params, psArray *x); 1074 1075 //XXX: What should these values be? 1076 #define PM_SOURCE_FIT_MODEL_NUM_ITERATIONS 100 1077 #define PM_SOURCE_FIT_MODEL_TOLERANCE 1.0 1078 /****************************************************************************** 1079 pmSourceFitModel(source, image): must create the appropiate arguments to the 1080 LM minimization routines for the various p_pmMinLM_XXXXXX_Vec() functions. 1081 1082 XXX: should there be a mask value? 1083 *****************************************************************************/ 1084 bool pmSourceFitModel(psSource *source, 1085 const psImage *image) 1086 { 1087 PS_PTR_CHECK_NULL(source, false); 1088 PS_PTR_CHECK_NULL(source->moments, false); 1089 PS_PTR_CHECK_NULL(source->peak, false); 1090 PS_PTR_CHECK_NULL(source->pixels, false); 1091 PS_PTR_CHECK_NULL(source->models, false); 1092 PS_IMAGE_CHECK_NULL(image, false); 1093 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); 1094 1095 psBool rc; 1096 psS32 count = 0; 1097 for (psS32 i = 0 ; i < source->pixels->numRows ; i++) { 1098 for (psS32 j = 0 ; j < source->pixels->numCols ; j++) { 1099 if (source->mask->data.U8[i][j] == 0) { 1100 count++; 1101 } 1102 } 1103 } 1104 psArray *x = psArrayAlloc(count); 1105 psVector *y = psVectorAlloc(count, PS_TYPE_F32); 1106 for (psS32 i = 0 ; i < source->pixels->numRows ; i++) { 1107 for (psS32 j = 0 ; j < source->pixels->numCols ; j++) { 1108 if (source->mask->data.U8[i][j] == 0) { 1109 psVector *coord = psVectorAlloc(2, PS_TYPE_F32); 1110 // XXX: Should we use the subimage offsets here, or does it not matter? 1111 coord->data.F32[0] = (psF32) (i); 1112 coord->data.F32[1] = (psF32) (j); 1113 x->data[count] = (psPtr *) coord; 1114 y->data.F32[count] = source->pixels->data.F32[i][j]; 1115 } 1116 } 1117 } 1118 1119 psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS, 1120 PM_SOURCE_FIT_MODEL_TOLERANCE); 1121 1122 psVector *params = psVectorAlloc(source->models->Nparams, PS_TYPE_F32); 1123 1124 switch (source->models->type) { 1125 case PS_MODEL_GAUSS: 1126 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1127 NULL, (psMinimizeLMChi2Func) p_pmMinLM_Gauss2D_Vec); 1128 break; 1129 case PS_MODEL_PGAUSS: 1130 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1131 NULL, (psMinimizeLMChi2Func) p_pmMinLM_PsuedoGauss2D_Vec); 1132 break; 1133 case PS_MODEL_TWIST_GAUSS: 1134 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1135 NULL, (psMinimizeLMChi2Func) p_pmMinLM_Wauss2D_Vec); 1136 break; 1137 case PS_MODEL_WAUSS: 1138 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1139 NULL, (psMinimizeLMChi2Func) p_pmMinLM_TwistGauss2D_Vec); 1140 break; 1141 case PS_MODEL_SERSIC: 1142 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1143 NULL, (psMinimizeLMChi2Func) p_pmMinLM_Sersic_Vec); 1144 break; 1145 case PS_MODEL_SERSIC_CORE: 1146 rc = psMinimizeLMChi2(myMin, NULL, params, NULL, x, y, 1147 NULL, (psMinimizeLMChi2Func) p_pmMinLM_SersicCore_Vec); 1148 break; 1149 default: 1150 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1151 rc = false; 1152 } 1153 1154 psFree(x); 1155 psFree(y); 1156 psFree(myMin); 1157 psFree(params); 1158 return(rc); 1159 } 1160 1161 bool p_pmSourceAddOrSubModel(psImage *image, 1162 psSource *source, 1163 bool center, 1164 psS32 flag) 1165 { 1166 PS_PTR_CHECK_NULL(source, false); 1167 PS_PTR_CHECK_NULL(source->moments, false); 1168 PS_PTR_CHECK_NULL(source->peak, false); 1169 PS_PTR_CHECK_NULL(source->pixels, false); 1170 PS_PTR_CHECK_NULL(source->models, false); 1171 PS_IMAGE_CHECK_NULL(image, false); 1172 PS_IMAGE_CHECK_TYPE(image, PS_TYPE_F32, false); 1173 1174 // XXX: make static, document. 1175 psVector *deriv = psVectorAlloc(100, PS_TYPE_F32); 1176 psVector *params = psVectorAlloc(100, PS_TYPE_F32); 1177 psVector *x = psVectorAlloc(2, PS_TYPE_F32); 1178 for (psS32 i = 0 ; i < source->models->Nparams ; i++) { 1179 params->data.F32[i] = source->models->params[i]; 1180 } 1181 1182 for (psS32 i = 0 ; i < source->pixels->numRows ; i++) { 1183 for (psS32 j = 0 ; j < source->pixels->numCols ; j++) { 1184 psF32 pixelValue; 1185 // XXX: Should you use offsets here? 1186 // XXX: Make sure you have col/row order correct. 1187 x->data.F32[0] = (float) j; 1188 x->data.F32[1] = (float) i; 1189 switch (source->models->type) { 1190 case PS_MODEL_GAUSS: 1191 pixelValue = pmMinLM_Gauss2D(deriv, params, x); 1192 break; 1193 case PS_MODEL_PGAUSS: 1194 pixelValue = pmMinLM_PsuedoGauss2D(deriv, params, x); 1195 break; 1196 case PS_MODEL_TWIST_GAUSS: 1197 pixelValue = pmMinLM_TwistGauss2D(deriv, params, x); 1198 break; 1199 case PS_MODEL_WAUSS: 1200 pixelValue = pmMinLM_Wauss2D(deriv, params, x); 1201 break; 1202 case PS_MODEL_SERSIC: 1203 pixelValue = pmMinLM_Sersic(deriv, params, x); 1204 break; 1205 case PS_MODEL_SERSIC_CORE: 1206 pixelValue = pmMinLM_SersicCore(deriv, params, x); 1207 break; 1208 default: 1209 psError(PS_ERR_UNKNOWN, true, "Undefined psModelType"); 1210 psFree(x); 1211 psFree(deriv); 1212 psFree(params); 1213 return(false); 1214 } 1215 if (flag == 1) { 1216 pixelValue = -pixelValue; 1217 } 1218 1219 // XXX: Must figure out how to calculate the image coordinates and 1220 // how to use the boolean "center" flag. 1221 psS32 imageRow = 0; 1222 psS32 imageCol = 0; 1223 image->data.F32[imageRow][imageCol]+= pixelValue; 1224 } 1225 } 1226 psFree(x); 1227 psFree(deriv); 1228 psFree(params); 1229 return(true); 1230 } 1231 1232 1233 1234 /****************************************************************************** 1235 *****************************************************************************/ 1236 bool pmSourceAddModel(psImage *image, 1237 psSource *source, 1238 bool center) 1239 { 1240 return(p_pmSourceAddOrSubModel(image, source, center, 0)); 1241 } 1242 1243 /****************************************************************************** 1244 *****************************************************************************/ 1245 bool pmSourceSubModel(psImage *image, 1246 psSource *source, 1247 bool center) 1248 { 1249 return(p_pmSourceAddOrSubModel(image, source, center, 1)); 1250 } 1251 1252 1253 // XXX: Put this is psConstants.h 1254 #define PS_VECTOR_CHECK_SIZE(VEC1, N, RVAL) \ 1255 if (VEC1->n != N) { \ 1256 psError(PS_ERR_BAD_PARAMETER_SIZE, true, \ 1257 "psVector %s has size %d, should be %d.", \ 1258 #VEC1, VEC1->n, N); \ 1259 return(RVAL); \ 1260 } 1261 1262 1263 /****************************************************************************** 1264 pmMinLM_Gauss2D(*deriv, *params, *x): the argument "x" contains a single "x, 1265 y" coordinate pair. This function computes the gaussian, specified by the 1266 parameters in "params" at that x,y point and returns the value. The 1267 derivatives are also caculated and returned in the "deriv" argument. 1268 1269 params->data.F32[0] = So; 1270 params->data.F32[1] = Zo; 1271 params->data.F32[2] = Xo; 1272 params->data.F32[3] = Yo; 1273 params->data.F32[4] = sqrt(2.0) / SigmaX; 1274 params->data.F32[5] = sqrt(2.0) / SigmaY; 1275 params->data.F32[6] = Sxy; 1276 1277 XXX: Consider getting rid of the parameter checks since this might consume 1278 a significant fraction of this function CPU time. 1279 *****************************************************************************/ 1280 psF32 pmMinLM_Gauss2D(psVector *deriv, 1281 psVector *params, 1282 psVector *x) 1283 { 1284 PS_VECTOR_CHECK_NULL(deriv, NAN); 1285 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1286 PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); 1287 PS_VECTOR_CHECK_NULL(params, NAN); 1288 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1289 PS_VECTOR_CHECK_SIZE(params, 7, NAN); 1290 PS_VECTOR_CHECK_NULL(x, NAN); 1291 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1292 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1293 1294 psF32 X = x->data.F32[0] - params->data.F32[2]; 1295 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1296 psF32 px = params->data.F32[4]*X; 1297 psF32 py = params->data.F32[5]*Y; 1298 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1299 psF32 r = exp(-z); 1300 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1301 1302 psF32 q = params->data.F32[1]*r; 1303 deriv->data.F32[0] = +1.0; 1304 deriv->data.F32[1] = +r; 1305 deriv->data.F32[2] = q*(2*px*params->data.F32[4] + params->data.F32[6]*Y); 1306 deriv->data.F32[3] = q*(2*py*params->data.F32[5] + params->data.F32[6]*X); 1307 deriv->data.F32[4] = -2.0*q*px*X; 1308 deriv->data.F32[5] = -2.0*q*py*Y; 1309 deriv->data.F32[6] = -q*X*Y; 1310 1311 return(f); 1312 } 1313 1314 /****************************************************************************** 1315 p_pmMinLM_Gauss2D_Vec(*deriv, *params, *x): this function wraps the above 1316 function in a form that is usable in the LM minimization routines. 1317 *****************************************************************************/ 1318 psVector *p_pmMinLM_Gauss2D_Vec(psImage *deriv, 1319 psVector *params, 1320 psArray *x) 1321 { 1322 PS_IMAGE_CHECK_NULL(deriv, NULL); 1323 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1324 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1325 PS_VECTOR_CHECK_NULL(params, NULL); 1326 PS_VECTOR_CHECK_EMPTY(params, NULL); 1327 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1328 PS_PTR_CHECK_NULL(x, NULL); 1329 if (deriv->numRows != x->n) { 1330 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1331 } 1332 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1333 // XXX: use static memory here. 1334 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1335 1336 for (psS32 i = 0 ; i < x->n ; i++) { 1337 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1338 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1339 } 1340 1341 tmpVec->data.F32[i] = pmMinLM_Gauss2D(tmpRow, 1342 params, 1343 (psVector *) x->data[i]); 1344 1345 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1346 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1347 } 1348 } 1349 1350 psFree(tmpRow); 1351 return(tmpVec); 1352 } 1353 1354 1355 /****************************************************************************** 1356 params->data.F32[0] = So; 1357 params->data.F32[1] = Zo; 1358 params->data.F32[2] = Xo; 1359 params->data.F32[3] = Yo; 1360 params->data.F32[4] = sqrt(2) / SigmaX; 1361 params->data.F32[5] = sqrt(2) / SigmaY; 1362 params->data.F32[6] = Sxy; 1363 *****************************************************************************/ 1364 psF32 pmMinLM_PsuedoGauss2D(psVector *deriv, 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, 1365 1503 psVector *params, 1366 1504 psVector *x) 1367 { 1368 PS_VECTOR_CHECK_NULL(deriv, NAN); 1369 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1370 PS_VECTOR_CHECK_SIZE(deriv, 7, NAN); 1371 PS_VECTOR_CHECK_NULL(params, NAN); 1372 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1373 PS_VECTOR_CHECK_SIZE(params, 7, NAN); 1374 PS_VECTOR_CHECK_NULL(x, NAN); 1375 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1376 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1377 1378 psF32 X = x->data.F32[0] - params->data.F32[2]; 1379 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1380 psF32 px = params->data.F32[4]*X; 1381 psF32 py = params->data.F32[5]*Y; 1382 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1383 psF32 t = 1 + z + 0.5*z*z; 1384 psF32 r = 1.0 / (t*(1 + z/3)); /* exp (-Z) */ 1385 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1386 1387 // 1388 // note difference from a pure gaussian: q = params->data.F32[1]*r 1389 // 1390 1391 psF32 q = params->data.F32[1]*r*r*t; 1392 deriv->data.F32[0] = +1.0; 1393 deriv->data.F32[1] = +r; 1394 deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); 1395 deriv->data.F32[3] = q * 1396 (2.0*py*params->data.F32[5] + params->data.F32[6]*X); 1397 deriv->data.F32[4] = -2.0*q*px*X; 1398 deriv->data.F32[5] = -2.0*q*py*Y; 1399 deriv->data.F32[6] = -q*X*Y; 1400 1401 return(f); 1402 } 1403 1404 /****************************************************************************** 1405 p_pmMinLM_PsuedoGauss2D_Vec(*deriv, *params, *x): this function wraps the 1406 above function in a form that is usable in the LM minimization routines. 1407 *****************************************************************************/ 1408 psVector *p_pmMinLM_PsuedoGauss2D_Vec(psImage *deriv, 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, 1409 1541 psVector *params, 1410 1542 psArray *x) 1411 { 1412 PS_IMAGE_CHECK_NULL(deriv, NULL); 1413 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1414 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1415 PS_VECTOR_CHECK_NULL(params, NULL); 1416 PS_VECTOR_CHECK_EMPTY(params, NULL); 1417 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1418 PS_PTR_CHECK_NULL(x, NULL); 1419 if (deriv->numRows != x->n) { 1420 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1421 } 1422 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1423 // XXX: use static memory here. 1424 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1425 1426 for (psS32 i = 0 ; i < x->n ; i++) { 1427 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1428 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1429 } 1430 1431 tmpVec->data.F32[i] = pmMinLM_PsuedoGauss2D(tmpRow, 1432 params, 1433 (psVector *) x->data[i]); 1434 1435 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1436 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1437 } 1438 } 1439 1440 psFree(tmpRow); 1441 return(tmpVec); 1442 } 1443 1444 1445 1446 1447 /****************************************************************************** 1448 params->data.F32[0] = So; 1449 params->data.F32[1] = Zo; 1450 params->data.F32[2] = Xo; 1451 params->data.F32[3] = Yo; 1452 params->data.F32[4] = Sx; 1453 params->data.F32[5] = Sy; 1454 params->data.F32[6] = Sxy; 1455 params->data.F32[7] = B2; 1456 params->data.F32[8] = B3; 1457 *****************************************************************************/ 1458 psF32 pmMinLM_Wauss2D(psVector *deriv, 1459 psVector *params, 1460 psVector *x) 1461 { 1462 PS_VECTOR_CHECK_NULL(deriv, NAN); 1463 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1464 PS_VECTOR_CHECK_SIZE(deriv, 9, NAN); 1465 PS_VECTOR_CHECK_NULL(params, NAN); 1466 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1467 PS_VECTOR_CHECK_SIZE(params, 9, NAN); 1468 PS_VECTOR_CHECK_NULL(x, NAN); 1469 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1470 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1471 1472 psF32 X = x->data.F32[0] - params->data.F32[2]; 1473 psF32 Y = x->data.F32[1] - params->data.F32[2]; 1474 psF32 px = params->data.F32[4]*X; 1475 psF32 py = params->data.F32[5]*Y; 1476 psF32 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + params->data.F32[6]*X*Y; 1477 psF32 t = 0.5*z*z*(1.0 + params->data.F32[8]*z/3.0); 1478 psF32 r = 1.0 / (1.0 + z + params->data.F32[7]*t); /* exp (-Z) */ 1479 psF32 f = params->data.F32[1]*r + params->data.F32[0]; 1480 1481 // 1482 // note difference from gaussian: q = params->data.F32[1]*r 1483 // 1484 psF32 q = params->data.F32[1]*r*r*(1.0 + params->data.F32[7]*z*(1.0 + params->data.F32[8]*z/2.0)); 1485 deriv->data.F32[0] = +1.0; 1486 deriv->data.F32[1] = +r; 1487 deriv->data.F32[2] = q*(2.0*px*params->data.F32[4] + params->data.F32[6]*Y); 1488 deriv->data.F32[3] = q*(2.0*py*params->data.F32[5] + params->data.F32[6]*X); 1489 deriv->data.F32[4] = -2.0*q*px*X; 1490 deriv->data.F32[5] = -2.0*q*py*Y; 1491 deriv->data.F32[6] = -q*X*Y; 1492 deriv->data.F32[7] = -100.0*params->data.F32[1]*r*r*t; 1493 deriv->data.F32[8] = -100.0*params->data.F32[1]*r*r*params->data.F32[7]*(z*z*z)/6.0; 1494 // 1495 // The values of 100 dampen the swing of params->data.F32[7,8] */ 1496 // 1497 1498 return(f); 1499 } 1500 1501 /****************************************************************************** 1502 p_pmMinLM_Wauss2D_Vec(*deriv, *params, *x): this function wraps the above 1503 function in a form that is usable in the LM minimization routines. 1504 *****************************************************************************/ 1505 psVector *p_pmMinLM_Wauss2D_Vec(psImage *deriv, 1506 psVector *params, 1507 psArray *x) 1508 { 1509 PS_IMAGE_CHECK_NULL(deriv, NULL); 1510 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1511 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1512 PS_VECTOR_CHECK_NULL(params, NULL); 1513 PS_VECTOR_CHECK_EMPTY(params, NULL); 1514 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1515 PS_PTR_CHECK_NULL(x, NULL); 1516 if (deriv->numRows != x->n) { 1517 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1518 } 1519 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1520 // XXX: use static memory here. 1521 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1522 1523 for (psS32 i = 0 ; i < x->n ; i++) { 1524 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1525 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1526 } 1527 1528 tmpVec->data.F32[i] = pmMinLM_Wauss2D(tmpRow, 1529 params, 1530 (psVector *) x->data[i]); 1531 1532 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1533 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1534 } 1535 } 1536 1537 psFree(tmpRow); 1538 return(tmpVec); 1539 } 1540 1541 1542 1543 1544 1545 1546 // XXX: What should these be? 1547 #define FFACTOR 1.0 1548 #define FSCALE 1.0 1549 /****************************************************************************** 1550 params->data.F32[0] = So; 1551 params->data.F32[1] = Zo; 1552 params->data.F32[2] = Xo; 1553 params->data.F32[3] = Yo; 1554 params->data.F32[4] = SxInner; 1555 params->data.F32[5] = SyInner; 1556 params->data.F32[6] = SxyInner; 1557 params->data.F32[7] = SxOuter; 1558 params->data.F32[8] = SyOuter; 1559 params->data.F32[9] = SxyOuter; 1560 params->data.F32[10] = N; 1561 *****************************************************************************/ 1562 psF32 pmMinLM_TwistGauss2D(psVector *deriv, 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, 1563 1895 psVector *params, 1564 1896 psVector *x) 1565 { 1566 PS_VECTOR_CHECK_NULL(deriv, NAN); 1567 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1568 PS_VECTOR_CHECK_SIZE(deriv, 11, NAN); 1569 PS_VECTOR_CHECK_NULL(params, NAN); 1570 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1571 PS_VECTOR_CHECK_SIZE(params, 11, NAN); 1572 PS_VECTOR_CHECK_NULL(x, NAN); 1573 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1574 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1575 1576 psF32 X = x->data.F32[0] - params->data.F32[2]; 1577 psF32 Y = x->data.F32[1] - params->data.F32[3]; 1578 psF32 px1 = params->data.F32[4]*X; 1579 psF32 py1 = params->data.F32[5]*Y; 1580 psF32 px2 = params->data.F32[7]*X; 1581 psF32 py2 = params->data.F32[8]*Y; 1582 psF32 z1 = 0.5*PS_SQR(px1) + 0.5*PS_SQR(py1) + params->data.F32[4]*X*Y; 1583 psF32 z2 = 0.5*PS_SQR(px2) + 0.5*PS_SQR(py2) + params->data.F32[9]*X*Y; 1584 psF32 r = 1.0 / (1.0 + z1 + pow(z2,params->data.F32[10])); 1585 1586 1587 psF32 f = params->data.F32[5]*r + params->data.F32[6]; 1588 psF32 q1 = params->data.F32[5]*PS_SQR(r); 1589 psF32 q2 = params->data.F32[5]*PS_SQR(r)*params->data.F32[10]*pow(z2,(params->data.F32[10]-1.0)); 1590 deriv->data.F32[0] = +1.0; 1591 deriv->data.F32[1] = +r; 1592 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); 1593 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); 1594 1595 // 1596 // These fudge factors impede the growth of params->data.F32[4] beyond 1597 // params->data.F32[7]. 1598 // 1599 psF32 f1 = fabs(params->data.F32[7]) / fabs(params->data.F32[4]); 1600 psF32 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; 1601 deriv->data.F32[4] = -2.0*q1*px1*X*f2; 1602 1603 // 1604 // These fudge factors impede the growth of params->data.F32[5] beyond 1605 // params->data.F32[8]. 1606 // 1607 f1 = fabs(params->data.F32[8]) / fabs(params->data.F32[5]); 1608 f2 = (f1 < FSCALE) ? 1.0 : FFACTOR*(f1 - FSCALE) + 1.0; 1609 deriv->data.F32[5] = -2.0*q1*py1*Y*f2; 1610 deriv->data.F32[6] = -q1*X*Y; 1611 deriv->data.F32[7] = -2.0*q2*px2*X; 1612 deriv->data.F32[8] = -2.0*q2*py2*Y; 1613 deriv->data.F32[9] = -q2*X*Y; 1614 deriv->data.F32[10] = -q1*log(z2); 1615 1616 return(f); 1617 } 1618 1619 /****************************************************************************** 1620 p_pmMinLM_TwistGauss2D_Vec(*deriv, *params, *x): this function wraps the above 1621 function in a form that is usable in the LM minimization routines. 1622 *****************************************************************************/ 1623 psVector *p_pmMinLM_TwistGauss2D_Vec(psImage *deriv, 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, 1624 1916 psVector *params, 1625 1917 psArray *x) 1626 { 1627 PS_IMAGE_CHECK_NULL(deriv, NULL); 1628 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1629 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1630 PS_VECTOR_CHECK_NULL(params, NULL); 1631 PS_VECTOR_CHECK_EMPTY(params, NULL); 1632 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1633 PS_PTR_CHECK_NULL(x, NULL); 1634 if (deriv->numRows != x->n) { 1635 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1636 } 1637 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1638 // XXX: use static memory here. 1639 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1640 1641 for (psS32 i = 0 ; i < x->n ; i++) { 1642 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1643 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1644 } 1645 1646 tmpVec->data.F32[i] = pmMinLM_TwistGauss2D(tmpRow, 1647 params, 1648 (psVector *) x->data[i]); 1649 1650 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1651 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1652 } 1653 } 1654 1655 psFree(tmpRow); 1656 return(tmpVec); 1657 } 1658 1659 1660 1661 /****************************************************************************** 1662 float Sersic() 1663 params->data.F32[0] = So; 1664 params->data.F32[1] = Zo; 1665 params->data.F32[2] = Xo; 1666 params->data.F32[3] = Yo; 1667 params->data.F32[4] = Sx; 1668 params->data.F32[5] = Sy; 1669 params->data.F32[6] = Sxy; 1670 params->data.F32[7] = Nexp; 1671 *****************************************************************************/ 1672 psF32 pmMinLM_Sersic(psVector *deriv, 1673 psVector *params, 1674 psVector *x) 1675 { 1676 PS_VECTOR_CHECK_NULL(deriv, NAN); 1677 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1678 PS_VECTOR_CHECK_SIZE(deriv, 8, NAN); 1679 PS_VECTOR_CHECK_NULL(params, NAN); 1680 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1681 PS_VECTOR_CHECK_SIZE(params, 8, NAN); 1682 PS_VECTOR_CHECK_NULL(x, NAN); 1683 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1684 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1685 1686 psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); 1687 return(0.0); 1688 } 1689 /****************************************************************************** 1690 p_pmMinLM_Sersic_Vec(*deriv, *params, *x): this function wraps the above 1691 function in a form that is usable in the LM minimization routines. 1692 *****************************************************************************/ 1693 psVector *p_pmMinLM_Sersic_Vec(psImage *deriv, 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, 1694 1967 psVector *params, 1695 psArray *x) 1696 { 1697 PS_IMAGE_CHECK_NULL(deriv, NULL); 1698 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1699 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1700 PS_VECTOR_CHECK_NULL(params, NULL); 1701 PS_VECTOR_CHECK_EMPTY(params, NULL); 1702 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1703 PS_PTR_CHECK_NULL(x, NULL); 1704 if (deriv->numRows != x->n) { 1705 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1706 } 1707 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1708 // XXX: use static memory here. 1709 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1710 1711 for (psS32 i = 0 ; i < x->n ; i++) { 1712 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1713 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1714 } 1715 1716 tmpVec->data.F32[i] = pmMinLM_Sersic(tmpRow, 1717 params, 1718 (psVector *) x->data[i]); 1719 1720 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1721 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1722 } 1723 } 1724 1725 psFree(tmpRow); 1726 return(tmpVec); 1727 } 1728 1729 /****************************************************************************** 1730 float SersicBulge() 1731 params->data.F32[0] So; 1732 params->data.F32[1] Zo; 1733 params->data.F32[2] Xo; 1734 params->data.F32[3] Yo; 1735 params->data.F32[4] SxInner; 1736 params->data.F32[5] SyInner; 1737 params->data.F32[6] SxyInner; 1738 params->data.F32[7] Zd; 1739 params->data.F32[8] SxOuter; 1740 params->data.F32[9] SyOuter; 1741 params->data.F32[10] = SxyOuter; 1742 params->data.F32[11] = Nexp; 1743 *****************************************************************************/ 1744 psF32 pmMinLM_SersicCore(psVector *deriv, 1745 psVector *params, 1746 psVector *x) 1747 { 1748 PS_VECTOR_CHECK_NULL(deriv, NAN); 1749 PS_VECTOR_CHECK_TYPE(deriv, PS_TYPE_F32, NAN); 1750 PS_VECTOR_CHECK_SIZE(deriv, 12, NAN); 1751 PS_VECTOR_CHECK_NULL(params, NAN); 1752 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NAN); 1753 PS_VECTOR_CHECK_SIZE(params, 12, NAN); 1754 PS_VECTOR_CHECK_NULL(x, NAN); 1755 PS_VECTOR_CHECK_TYPE(x, PS_TYPE_F32, NAN); 1756 PS_VECTOR_CHECK_SIZE(x, 2, NAN); 1757 1758 psError(PS_ERR_UNKNOWN, true, "This function is not implemented yet."); 1759 return(0.0); 1760 } 1761 /****************************************************************************** 1762 p_pmMinLM_SersicCore_Vec(*deriv, *params, *x): this function wraps the above 1763 function in a form that is usable in the LM minimization routines. 1764 *****************************************************************************/ 1765 psVector *p_pmMinLM_SersicCore_Vec(psImage *deriv, 1766 psVector *params, 1767 psArray *x) 1768 { 1769 PS_IMAGE_CHECK_NULL(deriv, NULL); 1770 PS_IMAGE_CHECK_EMPTY(deriv, NULL); 1771 PS_IMAGE_CHECK_TYPE(deriv, PS_TYPE_F32, NULL); 1772 PS_VECTOR_CHECK_NULL(params, NULL); 1773 PS_VECTOR_CHECK_EMPTY(params, NULL); 1774 PS_VECTOR_CHECK_TYPE(params, PS_TYPE_F32, NULL); 1775 PS_PTR_CHECK_NULL(x, NULL); 1776 if (deriv->numRows != x->n) { 1777 psError(PS_ERR_UNKNOWN, true, "deriv must have one row for each coordinate set in x."); 1778 } 1779 psVector *tmpVec = psVectorAlloc(x->n, PS_TYPE_F32); 1780 // XXX: use static memory here. 1781 psVector *tmpRow = psVectorAlloc(deriv->numCols, PS_TYPE_F32); 1782 1783 for (psS32 i = 0 ; i < x->n ; i++) { 1784 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1785 tmpRow->data.F32[j] = deriv->data.F32[i][j]; 1786 } 1787 1788 tmpVec->data.F32[i] = pmMinLM_SersicCore(tmpRow, 1789 params, 1790 (psVector *) x->data[i]); 1791 1792 for (psS32 j = 0 ; j < tmpRow->n ; j++) { 1793 deriv->data.F32[i][j] = tmpRow->data.F32[j]; 1794 } 1795 } 1796 1797 psFree(tmpRow); 1798 return(tmpVec); 1799 } 1800 1801 1802 1803 /****************************************************************************** 1804 *****************************************************************************/ 1805 psF32 pmMinLM_PsuedoSersic(psVector *deriv, 1806 psVector *params, 1807 psVector *x) 1808 { 1809 return(0.0); 1810 } 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 } -
trunk/psModules/src/pmObjects.h
r3498 r3510 5 5 * @author GLG, MHPCC 6 6 * 7 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $8 * @date $Date: 2005-03-2 4 22:36:52$7 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2005-03-25 23:18:35 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 65 65 } psModelType; 66 66 67 // XXX: It will be better if params, and dparams, were psVectors. 67 68 typedef struct 68 69 { … … 173 174 /****************************************************************************** 174 175 *****************************************************************************/ 176 typedef enum { 177 PS_CONTOUR_CRUDE, 178 } pmContourType; 179 175 180 psArray *pmSourceContour(psSource *source, 176 181 const psImage *image, 177 182 psF32 level, 178 p sS32mode);183 pmContourType mode); 179 184 180 185 /******************************************************************************
Note:
See TracChangeset
for help on using the changeset viewer.
