Changeset 35961 for branches/eam_branches/ipp-20130711/psModules
- Timestamp:
- Aug 15, 2013, 5:56:56 PM (13 years ago)
- Location:
- branches/eam_branches/ipp-20130711/psModules/src/objects
- Files:
-
- 10 edited
-
Makefile.am (modified) (3 diffs)
-
models/pmModel_DEV.c (modified) (7 diffs)
-
models/pmModel_EXP.c (modified) (6 diffs)
-
pmModel_CentralPixel.c (modified) (6 diffs)
-
pmPCM_MinimizeChisq.c (modified) (5 diffs)
-
pmPCMdata.c (modified) (2 diffs)
-
pmSourceFitModel.c (modified) (2 diffs)
-
pmSourceFitModel.h (modified) (1 diff)
-
pmSourceFitPCM.c (modified) (1 diff)
-
pmSourceFitSet.c (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
branches/eam_branches/ipp-20130711/psModules/src/objects/Makefile.am
r34823 r35961 20 20 pmModelClass.c \ 21 21 pmModelUtils.c \ 22 pmModel_CentralPixel.c \ 22 23 pmSource.c \ 23 24 pmPhotObj.c \ … … 97 98 pmModelClass.h \ 98 99 pmModelUtils.h \ 100 pmModel_CentralPixel.h \ 99 101 pmSource.h \ 100 102 pmPhotObj.h \ … … 110 112 pmSourceOutputs.h \ 111 113 pmSourceIO.h \ 112 pmSourceSatstar.h \ 114 pmSourceSatstar.h \ 113 115 pmSourcePlots.h \ 114 116 pmSourceVisual.h \ -
branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_DEV.c
r35876 r35961 16 16 * PM_PAR_SYY 5 - Y^2 term of elliptical contour (sqrt(2) / SigmaY) 17 17 * PM_PAR_SXY 6 - X*Y term of elliptical contour 18 * PM_PAR_7 7 - normalized dev parameter19 18 20 19 note that a standard dev model uses exp(-K*(z^(1/2n) - 1). the additional elements (K, … … 86 85 static float *paramsMinUse = paramsMinLax; 87 86 static float *paramsMaxUse = paramsMaxLax; 88 static float betaUse[] = { 1000, 3e6, 5, 5, 1.0, 1.0, 0.5 };87 static float betaUse[] = { 2, 3e6, 5, 5, 3.0, 3.0, 0.5 }; 89 88 90 89 static bool limitsApply = true; // Apply limits? 91 90 92 # include "pmModel_SERSIC.CP.h"91 // # include "pmModel_SERSIC.CP.h" 93 92 94 93 psF32 PM_MODEL_FUNC (psVector *deriv, … … 109 108 psAssert (z >= 0, "do not allow negative z values in model"); 110 109 111 float index = 0.5 / ALPHA; 112 float par7 = ALPHA; 113 float bn = 1.9992*index - 0.3271; 114 float Io = exp(bn); 115 116 psF32 f2 = bn*pow(z,ALPHA); 117 psF32 f1 = Io*exp(-f2); 118 110 // for DEV, we can hard-wire kappa(4): 111 // float index = 4.0; 112 float kappa = 7.670628; 113 114 // r = sqrt(z) 115 float q = kappa*pow(z,ALPHA); 116 psF32 f0 = exp(-q); 117 118 psF32 f1 = PAR[PM_PAR_I0]*f0; 119 psF32 f = PAR[PM_PAR_SKY] + f1; 120 121 assert (isfinite(q)); 122 assert (isfinite(f0)); 123 assert (isfinite(f1)); 124 assert (isfinite(f)); 125 126 // only worry about the central 4 pixels at most 127 // If I use DELTA = 0.2, I'm way off for the total flux 128 // If I use DELTA = 0.02, I'm totally good (but I am under on the total flux for R = 30 by 0.2 mags -- aperture failure) 129 // For DELTA = 0.02 & Rmin/Rmaj = 0.25, I'm over flux by 0.15 mags (due to the central pixel) 119 130 psF32 radius = hypot(X, Y); 120 131 if (radius < 1.0) { 121 122 // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center 123 124 // first, use Rmajor and index to find the central pixel flux (fraction of total flux) 125 psEllipseAxes axes; 126 pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); 127 128 // get the central pixel flux from the lookup table 129 float xPix = (axes.major - centralPixelXo) / centralPixeldX; 130 xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1); 131 float yPix = (index - centralPixelYo) / centralPixeldY; 132 yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1); 133 134 // the integral of a Sersic has an analytical form as follows: 135 float logGamma = lgamma(2.0*index); 136 float bnFactor = pow(bn, 2.0*index); 137 float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; 138 139 // XXX interpolate to get the value 140 // XXX for the moment, just integerize 141 // XXX I need to multiply by the integrated flux to get the flux in the central pixel 142 float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm; 143 144 float px1 = 1.0 / PAR[PM_PAR_SXX]; 145 float py1 = 1.0 / PAR[PM_PAR_SYY]; 146 float z10 = PS_SQR(px1); 147 float z01 = PS_SQR(py1); 148 149 // which pixels do we need for this interpolation? 150 // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...) 151 if ((X >= 0) && (Y >= 0)) { 152 float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive 153 float V00 = Vcenter; 154 float V10 = Io*exp(-bn*pow(z10,par7)); 155 float V01 = Io*exp(-bn*pow(z01,par7)); 156 float V11 = Io*exp(-bn*pow(z11,par7)); 157 f1 = interpolatePixels(V00, V10, V01, V11, X, Y); 132 // subdivide the central 2,3,4 pixels by Nx,Ny 133 float Npix = 0.0; 134 float Fpix = 0.0; 135 float Xpix = floor(pixcoord->data.F32[0]) - PAR[PM_PAR_XPOS]; 136 float Ypix = floor(pixcoord->data.F32[1]) - PAR[PM_PAR_YPOS]; 137 # define DELTA 0.02 138 for (float ix = 0.1; ix <= 0.9; ix += DELTA) { 139 for (float iy = 0.1; iy <= 0.9; iy += DELTA) { 140 psF32 X = Xpix + ix; 141 psF32 Y = Ypix + iy; 142 psF32 px = X / PAR[PM_PAR_SXX]; 143 psF32 py = Y / PAR[PM_PAR_SYY]; 144 psF32 z = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; 145 146 // sqrt(z) is r 147 float q = kappa*pow(z,ALPHA); 148 psF32 f0 = exp(-q); 149 150 psF32 f1 = PAR[PM_PAR_I0]*f0; 151 psF32 fx = PAR[PM_PAR_SKY] + f1; 152 Fpix += fx; 153 Npix += 1.0; 158 154 } 159 if ((X < 0) && (Y >= 0)) { 160 float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative 161 float V00 = Io*exp(-bn*pow(z10,par7)); 162 float V10 = Vcenter; 163 float V01 = Io*exp(-bn*pow(z11,par7)); 164 float V11 = Io*exp(-bn*pow(z01,par7)); 165 f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y); 166 } 167 if ((X >= 0) && (Y < 0)) { 168 float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative 169 float V00 = Io*exp(-bn*pow(z01,par7)); 170 float V10 = Io*exp(-bn*pow(z11,par7)); 171 float V01 = Vcenter; 172 float V11 = Io*exp(-bn*pow(z10,par7)); 173 f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y)); 174 } 175 if ((X < 0) && (Y < 0)) { 176 float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive 177 float V00 = Io*exp(-bn*pow(z11,par7)); 178 float V10 = Io*exp(-bn*pow(z10,par7)); 179 float V01 = Io*exp(-bn*pow(z01,par7)); 180 float V11 = Vcenter; 181 f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y)); 182 } 155 } 156 f = Fpix / Npix; 183 157 } 184 185 psF32 z0 = PAR[PM_PAR_I0]*f1;186 psF32 f0 = PAR[PM_PAR_SKY] + z0;187 188 assert (isfinite(f2));189 assert (isfinite(f1));190 assert (isfinite(z0));191 assert (isfinite(f0));192 158 193 159 if (deriv != NULL) { … … 195 161 196 162 dPAR[PM_PAR_SKY] = +1.0; 197 dPAR[PM_PAR_I0] = +2.0*f1; // XXX extra damping.. 198 199 // gradient is infinite for z = 0; saturate at z = 0.01 200 psF32 z1 = (z < 0.01) ? z0*bn*ALPHA*pow(0.01,ALPHA - 1.0) : z0*bn*ALPHA*pow(z,ALPHA - 1.0); 201 202 assert (isfinite(z1)); 203 204 dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); 205 dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); 206 dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; 207 dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; 208 dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; 209 } 210 return (f0); 163 dPAR[PM_PAR_I0] = +f0; 164 165 if (z > 0.01) { 166 float z1 = f1*kappa*ALPHA*pow(z,ALPHA-1.0); 167 dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]); 168 dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]); 169 dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; 170 dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; 171 dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; 172 } else { 173 // gradient -> 0 for z -> 0, but has undef form 174 float z1 = f1*kappa*ALPHA*pow(z,ALPHA); 175 dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]); 176 dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]); 177 dPAR[PM_PAR_SXX] = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX]; 178 dPAR[PM_PAR_SYY] = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY]; 179 dPAR[PM_PAR_SXY] = -1.0*z1; 180 } 181 } 182 return (f); 211 183 } 212 184 … … 302 274 } 303 275 304 // the normalization is modified by the slope305 float index = 0.5 / ALPHA;306 float bn = 1.9992*index - 0.3271;307 float Io = exp(0.5*bn);308 309 276 // set the model normalization 310 277 if (!pmModelSetNorm(&PAR[PM_PAR_I0], source)) { 311 278 return false; 312 279 } 313 PAR[PM_PAR_I0] /= Io;314 280 315 281 // set the model position … … 328 294 psEllipseAxes axes; 329 295 pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); 330 float AspectRatio = axes.minor / axes.major; 331 332 float index = 4.0; 333 float bn = 1.9992*index - 0.3271; 334 335 // the integral of a Sersic has an analytical form as follows: 336 float logGamma = lgamma(2.0*index); 337 float bnFactor = pow(bn, 2.0*index); 338 float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; 339 340 psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio; 341 342 return(Flux); 296 297 float norm = 0.00168012; 298 float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm; 299 300 return(flux); 343 301 } 344 302 … … 359 317 pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); 360 318 361 // f = Io exp(-z^n) -> z^n = ln(Io/f) 362 psF64 zn = log(PAR[PM_PAR_I0] / flux); 319 // static value for DEV: 320 float kappa = 7.670628; 321 322 // f = Io exp(-kappa*z^n) -> z^n = ln(Io/f) / kappa 323 psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa; 363 324 psF64 radius = axes.major * sqrt (2.0) * pow(zn, 0.5 / ALPHA); 364 325 -
branches/eam_branches/ipp-20130711/psModules/src/objects/models/pmModel_EXP.c
r35876 r35961 82 82 static bool limitsApply = true; // Apply limits? 83 83 84 # include "pmModel_SERSIC.CP.h" 84 // # include "pmModel_SERSIC.CP.h" 85 86 // the problems I'm having with the SERSIC-like functions are: 87 // 1) making sure I have the right functional form so that PAR[SXX,etc] represent R_eff (half-light radius) 88 // 2) getting the central pixel right 89 // 3) getting the derivaties right. 85 90 86 91 psF32 PM_MODEL_FUNC (psVector *deriv, … … 101 106 psAssert (z >= 0, "do not allow negative z values in model"); 102 107 103 float index = 1.0; 104 float par7 = 0.5; 105 float bn = 1.9992*index - 0.3271; 106 float Io = exp(bn); 107 108 psF32 f2 = bn*sqrt(z); 109 psF32 f1 = Io*exp(-f2); 110 108 // for EXP, we can hard-wire kappa(1): 109 // float index = 1.0; 110 float kappa = 1.70056; 111 112 // sqrt(z) is r 113 float q = kappa*sqrt(z); 114 psF32 f0 = exp(-q); 115 116 psF32 f1 = PAR[PM_PAR_I0]*f0; 117 psF32 f = PAR[PM_PAR_SKY] + f1; 118 119 assert (isfinite(q)); 120 assert (isfinite(f0)); 121 assert (isfinite(f1)); 122 assert (isfinite(f)); 123 124 // only worry about the central 4 pixels at most 111 125 psF32 radius = hypot(X, Y); 112 126 if (radius < 1.0) { 113 114 // ** use bilinear interpolation to the given location from the 4 surrounding pixels centered on the object center 115 116 // first, use Rmajor and index to find the central pixel flux (fraction of total flux) 117 psEllipseAxes axes; 118 pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); 119 120 // get the central pixel flux from the lookup table 121 float xPix = (axes.major - centralPixelXo) / centralPixeldX; 122 xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1); 123 float yPix = (index - centralPixelYo) / centralPixeldY; 124 yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1); 125 126 // the integral of a Sersic has an analytical form as follows: 127 float logGamma = lgamma(2.0*index); 128 float bnFactor = pow(bn, 2.0*index); 129 float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; 130 131 // XXX interpolate to get the value 132 // XXX for the moment, just integerize 133 // XXX I need to multiply by the integrated flux to get the flux in the central pixel 134 float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm; 135 136 float px1 = 1.0 / PAR[PM_PAR_SXX]; 137 float py1 = 1.0 / PAR[PM_PAR_SYY]; 138 float z10 = PS_SQR(px1); 139 float z01 = PS_SQR(py1); 140 141 // which pixels do we need for this interpolation? 142 // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...) 143 if ((X >= 0) && (Y >= 0)) { 144 float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive 145 float V00 = Vcenter; 146 float V10 = Io*exp(-bn*pow(z10,par7)); 147 float V01 = Io*exp(-bn*pow(z01,par7)); 148 float V11 = Io*exp(-bn*pow(z11,par7)); 149 f1 = interpolatePixels(V00, V10, V01, V11, X, Y); 127 // subdivide the central 2,3,4 pixels by Nx,Ny 128 float Npix = 0.0; 129 float Fpix = 0.0; 130 float Xpix = floor(pixcoord->data.F32[0]) - PAR[PM_PAR_XPOS]; 131 float Ypix = floor(pixcoord->data.F32[1]) - PAR[PM_PAR_YPOS]; 132 for (float ix = 0.1; ix < 1.0; ix += 0.2) { 133 for (float iy = 0.1; iy < 1.0; iy += 0.2) { 134 psF32 X = Xpix + ix; 135 psF32 Y = Ypix + iy; 136 psF32 px = X / PAR[PM_PAR_SXX]; 137 psF32 py = Y / PAR[PM_PAR_SYY]; 138 psF32 z = PS_SQR(px) + PS_SQR(py) + PAR[PM_PAR_SXY]*X*Y; 139 140 // sqrt(z) is r 141 float q = kappa*sqrt(z); 142 psF32 f0 = exp(-q); 143 144 psF32 f1 = PAR[PM_PAR_I0]*f0; 145 psF32 fx = PAR[PM_PAR_SKY] + f1; 146 Fpix += fx; 147 Npix += 1.0; 150 148 } 151 if ((X < 0) && (Y >= 0)) { 152 float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative 153 float V00 = Io*exp(-bn*pow(z10,par7)); 154 float V10 = Vcenter; 155 float V01 = Io*exp(-bn*pow(z11,par7)); 156 float V11 = Io*exp(-bn*pow(z01,par7)); 157 f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y); 158 } 159 if ((X >= 0) && (Y < 0)) { 160 float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative 161 float V00 = Io*exp(-bn*pow(z01,par7)); 162 float V10 = Io*exp(-bn*pow(z11,par7)); 163 float V01 = Vcenter; 164 float V11 = Io*exp(-bn*pow(z10,par7)); 165 f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y)); 166 } 167 if ((X < 0) && (Y < 0)) { 168 float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive 169 float V00 = Io*exp(-bn*pow(z11,par7)); 170 float V10 = Io*exp(-bn*pow(z10,par7)); 171 float V01 = Io*exp(-bn*pow(z01,par7)); 172 float V11 = Vcenter; 173 f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y)); 174 } 175 } 176 177 psF32 z0 = PAR[PM_PAR_I0]*f1; 178 psF32 f0 = PAR[PM_PAR_SKY] + z0; 179 180 assert (isfinite(f2)); 181 assert (isfinite(f1)); 182 assert (isfinite(z0)); 183 assert (isfinite(f0)); 149 } 150 f = Fpix / Npix; 151 } 184 152 185 153 if (deriv != NULL) { … … 187 155 188 156 dPAR[PM_PAR_SKY] = +1.0; 189 dPAR[PM_PAR_I0] = +f1; 190 191 // gradient is infinite for z = 0; saturate at z = 0.01 192 // z1 is -df/dz (the negative sign is canceled by most of dz/dPAR[i] 193 psF32 z1 = (z < 0.01) ? 0.5*bn*z0/sqrt(0.01) : 0.5*bn*z0/sqrt(z); 194 195 // XXX dampen SXX and SYY as in GAUSS? 196 dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px/PAR[PM_PAR_SXX] + Y*PAR[PM_PAR_SXY]); 197 dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py/PAR[PM_PAR_SYY] + X*PAR[PM_PAR_SXY]); 198 dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; 199 dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; 200 dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; 201 } 202 return (f0); 157 dPAR[PM_PAR_I0] = +f0; 158 159 if (z > 0.01) { 160 float z1 = 0.5*f1*kappa/sqrt(z); 161 dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0*px + Y*PAR[PM_PAR_SXY]); 162 dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0*py + X*PAR[PM_PAR_SXY]); 163 dPAR[PM_PAR_SXX] = +2.0*z1*px*px/PAR[PM_PAR_SXX]; 164 dPAR[PM_PAR_SYY] = +2.0*z1*py*py/PAR[PM_PAR_SYY]; 165 dPAR[PM_PAR_SXY] = -1.0*z1*X*Y; 166 } else { 167 // gradient -> 0 for z -> 0, but has undef form 168 float z1 = 0.5*f1*kappa; 169 dPAR[PM_PAR_XPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SXX] + PAR[PM_PAR_SXY]); 170 dPAR[PM_PAR_YPOS] = +1.0*z1*(2.0/PAR[PM_PAR_SYY] + PAR[PM_PAR_SXY]); 171 dPAR[PM_PAR_SXX] = +2.0*z1*px/PAR[PM_PAR_SXX]/PAR[PM_PAR_SXX]; 172 dPAR[PM_PAR_SYY] = +2.0*z1*py/PAR[PM_PAR_SYY]/PAR[PM_PAR_SYY]; 173 dPAR[PM_PAR_SXY] = -1.0*z1; 174 } 175 } 176 return (f); 203 177 } 204 178 … … 314 288 psEllipseAxes axes; 315 289 pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); 316 float AspectRatio = axes.minor / axes.major; 317 318 float index = 1.0; 319 float bn = 1.9992*index - 0.3271; 320 321 // the integral of a Sersic has an analytical form as follows: 322 float logGamma = lgamma(2.0*index); 323 float bnFactor = pow(bn, 2.0*index); 324 float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; 325 326 psF64 Flux = PAR[PM_PAR_I0] * norm * AspectRatio; 327 328 return(Flux); 290 291 // static value for EXP: 292 float norm = 0.34578; // \int exp(-kappa*sqrt(z)) r dr 293 294 float flux = PAR[PM_PAR_I0] * 2.0 * M_PI * axes.major * axes.minor * norm; 295 296 return(flux); 329 297 } 330 298 … … 345 313 pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); 346 314 347 // f = Io exp(-sqrt(z)) -> sqrt(z) = ln(Io/f) 348 psF64 zn = log(PAR[PM_PAR_I0] / flux); 315 // static value for EXP: 316 float kappa = 1.70056; 317 318 // f = Io exp(-kappa*sqrt(z)) -> sqrt(z) = ln(Io/f) / kappa 319 psF64 zn = log(PAR[PM_PAR_I0] / flux) / kappa; 349 320 psF64 radius = axes.major * sqrt (2.0) * zn; 350 321 … … 501 472 return; 502 473 } 474 475 # if (0) 476 void bilin_inter_function () { 477 // first, use Rmajor and index to find the central pixel flux (fraction of total flux) 478 psEllipseAxes axes; 479 pmModelParamsToAxes (&axes, PAR[PM_PAR_SXX], PAR[PM_PAR_SXY], PAR[PM_PAR_SYY], true); 480 481 // get the central pixel flux from the lookup table 482 float xPix = (axes.major - centralPixelXo) / centralPixeldX; 483 xPix = PS_MIN (PS_MAX(xPix, 0), centralPixelNX - 1); 484 float yPix = (index - centralPixelYo) / centralPixeldY; 485 yPix = PS_MIN (PS_MAX(yPix, 0), centralPixelNY - 1); 486 487 // the integral of a Sersic has an analytical form as follows: 488 float logGamma = lgamma(2.0*index); 489 float bnFactor = pow(bn, 2.0*index); 490 float norm = 2.0 * M_PI * PS_SQR(axes.major) * index * exp(bn) * exp(logGamma) / bnFactor; 491 492 // XXX interpolate to get the value 493 // XXX for the moment, just integerize 494 // XXX I need to multiply by the integrated flux to get the flux in the central pixel 495 float Vcenter = centralPixel[(int)yPix][(int)xPix] * norm; 496 497 float px1 = 1.0 / PAR[PM_PAR_SXX]; 498 float py1 = 1.0 / PAR[PM_PAR_SYY]; 499 float z10 = PS_SQR(px1); 500 float z01 = PS_SQR(py1); 501 502 // which pixels do we need for this interpolation? 503 // (I do not keep state information, so I don't know anything about other evaluations of nearby pixels...) 504 if ((X >= 0) && (Y >= 0)) { 505 float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive 506 float V00 = Vcenter; 507 float V10 = Io*exp(-bn*pow(z10,par7)); 508 float V01 = Io*exp(-bn*pow(z01,par7)); 509 float V11 = Io*exp(-bn*pow(z11,par7)); 510 f1 = interpolatePixels(V00, V10, V01, V11, X, Y); 511 } 512 if ((X < 0) && (Y >= 0)) { 513 float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative 514 float V00 = Io*exp(-bn*pow(z10,par7)); 515 float V10 = Vcenter; 516 float V01 = Io*exp(-bn*pow(z11,par7)); 517 float V11 = Io*exp(-bn*pow(z01,par7)); 518 f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), Y); 519 } 520 if ((X >= 0) && (Y < 0)) { 521 float z11 = z10 + z01 - PAR[PM_PAR_SXY]; // X * Y negative 522 float V00 = Io*exp(-bn*pow(z01,par7)); 523 float V10 = Io*exp(-bn*pow(z11,par7)); 524 float V01 = Vcenter; 525 float V11 = Io*exp(-bn*pow(z10,par7)); 526 f1 = interpolatePixels(V00, V10, V01, V11, X, (1.0 + Y)); 527 } 528 if ((X < 0) && (Y < 0)) { 529 float z11 = z10 + z01 + PAR[PM_PAR_SXY]; // X * Y positive 530 float V00 = Io*exp(-bn*pow(z11,par7)); 531 float V10 = Io*exp(-bn*pow(z10,par7)); 532 float V01 = Io*exp(-bn*pow(z01,par7)); 533 float V11 = Vcenter; 534 f1 = interpolatePixels(V00, V10, V01, V11, (1.0 + X), (1.0 + Y)); 535 } 536 } 537 # endif -
branches/eam_branches/ipp-20130711/psModules/src/objects/pmModel_CentralPixel.c
r35948 r35961 209 209 210 210 // XXX for test purposes only: 211 # define TEST_IMAGE 1211 # define TEST_IMAGE 0 212 212 # if (TEST_IMAGE) 213 213 static psImage *map = NULL; 214 214 # endif 215 215 216 float pmModelCP_GetFlux_RotSquare (pmModelCP *cp, float dx, float dy, float theta); 217 216 218 float pmModelCP_GetFlux (pmModelCP *cp, float dx, float dy, float theta) { 217 219 … … 221 223 # endif 222 224 223 float flux = pmModelCP_GetFlux_Bresen (cp, dx, dy, theta); 224 225 // float flux = pmModelCP_GetFlux_Bresen (cp, dx, dy, theta); 226 // float flux = pmModelCP_GetFlux_Old (cp, dx, dy, theta); 227 float flux = pmModelCP_GetFlux_RotSquare (cp, dx, dy, theta); 228 229 // RotSquare for theta = 0.0 & Bresen give the same answer 230 // if I count from x[0] <= ix < x[1] 231 225 232 # if (TEST_IMAGE) 226 233 psFits *fits = psFitsOpen ("map.fits", "w"); … … 309 316 } 310 317 311 float pmModelCP_GetFlux_RotSqaure (pmModelCP *cp, float dx, float dy, float theta) { 318 // *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors 319 # define COMPARE_INDEX(A,B) (y[A] < y[B]) 320 # define SWAP_INDEX(TYPE,A,B) { \ 321 int tmp; \ 322 if (A != B) { \ 323 tmp = x[A]; \ 324 x[A] = x[B]; \ 325 x[B] = tmp; \ 326 tmp = y[A]; \ 327 y[A] = y[B]; \ 328 y[B] = tmp; \ 329 } \ 330 } 331 332 bool pmModelCP_SortCorners (int *x, int *y, int Npar) { 333 334 if (Npar < 2) return true; 335 336 // sort the vector set by the radius 337 PSSORT (Npar, COMPARE_INDEX, SWAP_INDEX, NONE); 338 return true; 339 } 340 341 float pmModelCP_GetFlux_RotSquare (pmModelCP *cp, float dx, float dy, float theta) { 312 342 313 343 // the cp data is defined for the central 3x3 pixels. we allow dx,dy to have values of … … 324 354 float cs = cos(theta*PS_RAD_DEG); 325 355 float sn = sin(theta*PS_RAD_DEG); 326 327 356 float Nsub = 11.0; 328 int Xsub00 = ((dx - 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub; 329 int Ysub00 = ((dx - 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub; 330 int Xsub01 = ((dx - 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub; 331 int Ysub01 = ((dx - 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub; 332 int Xsub10 = ((dx + 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub; 333 int Ysub10 = ((dx + 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub; 334 int Xsub11 = ((dx + 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub; 335 int Ysub11 = ((dx + 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub; 336 337 /* generic rotated square: 338 339 */ 340 341 int Xmin, Xmax, Ymin, Ymax; 342 343 Xmin = PS_MIN(Xsub00,Xsub01); 344 Xmin = PS_MIN(Xsub10,Xmin); 345 Xmin = PS_MIN(Xsub11,Xmin); 346 Xmin = PS_MIN(Xmin, cp->flux->numCols - 1); 347 Xmin = PS_MAX(Xmin, 0); 348 Xmax = PS_MAX(Xsub00,Xsub01); 349 Xmax = PS_MAX(Xsub10,Xmax); 350 Xmax = PS_MAX(Xsub11,Xmax); 351 Xmax = PS_MIN(Xmax, cp->flux->numCols - 1); 352 Xmax = PS_MAX(Xmax, 0); 353 Ymin = PS_MIN(Ysub00,Ysub01); 354 Ymin = PS_MIN(Ysub10,Ymin); 355 Ymin = PS_MIN(Ysub11,Ymin); 356 Ymin = PS_MIN(Ymin, cp->flux->numRows - 1); 357 Ymin = PS_MAX(Ymin, 0); 358 Ymax = PS_MAX(Ysub00,Ysub01); 359 Ymax = PS_MAX(Ysub10,Ymax); 360 Ymax = PS_MAX(Ysub11,Ymax); 361 Ymax = PS_MIN(Ymax, cp->flux->numRows - 1); 362 Ymax = PS_MAX(Ymax, 0); 363 364 // integrate pixels from Xmin,Ymin to Xmax,Ymax, only include pixels contained in the 365 // target pixel 357 358 int Xsub[4], Ysub[4]; 359 360 Xsub[0] = ((dx - 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub; 361 Ysub[0] = ((dx - 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub; 362 Xsub[1] = ((dx - 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub; 363 Ysub[1] = ((dx - 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub; 364 Xsub[2] = ((dx + 0.5)*cs - (dy - 0.5)*sn + 1.5)*Nsub; 365 Ysub[2] = ((dx + 0.5)*sn + (dy - 0.5)*cs + 1.5)*Nsub; 366 Xsub[3] = ((dx + 0.5)*cs - (dy + 0.5)*sn + 1.5)*Nsub; 367 Ysub[3] = ((dx + 0.5)*sn + (dy + 0.5)*cs + 1.5)*Nsub; 368 369 // first, sort the corners in order of the Y coordinate: 370 pmModelCP_SortCorners (Xsub, Ysub, 4); 366 371 367 372 float flux = 0.0; 368 int npix = 0; 369 for (int i = Xmin; i < Xmax; i++) { 370 float dX = i / Nsub - 1.5; 371 for (int j = Ymin; j < Ymax; j++) { 372 float dY = j / Nsub - 1.5; 373 374 float Xim = dX*cs + dY*sn; 375 if (Xim < (dx - 0.5)) continue; 376 if (Xim > (dx + 0.5)) continue; 377 378 float Yim = -dX*sn + dY*cs; 379 if (Yim < (dy - 0.5)) continue; 380 if (Yim > (dy + 0.5)) continue; 381 382 flux += cp->flux->data.F32[j][i]; 383 npix ++; 384 } 385 } 386 387 float normFlux = flux / npix; 388 return normFlux; 373 float npix = 0.0; 374 375 // if (Ysub[0] == Ysub[1]), we have a simple square 376 if (Ysub[0] == Ysub[1]) { 377 psAssert (Ysub[2] == Ysub[3], "not square?"); 378 int Xmin = PS_MIN(Xsub[0], Xsub[1]); 379 int Xmax = PS_MAX(Xsub[0], Xsub[1]); 380 for (int iy = Ysub[0]; iy < Ysub[3]; iy++) { 381 for (int ix = Xmin; ix < Xmax; ix++) { 382 flux += cp->flux->data.F32[iy][ix]; 383 npix += 1.0; 384 # if (TEST_IMAGE) 385 fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]); 386 map->data.S32[iy][ix] ++; 387 # endif 388 } 389 } 390 float normFlux = flux / npix; 391 return normFlux; 392 } 393 394 // second case: Xsub[1] > Xsub[2]: 395 if (Xsub[1] > Xsub[2]) { 396 float dYdXp, dYdXm; 397 // first segment, Ysub[0] to Ysub[1]: 398 dYdXp = (Ysub[1] - Ysub[0]) / (float) (Xsub[1] - Xsub[0]); 399 dYdXm = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]); 400 for (int iy = Ysub[0]; iy < Ysub[1]; iy++) { 401 int Xs = (iy - Ysub[0]) / dYdXm + Xsub[0]; 402 int Xe = (iy - Ysub[0]) / dYdXp + Xsub[0]; 403 for (int ix = Xs; ix < Xe; ix ++) { 404 flux += cp->flux->data.F32[iy][ix]; 405 npix += 1.0; 406 # if (TEST_IMAGE) 407 fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]); 408 map->data.S32[iy][ix] ++; 409 # endif 410 } 411 } 412 // 2nd segment, Ysub[1] to Ysub[2]: 413 dYdXp = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]); 414 dYdXm = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]); 415 for (int iy = Ysub[1]; iy < Ysub[2]; iy++) { 416 int Xs = (iy - Ysub[0]) / dYdXm + Xsub[0]; 417 int Xe = (iy - Ysub[1]) / dYdXp + Xsub[1]; 418 for (int ix = Xs; ix < Xe; ix ++) { 419 flux += cp->flux->data.F32[iy][ix]; 420 npix += 1.0; 421 # if (TEST_IMAGE) 422 fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]); 423 map->data.S32[iy][ix] ++; 424 # endif 425 } 426 } 427 // first segment, Ysub[0] to Ysub[1]: 428 dYdXp = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]); 429 dYdXm = (Ysub[3] - Ysub[2]) / (float) (Xsub[3] - Xsub[2]); 430 for (int iy = Ysub[2]; iy < Ysub[3]; iy++) { 431 int Xs = (iy - Ysub[2]) / dYdXm + Xsub[2]; 432 int Xe = (iy - Ysub[1]) / dYdXp + Xsub[1]; 433 for (int ix = Xs; ix < Xe; ix ++) { 434 flux += cp->flux->data.F32[iy][ix]; 435 npix += 1.0; 436 # if (TEST_IMAGE) 437 fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]); 438 map->data.S32[iy][ix] ++; 439 # endif 440 } 441 } 442 float normFlux = flux / npix; 443 return normFlux; 444 } 445 446 // third case: Xsub[1] < Xsub[2]: 447 if (Xsub[2] > Xsub[1]) { 448 // first segment, Ysub[0] to Ysub[1]: 449 float dYdXp, dYdXm; 450 dYdXp = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]); 451 dYdXm = (Ysub[1] - Ysub[0]) / (float) (Xsub[1] - Xsub[0]); 452 for (int iy = Ysub[0]; iy < Ysub[1]; iy++) { 453 int Xs = (iy - Ysub[0]) / dYdXm + Xsub[0]; 454 int Xe = (iy - Ysub[0]) / dYdXp + Xsub[0]; 455 for (int ix = Xs; ix < Xe; ix ++) { 456 flux += cp->flux->data.F32[iy][ix]; 457 npix += 1.0; 458 # if (TEST_IMAGE) 459 fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]); 460 map->data.S32[iy][ix] ++; 461 # endif 462 } 463 } 464 // 2nd segment, Ysub[1] to Ysub[2]: 465 dYdXp = (Ysub[2] - Ysub[0]) / (float) (Xsub[2] - Xsub[0]); 466 dYdXm = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]); 467 for (int iy = Ysub[1]; iy < Ysub[2]; iy++) { 468 int Xs = (iy - Ysub[1]) / dYdXm + Xsub[1]; 469 int Xe = (iy - Ysub[0]) / dYdXp + Xsub[0]; 470 for (int ix = Xs; ix < Xe; ix ++) { 471 flux += cp->flux->data.F32[iy][ix]; 472 npix += 1.0; 473 # if (TEST_IMAGE) 474 fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]); 475 map->data.S32[iy][ix] ++; 476 # endif 477 } 478 } 479 // first segment, Ysub[0] to Ysub[1]: 480 dYdXp = (Ysub[3] - Ysub[2]) / (float) (Xsub[3] - Xsub[2]); 481 dYdXm = (Ysub[3] - Ysub[1]) / (float) (Xsub[3] - Xsub[1]); 482 for (int iy = Ysub[2]; iy < Ysub[3]; iy++) { 483 int Xs = (iy - Ysub[1]) / dYdXm + Xsub[1]; 484 int Xe = (iy - Ysub[2]) / dYdXp + Xsub[2]; 485 for (int ix = Xs; ix < Xe; ix ++) { 486 flux += cp->flux->data.F32[iy][ix]; 487 npix += 1.0; 488 # if (TEST_IMAGE) 489 fprintf (stderr, "%d %d | %f %f | %f\n", ix, iy, flux, npix, cp->flux->data.F32[iy][ix]); 490 map->data.S32[iy][ix] ++; 491 # endif 492 } 493 } 494 float normFlux = flux / npix; 495 return normFlux; 496 } 497 myAbort ("impossible case?"); 389 498 } 390 499 … … 478 587 } 479 588 float normFlux = flux / npix; 480 fprintf (stderr, "bres: %f %f %f\n", flux, (float) npix, normFlux);589 // fprintf (stderr, "bres: %f %f %f\n", flux, (float) npix, normFlux); 481 590 return normFlux; 482 591 } … … 574 683 } 575 684 float normFlux = flux / npix; 576 fprintf (stderr, "full : %f %f %f\n", flux, (float) npix, normFlux);685 // fprintf (stderr, "full : %f %f %f\n", flux, (float) npix, normFlux); 577 686 return normFlux; 578 687 } -
branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCM_MinimizeChisq.c
r35768 r35961 42 42 #include "pmPCMdata.h" 43 43 44 # define SAVE_IMAGES 0 45 # if (SAVE_IMAGES) 46 int psphotSaveImage (psMetadata *header, psImage *image, char *filename); 47 # endif 48 44 49 # define FACILITY "psModules.objects" 45 50 … … 130 135 } 131 136 137 char key[10]; // used for interactive responses 138 bool testValue = false; 139 132 140 // set a new guess for Alpha, Beta, Params 133 141 if (!psMinLM_GuessABP(Alpha, Beta, Params, alpha, beta, params, paramMask, checkLimits, lambda, &dLinear)) { 142 if (min->isInteractive) { 143 fprintf (stdout, "guess failed (singular matrix or NaN values), continue? [Y,n] "); 144 if (!fgets(key, 8, stdin)) { 145 psWarning("Unable to read option"); 146 } 147 switch (key[0]) { 148 case 'n': 149 case 'N': 150 done = true; 151 break; 152 case 'y': 153 case 'Y': 154 case '\n': 155 lambda *= 10.0; 156 continue; 157 default: 158 lambda *= 10.0; 159 continue; 160 } 161 if (done) break; 162 } 134 163 min->iter ++; 135 164 if (min->iter >= min->maxIter) break; … … 138 167 } 139 168 169 if (min->isInteractive) { 170 p_psVectorPrint(psTraceGetDestination(), Params, "current parameters: "); 171 fprintf (stdout, "last chisq : %f\n", min->value); 172 bool getOptions = true; 173 while (getOptions) { 174 fprintf (stdout, "options: (m)odify, (g)o, (q)uit: "); 175 if (!fgets(key, 8, stdin)) { 176 psWarning("Unable to read option"); 177 } 178 switch (key[0]) { 179 case 'm': 180 case 'M': 181 testValue = TRUE; 182 fprintf (stdout, "enter (Npar) (value): "); 183 int Npar = 0; 184 float value= 0; 185 fscanf (stdin, "%d %f", &Npar, &value); 186 Params->data.F32[Npar] = value; 187 break; 188 case 'g': 189 case 'G': 190 case '\n': 191 getOptions = false; 192 break; 193 default: 194 done = true; 195 break; 196 } 197 fprintf (stderr, "foo\n"); 198 } 199 if (done) break; 200 } 201 140 202 // dump some useful info if trace is defined 141 203 if (psTraceGetLevel(FACILITY) >= 6) { … … 202 264 // XXX : Madsen gives suggestion for better use of rho 203 265 // rho is positive if the new chisq is smaller 204 if ( rho >= -1e-6) {266 if (testValue || (rho >= -1e-6)) { 205 267 min->value = Chisq; 206 268 alpha = psImageCopy(alpha, Alpha, PS_TYPE_F32); … … 474 536 // XXX TEST : SAVE IMAGES 475 537 # if (SAVE_IMAGES) 476 psphotSaveImage (NULL, pcm->psf->image, "psf.fits"); 477 psphotSaveImage (NULL, pcm->modelFlux, "model.fits"); 478 psphotSaveImage (NULL, pcm->modelConvFlux, "modelConv.fits"); 479 psphotSaveImage (NULL, source->pixels, "obj.fits"); 480 psphotSaveImage (NULL, source->maskObj, "mask.fits"); 481 psphotSaveImage (NULL, source->variance, "variance.fits"); 538 static int Npass = 0; 539 char name[128]; 540 snprintf (name, 128, "psf.%03d.fits", Npass); psphotSaveImage (NULL, pcm->psf->image, name); 541 snprintf (name, 128, "mod.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelFlux, name); 542 snprintf (name, 128, "cnv.%03d.fits", Npass); psphotSaveImage (NULL, pcm->modelConvFlux, name); 543 snprintf (name, 128, "obj.%03d.fits", Npass); psphotSaveImage (NULL, source->pixels, name); 544 snprintf (name, 128, "msk.%03d.fits", Npass); psphotSaveImage (NULL, source->maskObj, name); 545 snprintf (name, 128, "var.%03d.fits", Npass); psphotSaveImage (NULL, source->variance, name); 546 for (int n = 0; n < pcm->dmodelsFlux->n; n++) { 547 psImage *dmodelConv = pcm->dmodelsConvFlux->data[n]; 548 if (!dmodelConv) continue; 549 snprintf (name, 128, "dpar.%01d.%03d.fits", n, Npass); psphotSaveImage (NULL, dmodelConv, name); 550 } 551 Npass ++; 482 552 # endif 483 553 -
branches/eam_branches/ipp-20130711/psModules/src/objects/pmPCMdata.c
r35768 r35961 242 242 constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1; 243 243 break; 244 case PM_SOURCE_FIT_EXT_AND_SKY: 245 // EXT model fits all params (including sky) 246 nParams = params->n; 247 psVectorInit (constraint->paramMask, 0); 248 break; 244 249 case PM_SOURCE_FIT_INDEX: 245 250 // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params … … 365 370 pcm->constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1; 366 371 break; 372 case PM_SOURCE_FIT_EXT_AND_SKY: 373 // EXT model fits all params (including sky) 374 nParams = model->params->n; 375 psVectorInit (pcm->constraint->paramMask, 0); 376 break; 367 377 case PM_SOURCE_FIT_INDEX: 368 378 // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params -
branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitModel.c
r35768 r35961 66 66 opt->gainFactorMode = 0; 67 67 opt->chisqConvergence = true; 68 opt->isInteractive = false; 68 69 69 70 return opt; … … 247 248 myMin->gainFactorMode = options->gainFactorMode; 248 249 myMin->chisqConvergence = options->chisqConvergence; 250 myMin->isInteractive = options->isInteractive; 249 251 250 252 psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); -
branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitModel.h
r35768 r35961 37 37 int gainFactorMode; 38 38 bool chisqConvergence; 39 bool isInteractive; 39 40 } pmSourceFitOptions; 40 41 -
branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitPCM.c
r35768 r35961 68 68 myMin->chisqConvergence = fitOptions->chisqConvergence; 69 69 myMin->gainFactorMode = fitOptions->gainFactorMode; 70 myMin->isInteractive = fitOptions->isInteractive; 70 71 71 72 psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32); -
branches/eam_branches/ipp-20130711/psModules/src/objects/pmSourceFitSet.c
r35768 r35961 570 570 myMin->gainFactorMode = options->gainFactorMode; 571 571 myMin->chisqConvergence = options->chisqConvergence; 572 myMin->isInteractive = options->isInteractive; 572 573 573 574 psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F32);
Note:
See TracChangeset
for help on using the changeset viewer.
