Changeset 15843
- Timestamp:
- Dec 14, 2007, 3:24:38 PM (18 years ago)
- Location:
- trunk
- Files:
-
- 5 edited
-
psModules/src/objects/pmModel.h (modified) (3 diffs)
-
psModules/src/objects/pmModelUtils.c (modified) (3 diffs)
-
psModules/src/objects/pmModelUtils.h (modified) (2 diffs)
-
psphot/src/models/pmModel_STRAIL.c (modified) (6 diffs)
-
psphot/src/models/pmModel_TEST1.c (modified) (10 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psModules/src/objects/pmModel.h
r15562 r15843 5 5 * @author EAM, IfA 6 6 * 7 * @version $Revision: 1.1 3$ $Name: not supported by cvs2svn $8 * @date $Date: 2007-1 1-10 01:09:20$7 * @version $Revision: 1.14 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2007-12-15 01:21:33 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 63 63 // This function constructs the PSF model for the given source based on the 64 64 // supplied psf and the EXT model for the object. 65 typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, pmPSF *psf);65 typedef bool (*pmModelFromPSFFunc)(pmModel *modelPSF, pmModel *modelEXT, const pmPSF *psf); 66 66 67 67 // This function sets the model parameters based on the PSF for a given coordinate and central 68 68 // intensity 69 typedef bool (*pmModelParamsFromPSF)(pmModel *model, pmPSF *psf, float Xo, float Yo, float Io);69 typedef bool (*pmModelParamsFromPSF)(pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io); 70 70 71 71 // This function returns the success / failure status of the given model fit … … 166 166 167 167 bool pmModelAddWithOffset(psImage *image, 168 169 170 171 172 173 168 psImage *mask, 169 pmModel *model, 170 pmModelOpMode mode, 171 psMaskType maskVal, 172 int dx, 173 int dy); 174 174 175 175 bool pmModelSubWithOffset(psImage *image, 176 177 178 179 180 181 176 psImage *mask, 177 pmModel *model, 178 pmModelOpMode mode, 179 psMaskType maskVal, 180 int dx, 181 int dy); 182 182 183 183 /** pmModelFitStatus() -
trunk/psModules/src/objects/pmModelUtils.c
r15562 r15843 5 5 * @author EAM, IfA 6 6 * 7 * @version $Revision: 1. 4$ $Name: not supported by cvs2svn $8 * @date $Date: 2007-1 1-10 01:09:20$7 * @version $Revision: 1.5 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2007-12-15 01:22:11 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 34 34 construct a realization of the PSF model at the object coordinates 35 35 *****************************************************************************/ 36 pmModel *pmModelFromPSF (pmModel *modelEXT, pmPSF *psf)36 pmModel *pmModelFromPSF (pmModel *modelEXT, const pmPSF *psf) 37 37 { 38 38 PS_ASSERT_PTR_NON_NULL(psf, NULL); … … 56 56 // instantiate a model for the PSF at this location with peak flux 57 57 // NOTE: psf and (Xo,Yo) are defined wrt chip coordinates 58 pmModel *pmModelFromPSFforXY ( pmPSF *psf, float Xo, float Yo, float Io)58 pmModel *pmModelFromPSFforXY (const pmPSF *psf, float Xo, float Yo, float Io) 59 59 { 60 60 PS_ASSERT_PTR_NON_NULL(psf, NULL); -
trunk/psModules/src/objects/pmModelUtils.h
r14652 r15843 5 5 * @author EAM, IfA 6 6 * 7 * @version $Revision: 1. 2$ $Name: not supported by cvs2svn $8 * @date $Date: 2007- 08-24 00:11:02$7 * @version $Revision: 1.3 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2007-12-15 01:22:11 $ 9 9 * Copyright 2007 IfA, University of Hawaii 10 10 */ … … 27 27 pmModel *pmModelFromPSF( 28 28 pmModel *model, ///< Add comment 29 pmPSF *psf///< Add comment29 const pmPSF *psf ///< Add comment 30 30 ); 31 31 32 32 pmModel *pmModelFromPSFforXY ( 33 pmPSF *psf,34 float Xo, 35 float Yo, 33 const pmPSF *psf, 34 float Xo, 35 float Yo, 36 36 float Io 37 37 ); 38 38 39 39 bool pmModelSetFlux ( 40 pmModel *model, 40 pmModel *model, 41 41 float flux 42 42 ); -
trunk/psphot/src/models/pmModel_STRAIL.c
r15000 r15843 1 2 /****************************************************************************** 3 params->data.F32[0] = So; 4 params->data.F32[1] = Zo; 5 params->data.F32[2] = Xo; 6 params->data.F32[3] = Yo; 7 params->data.F32[4] = 1 / SigmaX; 8 params->data.F32[5] = 1 / SigmaY; 9 params->data.F32[6] = Sxy; 10 params->data.F32[7] = length; 11 params->data.F32[8] = theta; 12 *****************************************************************************/ 13 14 # define PM_MODEL_FUNC 15 # define PM_MODEL_FLUX 16 # define PM_MODEL_GUESS 17 # define PM_MODEL_LIMITS 18 # define PM_MODEL_RADIUS 19 # define PM_MODEL_FROM_PSF 1 2 /****************************************************************************** 3 params->data.F32[0] = So; 4 params->data.F32[1] = Zo; 5 params->data.F32[2] = Xo; 6 params->data.F32[3] = Yo; 7 params->data.F32[4] = 1 / SigmaX; 8 params->data.F32[5] = 1 / SigmaY; 9 params->data.F32[6] = Sxy; 10 params->data.F32[7] = length; 11 params->data.F32[8] = theta; 12 *****************************************************************************/ 13 14 # define PM_MODEL_FUNC pmModelFunc_STRAIL 15 # define PM_MODEL_FLUX pmModelFlux_STRAIL 16 # define PM_MODEL_GUESS pmModelGuess_STRAIL 17 # define PM_MODEL_LIMITS pmModelLimits_STRAIL 18 # define PM_MODEL_RADIUS pmModelRadius_STRAIL 19 # define PM_MODEL_FROM_PSF pmModelFromPSF_STRAIL 20 20 # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_STRAIL 21 21 # define PM_MODEL_FIT_STATUS pmModelFitStatus_STRAIL 22 22 23 psF32 PM_MODEL_FUNC(psVector *deriv, 24 const psVector *params, 25 const psVector *x) 26 { 27 psF32 *PAR = params->data.F32; 28 29 psF32 trailLength = PAR[7]; 30 psF32 theta = PAR[8]; 31 32 psF32 x0 = PAR[2]; //streak center 33 psF32 y0 = PAR[3]; //streak center 34 35 //S values (1/sigma for x and y case, sigma for xy case) 36 psF32 sx = PAR[4]; 37 psF32 sy = PAR[5]; 38 psF32 sxy = PAR[6]; 39 40 psF32 sinT=sin(theta); 41 psF32 cosT=cos(theta); 42 psF32 sin2T=sin(2.0*theta); 43 psF32 cos2T=cos(2.0*theta); 44 45 // printf("Trying object at %4.1f,%4.1f with length %3.1f and angle %1.3f\r", x0, y0, length, theta); 46 47 //current location relative to trail center 48 psF32 X = x->data.F32[0] - x0; 49 psF32 Y = x->data.F32[1] - y0; 50 51 //x' and y' location (trail-orienter coords) 52 psF32 xs = X*cosT + Y*sinT; 53 psF32 ys = -1.0*X*sinT + Y*cosT; 54 55 //initialize variables to be changed below 56 psF32 x1 = 0; 57 psF32 y1 = 0; 58 psF32 px = 0; 59 psF32 py = 0; 60 psF32 z = 0; 61 psF32 zx = 0; 62 psF32 t = 0; 63 psF32 tx = 0; 64 psF32 r = 0; 65 psF32 rx = 0; 66 psF32 f = 0; 67 68 psF32 sxrot = 0; 69 psF32 syrot = 0; 70 psF32 sxyrot = 0; 71 psF32 dsxrot = 0; 72 psF32 dsyrot = 0; 73 psF32 dsxyrot = 0; 74 75 // psF32 Rx = 0; 76 // psF32 Ry = 0; 77 // psF32 Rxy = 0; 78 79 80 //calculate new S values (1/sigma) for rotated frame 81 psF32 sxrotsq = PS_SQR(cosT*sx) + PS_SQR(sinT*sy) + cosT*sinT*sxy; 82 psF32 syrotsq = PS_SQR(cosT*sy) + PS_SQR(sinT*sx) - cosT*sinT*sxy; 83 84 // psF32 testtwo=10.1; 85 // psF32 testone=fabsf(testtwo); 86 // fprintf (stderr, "Test: %f is the absolute value of %f?\n",testone,testtwo); 87 if (sxrotsq<0) { 88 sxrot = sqrt(-(sxrotsq)); 89 syrot = sqrt(syrotsq); 90 fprintf (stderr, "error in sxrotsq: Neg, sxrotsq=%f sx=%f sy=%f sxy=%f theta=%f\n",sxrotsq,sx,sy,sxy,theta); 91 } else if (syrotsq<0) { 92 sxrot = sqrt(sxrotsq); 93 syrot = sqrt(-(syrotsq)); 94 fprintf (stderr, "error in syrotsq: Neg, syrotsq=%f sx=%f sy=%f sxy=%f theta=%f\n",syrotsq,sx,sy,sxy,theta); 95 } else if (sxrotsq==0){ 96 sxrot = 0.01; 97 syrot = sqrt(syrotsq); 98 fprintf (stderr, "error in sxrotsq: Zero, sxrotsq=%f \n",sxrotsq); 99 } else if (syrotsq==0) { 100 syrot = 0.01; 101 sxrot = sqrt(sxrotsq); 102 fprintf (stderr, "error in syrotsq: Zero, syrotsq=%f \n",syrotsq); 103 // return(0); 104 }else { 105 sxrot = sqrt(sxrotsq); 106 syrot = sqrt(syrotsq); 107 } 108 109 sxyrot = sxy*cos2T + (PS_SQR(sy) - PS_SQR(sx))*sin2T; 110 111 if (isnan(sxrot)) { 112 fprintf (stderr, "error in sxrot syrot=%f sx=%f sy=%f sxy=%f cosT=%f sinT=%f\n",syrot,sx,sy,sxy,cosT,sinT); 113 } else if (isnan(syrot)) { 114 fprintf (stderr, "error in syrot sxrot=%f sx=%f sy=%f sxy=%f cosT=%f sinT=%f\n",sxrot,sx,sy,sxy,cosT,sinT); 115 } 116 117 //calculate length of pipe (not of trail motion) 118 psF32 length = trailLength - 2.0*2.0/sxrot; 119 120 121 if (xs < length/(-2.0)) { 122 x1 = (xs+length/2.0)*cosT - ys*sinT; //Endpoint1 123 y1 = (xs+length/2.0)*sinT + ys*cosT; //Endpoint1 124 //1.6 factor comes from by-eye testing of fits, sqrt from the later squaring of it 125 //1.6 ~= phi (golden mean)...coincidence? 126 px = sxrot*x1/sqrt(1.6); 127 py = syrot*y1; 128 129 //first find out what the falloff in the x direction is... 130 zx = 0.5*PS_SQR(px); 131 tx = 1 + zx + zx*zx/2.0 + zx*zx*zx/6.0; 132 rx = 1.0 / (tx + zx*zx*zx*zx/24.0); 133 134 //...and now in the y direction 135 z = 0.5*PS_SQR(py)+sxyrot*x1*y1; 136 t = 1 + z + z*z/2.0; 137 r = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */ 138 f = PAR[1]*rx*r + PAR[0]; 139 140 } else if (xs > length/2.0){ 141 x1 = (xs-length/2.0)*cosT - ys*sinT; //Endpoint2 142 y1 = (xs-length/2.0)*sinT + ys*cosT; //Endpoint2 143 px = sxrot*x1/sqrt(1.6); 144 py = syrot*y1; 145 zx = 0.5*PS_SQR(px); 146 tx = 1 + zx + zx*zx/2.0 + zx*zx*zx/6.0; 147 rx = 1.0 / (tx + zx*zx*zx*zx/24.0); 148 149 z = 0.5*PS_SQR(py)+sxyrot*x1*y1; 150 t = 1 + z + z*z/2.0; 151 r = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */ 152 f = PAR[1]*rx*r + PAR[0]; 153 154 if (isnan(r)) { 155 fprintf (stderr, "error in +r t=%f z=%f\n",t,z); 156 } 157 158 } else { 159 x1 = -ys*sinT; 160 y1 = ys*cosT; 161 px = sx*x1; 162 py = sy*y1; 163 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + sxy*x1*y1; 164 t = 1 + z + z*z/2.0; 165 r = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */ 166 rx = 1.0; //so that dF/dF0 can be generalized 167 f = PAR[1]*r + PAR[0]; 168 169 if (isnan(r)) { 170 fprintf (stderr, "error in midr t=%f z=%f\n",t,z); 171 } 172 } 173 174 175 //ok...so this is df/dPAR[X] 176 if (deriv != NULL) { 177 psF32 q = 1; 178 //stable 179 deriv->data.F32[0] = +1.0; 180 deriv->data.F32[1] = +r * rx; 181 182 if (isnan(deriv->data.F32[0])) { 183 fprintf (stderr, "error in deriv0\n"); 184 } else if (isnan(deriv->data.F32[1])) { 185 fprintf (stderr, "error in deriv1 r=%f rx=%f\n",r,rx); 186 } 187 188 dsxrot = 2.0*PS_SQR(sy)*sinT*cosT - 2.0*PS_SQR(sx)*sinT*cosT + sxy*(PS_SQR(cosT)-PS_SQR(sinT)); 189 dsyrot = 2.0*PS_SQR(sx)*sinT*cosT - 2.0*PS_SQR(sy)*sinT*cosT - sxy*(PS_SQR(cosT)-PS_SQR(sinT)); 190 dsxyrot = 2.0*cos2T*(PS_SQR(sy)-PS_SQR(sx)) - 2.0*sxy*sin2T; 191 192 if (isnan(dsxrot)) { 193 fprintf (stderr, "error in dsxrot\n"); 194 } else if (isnan(dsyrot)) { 195 fprintf (stderr, "error in dsyrot\n"); 196 } else if (isnan(dsxyrot)) { 197 fprintf (stderr, "error in dsxyrot\n"); 198 } 199 200 //initialize variables definied in the if statements 201 // psF32 XXX=0; 202 // psF32 YYY=0; 203 204 205 // variable over piecewise func 206 // change the endcaps to be 4th order gaussians with sxrot_fit=1.6*sxrot 207 // y' direction can ge adequately modelled by a 3rd order gaussian with syrot 208 if (xs < length/(-2.0)) { 209 210 q=PAR[1]*r*rx; 211 deriv->data.F32[2] = -q*(rx*tx/1.6*x1*PS_SQR(sxrot) + r*t*y1*sxyrot); 212 deriv->data.F32[3] = -q*r*t*(y1*PS_SQR(syrot) + x1*sxyrot); 213 deriv->data.F32[4] = -q*(rx*tx*x1*sx*PS_SQR(cosT)/1.6*(x1+2.0*cosT/sxrot) + r*t*y1*sx*sinT*(y1*sinT+2.0*PS_SQR(syrot)*PS_SQR(cosT)/(sxrot*sxrot*sxrot)) + 2.0*r*t*sx*(-1.0*x1*y1*sin2T+y1*sxyrot*(cosT*cosT*cosT)/(sxrot*sxrot*sxrot)+x1*sxyrot*(sinT*cosT*cosT)/(sxrot*sxrot*sxrot))); 214 deriv->data.F32[5] = -q*(rx*tx*x1*sy*PS_SQR(sinT)/1.6*(x1+2.0*cosT/sxrot) + r*t*y1*sy*(y1*PS_SQR(cosT)+2.0*PS_SQR(syrot)*(sinT*sinT*sinT)/(sxrot*sxrot*sxrot)) + 2.0*r*t*sy*(x1*y1*sin2T+y1*sxyrot*(sinT*sinT*cosT)/(sxrot*sxrot*sxrot)+x1*sxyrot*(sinT*sinT*sinT)/(sxrot*sxrot*sxrot))); 215 deriv->data.F32[6] = -q*(rx*tx*x1*cosT*sinT/3.2*(x1+2.0*cosT/sxrot) + r*t*y1*cosT*sinT/2.0*(-y1+2.0*PS_SQR(syrot)*sinT/(sxrot*sxrot*sxrot)) + r*t*(x1*y1*cos2T+y1*sxyrot*sinT*cosT*cosT/(sxrot*sxrot*sxrot)+x1*sxyrot*sinT*sinT*cosT/(sxrot*sxrot*sxrot))); 216 deriv->data.F32[7] = -q*(rx*tx*PS_SQR(sxrot)*x1*cosT/3.2 + r*t*PS_SQR(syrot)*y1*sinT/2.0 + r*t*sxyrot/2.0*(y1*cosT+x1*sinT)); 217 deriv->data.F32[8] = -q*(rx*tx*x1/3.2*(x1*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))+PS_SQR(sxrot)*(2.0*cosT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot) - sinT*length)) + r*t*y1/2.0*(y1*(2.0*sinT*cosT*(PS_SQR(sx)-PS_SQR(sy))-sxy*(PS_SQR(cosT)-PS_SQR(sinT)))+PS_SQR(syrot)*(2.0*sinT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)+cosT*length)) + r*t*(2.0*x1*y1*(cos2T*(PS_SQR(sy)-PS_SQR(sx))-sxy*sin2T)+y1*sxyrot/2.0*(2.0*cosT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)-length*sinT)+x1*sxyrot/2.0*(2.0*sinT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)+length*cosT))); 218 219 /* if (isnan(deriv->data.F32[2])) { 220 fprintf (stderr, "error in deriv2\n"); 221 } else if (isnan(deriv->data.F32[3])) { 222 fprintf (stderr, "error in deriv3\n"); 223 } else if (isnan(deriv->data.F32[4])) { 224 fprintf (stderr, "error in deriv4\n"); 225 } else if (isnan(deriv->data.F32[5])) { 226 fprintf (stderr, "error in deriv5\n"); 227 } else if (isnan(deriv->data.F32[6])) { 228 fprintf (stderr, "error in deriv6\n"); 229 } else if (isnan(deriv->data.F32[7])) { 230 fprintf (stderr, "error in deriv7\n"); 231 } else if (isnan(deriv->data.F32[8])) { 232 fprintf (stderr, "error in deriv8\n"); 233 } 234 */ 235 236 237 238 } else if (xs > length/2.0){ 239 240 q=PAR[1]*r*rx; 241 deriv->data.F32[2] = -q*(rx*tx/1.6*x1*PS_SQR(sxrot) + r*t*y1*sxyrot); 242 deriv->data.F32[3] = -q*r*t*(y1*PS_SQR(syrot) + x1*sxyrot); 243 deriv->data.F32[4] = -q*(rx*tx*x1*sx*PS_SQR(cosT)/1.6*(x1-2.0*cosT/sxrot) + r*t*y1*sx*sinT*(y1*sinT-2.0*PS_SQR(syrot)*PS_SQR(cosT)/(sxrot*sxrot*sxrot)) + 2.0*r*t*sx*(-1.0*x1*y1*sin2T-y1*sxyrot*(cosT*cosT*cosT)/(sxrot*sxrot*sxrot)-x1*sxyrot*(sinT*cosT*cosT)/(sxrot*sxrot*sxrot))); 244 deriv->data.F32[5] = -q*(rx*tx*x1*sy*PS_SQR(sinT)/1.6*(x1-2.0*cosT/sxrot) + r*t*y1*sy*(y1*PS_SQR(cosT)-2.0*PS_SQR(syrot)*(sinT*sinT*sinT)/(sxrot*sxrot*sxrot)) + 2.0*r*t*sy*(x1*y1*sin2T-y1*sxyrot*(sinT*sinT*cosT)/(sxrot*sxrot*sxrot)-x1*sxyrot*(sinT*sinT*sinT)/(sxrot*sxrot*sxrot))); 245 deriv->data.F32[6] = -q*(rx*tx*x1*cosT*sinT/3.2*(x1-2.0*cosT/sxrot) + r*t*y1*cosT*sinT/2.0*(-y1-2.0*PS_SQR(syrot)*sinT/(sxrot*sxrot*sxrot)) + r*t*(x1*y1*cos2T-y1*sxyrot*sinT*cosT*cosT/(sxrot*sxrot*sxrot)-x1*sxyrot*sinT*sinT*cosT/(sxrot*sxrot*sxrot))); 246 deriv->data.F32[7] = q*(rx*tx*PS_SQR(sxrot)*x1*cosT/3.2 + r*t*PS_SQR(syrot)*y1*sinT/2.0 + r*t*sxyrot/2.0*(y1*cosT+x1*sinT)); 247 deriv->data.F32[8] = -q*(rx*tx*x1/3.2*(x1*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))-PS_SQR(sxrot)*(2.0*cosT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot) - sinT*length)) + r*t*y1/2.0*(y1*(2.0*sinT*cosT*(PS_SQR(sx)-PS_SQR(sy))-sxy*(PS_SQR(cosT)-PS_SQR(sinT)))-PS_SQR(syrot)*(2.0*sinT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)+cosT*length)) + r*t*(2.0*x1*y1*(cos2T*(PS_SQR(sy)-PS_SQR(sx))-sxy*sin2T)-y1*sxyrot/2.0*(2.0*cosT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)-length*sinT)-x1*sxyrot/2.0*(2.0*sinT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)+length*cosT))); 248 249 250 /* if (isnan(deriv->data.F32[2])) { 251 fprintf (stderr, "error in deriv2\n"); 252 } else if (isnan(deriv->data.F32[3])) { 253 fprintf (stderr, "error in deriv3\n"); 254 } else if (isnan(deriv->data.F32[4])) { 255 fprintf (stderr, "error in deriv4\n"); 256 } else if (isnan(deriv->data.F32[5])) { 257 fprintf (stderr, "error in deriv5\n"); 258 } else if (isnan(deriv->data.F32[6])) { 259 fprintf (stderr, "error in deriv6\n"); 260 } else if (isnan(deriv->data.F32[7])) { 261 fprintf (stderr, "error in deriv7\n"); 262 } else if (isnan(deriv->data.F32[8])) { 263 fprintf (stderr, "error in deriv8\n"); 264 } 265 */ 266 267 } else { 268 // this does not change from before, as the y' falloff can be modelled by the standard 3rd order gaussian 269 // note difference from a pure gaussian: q = PAR[1]*r 270 q = PAR[1]*r*r*t; 271 deriv->data.F32[2] = q*(PS_SQR(sx*sinT)*x1 - PS_SQR(sy)*y1*cosT*sinT - sxy*x1*sinT*cosT + sxy*y1*PS_SQR(sinT)); 272 deriv->data.F32[3] = q*(-1*PS_SQR(sx)*x1*sinT*cosT + PS_SQR(sy*cosT)*y1 + sxy*x1*PS_SQR(cosT) - sxy*y1*sinT*cosT); 273 deriv->data.F32[4] = -q*sx*PS_SQR(x1); 274 deriv->data.F32[5] = -q*sy*PS_SQR(y1); 275 deriv->data.F32[6] = -q*x1*y1; 276 deriv->data.F32[7] = 0; 277 deriv->data.F32[8] = -q*( PS_SQR(sx)*x1*(xs*sinT - ys*cosT) + PS_SQR(sy)*y1*(xs*cosT - ys*sinT) + sxy*x1*(xs*cosT - ys*sinT) + sxy*y1*(xs*sinT - ys*cosT) ); 278 279 if (isnan(deriv->data.F32[2])) { 280 fprintf (stderr, "error in deriv2\n"); 281 } else if (isnan(deriv->data.F32[3])) { 282 fprintf (stderr, "error in deriv3\n"); 283 } else if (isnan(deriv->data.F32[4])) { 284 fprintf (stderr, "error in deriv4\n"); 285 } else if (isnan(deriv->data.F32[5])) { 286 fprintf (stderr, "error in deriv5\n"); 287 } else if (isnan(deriv->data.F32[6])) { 288 fprintf (stderr, "error in deriv6\n"); 289 } else if (isnan(deriv->data.F32[7])) { 290 fprintf (stderr, "error in deriv7\n"); 291 } else if (isnan(deriv->data.F32[8])) { 292 fprintf (stderr, "error in deriv8\n"); 293 } 294 295 296 } 297 } 298 return(f); 299 } 300 301 //fixed 23 psF32 PM_MODEL_FUNC(psVector *deriv, 24 const psVector *params, 25 const psVector *x) 26 { 27 psF32 *PAR = params->data.F32; 28 29 psF32 trailLength = PAR[7]; 30 psF32 theta = PAR[8]; 31 32 psF32 x0 = PAR[2]; //streak center 33 psF32 y0 = PAR[3]; //streak center 34 35 //S values (1/sigma for x and y case, sigma for xy case) 36 psF32 sx = PAR[4]; 37 psF32 sy = PAR[5]; 38 psF32 sxy = PAR[6]; 39 40 psF32 sinT=sin(theta); 41 psF32 cosT=cos(theta); 42 psF32 sin2T=sin(2.0*theta); 43 psF32 cos2T=cos(2.0*theta); 44 45 // printf("Trying object at %4.1f,%4.1f with length %3.1f and angle %1.3f\r", x0, y0, length, theta); 46 47 //current location relative to trail center 48 psF32 X = x->data.F32[0] - x0; 49 psF32 Y = x->data.F32[1] - y0; 50 51 //x' and y' location (trail-orienter coords) 52 psF32 xs = X*cosT + Y*sinT; 53 psF32 ys = -1.0*X*sinT + Y*cosT; 54 55 //initialize variables to be changed below 56 psF32 x1 = 0; 57 psF32 y1 = 0; 58 psF32 px = 0; 59 psF32 py = 0; 60 psF32 z = 0; 61 psF32 zx = 0; 62 psF32 t = 0; 63 psF32 tx = 0; 64 psF32 r = 0; 65 psF32 rx = 0; 66 psF32 f = 0; 67 68 psF32 sxrot = 0; 69 psF32 syrot = 0; 70 psF32 sxyrot = 0; 71 psF32 dsxrot = 0; 72 psF32 dsyrot = 0; 73 psF32 dsxyrot = 0; 74 75 // psF32 Rx = 0; 76 // psF32 Ry = 0; 77 // psF32 Rxy = 0; 78 79 80 //calculate new S values (1/sigma) for rotated frame 81 psF32 sxrotsq = PS_SQR(cosT*sx) + PS_SQR(sinT*sy) + cosT*sinT*sxy; 82 psF32 syrotsq = PS_SQR(cosT*sy) + PS_SQR(sinT*sx) - cosT*sinT*sxy; 83 84 // psF32 testtwo=10.1; 85 // psF32 testone=fabsf(testtwo); 86 // fprintf (stderr, "Test: %f is the absolute value of %f?\n",testone,testtwo); 87 if (sxrotsq<0) { 88 sxrot = sqrt(-(sxrotsq)); 89 syrot = sqrt(syrotsq); 90 fprintf (stderr, "error in sxrotsq: Neg, sxrotsq=%f sx=%f sy=%f sxy=%f theta=%f\n",sxrotsq,sx,sy,sxy,theta); 91 } else if (syrotsq<0) { 92 sxrot = sqrt(sxrotsq); 93 syrot = sqrt(-(syrotsq)); 94 fprintf (stderr, "error in syrotsq: Neg, syrotsq=%f sx=%f sy=%f sxy=%f theta=%f\n",syrotsq,sx,sy,sxy,theta); 95 } else if (sxrotsq==0){ 96 sxrot = 0.01; 97 syrot = sqrt(syrotsq); 98 fprintf (stderr, "error in sxrotsq: Zero, sxrotsq=%f \n",sxrotsq); 99 } else if (syrotsq==0) { 100 syrot = 0.01; 101 sxrot = sqrt(sxrotsq); 102 fprintf (stderr, "error in syrotsq: Zero, syrotsq=%f \n",syrotsq); 103 // return(0); 104 }else { 105 sxrot = sqrt(sxrotsq); 106 syrot = sqrt(syrotsq); 107 } 108 109 sxyrot = sxy*cos2T + (PS_SQR(sy) - PS_SQR(sx))*sin2T; 110 111 if (isnan(sxrot)) { 112 fprintf (stderr, "error in sxrot syrot=%f sx=%f sy=%f sxy=%f cosT=%f sinT=%f\n",syrot,sx,sy,sxy,cosT,sinT); 113 } else if (isnan(syrot)) { 114 fprintf (stderr, "error in syrot sxrot=%f sx=%f sy=%f sxy=%f cosT=%f sinT=%f\n",sxrot,sx,sy,sxy,cosT,sinT); 115 } 116 117 //calculate length of pipe (not of trail motion) 118 psF32 length = trailLength - 2.0*2.0/sxrot; 119 120 121 if (xs < length/(-2.0)) { 122 x1 = (xs+length/2.0)*cosT - ys*sinT; //Endpoint1 123 y1 = (xs+length/2.0)*sinT + ys*cosT; //Endpoint1 124 //1.6 factor comes from by-eye testing of fits, sqrt from the later squaring of it 125 //1.6 ~= phi (golden mean)...coincidence? 126 px = sxrot*x1/sqrt(1.6); 127 py = syrot*y1; 128 129 //first find out what the falloff in the x direction is... 130 zx = 0.5*PS_SQR(px); 131 tx = 1 + zx + zx*zx/2.0 + zx*zx*zx/6.0; 132 rx = 1.0 / (tx + zx*zx*zx*zx/24.0); 133 134 //...and now in the y direction 135 z = 0.5*PS_SQR(py)+sxyrot*x1*y1; 136 t = 1 + z + z*z/2.0; 137 r = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */ 138 f = PAR[1]*rx*r + PAR[0]; 139 140 } else if (xs > length/2.0){ 141 x1 = (xs-length/2.0)*cosT - ys*sinT; //Endpoint2 142 y1 = (xs-length/2.0)*sinT + ys*cosT; //Endpoint2 143 px = sxrot*x1/sqrt(1.6); 144 py = syrot*y1; 145 zx = 0.5*PS_SQR(px); 146 tx = 1 + zx + zx*zx/2.0 + zx*zx*zx/6.0; 147 rx = 1.0 / (tx + zx*zx*zx*zx/24.0); 148 149 z = 0.5*PS_SQR(py)+sxyrot*x1*y1; 150 t = 1 + z + z*z/2.0; 151 r = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */ 152 f = PAR[1]*rx*r + PAR[0]; 153 154 if (isnan(r)) { 155 fprintf (stderr, "error in +r t=%f z=%f\n",t,z); 156 } 157 158 } else { 159 x1 = -ys*sinT; 160 y1 = ys*cosT; 161 px = sx*x1; 162 py = sy*y1; 163 z = 0.5*PS_SQR(px) + 0.5*PS_SQR(py) + sxy*x1*y1; 164 t = 1 + z + z*z/2.0; 165 r = 1.0 / (t + z*z*z/6.0); /* exp (-Z) */ 166 rx = 1.0; //so that dF/dF0 can be generalized 167 f = PAR[1]*r + PAR[0]; 168 169 if (isnan(r)) { 170 fprintf (stderr, "error in midr t=%f z=%f\n",t,z); 171 } 172 } 173 174 175 //ok...so this is df/dPAR[X] 176 if (deriv != NULL) { 177 psF32 q = 1; 178 //stable 179 deriv->data.F32[0] = +1.0; 180 deriv->data.F32[1] = +r * rx; 181 182 if (isnan(deriv->data.F32[0])) { 183 fprintf (stderr, "error in deriv0\n"); 184 } else if (isnan(deriv->data.F32[1])) { 185 fprintf (stderr, "error in deriv1 r=%f rx=%f\n",r,rx); 186 } 187 188 dsxrot = 2.0*PS_SQR(sy)*sinT*cosT - 2.0*PS_SQR(sx)*sinT*cosT + sxy*(PS_SQR(cosT)-PS_SQR(sinT)); 189 dsyrot = 2.0*PS_SQR(sx)*sinT*cosT - 2.0*PS_SQR(sy)*sinT*cosT - sxy*(PS_SQR(cosT)-PS_SQR(sinT)); 190 dsxyrot = 2.0*cos2T*(PS_SQR(sy)-PS_SQR(sx)) - 2.0*sxy*sin2T; 191 192 if (isnan(dsxrot)) { 193 fprintf (stderr, "error in dsxrot\n"); 194 } else if (isnan(dsyrot)) { 195 fprintf (stderr, "error in dsyrot\n"); 196 } else if (isnan(dsxyrot)) { 197 fprintf (stderr, "error in dsxyrot\n"); 198 } 199 200 //initialize variables definied in the if statements 201 // psF32 XXX=0; 202 // psF32 YYY=0; 203 204 205 // variable over piecewise func 206 // change the endcaps to be 4th order gaussians with sxrot_fit=1.6*sxrot 207 // y' direction can ge adequately modelled by a 3rd order gaussian with syrot 208 if (xs < length/(-2.0)) { 209 210 q=PAR[1]*r*rx; 211 deriv->data.F32[2] = -q*(rx*tx/1.6*x1*PS_SQR(sxrot) + r*t*y1*sxyrot); 212 deriv->data.F32[3] = -q*r*t*(y1*PS_SQR(syrot) + x1*sxyrot); 213 deriv->data.F32[4] = -q*(rx*tx*x1*sx*PS_SQR(cosT)/1.6*(x1+2.0*cosT/sxrot) + r*t*y1*sx*sinT*(y1*sinT+2.0*PS_SQR(syrot)*PS_SQR(cosT)/(sxrot*sxrot*sxrot)) + 2.0*r*t*sx*(-1.0*x1*y1*sin2T+y1*sxyrot*(cosT*cosT*cosT)/(sxrot*sxrot*sxrot)+x1*sxyrot*(sinT*cosT*cosT)/(sxrot*sxrot*sxrot))); 214 deriv->data.F32[5] = -q*(rx*tx*x1*sy*PS_SQR(sinT)/1.6*(x1+2.0*cosT/sxrot) + r*t*y1*sy*(y1*PS_SQR(cosT)+2.0*PS_SQR(syrot)*(sinT*sinT*sinT)/(sxrot*sxrot*sxrot)) + 2.0*r*t*sy*(x1*y1*sin2T+y1*sxyrot*(sinT*sinT*cosT)/(sxrot*sxrot*sxrot)+x1*sxyrot*(sinT*sinT*sinT)/(sxrot*sxrot*sxrot))); 215 deriv->data.F32[6] = -q*(rx*tx*x1*cosT*sinT/3.2*(x1+2.0*cosT/sxrot) + r*t*y1*cosT*sinT/2.0*(-y1+2.0*PS_SQR(syrot)*sinT/(sxrot*sxrot*sxrot)) + r*t*(x1*y1*cos2T+y1*sxyrot*sinT*cosT*cosT/(sxrot*sxrot*sxrot)+x1*sxyrot*sinT*sinT*cosT/(sxrot*sxrot*sxrot))); 216 deriv->data.F32[7] = -q*(rx*tx*PS_SQR(sxrot)*x1*cosT/3.2 + r*t*PS_SQR(syrot)*y1*sinT/2.0 + r*t*sxyrot/2.0*(y1*cosT+x1*sinT)); 217 deriv->data.F32[8] = -q*(rx*tx*x1/3.2*(x1*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))+PS_SQR(sxrot)*(2.0*cosT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot) - sinT*length)) + r*t*y1/2.0*(y1*(2.0*sinT*cosT*(PS_SQR(sx)-PS_SQR(sy))-sxy*(PS_SQR(cosT)-PS_SQR(sinT)))+PS_SQR(syrot)*(2.0*sinT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)+cosT*length)) + r*t*(2.0*x1*y1*(cos2T*(PS_SQR(sy)-PS_SQR(sx))-sxy*sin2T)+y1*sxyrot/2.0*(2.0*cosT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)-length*sinT)+x1*sxyrot/2.0*(2.0*sinT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)+length*cosT))); 218 219 /* if (isnan(deriv->data.F32[2])) { 220 fprintf (stderr, "error in deriv2\n"); 221 } else if (isnan(deriv->data.F32[3])) { 222 fprintf (stderr, "error in deriv3\n"); 223 } else if (isnan(deriv->data.F32[4])) { 224 fprintf (stderr, "error in deriv4\n"); 225 } else if (isnan(deriv->data.F32[5])) { 226 fprintf (stderr, "error in deriv5\n"); 227 } else if (isnan(deriv->data.F32[6])) { 228 fprintf (stderr, "error in deriv6\n"); 229 } else if (isnan(deriv->data.F32[7])) { 230 fprintf (stderr, "error in deriv7\n"); 231 } else if (isnan(deriv->data.F32[8])) { 232 fprintf (stderr, "error in deriv8\n"); 233 } 234 */ 235 236 237 238 } else if (xs > length/2.0){ 239 240 q=PAR[1]*r*rx; 241 deriv->data.F32[2] = -q*(rx*tx/1.6*x1*PS_SQR(sxrot) + r*t*y1*sxyrot); 242 deriv->data.F32[3] = -q*r*t*(y1*PS_SQR(syrot) + x1*sxyrot); 243 deriv->data.F32[4] = -q*(rx*tx*x1*sx*PS_SQR(cosT)/1.6*(x1-2.0*cosT/sxrot) + r*t*y1*sx*sinT*(y1*sinT-2.0*PS_SQR(syrot)*PS_SQR(cosT)/(sxrot*sxrot*sxrot)) + 2.0*r*t*sx*(-1.0*x1*y1*sin2T-y1*sxyrot*(cosT*cosT*cosT)/(sxrot*sxrot*sxrot)-x1*sxyrot*(sinT*cosT*cosT)/(sxrot*sxrot*sxrot))); 244 deriv->data.F32[5] = -q*(rx*tx*x1*sy*PS_SQR(sinT)/1.6*(x1-2.0*cosT/sxrot) + r*t*y1*sy*(y1*PS_SQR(cosT)-2.0*PS_SQR(syrot)*(sinT*sinT*sinT)/(sxrot*sxrot*sxrot)) + 2.0*r*t*sy*(x1*y1*sin2T-y1*sxyrot*(sinT*sinT*cosT)/(sxrot*sxrot*sxrot)-x1*sxyrot*(sinT*sinT*sinT)/(sxrot*sxrot*sxrot))); 245 deriv->data.F32[6] = -q*(rx*tx*x1*cosT*sinT/3.2*(x1-2.0*cosT/sxrot) + r*t*y1*cosT*sinT/2.0*(-y1-2.0*PS_SQR(syrot)*sinT/(sxrot*sxrot*sxrot)) + r*t*(x1*y1*cos2T-y1*sxyrot*sinT*cosT*cosT/(sxrot*sxrot*sxrot)-x1*sxyrot*sinT*sinT*cosT/(sxrot*sxrot*sxrot))); 246 deriv->data.F32[7] = q*(rx*tx*PS_SQR(sxrot)*x1*cosT/3.2 + r*t*PS_SQR(syrot)*y1*sinT/2.0 + r*t*sxyrot/2.0*(y1*cosT+x1*sinT)); 247 deriv->data.F32[8] = -q*(rx*tx*x1/3.2*(x1*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))-PS_SQR(sxrot)*(2.0*cosT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot) - sinT*length)) + r*t*y1/2.0*(y1*(2.0*sinT*cosT*(PS_SQR(sx)-PS_SQR(sy))-sxy*(PS_SQR(cosT)-PS_SQR(sinT)))-PS_SQR(syrot)*(2.0*sinT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)+cosT*length)) + r*t*(2.0*x1*y1*(cos2T*(PS_SQR(sy)-PS_SQR(sx))-sxy*sin2T)-y1*sxyrot/2.0*(2.0*cosT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)-length*sinT)-x1*sxyrot/2.0*(2.0*sinT*(2.0*sinT*cosT*(PS_SQR(sy)-PS_SQR(sx))+sxy*(PS_SQR(cosT)-PS_SQR(sinT)))/(sxrot*sxrot*sxrot)+length*cosT))); 248 249 250 /* if (isnan(deriv->data.F32[2])) { 251 fprintf (stderr, "error in deriv2\n"); 252 } else if (isnan(deriv->data.F32[3])) { 253 fprintf (stderr, "error in deriv3\n"); 254 } else if (isnan(deriv->data.F32[4])) { 255 fprintf (stderr, "error in deriv4\n"); 256 } else if (isnan(deriv->data.F32[5])) { 257 fprintf (stderr, "error in deriv5\n"); 258 } else if (isnan(deriv->data.F32[6])) { 259 fprintf (stderr, "error in deriv6\n"); 260 } else if (isnan(deriv->data.F32[7])) { 261 fprintf (stderr, "error in deriv7\n"); 262 } else if (isnan(deriv->data.F32[8])) { 263 fprintf (stderr, "error in deriv8\n"); 264 } 265 */ 266 267 } else { 268 // this does not change from before, as the y' falloff can be modelled by the standard 3rd order gaussian 269 // note difference from a pure gaussian: q = PAR[1]*r 270 q = PAR[1]*r*r*t; 271 deriv->data.F32[2] = q*(PS_SQR(sx*sinT)*x1 - PS_SQR(sy)*y1*cosT*sinT - sxy*x1*sinT*cosT + sxy*y1*PS_SQR(sinT)); 272 deriv->data.F32[3] = q*(-1*PS_SQR(sx)*x1*sinT*cosT + PS_SQR(sy*cosT)*y1 + sxy*x1*PS_SQR(cosT) - sxy*y1*sinT*cosT); 273 deriv->data.F32[4] = -q*sx*PS_SQR(x1); 274 deriv->data.F32[5] = -q*sy*PS_SQR(y1); 275 deriv->data.F32[6] = -q*x1*y1; 276 deriv->data.F32[7] = 0; 277 deriv->data.F32[8] = -q*( PS_SQR(sx)*x1*(xs*sinT - ys*cosT) + PS_SQR(sy)*y1*(xs*cosT - ys*sinT) + sxy*x1*(xs*cosT - ys*sinT) + sxy*y1*(xs*sinT - ys*cosT) ); 278 279 if (isnan(deriv->data.F32[2])) { 280 fprintf (stderr, "error in deriv2\n"); 281 } else if (isnan(deriv->data.F32[3])) { 282 fprintf (stderr, "error in deriv3\n"); 283 } else if (isnan(deriv->data.F32[4])) { 284 fprintf (stderr, "error in deriv4\n"); 285 } else if (isnan(deriv->data.F32[5])) { 286 fprintf (stderr, "error in deriv5\n"); 287 } else if (isnan(deriv->data.F32[6])) { 288 fprintf (stderr, "error in deriv6\n"); 289 } else if (isnan(deriv->data.F32[7])) { 290 fprintf (stderr, "error in deriv7\n"); 291 } else if (isnan(deriv->data.F32[8])) { 292 fprintf (stderr, "error in deriv8\n"); 293 } 294 295 296 } 297 } 298 return(f); 299 } 300 301 //fixed 302 302 // XXX this needs to apply the axis ratio limits to prevent avoid solutions 303 303 # define AR_MAX 20.0 304 304 # define AR_RATIO 0.99 305 305 bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) { 306 306 307 307 float beta_lim = 0; 308 308 float params_min = 0; … … 310 310 float f1, f2, q1; 311 311 float q2 = 0; 312 312 313 313 // we need to calculate the limits for SXY specially 314 314 if (nParam == PM_PAR_SXY) { 315 316 317 318 319 315 f1 = 1.0 / PS_SQR(params[PM_PAR_SYY]) + 1.0 / PS_SQR(params[PM_PAR_SXX]); 316 f2 = 1.0 / PS_SQR(params[PM_PAR_SYY]) - 1.0 / PS_SQR(params[PM_PAR_SXX]); 317 q1 = PS_SQR(f1)*AR_RATIO - PS_SQR(f2); 318 assert (q1 > 0); 319 q2 = 0.5*sqrt (q1); 320 320 } 321 321 322 322 switch (mode) { 323 323 case PS_MINIMIZE_BETA_LIMIT: 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 324 switch (nParam) { 325 case PM_PAR_SKY: beta_lim = 1000; break; 326 case PM_PAR_I0: beta_lim = 10000; break; // too small? 327 case PM_PAR_XPOS: beta_lim = 50; break; 328 case PM_PAR_YPOS: beta_lim = 50; break; 329 case PM_PAR_SXX: beta_lim = 0.5; break; 330 case PM_PAR_SYY: beta_lim = 0.5; break; 331 case PM_PAR_SXY: beta_lim = 1.0; break; // set this to q2? 332 case 7: beta_lim = 10.0; break; 333 case 8: beta_lim = M_PI/6.0; break; 334 335 default: 336 psAbort("invalid parameter %d for beta test", nParam); 337 } 338 if (fabs(beta[nParam]) > fabs(beta_lim)) { 339 beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim); 340 return false; 341 } 342 return true; 343 343 case PS_MINIMIZE_PARAM_MIN: 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 344 switch (nParam) { 345 case PM_PAR_SKY: params_min = -1000; break; 346 case PM_PAR_I0: params_min = 0; break; 347 case PM_PAR_XPOS: params_min = -100; break; 348 case PM_PAR_YPOS: params_min = -100; break; 349 case PM_PAR_SXX: params_min = 0.5; break; 350 case PM_PAR_SYY: params_min = 0.5; break; 351 case PM_PAR_SXY: params_min = -5.0; break; // set this to -q2? 352 case 7: params_min = 0; break; 353 case 8: params_min = -1*M_PI; break; 354 355 default: 356 psAbort("invalid parameter %d for param min test", nParam); 357 } 358 if (params[nParam] < params_min) { 359 params[nParam] = params_min; 360 return false; 361 } 362 return true; 363 363 case PS_MINIMIZE_PARAM_MAX: 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 364 switch (nParam) { 365 case PM_PAR_SKY: params_max = 1e5; break; 366 case PM_PAR_I0: params_max = 1e8; break; 367 case PM_PAR_XPOS: params_max = 1e4; break; 368 case PM_PAR_YPOS: params_max = 1e4; break; 369 case PM_PAR_SXX: params_max = 100; break; 370 case PM_PAR_SYY: params_max = 100; break; 371 case PM_PAR_SXY: params_max = +q2; break; 372 case 7: params_max = 150; break; 373 case 8: params_max = M_PI; break; 374 default: 375 psAbort("invalid parameter %d for param max test", nParam); 376 } 377 if (params[nParam] > params_max) { 378 params[nParam] = params_max; 379 return false; 380 } 381 return true; 382 382 default: 383 383 psAbort("invalid choice for limits"); 384 384 } 385 385 psAbort("should not reach here"); … … 387 387 } 388 388 389 //fixed 390 psF64 PM_MODEL_FLUX(const psVector *params) 391 { 392 float f, norm, z; 393 394 psF32 *PAR = params->data.F32; 395 396 psF64 A1 = PS_SQR(PAR[4]); 397 psF64 A2 = PS_SQR(PAR[5]); 398 psF64 A3 = PS_SQR(PAR[6]); 399 psF32 Rx=2./PS_SQR(PAR[4]); 400 psF32 Ry=2./PS_SQR(PAR[5]); 401 psF32 Rxy=PAR[6]; 402 403 404 psF32 theta = PAR[8]; 405 psF32 sinT=sin(theta); 406 psF32 cosT=cos(theta); 407 408 psF32 Syrot = ( PS_SQR(sinT)/Rx + PS_SQR(cosT)/Ry - Rxy*sinT*cosT ); //rotated sigma y 409 410 psF64 A4 = Syrot*PAR[7]; 411 412 psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3) + A4; 413 // Area is equivalent to 2 pi sigma^2 + rectangle 414 415 // the area needs to be multiplied by the integral of f(z) 416 norm = 0.0; 417 for (z = 0.005; z < 50; z += 0.01) { 418 f = 1.0 / (1 + z + z*z/2 + z*z*z/6); 419 norm += f; 420 } 421 norm *= 0.01; 422 423 psF64 Flux = params->data.F32[1] * Area * norm; 424 425 return(Flux); 426 } 427 428 // define this function so it never returns Inf or NaN 429 // also prevent 0 returns, and just send a v small number 430 // return the radius which yields the requested flux 431 432 //fixed, but need to change how it is called to accomodate 2 radii 433 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) 434 { 435 if (flux <= 0) return (1.0); 436 if (params->data.F32[1] <= 0) return (1.0); 437 if (flux >= params->data.F32[1]) return (1.0); 438 439 psF32 *PAR = params->data.F32; 440 psF32 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); 441 442 psF32 theta = PAR[8]; 443 psF32 sinT=sin(theta); 444 psF32 cosT=cos(theta); 445 psF32 Rx=2./PS_SQR(PAR[4]); 446 psF32 Ry=2./PS_SQR(PAR[5]); 447 psF32 Rxy=PAR[6]; 448 psF32 length=PAR[7]; 449 450 psF32 Syrot = ( PS_SQR(sinT)/Rx + PS_SQR(cosT)/Ry - Rxy*sinT*cosT ); //rotated sigma y 451 452 psF64 radius = 0; 453 if (flux > 0){ 454 psF64 radius0 = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); 455 psF64 radius1 = Syrot * sqrt (2.0 * log(params->data.F32[1] / flux)); 456 457 if (radius0 > radius1) { 458 radius=radius0+length/2.0; 459 } else { 460 radius=radius1+length/2.0; 461 } 462 } else { 463 radius = 1000; 464 } 465 466 if (radius < 0.01){ 467 radius = 0.01; 468 } 469 470 if (isnan(radius)) { 471 fprintf (stderr, "error in code\n"); 472 } 473 return (radius); 474 } 475 476 //fixed I think...no good way of guessing as far as I can tell 389 //fixed 390 psF64 PM_MODEL_FLUX(const psVector *params) 391 { 392 float f, norm, z; 393 394 psF32 *PAR = params->data.F32; 395 396 psF64 A1 = PS_SQR(PAR[4]); 397 psF64 A2 = PS_SQR(PAR[5]); 398 psF64 A3 = PS_SQR(PAR[6]); 399 psF32 Rx=2./PS_SQR(PAR[4]); 400 psF32 Ry=2./PS_SQR(PAR[5]); 401 psF32 Rxy=PAR[6]; 402 403 404 psF32 theta = PAR[8]; 405 psF32 sinT=sin(theta); 406 psF32 cosT=cos(theta); 407 408 psF32 Syrot = ( PS_SQR(sinT)/Rx + PS_SQR(cosT)/Ry - Rxy*sinT*cosT ); //rotated sigma y 409 410 psF64 A4 = Syrot*PAR[7]; 411 412 psF64 Area = 2.0 * M_PI / sqrt(A1*A2 - A3) + A4; 413 // Area is equivalent to 2 pi sigma^2 + rectangle 414 415 // the area needs to be multiplied by the integral of f(z) 416 norm = 0.0; 417 for (z = 0.005; z < 50; z += 0.01) { 418 f = 1.0 / (1 + z + z*z/2 + z*z*z/6); 419 norm += f; 420 } 421 norm *= 0.01; 422 423 psF64 Flux = params->data.F32[1] * Area * norm; 424 425 return(Flux); 426 } 427 428 // define this function so it never returns Inf or NaN 429 // also prevent 0 returns, and just send a v small number 430 // return the radius which yields the requested flux 431 432 //fixed, but need to change how it is called to accomodate 2 radii 433 psF64 PM_MODEL_RADIUS (const psVector *params, psF64 flux) 434 { 435 if (flux <= 0) return (1.0); 436 if (params->data.F32[1] <= 0) return (1.0); 437 if (flux >= params->data.F32[1]) return (1.0); 438 439 psF32 *PAR = params->data.F32; 440 psF32 sigma = sqrt(2.0) * hypot (1.0 / params->data.F32[4], 1.0 / params->data.F32[5]); 441 442 psF32 theta = PAR[8]; 443 psF32 sinT=sin(theta); 444 psF32 cosT=cos(theta); 445 psF32 Rx=2./PS_SQR(PAR[4]); 446 psF32 Ry=2./PS_SQR(PAR[5]); 447 psF32 Rxy=PAR[6]; 448 psF32 length=PAR[7]; 449 450 psF32 Syrot = ( PS_SQR(sinT)/Rx + PS_SQR(cosT)/Ry - Rxy*sinT*cosT ); //rotated sigma y 451 452 psF64 radius = 0; 453 if (flux > 0){ 454 psF64 radius0 = sigma * sqrt (2.0 * log(params->data.F32[1] / flux)); 455 psF64 radius1 = Syrot * sqrt (2.0 * log(params->data.F32[1] / flux)); 456 457 if (radius0 > radius1) { 458 radius=radius0+length/2.0; 459 } else { 460 radius=radius1+length/2.0; 461 } 462 } else { 463 radius = 1000; 464 } 465 466 if (radius < 0.01){ 467 radius = 0.01; 468 } 469 470 if (isnan(radius)) { 471 fprintf (stderr, "error in code\n"); 472 } 473 return (radius); 474 } 475 476 //fixed I think...no good way of guessing as far as I can tell 477 477 bool PM_MODEL_GUESS (pmModel *model, pmSource *source) 478 { 479 pmMoments *Smoments = source->moments; 480 psF32 *params = model->params->data.F32; 481 482 psEllipseAxes axes; 483 psEllipseShape shape; 484 psEllipseMoments moments; 485 486 moments.x2 = PS_SQR(Smoments->Sx); 487 moments.y2 = PS_SQR(Smoments->Sy); 488 moments.xy = Smoments->Sxy; 489 //sometimes these moment inputs are zero...why? 490 491 // solve the math to go from Moments To Shape 478 { 479 pmMoments *Smoments = source->moments; 480 psF32 *params = model->params->data.F32; 481 482 psEllipseAxes axes; 483 psEllipseShape shape; 484 psEllipseMoments moments; 485 486 moments.x2 = PS_SQR(Smoments->Sx); 487 moments.y2 = PS_SQR(Smoments->Sy); 488 moments.xy = Smoments->Sxy; 489 //sometimes these moment inputs are zero...why? 490 491 // solve the math to go from Moments To Shape 492 492 // limit the axis ratio to 20 (a guess) 493 axes = psEllipseMomentsToAxes(moments, 20.0); 494 shape = psEllipseAxesToShape(axes); 495 496 params[0] = Smoments->Sky; 497 params[1] = Smoments->Peak - Smoments->Sky; 498 params[2] = Smoments->x; 499 params[3] = Smoments->y; 500 params[7] = 2 * axes.major; 501 params[8] = axes.theta; 502 503 if (moments.x2 == 0 || moments.y2 == 0){ 504 params[4] = 2.0; 505 params[5] = 2.0; 506 params[6] = 0.0; 507 } else { 508 params[4] = 1.0 / shape.sx; 509 params[5] = 1.0 / shape.sy; 510 params[6] = shape.sxy; 511 } 512 513 // printf("Who is NaN? momx: %4.3f momy: %4.3f momxy: %4.3f\n", moments.x2,moments.y2, moments.xy); 514 515 return(true); 516 } 517 518 //fixed 519 bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf)493 axes = psEllipseMomentsToAxes(moments, 20.0); 494 shape = psEllipseAxesToShape(axes); 495 496 params[0] = Smoments->Sky; 497 params[1] = Smoments->Peak - Smoments->Sky; 498 params[2] = Smoments->x; 499 params[3] = Smoments->y; 500 params[7] = 2 * axes.major; 501 params[8] = axes.theta; 502 503 if (moments.x2 == 0 || moments.y2 == 0){ 504 params[4] = 2.0; 505 params[5] = 2.0; 506 params[6] = 0.0; 507 } else { 508 params[4] = 1.0 / shape.sx; 509 params[5] = 1.0 / shape.sy; 510 params[6] = shape.sxy; 511 } 512 513 // printf("Who is NaN? momx: %4.3f momy: %4.3f momxy: %4.3f\n", moments.x2,moments.y2, moments.xy); 514 515 return(true); 516 } 517 518 //fixed 519 bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) 520 520 { 521 psF32 *out = modelPSF->params->data.F32; 522 psF32 *in = modelFLT->params->data.F32; 523 524 out[0] = in[0]; 525 out[1] = in[1]; 526 out[2] = in[2]; 527 out[3] = in[3]; 528 out[7] = in[7]; 529 out[8] = in[8]; 530 531 for (int i = 4; i < 7; i++) { 532 pmTrend2D *trend = psf->params->data[i-4]; 533 out[i] = pmTrend2DEval (trend, out[2], out[3]); 534 } 535 return(true); 536 } 537 521 psF32 *out = modelPSF->params->data.F32; 522 psF32 *in = modelFLT->params->data.F32; 523 524 out[0] = in[0]; 525 out[1] = in[1]; 526 out[2] = in[2]; 527 out[3] = in[3]; 528 out[7] = in[7]; 529 out[8] = in[8]; 530 531 for (int i = 4; i < 7; i++) { 532 pmTrend2D *trend = psf->params->data[i-4]; 533 out[i] = pmTrend2DEval (trend, out[2], out[3]); 534 } 535 return(true); 536 } 537 538 538 // construct the PSF model from the FLT model and the psf 539 539 // XXX is this sufficiently general do be a global function, not a pmModelClass function? 540 bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, pmPSF *psf, float Xo, float Yo, float Io)540 bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io) 541 541 { 542 542 psF32 *PAR = model->params->data.F32; … … 550 550 PAR[PM_PAR_XPOS] = Xo; 551 551 PAR[PM_PAR_YPOS] = Yo; 552 552 553 553 // supply the model-fitted parameters, or copy from the input 554 554 for (int i = 0; i < psf->params->n; i++) { 555 556 557 555 if (i == PM_PAR_SKY) continue; 556 pmTrend2D *trend = psf->params->data[i]; 557 PAR[i] = pmTrend2DEval(trend, Xo, Yo); 558 558 } 559 559 … … 562 562 // XXX user-defined value for limit? 563 563 if (!pmPSF_FitToModel (PAR, 0.1)) { 564 565 564 psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); 565 return false; 566 566 } 567 567 … … 573 573 continue; 574 574 575 575 bool status = true; 576 576 status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, PAR, NULL); 577 577 status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, PAR, NULL); 578 579 580 581 578 if (!status) { 579 psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", Xo, Yo); 580 model->flags |= PM_MODEL_STATUS_LIMITS; 581 } 582 582 } 583 583 return(true); 584 584 } 585 585 586 //done I think 586 //done I think 587 587 bool PM_MODEL_FIT_STATUS (pmModel *model) 588 588 { 589 psF32 dP; 590 bool status; 591 592 psF32 *PAR = model->params->data.F32; 593 psF32 *dPAR = model->dparams->data.F32; 594 595 dP = 0; 596 dP += PS_SQR(dPAR[4] / PAR[4]); 597 dP += PS_SQR(dPAR[5] / PAR[5]); 598 dP = sqrt (dP); 599 600 status = true; 601 status &= (dP < 0.5); 602 status &= (PAR[1] > 0); 603 status &= ((dPAR[1]/PAR[1]) < 0.5); 604 // status &= ((dPAR[7]/PAR[7]) < 0.5); 605 // status &= ((dPAR[8]/PAR[8]) < 0.5); 606 607 if (status) return true; 608 return false; 609 } 589 psF32 dP; 590 bool status; 591 592 psF32 *PAR = model->params->data.F32; 593 psF32 *dPAR = model->dparams->data.F32; 594 595 dP = 0; 596 dP += PS_SQR(dPAR[4] / PAR[4]); 597 dP += PS_SQR(dPAR[5] / PAR[5]); 598 dP = sqrt (dP); 599 600 status = true; 601 status &= (dP < 0.5); 602 status &= (PAR[1] > 0); 603 status &= ((dPAR[1]/PAR[1]) < 0.5); 604 // status &= ((dPAR[7]/PAR[7]) < 0.5); 605 // status &= ((dPAR[8]/PAR[8]) < 0.5); 606 607 if (status) return true; 608 return false; 609 } 610 610 611 611 # undef PM_MODEL_FUNC -
trunk/psphot/src/models/pmModel_TEST1.c
r15000 r15843 16 16 *****************************************************************************/ 17 17 18 # define PM_MODEL_FUNC 19 # define PM_MODEL_FLUX 20 # define PM_MODEL_GUESS 21 # define PM_MODEL_LIMITS 22 # define PM_MODEL_RADIUS 23 # define PM_MODEL_FROM_PSF 18 # define PM_MODEL_FUNC pmModelFunc_TEST1 19 # define PM_MODEL_FLUX pmModelFlux_TEST1 20 # define PM_MODEL_GUESS pmModelGuess_TEST1 21 # define PM_MODEL_LIMITS pmModelLimits_TEST1 22 # define PM_MODEL_RADIUS pmModelRadius_TEST1 23 # define PM_MODEL_FROM_PSF pmModelFromPSF_TEST1 24 24 # define PM_MODEL_PARAMS_FROM_PSF pmModelParamsFromPSF_TEST1 25 25 # define PM_MODEL_FIT_STATUS pmModelFitStatus_TEST1 … … 27 27 // the model is a function of the pixel coordinate (pixcoord[0,1] = x,y) 28 28 psF32 PM_MODEL_FUNC(psVector *deriv, 29 30 29 const psVector *params, 30 const psVector *pixcoord) 31 31 { 32 32 psF32 *PAR = params->data.F32; … … 57 57 58 58 // define the parameter limits 59 bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) 59 bool PM_MODEL_LIMITS (psMinConstraintMode mode, int nParam, float *params, float *beta) 60 60 { 61 61 float beta_lim = 0; … … 65 65 switch (mode) { 66 66 case PS_MINIMIZE_BETA_LIMIT: 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 67 switch (nParam) { 68 case PM_PAR_SKY: beta_lim = 1000; break; 69 case PM_PAR_I0: beta_lim = 3e6; break; 70 case PM_PAR_XPOS: beta_lim = 5; break; 71 case PM_PAR_YPOS: beta_lim = 5; break; 72 case PM_PAR_SXX: beta_lim = 0.5; break; 73 case PM_PAR_SYY: beta_lim = 0.5; break; 74 case PM_PAR_SXY: beta_lim = 0.5; break; 75 default: 76 psAbort("invalid parameter %d for beta test", nParam); 77 } 78 if (fabs(beta[nParam]) > fabs(beta_lim)) { 79 beta[nParam] = (beta[nParam] > 0) ? fabs(beta_lim) : -fabs(beta_lim); 80 return false; 81 } 82 return true; 83 83 case PS_MINIMIZE_PARAM_MIN: 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 84 switch (nParam) { 85 case PM_PAR_SKY: params_min = -1000; break; 86 case PM_PAR_I0: params_min = 0; break; 87 case PM_PAR_XPOS: params_min = -100; break; 88 case PM_PAR_YPOS: params_min = -100; break; 89 case PM_PAR_SXX: params_min = 0.5; break; 90 case PM_PAR_SYY: params_min = 0.5; break; 91 case PM_PAR_SXY: params_min = -5.0; break; 92 default: 93 psAbort("invalid parameter %d for param min test", nParam); 94 } 95 if (params[nParam] < params_min) { 96 params[nParam] = params_min; 97 return false; 98 } 99 return true; 100 100 case PS_MINIMIZE_PARAM_MAX: 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 101 switch (nParam) { 102 case PM_PAR_SKY: params_max = 1e5; break; 103 case PM_PAR_I0: params_max = 1e8; break; 104 case PM_PAR_XPOS: params_max = 1e4; break; 105 case PM_PAR_YPOS: params_max = 1e4; break; 106 case PM_PAR_SXX: params_max = 100; break; 107 case PM_PAR_SYY: params_max = 100; break; 108 case PM_PAR_SXY: params_max = +5.0; break; 109 default: 110 psAbort("invalid parameter %d for param max test", nParam); 111 } 112 if (params[nParam] > params_max) { 113 params[nParam] = params_max; 114 return false; 115 } 116 return true; 117 117 default: 118 118 psAbort("invalid choice for limits"); 119 119 } 120 120 psAbort("should not reach here"); … … 204 204 } 205 205 206 bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, pmPSF *psf)206 bool PM_MODEL_FROM_PSF (pmModel *modelPSF, pmModel *modelFLT, const pmPSF *psf) 207 207 { 208 208 psF32 *out = modelPSF->params->data.F32; … … 214 214 215 215 for (int i = 0; i < psf->params->n; i++) { 216 217 218 } else { 219 220 221 216 if (psf->params->data[i] == NULL) { 217 out[i] = in[i]; 218 } else { 219 pmTrend2D *trend = psf->params->data[i]; 220 out[i] = pmTrend2DEval(trend, in[PM_PAR_XPOS], in[PM_PAR_YPOS]); 221 } 222 222 } 223 223 … … 230 230 // construct the PSF model from the FLT model and the psf 231 231 // XXX is this sufficiently general do be a global function, not a pmModelClass function? 232 bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, pmPSF *psf, float Xo, float Yo, float Io)232 bool PM_MODEL_PARAMS_FROM_PSF (pmModel *model, const pmPSF *psf, float Xo, float Yo, float Io) 233 233 { 234 234 psF32 *PAR = model->params->data.F32; … … 242 242 PAR[PM_PAR_XPOS] = Xo; 243 243 PAR[PM_PAR_YPOS] = Yo; 244 244 245 245 // supply the model-fitted parameters, or copy from the input 246 246 for (int i = 0; i < psf->params->n; i++) { 247 248 249 247 if (i == PM_PAR_SKY) continue; 248 pmTrend2D *trend = psf->params->data[i]; 249 PAR[i] = pmTrend2DEval(trend, Xo, Yo); 250 250 } 251 251 … … 254 254 // XXX user-defined value for limit? 255 255 if (!pmPSF_FitToModel (PAR, 0.1)) { 256 257 256 psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo); 257 return false; 258 258 } 259 259 … … 265 265 continue; 266 266 267 267 bool status = true; 268 268 status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MIN, i, PAR, NULL); 269 269 status &= PM_MODEL_LIMITS (PS_MINIMIZE_PARAM_MAX, i, PAR, NULL); 270 271 272 273 270 if (!status) { 271 psTrace ("psModules.objects", 5, "Hitting parameter limits at (r,c) = (%.1f, %.1f)", Xo, Yo); 272 model->flags |= PM_MODEL_STATUS_LIMITS; 273 } 274 274 } 275 275 return(true);
Note:
See TracChangeset
for help on using the changeset viewer.
