IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Dec 28, 2004, 8:40:04 PM (22 years ago)
Author:
eugene
Message:

fixes for -Wall -Werror (mostly get_argument lines)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/cmd.astro/imfit.c

    r2598 r2843  
    11# include "astro.h"
    22
    3 # define SUBTRACT 0
    4 
    5 static float (*fitfunc)(float, float, float *, int, float *);
    6 
    73float fgaussTD (float, float, float *, int, float *);
    8 float fgalaxyTD (float, float, float *, int, float *);
    9 float fbarTD (float, float, float *, int, float *);
    10 void fix_ellipsegauss_pars (float *par, int Npar);
    11 
    12 /*** for some reason, fgaussTD fails occasionally if we use these two vars... */
    13 static float f1, f2;
    144
    155int imfit (int argc, char **argv) {
    166
    17   int i, j, N, done, Npts, Subtract;
    18   double Npar, DZ, DZ2;
    19   float dx, dy, chisq, nchisq;
    20   float **covar;
    21   float *V, *x, *y, *z, *dz, *par, tmpsky, DFact;
    22   int sx, sy, nx, ny;
    23   char *c, line[64];
     7  int i, j, N, Npts, Npar, Save, VERBOSE;
     8  int sx, sy, nx, ny, Nx, Ny;
     9  float chisq, ochisq, dchisq;
     10  float *x, *y, *z, *dz, *par, *V;
     11  char line[64];
    2412  Buffer *buf;
    2513
    26   Subtract = FALSE;
    27   if (N = get_argument (argc, argv, "-sub")) {
    28     remove_argument (N, &argc, argv);
    29     Subtract  = TRUE;
    30   }
    31 
    32   DFact = 1;
    33   if (N = get_argument (argc, argv, "-D")) {
    34     remove_argument (N, &argc, argv);
    35     DFact  = atof(argv[N]);
    36     remove_argument (N, &argc, argv);
    37   }
    38 
    39   fitfunc = fgaussTD;
    40   if (N = get_argument (argc, argv, "-gal")) {
    41     remove_argument (N, &argc, argv);
    42     fitfunc = fgalaxyTD;
    43   }
    44   if (N = get_argument (argc, argv, "-bar")) {
    45     remove_argument (N, &argc, argv);
    46     fitfunc = fbarTD;
     14  Save = FALSE;
     15  if ((N = get_argument (argc, argv, "-save"))) {
     16    remove_argument (N, &argc, argv);
     17    Save = TRUE;
     18  }
     19
     20  VERBOSE = FALSE;
     21  if ((N = get_argument (argc, argv, "-v"))) {
     22    remove_argument (N, &argc, argv);
     23    VERBOSE = TRUE;
    4724  }
    4825
    … …  
    5330
    5431  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    55 
    56   f1 = 1;
    57   if ((c = get_variable ("BETA1")) != (char *) NULL) f1 = atof (c);
    58 
    59   f2 = 1;
    60   if ((c = get_variable ("BETA2")) != (char *) NULL) f2 = atof (c);
    61 
    6232  sx = atof (argv[2]);
    6333  sy = atof (argv[3]);
    6434  nx = atof (argv[4]);
    6535  ny = atof (argv[5]);
    66 
    67   if ((sx < 0) || (sy < 0) ||
    68       (sx+nx > buf[0].matrix.Naxis[0]) ||
    69       (sy+ny > buf[0].matrix.Naxis[1])) {
    70     fprintf (stderr, "region out of range\n");
    71     return (FALSE);
    72   }
     36  Nx = buf[0].matrix.Naxis[0];
     37  Ny = buf[0].matrix.Naxis[1];
     38
     39  if (sx < 0) goto range;
     40  if (sy < 0) goto range;
     41  if (sx + nx >= Nx) goto range;
     42  if (sy + ny >= Ny) goto range;
    7343
    7444  Npar = 7;
    … …  
    7848  /* convert array z[x,y] to x[i], y[i], z[i] */
    7949  Npts = nx*ny;
    80   ALLOCATE (x, float, 2*Npts);
    81   ALLOCATE (y, float, 2*Npts);
    82   ALLOCATE (z, float, 2*Npts);
     50  ALLOCATE (x,  float, 2*Npts);
     51  ALLOCATE (y,  float, 2*Npts);
     52  ALLOCATE (z,  float, 2*Npts);
    8353  ALLOCATE (dz, float, 2*Npts);
    8454
    85   for (N = j = 0; j < ny; j++) {
     55  N = 0;
     56  for (j = 0; j < ny; j++) {
    8657    V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx;
    8758    for (i = 0; i < nx; i++, V++, N++) {
    88       dx = i + sx;
    89       dy = j + sy;
     59      x[N] = i + sx;
     60      y[N] = j + sy;
    9061      z[N] = *V;
    9162      dz[N] = 1.0 / *V;  /* temporarily no gain (or, gain = 1) */
    92       x[N] = dx;
    93       y[N] = dy;
    9463    }
    9564  }
    … …  
    9766  /* make initial guesses on parameters */
    9867  /* fitting X^2 / (2*Sx), Y^2 / (2*Sy) */
    99   get_aperture_stats (&buf[0].matrix, (int)(sx + 0.5*nx), (int)(sy + 0.5*ny), 15, 2, 65000);
    10068  par[0] = atof (get_variable ("Xg"));
    10169  par[1] = atof (get_variable ("Yg"));
    102   par[2] = SQ (atof (get_variable ("SXg")) / 2.355);
    103   par[3] = SQ (atof (get_variable ("SYg")) / 2.355);
    104   par[4] = 0;
    105   par[5] = atof (get_variable ("Zcg")) / 15.0;
    106   /* par[5] is peak intensity, Zcg is total flux, scale down some */
     70  par[2] = 2.35 * sqrt(2.0) / atof (get_variable ("SXg"));
     71  par[3] = 2.35 * sqrt(2.0) / atof (get_variable ("SYg"));
     72  par[4] = 0.0;
     73  par[5] = atof (get_variable ("Zpk"));
    10774  par[6] = atof (get_variable ("Sg"));
    10875
    109   /* soften the errors */
    110   for (i = 0; i < Npts; i++) {
    111     double tmp1;
    112     tmp1 = fitfunc (x[i], y[i], par, Npar, (float *) NULL);
    113     if (!finite (tmp1)) {
    114       fprintf (stderr, "nan for tmp @ %d (%f %f %f)\n", i, x[i], y[i], z[i]);
    115     } else {
    116       DZ = DFact * fabs (z[i] - tmp1);
    117       DZ2 = DFact * fabs (z[i] - tmp1);
    118       dz[i] *= 1 / (dz[i]*DZ + 1);  /* dz[i] = 1 / (DZ + 1/dz[i]) */
     76  /* run fit routine */
     77  ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fgaussTD, VERBOSE);
     78  dchisq = ochisq;
     79
     80  for (i = 0; (i < 20) && ((dchisq <= 0.0) || (dchisq > 0.1*(Npts - Npar))); i++) {
     81    chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fgaussTD, VERBOSE);
     82    dchisq = ochisq - chisq;
     83    ochisq = chisq;
     84    if (VERBOSE) fprintf (stderr, "dchisq: %f, Ndof: %d\n", dchisq, Npts - Npar);
     85  } 
     86
     87  /** create output image **/
     88  if (Save) {
     89    Buffer *out;
     90    float *Vi, *Vo, vr, vf, *dpar;
     91
     92    if ((out = SelectBuffer ("out",   ANYBUFFER, TRUE)) == NULL) return (FALSE);
     93    ALLOCATE (dpar, float, 2*Npar);
     94 
     95    free (out[0].header.buffer);
     96    free (out[0].matrix.buffer);
     97    strcpy (out[0].file, "(empty)");
     98    out[0].header.bitpix = -32;
     99    out[0].header.bscale = 1.0;
     100    out[0].header.bzero  = 0.0;
     101    out[0].header.unsign = FALSE;
     102    out[0].header.Naxes = 2;
     103    out[0].header.Naxis[0] = 2*nx;
     104    out[0].header.Naxis[1] = 2*ny;
     105    CreateBuffer (out);
     106
     107    /* four panels: 1) raw image. 2) fit  3) raw - fit   4) ?? */
     108    Vi = (float *)buf[0].matrix.buffer;
     109    Vo = (float *)out[0].matrix.buffer;
     110    for (j = 0; j < ny; j++) {
     111      for (i = 0; i < nx; i++) {
     112        vf = fgaussTD ((float)(i+sx), (float)(j+sy), par, Npar, dpar);
     113        vr = Vi[(i+sx)+(j+sy)*Nx];
     114        Vo[(i   )+(j   )*2*nx] = vr;
     115        Vo[(i+nx)+(j   )*2*nx] = vf;
     116        Vo[(i   )+(j+ny)*2*nx] = vr - vf;
     117        Vo[(i+nx)+(j+ny)*2*nx] = 0;
     118      }
    119119    }
    120   }
    121   /* run fit routine */
    122   chisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fitfunc);
    123   done = FALSE;
    124   for (i = 0; (i < 30) && !done; i++) {
    125     nchisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc);
    126     if (chisq - nchisq < 0.1) done = TRUE;
    127     if (nchisq > 0) chisq = nchisq;
    128   } 
    129   covar = mrq2dcovar (Npar);
     120    free (dpar);
     121  }
     122
     123  /* parameter renormalize */
     124  par[2] = 2.35*sqrt(2.0) / par[2];
     125  par[3] = 2.35*sqrt(2.0) / par[3];
     126  for (i = 0; i < Npar; i++) {
     127    sprintf (line, "PAR%d", i);
     128    set_variable (line, (double) par[i]);
     129    if (VERBOSE) fprintf (stderr, "%f ", par[i]);
     130  }
     131  if (VERBOSE) fprintf (stderr, "  %f\n", chisq/Npts);
     132
     133  free (par);
     134  free (x);
     135  free (y);
     136  free (z);
     137  free (dz);
     138
    130139  mrq2dfree (Npar);
     140  return (TRUE);
     141
     142range:
     143  fprintf (stderr, "region out of range\n");
     144  return (FALSE);
     145}
     146
     147/* pars: x, y, sx, sy, sxy, I, sky */
     148float fgaussTD (float x, float y, float *par, int Npar, float *dpar) {
     149
     150  float X, Y, px, py;
     151  float z, r, q, f;
     152
     153  X = x - par[0];
     154  Y = y - par[1];
     155 
     156  px = par[2]*X;
     157  py = par[3]*Y;
     158
     159  z = 0.5*SQ(px) + 0.5*SQ(py) + par[4]*X*Y;
     160  r = exp (-z);
     161  q = par[5]*r;
     162  f = q + par[6];
     163
     164  dpar[0] = q*(2*px*par[2] + par[4]*Y);
     165  dpar[1] = q*(2*py*par[3] + par[4]*X);
     166  dpar[2] = -2*q*px*X;
     167  dpar[3] = -2*q*py*Y;
     168  dpar[4] = -q*X*Y;
     169  dpar[5] = +r;
     170  dpar[6] = +1;
     171
     172  return (f);
     173}
     174/* alternate for exp(-z) */
     175/* E = 1.0 / (1 + f1*Z + f2*(0.5*Z*Z + 0.1666666*Z*Z*Z)); */  /* ~ exp (-Z) */
     176/* E = 1.0 / (1 +    Z +    (0.5*Z*Z + 0.1666666*Z*Z*Z)); */  /* ~ exp (-Z) */
     177
     178# if (0)
     179/* pars: x, y, sx, sy, sxy, I, sky */
     180float fgalaxyTD (float x, float y, float *par, int Npar, float *dpar) {
     181
     182  float X, Y, Z, E, F, q, R, f, p2, p3;
     183
     184  X = x - par[0];
     185  Y = y - par[1];
     186 
     187  p2 = X / par[2];
     188  p3 = Y / par[3];
     189
     190  Z = sqrt (0.5*p2*X + X*Y*par[4] + 0.5*p3*Y);                 /* R */
     191  E = 1.0 / (1 + Z);   
     192
     193  q = par[5] * E;
     194  R = q*E;
     195  F = 0.5 / Z;
     196 
     197  f = q + par[6];
     198
     199  dpar[0] = F*R*(p2 + par[4]*Y);
     200  dpar[1] = F*R*(p3 + par[4]*X);
     201  dpar[2] = F*0.5*R*p2*p2;
     202  dpar[3] = F*0.5*R*p3*p3;
     203  dpar[4] = -R*X*Y*F;
     204   
     205  dpar[5] = E;
     206  dpar[6] = 1;
     207  return (f);
     208}
     209
     210/* pars: x, y, sx, sy, sxy, I, sky */
     211float fbarTD (float x, float y, float *par, int Npar, float *dpar) {
     212
     213  float X, Y, Z, E, F, q, R, f, p2, p3;
     214
     215  X = x - par[0];
     216  Y = y - par[1];
     217 
     218  p2 = X / par[2];
     219  p3 = Y / par[3];
     220
     221  Z = 0.5*p2*X + X*Y*par[4] + 0.5*p3*Y;                 /* R */
     222  E = 1.0 / (1 + Z*Z*Z);   
     223
     224  q = par[5] * E;
     225  F = 3*Z*Z;
     226  R = q*E*F;
     227 
     228  f = q + par[6];
     229
     230  dpar[0] = R*(p2 + par[4]*Y);
     231  dpar[1] = R*(p3 + par[4]*X);
     232  dpar[2] = 0.5*R*p2*p2;
     233  dpar[3] = 0.5*R*p3*p3;
     234  dpar[4] = -R*X*Y;
     235   
     236  dpar[5] = E;
     237  dpar[6] = 1;
     238  return (f);
     239}
     240
     241/* convert from x,y to major,minor */
     242void fix_ellipsegauss_pars (float *par, int Npar) {
     243
     244  float p2, p4, angle, t1, t2, tmp, area;
     245
     246  /* par[0], par[1] = Xo, Yo - stay the same */
     247
     248  p2 = 1/par[2];
     249  p4 = 1/par[3];
     250
     251  angle = 0.5 * atan2 (-2*par[4], p4 - p2);
     252
     253  tmp = sqrt (SQ(p2 - p4) + 4*SQ(par[4]));
     254  t1 = (p2 + p4 + tmp) / 2;
     255  t2 = t1 - tmp;
     256
     257  par[2] = 2.35482*sqrt(1/t2);
     258  par[3] = 2.35482*sqrt(1/t1);
     259  par[4] = DEG_RAD * angle;
     260
     261  area = 2*M_PI/sqrt(t1*t2);
     262
     263  par[5] *= area;
     264
     265}
     266# endif
     267
     268/***  options for later
     269
     270  Subtract = FALSE;
     271  if ((N = get_argument (argc, argv, "-sub"))) {
     272    remove_argument (N, &argc, argv);
     273    Subtract  = TRUE;
     274  }
     275
     276  DFact = 1;
     277  if ((N = get_argument (argc, argv, "-D"))) {
     278    remove_argument (N, &argc, argv);
     279    DFact  = atof(argv[N]);
     280    remove_argument (N, &argc, argv);
     281  }
     282
     283  fitfunc = fgaussTD;
     284  if ((N = get_argument (argc, argv, "-gal"))) {
     285    remove_argument (N, &argc, argv);
     286    fitfunc = fgalaxyTD;
     287  }
     288  if ((N = get_argument (argc, argv, "-bar"))) {
     289    remove_argument (N, &argc, argv);
     290    fitfunc = fbarTD;
     291  }
     292
     293
     294  f1 = 1;
     295  if ((c = get_variable ("BETA1")) != (char *) NULL) f1 = atof (c);
     296
     297  f2 = 1;
     298  if ((c = get_variable ("BETA2")) != (char *) NULL) f2 = atof (c);
    131299
    132300  if (Subtract) {
    … …  
    144312  }
    145313
    146   fix_ellipsegauss_pars (par, Npar);
    147 
    148   for (i = 0; i < Npar; i++) {
    149     sprintf (line, "PAR%d", i);
    150     set_variable (line, (double) par[i]);
    151     fprintf (stderr, "%f ", par[i]);
    152   }
    153   fprintf (stderr, "  %f\n", nchisq/(1.0*Npts));
    154 
    155   free (par);
    156   free (x);
    157   free (y);
    158   free (z);
    159   free (dz);
    160 
    161   return (TRUE);
    162 
    163 }
    164 
    165 /* pars: x, y, sx, sy, sxy, I, sky */
    166 float fgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    167 
    168   float X, Y, Z, E, F, q, R, f, p2, p3;
    169 
    170   X = x - par[0];
    171   Y = y - par[1];
    172  
    173   p2 = X / par[2];
    174   p3 = Y / par[3];
    175 
    176   Z = 0.5*p2*X + X*Y*par[4] + 0.5*p3*Y;                 /* R^2 */
    177   /* E = 1.0 / (1 + f1*Z + f2*(0.5*Z*Z + 0.1666666*Z*Z*Z)); */  /* ~ exp (-Z) */
    178   E = 1.0 / (1 + Z + (0.5*Z*Z + 0.1666666*Z*Z*Z));   /* ~ exp (-Z) */
    179   F = (f1 + f2*(Z + 0.5*Z*Z));
    180 
    181   q = par[5] * E;
    182   R = q*E*F;
    183  
    184   f = q + par[6];
    185 
    186   if (dpar != (float *) NULL) {
    187     dpar[0] = R*(p2 + par[4]*Y);
    188     dpar[1] = R*(p3 + par[4]*X);
    189     dpar[2] = 0.5*R*p2*p2;
    190     dpar[3] = 0.5*R*p3*p3;
    191     dpar[4] = -R*X*Y;
    192    
    193     dpar[5] = E;
    194     dpar[6] = 1;
    195   }
    196 
    197   /*
    198   if ((fabs(X) < 2) && (fabs(Y) < 2)) {
    199     int i;
    200 
    201     for (i = 0; i < Npar; i++) {
    202       fprintf (stderr, "%d  %f  %f\n", i, par[i], dpar[i]);
    203     }
    204     fprintf (stderr, "func: (%f, %f)  %f\n", X, Y, f);
    205   }
    206   */
    207 
    208   if (!finite (f)) {
    209     fprintf (stderr, "nan for fgaussTD\n");
    210   }
    211 
    212   return (f);
    213 
    214 }
    215 
    216 /* pars: x, y, sx, sy, sxy, I, sky */
    217 float fgalaxyTD (float x, float y, float *par, int Npar, float *dpar) {
    218 
    219   float X, Y, Z, E, F, q, R, f, p2, p3;
    220 
    221   X = x - par[0];
    222   Y = y - par[1];
    223  
    224   p2 = X / par[2];
    225   p3 = Y / par[3];
    226 
    227   Z = sqrt (0.5*p2*X + X*Y*par[4] + 0.5*p3*Y);                 /* R */
    228   E = 1.0 / (1 + Z);   
    229 
    230   q = par[5] * E;
    231   R = q*E;
    232   F = 0.5 / Z;
    233  
    234   f = q + par[6];
    235 
    236   if (dpar != (float *) NULL) {
    237     dpar[0] = F*R*(p2 + par[4]*Y);
    238     dpar[1] = F*R*(p3 + par[4]*X);
    239     dpar[2] = F*0.5*R*p2*p2;
    240     dpar[3] = F*0.5*R*p3*p3;
    241     dpar[4] = -R*X*Y*F;
    242    
    243     dpar[5] = E;
    244     dpar[6] = 1;
    245   }
    246 
    247   /*
    248   if ((fabs(X) < 2) && (fabs(Y) < 2)) {
    249     int i;
    250 
    251     for (i = 0; i < Npar; i++) {
    252       fprintf (stderr, "%d  %f  %f\n", i, par[i], dpar[i]);
    253     }
    254     fprintf (stderr, "func: (%f, %f)  %f\n", X, Y, f);
    255   }
    256   */
    257   return (f);
    258 
    259 }
    260 
    261 
    262 /* pars: x, y, sx, sy, sxy, I, sky */
    263 float fbarTD (float x, float y, float *par, int Npar, float *dpar) {
    264 
    265   float X, Y, Z, E, F, q, R, f, p2, p3;
    266 
    267   X = x - par[0];
    268   Y = y - par[1];
    269  
    270   p2 = X / par[2];
    271   p3 = Y / par[3];
    272 
    273   Z = 0.5*p2*X + X*Y*par[4] + 0.5*p3*Y;                 /* R */
    274   E = 1.0 / (1 + Z*Z*Z);   
    275 
    276   q = par[5] * E;
    277   F = 3*Z*Z;
    278   R = q*E*F;
    279  
    280   f = q + par[6];
    281 
    282   if (dpar != (float *) NULL) {
    283     dpar[0] = R*(p2 + par[4]*Y);
    284     dpar[1] = R*(p3 + par[4]*X);
    285     dpar[2] = 0.5*R*p2*p2;
    286     dpar[3] = 0.5*R*p3*p3;
    287     dpar[4] = -R*X*Y;
    288    
    289     dpar[5] = E;
    290     dpar[6] = 1;
    291   }
    292 
    293   /*
    294   if ((fabs(X) < 2) && (fabs(Y) < 2)) {
    295     int i;
    296 
    297     for (i = 0; i < Npar; i++) {
    298       fprintf (stderr, "%d  %f  %f\n", i, par[i], dpar[i]);
    299     }
    300     fprintf (stderr, "func: (%f, %f)  %f\n", X, Y, f);
    301   }
    302   */
    303   return (f);
    304 
    305 }
    306 
    307 
    308 void fix_ellipsegauss_pars (float *par, int Npar) {
    309 
    310   float p2, p4, angle, t1, t2, tmp, area;
    311 
    312   /* par[0], par[1] = Xo, Yo - stay the same */
    313 
    314   p2 = 1/par[2];
    315   p4 = 1/par[3];
    316 
    317   angle = 0.5 * atan2 (-2*par[4], p4 - p2);
    318 
    319   tmp = sqrt (SQ(p2 - p4) + 4*SQ(par[4]));
    320   t1 = (p2 + p4 + tmp) / 2;
    321   t2 = t1 - tmp;
    322 
    323   par[2] = 2.35482*sqrt(1/t2);
    324   par[3] = 2.35482*sqrt(1/t1);
    325   par[4] = DEG_RAD * angle;
    326 
    327   area = 2*M_PI/sqrt(t1*t2);
    328 
    329   par[5] *= area;
    330 
    331 }
    332 
    333 
     314***/
Note: See TracChangeset for help on using the changeset viewer.