- Timestamp:
- Oct 2, 2009, 12:11:34 PM (17 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branches/20090715/psModules/src/objects/models/pmModel_RGAUSS.c
r25683 r25745 20 20 *****************************************************************************/ 21 21 22 #include <stdio.h> 23 #include <pslib.h> 24 25 #include "pmMoments.h" 26 #include "pmPeaks.h" 27 #include "pmSource.h" 28 #include "pmModel.h" 29 #include "pmModel_RGAUSS.h" 30 22 31 # define PM_MODEL_FUNC pmModelFunc_RGAUSS 23 32 # define PM_MODEL_FLUX pmModelFlux_RGAUSS … … 28 37 # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_RGAUSS 29 38 # define PM_MODEL_FIT_STATUS pmModelFitStatus_RGAUSS 39 # define PM_MODEL_SET_LIMITS pmModelSetLimits_RGAUSS 30 40 31 41 // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) 32 42 // 0.5 PIX: the parameters are defined in terms of pixel coords, so the incoming pixcoords 33 43 // values need to be pixel coords 44 45 // Lax parameter limits 46 static float paramsMinLax[] = { -1.0e3, 1.0e-2, -100, -100, 0.5, 0.5, -1.0, 1.25 }; 47 static float paramsMaxLax[] = { 1.0e5, 1.0e8, 1.0e4, 1.0e4, 100, 100, 1.0, 4.0 }; 48 49 // Strict parameter limits 50 static float *paramsMinStrict = paramsMinLax; 51 static float *paramsMaxStrict = paramsMaxLax; 52 53 // Parameter limits to use 54 static float *paramsMinUse = paramsMinLax; 55 static float *paramsMaxUse = paramsMaxLax; 56 static float betaUse[] = { 1000, 3e6, 5, 5, 0.5, 0.5, 0.5, 0.5 }; 57 58 static bool limitsApply = true; // Apply limits? 59 34 60 psF32 PM_MODEL_FUNC (psVector *deriv, 35 61 const psVector *params, … … 76 102 # define AR_MAX 20.0 77 103 # define AR_RATIO 0.99 104 78 105 bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) 79 106 { 80 float beta_lim = 0, params_min = 0, params_max = 0; 81 float f1 = 0, f2 = 0, q1 = 0, q2 = 0; 107 if (!limitsApply) { 108 return true; 109 } 110 psAssert(nParam >= 0 && nParam <= PM_PAR_7, "Parameter index is out of bounds"); 82 111 83 112 // we need to calculate the limits for SXY specially 113 float q2 = NAN; 84 114 if (nParam == PM_PAR_SXY) { 85 f 1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]);86 f 2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]);87 q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2);115 float f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]); 116 float f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]); 117 float q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); 88 118 q1 = (q1 < 0.0) ? 0.0 : q1; 89 119 // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg.. Saturate at that 90 120 // angle and let f2,f1 fight it out 91 q2 = 0.5*sqrt(q1);121 q2 = 0.5*sqrtf(q1); 92 122 } 93 123 94 124 switch (mode) { 95 case PS_MINIMIZE_BETA_LIMIT: 96 switch (nParam) { 97 case PM_PAR_SKY: 98 beta_lim = 1000; 99 break; 100 case PM_PAR_I0: 101 beta_lim = 3e6; 102 break; 103 case PM_PAR_XPOS: 104 beta_lim = 5; 105 break; 106 case PM_PAR_YPOS: 107 beta_lim = 5; 108 break; 109 case PM_PAR_SXX: 110 beta_lim = 0.5; 111 break; 112 case PM_PAR_SYY: 113 beta_lim = 0.5; 114 break; 115 case PM_PAR_SXY: 116 beta_lim = 0.5*q2; 117 break; 118 case PM_PAR_7: 119 beta_lim = 0.5; 120 break; 121 default: 122 psAbort("invalid parameter %d for beta test", nParam); 123 } 124 if (fabs(beta[nParam]) > fabs(beta_lim)) { 125 beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim); 126 psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g", 127 nParam, beta[nParam], beta_lim); 128 return false; 129 } 130 return true; 131 case PS_MINIMIZE_PARAM_MIN: 132 switch (nParam) { 133 case PM_PAR_SKY: 134 params_min = -1000; 135 break; 136 case PM_PAR_I0: 137 params_min = 0.01; 138 break; 139 case PM_PAR_XPOS: 140 params_min = -100; 141 break; 142 case PM_PAR_YPOS: 143 params_min = -100; 144 break; 145 case PM_PAR_SXX: 146 params_min = 0.5; 147 break; 148 case PM_PAR_SYY: 149 params_min = 0.5; 150 break; 151 case PM_PAR_SXY: 152 params_min = -q2; 153 break; 154 case PM_PAR_7: 155 params_min = 1.25; 156 break; 157 default: 158 psAbort("invalid parameter %d for param min test", nParam); 159 } 160 if (params[nParam] < params_min) { 161 params[nParam] = params_min; 162 psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g", 163 nParam, params[nParam], params_min); 164 return false; 165 } 166 return true; 167 case PS_MINIMIZE_PARAM_MAX: 168 switch (nParam) { 169 case PM_PAR_SKY: 170 params_max = 1e5; 171 break; 172 case PM_PAR_I0: 173 params_max = 1e8; 174 break; 175 case PM_PAR_XPOS: 176 params_max = 1e4; 177 break; 178 case PM_PAR_YPOS: 179 params_max = 1e4; 180 break; 181 case PM_PAR_SXX: 182 params_max = 100; 183 break; 184 case PM_PAR_SYY: 185 params_max = 100; 186 break; 187 case PM_PAR_SXY: 188 params_max = +q2; 189 break; 190 case PM_PAR_7: 191 params_max = 4.0; 192 break; 193 default: 194 psAbort("invalid parameter %d for param max test", nParam); 195 } 196 if (params[nParam] > params_max) { 197 params[nParam] = params_max; 198 psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g", 199 nParam, params[nParam], params_max); 200 return false; 201 } 202 return true; 203 default: 125 case PS_MINIMIZE_BETA_LIMIT: { 126 psAssert(beta, "Require beta to limit beta"); 127 float limit = betaUse[nParam]; 128 if (nParam == PM_PAR_SXY) { 129 limit *= q2; 130 } 131 if (fabs(beta[nParam]) > fabs(limit)) { 132 beta[nParam] = (beta[nParam] > 0) ? fabs(limit) : -fabs(limit); 133 psTrace("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g", 134 nParam, beta[nParam], limit); 135 return false; 136 } 137 return true; 138 } 139 case PS_MINIMIZE_PARAM_MIN: { 140 psAssert(params, "Require parameters to limit parameters"); 141 psAssert(paramsMinUse, "Require parameter limits to limit parameters"); 142 float limit = paramsMinUse[nParam]; 143 if (nParam == PM_PAR_SXY) { 144 limit *= q2; 145 } 146 if (params[nParam] < limit) { 147 params[nParam] = limit; 148 psTrace("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g", 149 nParam, params[nParam], limit); 150 return false; 151 } 152 return true; 153 } 154 case PS_MINIMIZE_PARAM_MAX: { 155 psAssert(params, "Require parameters to limit parameters"); 156 psAssert(paramsMaxUse, "Require parameter limits to limit parameters"); 157 float limit = paramsMaxUse[nParam]; 158 if (nParam == PM_PAR_SXY) { 159 limit *= q2; 160 } 161 if (params[nParam] > limit) { 162 params[nParam] = limit; 163 psTrace("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g", 164 nParam, params[nParam], limit); 165 return false; 166 } 167 return true; 168 } 169 default: 204 170 psAbort("invalid choice for limits"); 205 171 } … … 207 173 return false; 208 174 } 175 209 176 210 177 // make an initial guess for parameters … … 454 421 } 455 422 423 424 void PM_MODEL_SET_LIMITS(pmModelLimitsType type) 425 { 426 switch (type) { 427 case PM_MODEL_LIMITS_NONE: 428 paramsMinUse = NULL; 429 paramsMaxUse = NULL; 430 limitsApply = true; 431 break; 432 case PM_MODEL_LIMITS_IGNORE: 433 paramsMinUse = NULL; 434 paramsMaxUse = NULL; 435 limitsApply = false; 436 break; 437 case PM_MODEL_LIMITS_LAX: 438 paramsMinUse = paramsMinLax; 439 paramsMaxUse = paramsMaxLax; 440 limitsApply = true; 441 break; 442 case PM_MODEL_LIMITS_STRICT: 443 paramsMinUse = paramsMinStrict; 444 paramsMaxUse = paramsMaxStrict; 445 limitsApply = true; 446 break; 447 default: 448 psAbort("Unrecognised model limits type: %x", type); 449 } 450 return; 451 } 452 456 453 # undef PM_MODEL_FUNC 457 454 # undef PM_MODEL_FLUX … … 462 459 # undef PM_MODEL_PARAMS_FROM_PSF 463 460 # undef PM_MODEL_FIT_STATUS 461 # undef PM_MODEL_SET_LIMITS
Note:
See TracChangeset
for help on using the changeset viewer.
