IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 2938


Ignore:
Timestamp:
Jan 8, 2005, 2:45:18 PM (22 years ago)
Author:
eugene
Message:

minor changes

Location:
trunk/Ohana/src/opihi
Files:
7 edited

Legend:

Unmodified
Added
Removed

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

    r2886 r2938  
    11# include "astro.h"
     2# define FFACTOR 200
     3# define FSCALE 1.2
     4
     5float mrq2dchi (float *x, float *t, float *y, float *dy, int Npts,
     6                float *par, int Npar,
     7                float (funcs)(float, float, float *, int, float *));
    28
    39float fgaussTD (float, float, float *, int, float *);
    410float pgaussTD (float, float, float *, int, float *);
     11float sgaussTD (float, float, float *, int, float *);
    512float vgaussTD (float, float, float *, int, float *);
    6 float qgaussTD (float, float, float *, int, float *);
     13
     14float PgaussTD (float, float, float *, int, float *);
     15float SgaussTD (float, float, float *, int, float *);
     16
     17float tgaussTD (float, float, float *, int, float *);
    718float rgaussTD (float, float, float *, int, float *);
    8 float sgaussTD (float, float, float *, int, float *);
    919
    1020static float Npow;
     21static float Npin;
     22static int MissingVariable = FALSE;
     23static float *fpar;
     24static int Nfpar;
     25
     26int mrq2dlimits (float *pmin, float *pmax, int Npar);
     27
     28float get_variable_default (char *name, float dvalue) {
     29
     30  char *value;
     31  float fvalue;
     32
     33  value = get_variable (name);
     34  if (value == NULL) {
     35    MissingVariable = TRUE;
     36    return (dvalue);
     37  }
     38  fvalue = atof (value);
     39  return (fvalue);
     40}
    1141
    1242int imfit (int argc, char **argv) {
    1343
    14   int i, j, N, Npts, Npar, Save, VERBOSE;
     44  int i, j, N, Npts, Npar, Save, VERBOSE, ShapeVariation;
    1545  int sx, sy, nx, ny, Nx, Ny;
    1646  float chisq, ochisq, dchisq, Gain, RDnoise;
    17   float *x, *y, *z, *dz, *par, *V;
     47  float *x, *y, *z, *dz, *par, *V, *parmin, *parmax;
    1848  char line[64];
    1949  Buffer *buf;
    … …  
    2656  }
    2757
     58  ShapeVariation = FALSE;
     59  if ((N = get_argument (argc, argv, "-shapes"))) {
     60    remove_argument (N, &argc, argv);
     61    ShapeVariation = TRUE;
     62  }
     63
    2864  /* Gain in e/DN */
    2965  Gain = 1.0;
    … …  
    4985
    5086  /* set fitting function */
     87  Nfpar = 0;
    5188  Npar = 7;
    5289  fitfunc = fgaussTD;
    … …  
    6097      fitfunc = pgaussTD;
    6198    }
    62     if (!strcmp(argv[N], "qgauss")) {
    63       fitfunc = qgaussTD;
    64       Npar = 9;
    65     }
    6699    if (!strcmp(argv[N], "vgauss")) {
    67100      fitfunc = vgaussTD;
    68101      Npar = 9;
    69102    }
    70     if (!strcmp(argv[N], "rgauss")) {
    71       fitfunc = rgaussTD;
    72       Npar = 11;
    73     }
    74103    if (!strcmp(argv[N], "sgauss")) {
    75104      fitfunc = sgaussTD;
    76105      Npar = 10;
     106    }
     107    if (!strcmp(argv[N], "rgauss")) {
     108      fitfunc = rgaussTD;
     109      Npar = 10;
     110    }
     111    if (!strcmp(argv[N], "tgauss")) {
     112      fitfunc = tgaussTD;
     113      Npar = 10;
     114    }
     115    if (!strcmp(argv[N], "Pgauss")) {
     116      fitfunc = PgaussTD;
     117      Npar = 4;
     118      Nfpar = 3;
     119    }
     120    if (!strcmp(argv[N], "Sgauss")) {
     121      fitfunc = SgaussTD;
     122      Npar = 4;
     123      Nfpar = 6;
    77124    }
    78125    if (fitfunc == NULL) {
    … …  
    102149  if (sy + ny >= Ny) goto range;
    103150
    104   ALLOCATE (par, float, 2*Npar);
     151  ALLOCATE (par, float, Npar);
    105152  bzero (par, Npar*sizeof(float));
     153  if (Nfpar > 0) {
     154    ALLOCATE (fpar, float, Nfpar);
     155    bzero (fpar, Nfpar*sizeof(float));
     156  }   
    106157
    107158  /* convert array z[x,y] to x[i], y[i], z[i] */
    … …  
    127178
    128179  /* make initial guesses on parameters */
    129   par[0] = atof (get_variable ("Xg"));
    130   par[1] = atof (get_variable ("Yg"));
    131   par[2] = 2.35 * sqrt(2.0) / atof (get_variable ("SXg"));
    132   par[3] = 2.35 * sqrt(2.0) / atof (get_variable ("SYg"));
    133   par[4] = 0.0;
    134   par[5] = atof (get_variable ("Zpk"));
    135   par[6] = atof (get_variable ("Sg"));
    136   if (Npar == 9) {
     180  if (Npar >= 7) {
     181    par[0] = get_variable_default ("Xg", 0);
     182    par[1] = get_variable_default ("Yg", 0);
     183    par[2] = 2.35 * sqrt(2.0) / get_variable_default ("SXg", 2.0);
     184    par[3] = 2.35 * sqrt(2.0) / get_variable_default ("SYg", 2.0);
     185    par[4] = 0.0;
     186    par[5] = get_variable_default ("Zpk", 10000);
     187    par[6] = get_variable_default ("Sg", 0.0);
     188  }
     189  if (Npar == 4) {
     190    par[0] = get_variable_default ("Xg", 0);
     191    par[1] = get_variable_default ("Yg", 0);
     192    par[2] = get_variable_default ("Zpk", 10000);
     193    par[3] = get_variable_default ("Sg", 0.0);
     194  }
     195  if (Npar == 9) { /** vgauss **/
    137196    par[7] = par[8] = 1;
    138197  }
    139   if (Npar == 11) {
    140     par[7] = 0.2*par[5];
    141     par[8] = 0.25*par[2];
    142     par[9] = 0.25*par[3];
    143     par[10] = 0.0;
    144   }
    145   if (Npar == 10) {
    146     Npow = atof (get_variable ("Npow"));
    147     par[7] = 2.35 * sqrt(2.0) / atof (get_variable ("SXf"));
    148     par[8] = 2.35 * sqrt(2.0) / atof (get_variable ("SYf"));
     198  if (Npar == 10) { /** sgauss, tgauss, rgauss **/
     199    Npow = get_variable_default ("Npow", 2.25);
     200    Npin = get_variable_default ("Npin", 1.00);
     201    par[7] = 2.35 * sqrt(2.0) / get_variable_default ("SXf", 15.0);
     202    par[8] = 2.35 * sqrt(2.0) / get_variable_default ("SYf", 15.0);
    149203    par[9] = 0.0;
    150204  }
    151 
     205  if (fitfunc == PgaussTD) {
     206    fpar[0] = 2.35 * sqrt(2.0) / get_variable_default ("SXg", 15.0);
     207    fpar[1] = 2.35 * sqrt(2.0) / get_variable_default ("SYg", 15.0);
     208    fpar[2] = get_variable_default ("SXYg", 15.0);
     209  }   
     210  if (fitfunc == SgaussTD) {
     211    fpar[0] = 2.35 * sqrt(2.0) / get_variable_default ("SXg", 15.0);
     212    fpar[1] = 2.35 * sqrt(2.0) / get_variable_default ("SYg", 15.0);
     213    fpar[2] = get_variable_default ("SXYg", 15.0);
     214    fpar[3] = 2.35 * sqrt(2.0) / get_variable_default ("SXf", 15.0);
     215    fpar[4] = 2.35 * sqrt(2.0) / get_variable_default ("SYf", 15.0);
     216    fpar[5] = get_variable_default ("SXYf", 15.0);
     217  }   
    152218
    153219  /* run fit routine */
    … …  
    159225    dchisq = ochisq - chisq;
    160226    ochisq = chisq;
    161     if (VERBOSE) fprintf (stderr, "dchisq: %f, Ndof: %d\n", dchisq, Npts - Npar);
    162227  } 
    163228
    164   /** create output image (keep in sky : par[6]?) **/
     229  if (ShapeVariation) {
     230    /* find dChi/dSx and dChi/dSy given by increasing shape terms by 5% */
     231    float tp1, tp2, chix, chiy;
     232    chix = chiy = 0;
     233    if (fitfunc == sgaussTD) {
     234      tp1 = par[2];
     235      tp2 = par[7];
     236      par[2] = par[2]*1.05;
     237      par[7] = par[7]*1.05;
     238      chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
     239      par[2] = tp1;
     240      par[7] = tp2;
     241
     242      tp1 = par[3];
     243      tp2 = par[8];
     244      par[3] = par[3]*1.05;
     245      par[8] = par[8]*1.05;
     246      chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
     247      par[3] = tp1;
     248      par[8] = tp2;
     249    }
     250    if (fitfunc == pgaussTD) {
     251      tp1 = par[2];
     252      par[2] = par[2]*1.05;
     253      chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
     254      par[2] = tp1;
     255
     256      tp1 = par[3];
     257      par[3] = par[3]*1.05;
     258      chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
     259      par[3] = tp1;
     260    }
     261    if (fitfunc == SgaussTD) {
     262      tp1 = par[0];
     263      tp2 = par[3];
     264      par[0] = par[0]*1.05;
     265      par[3] = par[3]*1.05;
     266      chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
     267      par[0] = tp1;
     268      par[3] = tp2;
     269
     270      tp1 = par[1];
     271      tp2 = par[4];
     272      par[1] = par[1]*1.05;
     273      par[4] = par[4]*1.05;
     274      chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
     275      par[1] = tp1;
     276      par[4] = tp2;
     277    }
     278    if (fitfunc == PgaussTD) {
     279      tp1 = par[0];
     280      par[0] = par[0]*1.05;
     281      chix = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
     282      par[0] = tp1;
     283
     284      tp1 = par[1];
     285      par[1] = par[1]*1.05;
     286      chiy = mrq2dchi (x, y, z, dz, Npts, par, Npar, fitfunc) - chisq;
     287      par[1] = tp1;
     288    }
     289    set_variable ("dChiX", chix/chisq);
     290    set_variable ("dChiY", chiy/chisq);
     291  }
     292
     293  /** create output image (keep in sky) **/
    165294  if (Save) {
    166295    Buffer *out;
    167     float *Vi, *Vo, vr, vf, *dpar;
     296    float *Vi, *Vo, vr, vf, sky;
    168297
    169298    if ((out = SelectBuffer ("out",   ANYBUFFER, TRUE)) == NULL) return (FALSE);
    170     ALLOCATE (dpar, float, 2*Npar);
    171  
    172299    free (out[0].header.buffer);
    173300    free (out[0].matrix.buffer);
    … …  
    181308    out[0].header.Naxis[1] = 2*ny;
    182309    CreateBuffer (out);
     310    if (Npar == 4) {
     311      sky = par[3];
     312    } else {
     313      sky = par[6];
     314    }
    183315
    184316    /* four panels: 1) raw image. 2) fit  3) raw - fit   4) ?? */
    … …  
    187319    for (j = 0; j < ny; j++) {
    188320      for (i = 0; i < nx; i++) {
    189         vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, dpar);
     321        vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, NULL);
    190322        vr = Vi[(i+sx)+(j+sy)*Nx];
    191323        Vo[(i   )+(j   )*2*nx] = vr;
    192324        Vo[(i+nx)+(j   )*2*nx] = vf;
    193         Vo[(i   )+(j+ny)*2*nx] = vr - vf + par[6];
    194         Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf) + par[6];
     325        Vo[(i   )+(j+ny)*2*nx] = vr - vf + sky;
     326        Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf) + sky;
    195327      }
    196328    }
    197     free (dpar);
    198329  }
    199330
    200331  /* parameter renormalize */
    201   par[2] = 2.35*sqrt(2.0) / par[2];
    202   par[3] = 2.35*sqrt(2.0) / par[3];
    203   if (Npar == 11) {
    204     par[8] = 2.35*sqrt(2.0) / par[8];
    205     par[9] = 2.35*sqrt(2.0) / par[9];
    206   }
    207   if (Npar == 10
    208 ) {
    209     par[7] = 2.35*sqrt(2.0) / par[7];
    210     par[8] = 2.35*sqrt(2.0) / par[8];
    211   }
    212   for (i = 0; i < Npar; i++) {
    213     sprintf (line, "PAR%d", i);
    214     set_variable (line, (double) par[i]);
    215     if (VERBOSE) fprintf (stderr, "%f ", par[i]);
    216   }
    217   if (VERBOSE) fprintf (stderr, "  %f\n", chisq/Npts);
    218   set_variable ("ChiSq", (double) (chisq/Npts));
    219 
    220   free (par);
     332  if (Npar >= 7) {
     333    set_variable ("Xg",  par[0]);
     334    set_variable ("Yg",  par[1]);
     335    set_variable ("SXg", 2.35 * sqrt(2.0) / par[2]);
     336    set_variable ("SYg", 2.35 * sqrt(2.0) / par[3]);
     337    set_variable ("SXYg", par[4]);
     338    set_variable ("Zpk", par[5]);
     339    set_variable ("Sg",  par[6]);
     340  }
     341  if (Npar == 9) {
     342    set_variable ("SXf", par[7]);
     343    set_variable ("SYf", par[8]);
     344  }
     345  if (Npar == 10) {
     346    set_variable ("SXf", 2.35 * sqrt(2.0) / par[7]);
     347    set_variable ("SYf", 2.35 * sqrt(2.0) / par[8]);
     348    set_variable ("SXYf", par[9]);
     349  }
     350  if (fitfunc == PgaussTD) {
     351    set_variable ("Xg",  par[0]);
     352    set_variable ("Yg",  par[1]);
     353    set_variable ("Zpk", par[2]);
     354    set_variable ("Sg",  par[3]);
     355    set_variable ("SXg", 2.35 * sqrt(2.0) / fpar[0]);
     356    set_variable ("SYg", 2.35 * sqrt(2.0) / fpar[1]);
     357    set_variable ("SXYg", fpar[2]);
     358  }
     359  if (fitfunc == SgaussTD) {
     360    set_variable ("Xg",  par[0]);
     361    set_variable ("Yg",  par[1]);
     362    set_variable ("Zpk", par[2]);
     363    set_variable ("Sg",  par[3]);
     364    set_variable ("SXg", 2.35 * sqrt(2.0) / fpar[0]);
     365    set_variable ("SYg", 2.35 * sqrt(2.0) / fpar[1]);
     366    set_variable ("SXYg", fpar[2]);
     367    set_variable ("SXf", 2.35 * sqrt(2.0) / fpar[3]);
     368    set_variable ("SYf", 2.35 * sqrt(2.0) / fpar[4]);
     369    set_variable ("SXYf", fpar[5]);
     370  }   
     371  set_variable ("ChiSq", chisq/(Npts - Npar));
     372
    221373  free (x);
    222374  free (y);
    223375  free (z);
    224376  free (dz);
     377  free (par);
     378  if (Nfpar > 0) free (fpar);
    225379
    226380  mrq2dfree (Npar);
    … …  
    232386}
    233387
    234 /* pars: x, y, sx, sy, sxy, I, sky */
     388/* real 2D gaussian -- x, y, sx, sy, sxy, I, sky */
    235389float fgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    236390
    … …  
    249403  f = q + par[6];
    250404
    251   dpar[0] = q*(2*px*par[2] + par[4]*Y);
    252   dpar[1] = q*(2*py*par[3] + par[4]*X);
    253   dpar[2] = -2*q*px*X;
    254   dpar[3] = -2*q*py*Y;
    255   dpar[4] = -q*X*Y;
    256   dpar[5] = +r;
    257   dpar[6] = +1;
    258 
    259   return (f);
    260 }
    261 
    262 /***** fix the derivatives.  q should not be in function */
    263 /* pars: x, y, sx, sy, sxy, I, sky */
     405  if (dpar != NULL) {
     406    dpar[0] = q*(2*px*par[2] + par[4]*Y);
     407    dpar[1] = q*(2*py*par[3] + par[4]*X);
     408    dpar[2] = -2*q*px*X;
     409    dpar[3] = -2*q*py*Y;
     410    dpar[4] = -q*X*Y;
     411    dpar[5] = +r;
     412    dpar[6] = +1;
     413  }
     414  return (f);
     415}
     416
     417/* pseudo 2D gaussian -- x, y, sx, sy, sxy, I, sky */
    264418float pgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    265419
    … …  
    279433  /* note difference from gaussian: q = par[5]*r */
    280434
    281   dpar[0] = q*(2*px*par[2] + par[4]*Y);
    282   dpar[1] = q*(2*py*par[3] + par[4]*X);
    283   dpar[2] = -2*q*px*X;
    284   dpar[3] = -2*q*py*Y;
    285   dpar[4] = -q*X*Y;
    286   dpar[5] = +r;
    287   dpar[6] = +1;
    288 
    289   return (f);
    290 }
    291 
    292 /* pars: x, y, sx, sy, sxy, I, sky, f1, f2 */
     435  if (dpar != NULL) {
     436    dpar[0] = q*(2*px*par[2] + par[4]*Y);
     437    dpar[1] = q*(2*py*par[3] + par[4]*X);
     438    dpar[2] = -2*q*px*X;
     439    dpar[3] = -2*q*py*Y;
     440    dpar[4] = -q*X*Y;
     441    dpar[5] = +r;
     442    dpar[6] = +1;
     443  }
     444  return (f);
     445}
     446
     447/* pseudo 2D gaussian -- x, y, (sx), (sy), (sxy), I, sky */
     448float PgaussTD (float x, float y, float *par, int Npar, float *dpar) {
     449
     450  float X, Y, px, py;
     451  float z, r, q, f;
     452
     453  /* par -> fpar: (2,0), (3,1), (4,2) */
     454
     455  X = x - par[0];
     456  Y = y - par[1];
     457 
     458  px = fpar[0]*X;
     459  py = fpar[1]*Y;
     460
     461  z = 0.5*SQ(px) + 0.5*SQ(py) + fpar[2]*X*Y;
     462  r = 1.0 / (1 + z + 0.5*z*z*(1 + z/3)); /* ~ exp (-Z) */
     463  f = par[2]*r + par[3];
     464  q = par[2]*r*r*(1 + z + 0.5*z*z);
     465  /* note difference from gaussian: q = par[5]*r */
     466
     467  if (dpar != NULL) {
     468    dpar[0] = q*(2*px*fpar[0] + fpar[2]*Y);
     469    dpar[1] = q*(2*py*fpar[1] + fpar[2]*X);
     470    dpar[2] = +r;
     471    dpar[3] = +1;
     472  }
     473  return (f);
     474}
     475
     476/* pseudo 2D gaussian with floating 2nd and 3rd order terms -- x, y, sx, sy, sxy, I, sky, f1, f2 */
    293477float vgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    294478
    … …  
    309493  /* note difference from gaussian: q = par[5]*r */
    310494
    311   dpar[0] = q*(2*px*par[2] + par[4]*Y);
    312   dpar[1] = q*(2*py*par[3] + par[4]*X);
    313   dpar[2] = -2*q*px*X;
    314   dpar[3] = -2*q*py*Y;
    315   dpar[4] = -q*X*Y;
    316   dpar[5] = +r;
    317   dpar[6] = +1;
    318   dpar[7] = -100*par[5]*r*r*k;
    319   dpar[8] = -100*par[5]*r*r*par[7]*(z*z*z)/6;
    320 
    321   return (f);
    322 }
    323 
     495  if (dpar != NULL) {
     496    dpar[0] = q*(2*px*par[2] + par[4]*Y);
     497    dpar[1] = q*(2*py*par[3] + par[4]*X);
     498    dpar[2] = -2*q*px*X;
     499    dpar[3] = -2*q*py*Y;
     500    dpar[4] = -q*X*Y;
     501    dpar[5] = +r;
     502    dpar[6] = +1;
     503    dpar[7] = -100*par[5]*r*r*k;
     504    dpar[8] = -100*par[5]*r*r*par[7]*(z*z*z)/6;
     505  }
     506  return (f);
     507}
     508
     509/* two components: (1 + z_1 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */
     510float sgaussTD (float x, float y, float *par, int Npar, float *dpar) {
     511
     512  float X, Y, px1, py1, px2, py2;
     513  float z1, z2, r, q1, q2, f, f1, f2;
     514
     515  X = x - par[0];
     516  Y = y - par[1];
     517 
     518  px1 = par[2]*X;
     519  py1 = par[3]*Y;
     520  px2 = par[7]*X;
     521  py2 = par[8]*Y;
     522
     523  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
     524  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y;
     525
     526  r = 1.0 / (1 + z1 + pow(z2,Npow));
     527  f = par[5]*r + par[6];
     528
     529  q1 = par[5]*SQ(r);
     530  q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-1));
     531
     532  if (dpar != NULL) {
     533    dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y);
     534    dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X);
     535
     536    /* these fudge factors impede the growth of par[2] beyond par[7] */
     537    f1 = fabs(par[7]) / fabs(par[2]);
     538    f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1;
     539    dpar[2] = -2*q1*px1*X*f2;
     540
     541    f1 = fabs(par[8]) / fabs(par[3]);
     542    f2 = (f1 < FSCALE) ? 1 : FFACTOR*(f1 - FSCALE) + 1;
     543    dpar[3] = -2*q1*py1*Y*f2;
     544
     545    dpar[4] = -q1*X*Y;
     546    dpar[5] = +r;
     547    dpar[6] = +1;
     548    dpar[7] = -2*q2*px2*X;
     549    dpar[8] = -2*q2*py2*Y;
     550    dpar[9] = -q2*X*Y;
     551  }
     552  return (f);
     553}
     554
     555/* two components: (1 + z_1 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */
     556float SgaussTD (float x, float y, float *par, int Npar, float *dpar) {
     557
     558  float X, Y, px1, py1, px2, py2;
     559  float z1, z2, r, q1, q2, f, f1, f2;
     560
     561  X = x - par[0];
     562  Y = y - par[1];
     563 
     564  px1 = fpar[0]*X;
     565  py1 = fpar[1]*Y;
     566  px2 = fpar[3]*X;
     567  py2 = fpar[4]*Y;
     568
     569  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + fpar[2]*X*Y;
     570  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + fpar[5]*X*Y;
     571
     572  r = 1.0 / (1 + z1 + pow(z2,Npow));
     573  f = par[2]*r + par[3];
     574
     575  q1 = par[2]*SQ(r);
     576  q2 = par[2]*SQ(r)*Npow*pow(z2,(Npow-1));
     577
     578  if (dpar != NULL) {
     579    dpar[0] = q1*(2*px1*fpar[0] + fpar[2]*Y) + q2*(2*px2*fpar[3] + fpar[5]*Y);
     580    dpar[1] = q1*(2*py1*fpar[1] + fpar[2]*X) + q2*(2*py2*fpar[4] + fpar[5]*X);
     581    dpar[2] = +r;
     582    dpar[3] = +1;
     583  }
     584  return (f);
     585}
     586
     587/* two components: (1 + z_1 + 0.5*z_1^2 + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */
     588float rgaussTD (float x, float y, float *par, int Npar, float *dpar) {
     589
     590  float X, Y, px1, py1, px2, py2;
     591  float z1, z2, r, q1, q2, f;
     592
     593  X = x - par[0];
     594  Y = y - par[1];
     595 
     596  px1 = par[2]*X;
     597  py1 = par[3]*Y;
     598  px2 = par[7]*X;
     599  py2 = par[8]*Y;
     600
     601  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
     602  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y;
     603
     604  r = 1.0 / (1 + z1 + 0.5*SQ(z1)+ pow(z2,Npow));
     605  f = par[5]*r + par[6];
     606
     607  q1 = par[5]*SQ(r)*(1 + z1);
     608  q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-1));
     609
     610  if (dpar != NULL) {
     611    dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y);
     612    dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X);
     613    dpar[2] = -2*q1*px1*X;
     614    dpar[3] = -2*q1*py1*Y;
     615    dpar[4] = -q1*X*Y;
     616    dpar[5] = +r;
     617    dpar[6] = +1;
     618    dpar[7] = -2*q2*px2*X;
     619    dpar[8] = -2*q2*py2*Y;
     620    dpar[9] = -q2*X*Y;
     621  }
     622  return (f);
     623}
     624
     625/* two components: (1 + z_1^M + z_2^N)^(-1) -- x, y, sx1, sy1, sxy1, I, sky, sx2, sy2, sxy2 */
     626float tgaussTD (float x, float y, float *par, int Npar, float *dpar) {
     627
     628  float X, Y, px1, py1, px2, py2;
     629  float z1, z2, r, q1, q2, f;
     630
     631  X = x - par[0];
     632  Y = y - par[1];
     633 
     634  px1 = par[2]*X;
     635  py1 = par[3]*Y;
     636  px2 = par[7]*X;
     637  py2 = par[8]*Y;
     638
     639  z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
     640  z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y;
     641
     642  r = 1.0 / (1 + pow(z1,Npin) + pow(z2,Npow));
     643  f = par[5]*r + par[6];
     644
     645  q1 = par[5]*SQ(r)*Npin*pow(z1,(Npin-1));
     646  q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-1));
     647
     648  if (dpar != NULL) {
     649    dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y);
     650    dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X);
     651    dpar[2] = -2*q1*px1*X*2;
     652    dpar[3] = -2*q1*py1*Y*2;
     653    dpar[4] = -q1*X*Y;
     654    dpar[5] = +r;
     655    dpar[6] = +1;
     656    dpar[7] = -2*q2*px2*X;
     657    dpar[8] = -2*q2*py2*Y;
     658    dpar[9] = -q2*X*Y;
     659  }
     660  return (f);
     661}
     662
     663/* alternate for exp(-z) */
     664/* E = 1.0 / (1 + f1*Z + f2*(0.5*Z*Z + 0.1666666*Z*Z*Z)); */  /* ~ exp (-Z) */
     665/* E = 1.0 / (1 +    Z +    (0.5*Z*Z + 0.1666666*Z*Z*Z)); */  /* ~ exp (-Z) */
     666
     667# if (0)
    324668/* pars: x, y, sx, sy, sxy, I, sky */
     669float fgalaxyTD (float x, float y, float *par, int Npar, float *dpar) {
     670
     671  float X, Y, Z, E, F, q, R, f, p2, p3;
     672
     673  X = x - par[0];
     674  Y = y - par[1];
     675 
     676  p2 = X / par[2];
     677  p3 = Y / par[3];
     678
     679  Z = sqrt (0.5*p2*X + X*Y*par[4] + 0.5*p3*Y);                 /* R */
     680  E = 1.0 / (1 + Z);   
     681
     682  q = par[5] * E;
     683  R = q*E;
     684  F = 0.5 / Z;
     685 
     686  f = q + par[6];
     687
     688  dpar[0] = F*R*(p2 + par[4]*Y);
     689  dpar[1] = F*R*(p3 + par[4]*X);
     690  dpar[2] = F*0.5*R*p2*p2;
     691  dpar[3] = F*0.5*R*p3*p3;
     692  dpar[4] = -R*X*Y*F;
     693   
     694  dpar[5] = E;
     695  dpar[6] = 1;
     696  return (f);
     697}
     698
     699/* pars: x, y, sx, sy, sxy, I, sky */
     700float fbarTD (float x, float y, float *par, int Npar, float *dpar) {
     701
     702  float X, Y, Z, E, F, q, R, f, p2, p3;
     703
     704  X = x - par[0];
     705  Y = y - par[1];
     706 
     707  p2 = X / par[2];
     708  p3 = Y / par[3];
     709
     710  Z = 0.5*p2*X + X*Y*par[4] + 0.5*p3*Y;                 /* R */
     711  E = 1.0 / (1 + Z*Z*Z);   
     712
     713  q = par[5] * E;
     714  F = 3*Z*Z;
     715  R = q*E*F;
     716 
     717  f = q + par[6];
     718
     719  dpar[0] = R*(p2 + par[4]*Y);
     720  dpar[1] = R*(p3 + par[4]*X);
     721  dpar[2] = 0.5*R*p2*p2;
     722  dpar[3] = 0.5*R*p3*p3;
     723  dpar[4] = -R*X*Y;
     724   
     725  dpar[5] = E;
     726  dpar[6] = 1;
     727  return (f);
     728}
     729
     730/* convert from x,y to major,minor */
     731void fix_ellipsegauss_pars (float *par, int Npar) {
     732
     733  float p2, p4, angle, t1, t2, tmp, area;
     734
     735  /* par[0], par[1] = Xo, Yo - stay the same */
     736
     737  p2 = 1/par[2];
     738  p4 = 1/par[3];
     739
     740  angle = 0.5 * atan2 (-2*par[4], p4 - p2);
     741
     742  tmp = sqrt (SQ(p2 - p4) + 4*SQ(par[4]));
     743  t1 = (p2 + p4 + tmp) / 2;
     744  t2 = t1 - tmp;
     745
     746  par[2] = 2.35482*sqrt(1/t2);
     747  par[3] = 2.35482*sqrt(1/t1);
     748  par[4] = DEG_RAD * angle;
     749
     750  area = 2*M_PI/sqrt(t1*t2);
     751
     752  par[5] *= area;
     753
     754}
     755# endif
     756
     757/***  options for later
     758
     759  Subtract = FALSE;
     760  if ((N = get_argument (argc, argv, "-sub"))) {
     761    remove_argument (N, &argc, argv);
     762    Subtract  = TRUE;
     763  }
     764
     765  DFact = 1;
     766  if ((N = get_argument (argc, argv, "-D"))) {
     767    remove_argument (N, &argc, argv);
     768    DFact  = atof(argv[N]);
     769    remove_argument (N, &argc, argv);
     770  }
     771
     772  fitfunc = fgaussTD;
     773  if ((N = get_argument (argc, argv, "-gal"))) {
     774    remove_argument (N, &argc, argv);
     775    fitfunc = fgalaxyTD;
     776  }
     777  if ((N = get_argument (argc, argv, "-bar"))) {
     778    remove_argument (N, &argc, argv);
     779    fitfunc = fbarTD;
     780  }
     781
     782
     783  f1 = 1;
     784  if ((c = get_variable ("BETA1")) != (char *) NULL) f1 = atof (c);
     785
     786  f2 = 1;
     787  if ((c = get_variable ("BETA2")) != (char *) NULL) f2 = atof (c);
     788
     789  if (Subtract) {
     790    tmpsky = par[6];
     791    par[6] = 0;
     792    for (N = j = 0; j < ny; j++) {
     793      V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx;
     794      for (i = 0; i < nx; i++, V++, N++) {
     795        dx = i + sx;
     796        dy = j + sy;
     797        *V -= fitfunc (dx, dy, par, Npar, (float *) NULL);
     798      }
     799    }
     800    par[6] = tmpsky;
     801  }
     802
     803***/
     804
     805# if (0)
     806
     807/* these two tests were not very succcessful.  the first did not model the shape well because
     808   it could not match the change in roundness with radius.  the second did not work because the
     809   parameters were degenerate (amplitude and slope of second component) */
     810
     811/* test: fixed, non-integer higher-order term -- x, y, sx, sy, sxy, I, sky, f1, f2 */
    325812float qgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    326813
    … …  
    354841}
    355842
    356 /* pars: x, y, sx, sy, sxy, I, sky */
     843/* test: two component model: inner pseudo gaussian with outer z^1.75 x, y, sx, sy, sxy, I, sky */
    357844float rgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    358845
    … …  
    394881}
    395882
    396 /* pars: x, y, sx, sy, sxy, I, sky */
    397 float sgaussTD (float x, float y, float *par, int Npar, float *dpar) {
    398 
    399   float X, Y, px1, py1, px2, py2;
    400   float z1, z2, r, q1, q2, f;
    401 
    402   X = x - par[0];
    403   Y = y - par[1];
    404  
    405   px1 = par[2]*X;
    406   py1 = par[3]*Y;
    407   px2 = par[7]*X;
    408   py2 = par[8]*Y;
    409 
    410   z1 = 0.5*SQ(px1) + 0.5*SQ(py1) + par[4]*X*Y;
    411   z2 = 0.5*SQ(px2) + 0.5*SQ(py2) + par[9]*X*Y;
    412 
    413   r = 1.0 / (1 + z1 + pow(z2,Npow));
    414   f = par[5]*r + par[6];
    415 
    416   q1 = par[5]*SQ(r);
    417   q2 = par[5]*SQ(r)*Npow*pow(z2,(Npow-1));
    418 
    419   dpar[0] = q1*(2*px1*par[2] + par[4]*Y) + q2*(2*px2*par[7] + par[9]*Y);
    420   dpar[1] = q1*(2*py1*par[3] + par[4]*X) + q2*(2*py2*par[8] + par[9]*X);
    421   dpar[2] = -2*q1*px1*X;
    422   dpar[3] = -2*q1*py1*Y;
    423   dpar[4] = -q1*X*Y;
    424   dpar[5] = +r;
    425   dpar[6] = +1;
    426   dpar[7] = -2*q2*px2*X;
    427   dpar[8] = -2*q2*py2*Y;
    428   dpar[9] = -q2*X*Y;
    429 
    430   return (f);
    431 }
    432 
    433 /* alternate for exp(-z) */
    434 /* E = 1.0 / (1 + f1*Z + f2*(0.5*Z*Z + 0.1666666*Z*Z*Z)); */  /* ~ exp (-Z) */
    435 /* E = 1.0 / (1 +    Z +    (0.5*Z*Z + 0.1666666*Z*Z*Z)); */  /* ~ exp (-Z) */
    436 
    437 # if (0)
    438 /* pars: x, y, sx, sy, sxy, I, sky */
    439 float fgalaxyTD (float x, float y, float *par, int Npar, float *dpar) {
    440 
    441   float X, Y, Z, E, F, q, R, f, p2, p3;
    442 
    443   X = x - par[0];
    444   Y = y - par[1];
    445  
    446   p2 = X / par[2];
    447   p3 = Y / par[3];
    448 
    449   Z = sqrt (0.5*p2*X + X*Y*par[4] + 0.5*p3*Y);                 /* R */
    450   E = 1.0 / (1 + Z);   
    451 
    452   q = par[5] * E;
    453   R = q*E;
    454   F = 0.5 / Z;
    455  
    456   f = q + par[6];
    457 
    458   dpar[0] = F*R*(p2 + par[4]*Y);
    459   dpar[1] = F*R*(p3 + par[4]*X);
    460   dpar[2] = F*0.5*R*p2*p2;
    461   dpar[3] = F*0.5*R*p3*p3;
    462   dpar[4] = -R*X*Y*F;
    463    
    464   dpar[5] = E;
    465   dpar[6] = 1;
    466   return (f);
    467 }
    468 
    469 /* pars: x, y, sx, sy, sxy, I, sky */
    470 float fbarTD (float x, float y, float *par, int Npar, float *dpar) {
    471 
    472   float X, Y, Z, E, F, q, R, f, p2, p3;
    473 
    474   X = x - par[0];
    475   Y = y - par[1];
    476  
    477   p2 = X / par[2];
    478   p3 = Y / par[3];
    479 
    480   Z = 0.5*p2*X + X*Y*par[4] + 0.5*p3*Y;                 /* R */
    481   E = 1.0 / (1 + Z*Z*Z);   
    482 
    483   q = par[5] * E;
    484   F = 3*Z*Z;
    485   R = q*E*F;
    486  
    487   f = q + par[6];
    488 
    489   dpar[0] = R*(p2 + par[4]*Y);
    490   dpar[1] = R*(p3 + par[4]*X);
    491   dpar[2] = 0.5*R*p2*p2;
    492   dpar[3] = 0.5*R*p3*p3;
    493   dpar[4] = -R*X*Y;
    494    
    495   dpar[5] = E;
    496   dpar[6] = 1;
    497   return (f);
    498 }
    499 
    500 /* convert from x,y to major,minor */
    501 void fix_ellipsegauss_pars (float *par, int Npar) {
    502 
    503   float p2, p4, angle, t1, t2, tmp, area;
    504 
    505   /* par[0], par[1] = Xo, Yo - stay the same */
    506 
    507   p2 = 1/par[2];
    508   p4 = 1/par[3];
    509 
    510   angle = 0.5 * atan2 (-2*par[4], p4 - p2);
    511 
    512   tmp = sqrt (SQ(p2 - p4) + 4*SQ(par[4]));
    513   t1 = (p2 + p4 + tmp) / 2;
    514   t2 = t1 - tmp;
    515 
    516   par[2] = 2.35482*sqrt(1/t2);
    517   par[3] = 2.35482*sqrt(1/t1);
    518   par[4] = DEG_RAD * angle;
    519 
    520   area = 2*M_PI/sqrt(t1*t2);
    521 
    522   par[5] *= area;
    523 
    524 }
    525883# endif
    526884
    527 /***  options for later
    528 
    529   Subtract = FALSE;
    530   if ((N = get_argument (argc, argv, "-sub"))) {
    531     remove_argument (N, &argc, argv);
    532     Subtract  = TRUE;
    533   }
    534 
    535   DFact = 1;
    536   if ((N = get_argument (argc, argv, "-D"))) {
    537     remove_argument (N, &argc, argv);
    538     DFact  = atof(argv[N]);
    539     remove_argument (N, &argc, argv);
    540   }
    541 
    542   fitfunc = fgaussTD;
    543   if ((N = get_argument (argc, argv, "-gal"))) {
    544     remove_argument (N, &argc, argv);
    545     fitfunc = fgalaxyTD;
    546   }
    547   if ((N = get_argument (argc, argv, "-bar"))) {
    548     remove_argument (N, &argc, argv);
    549     fitfunc = fbarTD;
    550   }
    551 
    552 
    553   f1 = 1;
    554   if ((c = get_variable ("BETA1")) != (char *) NULL) f1 = atof (c);
    555 
    556   f2 = 1;
    557   if ((c = get_variable ("BETA2")) != (char *) NULL) f2 = atof (c);
    558 
    559   if (Subtract) {
    560     tmpsky = par[6];
    561     par[6] = 0;
    562     for (N = j = 0; j < ny; j++) {
    563       V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx;
    564       for (i = 0; i < nx; i++, V++, N++) {
    565         dx = i + sx;
    566         dy = j + sy;
    567         *V -= fitfunc (dx, dy, par, Npar, (float *) NULL);
    568       }
    569     }
    570     par[6] = tmpsky;
    571   }
    572 
    573 ***/
     885  /* forcing values to have a rational range
     886  ALLOCATE (parmin, float, Npar);
     887  ALLOCATE (parmax, float, Npar);
     888  bzero (parmin, Npar*sizeof(float));
     889  bzero (parmax, Npar*sizeof(float));
     890  parmin[0] = parmin[1] = 0;
     891  parmax[0] = buf[0].matrix.Naxis[0];
     892  parmax[1] = buf[0].matrix.Naxis[1];
     893
     894  parmin[2] = parmin[3] = 0.01;
     895  parmax[2] = parmax[3] = 100.0;
     896  parmin[4] = -1000;
     897  parmax[4] = -1000;
     898 
     899  parmin[5] = 1;
     900  parmax[5] = 1e5;
     901
     902  parmin[6] = 0.0;
     903  parmax[6] = 1e5;
     904
     905  if (Npar == 9) {
     906    parmin[7] = parmin[8] = 0.01;
     907    parmax[7] = parmax[8] = 10.0;
     908  }
     909  if (Npar == 10) {
     910    parmin[7] = parmin[8] = 0.01;
     911    parmax[7] = parmax[8] = 10.0;
     912    parmin[9] = -1000;
     913    parmax[9] = -1000;
     914  }
     915  */
     916

  • trunk/Ohana/src/opihi/cmd.data/gridify.c

    r2843 r2938  
    33int gridify (int argc, char **argv) {
    44
    5   int i, Nx, Ny, Xb, Yb;
     5  int i, Nx, Ny, Xb, Yb, Normalize, N;
    66  float Xmin, Xmax, dX, Ymin, Ymax, dY;
    77  float *buf, *val, *x, *y, *z;
    … …  
    99  Buffer *bf;
    1010  Vector *vx, *vy, *vz;
     11
     12  Normalize = TRUE;
     13  if ((N = get_argument (argc, argv, "-raw"))) {
     14    remove_argument (N, &argc, argv);
     15    Normalize = FALSE;
     16  }
    1117
    1218  if (argc != 11) {
    … …  
    6672    if (Xb >= Nx) continue;
    6773    if (Yb >= Ny) continue;
    68     val[Xb + Yb*Nx] = *z;
     74    val[Xb + Yb*Nx] += *z;
    6975    Nval[Xb + Yb*Nx]++;
    7076  }
    … …  
    7278  buf = (float *) bf[0].matrix.buffer;
    7379  for (i = 0; i < Nx*Ny; i++) {
    74     if (Nval[i] == 0) {
    75       buf[i] = 0;
    76       continue;
     80    if (Normalize) {
     81      if (Nval[i] == 0) {
     82        buf[i] = 0;
     83        continue;
     84      }
     85      buf[i] = val[i] / Nval[i];
     86    } else {
     87      buf[i] = val[i];
    7788    }
    78     buf[i] = val[i] / Nval[i];
    7989  }
    8090

  • trunk/Ohana/src/opihi/cmd.data/imsmooth.c

    r2843 r2938  
    11# include "data.h"
    2 # define NSIGMA 3
    32
    43int imsmooth (int argc, char **argv) {
    54 
    6   int i, j, n, Nx, Ny, Ns, Ngauss;
     5  int i, j, n, N, Nx, Ny, Ns, Ngauss;
    76  float *vi, *vo, *gauss, *gaussnorm;
    8   float g, s, sigma;
     7  float g, s, sigma, Nsigma;
    98  Buffer *in;
    109  float *temp;
     10
     11  Nsigma = 3;
     12  if ((N = get_argument (argc, argv, "-Nsigma"))) {
     13    remove_argument (N, &argc, argv);
     14    Nsigma = atof (argv[N]);
     15    remove_argument (N, &argc, argv);
     16  }
    1117
    1218  if (argc != 3) {
    … …  
    2329
    2430  /* build a 1D gaussian */
    25   Ns = (int) (NSIGMA*sigma + 0.5);
     31  Ns = (int) (Nsigma*sigma + 0.5);
    2632  Ngauss = 2*Ns + 1;
    2733  ALLOCATE (gaussnorm, float, Ngauss);

  • trunk/Ohana/src/opihi/lib.data/mrq2dmin.c

    r2843 r2938  
    77*/
    88
     9# define VERY_VERBOSE 0
     10
    911static float **alpha, **talpha;
    1012static float **beta, **tbeta;
    1113static float *partry, *dyda;
    1214static float ochisq, lambda;
     15
     16static float *parmin = NULL;
     17static float *parmax = NULL;
    1318
    1419float mrq2dcof (float *x, float *t, float *y, float *dy, int Npts,
    … …  
    4146      ta[k][j] = ta[j][k];
    4247     
     48# if (VERY_VERBOSE)
     49  for (j = 0; j < Npar; j++) {
     50    for (k = 0; k < Npar; k++) {
     51      fprintf (stderr, "%9.3e  ", ta[j][k]);
     52    }
     53    fprintf (stderr, "    :   %9.3e  ", tb[j][0]);
     54    fprintf (stderr, "\n");
     55  }
     56# endif
     57
    4358  return (chisq);
    4459
     60}
     61
     62float mrq2dchi (float *x, float *t, float *y, float *dy, int Npts,
     63                float *par, int Npar,
     64                float (funcs)(float, float, float *, int, float *)) {
     65
     66  int k, j, i;
     67  float ydiff, wt, chisq;
     68
     69  chisq = 0.0;
     70  for (i = 0; i < Npts; i++) {
     71    ydiff = funcs (x[i], t[i], par, Npar, dyda) - y[i];
     72    chisq += SQ(ydiff) * dy[i];
     73  }
     74  return (chisq);
    4575}
    4676
    … …  
    6595  }
    6696
    67   for (j = 0; j < Npar; j++) partry[j] = par[j] - tbeta[j][0];
     97  for (j = 0; j < Npar; j++) {
     98    partry[j] = par[j] - tbeta[j][0];
     99    /*
     100    if (parmin != NULL) partry[j] = MAX (parmin[j], partry[j]);
     101    if (parmax != NULL) partry[j] = MIN (parmax[j], partry[j]);
     102    */
     103  }
    68104
    69105  chisq = mrq2dcof (x, t, y, dy, Npts, partry, Npar, talpha, tbeta, funcs);
    … …  
    93129  return (chisq);
    94130
     131}
     132
     133int mrq2dlimits (float *pmin, float *pmax, int Npar) {
     134
     135  int i;
     136
     137  ALLOCATE (parmin, float, Npar);
     138  ALLOCATE (parmax, float, Npar);
     139  for (i = 0; i < Npar; i++) {
     140    parmin[i] = pmin[i];
     141    parmax[i] = pmax[i];
     142  }
     143  return (TRUE);
    95144}
    96145
    … …  
    156205  free (dyda);
    157206
    158 }
     207  if (parmin != NULL) free (parmin);
     208  if (parmax != NULL) free (parmax);
     209
     210}

  • trunk/Ohana/src/opihi/mana/findpeaks.c

    r2826 r2938  
    77  int Npeak, NPEAK, Npeaks;
    88  float *v;
    9   int *peaks, *keep, *xp, *yp;
     9  int *peaks, *keep, *xp, *yp, *zp;
    1010  float threshold, vt, vo;
    11   Vector *vecx, *vecy;
     11  Vector *vecx, *vecy, *vecz;
    1212  Buffer *buf;
    1313
    … …  
    1919  if ((vecx = SelectVector ("xp", ANYVECTOR, TRUE)) == NULL) return (FALSE);
    2020  if ((vecy = SelectVector ("yp", ANYVECTOR, TRUE)) == NULL) return (FALSE);
     21  if ((vecz = SelectVector ("zp", ANYVECTOR, TRUE)) == NULL) return (FALSE);
    2122
    2223  Nx = buf[0].matrix.Naxis[0];
    … …  
    2728  ALLOCATE (xp, int, NPEAK);
    2829  ALLOCATE (yp, int, NPEAK);
     30  ALLOCATE (zp, int, NPEAK);
    2931
    3032  /* find peaks for each row */
    … …  
    3739      REALLOCATE (xp, int, NPEAK);
    3840      REALLOCATE (yp, int, NPEAK);
     41      REALLOCATE (zp, int, NPEAK);
    3942    }
    4043    for (i = 0; i < Npeaks; i++) {
    4144      xp[Npeak + i] = peaks[i];
    4245      yp[Npeak + i] = j;
     46      zp[Npeak + i] = v[peaks[i] + j*Nx];
    4347    }
    4448    Npeak += Npeaks;
    … …  
    7478  REALLOCATE (vecx[0].elements, float, MAX (Npeak, 1));
    7579  REALLOCATE (vecy[0].elements, float, MAX (Npeak, 1));
     80  REALLOCATE (vecz[0].elements, float, MAX (Npeak, 1));
    7681  /* eliminate non-local peaks */
    7782  for (N = n = 0; n < Npeak; n++) {
    … …  
    7984    vecx[0].elements[N] = xp[n];
    8085    vecy[0].elements[N] = yp[n];
     86    vecz[0].elements[N] = zp[n];
    8187    N ++;
    8288  }
     89  free (xp);
     90  free (yp);
     91  free (zp);
     92  free (keep);
    8393
    8494  REALLOCATE (vecx[0].elements, float, MAX (N, 1));
    8595  REALLOCATE (vecy[0].elements, float, MAX (N, 1));
    86   vecx[0].Nelements = vecy[0].Nelements = N;
     96  REALLOCATE (vecz[0].elements, float, MAX (N, 1));
     97  vecx[0].Nelements = vecy[0].Nelements = vecz[0].Nelements = N;
    8798
    8899  return (TRUE);
    … …  
    93104}
    94105
     106

  • trunk/Ohana/src/opihi/mana/findrowpeaks.c

    r2826 r2938  
    33int *findrowpeaks (float *row, int Nrow, float threshold, int *npeaks) {
    44 
    5   int i, rising;
     5  int i;
    66  int Npeaks, NPEAKS;
    77  int *peaks;
    8 
    9 # if (0)
    10   ALLOCATE (peaks, int, 3);
    11   peaks[0] = 20;
    12   peaks[1] = 30;
    13   peaks[2] = 40;
    14   *npeaks = 3;
    15   return (peaks);
    16 # endif
    178
    189  Npeaks = 0;
    … …  
    2516  }
    2617   
    27   rising = (row[1] > row[0]);
     18  /* special case for first pixel in row */
     19  if ((row[0] > row[1]) && (row[0] > threshold)) {
     20    peaks[Npeaks] = 0;
     21    Npeaks ++;
     22  }   
    2823
    29   for (i = 2; i < Nrow; i++) {
    30     if (row[i] < threshold) {
    31       rising = TRUE;
    32       continue;
    33     }
    34     if (rising) {
    35       if (row[i] < row[i-1]) {
    36         rising = FALSE;
    37         peaks[Npeaks] = i-1;
    38         Npeaks ++;
    39         if (Npeaks >= NPEAKS) {
    40           NPEAKS += 100;
    41           REALLOCATE (peaks, int, NPEAKS);
    42           fprintf (stderr, "#");
    43         }
    44       }
    45       continue;
    46     }
    47     if (!rising) {
    48       if (row[i] > row[i-1]) {
    49         rising = TRUE;
    50       }
    51       continue;
     24  for (i = 1; i < Nrow - 1; i++) {
     25    if (row[i] < threshold) continue;
     26    if (row[i] < row[i-1]) continue;
     27    if (row[i] <= row[i+1]) continue;
     28
     29    peaks[Npeaks] = i;
     30    Npeaks ++;
     31    if (Npeaks >= NPEAKS) {
     32      NPEAKS += 100;
     33      REALLOCATE (peaks, int, NPEAKS);
    5234    }
    5335  }     
     36
     37  /* special case for last pixel in row */
     38  if ((row[Nrow-1] >= row[Nrow-2]) && (row[Nrow-1] > threshold)) {
     39    peaks[Npeaks] = Nrow-1;
     40    Npeaks ++;
     41  }   
     42
    5443  *npeaks = Npeaks;
    5544  return (peaks);

  • trunk/Ohana/src/opihi/mana/rawstars.c

    r2843 r2938  
    77  double Raper, Rinner, Router;
    88  Vector *xp, *yp;
    9   Vector *xc, *yc, *sx, *sy, *sxy, *zs, *zp, *sk;
     9  Vector *xc, *yc, *sx, *sy, *sxy, *zs, *zc, *sk;
    1010  Buffer *buff;
    1111
    … …  
    4444  if ((xc = SelectVector ("xc", ANYVECTOR, TRUE)) == NULL) return (FALSE);
    4545  if ((yc = SelectVector ("yc", ANYVECTOR, TRUE)) == NULL) return (FALSE);
     46  if ((zc = SelectVector ("zc", ANYVECTOR, TRUE)) == NULL) return (FALSE);
     47  if ((zs = SelectVector ("zs", ANYVECTOR, TRUE)) == NULL) return (FALSE);
     48  if ((sk = SelectVector ("sk", ANYVECTOR, TRUE)) == NULL) return (FALSE);
    4649  if ((sx = SelectVector ("sx", ANYVECTOR, TRUE)) == NULL) return (FALSE);
    4750  if ((sy = SelectVector ("sy", ANYVECTOR, TRUE)) == NULL) return (FALSE);
    4851  if ((sxy = SelectVector ("sxy", ANYVECTOR, TRUE)) == NULL) return (FALSE);
    49   if ((zs = SelectVector ("zs", ANYVECTOR, TRUE)) == NULL) return (FALSE);
    50   if ((zp = SelectVector ("zp", ANYVECTOR, TRUE)) == NULL) return (FALSE);
    51   if ((sk = SelectVector ("sk", ANYVECTOR, TRUE)) == NULL) return (FALSE);
    5252
    5353  Nx = buff[0].matrix.Naxis[0];
    … …  
    6161  REALLOCATE (sxy[0].elements, float, Np);
    6262  REALLOCATE (zs[0].elements, float, Np);
    63   REALLOCATE (zp[0].elements, float, Np);
     63  REALLOCATE (zc[0].elements, float, Np);
    6464  REALLOCATE (sk[0].elements, float, Np);
    6565  xc[0].Nelements = yc[0].Nelements = sx[0].Nelements = Np;
    66   sy[0].Nelements = zs[0].Nelements = zp[0].Nelements = Np;
     66  sy[0].Nelements = zs[0].Nelements = zc[0].Nelements = Np;
    6767  sxy[0].Nelements = sk[0].Nelements = Np;
    6868
    … …  
    8383                    &sxy[0].elements[i],
    8484                    &zs[0].elements[i],
    85                     &zp[0].elements[i],
     85                    &zc[0].elements[i],
    8686                    &sk[0].elements[i]);
    8787  }
Note: See TracChangeset for help on using the changeset viewer.