Changeset 14323
- Timestamp:
- Jul 19, 2007, 2:27:31 PM (19 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/objects/models/pmModel_RGAUSS.c
r11687 r14323 28 28 # define PM_MODEL_FIT_STATUS pmModelFitStatus_RGAUSS 29 29 30 psF 64 PS_MODEL_FUNC (psVector *deriv,30 psF32 PM_MODEL_FUNC (psVector *deriv, 31 31 const psVector *params, 32 32 const psVector *pixcoord) … … 39 39 psF32 py = Y / PAR[PM_PAR_SYY]; 40 40 psF32 z = PS_SQR(px) + PS_SQR(py) + X*Y*PAR[PM_PAR_SXY]; 41 42 assert (z >= 0); 41 43 42 44 psF32 p = pow(z, PAR[PM_PAR_7] - 1.0); … … 58 60 dPAR[PM_PAR_SYY] = +2.0*q*py*py/PAR[PM_PAR_SYY]; 59 61 dPAR[PM_PAR_SXY] = -q*X*Y; 60 dPAR[PM_PAR_7] = -5.0*t*log(z)*p*z; 62 63 // this model derivative is undefined at z = 0.0, but is actually 0.0 64 dPAR[PM_PAR_7] = (z == 0.0) ? 0.0 : -5.0*t*log(z)*p*z; 61 65 } 62 66 return(f); 63 67 } 64 68 65 bool PM_MODEL_LIMITS (psVector **beta_lim, psVector **params_min, psVector **params_max) 66 { 67 68 *beta_lim = psVectorAlloc (8, PS_TYPE_F32); 69 *params_min = psVectorAlloc (8, PS_TYPE_F32); 70 *params_max = psVectorAlloc (8, PS_TYPE_F32); 71 72 beta_lim[0][0].data.F32[PM_PAR_SKY] = 1000; 73 beta_lim[0][0].data.F32[PM_PAR_I0] = 3e6; 74 beta_lim[0][0].data.F32[PM_PAR_XPOS] = 5; 75 beta_lim[0][0].data.F32[PM_PAR_YPOS] = 5; 76 beta_lim[0][0].data.F32[PM_PAR_SXX] = 0.5; 77 beta_lim[0][0].data.F32[PM_PAR_SYY] = 0.5; 78 beta_lim[0][0].data.F32[PM_PAR_SXY] = 0.5; 79 beta_lim[0][0].data.F32[PM_PAR_7] = 0.5; 80 81 params_min[0][0].data.F32[PM_PAR_SKY] = -1000; 82 params_min[0][0].data.F32[PM_PAR_I0] = 0; 83 params_min[0][0].data.F32[PM_PAR_XPOS] = -100; 84 params_min[0][0].data.F32[PM_PAR_YPOS] = -100; 85 params_min[0][0].data.F32[PM_PAR_SXX] = 0.5; 86 params_min[0][0].data.F32[PM_PAR_SYY] = 0.5; 87 params_min[0][0].data.F32[PM_PAR_SXY] = -5.0; 88 params_min[0][0].data.F32[PM_PAR_7] = 1.25; 89 90 params_max[0][0].data.F32[PM_PAR_SKY] = 1e5; 91 params_max[0][0].data.F32[PM_PAR_I0] = 1e8; 92 params_max[0][0].data.F32[PM_PAR_XPOS] = 1e4; // this should be set by image dimensions! 93 params_max[0][0].data.F32[PM_PAR_YPOS] = 1e4; // this should be set by image dimensions! 94 params_max[0][0].data.F32[PM_PAR_SXX] = 100.0; 95 params_max[0][0].data.F32[PM_PAR_SYY] = 100.0; 96 params_max[0][0].data.F32[PM_PAR_SXY] = +5.0; 97 params_max[0][0].data.F32[PM_PAR_7] = 4.0; 98 99 return (TRUE); 100 } 101 102 bool PS_MODEL_GUESS (psModel *model, psSource *source) 69 // define the parameter limits 70 // AR_MAX is the maximum allowed axis ratio 71 // AR_RATIO is ((1-R)/(1+R))^2 where R = AR_MAX^(-2) 72 # define AR_MAX 20.0 73 # define AR_RATIO 0.99 74 bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) 75 { 76 float beta_lim = 0, params_min = 0, params_max = 0; 77 float f1 = 0, f2 = 0, q1 = 0, q2 = 0; 78 79 // we need to calculate the limits for SXY specially 80 if (nParam == PM_PAR_SXY) { 81 f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]); 82 f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]); 83 q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); 84 q1 = (q1 < 0.0) ? 0.0 : q1; 85 // if q1 < 0.0, f2 ~ f1, we have a very large axis ratio near 45deg.. Saturate at that 86 // angle and let f2,f1 fight it out 87 q2 = 0.5*sqrt (q1); 88 } 89 90 switch (mode) { 91 case PS_MINIMIZE_BETA_LIMIT: 92 switch (nParam) { 93 case PM_PAR_SKY: 94 beta_lim = 1000; 95 break; 96 case PM_PAR_I0: 97 beta_lim = 3e6; 98 break; 99 case PM_PAR_XPOS: 100 beta_lim = 5; 101 break; 102 case PM_PAR_YPOS: 103 beta_lim = 5; 104 break; 105 case PM_PAR_SXX: 106 beta_lim = 0.5; 107 break; 108 case PM_PAR_SYY: 109 beta_lim = 0.5; 110 break; 111 case PM_PAR_SXY: 112 beta_lim = 0.5*q2; 113 break; 114 case PM_PAR_7: 115 beta_lim = 0.5; 116 break; 117 default: 118 psAbort("invalid parameter %d for beta test", nParam); 119 } 120 if (fabs(beta[nParam]) > fabs(beta_lim)) { 121 beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim); 122 psTrace ("psModules.objects", 5, "|beta[nParam==%d]| > |beta_lim|; %g v. %g", 123 nParam, beta[nParam], beta_lim); 124 return false; 125 } 126 return true; 127 case PS_MINIMIZE_PARAM_MIN: 128 switch (nParam) { 129 case PM_PAR_SKY: 130 params_min = -1000; 131 break; 132 case PM_PAR_I0: 133 params_min = 0; 134 break; 135 case PM_PAR_XPOS: 136 params_min = -100; 137 break; 138 case PM_PAR_YPOS: 139 params_min = -100; 140 break; 141 case PM_PAR_SXX: 142 params_min = 0.5; 143 break; 144 case PM_PAR_SYY: 145 params_min = 0.5; 146 break; 147 case PM_PAR_SXY: 148 params_min = -q2; 149 break; 150 case PM_PAR_7: 151 params_min = 1.25; 152 break; 153 default: 154 psAbort("invalid parameter %d for param min test", nParam); 155 } 156 if (params[nParam] < params_min) { 157 params[nParam] = params_min; 158 psTrace ("psModules.objects", 5, "params[nParam==%d] < params_min; %g v. %g", 159 nParam, params[nParam], params_min); 160 return false; 161 } 162 return true; 163 case PS_MINIMIZE_PARAM_MAX: 164 switch (nParam) { 165 case PM_PAR_SKY: 166 params_max = 1e5; 167 break; 168 case PM_PAR_I0: 169 params_max = 1e8; 170 break; 171 case PM_PAR_XPOS: 172 params_max = 1e4; 173 break; 174 case PM_PAR_YPOS: 175 params_max = 1e4; 176 break; 177 case PM_PAR_SXX: 178 params_max = 100; 179 break; 180 case PM_PAR_SYY: 181 params_max = 100; 182 break; 183 case PM_PAR_SXY: 184 params_max = +q2; 185 break; 186 case PM_PAR_7: 187 params_max = 4.0; 188 break; 189 default: 190 psAbort("invalid parameter %d for param max test", nParam); 191 } 192 if (params[nParam] > params_max) { 193 params[nParam] = params_max; 194 psTrace ("psModules.objects", 5, "params[nParam==%d] > params_max; %g v. %g", 195 nParam, params[nParam], params_max); 196 return false; 197 } 198 return true; 199 default: 200 psAbort("invalid choice for limits"); 201 } 202 psAbort("should not reach here"); 203 return false; 204 } 205 206 // make an initial guess for parameters 207 bool PM_MODEL_GUESS (pmModel *model, pmSource *source) 103 208 { 104 209 pmMoments *moments = source->moments; 105 210 pmPeak *peak = source->peak; 106 psF32 *PAR = model->params->data.F32; 107 108 psEllipseAxes axes; 109 psEllipseShape shape; 110 psEllipseMoments tmpMoments; 111 112 // XXX fix this stuff : should be using correct ellipse relationships... 113 tmpMoments.x2 = PS_SQR(source->moments->Sx); 114 tmpMoments.y2 = PS_SQR(source->moments->Sy); 115 tmpMoments.xy = source->moments->Sxy; 116 117 axes = psEllipseMomentsToAxes(tmpMoments); 118 shape = psEllipseAxesToShape(axes); 119 120 PAR[PM_PAR_SKY] = moments->Sky; 121 PAR[PM_PAR_I0] = moments->Peak - moments->Sky; 211 psF32 *PAR = model->params->data.F32; 212 213 psEllipseMoments emoments; 214 emoments.x2 = moments->Sx; 215 emoments.y2 = moments->Sy; 216 emoments.xy = moments->Sxy; 217 218 // force the axis ratio to be < 20.0 219 psEllipseAxes axes = psEllipseMomentsToAxes (emoments, 20.0); 220 221 if (!isfinite(axes.major)) return false; 222 if (!isfinite(axes.minor)) return false; 223 if (!isfinite(axes.theta)) return false; 224 225 psEllipseShape shape = psEllipseAxesToShape (axes); 226 227 if (!isfinite(shape.sx)) return false; 228 if (!isfinite(shape.sy)) return false; 229 if (!isfinite(shape.sxy)) return false; 230 231 PAR[PM_PAR_SKY] = moments->Sky; 232 PAR[PM_PAR_I0] = moments->Peak - moments->Sky; 122 233 PAR[PM_PAR_XPOS] = peak->x; 123 234 PAR[PM_PAR_YPOS] = peak->y; 124 PAR[PM_PAR_SXX] = PS_MAX(0.5, moments->Sx);125 PAR[PM_PAR_SYY] = PS_MAX(0.5, moments->Sy);235 PAR[PM_PAR_SXX] = PS_MAX(0.5, M_SQRT2*shape.sx); 236 PAR[PM_PAR_SYY] = PS_MAX(0.5, M_SQRT2*shape.sy); 126 237 PAR[PM_PAR_SXY] = shape.sxy; 127 PAR[PM_PAR_7] = 2.0; 238 PAR[PM_PAR_7] = 2.25; 239 128 240 return(true); 129 241 } 130 242 131 psF64 P S_MODEL_FLUX (const psVector *params)243 psF64 PM_MODEL_FLUX (const psVector *params) 132 244 { 133 245 float norm, z; … … 136 248 psF32 *PAR = params->data.F32; 137 249 138 shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0);139 shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0);250 shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; 251 shape.sy = PAR[PM_PAR_SYY] / M_SQRT2; 140 252 shape.sxy = PAR[PM_PAR_SXY]; 141 253 142 254 // Area is equivalent to 2 pi sigma^2 143 psEllipseAxes axes = psEllipseShapeToAxes (shape );255 psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); 144 256 psF64 Area = 2.0 * M_PI * axes.major * axes.minor; 145 257 … … 167 279 // define this function so it never returns Inf or NaN 168 280 // return the radius which yields the requested flux 169 psF64 PS_MODEL_RADIUS (const psVector *params, psF64 flux) 170 { 171 psF64 z, f, p; 281 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) 282 { 283 psF64 z, f; 284 int Nstep = 0; 172 285 psEllipseShape shape; 173 286 … … 181 294 return (1.0); 182 295 183 shape.sx = PAR[PM_PAR_SXX] / sqrt(2.0);184 shape.sy = PAR[PM_PAR_SYY] / sqrt(2.0);296 shape.sx = PAR[PM_PAR_SXX] / M_SQRT2; 297 shape.sy = PAR[PM_PAR_SYY] / M_SQRT2; 185 298 shape.sxy = PAR[PM_PAR_SXY]; 186 299 187 psEllipseAxes axes = psEllipseShapeToAxes (shape );300 psEllipseAxes axes = psEllipseShapeToAxes (shape, 20.0); 188 301 psF64 sigma = axes.major; 189 302 190 psF64 dz = 1.0 / (2.0 * sigma*sigma);191 303 psF64 limit = flux / PAR[PM_PAR_I0]; 192 304 … … 196 308 197 309 // choose a z value guaranteed to be beyond our limit 198 float z0 = pow((1.0 / limit), (1.0 / 2.25));199 float z1 = (1.0 / limit) / PAR[PM_PAR_7];310 float z0 = pow((1.0 / limit), (1.0 / PAR[PM_PAR_7])); 311 float z1 = (1.0 / limit); 200 312 z1 = PS_MAX (z0, z1); 201 313 z0 = 0.0; … … 224 336 } 225 337 226 bool P S_MODEL_FROM_PSF (psModel *modelPSF, psModel *modelFLT, pmPSF *psf)338 bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf) 227 339 { 228 340 … … 243 355 } 244 356 245 // the 2D model for SXY actually fits SXY / (SXX^-2 + SYY^-2); correct here 246 out[PM_PAR_SXY] = pmPSF_SXYtoModel (out); 247 248 return(true); 357 // the 2D PSF model fits polarization terms (E0,E1,E2) 358 // convert to shape terms (SXX,SYY,SXY) 359 if (!pmPSF_FitToModel (out, 0.1)) { 360 psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", 361 in[PM_PAR_YPOS], in[PM_PAR_XPOS]); 362 return false; 363 } 364 365 // apply the model limits here: this truncates excessive extrapolation 366 // XXX do we need to do this still? should we put in asserts to test? 367 for (int i = 0; i < psf->params_NEW->n; i++) { 368 // apply the limits to all components or just the psf-model parameters? 369 if (psf->params_NEW->data[i] == NULL) 370 continue; 371 372 bool status = true; 373 status &= PM_MODEL_LIMITS(PS_MINIMIZE_PARAM_MIN, i, out, NULL); 374 status &= PM_MODEL_LIMITS(PS_MINIMIZE_PARAM_MAX, i, out, NULL); 375 if (!status) { 376 psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", 377 in[PM_PAR_XPOS], in[PM_PAR_YPOS]); 378 modelPSF->flags |= PM_MODEL_STATUS_LIMITS; 379 } 380 } 381 382 return true; 249 383 } 250 384
Note:
See TracChangeset
for help on using the changeset viewer.
