IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 27, 2019, 11:39:38 AM (7 years ago)
Author:
eugene
Message:

merge changes from trunk (mostly Ohana & ippMonitor improvements, some work on ippTools)

Location:
branches/eam_branches/ipp-20191011
Files:
1 deleted
97 edited
37 copied

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20191011

  • branches/eam_branches/ipp-20191011/Ohana

  • branches/eam_branches/ipp-20191011/Ohana/src/addstar/Makefile

    r40543 r41170  
    2727loadgaia     : $(BIN)/loadgaia.$(ARCH)
    2828loadgaia_dr2 : $(BIN)/loadgaia_dr2.$(ARCH)
     29loadukirt_uhs : $(BIN)/loadukirt_uhs.$(ARCH)
    2930loadstarpar  : $(BIN)/loadstarpar.$(ARCH)
    3031loadstarpar_client : $(BIN)/loadstarpar_client.$(ARCH)
    … …  
    4445# programs in 'SERVER' use the client-server concept and are out of date
    4546
    46 INSTALL = addstar addstar_client sedstar loadgalphot loadgaia loadgaia_dr2 loadstarpar loadstarpar_client setobjflags setobjflags_client loadICRF loadICRF_client skycells mkcmf dumpskycells findskycell load2mass loadwise loadtycho loadbsc loadsupercos
     47INSTALL = addstar addstar_client sedstar loadgalphot loadgaia loadgaia_dr2 loadukirt_uhs loadstarpar loadstarpar_client setobjflags setobjflags_client loadICRF loadICRF_client skycells mkcmf dumpskycells findskycell load2mass loadwise loadtycho loadbsc loadsupercos
    4748SERVER  = addstarc addstard addstart
    4849
    … …  
    305306$(SRC)/psps_ids.$(ARCH).o
    306307
     308LOAD-UKIRT_UHS = \
     309$(SRC)/loadukirt_uhs.$(ARCH).o \
     310$(SRC)/ConfigInit.$(ARCH).o \
     311$(SRC)/SetSignals.$(ARCH).o \
     312$(SRC)/Shutdown.$(ARCH).o \
     313$(SRC)/args_loadukirt_uhs.$(ARCH).o \
     314$(SRC)/find_matches_ukirt_uhs.$(ARCH).o \
     315$(SRC)/loadukirt_uhs_catalog.$(ARCH).o \
     316$(SRC)/loadukirt_uhs_make_subset.$(ARCH).o \
     317$(SRC)/loadukirt_uhs_readstars.$(ARCH).o \
     318$(SRC)/loadukirt_uhs_readstars_uhs.$(ARCH).o \
     319$(SRC)/loadukirt_uhs_readstars_ugcs.$(ARCH).o \
     320$(SRC)/loadukirt_uhs_readstars_ugps.$(ARCH).o \
     321$(SRC)/loadukirt_uhs_readstars_ulas.$(ARCH).o \
     322$(SRC)/loadukirt_uhs_table.$(ARCH).o \
     323$(SRC)/parse_csv.$(ARCH).o \
     324$(SRC)/resort_catalog.$(ARCH).o \
     325$(SRC)/build_links.$(ARCH).o \
     326$(SRC)/strhash.$(ARCH).o \
     327$(SRC)/sortIDs.$(ARCH).o \
     328$(SRC)/psps_ids.$(ARCH).o
     329
    307330LOAD-STARPAR = \
    308331$(SRC)/loadstarpar.$(ARCH).o \
    … …  
    528551$(LOAD-GAIA)           : $(INC)/addstar.h $(INC)/gaia.h
    529552$(LOAD-GAIA_DR2)       : $(INC)/addstar.h $(INC)/gaia_dr2.h
     553$(LOAD-UKIRT_UHS)      : $(INC)/addstar.h $(INC)/ukirt_uhs.h
    530554$(LOAD-STARPAR)        : $(INC)/addstar.h $(INC)/loadstarpar.h
    531555$(LOAD-STARPAR-CLIENT) : $(INC)/addstar.h $(INC)/loadstarpar.h
    … …  
    551575$(BIN)/loadgaia.$(ARCH)       : $(LOAD-GAIA)
    552576$(BIN)/loadgaia_dr2.$(ARCH)   : $(LOAD-GAIA_DR2)
     577$(BIN)/loadukirt_uhs.$(ARCH)  : $(LOAD-UKIRT_UHS)
    553578$(BIN)/loadstarpar.$(ARCH)    : $(LOAD-STARPAR)
    554579$(BIN)/loadstarpar_client.$(ARCH) : $(LOAD-STARPAR-CLIENT)

  • branches/eam_branches/ipp-20191011/Ohana/src/addstar/src/ReadImageHeader.c

    r40291 r41170  
    118118  if (photcode == 0) {
    119119    if (!gfits_scan (header, "PHOTCODE", "%s", 1, photname)) {
    120       fprintf (stderr, "photcode not supplied in header\n");
     120      fprintf (stderr, "failure: photcode not supplied in header\n");
    121121      return (FALSE);
    122122    }

  • branches/eam_branches/ipp-20191011/Ohana/src/dvomerge/include/dvomerge.h

    r40805 r41170  
    3131int    IMAGES_ONLY;
    3232int    ACCEPT_MOTION;
     33
     34int    ACCEPT_ASTROM;
    3335int    RETAIN_AVE_PHOTOMETRY;
     36
    3437char   CATDIR[DVO_MAX_PATH];
    3538char   GSCFILE[DVO_MAX_PATH];

  • branches/eam_branches/ipp-20191011/Ohana/src/dvomerge/src/args.c

    r40805 r41170  
    112112  }
    113113
    114   /* limit the impact of a dvomerge -parallel */
     114  /* accept input database average astrometry motions */
    115115  ACCEPT_MOTION = FALSE;
    116116  if ((N = get_argument (*argc, argv, "-accept-motion"))) {
    117117    remove_argument (N, argc, argv);
    118118    ACCEPT_MOTION = TRUE;
     119  }
     120  /* accept input database average astrometry information */
     121  ACCEPT_ASTROM = FALSE;
     122  if ((N = get_argument (*argc, argv, "-accept-astrom"))) {
     123    remove_argument (N, argc, argv);
     124    ACCEPT_ASTROM = TRUE;
    119125  }
    120126
    … …  
    337343  }
    338344
    339   /* limit the impact of a dvomerge -parallel */
     345  /* accept input database average astrometry motions */
    340346  ACCEPT_MOTION = FALSE;
    341347  if ((N = get_argument (*argc, argv, "-accept-motion"))) {
    342348    remove_argument (N, argc, argv);
    343349    ACCEPT_MOTION = TRUE;
     350  }
     351  /* accept input database average astrometry information */
     352  ACCEPT_ASTROM = FALSE;
     353  if ((N = get_argument (*argc, argv, "-accept-astrom"))) {
     354    remove_argument (N, argc, argv);
     355    ACCEPT_ASTROM = TRUE;
    344356  }
    345357

  • branches/eam_branches/ipp-20191011/Ohana/src/dvomerge/src/dvomergeUpdate_catalogs.c

    r40805 r41170  
    328328    if (FORCE_MERGE)           { strextend (&command, "-force-merge"); }
    329329    if (ACCEPT_MOTION)         { strextend (&command, "-accept-motion"); }
     330    if (ACCEPT_ASTROM)         { strextend (&command, "-accept-astrom"); }
    330331    if (RETAIN_AVE_PHOTOMETRY) { strextend (&command, "-retain-ave-photometry"); }
    331332    if (MATCHED_TABLES)        { strextend (&command, "-matched-tables"); }

  • branches/eam_branches/ipp-20191011/Ohana/src/dvomerge/src/merge_catalogs_old.c

    r40805 r41170  
    363363    }
    364364
     365    // XXX: add choice of secfilt
    365366    // update the average properties to reflect the incoming entries:
    366367    // if RETAIN_AVE_PHOTOMETRY is true and the original value is NAN, but the input value is not, accept the input:
    … …  
    389390
    390391    // we can choose to accept the proper-motion and parallax from the reference tgtcat
    391     if (ACCEPT_MOTION) {
     392    if (ACCEPT_MOTION || ACCEPT_ASTROM) {
    392393      output[0].average[n].dR         = input[0].average[N].dR;
    393394      output[0].average[n].dD         = input[0].average[N].dD;
    … …  
    400401      output[0].average[n].Tmean      = input[0].average[N].Tmean;
    401402    }
     403    if (ACCEPT_ASTROM) {
     404      output[0].average[n].R          = input[0].average[N].R;
     405      output[0].average[n].D          = input[0].average[N].D;
     406    }
    402407
    403408    /* Nm is updated, but not written out in -update mode (for existing entries)
    … …  
    461466
    462467    // we can choose to accept the proper-motion and parallax from the reference tgtcat
    463     if (ACCEPT_MOTION) {
     468    if (ACCEPT_MOTION || ACCEPT_ASTROM) {
    464469      output[0].average[Nave].dR         = input[0].average[N].dR;
    465470      output[0].average[Nave].dD         = input[0].average[N].dD;
    … …  
    471476      output[0].average[Nave].dP         = input[0].average[N].dP;
    472477      output[0].average[Nave].Tmean      = input[0].average[N].Tmean;
     478    }
     479    if (ACCEPT_ASTROM) {
     480      output[0].average[Nave].R          = input[0].average[N].R;
     481      output[0].average[Nave].D          = input[0].average[N].D;
    473482    }
    474483

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/Makefile

    r39401 r41170  
    4242$(SRC)/bDrawObjects.$(ARCH).o             $(SRC)/bDrawFrame.$(ARCH).o         \
    4343$(SRC)/bDrawLabels.$(ARCH).o              $(SRC)/bDrawIt.$(ARCH).o            \
    44 $(SRC)/bDrawImage.$(ARCH).o               \
     44$(SRC)/bDrawImage.$(ARCH).o               $(SRC)/MemoryDump.$(ARCH).o         \
    4545$(SRC)/PNGit.$(ARCH).o                    $(SRC)/PPMit.$(ARCH).o              \
    4646$(SRC)/PSit.$(ARCH).o                     $(SRC)/CrossHairs.$(ARCH).o         \

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/include/prototypes.h

    r40165 r41170  
    33int           args                PROTO((int *argc, char **argv));
    44void          SetUpGraphic        PROTO((int *argc, char **argv));
     5void          FreeGraphic         PROTO((void));
    56void          DefineLayout        PROTO((int, char **));
    67int           EventLoop           PROTO((void));
    78void          CloseDisplay        PROTO((void));
     9int           MemoryDump          PROTO((int sock));
     10int           MemoryDumpLines     PROTO((int sock));
     11int           MemoryDumpOnExit    PROTO((int sock));
     12int           MemoryDumpSetOnExit PROTO((int state));
     13int           MemoryDumpAndExit   PROTO((void));
    814
    915/* SetUpGraphic */
    … …  
    3339void          DrawTextlines       PROTO((KapaGraphWidget *graph));
    3440void          DrawConnect         PROTO((Graphic *graphic, KapaGraphWidget *graph, Gobjects *objects));
     41void          DrawPolygon         PROTO((Graphic *graphic, KapaGraphWidget *graph, Gobjects *objects));
     42void          DrawPolyfill        PROTO((Graphic *graphic, KapaGraphWidget *graph, Gobjects *objects));
    3543void          DrawHistogram       PROTO((Graphic *graphic, KapaGraphWidget *graph, Gobjects *objects));
    3644int           DrawObjectN         PROTO((Graphic *graphic, KapaGraphWidget *graph, Gobjects *objects));
    … …  
    8593int           SetColormapFromPipe PROTO((int sock));
    8694int           SetNanColorFromPipe PROTO((int sock));
     95int           SetSmoothSigma      PROTO((int sock));
    8796
    8897int           LoadVectorData      PROTO((int sock, KapaGraphWidget *graph, int N, char *type));
    … …  
    109118
    110119void          InitLayout          PROTO((int argc, char **argv));
     120void          FreeLayout          PROTO((void));
    111121
    112122/* PS drawing utilities */
    … …  
    142152void          bDrawBars           PROTO((bDrawBuffer *buffer, KapaGraphWidget *graph, Gobjects *object, int mode));
    143153void          bDrawPoints         PROTO((bDrawBuffer *buffer, KapaGraphWidget *graph, Gobjects *object));
     154void          bDrawPolygons       PROTO((bDrawBuffer *buffer, KapaGraphWidget *graph, Gobjects *object));
     155void          bDrawFillPolygons   PROTO((bDrawBuffer *buffer, KapaGraphWidget *graph, Gobjects *object));
    144156void          bDrawXErrors        PROTO((bDrawBuffer *buffer, KapaGraphWidget *graph, Gobjects *object));
    145157void          bDrawYErrors        PROTO((bDrawBuffer *buffer, KapaGraphWidget *graph, Gobjects *object));
    … …  
    179191void          InitButtonFunc      PROTO((Button *button, int (*function)(Graphic *graphic, KapaImageWidget *image)));
    180192void          DrawImage           PROTO((KapaImageWidget *image));
     193void          DrawImageTool       PROTO((KapaImageWidget *image));
    181194void          DrawButton          PROTO((Graphic *graphic, Button *button));
    182195void          DrawBitmap          PROTO((Graphic *graphic, int x, int y, int dx, int dy, unsigned char *bitmap, int mode));

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/include/structures.h

    r40570 r41170  
    4848  unsigned long  back;        // basic background color
    4949
     50  float smooth_sigma; // anti-aliasing smoothing scale
     51
    5052  unsigned long  overlay_color[NOVERLAYS]; // image plotting colors
    5153} Graphic;
    … …  
    131133} Overlay;
    132134 
    133 typedef struct {
    134   float *x, *y, *z, *dxp, *dxm, *dyp, *dym;
    135   int Npts;
    136   int style, ptype, ltype, color, etype, ebar;
    137   double lweight, size;
     135// a set of objects all have the same basic properties
     136typedef struct {
     137  float *x; // x-coordinates of the points
     138  float *y; // y-coordinates of the points
     139  float *z; // size/colorscale for points
     140  float *dxp; // lower-errorbar in x
     141  float *dxm; // upper-errorbar in x
     142  float *dyp; // lower-errorbar in y
     143  float *dym; // upper-errorbar in y
     144  int Npts;   // number of points in this set
     145  int style;  // how are the object draw: CONNECT, HISTOGRAM,
     146  int ptype;  // shape of object at each point
     147  // ptype is overloaded for NPOLYGON to be the number of points / polygon
     148  int ltype;  // style of line (solid, dot, dash, etc)
     149  int color;  // color for point (if not colorscaled)
     150  int etype;  // errorbars to draw (0x01 = y, 0x02 = x)
     151  int ebar;   // draw a cap on the error bar
     152  double lweight; // line thickness
     153  double size; // size of the object
    138154  double x0, x1, y0, y1;  /* limits for this object */
    139155  double alpha;

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/CheckPipe.c

    r38465 r41170  
    352352  }
    353353
     354  if (!strcmp (word, "SIGM")) {
     355    SetSmoothSigma (sock);
     356    KiiSendCommand (sock, 4, "DONE");
     357    FINISHED (TRUE);
     358  }
     359
     360  if (!strcmp (word, "MEMD")) {
     361    MemoryDump (sock);
     362    KiiSendCommand (sock, 4, "DONE");
     363    FINISHED (TRUE);
     364  }
     365
     366  if (!strcmp (word, "MEML")) {
     367    MemoryDumpLines (sock);
     368    KiiSendCommand (sock, 4, "DONE");
     369    FINISHED (TRUE);
     370  }
     371
     372  if (!strcmp (word, "MEMX")) {
     373    MemoryDumpOnExit (sock);
     374    KiiSendCommand (sock, 4, "DONE");
     375    FINISHED (TRUE);
     376  }
     377
    354378  fprintf (stderr, "unknown signal %s\n", word);
    355379  KiiSendCommand (sock, 4, "DONE");

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/DrawObjects.c

    r40291 r41170  
    9797  switch (object[0].style) {
    9898    case KAPA_PLOT_CONNECT:
    99       // fprintf (stderr, "plot KAPA_PLOT_CONNECT: \n");
    10099      DrawConnect (graphic, graph, object);
    101100      break;
     101    case KAPA_PLOT_POLYGON:
     102      DrawPolygon (graphic, graph, object);
     103      break;
     104    case KAPA_PLOT_POLYFILL:
     105      DrawPolyfill (graphic, graph, object);
     106      break;
    102107    case KAPA_PLOT_HISTOGRAM:
    103       // fprintf (stderr, "plot KAPA_PLOT_HISTOGRAM:\n");
    104108      DrawHistogram (graphic, graph, object);
    105109      break;
    106110    case KAPA_PLOT_BARS_SOLID:
    107       // fprintf (stderr, "plot KAPA_PLOT_BARS_SOLID:\n");
    108111      DrawBars (graphic, graph, object, KAPA_PLOT_BARS_SOLID);
    109112      break;
    110113    case KAPA_PLOT_BARS_OUTLINE:
    111       // fprintf (stderr, "plot KAPA_PLOT_BARS_OUTLINE:\n");
    112114      DrawBars (graphic, graph, object, KAPA_PLOT_BARS_OUTLINE);
    113115      break;
    114116    case KAPA_PLOT_BARS_OUTFILL:
    115       // fprintf (stderr, "plot KAPA_PLOT_BARS_OUTFILL:\n");
    116117      DrawBars (graphic, graph, object, KAPA_PLOT_BARS_OUTFILL);
    117118      break;
    118119    case KAPA_PLOT_POINTS:
    119120    default:
    120       // fprintf (stderr, "plot KAPA_PLOT_POINTS:\n");
    121121      DrawPoints (graphic, graph, object);
    122122      break;
    … …  
    172172    sx0 = sx1; sy0 = sy1;
    173173  }
    174  
     174}
     175
     176/******/
     177void DrawPolygon (Graphic *graphic, KapaGraphWidget *graph, Gobjects *object) {
     178
     179  float *x, *y;
     180  double mxi, mxj, myi, myj, bxi, bxj, byi, byj, bx, by;
     181  double sx0, sy0, sx1, sy1;
     182  double X0, X1, Y0, Y1;
     183
     184  mxi = graph[0].axis[0].dfx / (object[0].x1 - object[0].x0);
     185  mxj = graph[0].axis[1].dfx / (object[0].y1 - object[0].y0);
     186  myi = graph[0].axis[0].dfy / (object[0].x1 - object[0].x0);
     187  myj = graph[0].axis[1].dfy / (object[0].y1 - object[0].y0);
     188 
     189  bxi  =  graph[0].axis[0].fx - object[0].x0*graph[0].axis[0].dfx/(object[0].x1 - object[0].x0);
     190  bxj  =  -object[0].y0*graph[0].axis[1].dfx/(object[0].y1 - object[0].y0);
     191  byi  =  -object[0].x0*graph[0].axis[0].dfy/(object[0].x1 - object[0].x0);
     192  byj  =  graph[0].axis[1].fy - object[0].y0*graph[0].axis[1].dfy/(object[0].y1 - object[0].y0);
     193 
     194  bx = bxi + bxj;
     195  by = byi + byj;
     196 
     197  X0 = graph[0].axis[0].fx;
     198  X1 = graph[0].axis[0].fx + graph[0].axis[0].dfx;
     199  Y0 = graph[0].axis[1].fy;
     200  Y1 = graph[0].axis[1].fy + graph[0].axis[1].dfy;
     201
     202  x = object[0].x; y = object[0].y;
     203
     204  // ptype must > 2
     205  // Npts % ptype must be 0
     206  // who must validate that?
     207
     208  // each polygon is made of (N = ptype) points
     209  // we connect each point and the last one back
     210  for (int i = 0; i < object[0].Npts; i += object[0].ptype) {
     211    // first check for any invalid values for this polygon
     212    int skipObject = FALSE;
     213    for (int j = 0; (j < object[0].ptype) && !skipObject; j++) {
     214      int k = i + j;
     215      if (!(finite(x[k]) && finite(y[k]))) skipObject = TRUE;
     216    }
     217    if (skipObject) continue;
     218
     219    for (int j = 0; (j < object[0].ptype); j++) {
     220      int k = i + j;
     221      sx0 = x[k]*mxi + y[k]*mxj + bx + XCENTER;
     222      sy0 = x[k]*myi + y[k]*myj + by + YCENTER;
     223
     224      // last point connects to first
     225      if (j == object[0].ptype - 1) {
     226        sx1 = x[i]*mxi + y[i]*mxj + bx + XCENTER;
     227        sy1 = x[i]*myi + y[i]*myj + by + YCENTER;
     228      } else {
     229        sx1 = x[k+1]*mxi + y[k+1]*mxj + bx + XCENTER;
     230        sy1 = x[k+1]*myi + y[k+1]*myj + by + YCENTER;
     231      }
     232      ClipLine (graphic, sx0, sy0, sx1, sy1, X0, Y0, X1, Y1);
     233    }
     234  }
     235}
     236
     237/******/
     238void DrawPolyfill (Graphic *graphic, KapaGraphWidget *graph, Gobjects *object) {
     239
     240  float *x, *y, *z;
     241  double mxi, mxj, myi, myj, bxi, bxj, byi, byj, bx, by;
     242  double sx0, sy0;
     243  // double X0, X1, Y0, Y1;
     244
     245  mxi = graph[0].axis[0].dfx / (object[0].x1 - object[0].x0);
     246  mxj = graph[0].axis[1].dfx / (object[0].y1 - object[0].y0);
     247  myi = graph[0].axis[0].dfy / (object[0].x1 - object[0].x0);
     248  myj = graph[0].axis[1].dfy / (object[0].y1 - object[0].y0);
     249 
     250  bxi  =  graph[0].axis[0].fx - object[0].x0*graph[0].axis[0].dfx/(object[0].x1 - object[0].x0);
     251  bxj  =  -object[0].y0*graph[0].axis[1].dfx/(object[0].y1 - object[0].y0);
     252  byi  =  -object[0].x0*graph[0].axis[0].dfy/(object[0].x1 - object[0].x0);
     253  byj  =  graph[0].axis[1].fy - object[0].y0*graph[0].axis[1].dfy/(object[0].y1 - object[0].y0);
     254 
     255  bx = bxi + bxj;
     256  by = byi + byj;
     257 
     258  /*
     259    window boundary so we can clip objects
     260  X0 = graph[0].axis[0].fx;
     261  X1 = graph[0].axis[0].fx + graph[0].axis[0].dfx;
     262  Y0 = graph[0].axis[1].fy;
     263  Y1 = graph[0].axis[1].fy + graph[0].axis[1].dfy;
     264  */
     265
     266  x = object[0].x; y = object[0].y; z = object[0].z;
     267
     268  // ptype must > 2
     269  // Npts % ptype must be 0
     270  // who must validate that?
     271
     272  ALLOCATE_PTR (points, XPoint, object[0].ptype);
     273
     274  // NOTE that LoadObject.c:45 limits the allow
     275  // object styles which may have negative (scaled) colors
     276  int scaleColor = (object[0].color < 0);
     277
     278  // each polygon is made of (N = ptype) points
     279  // we connect each point and the last one back
     280  for (int i = 0; i < object[0].Npts; i += object[0].ptype) {
     281    // first check for any invalid values for this polygon
     282    int skipObject = FALSE;
     283    for (int j = 0; (j < object[0].ptype) && !skipObject; j++) {
     284      int k = i + j;
     285      if (!(finite(x[k]) && finite(y[k]))) skipObject = TRUE;
     286    }
     287    if (skipObject) continue;
     288
     289    for (int j = 0; (j < object[0].ptype); j++) {
     290      int k = i + j;
     291      sx0 = x[k]*mxi + y[k]*mxj + bx + XCENTER;
     292      sy0 = x[k]*myi + y[k]*myj + by + YCENTER;
     293
     294      points[j].x = sx0;
     295      points[j].y = sy0;
     296    }
     297
     298    if (scaleColor) {
     299      if (!finite(z[i])) continue;
     300      int pixel = MIN (graphic->Npixels - 2, MAX (0, z[i]*(graphic->Npixels - 1)));
     301      XSetForeground (graphic->display, graphic->gc, graphic->cmap[pixel].pixel);
     302    }
     303    XFillPolygon (graphic->display, graphic->window, graphic->gc, points, object[0].ptype, Convex, CoordModeOrigin);
     304  }
     305
     306  free (points);
    175307}
    176308
    … …  
    580712  int scaleColor = (object[0].color < 0);
    581713
     714  // NOTE that LoadObject.c:45 limits the allow
     715  // object styles which may have negative (scaled) colors
     716
    582717  ds = 0.5 * (graphic->dx + graphic->dy) * 0.003 * object[0].size;
    583718  dz = 0.5 * (graphic->dx + graphic->dy) * 0.010;
    … …  
    9161051    // for open circles, only go to the outer radius
    9171052    D = 0;
    918     if (object[0].ptype ==  7) { D = scaleSize ? dz*z[i] : ds; }
    919     if (object[0].ptype ==  1) { D = scaleSize ? dz*z[i] : ds; }
    920     if (object[0].ptype ==  5) { D = scaleSize ? 0.66*dz*z[i] : 0.66*ds; }
    921     if (object[0].ptype == 15) { D = scaleSize ? 0.66*dz*z[i] : 0.66*ds; }
     1053    if (object[0].ptype == KAPA_POINT_CIRCLE_OPEN       ) { D = scaleSize ? dz*z[i] : ds; }
     1054    if (object[0].ptype == KAPA_POINT_BOX_OPEN          ) { D = scaleSize ? dz*z[i] : ds; }
     1055    if (object[0].ptype == KAPA_POINT_TRIANGLE_OPEN     ) { D = scaleSize ? 0.66*dz*z[i] : 0.66*ds; }
     1056    if (object[0].ptype == KAPA_POINT_TRIANGLE_OPEN_DOWN) { D = scaleSize ? 0.66*dz*z[i] : 0.66*ds; }
    9221057    if (!(finite(x[i]) && finite(y[i]) && finite(dxp[i]))) goto skip_dxp;
    9231058    if (D > fabs(dxp[i]*mxi)) goto skip_dxp;
    … …  
    9991134    // for open circles, only go to the outer radius
    10001135    D = 0;
    1001     if (object[0].ptype ==  7) { D = scaleSize ? dz*z[i] : ds; }
    1002     if (object[0].ptype ==  1) { D = scaleSize ? dz*z[i] : ds; }
    1003     if (object[0].ptype ==  5) { D = scaleSize ? 1.15*dz*z[i] : 1.15*ds; }
    1004     if (object[0].ptype == 15) { D = scaleSize ? 0.58*dz*z[i] : 0.58*ds; }
     1136    if (object[0].ptype == KAPA_POINT_CIRCLE_OPEN       ) { D = scaleSize ? dz*z[i] : ds; }
     1137    if (object[0].ptype == KAPA_POINT_BOX_OPEN          ) { D = scaleSize ? dz*z[i] : ds; }
     1138    if (object[0].ptype == KAPA_POINT_TRIANGLE_OPEN     ) { D = scaleSize ? 1.15*dz*z[i] : 1.15*ds; }
     1139    if (object[0].ptype == KAPA_POINT_TRIANGLE_OPEN_DOWN) { D = scaleSize ? 0.58*dz*z[i] : 0.58*ds; }
    10051140    if (!(finite(x[i]) && finite(y[i]) && finite(dyp[i]))) goto skip_dyp;
    10061141    if (D > fabs(dyp[i]*myj)) goto skip_dyp;
    … …  
    10231158  skip_dyp:
    10241159    if (!(finite(x[i]) && finite(y[i]) && finite(dym[i]))) continue;
    1025     if (object[0].ptype ==  5) { D = scaleSize ? 0.58*dz*z[i] : 0.58*ds; }
    1026     if (object[0].ptype == 15) { D = scaleSize ? 1.15*dz*z[i] : 1.15*ds; }
     1160    if (object[0].ptype == KAPA_POINT_TRIANGLE_OPEN     ) { D = scaleSize ? 0.58*dz*z[i] : 0.58*ds; }
     1161    if (object[0].ptype == KAPA_POINT_TRIANGLE_OPEN_DOWN) { D = scaleSize ? 1.15*dz*z[i] : 1.15*ds; }
    10271162    if (D > fabs(dym[i]*myj)) continue;
    10281163    sx0 = x[i]*mxi + y[i]*mxj + bx + XCENTER;

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/Image.c

    r39457 r41170  
    148148void DrawImage (KapaImageWidget *image) {
    149149
     150  Graphic *graphic;
     151
     152  if (image == NULL) return;
     153
     154  graphic = GetGraphic ();
     155
     156  if (image[0].picture.pix) {
     157    XPutImage (graphic[0].display, graphic[0].window, graphic[0].gc,
     158               image[0].picture.pix, 0, 0,
     159               image[0].picture.x + 1, image[0].picture.y + 1,
     160               image[0].picture.dx, image[0].picture.dy);
     161  }
     162}
     163
     164// add the zoom, pan, crosshairs, status box, buttons
     165void DrawImageTool (KapaImageWidget *image) {
     166
    150167  int i;
    151168  Graphic *graphic;
    … …  
    160177                  image[0].picture.dx+1, image[0].picture.dy+1);
    161178 
    162   if (image[0].picture.pix) {
    163     XPutImage (graphic[0].display, graphic[0].window, graphic[0].gc,
    164                image[0].picture.pix, 0, 0,
    165                image[0].picture.x + 1, image[0].picture.y + 1,
    166                image[0].picture.dx, image[0].picture.dy);
    167   }
    168 
    169   for (i = 0; i < NOVERLAYS; i++) {
    170     if (image[0].overlay[i].active) {
    171       PaintOverlay (graphic, image, i);
    172     }
    173   }
    174 
    175179  if (image[0].location) {
    176180    if (image[0].cmapbar.pix) {
    … …  
    211215    StatusBox (graphic, image);
    212216  }
    213 
    214   FlushDisplay ();
    215217}
    216218

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/InterpretKeys.c

    r38465 r41170  
    160160      break;
    161161
     162    case XK_plus:
     163    case XK_equal:
    162164    case XK_KP_Add:
    163165      if (modstate & ControlMask) {
    … …  
    171173      Reorient (graphic, image, image[0].picture.Xc, image[0].picture.Yc, 0);
    172174      break;
     175
     176    case XK_minus:
     177    case XK_underscore:
    173178    case XK_KP_Subtract:
    174179      if (modstate & ControlMask) {

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/JPEGit24.c

    r34088 r41170  
    212212    }
    213213    bDrawBufferFree (buffer);
     214    free (palette);
    214215  }
    215216

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/Layout.c

    r35237 r41170  
    4545  AddSection ("default", 0.0, 0.0, 1.0, 1.0, -1);
    4646}
     47
     48void FreeLayout (void) {
     49  FreeSections ();
     50  FreeRotFonts ();
     51}

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/LoadFrame.c

    r38986 r41170  
    9393  SetSectionSizes (section);
    9494  Refresh();
    95   // if (USE_XWINDOW) DrawFrame (graph);
    96   // FlushDisplay ();
    97 
    98   /* XXX why did I do this??? */
    99 # if (0)
    100   status = TRUE;
    101   if (status) {
    102     for (i = 0; i < Nsection; i++) {
    103       PositionPicture (&section[i]);
    104     }
    105     if (USE_XWINDOW) XClearWindow (graphic.display, graphic.window);
    106     Refresh ();
    107   } else {
    108     FlushDisplay ();
    109   }
    110 # endif
    11195
    11296  return (TRUE);

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/LoadLabels.c

    r39577 r41170  
    2626  label[Nbytes] = 0;
    2727 
     28  FREE (label);
     29
    2830  c = GetRotFont (&size);
    2931

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/LoadObject.c

    r40570 r41170  
    4242      graph[0].objects[N].color = KapaColormapSize() - 1;
    4343  }
    44   if ((graph[0].objects[N].style != KAPA_PLOT_POINTS) && (graph[0].objects[N].color < 0)) {
     44
     45  // XXX watch out for this restriction in DrawObjects / bDrawObjects
     46  int canScaleColor = FALSE;
     47  canScaleColor |= (graph[0].objects[N].style == KAPA_PLOT_POINTS);
     48  canScaleColor |= (graph[0].objects[N].style == KAPA_PLOT_POLYFILL);
     49
     50  if (!canScaleColor && (graph[0].objects[N].color < 0)) {
    4551      graph[0].objects[N].color = 0;
    4652  }
    … …  
    103109
    104110  if (USE_XWINDOW) {
    105     Graphic *graphic = GetGraphic();
    106     DrawObjectN (graphic, graph, &graph[0].objects[graph[0].Nobjects-1]);
    107   }
    108   FlushDisplay ();
     111    if (0) {
     112      // use this if we are not using buffered plotting
     113      Graphic *graphic = GetGraphic();
     114      DrawObjectN (graphic, graph, &graph[0].objects[graph[0].Nobjects-1]);
     115      FlushDisplay ();
     116    } else {
     117      Refresh ();
     118    }
     119  }
    109120
    110121  return (TRUE);

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/PNGit.c

    r39457 r41170  
    9898  png_write_info (png_ptr, info_ptr);
    9999
     100  if (graphic->smooth_sigma == 0.0) {
     101    fprintf (stderr, "making PNG without antialias smoothing\n");
     102    fprintf (stderr, "use 'antialias (value)' to set smoothing scale (0.4 - 0.7 recommended)\n");
     103  }
     104   
    100105  if (haveImage) {
    101106    buffer = bDrawIt (palette, Npalette, 3);
    … …  
    109114 
    110115  bDrawBufferFree (buffer);
     116  free (palette);
     117 
    111118  fclose (f);
    112119  return (TRUE);

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/PPMit.c

    r29938 r41170  
    4646 
    4747  bDrawBufferFree (buffer);
     48  free (palette);
     49
    4850  return (TRUE);
    4951}

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/Refresh.c

    r36833 r41170  
    11# include "Ximage.h"
    22
     3# define USE_BUFFERED_DRAW 1
     4
     5void bDrawXimage (bDrawBuffer *buffer);
     6void Refresh_Buffered (void);
     7void Refresh_Unbuffered (void);
     8
    39void Refresh (void) {
     10  if (USE_BUFFERED_DRAW) {
     11    Refresh_Buffered();
     12  } else {
     13    Refresh_Unbuffered();
     14  }
     15}
     16
     17void Refresh_Buffered (void) {
     18
     19  int Npalette;
     20  Graphic *graphic;
     21
     22  if (!USE_XWINDOW) return;
     23  // if (HAVE_BACKING) return;
     24
     25  graphic = GetGraphic();
     26 
     27  // limit the png window to the min needed to contain the active graphic regions
     28  SectionMinBoundary (graphic);
     29
     30  // is the palette reasonable in modern context?
     31  png_color *palette = KapaPNGPalette (&Npalette);
     32
     33  bDrawBuffer *buffer = bDrawIt (palette, Npalette, 3);
     34 
     35  /* XClearWindow   (graphic.display, graphic.window); */
     36  XSetForeground (graphic->display, graphic->gc, graphic->back);
     37  XFillRectangle (graphic->display, graphic->window, graphic->gc, 0, 0, graphic->dx, graphic->dy);
     38  XSetForeground (graphic->display, graphic->gc, graphic->fore);
     39 
     40  // copy buffer to Xwindow as image?
     41  bDrawXimage (buffer);
     42  bDrawBufferFree (buffer);
     43  free (palette);
     44
     45  // draw image tool for all sections
     46  int Nsection = GetNumberOfSections ();
     47  for (int i = 0; i < Nsection; i++) {
     48    Section *section = GetSectionByNumber (i);
     49
     50    KapaImageWidget *image = section->image;
     51    DrawImageTool (image);
     52
     53    /*** PaintOverlay is called in DrawImage ***
     54    if (!image) continue;
     55    for (int j = 0; j < NOVERLAYS; j++) {
     56      if (image[0].overlay[j].active) {
     57        PaintOverlay (graphic, image, j);
     58      }
     59    }
     60    */
     61  }
     62
     63  FlushDisplay ();
     64}
     65
     66void Refresh_Unbuffered (void) {
    467
    568  int i, Nsection;
    … …  
    2083  Nsection = GetNumberOfSections ();
    2184  for (i = 0; i < Nsection; i++) {
    22       section = GetSectionByNumber (i);
    23       DrawSectionBG (graphic, section);
    24       DrawImage (section->image);
    25       DrawGraph (section->graph);
     85    section = GetSectionByNumber (i);
     86    DrawSectionBG (graphic, section);
     87
     88    KapaImageWidget *image = section->image;
     89    DrawImage (image);
     90    DrawImageTool (image);
     91
     92    /*** the overlay is added in DrawImage  ***/
     93    if (image) {
     94      for (int j = 0; j < NOVERLAYS; j++) {
     95        if (image[0].overlay[j].active) {
     96          PaintOverlay (graphic, image, j);
     97        }
     98      }
     99    }
     100
     101    DrawGraph (section->graph);
    26102  }
    27103

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/SetColormap.c

    r40627 r41170  
    7272          }
    7373      }
    74       fprintf (stderr, "ruff Npix: %d vs %d\n", i, graphic[0].Npixels);
     74      // fprintf (stderr, "ruff Npix: %d vs %d\n", i, graphic[0].Npixels);
    7575
    7676      // all other modes are 1D: set the flag:
    … …  
    217217      !strncmp (graphic->colormapName, "lgcy:", 5) ||
    218218      !strncmp (graphic->colormapName, "cetf:", 5) ||
     219      !strncmp (graphic->colormapName, "cetr:", 5) ||
    219220      !strncmp (graphic->colormapName, "csvf:", 5)) {
    220221
    … …  
    233234    int isCSV = !strncmp (graphic->colormapName, "csvf:", 5);
    234235    int isCET = !strncmp (graphic->colormapName, "cetf:", 5);
     236    int isCETRev = !strncmp (graphic->colormapName, "cetr:", 5);
    235237    int isLegacy = !strncmp (graphic->colormapName, "lgcy:", 5);
    236238
    237239    float fracIndex, fracRed, fracBlue, fracGreen;
    238240
    239     if (isCET) fracIndex = 0.0;
     241    if (isCET || isCETRev) fracIndex = 0.0;
    240242
    241243    while (scan_line_maxlen (f, line, 1024) != EOF) {
    … …  
    247249        if (Nscan != 4) continue;
    248250      }
    249       if (isCET) {
     251      if (isCET || isCETRev) {
    250252        Nscan = sscanf (line, "%f,%f,%f", &fracRed, &fracGreen, &fracBlue);
    251253        fracIndex += 1.0 / 256.0;
    … …  
    256258        if (Nscan != 4) continue;
    257259      }
    258       if (!isCSV && !isLegacy && !isCET) {
     260      if (!isCSV && !isLegacy && !isCET && !isCETRev) {
    259261        Nscan = sscanf (line, "%f %f %f %f", &fracIndex, &fracRed, &fracGreen, &fracBlue);
    260262        if (Nscan != 4) continue;
    … …  
    299301      }
    300302    }
     303
     304    // reverse the color sequence:
     305    if (isCETRev) {
     306      for (i = 0; i < MaxValue / 2; i++) { 
     307        unsigned short tmp;
     308        // fprintf (stderr, "swap: %d %d\n", graphic[0].cmap[i].red, graphic[0].cmap[MaxValue - 1 - i].red);
     309        tmp = graphic[0].cmap[i].red;   graphic[0].cmap[i].red   = graphic[0].cmap[MaxValue - 1 - i].red;   graphic[0].cmap[MaxValue - 1 - i].red   = tmp;
     310        tmp = graphic[0].cmap[i].blue;  graphic[0].cmap[i].blue  = graphic[0].cmap[MaxValue - 1 - i].blue;  graphic[0].cmap[MaxValue - 1 - i].blue  = tmp;
     311        tmp = graphic[0].cmap[i].green; graphic[0].cmap[i].green = graphic[0].cmap[MaxValue - 1 - i].green; graphic[0].cmap[MaxValue - 1 - i].green = tmp;
     312      }
     313    }
     314
    301315    goto store_colors;
    302316  }

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/SetUpGraphic.c

    r39532 r41170  
    3232  graphic->pixels       = NULL;
    3333
     34  graphic->smooth_sigma = 0.0;
     35
    3436  if (!USE_XWINDOW) {
    3537    ALLOCATE (graphic[0].pixels, unsigned long, NPIXELS_STATIC);
    … …  
    5052  CheckVisual (graphic, argc, argv);
    5153  CheckColors (graphic, argc, argv);
    52 
    5354
    5455  icon.width = icon_width;
    … …  
    8485}
    8586
     87int SetSmoothSigma (int sock) {
     88
     89  float sigma;
     90  KiiScanMessage (sock, "%f", &sigma);
     91
     92  if (isfinite(sigma)) {
     93    if ((sigma <= 1.1) && (sigma >= 0.0)) {
     94      graphic->smooth_sigma = sigma;
     95    }
     96  }
     97  return TRUE;
     98}
     99
     100void FreeGraphic () {
     101  free (graphic->pixels);
     102  free (graphic->cmap);
     103  free (graphic->color);
     104  free (graphic->colormapName);
     105  free (graphic);
     106}

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/args.c

    r39225 r41170  
    2323    NAME_WINDOW = strcreate (argv[N]);
    2424    remove_argument(N, argc, argv);
     25  }
     26
     27  if ((N = get_argument (*argc, argv, "-memdump"))) {
     28    remove_argument(N, argc, argv);
     29    MemoryDumpSetOnExit (TRUE);
    2530  }
    2631

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/bDrawImage.c

    r34088 r41170  
    11# include "Ximage.h"
     2
     3// XXX for the moment, this function does NOT set the mask bits.
     4// since we lay graphics on top of images, this is OK for now
    25
    36# define WHITE_R 255
    … …  
    5457
    5558  // the created buffer is supposed to contain the output windows
     59  if (Xs < 0) {
     60    fprintf (stderr, "image display boundaries out of range of window (invalid condition) : fix Kapa\n");
     61    abort();
     62  }
    5663  if (buffer[0].Nx < Xs + dx) {
    57     fprintf (stderr, "invalid condition\n");
     64    fprintf (stderr, "image display boundaries out of range of window (invalid condition) : fix Kapa\n");
     65    abort();
     66  }
     67  if (Ys < 0) {
     68    fprintf (stderr, "image display boundaries out of range of window (invalid condition) : fix Kapa\n");
    5869    abort();
    5970  }
    6071  if (buffer[0].Ny < Ys + dy) {
    61     fprintf (stderr, "invalid condition\n");
     72    fprintf (stderr, "image display boundaries out of range of window (invalid condition) : fix Kapa\n");
    6273    abort();
    6374  }
    … …  
    153164
    154165  return (TRUE);
     166}
     167
     168void bDrawXimage (bDrawBuffer *buffer) {
     169
     170  Graphic *graphic = GetGraphic ();
     171
     172  ALLOCATE_PTR (data, char, 4*buffer->Nx*buffer->Ny);
     173
     174  for (int iy = 0; iy < buffer->Ny; iy++) {
     175    for (int ix = 0; ix < buffer->Nx; ix++) {
     176      data[4*(iy*buffer->Nx + ix) + 0] = buffer->pixels[iy][3*ix + 2];
     177      data[4*(iy*buffer->Nx + ix) + 1] = buffer->pixels[iy][3*ix + 1];
     178      data[4*(iy*buffer->Nx + ix) + 2] = buffer->pixels[iy][3*ix + 0];
     179      data[4*(iy*buffer->Nx + ix) + 3] = 0;
     180    }
     181  }
     182
     183  XImage *pix = XCreateImage (graphic[0].display, graphic[0].visual, graphic[0].depth, ZPixmap, 0,
     184                              data, buffer->Nx, buffer->Ny, 32, 0);
     185
     186  XPutImage (graphic[0].display, graphic[0].window, graphic[0].gc, pix, 0, 0, 1, 1, buffer->Nx, buffer->Ny);
     187
     188  free (data);
    155189}
    156190
    … …  
    186220    }
    187221    bDrawBufferFree (buffer);
     222    free (palette);
    188223  }
    189224

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/bDrawIt.c

    r40558 r41170  
    33bDrawBuffer *bDrawIt (png_color *palette, int Npalette, int Nbyte) {
    44
    5   int i, j, Nsection;
    6   bDrawBuffer *buffer;
    7   bDrawColor black;
    8   Graphic *graphic;
    9   Section *section;
     5  Graphic *graphic = GetGraphic();
     6  bDrawColor black = KapaColorByName ("black");
    107
    11   graphic = GetGraphic();
     8  // get the number of sections
     9  int Nsection = GetNumberOfSections ();
    1210
    13   black = KapaColorByName ("black");
     11  // in order to allow the anti-aliasing to affect the text & graphs but not the images
     12  // I need to generate the images in one buffer and the graphs in a second buffer
     13  // then merge the two buffers.
    1414
    15   // if we want to trim, we'll need to carry about the start in graphic coords and
    16   // the dx,dy size. 
    17   buffer = bDrawBufferCreate (graphic->dxwin, graphic->dywin, Nbyte, palette, Npalette);
    18   bDrawSetStyle (buffer, black, 0, 0, 1.0);
     15  // base will hold the images
     16  bDrawBuffer *base = bDrawBufferCreate (graphic->dxwin, graphic->dywin, Nbyte, palette, Npalette);
     17  bDrawSetStyle (base, black, 0, 0, 1.0);
    1918 
    20   // reset the sizes for all sections
    21   Nsection = GetNumberOfSections ();
    22   for (i = 0; i < Nsection; i++) {
    23       section = GetSectionByNumber (i);
    24       bDrawImage (buffer, section->image, graphic);
    25       for (j = 0; section->image && (j < NOVERLAYS); j++) {
    26         if (section->image->overlay[j].active) bDrawOverlay (buffer, section->image, j);
    27       }
    28       bDrawGraph (buffer, section->graph);
     19  for (int i = 0; i < Nsection; i++) {
     20    Section *section = GetSectionByNumber (i);
     21    bDrawImage (base, section->image, graphic);
    2922  }
    3023
    31   return (buffer);
     24  // graph will hold the graphic overlay
     25  bDrawBuffer *graph = bDrawBufferCreate (graphic->dxwin, graphic->dywin, Nbyte, palette, Npalette);
     26  bDrawSetStyle (graph, black, 0, 0, 1.0);
     27
     28  for (int i = 0; i < Nsection; i++) {
     29    Section *section = GetSectionByNumber (i);
     30    for (int j = 0; section->image && (j < NOVERLAYS); j++) {
     31      if (section->image->overlay[j].active) bDrawOverlay (graph, section->image, j);
     32    }
     33    bDrawGraph (graph, section->graph);
     34  }
     35
     36  // apply anti-aliasing only to the graph
     37  if (graphic->smooth_sigma > 0.0) {
     38    // anything > 1.1 blurs the image too much
     39    graphic->smooth_sigma = MIN (graphic->smooth_sigma, 1.1);
     40    bDrawSmooth (graph, graphic->smooth_sigma);
     41  }
     42 
     43  // place graph on base
     44  bDrawMerge (base, graph);
     45  bDrawBufferFree (graph);
     46
     47  return (base);
    3248}
    3349

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/bDrawObjects.c

    r40570 r41170  
    99# define OpenTriangle(BUF,X1,Y1,X2,Y2,X3,Y3) (bDrawTriOpen (BUF, (X1), (Y1), (X2), (Y2), (X3), (Y3)))
    1010
     11// I should not have a local static variable for this:
     12// I should just pass the graphic structure to the called functions below
    1113static Graphic *graphic;
    1214
    … …  
    3537    case KAPA_PLOT_CONNECT:
    3638      bDrawConnect (buffer, graph, object);
     39      break;
     40    case KAPA_PLOT_POLYGON:
     41      bDrawPolygons (buffer, graph, object);
     42      break;
     43    case KAPA_PLOT_POLYFILL:
     44      bDrawFillPolygons (buffer, graph, object);
    3745      break;
    3846    case KAPA_PLOT_HISTOGRAM:
    … …  
    103111    sx0 = sx1; sy0 = sy1;
    104112  }
     113}
     114
     115void bDrawPolygons (bDrawBuffer *buffer, KapaGraphWidget *graph, Gobjects *object) {
     116
     117  float *x, *y;
     118  double mxi, mxj, myi, myj, bxi, bxj, byi, byj, bx, by;
     119  double sx0, sy0, sx1, sy1;
     120  double X0, X1, Y0, Y1;
     121
     122  mxi = graph[0].axis[0].dfx / (object[0].x1 - object[0].x0);
     123  mxj = graph[0].axis[1].dfx / (object[0].y1 - object[0].y0);
     124  myi = graph[0].axis[0].dfy / (object[0].x1 - object[0].x0);
     125  myj = graph[0].axis[1].dfy / (object[0].y1 - object[0].y0);
     126 
     127  bxi  =  graph[0].axis[0].fx - object[0].x0*graph[0].axis[0].dfx/(object[0].x1 - object[0].x0);
     128  bxj  =  -object[0].y0*graph[0].axis[1].dfx/(object[0].y1 - object[0].y0);
     129  byi  =  -object[0].x0*graph[0].axis[0].dfy/(object[0].x1 - object[0].x0);
     130  byj  =  graph[0].axis[1].fy - object[0].y0*graph[0].axis[1].dfy/(object[0].y1 - object[0].y0);
     131 
     132  bx = bxi + bxj;
     133  by = byi + byj;
     134 
     135  X0 = graph[0].axis[0].fx;
     136  X1 = graph[0].axis[0].fx + graph[0].axis[0].dfx;
     137  Y0 = graph[0].axis[1].fy;
     138  Y1 = graph[0].axis[1].fy + graph[0].axis[1].dfy;
     139
     140  x = object[0].x; y = object[0].y;
     141
     142  // ptype must > 2
     143  // Npts % ptype must be 0
     144  // who must validate that?
     145
     146  // each polygon is made of (N = ptype) points
     147  // we connect each point and the last one back
     148  for (int i = 0; i < object[0].Npts; i += object[0].ptype) {
     149    // first check for any invalid values for this polygon
     150    int skipObject = FALSE;
     151    for (int j = 0; (j < object[0].ptype) && !skipObject; j++) {
     152      int k = i + j;
     153      if (!(finite(x[k]) && finite(y[k]))) skipObject = TRUE;
     154    }
     155    if (skipObject) continue;
     156
     157    for (int j = 0; (j < object[0].ptype); j++) {
     158      int k = i + j;
     159      sx0 = x[k]*mxi + y[k]*mxj + bx;
     160      sy0 = x[k]*myi + y[k]*myj + by;
     161
     162      // last point connects to first
     163      if (j == object[0].ptype - 1) {
     164        sx1 = x[i]*mxi + y[i]*mxj + bx;
     165        sy1 = x[i]*myi + y[i]*myj + by;
     166      } else {
     167        sx1 = x[k+1]*mxi + y[k+1]*mxj + bx;
     168        sy1 = x[k+1]*myi + y[k+1]*myj + by;
     169      }
     170      bDrawClipLine (buffer, sx0, sy0, sx1, sy1, X0, Y0, X1, Y1);
     171    }
     172  }
     173}
     174
     175void bDrawFillPolygons (bDrawBuffer *buffer, KapaGraphWidget *graph, Gobjects *object) {
     176
     177  float *x, *y;
     178  double mxi, mxj, myi, myj, bxi, bxj, byi, byj, bx, by;
     179  double sx0, sy0;
     180  // double X0, X1, Y0, Y1;
     181
     182  mxi = graph[0].axis[0].dfx / (object[0].x1 - object[0].x0);
     183  mxj = graph[0].axis[1].dfx / (object[0].y1 - object[0].y0);
     184  myi = graph[0].axis[0].dfy / (object[0].x1 - object[0].x0);
     185  myj = graph[0].axis[1].dfy / (object[0].y1 - object[0].y0);
     186 
     187  bxi  =  graph[0].axis[0].fx - object[0].x0*graph[0].axis[0].dfx/(object[0].x1 - object[0].x0);
     188  bxj  =  -object[0].y0*graph[0].axis[1].dfx/(object[0].y1 - object[0].y0);
     189  byi  =  -object[0].x0*graph[0].axis[0].dfy/(object[0].x1 - object[0].x0);
     190  byj  =  graph[0].axis[1].fy - object[0].y0*graph[0].axis[1].dfy/(object[0].y1 - object[0].y0);
     191 
     192  bx = bxi + bxj;
     193  by = byi + byj;
     194 
     195  /*
     196    window boundary so we can clip objects
     197  X0 = graph[0].axis[0].fx;
     198  X1 = graph[0].axis[0].fx + graph[0].axis[0].dfx;
     199  Y0 = graph[0].axis[1].fy;
     200  Y1 = graph[0].axis[1].fy + graph[0].axis[1].dfy;
     201  */
     202
     203  x = object[0].x; y = object[0].y;
     204
     205  // ptype must > 2
     206  // Npts % ptype must be 0
     207  // who must validate that?
     208
     209  ALLOCATE_PTR (xpts, double, object[0].ptype);
     210  ALLOCATE_PTR (ypts, double, object[0].ptype);
     211
     212  int scaleColor = (object[0].color < 0);
     213  unsigned char *pixel1, *pixel2, *pixel3;
     214  float *z = object[0].z;
     215
     216  if (scaleColor) {
     217    // scaled colors use the colormap defined for the graphic
     218    ALLOCATE (pixel1, unsigned char, graphic[0].Npixels);
     219    ALLOCATE (pixel2, unsigned char, graphic[0].Npixels);
     220    ALLOCATE (pixel3, unsigned char, graphic[0].Npixels);
     221
     222    /** cmap[i].pixel must be defined even if X is not used **/
     223    for (int i = 0; i < graphic[0].Npixels; i++) { /* set up pixel array */
     224      pixel1[i] = graphic[0].cmap[i].red >> 8;
     225      pixel2[i] = graphic[0].cmap[i].green >> 8;
     226      pixel3[i] = graphic[0].cmap[i].blue >> 8;
     227    }
     228
     229  }
     230
     231  // each polygon is made of (N = ptype) points
     232  // we connect each point and the last one back
     233  for (int i = 0; i < object[0].Npts; i += object[0].ptype) {
     234    // first check for any invalid values for this polygon
     235    int skipObject = FALSE;
     236    for (int j = 0; (j < object[0].ptype) && !skipObject; j++) {
     237      int k = i + j;
     238      if (!(finite(x[k]) && finite(y[k]))) skipObject = TRUE;
     239    }
     240    if (skipObject) continue;
     241
     242    for (int j = 0; (j < object[0].ptype); j++) {
     243      int k = i + j;
     244      sx0 = x[k]*mxi + y[k]*mxj + bx;
     245      sy0 = x[k]*myi + y[k]*myj + by;
     246
     247      xpts[j] = sx0;
     248      ypts[j] = sy0;
     249    }
     250   
     251    if (scaleColor) {
     252      if (!finite(z[i])) continue;
     253      int pixel = MIN (graphic->Npixels - 2, MAX (0, z[i]*(graphic->Npixels - 1)));
     254      buffer->bColor_R = pixel1[pixel];
     255      buffer->bColor_G = pixel2[pixel];
     256      buffer->bColor_B = pixel3[pixel];
     257    }
     258    bDrawPolyFill (buffer, xpts, ypts, object[0].ptype);
     259  }
     260  free (xpts);
     261  free (ypts);
    105262}
    106263
    … …  
    8861043    // for open circles, only go to the outer radius
    8871044    D = 0;
    888     if (object[0].ptype ==  7) { D = scaleSize ? dz*z[i] : ds; }
    889     if (object[0].ptype ==  1) { D = scaleSize ? dz*z[i] : ds; }
    890     if (object[0].ptype ==  5) { D = scaleSize ? 0.66*dz*z[i] : 0.66*ds; }
    891     if (object[0].ptype == 15) { D = scaleSize ? 0.66*dz*z[i] : 0.66*ds; }
     1045    if (object[0].ptype == KAPA_POINT_CIRCLE_OPEN       ) { D = scaleSize ? dz*z[i] : ds; }
     1046    if (object[0].ptype == KAPA_POINT_BOX_OPEN          ) { D = scaleSize ? dz*z[i] : ds; }
     1047    if (object[0].ptype == KAPA_POINT_TRIANGLE_OPEN     ) { D = scaleSize ? 0.66*dz*z[i] : 0.66*ds; }
     1048    if (object[0].ptype == KAPA_POINT_TRIANGLE_OPEN_DOWN) { D = scaleSize ? 0.66*dz*z[i] : 0.66*ds; }
    8921049    if (!(finite(x[i]) && finite(y[i]) && finite(dxp[i]))) goto skip_dxp;
    8931050    if (D > fabs(dxp[i]*mxi)) goto skip_dxp;
    … …  
    9691126    // for open circles, only go to the outer radius
    9701127    D = 0;
    971     if (object[0].ptype ==  7) { D = scaleSize ? dz*z[i] : ds; }
    972     if (object[0].ptype ==  1) { D = scaleSize ? dz*z[i] : ds; }
    973     if (object[0].ptype ==  5) { D = scaleSize ? 1.15*dz*z[i] : 1.15*ds; }
    974     if (object[0].ptype == 15) { D = scaleSize ? 0.58*dz*z[i] : 0.58*ds; }
     1128    if (object[0].ptype == KAPA_POINT_CIRCLE_OPEN       ) { D = scaleSize ? dz*z[i] : ds; }
     1129    if (object[0].ptype == KAPA_POINT_BOX_OPEN          ) { D = scaleSize ? dz*z[i] : ds; }
     1130    if (object[0].ptype == KAPA_POINT_TRIANGLE_OPEN     ) { D = scaleSize ? 1.15*dz*z[i] : 1.15*ds; }
     1131    if (object[0].ptype == KAPA_POINT_TRIANGLE_OPEN_DOWN) { D = scaleSize ? 0.58*dz*z[i] : 0.58*ds; }
    9751132    if (!(finite(x[i]) && finite(y[i]) && finite(dyp[i]))) goto skip_dyp;
    9761133    if (D > fabs(dyp[i]*myj)) goto skip_dyp;
    … …  
    9931150  skip_dyp:
    9941151    if (!(finite(x[i]) && finite(y[i]) && finite(dym[i]))) continue;
    995     if (object[0].ptype ==  5) { D = scaleSize ? 0.58*dz*z[i] : 0.58*ds; }
    996     if (object[0].ptype == 15) { D = scaleSize ? 1.15*dz*z[i] : 1.15*ds; }
     1152    if (object[0].ptype == KAPA_POINT_TRIANGLE_OPEN     ) { D = scaleSize ? 0.58*dz*z[i] : 0.58*ds; }
     1153    if (object[0].ptype == KAPA_POINT_TRIANGLE_OPEN_DOWN) { D = scaleSize ? 1.15*dz*z[i] : 1.15*ds; }
    9971154    if (D > fabs(dym[i]*myj)) continue;
    9981155    sx0 = x[i]*mxi + y[i]*mxj + bx;

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/bDrawOverlay.c

    r40558 r41170  
    2424  }
    2525
    26   Xmin = 0;
    27   Ymin = 0;
    28   Xmax = image[0].picture.dx;
    29   Ymax = image[0].picture.dy;
     26  // Xmin = 0;
     27  // Ymin = 0;
     28  // Xmax = image[0].picture.dx;
     29  // Ymax = image[0].picture.dy;
     30
     31  Xmin = image[0].picture.x;
     32  Ymin = image[0].picture.y;
     33  Xmax = image[0].picture.x + image[0].picture.dx; // maybe this should be just dx?
     34  Ymax = image[0].picture.y + image[0].picture.dy;
    3035
    3136  if (N == INFRONT) {
    … …  
    3540
    3641  for (i = 0; i < image[0].overlay[N].Nobjects; i++) {
    37     Image_to_Picture (&X, &Y, image[0].overlay[N].objects[i].x, image[0].overlay[N].objects[i].y, &image[0].picture);
     42    // XXX the 0.5,0.5 offset is apparently needed here.
     43    // work on rationalizing these functions in the context of their different plotting types
     44    Image_to_Screen (&X, &Y, image[0].overlay[N].objects[i].x - 0.5, image[0].overlay[N].objects[i].y - 0.5, &image[0].picture);
    3845    dX = image[0].overlay[N].objects[i].dx * expand;
    3946    dY = image[0].overlay[N].objects[i].dy * expand;

  • branches/eam_branches/ipp-20191011/Ohana/src/kapa2/src/kapa.c

    r13479 r41170  
    1111
    1212  CloseDisplay ();
    13   exit (0);
     13
     14  // free things
     15  FreeLayout();
     16  FreeGraphic();
     17  FREE (NAME_WINDOW);
     18
     19  MemoryDumpAndExit ();
    1420}

  • branches/eam_branches/ipp-20191011/Ohana/src/libdvo/src/dvosorts.c

    r40397 r41170  
    7676}
    7777
    78 /* sort a coordinate pair (X,Y) and the associated index (S) */
     78/* sort the index of a coordinate pair (X,Y) on X (vector pair stays unsorted) */
     79// XXX isn't this function equivalent to dsort_indexonly?
    7980void sort_coords_indexonly (double *X, double *Y, off_t *S, off_t N) {
    8081  OHANA_UNUSED_PARAM(Y);

  • branches/eam_branches/ipp-20191011/Ohana/src/libkapa/include/kapa.h

    r40558 r41170  
    4545  KAPA_PLOT_BARS_OUTLINE =  4,
    4646  KAPA_PLOT_BARS_OUTFILL =  5,
    47   KAPA_PLOT_INVALID_MAX  =  6,
     47  KAPA_PLOT_POLYGON      =  6,
     48  KAPA_PLOT_POLYFILL     =  7,
     49  KAPA_PLOT_INVALID_MAX  =  8,
    4850} KapaPlotStyle;
    4951
    … …  
    172174  int Nx, Ny, Nbyte;
    173175  bDrawColor **pixels;
     176  char **mask;
    174177  png_color *palette;
    175178  int Npalette;
     179
    176180  // current drawing values:
    177181  int bWeight;
    … …  
    256260int KapaGetImageData (int fd, KapaImageData *graphmode);
    257261int KapaSetToolbox (int fd, int location);
     262int KapaSetSmoothSigma (int fd, float sigma);
     263int KapaMemoryDump (int fd);
     264int KapaMemoryDumpLines (int fd, int Nlines);
     265int KapaMemoryDumpOnExit (int fd, int state);
    258266
    259267/* KapaColors */
    … …  
    272280/* RotFont.c */
    273281void InitRotFonts PROTO((void));
     282void FreeRotFonts PROTO((void));
    274283int SetRotFont PROTO((char *name, int size));
    275284char *GetRotFont PROTO((int *size));
    … …  
    292301bDrawBuffer *bDrawBufferCreate (int Nx, int Ny, int Nbyte, png_color *palette, int Npalette);
    293302void bDrawBufferFree (bDrawBuffer *buffer);
     303int bDrawMerge (bDrawBuffer *base, bDrawBuffer *layer);
    294304void bDrawSetBuffer (bDrawBuffer *buffer);
    295305void bDrawSetColor (bDrawBuffer *buffer, bDrawColor color);
    … …  
    312322void bDrawTriOpen  (bDrawBuffer *buffer, double x1, double y1, double x2, double y2, double x3, double y3);
    313323void bDrawTriFill  (bDrawBuffer *buffer, double x1, double y1, double dx, double dy);
     324void bDrawPolyFill (bDrawBuffer *buffer, double *x, double *y, int Npoints);
     325
     326void bDrawSmooth (bDrawBuffer *buffer, float sigma);
    314327
    315328/* bDrawRotFont.c */

  • branches/eam_branches/ipp-20191011/Ohana/src/libkapa/src/KapaStyles.c

    r40106 r41170  
    5454  if (!strcasecmp (string,  "outline")) return KAPA_PLOT_BARS_OUTLINE;
    5555  if (!strcasecmp (string,  "outfill")) return KAPA_PLOT_BARS_OUTFILL;
     56
     57  if (!strcasecmp (string,  "polygon"))   return KAPA_PLOT_POLYGON;
     58  if (!strcasecmp (string,  "polyfill"))  return KAPA_PLOT_POLYFILL;
    5659
    5760  if (strlen(string) > 2) {

  • branches/eam_branches/ipp-20191011/Ohana/src/libkapa/src/KapaWindow.c

    r40570 r41170  
    455455}
    456456
     457int KapaSetSmoothSigma (int fd, float sigma) {
     458
     459  KiiSendCommand (fd, 4, "SIGM");
     460  KiiSendMessage (fd, "%4.1f ", sigma);
     461  KiiWaitAnswer (fd, "DONE");
     462  return (TRUE);
     463}
     464
     465int KapaMemoryDump (int fd) {
     466
     467  KiiSendCommand (fd, 4, "MEMD");
     468  KiiWaitAnswer (fd, "DONE");
     469  return (TRUE);
     470}
     471
     472int KapaMemoryDumpLines (int fd, int Nlines) {
     473
     474  KiiSendCommand (fd, 4, "MEML");
     475  KiiSendMessage (fd, "%d ", Nlines);
     476  KiiWaitAnswer (fd, "DONE");
     477  return (TRUE);
     478}
     479
     480int KapaMemoryDumpOnExit (int fd, int state) {
     481
     482  KiiSendCommand (fd, 4, "MEMX");
     483  KiiSendMessage (fd, "%d ", state);
     484  KiiWaitAnswer (fd, "DONE");
     485  return (TRUE);
     486}
     487

  • branches/eam_branches/ipp-20191011/Ohana/src/libkapa/src/RotFont.c

    r39576 r41170  
    3232  strncpy (currentname, RotFonts[DEFFONT].name, 63); currentname[63] = 0;
    3333  currentsize = RotFonts[DEFFONT].size;
     34}
     35
     36void FreeRotFonts (void) {
     37  free (RotFonts);
     38  RotFonts = FALSE;
     39  RotFontInited = FALSE;
    3440}
    3541

  • branches/eam_branches/ipp-20191011/Ohana/src/libkapa/src/bDrawFuncs.c

    r40572 r41170  
    11# include <kapa_internal.h>
    22
    3 // move these to the bDrawBuffer type
    4 // static int bWeight;
    5 // static int bType;
    6 // static bDrawColor bColor;
    7 // static bDrawColor bColor_R;
    8 // static bDrawColor bColor_G;
    9 // static bDrawColor bColor_B;
    10 // static bDrawBuffer *bBuffer;
     3// buffer->pixels carries the plot image
     4// buffer->mask is 0 if the pixel is untouched, 1 if the pixel has data
    115
    126void bDrawCircleSingle (bDrawBuffer *buffer, double xc, double yc, double radius);
     7
     8int bDrawMerge (bDrawBuffer *base, bDrawBuffer *layer) {
     9
     10  // the two bDrawBuffers must match in size and depth
     11
     12  if (base->Nx != layer->Nx) return FALSE;
     13  if (base->Ny != layer->Ny) return FALSE;
     14  if (base->Nbyte != layer->Nbyte) return FALSE;
     15
     16  for (int j = 0; j < base->Ny; j++) {
     17    for (int i = 0; i < base->Nx; i++) {
     18
     19      if (!layer->mask[j][i]) continue;
     20     
     21      // completely opaque top layer:
     22      if (base->Nbyte == 1) {
     23        base->pixels[j][i] = layer->pixels[j][i];
     24      } else {
     25        base->pixels[j][3*i+0] = layer->pixels[j][3*i+0];
     26        base->pixels[j][3*i+1] = layer->pixels[j][3*i+1];
     27        base->pixels[j][3*i+2] = layer->pixels[j][3*i+2];
     28      }
     29    }
     30  }
     31  return TRUE;
     32}
    1333
    1434// create a drawing buffer with either 1 or 3 byte colors
    … …  
    3555
    3656  ALLOCATE (buffer[0].pixels, bDrawColor *, Ny);
     57  ALLOCATE (buffer[0].mask, char *, Ny);
    3758  for (i = 0; i < Ny; i++) {
    3859    ALLOCATE (buffer[0].pixels[i], bDrawColor, Nbyte*Nx);
     60    ALLOCATE (buffer[0].mask[i], char, Nx);
    3961    for (j = 0; j < Nx; j++) {
    4062      if (Nbyte == 1) {
    … …  
    4567        buffer[0].pixels[i][3*j+2] = white_B;
    4668      }
     69      buffer[0].mask[i][j] = 0; // for now: 0 = no data, 1 = data
    4770    }
    4871  }
    … …  
    6386  for (i = 0; i < buffer[0].Ny; i++) {
    6487    free (buffer[0].pixels[i]);
     88    free (buffer[0].mask[i]);
    6589  }
    6690  free (buffer[0].pixels);
    67   free (buffer[0].palette);
     91  free (buffer[0].mask);
     92  // free (buffer[0].palette);
    6893  free (buffer);
    6994  return;
    … …  
    121146    buffer[0].pixels[y][x] = buffer->bColor;
    122147  } else {
    123     // buffer[0].pixels[y][3*x+0] = buffer->bColor_R;
    124     // buffer[0].pixels[y][3*x+1] = buffer->bColor_G;
    125     // buffer[0].pixels[y][3*x+2] = buffer->bColor_B;
    126     buffer[0].pixels[y][3*x+0] = MAX (0x00, MIN(0xff, beta * buffer[0].pixels[y][3*x+0] + alpha * buffer->bColor_R));
     148    buffer[0].pixels[y][3*x+0] = MAX (0x00, MIN(0xff, beta * buffer[0].pixels[y][3*x+0] + alpha * buffer->bColor_R)); // was just buffer->bColor_R, etc
    127149    buffer[0].pixels[y][3*x+1] = MAX (0x00, MIN(0xff, beta * buffer[0].pixels[y][3*x+1] + alpha * buffer->bColor_G));
    128150    buffer[0].pixels[y][3*x+2] = MAX (0x00, MIN(0xff, beta * buffer[0].pixels[y][3*x+2] + alpha * buffer->bColor_B));
    129151  }
     152  buffer[0].mask[y][x] = 1;
    130153  return;
    131154}
    … …  
    150173void bDrawRectOpen (bDrawBuffer *buffer, double x1, double y1, double x2, double y2) {
    151174
    152   int X1, Y1, X2, Y2;
    153 
    154175  if (x1 > x2) SWAP (x1, x2);
    155176  if (y1 > y2) SWAP (y1, y2);
    156177
     178  int X1 = MIN (MAX (ROUND (x1), 0), buffer[0].Nx - 1);
     179  int X2 = MIN (MAX (ROUND (x2), 1), buffer[0].Nx - 1);
     180
     181  int Y1 = MIN (MAX (ROUND (y1), 0), buffer[0].Ny - 1);
     182  int Y2 = MIN (MAX (ROUND (y2), 1), buffer[0].Ny - 1);
     183
     184  int dNs = -0.5*(buffer->bWeight - 1);
     185  /* 0, 0, 0, -1, -1, -2, -2 */
     186
     187  int dNe = +0.5*buffer->bWeight + 1;
     188  /* 1, 1, 2, 2, 2, 3, 3 */
     189
     190  for (int dN = dNs; dN < dNe; dN ++) {
     191    // line on the bottom needs to run longer for the negative dN values
     192    bDrawLineHorizontal (buffer, X1 + dN, X2 + 1 - dN, Y1 + dN);
     193
     194    // line on the top needs to run longer for the positive dN values
     195    bDrawLineHorizontal (buffer, X1 + dN, X2 + 1 - dN, Y2 - dN);
     196
     197    // line on the left needs to run longer for the negative dN values
     198    bDrawLineVertical   (buffer, X1 + dN, Y1 + dN, Y2 - dN);
     199
     200    // line on the right needs to run longer for the positive dN values
     201    bDrawLineVertical   (buffer, X2 - dN, Y1 + dN, Y2 - dN);
     202  }
     203  return;
     204}
     205
     206void bDrawRectFill (bDrawBuffer *buffer, double x1, double y1, double x2, double y2) {
     207
     208  int i;
     209  int X1, Y1, X2, Y2;
     210
     211  if (x1 > x2) SWAP (x1, x2);
     212  if (y1 > y2) SWAP (y1, y2);
     213
    157214  X1 = MIN (MAX (ROUND (x1), 0), buffer[0].Nx - 1);
    158215  X2 = MIN (MAX (ROUND (x2), 1), buffer[0].Nx - 1);
    159216
    160217  Y1 = MIN (MAX (ROUND (y1), 0), buffer[0].Ny - 1);
    161   Y2 = MIN (MAX (ROUND (y2), 1), buffer[0].Ny - 1);
    162 
    163   bDrawLineHorizontal (buffer, X1, X2 + 1, Y1);
    164   bDrawLineHorizontal (buffer, X1, X2 + 1, Y2);
    165   bDrawLineVertical   (buffer, X1, Y1, Y2);
    166   bDrawLineVertical   (buffer, X2, Y1, Y2);
    167   return;
    168 }
    169 
    170 void bDrawRectFill (bDrawBuffer *buffer, double x1, double y1, double x2, double y2) {
    171 
    172   int i;
    173   int X1, Y1, X2, Y2;
    174 
    175   if (x1 > x2) SWAP (x1, x2);
    176   if (y1 > y2) SWAP (y1, y2);
    177 
    178   X1 = MIN (MAX (ROUND (x1), 0), buffer[0].Nx - 1);
    179   X2 = MIN (MAX (ROUND (x2), 1), buffer[0].Nx - 1);
    180 
    181   Y1 = MIN (MAX (ROUND (y1), 0), buffer[0].Ny - 1);
    182218  Y2 = MIN (MAX (ROUND (y2), 0), buffer[0].Ny - 1);
    183219
    184220  for (i = Y1; i < Y2; i++) {
     221    // this should be a line of width 1 or we duplicate pixels
    185222    bDrawLineHorizontal (buffer, X1, X2, i);
    186223  }
    … …  
    235272// use the Bresenham line drawing technique
    236273// integer-only Bresenham line-draw version which is fast
     274// bresenham assumes x1 < x2 and (y2 - y1) < (x2 - x1)
    237275void bDrawLineBresen (bDrawBuffer *buffer, int X1, int Y1, int X2, int Y2, int swapcoords) {
    238276
    … …  
    299337      buffer[0].pixels[Y][3*i+1] = MAX (0x00, MIN(0xff, beta * buffer[0].pixels[Y][3*i+1] + alpha * buffer->bColor_G));
    300338      buffer[0].pixels[Y][3*i+2] = MAX (0x00, MIN(0xff, beta * buffer[0].pixels[Y][3*i+2] + alpha * buffer->bColor_B));
    301       // buffer[0].pixels[Y][3*i+0] = buffer->bColor_R;
    302       // buffer[0].pixels[Y][3*i+1] = buffer->bColor_G;
    303       // buffer[0].pixels[Y][3*i+2] = buffer->bColor_B;
    304     }
     339    }
     340    buffer[0].mask[Y][i] = 1;
    305341  }
    306342  return;
    … …  
    326362      buffer[0].pixels[i][3*X+1] = MAX (0x00, MIN(0xff, beta * buffer[0].pixels[i][3*X+1] + alpha * buffer->bColor_G));
    327363      buffer[0].pixels[i][3*X+2] = MAX (0x00, MIN(0xff, beta * buffer[0].pixels[i][3*X+2] + alpha * buffer->bColor_B));
    328       // buffer[0].pixels[i][3*X+0] = buffer->bColor_R;
    329       // buffer[0].pixels[i][3*X+1] = buffer->bColor_G;
    330       // buffer[0].pixels[i][3*X+2] = buffer->bColor_B;
    331     }
     364    }
     365    buffer[0].mask[i][X] = 1;
    332366  }
    333367  return;
    … …  
    384418 
    385419  return;
     420}
     421
     422# define BOT_LEFT 0
     423# define BOT_RGHT 1
     424# define TOP_LEFT 2
     425# define TOP_RGHT 3
     426
     427// I should probably look up an appropriate recipe for this
     428// source code for XDrawPolyFill?
     429
     430int GetNextSeqNumber (int seq, int Npoint, int Clockwise, int LeftSide);
     431int bDrawFillBetweenSegments (bDrawBuffer *buffer, int *x, int *y, int ystart);
     432
     433void bDrawPolyFill (bDrawBuffer *buffer, double *x, double *y, int Npoints) {
     434
     435  // The coord list which is passed in (x,y) must be in order around the contour
     436  for (int i = 0; i < Npoints; i++) {
     437    fprintf (stderr, "%d : %f,%f\n", i, x[i], y[i]);
     438  }
     439
     440  // First, find the lowest point and start at that sequence number
     441  int iMin = 0;
     442  double yMin = y[iMin];
     443  for (int i = 1; i < Npoints; i++) {
     444    if (y[i] < yMin) { iMin = i; yMin = y[iMin]; }
     445  }
     446
     447  // we generate two line segments and will fill between them
     448  int xval[4];
     449  int yval[4];
     450
     451  // the starting point is a vertex
     452  xval[BOT_LEFT] = x[iMin];  yval[BOT_LEFT] = y[iMin];
     453  xval[BOT_RGHT] = x[iMin];  yval[BOT_RGHT] = y[iMin];
     454
     455  // get iNext, iPrev assuming ClockWise, then test
     456  int isCW = TRUE;
     457  int iNextLeft = GetNextSeqNumber (iMin, Npoints, isCW, TRUE); // TRUE => left-side
     458  int iNextRght = GetNextSeqNumber (iMin, Npoints, isCW, FALSE);
     459
     460  // if this is true, the points are not clockwise
     461  if (x[iNextLeft] > x[iNextRght]) {
     462    isCW = FALSE;
     463    int tmp = iNextLeft;
     464    iNextLeft = iNextRght;
     465    iNextRght = tmp;
     466  }
     467
     468  int iTL = iNextLeft; // sequence number of the top,left point
     469  int iTR = iNextRght; // sequence number of the top,right point
     470  xval[TOP_LEFT] = x[iTL]; yval[TOP_LEFT] = y[iTL];
     471  xval[TOP_RGHT] = x[iTR]; yval[TOP_RGHT] = y[iTR];
     472
     473  // we need to track the starting row, start at the bottom
     474  int ystart = yval[BOT_LEFT];
     475
     476  while (TRUE) {
     477    int isLeft = bDrawFillBetweenSegments (buffer, xval, yval, ystart);
     478
     479    // the last pair of segments in the sequence end at the same top point:
     480    if (iTL == iTR) break;
     481
     482    // after one pass, if isLeft is true, then cycle the points on the left
     483    // segment, otherwise cycle the right segment:
     484
     485    if (isLeft) {
     486      iTL = GetNextSeqNumber (iTL, Npoints, isCW, TRUE);
     487      xval[BOT_LEFT] = xval[TOP_LEFT]; yval[BOT_LEFT] = yval[TOP_LEFT];
     488      xval[TOP_LEFT] = x[iTL];         yval[TOP_LEFT] = y[iTL];
     489      ystart = yval[BOT_LEFT];
     490    } else {
     491      iTR = GetNextSeqNumber (iTR, Npoints, isCW, FALSE);
     492      xval[BOT_RGHT] = xval[TOP_RGHT]; yval[BOT_RGHT] = yval[TOP_RGHT];
     493      xval[TOP_RGHT] = x[iTR];         yval[TOP_RGHT] = y[iTR];
     494      ystart = yval[BOT_RGHT];
     495    }
     496  }
     497  return;
     498}
     499
     500// we have two line segments defined by x[],y[]. 
     501
     502// [0],[1] is the left-side segment
     503// [2],[3] is the right-side segment
     504
     505// Fill the region between the two segments until
     506// the lower of the two y coordinates y[1], y[3]
     507
     508// return which of y[1] or y[3] was the ending row
     509
     510int bDrawFillBetweenSegments (bDrawBuffer *buffer, int *x, int *y, int ystart) {
     511
     512  int dY_L = y[TOP_LEFT] - y[BOT_LEFT];
     513  int dX_L = x[TOP_LEFT] - x[BOT_LEFT];
     514
     515  int dY_R = y[TOP_RGHT] - y[BOT_RGHT];
     516  int dX_R = x[TOP_RGHT] - x[BOT_RGHT];
     517   
     518  double Slope_L = (dY_L == 0) ? NAN : dX_L / (double) dY_L;
     519  double Slope_R = (dY_R == 0) ? NAN : dX_R / (double) dY_R;
     520
     521  // stop at whichever comes first
     522  // XXX not sure if I need to use < or <= here:
     523  // with <=, the value of Y after the loop will be one too high.
     524  int Y;
     525  for (Y = ystart; (Y < y[TOP_LEFT]) && (Y < y[TOP_RGHT]); Y++) {
     526
     527    // calculating X_L, X_R based on the slope and current y point
     528    int X_L = (dY_L == 0) ? x[BOT_LEFT] : (Y - y[BOT_LEFT]) * Slope_L + x[BOT_LEFT];
     529    int X_R = (dY_R == 0) ? x[BOT_RGHT] : (Y - y[BOT_RGHT]) * Slope_R + x[BOT_RGHT];
     530     
     531    // draw horizontal line from X_L to X_R at Y
     532    bDrawLineHorizontal (buffer, X_L, X_R, Y);
     533  }
     534
     535  if (Y == y[TOP_LEFT]) {
     536    return TRUE; // TRUE is LEFT SIDE
     537  } else {
     538    return FALSE;
     539  }
     540}
     541
     542int GetNextSeqNumber (int seq, int Npoints, int Clockwise, int LeftSide) {
     543
     544  int iNext;
     545
     546  if (Clockwise) {
     547    if (LeftSide) {
     548      // next point is next in sequence
     549      iNext = (seq < Npoints - 1) ? seq + 1 : 0;
     550    } else {
     551      // next point is prev in sequence
     552      iNext = (seq > 0) ? seq - 1 : Npoints - 1;
     553    }
     554  } else {
     555    if (LeftSide) {
     556      // next point is prev in sequence
     557      iNext = (seq > 0) ? seq - 1 : Npoints - 1;
     558    } else {
     559      // next point is next in sequence
     560      iNext = (seq < Npoints - 1) ? seq + 1 : 0;
     561    }
     562  }
     563  return iNext;
    386564}
    387565
    … …  
    585763  if (x == 0) {
    586764    if (decrementY) {
     765      // this should must be a line or width 1 or we duplicate pixels
    587766      bDrawLineHorizontal (buffer, Xc    , Xc + 1, Yc + y);
    588767      bDrawLineHorizontal (buffer, Xc    , Xc + 1, Yc - y);
    … …  
    590769
    591770    // center line
     771    // this should must be a line or width 1 or we duplicate pixels
    592772    bDrawLineHorizontal (buffer, Xc - y, Xc + y + 1, Yc    );
    593773    return;
    … …  
    595775
    596776  if (x == y) {
     777    // this should must be a line or width 1 or we duplicate pixels
    597778    bDrawLineHorizontal (buffer, Xc - x, Xc + x + 1, Yc + y);
    598779    bDrawLineHorizontal (buffer, Xc - x, Xc + x + 1, Yc - y);
    … …  
    602783  // only draw these two lines if we decrement y
    603784  if (decrementY) {
     785    // this must be a line or width 1 or we duplicate pixels
    604786    bDrawLineHorizontal (buffer, Xc - x, Xc + x + 1, Yc + y);
    605787    bDrawLineHorizontal (buffer, Xc - x, Xc + x + 1, Yc - y);
    … …  
    607789
    608790  // always draw these two lines:
     791  // this must be a line or width 1 or we duplicate pixels
    609792  bDrawLineHorizontal (buffer, Xc - y, Xc + y + 1, Yc + x);
    610793  bDrawLineHorizontal (buffer, Xc - y, Xc + y + 1, Yc - x);
    … …  
    648831}
    649832
     833// run a smoothing kernel on each channel for anti-aliasing
     834// XXX I need to consider the mask as well
     835void bDrawSmooth (bDrawBuffer *buffer, float sigma) {
     836
     837  // generate a temp buffer for storage of the smoothed result
     838  int Nx = buffer[0].Nx;
     839  int Ny = buffer[0].Ny;
     840  ALLOCATE_PTR (temp, bDrawColor, Nx*Ny);
     841
     842  //* build a 1D gaussian
     843  int Nsigma = 2;
     844  int Ns = (int) (Nsigma*sigma + 0.5);
     845  int Ngauss = 2*Ns + 1;
     846  ALLOCATE_PTR (gaussnorm, float, Ngauss);
     847  float *gauss = &gaussnorm[Ns];
     848  for (int i = -Ns; i < Ns + 1; i++) {
     849    gauss[i] = exp ((i*i)/(-2*sigma*sigma));
     850  }
     851
     852  // NOTE XXX: need to handle Nbyte = 1 or 3
     853  int Nchannel = 3;
     854
     855  // we need to smooth the 3 channels independently, do this as three passes
     856  for (int ch = 0; ch < Nchannel; ch++) {
     857
     858    // smooth in X direction
     859    for (int j = 0; j < Ny; j++) {
     860      bDrawColor *vi = (bDrawColor *) &buffer[0].pixels[j][ch];
     861      bDrawColor *vo = &temp[j*Nx];
     862      for (int i = 0; i < Nx; i++, vo++, vi += Nchannel) {
     863        float g = 0.0, s = 0.0;
     864        for (int n = -Ns; n <= Ns; n++) {
     865          if (i+n < 0) continue;
     866          if (i+n >= Nx) continue;
     867          s += gauss[n]*vi[Nchannel*n];
     868          g += gauss[n];
     869        }
     870        *vo = s / g;
     871      }
     872    }
     873
     874    // NOTE: output maps back into original image
     875    // careful on the counting
     876
     877    /* smooth in Y direction */
     878    for (int i = 0; i < Nx; i++) {
     879      bDrawColor *vi = &temp[i];
     880      for (int j = 0; j < Ny; j++) {
     881        float g = 0.0, s = 0.0;
     882        for (int n = -Ns; n <= Ns; n++) {
     883          if (j+n < 0) continue;
     884          if (j+n >= Ny) continue;
     885          s += gauss[n]*vi[(n+j)*Nx]; // vi is a single-index array, so this is stepping by rows
     886          g += gauss[n];
     887        }
     888        buffer[0].pixels[j][Nchannel*i + ch] = s / g;
     889      }
     890    }
     891  }
     892
     893  // smooth the mask
     894  ALLOCATE_PTR (temp_mask, float, Nx*Ny);
     895
     896  // smooth in X direction
     897  for (int j = 0; j < Ny; j++) {
     898    char *vi = (char *) &buffer[0].mask[j][0];
     899    float *vo = &temp_mask[j*Nx];
     900    for (int i = 0; i < Nx; i++, vo++, vi++) {
     901        float g = 0.0, s = 0.0;
     902        for (int n = -Ns; n <= Ns; n++) {
     903          if (i+n < 0) continue;
     904          if (i+n >= Nx) continue;
     905          s += gauss[n]*vi[n];
     906          g += gauss[n];
     907        }
     908        *vo = s / g;
     909    }
     910  }
     911
     912  // NOTE: output maps back into original image
     913  // careful on the counting
     914
     915  /* smooth in Y direction */
     916  for (int i = 0; i < Nx; i++) {
     917    float *vi = &temp_mask[i];
     918    for (int j = 0; j < Ny; j++) {
     919        float g = 0.0, s = 0.0;
     920        for (int n = -Ns; n <= Ns; n++) {
     921          if (j+n < 0) continue;
     922          if (j+n >= Ny) continue;
     923          s += gauss[n]*vi[(n+j)*Nx]; // vi is a single-index array, so this is stepping by rows
     924          g += gauss[n];
     925        }
     926        buffer[0].mask[j][i] = (s / g > 0.05) ? 1 : 0;
     927    }
     928  }
     929
     930  free (temp);
     931  free (temp_mask);
     932  free (gaussnorm);
     933
     934  return;
     935}
     936
    650937/*
     938NOTES on Bresenham-style circles:
    651939the discriminant of inside or outside the circle is:
    652940

  • branches/eam_branches/ipp-20191011/Ohana/src/libohana/include/ohana_allocate.h

    r39457 r41170  
    4444void ohana_memcheck_block_func (const char *myFile, int myLine, const char *myFunc, void *in);
    4545void  real_free (void *in);
     46
     47void ohana_memdump_strings_file (FILE *f, int VERBOSE);
     48void ohana_memdump_set_maxlines (int N);
    4649
    4750# define ohana_memcheck(X) ohana_memcheck_func (__FILE__, __LINE__, __func__, X)

  • branches/eam_branches/ipp-20191011/Ohana/src/libohana/include/ohana_sort.h

    r40291 r41170  
    8383void isortfour (int *X, int *Y, int *Z, int *W, int N);
    8484
     85void dsort_indexonly (double *X, off_t *S, off_t N);
     86void dsort_int_indexonly (double *X, int *S, int N);
     87
    8588#endif

  • branches/eam_branches/ipp-20191011/Ohana/src/libohana/src/ohana_allocate.c

    r40117 r41170  
    4545
    4646static pthread_mutex_t memBlockListMutex = PTHREAD_MUTEX_INITIALIZER;
     47
     48static int NblockMaxDump = 100;
    4749
    4850void ohana_meminit () {
    … …  
    435437        fprintf (stderr, "memory corruption\n");
    436438      }
    437       if (Nbad < 100) {
     439      if (Nbad < NblockMaxDump) {
    438440        fprintf (stderr, "  file: %s, line: %d, func: %s\n", thisBlock->file, thisBlock->line, thisBlock->func);
    439441      }
    … …  
    613615    Nbytes += thisBlock->size;
    614616
    615     if (Ntotal < 100) {
     617    if (Ntotal < NblockMaxDump) {
    616618      if (good) {
    617619        fprintf (f, "  %zd  %zd  %zd  GOOD  %s %d, func: %s\n", Ntotal, thisBlock->size, Nbytes, thisBlock->file, thisBlock->line, thisBlock->func);
    … …  
    677679  return memstats;
    678680}
     681
     682void ohana_memdump_strings_file (FILE *f, int VERBOSE) {
     683
     684  if (!lastBlock) {
     685    if (VERBOSE) fprintf (f, "no memory allocated\n");
     686    return;
     687  }
     688
     689  OhanaMemblock *thisBlock = lastBlock;
     690
     691  size_t Ntotal = 0;
     692  size_t Nbytes = 0;
     693  size_t Nstring = 0;
     694
     695  fprintf (f, "    entry |    bytes | cumulative : line : string\n");
     696
     697  while (thisBlock) {
     698
     699    // pointer to the start of the user memory
     700    char *ptr = (char *)(thisBlock + 1);
     701
     702    Ntotal ++;
     703    Nbytes += thisBlock->size;
     704
     705    int N = strlen(thisBlock->file);
     706    if (N > 22) {
     707      int found = !strcmp("libohana/src/string.c", &thisBlock->file[N - 21]);
     708      if (found) {
     709        if (Nstring < NblockMaxDump) {
     710          fprintf (f, " %8zd | %8zd | %10zd : %4d : %s\n", Ntotal, thisBlock->size, Nbytes, thisBlock->line, ptr);
     711        }
     712        Nstring ++;
     713      }
     714      thisBlock = thisBlock->prevBlock;
     715    }
     716  }
     717
     718  if (Nstring || VERBOSE) {
     719    fprintf (f, "%zd strings allocated\n", Nstring);
     720  }
     721
     722  return;
     723}
     724
     725void ohana_memdump_set_maxlines (int N) {
     726  NblockMaxDump = N;
     727}

  • branches/eam_branches/ipp-20191011/Ohana/src/libohana/src/sorts.c

    r40291 r41170  
    196196}
    197197
     198/* sort the index of a vector (vector stays unsorted) */
     199void dsort_indexonly (double *X, off_t *S, off_t N) {
     200 
     201# define SWAPFUNC(A,B){ off_t itmp; \
     202  itmp = S[A]; S[A] = S[B]; S[B] = itmp; \
     203}
     204# define COMPARE(A,B)((!isfinite(X[S[A]]) && isfinite(X[S[B]])) || (X[S[A]] < X[S[B]]))
     205
     206  OHANA_SORT (N, COMPARE, SWAPFUNC);
     207
     208# undef SWAPFUNC
     209# undef COMPARE
     210
     211}
     212
     213/* sort the index of a vector (vector stays unsorted) */
     214void dsort_int_indexonly (double *X, int *S, int N) {
     215 
     216# define SWAPFUNC(A,B){ int itmp; \
     217  itmp = S[A]; S[A] = S[B]; S[B] = itmp; \
     218}
     219# define COMPARE(A,B)((!isfinite(X[S[A]]) && isfinite(X[S[B]])) || (X[S[A]] < X[S[B]]))
     220
     221  OHANA_SORT (N, COMPARE, SWAPFUNC);
     222
     223# undef SWAPFUNC
     224# undef COMPARE
     225
     226}
     227
     228

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.astro/cdensify.c

    r39294 r41170  
    3939  }
    4040
     41  int binning = 1;
     42  if ((N = get_argument (argc, argv, "-binning"))) {
     43    remove_argument (N, &argc, argv);
     44    binning = atoi(argv[N]);
     45    remove_argument (N, &argc, argv);
     46  }
     47
    4148  int PSFTYPE = IS_DOT;
    4249  if ((N = get_argument (argc, argv, "-psf"))) {
    … …  
    8087  Ymax = graphmode.ymax;
    8188  Ymin = graphmode.ymin;
    82   dX = (Xmax - Xmin) / (Xpix - 1);
    83   dY = (Ymax - Ymin) / (Ypix - 1);
     89
     90  // (dX, dY) are the pixel scale of the output image (output pixels / input pixels)
     91
     92  dX = binning * (Xmax - Xmin) / (Xpix - 1);
     93  dY = binning * (Ymax - Ymin) / (Ypix - 1);
    8494
    8595  CHECKVAL(Xmin);

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.astro/cdhistogram.c

    r39233 r41170  
    1313  int kapa;
    1414  Graphdata graphmode;
     15
     16  int binning = 1;
     17  if ((N = get_argument (argc, argv, "-binning"))) {
     18    remove_argument (N, &argc, argv);
     19    binning = atoi(argv[N]);
     20    remove_argument (N, &argc, argv);
     21  }
    1522
    1623  range = NULL;
    … …  
    7986  Ymax = graphmode.ymax;
    8087  Ymin = graphmode.ymin;
    81   dX = (Xmax - Xmin) / (Xpix - 1);
    82   dY = (Ymax - Ymin) / (Ypix - 1);
     88
     89  // (dX, dY) are the pixel scale of the output image (output pixels / input pixels)
     90
     91  dX = binning * (Xmax - Xmin) / (Xpix - 1);
     92  dY = binning * (Ymax - Ymin) / (Ypix - 1);
    8393
    8494  CHECKVAL(Xmin);
    … …  
    90100  CHECKVAL(dY);
    91101
    92   int Nx = (Xmax - Xmin) / dX + 1;
    93   int Ny = (Ymax - Ymin) / dY + 1;
     102  int Nx = abs((Xmax - Xmin) / dX) + 1;
     103  int Ny = abs((Ymax - Ymin) / dY) + 1;
    94104 
    95105  Coords newcoords = graphmode.coords;

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.astro/region.c

    r39225 r41170  
    33int region (int argc, char **argv) {
    44 
    5   double Ra, Dec, Radius;
     5  double Ra, Dec;
    66  float dx, dy;
    77  int N, kapa, NoClear, dXpix, dYpix;
    … …  
    5454  }
    5555
     56  float XSIZE = NAN;
     57  if ((N = get_argument (argc, argv, "-xsize"))) {
     58    remove_argument (N, &argc, argv);
     59    XSIZE = atof (argv[N]);
     60    remove_argument (N, &argc, argv);
     61    if (XSIZE <= 0) {
     62      gprint (GP_ERR, "ERROR: xsize <= 0\n");
     63      return FALSE;
     64    }
     65  }
     66
     67  float YSIZE = NAN;
     68  if ((N = get_argument (argc, argv, "-ysize"))) {
     69    remove_argument (N, &argc, argv);
     70    YSIZE = atof (argv[N]);
     71    remove_argument (N, &argc, argv);
     72    if (YSIZE <= 0) {
     73      gprint (GP_ERR, "ERROR: xsize <= 0\n");
     74      return FALSE;
     75    }
     76  }
     77
    5678  if ((N = get_argument (argc, argv, "-ew"))) {
    5779    remove_argument (N, &argc, argv);
    … …  
    81103  }
    82104
    83   if ((argc != 4) && (argc != 5)) {
    84     gprint (GP_ERR, "USAGE: region Ra Dec Radius [projection] [orientation]\n");
     105  int RadiusArg, CtypeArg;
     106  if (!isnan(XSIZE) || !isnan(YSIZE)) {
     107    RadiusArg = -1;
     108    CtypeArg = 3;
     109    if ((argc != 3) && (argc != 4)) {
     110    gprint (GP_ERR, "USAGE: region Ra Dec [projection] [-xsize deg] [-ysize deg]\n");
     111    gprint (GP_ERR, "   OR: region Ra Dec Radius [projection]\n");
    85112    gprint (GP_ERR, "  [-image] [-ew] [+ew] [-ns] [+ns] [-no-clear] [-angle theta]\n");
    86113    gprint (GP_ERR, " current: %f %f (%f x %f) (%s)\n",
    … …  
    90117             &graphmode.coords.ctype[5]);
    91118    return (FALSE);
    92   }
    93  
     119    }
     120  } else {
     121    RadiusArg = 3;
     122    CtypeArg = 4;
     123    if ((argc != 4) && (argc != 5)) {
     124      gprint (GP_ERR, "USAGE: region Ra Dec Radius [projection]\n");
     125      gprint (GP_ERR, "   OR: region Ra Dec [projection] [-xsize deg] [-ysize deg]\n");
     126      gprint (GP_ERR, "  [-image] [-ew] [+ew] [-ns] [+ns] [-no-clear] [-angle theta]\n");
     127      gprint (GP_ERR, " current: %f %f (%f x %f) (%s)\n",
     128              graphmode.coords.crval1, graphmode.coords.crval2,
     129              fabs(graphmode.xmax - graphmode.xmin),
     130              fabs(graphmode.ymax - graphmode.ymin),
     131              &graphmode.coords.ctype[5]);
     132      return (FALSE);
     133    }
     134  } 
    94135  if (!ohana_str_to_radec (&Ra, &Dec, argv[1], argv[2])) return (FALSE);
    95   Radius = atof (argv[3]);
     136
     137  // region 0 0 sin -- should raise an error (radius = 0 or non-numeric)
     138
     139  // I want to be able to support the old style call in which both of these were valid:
     140  // region 0 0 90 ait
     141  // region 0 0 90      <- uses existing, sticky projection type
     142  // but I also want to be able to use
     143  // region 0 0 ait -xsize 5
     144  // region 0 0 ait -ysize 5
     145
     146  double Radius = NAN;
     147  if (RadiusArg >= 0) {
     148    Radius = atof (argv[RadiusArg]);
     149  }
     150
    96151  InitCoords (&graphmode.coords, "DEC--TAN");
    97 
    98   if (argc == 5) {
    99     if (!strcasecmp (argv[4], "TAN"))
    100       strcpy (graphmode.coords.ctype, "DEC--TAN");
    101     if (!strcasecmp (argv[4], "SIN"))
    102       strcpy (graphmode.coords.ctype, "DEC--SIN");
    103     if (!strcasecmp (argv[4], "ARC"))
    104       strcpy (graphmode.coords.ctype, "DEC--ARC");
    105     if (!strcasecmp (argv[4], "STG"))
    106       strcpy (graphmode.coords.ctype, "DEC--STG");
    107     if (!strcasecmp (argv[4], "ZEA"))
    108       strcpy (graphmode.coords.ctype, "DEC--ZEA");
    109     if (!strcasecmp (argv[4], "AIT"))
    110       strcpy (graphmode.coords.ctype, "DEC--AIT");
    111     if (!strcasecmp (argv[4], "GLS"))
    112       strcpy (graphmode.coords.ctype, "DEC--GLS");
    113     if (!strcasecmp (argv[4], "PAR"))
    114       strcpy (graphmode.coords.ctype, "DEC--PAR");
    115   }
     152  if (argc == CtypeArg + 1) {
     153    if (!strcasecmp (argv[CtypeArg], "TAN")) { strcpy (graphmode.coords.ctype, "DEC--TAN"); goto got_ctype; }
     154    if (!strcasecmp (argv[CtypeArg], "SIN")) { strcpy (graphmode.coords.ctype, "DEC--SIN"); goto got_ctype; }
     155    if (!strcasecmp (argv[CtypeArg], "ARC")) { strcpy (graphmode.coords.ctype, "DEC--ARC"); goto got_ctype; }
     156    if (!strcasecmp (argv[CtypeArg], "STG")) { strcpy (graphmode.coords.ctype, "DEC--STG"); goto got_ctype; }
     157    if (!strcasecmp (argv[CtypeArg], "ZEA")) { strcpy (graphmode.coords.ctype, "DEC--ZEA"); goto got_ctype; }
     158    if (!strcasecmp (argv[CtypeArg], "AIT")) { strcpy (graphmode.coords.ctype, "DEC--AIT"); goto got_ctype; }
     159    if (!strcasecmp (argv[CtypeArg], "GLS")) { strcpy (graphmode.coords.ctype, "DEC--GLS"); goto got_ctype; }
     160    if (!strcasecmp (argv[CtypeArg], "PAR")) { strcpy (graphmode.coords.ctype, "DEC--PAR"); goto got_ctype; }
     161    gprint (GP_ERR, "ERROR: invalid projection type %s\n", argv[CtypeArg]);
     162    gprint (GP_ERR, "allowed values: TAN, SIN, ARC, STG, ZEA, AIT, GLS, PAR\n");
     163    return FALSE;
     164  }
     165got_ctype:
    116166 
    117167  graphmode.coords.crval1 = Ra;
    … …  
    126176  graphmode.coords.pc2_2 =  cos(Angle*RAD_DEG)*pc2_2;
    127177
    128   /* ask kapa for coordinate limits, to get the right aspect ratio */
     178  // ask kapa for coordinate limits, to get the right aspect ratio
     179  // dx, dy are the size of the graph region in pixels
    129180  KapaGetLimits (kapa, &dx, &dy);
    130181  dx = fabs (dx);
    … …  
    132183
    133184  /* define limits for Ra, Dec at center, grid in degrees */
    134   if (dy < dx) {
    135     graphmode.xmin = -(dx/dy)*Radius;
    136     graphmode.ymin = -Radius;
    137     graphmode.xmax = (dx/dy)*Radius;
    138     graphmode.ymax = Radius;
     185  if (RadiusArg >= 0) {
     186    // force non-anamorphic projection with Radius set to smaller axis
     187    if (dy < dx) {
     188      graphmode.xmin = -(dx/dy)*Radius;
     189      graphmode.ymin = -Radius;
     190      graphmode.xmax = (dx/dy)*Radius;
     191      graphmode.ymax = Radius;
     192    } else {
     193      graphmode.xmin = -Radius;
     194      graphmode.ymin = -(dy/dx)*Radius;
     195      graphmode.xmax = Radius;
     196      graphmode.ymax = (dy/dx)*Radius;
     197    }
    139198  } else {
    140     graphmode.xmin = -Radius;
    141     graphmode.ymin = -(dy/dx)*Radius;
    142     graphmode.xmax = Radius;
    143     graphmode.ymax = (dy/dx)*Radius;
    144   }
     199    if (isnan(XSIZE)) {
     200      graphmode.xmin = -(dx/dy)*YSIZE/2.0;
     201      graphmode.ymin = -YSIZE/2.0;
     202      graphmode.xmax = (dx/dy)*YSIZE/2.0;
     203      graphmode.ymax = YSIZE/2.0;
     204    }
     205    if (isnan(YSIZE)) {
     206      graphmode.xmin = -XSIZE/2.0;
     207      graphmode.ymin = -(dy/dx)*XSIZE/2.0;
     208      graphmode.xmax = XSIZE/2.0;
     209      graphmode.ymax = (dy/dx)*XSIZE/2.0;
     210    }
     211    // anamorphic projection:
     212    if (!isnan(XSIZE) && !isnan(YSIZE)) {
     213      graphmode.xmin = -XSIZE/2.0;
     214      graphmode.ymin = -YSIZE/2.0;
     215      graphmode.xmax =  XSIZE/2.0;
     216      graphmode.ymax =  YSIZE/2.0;
     217    }
     218  }
    145219
    146220  set_variable ("XMIN", graphmode.xmin);

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.basic/init.c

    r40522 r41170  
    7373  {1, "??",            list_vars,          "list variables *"},
    7474  {1, "#",             nop,                "a NOP function"},
     75  {1, "##",            nop,                "a NOP function"},
     76  {1, "###",           nop,                "a NOP function"},
    7577  {1, "local",         local,              "define local variables"},
    7678  {1, "macro",         macro,              "deal with the macros *"},

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.basic/memory.c

    r14177 r41170  
    55 
    66  if (argc < 2) goto usage;
     7
     8  if (!strcasecmp (argv[1], "max-lines")) {
     9    if (argc != 3) goto usage;
     10    int MaxLines = atoi(argv[2]);
     11    ohana_memdump_set_maxlines (MaxLines);
     12    return TRUE;
     13  }
    714
    815  if (!strcasecmp (argv[1], "all")) {
    … …  
    4350  }
    4451
     52  if (!strncasecmp ("strings", argv[1], strlen(argv[1]))) {
     53    ohana_memdump_strings_file (stderr, FALSE);
     54    return (TRUE);
     55  }
     56
    4557usage:
    46   gprint (GP_ERR, "USAGE: memory (all/leaks)\n");
     58  gprint (GP_ERR, "USAGE: memory (mode) [-max-lines]\n");
     59  gprint (GP_ERR, "USAGE: memory (max-lines) (N)\n");
     60  gprint (GP_ERR, " mode options: all, leaks, check, checkfree, variables, vectors, buffers, macros, commands, strings\n");
    4761  return (FALSE);
    4862}

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.basic/test/if.sh

    r16894 r41170  
    320320 end
    321321end
     322
     323# check memleaks
     324macro memtest1
     325
     326 local a b i N
     327
     328 $i = 0
     329 $a = 1
     330
     331 memory check
     332 output /dev/null
     333 for i 0 1000
     334   if ($a == 1)
     335    echo "run"   
     336   end
     337 end   
     338 output stdout
     339 memory check
     340end
     341
     342# check memleaks
     343macro memtest2
     344
     345 local a b i N
     346
     347 $i = 0
     348 $a = 1
     349 $b = 2
     350
     351 memory check
     352 output /dev/null
     353 for i 0 1000
     354   if (($a == 1) && ($b == 2))
     355    echo "run"   
     356   end
     357 end   
     358 output stdout
     359 memory check
     360end
     361
     362# memory test using continue
     363macro memtest3
     364
     365 local a b i N
     366
     367 $i = 0
     368 $a = 1
     369 $b = 2
     370
     371 memory check
     372 output /dev/null
     373 for i 0 1000
     374   if (($a == 1) && ($b == 2))
     375    continue
     376    echo "run"   
     377   end
     378 end   
     379 output stdout
     380 memory check
     381end
     382
     383# memory test using continue
     384macro memtest4
     385
     386 local a b i N bool
     387
     388
     389 $i = 0
     390 $a = 1
     391 $b = 2
     392 $bool = 0
     393
     394 memory check
     395 output /dev/null
     396 for i 0 1000
     397   $bool = (($a == 1) && ($b == 2))
     398 end   
     399 output stdout
     400 memory check
     401end

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.basic/test/macro.sh

    r10311 r41170  
    77# Is the macro working?
    88macro test_prep
     9 $test_var1 = check1
     10 $test_var2 = check2
     11 $var_count = $0
     12end
     13
     14# Is the macro working?
     15macro test_local
     16 local test_var test_var2 var_count
     17
    918 $test_var1 = check1
    1019 $test_var2 = check2
    … …  
    4756 end
    4857end
     58
     59# Memory Test for macro
     60macro memtest2
     61
     62 $i = 0
     63 test_prep
     64
     65 memory check
     66 for i 0 10000
     67  test_prep
     68 end   
     69 memory check
     70 
     71end
     72
     73# Memory Test for macro
     74macro memtest3
     75
     76 $i = 0
     77 test_local
     78
     79 memory check
     80 for i 0 10000
     81  test_local
     82 end   
     83 memory check
     84 
     85end

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/Makefile

    r40637 r41170  
    1919$(SRC)/init.$(ARCH).o          \
    2020$(SRC)/accum.$(ARCH).o          \
     21$(SRC)/antialias.$(ARCH).o              \
    2122$(SRC)/applyfit1d.$(ARCH).o     \
    2223$(SRC)/applyfit2d.$(ARCH).o     \
    … …  
    7879$(SRC)/join.$(ARCH).o           \
    7980$(SRC)/jpeg.$(ARCH).o           \
     81$(SRC)/kapamemory.$(ARCH).o             \
    8082$(SRC)/kern.$(ARCH).o           \
    8183$(SRC)/keyword.$(ARCH).o        \
    … …  
    8991$(SRC)/lookup.$(ARCH).o \
    9092$(SRC)/matrix.$(ARCH).o \
     93$(SRC)/match1d.$(ARCH).o        \
    9194$(SRC)/match2d.$(ARCH).o        \
    9295$(SRC)/mkrgb.$(ARCH).o  \
    9396$(SRC)/mcreate.$(ARCH).o        \
    9497$(SRC)/medacc.$(ARCH).o \
     98$(SRC)/mgaussdev.$(ARCH).o      \
    9599$(SRC)/mget.$(ARCH).o           \
    96100$(SRC)/mget3d.$(ARCH).o         \
    … …  
    99103$(SRC)/medimage.$(ARCH).o       \
    100104$(SRC)/medimage_commands.$(ARCH).o \
     105$(SRC)/medimage_calc.$(ARCH).o \
    101106$(SRC)/mset.$(ARCH).o           \
    102107$(SRC)/needles.$(ARCH).o        \
    … …  
    169174$(SRC)/vgauss.$(ARCH).o            \
    170175$(SRC)/vlorentz.$(ARCH).o          \
     176$(SRC)/vsigmoid.$(ARCH).o          \
    171177$(SRC)/vellipse.$(ARCH).o          \
    172178$(SRC)/vmaxwell.$(ARCH).o          \
    … …  
    176182$(SRC)/vzload.$(ARCH).o            \
    177183$(SRC)/vpeaks.$(ARCH).o            \
     184$(SRC)/vtransitions.$(ARCH).o     \
    178185$(SRC)/vpop.$(ARCH).o              \
    179186$(SRC)/vroll.$(ARCH).o             \
    … …  
    182189$(SRC)/vsmooth.$(ARCH).o           \
    183190$(SRC)/vstats.$(ARCH).o            \
     191$(SRC)/virls.$(ARCH).o             \
     192$(SRC)/vwtmean.$(ARCH).o \
    184193$(SRC)/xsection.$(ARCH).o          \
    185194$(SRC)/vsh.$(ARCH).o               \

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/cut.c

    r33662 r41170  
    11# include "data.h"
    22
    3 enum {SUM, MEAN, MEDIAN};
     3enum {SUM, MEAN, MEDIAN, INNER, NGOOD};
     4
     5double cutstat (float *value, int Nvalue, int mode) {
     6
     7  double output = 0.0;
     8
     9  if (Nvalue == 0) return NAN;
     10
     11  if (mode == MEDIAN) {
     12    fsort (value, Nvalue);
     13    if (Nvalue % 2) {
     14      int Ncenter = Nvalue / 2;
     15      myAssert ((Ncenter >= 0) && (Ncenter < Nvalue), "cutstat");
     16      output = value[Ncenter];
     17    } else {
     18      int Ncenter = Nvalue / 2 - 1;
     19      myAssert ((Ncenter >= 0) && (Ncenter < Nvalue - 1), "cutstat");
     20      output = 0.5*(value[Ncenter] + value[Ncenter + 1]);
     21    }
     22    return output;
     23  }
     24  if (mode == INNER) {
     25    fsort (value, Nvalue);
     26
     27    int Ns = 0, Ne = 0;
     28    if (Nvalue % 2) {
     29      // for an odd number of points take the same number below
     30      // and above the center value
     31      int Ncenter = Nvalue / 2;
     32      int Nquarter = 0.25*Nvalue;
     33      Ns = Ncenter - Nquarter;
     34      Ne = Ncenter + Nquarter;
     35    } else {
     36      // for an even number, the middle lies between two points
     37      // take the same number below as above
     38      int Ncenter = (int)(Nvalue / 2) - 1;
     39      int Nquarter = 0.25*Nvalue;
     40      Ns = Ncenter     - Nquarter;
     41      Ne = Ncenter + 1 + Nquarter;
     42    }
     43    int Nv = 0;
     44    for (int i = Ns; i <= Ne; i++) {
     45      myAssert ((i >= 0) && (i < Nvalue), "cutstat");
     46      output += value[i];
     47      Nv ++;
     48    }
     49    output /= Nv;
     50    return output;
     51  }
     52  if (mode == NGOOD) {
     53    return Nvalue;
     54  }
     55  if ((mode == MEAN) || (mode == SUM)) {
     56    for (int i = 0; i < Nvalue; i++) {
     57      myAssert ((i >= 0) && (i < Nvalue), "cutstat");
     58      output += value[i];
     59    }
     60    if (mode == MEAN) { output /= Nvalue; }
     61    return output;
     62  }
     63  return NAN;
     64}
    465
    566int cut (int argc, char **argv) {
    667 
    768  int i, j, N, Nx, Ny, Mode;
    8   float *Vin, *Vbuf, value;
     69  float *Vin, *Vbuf;
    970  int sx, sy, nx, ny;
    1071  Vector *xvec, *yvec;
    … …  
    1273
    1374  Mode = SUM;
     75  if ((N = get_argument (argc, argv, "-sum"))) {
     76    remove_argument (N, &argc, argv);
     77    Mode = SUM;
     78  }
    1479  if ((N = get_argument (argc, argv, "-median"))) {
    1580    remove_argument (N, &argc, argv);
    … …  
    1984    remove_argument (N, &argc, argv);
    2085    Mode = MEAN;
     86  }
     87  if ((N = get_argument (argc, argv, "-inner"))) {
     88    remove_argument (N, &argc, argv);
     89    Mode = INNER;
     90  }
     91  if ((N = get_argument (argc, argv, "-ngood"))) {
     92    remove_argument (N, &argc, argv);
     93    Mode = NGOOD;
    2194  }
    2295
    … …  
    36109  Ny = buf[0].matrix.Naxis[1];
    37110
    38   if ((sx < 0) || (sy < 0) || (sx+nx > Nx) || (sy+ny > Ny)) {
    39     gprint (GP_ERR, "region out of range\n");
     111  // ny & nx do not need to be constrained by the buffer, but they do need to be sensible
     112  if ((nx < 0) || (ny < 0) || (nx > 1e8) || (ny > 1e8)) {
     113    gprint (GP_ERR, "warning : extraction size is crazy: %d,%d\n", nx, ny);
    40114    return (FALSE);
    41115  }
    … …  
    57131
    58132    for (i = 0; i < nx; i++) {
     133      // for out-of-range areas, set the output pixels to NAN
     134      if (i + sx < 0) {
     135        yvec[0].elements.Flt[i] = NAN;
     136        continue;
     137      }
     138      if (i + sx >= Nx) {
     139        yvec[0].elements.Flt[i] = NAN;
     140        continue;
     141      }
     142
    59143      /* accumulate values */
     144      // skip out-of-range areas in y
     145      if (sy < 0) {
     146        ny += sy;
     147        sy = 0;
     148      }
     149      if (sy > Ny) sy = Ny;
     150      if (sy + ny >= Ny) ny = Ny - sy;
     151
    60152      Vin = (float *)(buf[0].matrix.buffer) + sy*Nx + sx + i;
     153      int Npix = 0;
    61154      for (j = 0; j < ny; j++, Vin += Nx) {
    62         Vbuf[j] = *Vin;
    63       }
    64       /* apply stat of choice */
    65       if (Mode == MEDIAN) {
    66         fsort (Vbuf, ny);
    67         value = Vbuf[(int)(0.5*ny)];
    68       } else {
    69         value = 0;
    70         for (j = 0; j < ny; j++) {
    71           value += Vbuf[j];
    72         }
    73         if (Mode == MEAN) { value /= ny; }
    74       }
    75       yvec[0].elements.Flt[i] = value;
     155        if (!isfinite(*Vin)) continue;
     156        Vbuf[Npix] = *Vin;
     157        Npix ++;
     158      }
     159      yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode);
    76160    }
    77161    free (Vbuf);
    … …  
    89173
    90174    for (i = 0; i < ny; i++) {
     175      // for out-of-range areas, set the output pixels to NAN
     176      if (i + sy < 0) {
     177        yvec[0].elements.Flt[i] = NAN;
     178        continue;
     179      }
     180      if (i + sy >= Ny) {
     181        yvec[0].elements.Flt[i] = NAN;
     182        continue;
     183      }
     184
    91185      /* accumulate values */
     186      // skip out-of-range areas in x
     187      if (sx < 0) {
     188        nx += sx;
     189        sx = 0;
     190      }
     191      if (sx > Nx) sx = Nx;
     192      if (sx + nx >= Nx) nx = Nx - sx;
     193
    92194      Vin = (float *)(buf[0].matrix.buffer) + (sy + i)*Nx + sx;
     195      int Npix = 0;
    93196      for (j = 0; j < nx; j++, Vin ++) {
    94         Vbuf[j] = *Vin;
    95       }
    96       /* apply stat of choice */
    97       if (Mode == MEDIAN) {
    98         fsort (Vbuf, nx);
    99         value = Vbuf[(int)(0.5*nx)];
    100       } else {
    101         value = 0;
    102         for (j = 0; j < nx; j++) {
    103           value += Vbuf[j];
    104         }
    105         if (Mode == MEAN) { value /= nx; }
    106       }
    107       yvec[0].elements.Flt[i] = value;
     197        if (!isfinite(*Vin)) continue;
     198        Vbuf[Npix] = *Vin;
     199        Npix ++;
     200      }
     201      yvec[0].elements.Flt[i] = cutstat (Vbuf, Npix, Mode);
    108202    }
    109203    free (Vbuf);

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/extract.c

    r28241 r41170  
    11# include "data.h"
     2
     3/* <from> : source image, must exist
     4   <to>   : target image -- if it does not exist, it is created of size (Nx,Ny)
     5   sx, sy : source starting coordinate -- need not be on a valid image pixel
     6   nx, ny : number of source pixels -- currently must not go out of bounds -> allow out-of-bounds
     7   Sx, Sy : target starting coordinate -- need not be on a valid image pixel
     8   Nx, Ny : redundant information UNLESS <to> does exist (in which case it is required information) -> make Nx,Ny optional if <to> exists?
     9 */
    210
    311int extract (int argc, char **argv) {
    412 
    5   int i, j;
    6   float *Vin, *Vout;
    7   int sx, sy, nx, ny, NX, NY;
    8   int Sx, Sy, Nx, Ny;
     13  int N;
    914  Buffer *in, *out;
     15
     16  float initValue = 0.0;
     17  int initOutput = FALSE;
     18  if ((N = get_argument (argc, argv, "-init"))) {
     19    remove_argument (N, &argc, argv);
     20    initValue = atof (argv[N]);
     21    remove_argument (N, &argc, argv);
     22    initOutput = TRUE;
     23  }
    1024
    1125  if (argc != 11) {
    … …  
    1529
    1630  if ((in = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    17   NX = in[0].matrix.Naxis[0];
    18   NY = in[0].matrix.Naxis[1];
     31  int NX = in[0].matrix.Naxis[0];
     32  int NY = in[0].matrix.Naxis[1];
    1933
    20   sx = atof (argv[3]);
    21   sy = atof (argv[4]);
    22   nx = atof (argv[5]);
    23   ny = atof (argv[6]);
     34  int sx = atof (argv[3]);
     35  int sy = atof (argv[4]);
     36  int nx = atof (argv[5]);
     37  int ny = atof (argv[6]);
    2438
    25   Sx = atof (argv[7]);
    26   Sy = atof (argv[8]);
    27   Nx = atof (argv[9]);
    28   Ny = atof (argv[10]);
     39  int Sx = atof (argv[7]);
     40  int Sy = atof (argv[8]);
     41  int Nx = atof (argv[9]);
     42  int Ny = atof (argv[10]);
    2943
    3044  if ((Sy + ny > Ny) || (Sx + nx > Nx)) {
    31     gprint (GP_ERR, "mismatch between source and dest regions\n");
     45    gprint (GP_ERR, "source pixels extend beyond target pixels\n");
    3246    gprint (GP_ERR, "%d + %d > %d or %d + %d > %d\n", Sy, ny, Ny, Sx, nx, Nx);
    33     return (FALSE);
     47    // return (FALSE);
    3448  }
     49
     50  // XXX : allow source region to fall outside source image
    3551
    3652  /* region is not on first image */
    … …  
    3955      (sy > in[0].matrix.Naxis[1])) {
    4056    gprint (GP_ERR, "region outside of source image\n");
    41     return (FALSE);
     57    // return (FALSE);
    4258  }
    4359
    44   if ((Sx + nx > Nx) || (Sy + ny > Ny)) {
    45     gprint (GP_ERR, "source region larger than dest region\n");
    46     return (FALSE);
    47   }
    4860  if ((Sx < 0) || (Sy < 0)) {
    4961    gprint (GP_ERR, "dest region out of range\n");
    50     return (FALSE);
     62    // return (FALSE);
    5163  }
    5264
    … …  
    7688  }
    7789
    78   for (j = 0; j < ny; j++) {
    79     if (j + sy < 0) continue;
    80     if (j + sy >= NY) continue;
    81     Vin = (float *)(in[0].matrix.buffer) + (j + sy)*in[0].matrix.Naxis[0] + sx; 
    82     Vout = (float *)(out[0].matrix.buffer) + (j + Sy)*out[0].matrix.Naxis[0] + Sx;   
    83     for (i = 0; i < nx; i++, Vin++, Vout++) {
     90  // (NX,NY) : source image dimensions
     91  // (Nx,Ny) : target image dimensions
     92
     93  // allow (sx, sy), (Sx, Sy) to be off image
     94  // allow (sx+nx, sy+ny) to be off image
     95
     96  // if (sx < 0) {
     97  //   nx += sx;
     98  //   sx = 0;
     99  // }
     100  // if (sx > NX) sx = NX;
     101  //
     102  // if (Sx < 0) {
     103  //   Nx += sx;
     104  //   sx = 0;
     105  // }
     106  // if (sx > NX) sx = NX;
     107 
     108  float *Vin = (float *)(in[0].matrix.buffer);
     109  float *Vout = (float *)(out[0].matrix.buffer);
     110
     111  if (initOutput) {
     112    for (int j = 0; j < Ny; j++) {
     113      for (int i = 0; i < Nx; i++) {
     114        Vout[i + Nx*j] = initValue;
     115      }
     116    }
     117  }
     118
     119  int Nps = NX*NY;
     120  int Npt = Nx*Ny;
     121
     122  for (int j = 0; j < ny; j++) {
     123    if (j + sy < 0) continue; // not yet on source image
     124    if (j + Sy < 0) continue; // not yet on target image
     125    if (j + sy >= NY) continue; // past edge of source image
     126    if (j + Sy >= Ny) continue; // past edge of target image
     127
     128    // Vin = (float *)(in[0].matrix.buffer) + (j + sy)*in[0].matrix.Naxis[0] + sx; 
     129    // Vout = (float *)(out[0].matrix.buffer) + (j + Sy)*out[0].matrix.Naxis[0] + Sx;   
     130
     131    for (int i = 0; i < nx; i++) {
    84132      if (i + sx < 0) continue;
    85133      if (i + sx >= NX) continue;
    86       *Vout = *Vin;
     134      if (i + Sx < 0) continue;
     135      if (i + Sx >= Nx) continue;
     136
     137      int ps = i + sx + (j + sy)*NX;
     138      int pt = i + Sx + (j + Sy)*Nx;
     139
     140      myAssert (pt >= 0, "oops 1");
     141      myAssert (pt < Npt, "oops 2");
     142
     143      myAssert (ps >= 0, "oops 3");
     144      myAssert (ps < Nps, "oops 4");
     145
     146      Vout[pt] = Vin[ps];
    87147    }
    88148  }
    … …  
    91151
    92152}
    93 

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/init.c

    r40637 r41170  
    33
    44int accum            PROTO((int, char **));
     5int antialias        PROTO((int, char **));
    56int applyfit         PROTO((int, char **));
    67int applyfit1d       PROTO((int, char **));
    … …  
    6768int join             PROTO((int, char **));
    6869int jpeg             PROTO((int, char **));
     70int kapamemory       PROTO((int, char **));
    6971int kern             PROTO((int, char **));
    7072int keyword          PROTO((int, char **));
    … …  
    7981int lookup           PROTO((int, char **));
    8082int matrix           PROTO((int, char **));
     83int match1d          PROTO((int, char **));
    8184int match2d          PROTO((int, char **));
    8285int mkrgb            PROTO((int, char **));
    8386int mcreate          PROTO((int, char **));
    8487int medacc           PROTO((int, char **));
     88int mgaussdev        PROTO((int, char **));
    8589int mget             PROTO((int, char **));
    8690int mget3d           PROTO((int, char **));
    … …  
    156160int vgauss           PROTO((int, char **));
    157161int vlorentz         PROTO((int, char **));
     162int vsigmoid         PROTO((int, char **));
    158163int vellipse         PROTO((int, char **));
    159164int vmaxwell         PROTO((int, char **));
    … …  
    163168int vzload           PROTO((int, char **));
    164169int vstats           PROTO((int, char **));
     170int vstats           PROTO((int, char **));
     171int virls            PROTO((int, char **));
     172int vwtmean          PROTO((int, char **));
    165173int vroll            PROTO((int, char **));
    166174int vshift           PROTO((int, char **));
    167175int vpeaks           PROTO((int, char **));
     176int vtransitions     PROTO((int, char **));
    168177int vpop             PROTO((int, char **));
    169178int vsmooth          PROTO((int, char **));
    … …  
    187196static Command cmds[] = { 
    188197  {1, "accum",        accum,            "accumulate vector values in another vector"},
     198  {1, "antialias",    antialias,        "set anti-alias sigma value for display"},
    189199  {1, "applyfit",     applyfit1d,       "apply 1-d fit to new vector"},
    190200  {1, "applyfit1d",   applyfit1d,       "apply 1-d fit to new vector"},
    … …  
    255265  {1, "join",         join,             "find the join of two ID vectors"},
    256266  {1, "jpeg",         jpeg,             "convert display image to JPEG"},
     267  {1, "kapamemory",   kapamemory,       "manage kapa memory dump options"},
    257268  {1, "kern",         kern,             "convolve with 3x3 kernel"},
    258269  {1, "keyword",      keyword,          "extract a FITS keyword from image header"},
    … …  
    265276  {1, "imcreate",     mcreate,          "create an image"},
    266277  {1, "medacc",       medacc,           "accumulate vector values in another vector"},
     278  {1, "mgaussdev",    mgaussdev,        "generate a gaussian deviate image"},
    267279  {1, "mget",         mget,             "extract a vector from an image"},
    268280  {1, "mget3d",       mget3d,           "extract a vector from a 3D image"},
    … …  
    273285  {1, "medimage",     medimage_command, "median image manipulation"},
    274286  {1, "matrix",       matrix,           "matrix math operations"},
     287  {1, "match1d",      match1d,          "match 2 vectors and return matched indexes"},
    275288  {1, "match2d",      match2d,          "match 2 pairs of X,Y vectors and return matched indexes"},
    276289  {1, "mkrgb",        mkrgb,            "convert 3 images to rgb jpeg (use Kapa for better control)"},
    … …  
    347360  {1, "vgauss",       vgauss,           "fit a Gaussian to a vector"},
    348361  {1, "vlorentz",     vlorentz,         "fit a Lorentzian to a vector"},
     362  {1, "vsigmoid",     vsigmoid,         "fit a Sigmoid to a vector"},
    349363  {1, "vellipse",     vellipse,         "fit a Ellipse to a vector pair"},
    350364  {1, "vgrid",        vgrid,            "generate an image from a triplet of vectors"},
    … …  
    354368  {1, "vload",        vload,            "load vectors as overlay on image display"},
    355369  {1, "vmaxwell",     vmaxwell,         "fit a Maxwellian to a vector"},
    356   {1, "vpeaks",       vpeaks,           "fine coord and flux of peaks in vector"},
     370  {1, "vpeaks",       vpeaks,           "find coord and flux of peaks in vector"},
     371  {1, "vtransitions", vtransitions,     "find points in vector that cross the transition value"},
    357372  {1, "vpop",         vpop,             "remove first element of a vector"},
    358373  {1, "vroll",        vroll,            "roll vector elements by 1 entry"},
    … …  
    360375  {1, "vsmooth",      vsmooth,          "Gaussian smooth of a vector"},
    361376  {1, "vstats",       vstats,           "statistics on a vector"},
     377  {1, "vwtmean",      vwtmean,          "weighted mean of a vector"},
     378  {1, "virls",        virls,            "IRLS mean of a vector"},
    362379  {1, "vzload",       vzload,           "load vectors as overlay on image display (scaled points)"},
    363380  {1, "vsh",          vsh,              "Vector Spherical Harmonics"},

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/interpolate_presort.c

    r40007 r41170  
    77  double x0, dx, dy, y0;
    88  Vector *xout, *yout, *xin, *yin;
     9
     10  int N;
     11  int FillEnds = FALSE;
     12  if ((N = get_argument (argc, argv, "-fill-ends"))) {
     13    FillEnds = TRUE;
     14    remove_argument (N, &argc, argv);
     15  }
    916
    1017  /** check basic syntax **/
    … …  
    3643  y0 = yin[0].elements.Flt[0];
    3744 
    38   /* in and out vectors are sorted */
    39   j = 0;
     45  // fill in the start with NANs
     46  for (i = 0; (i < xout[0].Nelements) && (xout[0].elements.Flt[i] < xin[0].elements.Flt[0]); i++) {
     47    yout[0].elements.Flt[i] = FillEnds ? yin[0].elements.Flt[0] : NAN;
     48  }
    4049
    4150  // every output pixel should get a value unless
    4251  // we are below the lowest in or above the highest in
    43   for (i = 0, j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) {
     52  for (j = 0; (i < xout[0].Nelements) && (j < xin[0].Nelements - 1); ) {
    4453
    4554    yout[0].elements.Flt[i] = NAN;
    … …  
    6675
    6776  // fill in the rest with NANs
     77  int NinLast = xin[0].Nelements - 1;
    6878  while (i < yout[0].Nelements) {
    69     yout[0].elements.Flt[i] = NAN;
     79    yout[0].elements.Flt[i] = FillEnds ? yin[0].elements.Flt[NinLast] : NAN;
    7080    i++;
    7181  }
    7282
    7383  return (TRUE);
    74    
    7584}

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/medimage_commands.c

    r40654 r41170  
    1515int medimage_add (int argc, char **argv) {
    1616
     17  int N;
    1718  Buffer *image;
     19  Buffer *var = NULL;
     20
     21  if ((N = get_argument (argc, argv, "-variance"))) {
     22    remove_argument (N, &argc, argv);
     23    if ((var = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     24    remove_argument (N, &argc, argv);
     25  }
    1826
    1927  if (argc != 3) {
    20     gprint (GP_ERR, "USAGE: medimage add (name) (image)\n");
     28    gprint (GP_ERR, "USAGE: medimage add (name) (image) [-variance variance]\n");
    2129    gprint (GP_ERR, "       add the given image to the set of images to be medianed\n");
     30    gprint (GP_ERR, "       optionally supply variance image (for weighted calculations)\n");
    2231    return FALSE;
    2332  }
    2433
    2534  if ((image = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     35
     36  if (var) {
     37    if ((var->matrix.Naxis[0] != image->matrix.Naxis[0]) ||
     38        (var->matrix.Naxis[1] != image->matrix.Naxis[1])) {
     39      gprint (GP_ERR, "variance buffer does not match image buffer dimensions\n");
     40      return FALSE;
     41    }
     42  }
    2643
    2744  MedImageType *median = FindMedImage (argv[1]);
    … …  
    4663  // new image should match existing medimage dimensions
    4764
     65  // AddMedImage (median, image, var);
    4866  int Ninput = median->Ninput;
    4967  median->Ninput ++;
    50   REALLOCATE (median->buffers, float *, median->Ninput);
     68  REALLOCATE (median->flx, float *, median->Ninput);
     69  REALLOCATE (median->var, float *, median->Ninput);
    5170
    52   ALLOCATE (median->buffers[Ninput], float, median->Nx*median->Ny);
    53   memcpy (median->buffers[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
     71  ALLOCATE (median->flx[Ninput], float, median->Nx*median->Ny);
     72  memcpy (median->flx[Ninput], image->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
    5473
    55   return TRUE;
    56 }
    57 
    58 int medimage_calc (int argc, char **argv) {
    59 
    60   int ix, iy, n, N;
    61   Buffer *output;
    62 
    63   int CALC_MEAN = FALSE;
    64   if ((N = get_argument (argc, argv, "-mean"))) {
    65     CALC_MEAN = TRUE;
    66     remove_argument (N, &argc, argv);
    67   }
    68 
    69   if (argc != 3) {
    70     gprint (GP_ERR, "USAGE: medimage calc (name) (output) [-mean]\n");
    71     gprint (GP_ERR, "       calculate the median image for the median image set\n");
    72     return FALSE;
    73   }
    74 
    75   MedImageType *median = FindMedImage (argv[1]);
    76   if (!median) {
    77     gprint (GP_ERR, "median image %s not found\n", argv[1]);
    78     return FALSE;
    79   }
    80 
    81   if ((output = SelectBuffer (argv[2], ANYBUFFER, TRUE)) == NULL) return (FALSE);
    82 
    83   int Ninput = median->Ninput;
    84   int Nx = median->Nx;
    85   int Ny = median->Ny;
    86 
    87   float *value = NULL;
    88   ALLOCATE (value, float, Ninput);
    89 
    90   gfits_free_matrix (&output->matrix);
    91   gfits_free_header (&output->header);
    92   if (!CreateBuffer (output, Nx, Ny, -32, 0.0, 1.0)) return FALSE;
    93 
    94   float *outvalue = (float *) output->matrix.buffer;
    95 
    96   for (iy = 0; iy < Ny; iy++) {
    97     for (ix = 0; ix < Nx; ix++) {
    98 
    99       int N = 0;
    100       int Npix = ix + Nx*iy;
    101       for (n = 0; n < Ninput; n++) {
    102         float v = median->buffers[n][Npix];
    103         if (!isfinite(v)) continue;
    104         value[N] = v;
    105         N++;
    106       }
    107       if (N == 0) continue;
    108 
    109       if (CALC_MEAN) {
    110         float sum = 0.0;
    111         for (n = 0; n < N; n++) {
    112           sum += value[n];
    113         }
    114         outvalue[Npix] = sum / (float) N;
    115       } else {
    116         fsort (value, N);
    117         outvalue[Npix] = value[(int)(0.5*N)];
    118       }
    119     }
     74  median->var[Ninput] = NULL;
     75  if (var) {
     76    ALLOCATE (median->var[Ninput], float, median->Nx*median->Ny);
     77    memcpy (median->var[Ninput], var->matrix.buffer, sizeof(float)*median->Nx*median->Ny);
    12078  }
    12179
    … …  
    12482
    12583/*
    126 int medimage_save (int argc, char **argv) {
     84   int medimage_save (int argc, char **argv) {
    12785
    128   int N;
     86   int N;
    12987
    130   int APPEND = FALSE;
    131   if ((N = get_argument (argc, argv, "-append"))) {
    132     APPEND = TRUE;
    133     remove_argument (N, &argc, argv);
    134   }
     88   int APPEND = FALSE;
     89   if ((N = get_argument (argc, argv, "-append"))) {
     90   APPEND = TRUE;
     91   remove_argument (N, &argc, argv);
     92   }
    13593
    136   if (argc != 3) {
    137     gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n");
    138     return FALSE;
    139   }
     94   if (argc != 3) {
     95   gprint (GP_ERR, "USAGE: medimage save (name) (filename) [-append]\n");
     96   return FALSE;
     97   }
    14098
    141   if (!SaveMedImage(argv[2], argv[1], APPEND)) {
    142     gprint (GP_ERR, "failed to save medimage %s\n", argv[1]);
    143     return (FALSE);
    144   }
    145   return TRUE;
    146 }
     99   if (!SaveMedImage(argv[2], argv[1], APPEND)) {
     100   gprint (GP_ERR, "failed to save medimage %s\n", argv[1]);
     101   return (FALSE);
     102   }
     103   return TRUE;
     104   }
    147105
    148 int medimage_load (int argc, char **argv) {
     106   int medimage_load (int argc, char **argv) {
    149107
    150   if (argc != 3) {
    151     gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n");
    152     return FALSE;
    153   }
     108   if (argc != 3) {
     109   gprint (GP_ERR, "USAGE: medimage load (name) (filename)\n");
     110   return FALSE;
     111   }
    154112
    155   if (!LoadMedImage(argv[2], argv[1])) {
    156     gprint (GP_ERR, "failed to load medimage %s\n", argv[1]);
    157     return (FALSE);
    158   }
    159   return TRUE;
    160 }
     113   if (!LoadMedImage(argv[2], argv[1])) {
     114   gprint (GP_ERR, "failed to load medimage %s\n", argv[1]);
     115   return (FALSE);
     116   }
     117   return TRUE;
     118   }
    161119*/
    162120

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/peak.c

    r20936 r41170  
    3939      if (*X < start) continue;
    4040      if (*X > end) continue;
     41      if (!isfinite(*Y)) continue;
    4142      if (*Y < ymax) continue;
    4243      xmax = *X;
    … …  
    6364      if (*X < start) continue;
    6465      if (*X > end) continue;
     66      if (!isfinite(*Y)) continue;
    6567      if (*Y < ymax) continue;
    6668      xmax = *X;

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/roll.c

    r7917 r41170  
    33int roll (int argc, char **argv) {
    44 
    5   int Nbytes, Nextra;
    6   int dX, dx, dy, nx, ny;
    75  Buffer *buf;
     6
     7  // this function is probably not needed
     8  gprint (GP_ERR, "ERROR: use 'shift' instead of 'roll'\n");
     9  return (FALSE);
    810
    911  if (argc != 4) {
    … …  
    1416  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    1517
    16   dx = atof (argv[2]);
    17   dy = atof (argv[3]);
    18   if (dy != 0) {
    19     gprint (GP_ERR, "only x rolls implemented for now\n");
     18  int dx = atoi (argv[2]);
     19  int dy = atoi (argv[3]);
     20
     21  /* if (dx,dy < 0), we are moving the start position backwards,
     22     if (dx,dy > 0), we are moving the start position forward */
     23 
     24  int Nx = buf[0].matrix.Naxis[0];
     25  int Ny = buf[0].matrix.Naxis[1];
     26
     27  if (dy == 0) {
     28    // moves in just dx can be done row-by-row
     29
     30    int dX = abs(dx);
     31    int Nbytes = (Nx * Ny - dX) * sizeof (float);
     32    int Nextra = dX * sizeof (float);
     33
     34    if (dx < 0) {
     35      memmove (buf[0].matrix.buffer, &buf[0].matrix.buffer[dX*sizeof(float)], Nbytes);
     36      memset (&buf[0].matrix.buffer[Nbytes], 0, Nextra);
     37    } else {
     38      memmove (&buf[0].matrix.buffer[dX*sizeof(float)], buf[0].matrix.buffer, Nbytes);
     39      memset (buf[0].matrix.buffer, 0, Nextra);
     40    }
     41
     42    return TRUE;
    2043  }
    2144
    22   /* if (dx < 0), we are moving the start position back by dx pixels,
    23      if (dx > 0), we are moving the start position forward by dx pixels */
    24  
    25   nx = buf[0].matrix.Naxis[0];
    26   ny = buf[0].matrix.Naxis[1];
     45  if (dx == 0) {
     46    // moves in just dy can be done row-by-row
    2747
    28   dX = abs(dx);
    29   Nbytes = nx * ny * sizeof (float);
    30   Nextra = (nx * ny + dX) * sizeof (float);
     48    ALLOCATE_PTR (output, float, Nx*Ny);
    3149
    32   REALLOCATE (buf[0].matrix.buffer, char, Nextra);
     50    float *input = (float *) buf[0].matrix.buffer;
    3351
    34   if (dx < 0) {
    35     memmove (buf[0].matrix.buffer, &buf[0].matrix.buffer[dX*sizeof(float)], Nbytes);
    36   } else {
    37     memmove (&buf[0].matrix.buffer[dX*sizeof(float)], buf[0].matrix.buffer, Nbytes);
     52    for (int iy = 0; iy < Ny; iy++) {
     53      if (iy + dy < 0) continue;
     54      if (iy + dy >= Ny) continue;
     55      memcpy (&output[iy*Nx], &input[(iy + dy)*Nx], Nx);
     56    }
     57    if (dy > 0) {
     58      memset (output, 0, dy*Nx*sizeof(float));
     59    } else {
     60      int Nbytes = Nx * (Ny - abs(dy)) * sizeof (float);
     61      int Nextra = abs(dy) * Nx * sizeof (float);
     62      memset (&output[Nbytes], 0, Nextra);
     63    }
     64    free (buf[0].matrix.buffer);
     65    buf[0].matrix.buffer = (char *) output;
     66    return TRUE;
    3867  }
    3968
    40   REALLOCATE (buf[0].matrix.buffer, char, Nbytes);
     69  ALLOCATE_PTR (output, float, Nx*Ny);
     70  float *input = (float *) buf[0].matrix.buffer;
    4171
    42   return (TRUE);
    43 
     72  for (int iy = 0; iy < Ny; iy++) {
     73    if (iy + dy < 0) continue;
     74    if (iy + dy >= Ny) continue;
     75    for (int ix = 0; ix < Nx; ix++) {
     76      if (ix + dx < 0) continue;
     77      if (ix + dx >= Nx) continue;
     78      output[ix + iy*Nx] = input[(ix + dx) + (iy + dy)*Nx];
     79    }
     80  }
     81  free (buf[0].matrix.buffer);
     82  buf[0].matrix.buffer = (char *) output;
     83  return TRUE;
    4484}
    4585
     86

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/rotate.c

    r34088 r41170  
    33int rotate (int argc, char **argv) {
    44 
    5   int i, j, NX, NY, X, Y, Lx, Ly, N;
     5  int i, j, NX, NY;
    66  float *in_buff, *out_buff, *c;
    7   double angle, CosAngle, SinAngle, Xo, Yo, dX, dY, fx, fy, x, y, X1, Y1;
     7  double Xo, Yo, x, y, X1, Y1;
    88  double pc11, pc12, pc21, pc22, PC11, PC12, PC21, PC22;
    99  Buffer *buf;
    1010
    11   Xo = 0;
    12   Yo = 0;
    13   if ((N = get_argument (argc, argv, "-center"))) {
    14     remove_argument (N, &argc, argv);
    15     Xo  = atof(argv[N]);
    16     remove_argument (N, &argc, argv);
    17     Yo  = atof(argv[N]);
    18     remove_argument (N, &argc, argv);
    19   }
     11//  Xo = 0;
     12//  Yo = 0;
     13//  if ((N = get_argument (argc, argv, "-center"))) {
     14//    remove_argument (N, &argc, argv);
     15//    Xo  = atof(argv[N]);
     16//    remove_argument (N, &argc, argv);
     17//    Yo  = atof(argv[N]);
     18//    remove_argument (N, &argc, argv);
     19//  }
    2020
    2121  if (argc != 3) {
    … …  
    221221  }
    222222
    223   angle = atof (argv[2]);
    224   CosAngle = cos (angle*RAD_DEG);
    225   SinAngle = sin (angle*RAD_DEG);
     223  double angle = atof (argv[2]);
     224  double CosAngle = cos (angle*RAD_DEG);
     225  double SinAngle = sin (angle*RAD_DEG);
    226226 
    227227  gprint (GP_ERR, "rotating: %f %f %f\n", angle, CosAngle, SinAngle);
    228228
    229   Lx = NX*fabs(CosAngle) + NY*fabs(SinAngle);
    230   Ly = NX*fabs(SinAngle) + NY*fabs(CosAngle);
    231   dX = MAX(0,NY*SinAngle);
    232   dY = MAX(0,-NX*SinAngle);
    233   /*
    234   gprint (GP_ERR, "%f %f  -->  ", Xo, Yo);
    235   X1 = Xo*CosAngle - Yo*SinAngle + dX;
    236   Y1 =  Xo*SinAngle + Yo*CosAngle + dY;
    237   gprint (GP_ERR, "%f %f\n", X1, Y1);
    238   */
     229  // we are rotating about the center pixel, (NX/2, NY/2),
     230  // but we are then putting the result in a new image
     231  // of size (Lx,Ly).
     232
     233  // (x,y) = (i - Nx/2),(j - Ny/2)
     234  // (x',y') = R(theta) (x,y)
     235  // (I,J) = (x' + Lx/2),(y' + Ly/2)
     236
     237  int Lx = NX*fabs(CosAngle) + NY*fabs(SinAngle);
     238  int Ly = NX*fabs(SinAngle) + NY*fabs(CosAngle);
    239239
    240240  /* fix reference pixel */
    241241  gfits_scan (&buf[0].header, "CRPIX1", "%lf", 1, &Xo);
    242242  gfits_scan (&buf[0].header, "CRPIX2", "%lf", 1, &Yo);
    243   /*
    244   X1 = (Xo - dX)*CosAngle - (Yo - dY)*SinAngle;
    245   Y1 = (dX - Xo)*SinAngle + (Yo - dY)*CosAngle;
    246   */
    247   X1 = Xo*CosAngle - Yo*SinAngle + dX;
    248   Y1 =  Xo*SinAngle + Yo*CosAngle + dY;
     243  X1 = (Xo - NX/2)*CosAngle + (Yo - NY/2)*SinAngle + Lx/2;
     244  Y1 = (NX/2 - Xo)*SinAngle + (Yo - NY/2)*CosAngle + Ly/2;
    249245  gfits_modify (&buf[0].header, "CRPIX1", "%lf", 1, X1);
    250246  gfits_modify (&buf[0].header, "CRPIX2", "%lf", 1, Y1);
    … …  
    255251  gfits_scan (&buf[0].header, "PC002001", "%lf", 1, &pc21);
    256252  gfits_scan (&buf[0].header, "PC002002", "%lf", 1, &pc22);
    257   PC11 = pc11*CosAngle - pc21*SinAngle;
    258   PC12 = pc12*CosAngle - pc22*SinAngle;
    259   PC21 = pc21*CosAngle + pc11*SinAngle;
    260   PC22 = pc22*CosAngle + pc12*SinAngle;
     253  PC11 = pc11*CosAngle + pc21*SinAngle;
     254  PC12 = pc12*CosAngle + pc22*SinAngle;
     255  PC21 = pc21*CosAngle - pc11*SinAngle;
     256  PC22 = pc22*CosAngle - pc12*SinAngle;
    261257  gfits_modify (&buf[0].header, "PC001001", "%le", 1, PC11);
    262258  gfits_modify (&buf[0].header, "PC001002", "%le", 1, PC12);
    … …  
    270266  gfits_create_matrix (&buf[0].header, &buf[0].matrix);
    271267  gfits_print_alt (&buf[0].header, "HISTORY", "%S", 1, "WARNING: rotated image!");
     268
    272269  out_buff = (float *)buf[0].matrix.buffer;
    273270  for (j = 0; j < Ly; j++) {
    274271    for (i = 0; i < Lx; i++, out_buff++) {
    275272
    276       x = (i - dX)*CosAngle + (j - dY)*SinAngle;
    277       y = (dX - i)*SinAngle + (j - dY)*CosAngle;
    278       X = (int) x;
    279       Y = (int) y;
    280 
    281       if (X < 0) continue;
    282       if (X >= NX - 1) continue;
    283       if (Y < 0) continue;
    284       if (Y >= NY - 1) continue;
    285 
    286       c = &in_buff[X + NX*Y];
    287       fx = x - X;
    288       fy = y - Y;
     273      float xo = (i - Lx/2);
     274      float yo = (j - Ly/2);
     275
     276      x =  xo*CosAngle + yo*SinAngle;
     277      y = -xo*SinAngle + yo*CosAngle;
     278
     279      int I = x + NX/2;
     280      int J = y + NY/2;
     281
     282      if (I < 0) continue;
     283      if (I >= NX - 1) continue;
     284      if (J < 0) continue;
     285      if (J >= NY - 1) continue;
     286
     287      c = &in_buff[I + NX*J];
     288
     289      int X = (int) x;
     290      int Y = (int) y;
     291      double fx = x - X;
     292      double fy = y - Y;
     293
    289294      *out_buff = (c[0]*(1-fx) + c[1]*fx)*(1-fy) + (c[NX+1]*fx + c[NX]*(1-fx))*fy;
     295//    *out_buff = c[0];
    290296    }
    291297  }

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/spline.c

    r39457 r41170  
    77int spline_load (int argc, char **argv);
    88int spline_save (int argc, char **argv);
     9int spline_print (int argc, char **argv);
    910int spline_delete (int argc, char **argv);
    1011int spline_rename (int argc, char **argv);
    … …  
    1819  {1, "load",       spline_load,       "write a spline to a FITS file"},
    1920  {1, "save",       spline_save,       "read a spline from a FITS file"},
     21  {1, "print",      spline_print,      "print a spline to stdout"},
    2022  {1, "delete",     spline_delete,     "delete a spline"},
    2123  {1, "rename",     spline_rename,     "rename a spline"},
    … …  
    3537  gprint (GP_ERR, "    spline load   (spline) (filename)            : load a spline from a FITS file\n");
    3638  gprint (GP_ERR, "    spline save   (spline) (filename) [-append]  : save a spline in FITS format\n");
     39  gprint (GP_ERR, "    spline print  (spline)                       : print a spline to stdout\n");
    3740
    3841  return FALSE;

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/spline_commands.c

    r39457 r41170  
    2828  if (xvec->Nelements != yvec->Nelements) {
    2929    gprint (GP_ERR, "x and y vectors must be the same length\n");
     30    return FALSE;
     31  }
     32
     33  if (xvec->Nelements < 3) {
     34    gprint (GP_ERR, "cannot make a spline with fewer than 3 knots\n");
    3035    return FALSE;
    3136  }
    … …  
    9398}
    9499
     100int spline_print (int argc, char **argv) {
     101
     102  int i;
     103
     104  if (argc != 2) {
     105    gprint (GP_ERR, "USAGE: spline print (name)\n");
     106    return FALSE;
     107  }
     108
     109  Spline *myspline = FindSpline (argv[1]);
     110  if (!myspline) {
     111    gprint (GP_ERR, "spline %s not found\n", argv[1]);
     112    return (FALSE);
     113  }
     114
     115  for (i = 0; i < myspline->Nknots; i++) {
     116    gprint (GP_LOG, "%lf : %lf : %lf\n", myspline->xk[i], myspline->yk[i], myspline->y2[i]);
     117  }   
     118
     119  return TRUE;
     120}
     121
    95122int spline_load (int argc, char **argv) {
    96123
    … …  
    109136int spline_delete (int argc, char **argv) {
    110137
    111   int status;
    112   Spline *spline;
     138  int N, status;
     139
     140  int QUIET = FALSE;
     141  if ((N = get_argument (argc, argv, "-q"))) {
     142    QUIET = TRUE;
     143    remove_argument (N, &argc, argv);
     144  }
    113145
    114146  if (argc != 2) {
    … …  
    117149  }
    118150
    119   spline = FindSpline (argv[1]);
     151  Spline *spline = FindSpline (argv[1]);
    120152  if (spline == NULL) {
     153    if (QUIET) return TRUE;
    121154    gprint (GP_ERR, "spline %s not found\n", argv[1]);
    122155    return FALSE;
    … …  
    130163int spline_rename (int argc, char **argv) {
    131164
    132   Spline *spline;
    133 
    134165  if (argc != 3) {
    135166    gprint (GP_ERR, "USAGE: spline rename (spline) (newname)\n");
    … …  
    137168  }
    138169
    139   spline = FindSpline (argv[1]);
     170  Spline *spline = FindSpline (argv[1]);
    140171  if (spline == NULL) {
    141172    gprint (GP_ERR, "spline %s not found\n", argv[1]);
    142173    return FALSE;
     174  }
     175
     176  Spline *newname = FindSpline (argv[2]);
     177  if (newname != NULL) {
     178    DeleteSpline (newname);
    143179  }
    144180
    … …  
    147183  return TRUE;
    148184}
    149 
    150 /*
    151 int spline_listspline (int argc, char **argv) {
    152 
    153   Spline *spline;
    154 
    155   if (argc != 2) {
    156     gprint (GP_ERR, "USAGE: spline listspline (spline)\n");
    157     return FALSE;
    158   }
    159 
    160   spline = FindSpline (argv[1]);
    161   if (spline == NULL) {
    162     gprint (GP_ERR, "spline %s not found\n", argv[1]);
    163     return FALSE;
    164   }
    165 
    166   ListPages (spline);
    167   return TRUE;
    168 }
    169 */
    170 

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/test/cut.sh

    r16056 r41170  
    44 memtest1
    55end
     6
     7macro test2
     8
     9  mgaussdev z 50 100 0.0 1.0
     10
     11  cut z x y_sum    x 0 0 z[][0] z[0][]
     12  cut z x y_median x 0 0 z[][0] z[0][] -median
     13  cut z x y_mean   x 0 0 z[][0] z[0][] -mean
     14  cut z x y_inner  x 0 0 z[][0] z[0][] -inner
     15 
     16end
     17
     18macro checkrange
     19  if ($0 != 2)
     20    echo "USAGE: checkrange (Npts)"
     21    break
     22  end
     23
     24  $Npts = $1
     25  if ($Npts % 2)
     26    $Ncenter = int($Npts / 2)
     27    $Nquarter = int(0.25 * $Npts)
     28    $Ns = $Ncenter - $Nquarter
     29    $Ne = $Ncenter + $Nquarter
     30  else
     31    $Ncenter = int($Npts / 2) - 1
     32    $Nquarter = int(0.25 * $Npts)
     33    $Ns = $Ncenter - $Nquarter
     34    $Ne = $Ncenter + $Nquarter + 1
     35  end
     36  if ($Ns < 0)
     37    echo "error: Ns < 0: $Ns $Ncenter $Nquarter"
     38  end
     39  if ($Ne >= $Npts)
     40    echo "error: Ne >= Npts: $Ne $Ncenter $Nquarter"
     41  end
     42  echo "$Npts : $Ncenter : $Nquarter : $Ns $Ne"
     43  for i 0 $Npts
     44    echo $i {($i >= $Ns) && ($i <= $Ne)}
     45  end
     46end
     47
    648
    749# Test if cut works
    … …  
    3779 end
    3880end
     81
     82macro memtest1
     83
     84 $i = 0
     85 mcreate tim 100 100
     86 cut tim xdir imx X 40 4 60 6
     87 # do one to set up memory that should stay used
     88
     89 memory check
     90 for i 0 100
     91   cut tim xdir imx X 40 4 60 6
     92 end
     93 memory check
     94   
     95 for i 0 100
     96   cut tim xdir imx y 40 4 60 6
     97 end
     98 memory check
     99   
     100 for i 0 100
     101   cut -median tim xdir imx X 40 4 60 6
     102 end
     103 memory check
     104   
     105 for i 0 100
     106   cut -median tim xdir imx y 40 4 60 6
     107 end
     108 memory check
     109   
     110 for i 0 100
     111   cut -mean tim xdir imx X 40 4 60 6
     112 end
     113 memory check
     114   
     115 for i 0 100
     116   cut -mean tim xdir imx y 40 4 60 6
     117 end
     118 memory check
     119   
     120end

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/test/delete.sh

    r14176 r41170  
    4949
    5050end
     51
     52# Test delete
     53macro test2
     54
     55 $PASS = 1
     56
     57 for i 0 1000
     58   create v 0 200
     59   delete v
     60 end
     61end

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/test/medimage.sh

    r36679 r41170  
    11
    2 macro go
    3  mcreate a 30 30
    4  for i 0 40
    5   set a$i = zero(a) + $i
    6  end
    7 
    8  for i 0 40
    9   medimage add t1 a$i
    10  end
    11 end
     2macro test.mean
     3
     4 $Nsample = 16
     5 medimage delete -q t1
     6 for i 0 $Nsample
     7  mgaussdev t 50 50 0.0 1.0
     8  medimage add t1 t
     9  imhist -q t x y -range -10 10 -delta 0.1
     10
     11  $C0 = 0
     12  $C1 = 1.5
     13  $C2 = 400
     14  $C3 = 0
     15  vgauss -q x y con yf
     16  # echo $C1
     17 end
     18
     19 medimage calc t1 T -mean
     20
     21 imhist -q T x y -range -10 10 -delta 0.1
     22 lim -n 1 x y; clear; box; plot -x hist x y
     23 $C0 = 0
     24 $C1 = 1.5
     25 $C2 = 400
     26 $C3 = 0
     27 vgauss -q x y con yf
     28 echo "expect {1/sqrt($Nsample)} : $C1"
     29 plot -c red -x line x yf
     30end
     31
     32macro test.median
     33
     34 $Nsample = 16
     35 medimage delete -q t1
     36 for i 0 $Nsample
     37  mgaussdev t 50 50 0.0 1.0
     38  medimage add t1 t
     39  imhist -q t x y -range -10 10 -delta 0.1
     40
     41  $C0 = 0
     42  $C1 = 1.5
     43  $C2 = 400
     44  $C3 = 0
     45  vgauss -q x y con yf
     46 end
     47
     48 # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
     49 # (somewhat higher scatter)
     50 medimage calc t1 T
     51
     52 imhist -q T x y -range -10 10 -delta 0.1
     53 lim -n 1 x y; clear; box; plot -x hist x y
     54 $C0 = 0
     55 $C1 = 1.5
     56 $C2 = 400
     57 $C3 = 0
     58 vgauss -q x y con yf
     59 echo "expect {1/sqrt($Nsample)} : $C1 (actually should be a bit higher)"
     60 plot -c red -x line x yf
     61end
     62 
     63macro test.wtmean
     64
     65 $Nsample = 8
     66 $sig1 = 1.0
     67 $sig2 = 3.0
     68
     69 medimage delete -q t1
     70 for i 0 $Nsample
     71  mgaussdev t 50 50 0.0 $sig1
     72  set v = $sig1^2 + zero(t)     
     73  medimage add t1 t -variance v
     74  imhist -q t x y -range -10 10 -delta 0.1
     75
     76  $C0 = 0
     77  $C1 = 1.5
     78  $C2 = 400
     79  $C3 = 0
     80  vgauss -q x y con yf
     81  # echo $C1
     82 end
     83
     84 for i 0 $Nsample
     85  mgaussdev t 50 50 0.0 $sig2
     86  set v = $sig2^2 + zero(t)     
     87  medimage add t1 t -variance v
     88  imhist -q t x y -range -10 10 -delta 0.1
     89
     90  $C0 = 0
     91  $C1 = 1.5
     92  $C2 = 400
     93  $C3 = 0
     94  vgauss -q x y con yf
     95  # echo $C1
     96 end
     97
     98 # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
     99 # (somewhat higher scatter)
     100 medimage calc t1 T -wtmean
     101
     102 imhist -q T x y -range -10 10 -delta 0.1
     103 lim -n 1 x y; clear; box; plot -x hist x y
     104 $C0 = 0
     105 $C1 = 1.5
     106 $C2 = 400
     107 $C3 = 0
     108 vgauss -q x y con yf
     109 $S1 = $Nsample / $sig1^2 + $Nsample / $sig2^2
     110 echo "expect {1/sqrt($S1)} : $C1"
     111 plot -c red -x line x yf
     112end
     113
     114macro test.irls
     115 medimage delete -q t1
     116 $Nsample = 16
     117 $sig = 1.0
     118 for i 0 $Nsample
     119  mgaussdev t 50 50 0.0 $sig
     120  set v = $sig^2 + zero(t)     
     121
     122  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
     123  set ts = t + bad
     124
     125  medimage add t1 ts -variance v
     126  imhist -q t x y -range -10 10 -delta 0.1
     127
     128  $C0 = 0
     129  $C1 = 1.5
     130  $C2 = 400
     131  $C3 = 0
     132  vgauss -q x y con yf
     133  # echo $C1
     134 end
     135
     136 # get stats for straight mean:
     137 medimage calc t1 Tm -mean
     138
     139 imhist -q Tm x y -range -10 10 -delta 0.1
     140 lim -n 1 x y; clear; box; plot -x hist x y
     141 $C0 = 0
     142 $C1 = 1.5
     143 $C2 = 400
     144 $C3 = 0
     145 vgauss -q x y con yf
     146 echo "sigma from straight stdev: $C1"
     147 # stats Tm
     148 
     149 plot -c red -x line x yf
     150
     151 # get stats for irls
     152 medimage calc t1 Ti -irls
     153
     154 imhist -q Ti x y -range -10 10 -delta 0.1
     155 lim -n 2 x y; clear; box; plot -x hist x y
     156 $C0 = 0
     157 $C1 = 1.5
     158 $C2 = 400
     159 $C3 = 0
     160 vgauss -q x y con yf
     161 echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"
     162 # stats Ti
     163 
     164 plot -c red -x line x yf
     165end
     166
     167
     168###################33
     169
     170
     171macro test.mean.var
     172
     173 $Nsample = 64
     174 $sig = 2.0
     175 medimage delete -q t1
     176 for i 0 $Nsample
     177  mgaussdev t 100 100 0.0 $sig
     178  medimage add t1 t
     179  imhist -q t x y -range -10 10 -delta 0.1
     180
     181  $C0 = 0
     182  $C1 = 1.5
     183  $C2 = 400
     184  $C3 = 0
     185  vgauss -q x y con yf
     186  # echo $C1
     187 end
     188
     189 medimage calc t1 T -mean -variance Tv
     190
     191 imhist -q T x y -range -10 10 -delta 0.1
     192 lim -n 1 x y; clear; box; plot -x hist x y
     193 $C0 = 0
     194 $C1 = 1.5
     195 $C2 = 400
     196 $C3 = 0
     197 vgauss -q x y con yf
     198 plot -c red -x line x yf
     199
     200 imhist Tv xv yv -range -1 4 -delta 0.1
     201 lim -n 2 xv yv; clear; box; plot xv yv -x hist
     202
     203 stat Tv
     204 echo "$C1 vs {sqrt($MEDIAN)} : expect {$sig/sqrt($Nsample)}"
     205end
     206
     207macro test.median.var
     208
     209 $Nsample = 64
     210 $sig = 2.0
     211 medimage delete -q t1
     212 for i 0 $Nsample
     213  mgaussdev t 50 50 0.0 $sig
     214  medimage add t1 t
     215  imhist -q t x y -range -10 10 -delta 0.1
     216
     217  $C0 = 0
     218  $C1 = 1.5
     219  $C2 = 400
     220  $C3 = 0
     221  vgauss -q x y con yf
     222  # echo $C1
     223 end
     224
     225 # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
     226 # (somewhat higher scatter)
     227 medimage calc t1 T -variance Tv
     228
     229 imhist -q T x y -range -10 10 -delta 0.1
     230 lim -n 1 x y; clear; box; plot -x hist x y
     231 $C0 = 0
     232 $C1 = 1.5
     233 $C2 = 400
     234 $C3 = 0
     235 vgauss -q x y con yf
     236 plot -c red -x line x yf
     237
     238 imhist Tv xv yv -range -1 4 -delta 0.1
     239 lim -n 2 xv yv; clear; box; plot xv yv -x hist
     240
     241 stat Tv
     242 echo "$C1 vs {sqrt($MEDIAN)} : expect {$sig/sqrt($Nsample)}"
     243end
     244 
     245macro test.wtmean.var
     246
     247 $Nsample = 32
     248 $sig1 = 1.0
     249 $sig2 = 1.0
     250
     251 medimage delete -q t1
     252 for i 0 $Nsample
     253  mgaussdev t 50 50 0.0 $sig1
     254  set v = $sig1^2 + zero(t)     
     255  medimage add t1 t -variance v
     256  imhist -q t x y -range -10 10 -delta 0.1
     257
     258  $C0 = 0
     259  $C1 = 1.5
     260  $C2 = 400
     261  $C3 = 0
     262  vgauss -q x y con yf
     263  # echo $C1
     264 end
     265
     266 for i 0 $Nsample
     267  mgaussdev t 50 50 0.0 $sig2
     268  set v = $sig2^2 + zero(t)     
     269  medimage add t1 t -variance v
     270  imhist -q t x y -range -10 10 -delta 0.1
     271
     272  $C0 = 0
     273  $C1 = 1.5
     274  $C2 = 400
     275  $C3 = 0
     276  vgauss -q x y con yf
     277  # echo $C1
     278 end
     279
     280 # note that median of a gaussian distributed variable is not distributed with sigma' = sigma / sqrt(N)
     281 # (somewhat higher scatter)
     282 medimage calc t1 T -wtmean -variance Tv
     283
     284 imhist -q T x y -range -10 10 -delta 0.1
     285 lim -n 1 x y; clear; box; plot -x hist x y
     286 $C0 = 0
     287 $C1 = 1.5
     288 $C2 = 400
     289 $C3 = 0
     290 vgauss -q x y con yf
     291 plot -c red -x line x yf
     292
     293 stat -q Tv
     294 $S1 = $Nsample / $sig1^2 + $Nsample / $sig2^2
     295 echo $C1 vs {sqrt($MEDIAN)} : expect {1/sqrt($S1)}
     296end
     297
     298macro test.irls.var
     299
     300 $Nsample = 16
     301 $sig = 1.0
     302
     303 medimage delete -q t1
     304 for i 0 $Nsample
     305  mgaussdev t 50 50 0.0 $sig
     306  set v = $sig^2 + zero(t)     
     307
     308  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
     309  set ts = t + bad
     310
     311  medimage add t1 ts -variance v
     312  imhist -q t x y -range -10 10 -delta 0.1
     313
     314  $C0 = 0
     315  $C1 = 1.5
     316  $C2 = 400
     317  $C3 = 0
     318  vgauss -q x y con yf
     319  # echo $C1
     320 end
     321
     322 # get stats for straight mean:
     323 medimage calc t1 Tm -mean -variance Tv
     324
     325 imhist -q Tm x y -range -10 10 -delta 0.1
     326 lim -n 1 x y; clear; box; plot -x hist x y
     327 $C0 = 0
     328 $C1 = 1.5
     329 $C2 = 400
     330 $C3 = 0
     331 vgauss -q x y con yf
     332 echo "sigma from straight stdev: $C1"
     333 # stats Tm
     334 
     335 plot -c red -x line x yf
     336
     337 # get stats for irls
     338 medimage calc t1 Ti -irls -variance Tv
     339
     340 imhist -q Ti x y -range -10 10 -delta 0.1
     341 lim -n 2 x y; clear; box; plot -x hist x y
     342 $C0 = 0
     343 $C1 = 1.5
     344 $C2 = 400
     345 $C3 = 0
     346 vgauss -q x y con yf
     347 echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"
     348 # stats Ti
     349 
     350 plot -c red -x line x yf
     351
     352 set dTv = sqrt(Tv)
     353 imhist dTv xv yv -range -1 4 -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist
     354
     355 stat -q Tv
     356 echo $C1 vs {sqrt($MEDIAN)} (ideal is {$sig/sqrt($Nsample)})"
     357end
     358
     359macro test.irls.boot.var
     360
     361 $Nsample = 64
     362 $sig = 1.0
     363
     364 medimage delete -q t1
     365 for i 0 $Nsample
     366  mgaussdev t 200 200 0.0 $sig
     367  set v = $sig^2 + zero(t)     
     368
     369  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
     370  set ts = t + bad
     371
     372  mgaussdev noise 200 200 0.0 0.5
     373  set ts = ts + noise
     374
     375  medimage add t1 ts -variance v
     376  imhist -q t x y -range -10 10 -delta 0.1
     377
     378  $C0 = 0
     379  $C1 = 1.5
     380  $C2 = 400
     381  $C3 = 0
     382  vgauss -q x y con yf
     383  # echo $C1
     384 end
     385
     386 # get stats for straight mean:
     387 medimage calc t1 Tm -mean -variance Tv
     388
     389 imhist -q Tm x y -range -10 10 -delta 0.02
     390 lim -n 1 x y; clear; box; plot -x hist x y
     391 peak -q x y
     392 $C0 = $peakpos
     393 $C1 = 1.5*$sig / sqrt($Nsample)
     394 $C2 = $peakval
     395 $C3 = 0
     396 vgauss -q x y con yf
     397 echo "sigma from straight stdev: $C1"
     398 # stats Tm
     399 
     400 plot -c red -x line x yf
     401
     402 # get stats for irls
     403 date
     404 medimage calc t1 Ti -irls -variance Tv -bootstrap-iter 100
     405 date
     406
     407 imhist -q Ti x y -range -10 10 -delta 0.02
     408 lim -n 2 x y; clear; box; plot -x hist x y
     409 peak -q x y
     410 $C0 = $peakpos
     411 $C1 = 1.5*$sig / sqrt($Nsample)
     412 $C2 = $peakval
     413 $C3 = 0
     414 vgauss -q x y con yf
     415 echo "sigma from irls: $C1 (ideal is {$sig/sqrt($Nsample)})"
     416 # stats Ti
     417 
     418 plot -c red -x line x yf
     419
     420 set dTv = sqrt(Tv)
     421 imhist dTv xv yv -range 0 {5*$sig/sqrt($Nsample)} -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist
     422
     423 stat -q Tv
     424 echo $C1 vs {sqrt($MEDIAN)} (ideal is {$sig/sqrt($Nsample)})"
     425end
     426
     427##############################
     428macro test.irls.boot.test
     429
     430 $Nsample = 100
     431 $sig1 = 1.0
     432
     433 medimage delete -q t1
     434 for i 0 $Nsample
     435  mgaussdev t 100 100 0.0 $sig1
     436  set v = $sig1^2 + zero(t)     
     437
     438  medimage add t1 t -variance v
     439 end
     440
     441 # get stats for irls
     442 medimage calc t1 Ti -irls -variance Tv -bootstrap
     443
     444 imhist -q Ti x y -range {-10*$sig1/sqrt($Nsample)} {10*$sig1/sqrt($Nsample)} -delta 0.01
     445 lim -n 2 x y; clear; box; plot -x hist x y
     446 peak -q x y
     447 $C0 = $peakpos
     448 $C1 = 1.5*$sig1/sqrt($Nsample)
     449 $C2 = $peakval
     450 $C3 = 0
     451 vgauss x y con yf
     452 echo "sigma from irls: $C1 (ideal is {$sig1/sqrt($Nsample)})"
     453 # stats Ti
     454 
     455 plot -c red -x line x yf
     456
     457 set dTv = sqrt(Tv)
     458 imhist dTv xv yv -range 0 {5*$sig1/sqrt($Nsample)} -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist
     459
     460 stat -q irls_npt
     461 $Npix = $MEAN
     462
     463 stat -q Tv
     464 echo "sigma of irls average: $C1, sqrt(mean) of irls variance: {sqrt($MEAN)}, (ideal is {$sig1/sqrt($Npix)})"
     465end
     466
     467macro test.irls.range.var
     468 medimage delete -q t1
     469 for i 0 8
     470  mgaussdev t 50 50 0.0 1.0
     471  set v = 1.0 + zero(t)     
     472
     473  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
     474  set ts = t + bad
     475
     476  medimage add t1 ts -variance v
     477  imhist -q t x y -range -10 10 -delta 0.1
     478
     479  $C0 = 0
     480  $C1 = 1.5
     481  $C2 = 400
     482  $C3 = 0
     483  vgauss -q x y con yf
     484  echo $C1
     485 end
     486
     487 for i 0 8
     488  mgaussdev t 50 50 0.0 3.0
     489  set v = 3.0 + zero(t)     
     490
     491  set bad = (rnd(t) < 0.05) ? 10*rnd(t) + 5 : zero(t)
     492  set ts = t + bad
     493
     494  medimage add t1 ts -variance v
     495  imhist -q t x y -range -10 10 -delta 0.1
     496
     497  $C0 = 0
     498  $C1 = 1.5
     499  $C2 = 400
     500  $C3 = 0
     501  vgauss -q x y con yf
     502  echo $C1
     503 end
     504
     505 # get stats for straight mean:
     506 medimage calc t1 Tm -mean -variance Tv
     507
     508 imhist -q Tm x y -range -10 10 -delta 0.1
     509 lim -n 1 x y; clear; box; plot -x hist x y
     510 $C0 = 0
     511 $C1 = 1.5
     512 $C2 = 400
     513 $C3 = 0
     514 vgauss -q x y con yf
     515 echo $C1
     516 stats Tm
     517 
     518 plot -c red -x line x yf
     519
     520 # get stats for irls
     521 medimage calc t1 Ti -irls -variance Tv
     522
     523 imhist -q Ti x y -range -10 10 -delta 0.1
     524 lim -n 2 x y; clear; box; plot -x hist x y
     525 $C0 = 0
     526 $C1 = 1.5
     527 $C2 = 400
     528 $C3 = 0
     529 vgauss -q x y con yf
     530 echo $C1
     531 stats Ti
     532 
     533 plot -c red -x line x yf
     534
     535 stat -q Tv
     536 echo $C1 vs {sqrt($MEDIAN)}
     537
     538 set dTv = sqrt(Tv)
     539 imhist dTv xv yv -range -1 4 -delta 0.02; lim -n 3 xv yv; clear; box; plot xv yv -x hist
     540end

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/threshold.c

    r40623 r41170  
    2525    remove_argument (N, &argc, argv);
    2626    BinMax = atoi (argv[N]);
     27    remove_argument (N, &argc, argv);
     28  }
     29
     30  float ValMin = NAN;
     31  float ValMax = NAN;
     32  if ((N = get_argument (argc, argv, "-vrange"))) {
     33    remove_argument (N, &argc, argv);
     34    ValMin = atof (argv[N]);
     35    remove_argument (N, &argc, argv);
     36    ValMax = atof (argv[N]);
    2737    remove_argument (N, &argc, argv);
    2838  }
    … …  
    6474  }
    6575
     76  // ValMin, ValMax not compatible with BinMin, BinMax
     77  if (isfinite(ValMin)) {
     78    int binTest = ohana_bisection_double (vecx[0].elements.Flt, vecx[0].Nelements, ValMin);
     79    BinMin = (vecx[0].elements.Flt[binTest+1] == ValMin) ? binTest : binTest + 1;
     80    // XXX dangerous: check that this does not run off the end, etc
     81  }
     82  if (isfinite(ValMax)) {
     83    int binTest = ohana_bisection_double (vecx[0].elements.Flt, vecx[0].Nelements, ValMax);
     84    BinMax = (vecx[0].elements.Flt[binTest+1] == ValMax) ? binTest : binTest + 1;
     85    // XXX dangerous: check that this does not run off the end, etc
     86  }
     87
    6688  // use bisection to find the value
    6789
    … …  
    7193
    7294  if (!REVERSE) {
     95    // this algorithm assumes vecy[BinMax] > threshold, vecy[BinMin] < threshold
     96    if (vecy[0].elements.Flt[BinMin] > value) {
     97      gprint (GP_ERR, "ERROR: all values above threshold\n");
     98      return FALSE;
     99    }
     100    if (vecy[0].elements.Flt[BinMax] < value) {
     101      gprint (GP_ERR, "ERROR: all values below threshold\n");
     102      return FALSE;
     103    }
    73104    while (Nhi - Nlo > 10) {
    74105      N = 0.5*(Nlo + Nhi);
    … …  
    89120    }
    90121  } else {
     122    // this algorithm assumes vecy[BinMin] > threshold, vecy[BinMax] < threshold
     123    if (vecy[0].elements.Flt[BinMin] < value) {
     124      gprint (GP_ERR, "ERROR: all values below threshold\n");
     125      return FALSE;
     126    }
     127    if (vecy[0].elements.Flt[BinMax] > value) {
     128      gprint (GP_ERR, "ERROR: all values above threshold\n");
     129      return FALSE;
     130    }
    91131    while (Nhi - Nlo > 10) {
    92132      N = 0.5*(Nlo + Nhi);

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/vgauss.c

    r39457 r41170  
    33/* local private functions */
    44opihi_flt fgaussOD (opihi_flt, opihi_flt *, int, opihi_flt *);
     5int vgauss_apply (int argc, char **argv);
    56
    67# define GET_VAR(V,A) \
    … …  
    2021  Vector *xvec, *yvec, *svec, *ovec;
    2122  char *c, name[16];
     23
     24  if ((N = get_argument (argc, argv, "-apply"))) {
     25    remove_argument (N, &argc, argv);
     26    int status = vgauss_apply (argc, argv);
     27    return status;
     28  }
    2229
    2330  Quiet = FALSE;
    … …  
    104111  } 
    105112  if (!Quiet) gprint (GP_ERR, "%d iterations\n", i);
     113  set_variable ("Chisq", chisq);
     114  set_variable ("ChisqNu", chisq / (float) (Npts - Npar));
    106115
    107116  for (i = 0; i < Npts; i++) {
    … …  
    125134}
    126135
     136int vgauss_apply (int argc, char **argv) {
     137
     138  opihi_flt par[4];
     139  Vector *xvec, *ovec;
     140  char *c;
     141
     142  // NOTE: -apply has already been stripped by vgauss call
     143  if (argc != 3) {
     144    gprint (GP_ERR, "USAGE: vgauss -apply <x> <out>\n");
     145    return (FALSE);
     146  }
     147 
     148  if ((xvec = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
     149  if ((ovec = SelectVector (argv[2], ANYVECTOR, TRUE)) == NULL) return (FALSE);
     150
     151  CastVector (xvec, OPIHI_FLT);
     152
     153  GET_VAR (par[0], "C0");
     154  GET_VAR (par[1], "C1");
     155  GET_VAR (par[2], "C2");
     156  GET_VAR (par[3], "C3");
     157  int Npar = 4;
     158
     159  int Npts = xvec[0].Nelements;
     160  ResetVector (ovec, OPIHI_FLT, Npts);
     161
     162  for (int i = 0; i < Npts; i++) {
     163    ovec[0].elements.Flt[i] = fgaussOD (xvec[0].elements.Flt[i], par, Npar, NULL);
     164  }
     165  ovec[0].Nelements = Npts;
     166
     167  return TRUE;
     168}
     169
    127170/* pars: x_o, sigma, I, back */
    128171opihi_flt fgaussOD (opihi_flt x, opihi_flt *par, int Npar, opihi_flt *dpar) {
    … …  
    135178  f = par[2]*r + par[3];
    136179
    137   dpar[0] = par[2]*r*z/par[1];
    138   dpar[1] = par[2]*r*z*z/par[1];
    139   dpar[2] = r;
    140   dpar[3] = 1;
     180  if (dpar) {
     181    dpar[0] = par[2]*r*z/par[1];
     182    dpar[1] = par[2]*r*z*z/par[1];
     183    dpar[2] = r;
     184    dpar[3] = 1;
     185  }
    141186 
    142187  return (f);

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/cmd.data/write_vectors.c

    r40574 r41170  
    3232  }
    3333
    34   /* option generate a FITS output table */
     34  // option generate a FITS output table (FITS holds the filename)
     35  // in FITS output context, -header is interpretted as a buffer containing
     36  // header keywords to supplement the FITS table header
    3537  char *FITS = NULL;
     38  Header *fitsheader = NULL;
     39  Buffer *headbuffer = NULL;
    3640  if ((N = get_argument (argc, argv, "-fits"))) {
    3741    remove_argument (N, &argc, argv);
    3842    FITS = strcreate (argv[N]);
    3943    remove_argument (N, &argc, argv);
    40   }
    41 
     44
     45    if ((N = get_argument (argc, argv, "-header"))) {
     46      remove_argument (N, &argc, argv);
     47      if ((headbuffer = SelectBuffer (argv[N], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     48      fitsheader = &headbuffer->header;
     49      remove_argument (N, &argc, argv);
     50    }
     51  }
     52 
    4253  /* option generate a FITS output table */
    4354  int CSV = FALSE;
    … …  
    124135
    125136  if (FITS) {
    126     int status = WriteVectorTableFITS (argv[1], FITS, vec, Nvec, append, compress, format, Ntile);
     137    int status = WriteVectorTableFITS (argv[1], FITS, fitsheader, vec, Nvec, append, compress, format, Ntile);
    127138    free (vec);
    128139    return status;
    … …  
    260271    gprint (GP_ERR, "  -append : write to the end of the existing file\n");
    261272    gprint (GP_ERR, "  -fits NAME : write a fits table (extention name is NAME, column names match vector names)\n");
     273    gprint (GP_ERR, "     in FITS output context, -header takes an additional argument which is interpretted\n");
     274    gprint (GP_ERR, "     as a buffer containing header keywords to supplement the FITS table header\n");
    262275    gprint (GP_ERR, "  -csv : write a comma-separated values file (eg, to read in excel)\n");
    263276    gprint (GP_ERR, "  -f \"format\" : provide formatting codes for output:\n");

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/dvo/Makefile

    r40523 r41170  
    6767$(SRC)/gstar.$(ARCH).o                  \
    6868$(SRC)/hosts.$(ARCH).o          \
     69$(SRC)/cimages.$(ARCH).o                \
    6970$(SRC)/images.$(ARCH).o         \
    7071$(SRC)/imbox.$(ARCH).o                  \

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/dvo/avextract.c

    r40574 r41170  
    272272  // write vectors to a table (this is used by parallel dvo operations, but can be used elsewhere)
    273273  if (RESULT_FILE && !SKIP_RESULTS) {
    274     int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, Nreturn, FALSE, FALSE, NULL, 0);
     274    int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", NULL, vec, Nreturn, FALSE, FALSE, NULL, 0);
    275275    if (!status) {
    276276      goto escape;

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/dvo/avmatch.c

    r40523 r41170  
    114114
    115115      CoordsFile = abspath("coords.fits", 1024);
    116       int status = WriteVectorTableFITS (CoordsFile, "COORDS", vec, 2, FALSE, FALSE, NULL, 0);
     116      int status = WriteVectorTableFITS (CoordsFile, "COORDS", NULL, vec, 2, FALSE, FALSE, NULL, 0);
    117117      if (!status) goto escape;
    118118    }
    … …  
    301301      vec[i][0].Nelements = Nfound;
    302302    }
    303     int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, Nfields + 1, FALSE, FALSE, NULL, 0);
     303    int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", NULL, vec, Nfields + 1, FALSE, FALSE, NULL, 0);
    304304    free (vec[Nfields]->elements.Int);
    305305    free (vec[Nfields]);

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/dvo/avperiodomatch.c

    r40523 r41170  
    9999
    100100      CoordsFile = abspath("coords.fits", 1024);
    101       int status = WriteVectorTableFITS (CoordsFile, "COORDS", vec, 2, FALSE, FALSE, NULL, 0);
     101      int status = WriteVectorTableFITS (CoordsFile, "COORDS", NULL, vec, 2, FALSE, FALSE, NULL, 0);
    102102      if (!status) goto escape;
    103103      free (vec);

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/dvo/avselect.c

    r40575 r41170  
    132132
    133133      CoordsFile = abspath("coords.fits", 1024);
    134       int status = WriteVectorTableFITS (CoordsFile, "COORDS", vec, 2, FALSE, FALSE, NULL, 0);
     134      int status = WriteVectorTableFITS (CoordsFile, "COORDS", NULL, vec, 2, FALSE, FALSE, NULL, 0);
    135135      if (!status) goto escape;
    136136    }
    … …  
    327327      vec[NfieldsOut-1] = RADvec;
    328328    }
    329     int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, NfieldsOut, FALSE, FALSE, NULL, 0);
     329    int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", NULL, vec, NfieldsOut, FALSE, FALSE, NULL, 0);
    330330    if (!status) goto escape;
    331331  }

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/dvo/dvo_host_utils.c

    r40549 r41170  
    233233  // write vectors to a table (this is used by parallel dvo operations, but can be used elsewhere)
    234234  if (ResultFile) {
    235     int status = WriteVectorTableFITS (ResultFile, "RESULT", vec, Nvec, FALSE, FALSE, NULL, 0);
     235    int status = WriteVectorTableFITS (ResultFile, "RESULT", NULL, vec, Nvec, FALSE, FALSE, NULL, 0);
    236236    if (!status) {
    237237      gprint (GP_ERR, "failed to write result file %s\n", ResultFile);

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/dvo/images.c

    r40165 r41170  
    1212
    1313  off_t i, Nimage, Nmosaic;
    14   int j, status, InPic, leftside, *plist, TimeSelect, ByName;
     14  int j, status, InPic, leftside, TimeSelect, ByName;
     15  // int *plist, n, npts; -- what is/was plist?
    1516  int WITH_MOSAIC, SOLO_MOSAIC;
    1617  time_t tzero, tend;
    17   int N, NPTS, n, npts, Npts, kapa, *foundMosaic;
     18  int n, N, NPTS, Npts, kapa, *foundMosaic;
    1819  Vector Xvec, Yvec;
    1920  double r[8], d[8], x[8], y[8], Rmin, Rmax, Rmid, trange, Radius;
    … …  
    148149  int TrianglePts  = wordhash ("TRI-");
    149150
    150   npts = NPTS = 200;
     151  // npts = 200
     152  NPTS = 200;
    151153  SetVector (&Xvec, OPIHI_FLT, NPTS);
    152154  SetVector (&Yvec, OPIHI_FLT, NPTS);
    … …  
    154156  Image *mosaic = NULL;
    155157
    156   ALLOCATE (plist, int, NPTS);
    157   n = N = 0;
     158  /******** stage 1 : find the images to plot *********/
     159
     160  // ALLOCATE (plist, int, NPTS);
     161  n = 0; // tracks number of images
     162  N = 0; // tracks number of vertices
    158163  for (i = 0; i < Nimage; i++) {
    159164    if (ByName && strncmp (image[i].name, name, strlen(name))) continue;
    … …  
    335340    if (!InPic) continue;
    336341
    337     plist[n] = i;
    338     n++;
    339     if (n > npts - 1) {
    340       npts += 200;
    341       REALLOCATE (plist, int, npts);
    342     }
    343     N+=2*Npts;
     342    // plist[n] = i;
     343    // n++;
     344    // if (n > npts - 1) {
     345    //   npts += 200;
     346    //   REALLOCATE (plist, int, npts);
     347    // }
     348    n ++;
     349    N += 2*Npts;
    344350    if (N + 16 >= NPTS) {  /* need to leave room for 8 point image */
    345351      NPTS += 400;

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/dvo/init.c

    r40523 r41170  
    1414int catdir_define   PROTO((int, char **));
    1515int ccd             PROTO((int, char **));
     16int cimages         PROTO((int, char **));
    1617int cmatch          PROTO((int, char **));
    1718int cmd             PROTO((int, char **));
    … …  
    8081  {1, "catdir",      catdir_define,"re-define CATDIR"},
    8182//  {1, "ccd",         ccd,          "plot color-color diagram"},
     83  {1, "cimages",     cimages,      "fill image boxes with a color"},
    8284  {1, "cmatch",      cmatch,       "match two catalogs"},
    8385//  {1, "cmd",         cmd,          "plot cmd of stars in current region"},

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/dvo/mextract.c

    r40574 r41170  
    352352  // write vectors to a table (this is used by parallel dvo operations, but can be used elsewhere)
    353353  if (RESULT_FILE && !SKIP_RESULTS) {
    354     int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, Nreturn, FALSE, FALSE, NULL, 0);
     354    int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", NULL, vec, Nreturn, FALSE, FALSE, NULL, 0);
    355355    if (!status) goto escape;
    356356  }

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/dvo/mmatch.c

    r40523 r41170  
    148148      // XXX this is now set for both cases...
    149149      CoordsFile = abspath("coords.fits", 1024);
    150       int status = WriteVectorTableFITS (CoordsFile, "COORDS", vec, 2, FALSE, FALSE, NULL, 0);
     150      int status = WriteVectorTableFITS (CoordsFile, "COORDS", NULL, vec, 2, FALSE, FALSE, NULL, 0);
    151151      if (!status) goto escape;
    152152    }
    … …  
    349349      vec[Nfields-1] = IDXvec;
    350350    }
    351     int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", vec, Nfields, FALSE, FALSE, NULL, 0);
     351    int status = WriteVectorTableFITS (RESULT_FILE, "RESULT", NULL, vec, Nfields, FALSE, FALSE, NULL, 0);
    352352    if (!status) goto escape;
    353353  }

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/include/data.h

    r40642 r41170  
    4848  int Nx;
    4949  int Ny;
    50   float **buffers;
     50  float **flx;
     51  float **var;
    5152} MedImageType;
    5253

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/include/deimos.h

    r40642 r41170  
    77} DeimosResult;
    88
    9 float     *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *trace, int Nx, int Ny, int row);
     9float     *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, Spline *psf_trace, int Nx, int Ny, int row);
    1010void       deimos_make_kernel (float stilt, int Nx);
    1111void       deimos_free_kernel ();
    … …  
    1313
    1414// internal to make_model:
    15 float     *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, int Nx, int Ny, int row);
     15float     *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, Spline *psf_trace, int Nx, int Ny, int row);
    1616float     *deimos_apply_tilt (float *input, int Nx, int Ny);
    1717void       deimos_apply_profile (Vector *profile, float *out, int Nx, int Ny);
    1818void       deimos_add_background (opihi_flt *backgnd, float *out, int Nx, int Ny, int row);
    19 float     *deimos_apply_trace (Spline *trace, float *input, int Nx, int Ny, int row);
     19float     *deimos_apply_trace (Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, float *input, int Nx, int Ny, int row);
    2020

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/include/dvomath.h

    r39640 r41170  
    169169
    170170/* vector IO functions */
    171 int           WriteVectorTableFITS  PROTO((char *filename, char *extname, Vector **vec, int Nvec, int append, char *compress, char *format, int Ntile));
     171int           WriteVectorTableFITS  PROTO((char *filename, char *extname, Header *extraheader, Vector **vec, int Nvec, int append, char *compress, char *format, int Ntile));
    172172Vector      **ReadVectorTableFITS   PROTO((char *filename, char *extname, int *Nvec));
    173173

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/lib.data/MedImageOps.c

    r36679 r41170  
    3232  free (medimage[0].name);
    3333  for (i = 0; i < medimage[0].Ninput; i++) {
    34     free (medimage[0].buffers[i]);
     34    free (medimage[0].flx[i]);
     35    FREE (medimage[0].var[i]);
    3536  }
    36   free (medimage[0].buffers);
     37  free (medimage[0].flx);
     38  free (medimage[0].var);
    3739  free (medimage);
    3840}
    … …  
    7577  medimage->Nx = Nx;
    7678  medimage->Ny = Ny;
    77   ALLOCATE (medimage->buffers, float *, 1);
     79  ALLOCATE (medimage->flx, float *, 1);
     80  ALLOCATE (medimage->var, float *, 1);
    7881
    7982  medimages[N] = medimage;

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/lib.data/SplineOps.c

    r41136 r41170  
    4141    ALLOCATE (spline[0].y2, opihi_flt, spline[0].Nknots);
    4242  }
    43   memset (spline[0].xk, 0, spline[0].Nknots * sizeof(opihi_flt));
    44   memset (spline[0].yk, 0, spline[0].Nknots * sizeof(opihi_flt));
    45   memset (spline[0].y2, 0, spline[0].Nknots * sizeof(opihi_flt));
     43  if (spline[0].Nknots) {
     44    memset (spline[0].xk, 0, spline[0].Nknots * sizeof(opihi_flt));
     45    memset (spline[0].yk, 0, spline[0].Nknots * sizeof(opihi_flt));
     46    memset (spline[0].y2, 0, spline[0].Nknots * sizeof(opihi_flt));
     47  }
    4648}
    4749

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/lib.data/open_kapa.c

    r31635 r41170  
    266266int close_kapa (char *name) {
    267267
     268  struct timespec request, remain;
    268269  int N;
    269270
    … …  
    274275  }
    275276  DelKapaDevice (name);
     277
     278  // avoid blocking on waitpid, test every 1 msec, up to 500 msec
     279  request.tv_sec = 0;
     280  request.tv_nsec = 1000000;
     281
     282  // try to harvest the child PID
     283  int waitstatus = 0;
     284  int result = waitpid (-1, &waitstatus, WNOHANG);
     285  for (int i = 0; (i < 500) && (result == 0); i++) {
     286    nanosleep (&request, &remain);
     287    result = waitpid (-1, &waitstatus, WNOHANG);
     288  }
     289
     290  if ((result == -1) && (errno != ECHILD)) {
     291    fprintf (stderr, "unexpected error from waitpid (%d): programming error\n", errno);
     292  }
     293  if (result == 0) {
     294    fprintf (stderr, "child did not exit (close_kapa), timeout");
     295  }
     296
    276297  return (TRUE);
    277298}

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/lib.data/spline.c

    r39231 r41170  
    77  float dy, dx, *tmp;
    88 
     9  // spline is not valid with < 3 points
     10  if (N < 3) return;
     11
    912  ALLOCATE (tmp, float, N);
    1013
    … …  
    6972  opihi_flt dy, dx, *tmp;
    7073 
     74  // spline is not valid with < 3 points
     75  if (N < 3) return;
     76
    7177  ALLOCATE (tmp, opihi_flt, N);
    7278

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/lib.data/starfuncs.c

    r40398 r41170  
    7575  FWHMx = 2.355*sqrt (fabs(x2 / I - x*x));
    7676  FWHMy = 2.355*sqrt (fabs(y2 / I - y*y));
     77
     78  fprintf (stderr, "Mxx, Myy: %f, %f\n", x2/I - x*x, y2/I - y*y);
    7779  Sxy   = xy / I - x*y;
    7880  mag = -2.5*log10(I);

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/lib.data/style_args.c

    r40570 r41170  
    7272  }
    7373
     74  if ((graphmode[0].style == KAPA_PLOT_POLYGON) || (graphmode[0].style == KAPA_PLOT_POLYFILL)) {
     75    if ((N = get_argument (*argc, argv, "-npoint"))) {
     76      remove_argument (N, argc, argv);
     77      graphmode[0].ptype = atoi (argv[N]);
     78      remove_argument (N, argc, argv);
     79    } else {
     80      gprint (GP_ERR, "polygon & polyfill styles require number of points argument: -npoint N\n");
     81      return FALSE;
     82    }
     83  }
    7484  return (TRUE);
    7585}

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/lib.shell/VectorIO.c

    r41136 r41170  
    162162}
    163163
     164// insert a new header line, return end pointer
     165char *gfits_insert_header_line (Header *header, char *endptr, char *newline) {
     166
     167  if (0) {
     168    char temp1[32], temp2[32];
     169    memset (temp1, 0, 32);
     170    memset (temp2, 0, 32);
     171
     172    if (endptr) {
     173      memcpy (temp1, endptr, 8);
     174    }
     175    memcpy (temp2, newline, 8);
     176    fprintf (stderr, "insert: %s : %s\n", temp1, temp2);
     177  }
     178
     179  /* find the END of the header, if not supplied */
     180  if (!endptr) {
     181    endptr = gfits_header_field (header, "END", 1);
     182    if (endptr == NULL) return NULL;
     183  }
     184
     185  /* is there enough space for 1 more line? */
     186  if (header[0].datasize - (endptr - (header[0].buffer)) < 2*FT_LINE_LENGTH) {
     187    // no, so expand by a full header block (FT_RECORD_SIZE = 32*80 = 2880)
     188    header[0].datasize += FT_RECORD_SIZE;
     189    REALLOCATE (header[0].buffer, char, header[0].datasize);
     190    // re-find the "END" marker, in case new memory block is used
     191    endptr = gfits_header_field (header, "END", 1);
     192    if (endptr == NULL) return NULL;
     193    memset (endptr + FT_LINE_LENGTH, ' ', FT_RECORD_SIZE);
     194  }
     195
     196  /* push END line back 1 */
     197  memmove ((endptr + FT_LINE_LENGTH), endptr, FT_LINE_LENGTH);
     198  memset (endptr, ' ', FT_LINE_LENGTH);
     199
     200  strncpy (endptr, newline, 80);
     201  endptr += FT_LINE_LENGTH;
     202
     203  return endptr;
     204}
     205 
     206static char *rawkeywords[] = {"SIMPLE", "BITPIX", "NAXIS", "PCOUNT", "GCOUNT", "EXTEND", "EXTNAME", "BSCALE", "BZERO", "TFIELDS", "TFORM", "TZERO", "TSCAL", "        ", NULL};
     207
    164208// write a set of vectors to a FITS file (vectors names become fits column names)
    165 int WriteVectorTableFITS (char *filename, char *extname, Vector **vec, int Nvec, int append, char *compress, char *format, int Ntile) {
     209int WriteVectorTableFITS (char *filename, char *extname, Header *extraheader, Vector **vec, int Nvec, int append, char *compress, char *format, int Ntile) {
    166210 
    167211  Header rawheader;
    … …  
    191235  rawtable.header = &rawheader;
    192236  if (!WriteVectorTable (&rawtable, extname, vec, Nvec, format, (compress != NULL))) goto escape;
    193   // NOTE : for compression, the table is constructed in native by order
     237  // NOTE : for compression, the table is constructed in native byte order
    194238
    195239  FTable *outtable = &rawtable;
    … …  
    215259
    216260  if (!append) {
     261    // generate a blank PHU header
    217262    Header header;
    218263    Matrix matrix;
    … …  
    228273  }
    229274
     275  if (extraheader) {
     276    // copy keywords which are not the standard or table keywords
     277    char *buf = extraheader->buffer;
     278    char *endptr = NULL;
     279    for (int i = 0; i < extraheader->datasize; i+= FT_LINE_LENGTH, buf += FT_LINE_LENGTH) {
     280
     281      if (0) {
     282        char temp1[32];
     283        memset (temp1, 0, 32);
     284        memcpy (temp1, buf, 8);
     285        fprintf (stderr, "buffer: %s\n", temp1);
     286      }
     287
     288      for (int j = 0; rawkeywords[j] != NULL; j++) {
     289        if (!strncmp (buf, rawkeywords[j], strlen(rawkeywords[j]))) goto skip_insert;
     290      }
     291      endptr = gfits_insert_header_line (outheader, endptr, buf);
     292      if (!endptr) {
     293        gprint (GP_ERR, "failed to update FITS header with extra keywords\n");
     294        return (FALSE);
     295      }
     296    skip_insert:
     297      continue;
     298    }
     299  }
     300
    230301  // write the actual table data
    231302  gfits_fwrite_Theader (f, outheader);

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/lib.shell/multicommand.c

    r39457 r41170  
    6363
    6464  p = line;
     65  int Nline = strlen(line);
     66 
    6567  done = FALSE;
    6668  status = TRUE;
    … …  
    7476    tmpline = strncreate (p, q - p);
    7577    stripwhite (tmpline);
    76     if (*tmpline) {
     78    myAssert (tmpline, "oops");
    7779
    78       if (bufferPending) {
    79         // flush old messages
    80         ExpectMessage (server, 0.2, &message);
    81         bufferPending = FALSE;
     80    // empty command, free and continue
     81    if (*tmpline == 0) {
     82      free (tmpline);
     83      if (q == line + Nline) {
     84        done = TRUE;
     85      } else {
     86        p = q + 1;
     87        myAssert (p - line <= Nline, "oops");
    8288      }
     89      continue;
     90    }
    8391
    84       status = command (tmpline, &outline, verbose);
     92    if (bufferPending) {
     93      // flush old messages
     94      ExpectMessage (server, 0.2, &message);
     95      bufferPending = FALSE;
     96    }
    8597
    86       if (status == -1) {
    87         if (server) {
    88           // send the command to the server instead
    89           if (!SendMessage (server, "%s", outline)) {
    90             switch (errno) {
    91               case EPIPE:
    92                 gprint (GP_ERR, "server connection has died\n");
    93                 exit (32);
    94               default:
    95                 gprint (GP_ERR, "server is busy...32\n");
    96                 bufferPending = TRUE;
    97                 goto escape;
    98             }
     98    // input tmpline is freed by command
     99    status = command (tmpline, &outline, verbose);
     100
     101    if (status == -1) {
     102      if (server) {
     103        // send the command to the server instead
     104        if (!SendMessage (server, "%s", outline)) {
     105          switch (errno) {
     106            case EPIPE:
     107              gprint (GP_ERR, "server connection has died\n");
     108              exit (32);
     109            default:
     110              gprint (GP_ERR, "server is busy...32\n");
     111              bufferPending = TRUE;
     112              goto escape;
    99113          }
     114        }
    100115
    101           // receive the command exit status
    102           if (ExpectMessage (server, MSG_TIMEOUT, &message)) {
    103             switch (errno) {
    104               case EPIPE:
    105                 gprint (GP_ERR, "server connection has died\n");
    106                 exit (33);
    107               default:
    108                 gprint (GP_ERR, "server is busy...33\n");
    109                 bufferPending = TRUE;
    110                 goto escape;
    111             }
    112           } else {
    113             sscanf (message.buffer, "STATUS %d", &status);
    114           }
    115 
    116           // receive the resulting stderr
    117           if (ExpectMessage (server, MSG_TIMEOUT, &message)) {
    118             switch (errno) {
    119               case EPIPE:
    120                 gprint (GP_ERR, "server connection has died\n");
    121                 exit (34);
    122               default:
    123                 gprint (GP_ERR, "server is busy...34\n");
    124                 bufferPending = TRUE;
    125                 goto escape;
    126             }
    127           } else {
    128             gwrite (message.buffer, 1, message.Nbuffer, GP_ERR);
    129           }
    130 
    131           // receive the resulting stdout
    132           if (ExpectMessage (server, MSG_TIMEOUT, &message)) {
    133             switch (errno) {
    134               case EPIPE:
    135                 gprint (GP_ERR, "server connection has died\n");
    136                 exit (35);
    137               default:
    138                 gprint (GP_ERR, "server is busy...35\n");
    139                 bufferPending = TRUE;
    140                 goto escape;
    141             }
    142           } else {
    143             gwrite (message.buffer, 1, message.Nbuffer, GP_LOG);
     116        // receive the command exit status
     117        if (ExpectMessage (server, MSG_TIMEOUT, &message)) {
     118          switch (errno) {
     119            case EPIPE:
     120              gprint (GP_ERR, "server connection has died\n");
     121              exit (33);
     122            default:
     123              gprint (GP_ERR, "server is busy...33\n");
     124              bufferPending = TRUE;
     125              goto escape;
    144126          }
    145127        } else {
    146           // if no server is defined, we treat an unknown command as an error
    147           status = FALSE;
    148         }       
    149       }
    150     escape:
    151       if (outline != NULL) free (outline);
    152       if (!status && auto_break) done = TRUE;
     128          sscanf (message.buffer, "STATUS %d", &status);
     129        }
     130
     131        // receive the resulting stderr
     132        if (ExpectMessage (server, MSG_TIMEOUT, &message)) {
     133          switch (errno) {
     134            case EPIPE:
     135              gprint (GP_ERR, "server connection has died\n");
     136              exit (34);
     137            default:
     138              gprint (GP_ERR, "server is busy...34\n");
     139              bufferPending = TRUE;
     140              goto escape;
     141          }
     142        } else {
     143          gwrite (message.buffer, 1, message.Nbuffer, GP_ERR);
     144        }
     145
     146        // receive the resulting stdout
     147        if (ExpectMessage (server, MSG_TIMEOUT, &message)) {
     148          switch (errno) {
     149            case EPIPE:
     150              gprint (GP_ERR, "server connection has died\n");
     151              exit (35);
     152            default:
     153              gprint (GP_ERR, "server is busy...35\n");
     154              bufferPending = TRUE;
     155              goto escape;
     156          }
     157        } else {
     158          gwrite (message.buffer, 1, message.Nbuffer, GP_LOG);
     159        }
     160      } else {
     161        // if no server is defined, we treat an unknown command as an error
     162        status = FALSE;
     163      }
    153164    }
     165  escape:
     166    if (outline != NULL) free (outline);
     167    if (!status && auto_break) done = TRUE;
     168
    154169    p = q + 1;
    155170  }

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/lib.shell/startup.c

    r39926 r41170  
    107107      LOAD_RC = FALSE;
    108108    }
     109    if ((N = get_argument (*argc, argv, "--no-rc"))) {
     110      remove_argument (N, argc, argv);
     111      LOAD_RC = FALSE;
     112    }
     113    if ((N = get_argument (*argc, argv, "-norc"))) {
     114      remove_argument (N, argc, argv);
     115      LOAD_RC = FALSE;
     116    }
     117    if ((N = get_argument (*argc, argv, "-no-rc"))) {
     118      remove_argument (N, argc, argv);
     119      LOAD_RC = FALSE;
     120    }
    109121   
    110122    ONLY_INPUT = FALSE;

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/mana/Makefile

    r40655 r41170  
    2626$(SRC)/deimos_getobj.$(ARCH).o \
    2727$(SRC)/deimos_fitobj.$(ARCH).o \
    28 $(SRC)/deimos_fitalt.$(ARCH).o \
    29 $(SRC)/deimos_mkalt.$(ARCH).o \
    3028$(SRC)/deimos_mkmodel.$(ARCH).o \
    3129$(SRC)/deimos_arclines.$(ARCH).o \
     30$(SRC)/deimos_fitarc.$(ARCH).o \
     31$(SRC)/deimos_fitprofile.$(ARCH).o \
    3232$(SRC)/findrowpeaks.$(ARCH).o
    3333

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/mana/deimos.c

    r40655 r41170  
    44int deimos_fitslit (int argc, char **argv);
    55int deimos_mkobj (int argc, char **argv);
    6 int deimos_mkalt (int argc, char **argv);
    76int deimos_getobj (int argc, char **argv);
     7// int deimos_getalt (int argc, char **argv);
    88int deimos_fitobj (int argc, char **argv);
    9 int deimos_fitalt (int argc, char **argv);
     9// int deimos_fitalt (int argc, char **argv);
    1010int deimos_arclines (int argc, char **argv);
     11int deimos_fitarc (int argc, char **argv);
     12int deimos_fitprofile (int argc, char **argv);
    1113
    1214static Command deimos_commands[] = {
    13   {1, "fitobj", deimos_fitobj, "fit for object parameters"},
    14   {1, "fitalt", deimos_fitalt, "fit for object parameters using LMM"},
    15   {1, "getobj", deimos_getobj, "determine crude object parameters"},
    16   {1, "mkobj", deimos_mkobj, "make a full object image"},
    17   {1, "mkalt", deimos_mkalt, "make a full object image (uses deimos_make_object, for a test)"},
    18   {1, "mkslit", deimos_mkslit, "make a slit image"},
    19   {1, "fitslit", deimos_fitslit, "fit slit image to observed slit flux"},
    20   {1, "arclines", deimos_arclines, "detect arclines using LSF and STILT"},
     15  {1, "fitobj",     deimos_fitobj,     "fit for object parameters using LMM"},
     16  {1, "fitalt",     deimos_fitobj,     "fit for object parameters using LMM"},
     17  {1, "getobj",     deimos_getobj,     "determine crude object parameters"},
     18  {1, "getalt",     deimos_getobj,     "determine crude object parameters"},
     19  {1, "mkobj",      deimos_mkobj,      "make a full object image"},
     20  {1, "mkslit",     deimos_mkslit,     "make a slit image"},
     21  {1, "fitslit",    deimos_fitslit,    "fit slit image to observed slit flux"},
     22  {1, "arclines",   deimos_arclines,   "detect arclines using LSF and STILT"},
     23  {1, "fitarc",     deimos_fitarc,     "fit arclamp lines"},
     24  {1, "fitprofile", deimos_fitprofile, "fit slit profile"},
    2125};
    2226

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/mana/deimos_arclines.c

    r40655 r41170  
    5353  deimos_make_LSF_kernel (LSF, stilt, Nx);
    5454
     55  // FILE *Fout = NULL;
     56
    5557  // loop over the rows
    5658  for (int iy = 0; iy < Ny; iy++) {
    5759    float Fsum = 0;
     60    float Ksum = 0;
     61
     62    // if (iy == 5802) Fout = fopen ("test.5802.dat", "w");
     63    // if (iy == 5815) Fout = fopen ("test.5815.dat", "w");
     64
    5865    // cross-correlation of buffer values and kernel values
    5966    for (int ix = 0; ix < Nx; ix++) {
    … …  
    6269        if ((iy + ky) < 0) continue;
    6370        if ((iy + ky) >= Ny) continue;
     71        if (!isfinite(Fin[ix + (iy + ky)*Nx])) continue;
     72        if (!isfinite(kernel[ix + ko*Nx])) continue;
    6473        Fsum += Fin[ix + (iy + ky)*Nx] * kernel[ix + ko*Nx];
     74        Ksum += kernel[ix + ko*Nx];
     75        // if ((iy == 5815) || (iy == 5802)) {
     76        //   fprintf (Fout, "%d %d %d : %f %f\n", iy, ix, ky, Fin[ix + (iy + ky)*Nx], kernel[ix + ko*Nx]);
     77        // }
    6578      }
    6679    }
     80
     81    // if ((iy == 5815) || (iy == 5802)) fclose (Fout);
     82
    6783    coord->elements.Flt[iy] = iy;
    6884    flux->elements.Flt[iy] = Fsum;

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/mana/deimos_fitslit.c

    r40659 r41170  
    11# include "data.h"
     2# define IRLS_TOLERANCE 1e-4
     3
     4float weight_cauchy_square_flt (float x2);
     5static void fitflux (float *Fwind, float *Fprof, float *weight, int Nx, double *flux, double *sky);
    26
    37int deimos_fitslit (int argc, char **argv) {
    48
    5   // fitslit (window) (slit) (flux) (sky)
    6 
    7   // int N;
     9  // fitslit (window) (vari) (slit) (flux) (sky)
    810
    911  Vector *flux    = NULL;
    … …  
    1113
    1214  Buffer *wind    = NULL;
     15  Buffer *vari    = NULL;
    1316  Buffer *slit    = NULL;
    1417
    15   if (argc != 5) {
    16     gprint (GP_ERR, "USAGE: deimos fitslit (window) (slit) (flux) (sky)\n");
    17     gprint (GP_ERR, "  inputs:  window (observed 2D flux), slit (model 2D flux)\n");
     18  int N;
     19  int Niter = 10;
     20  if ((N = get_argument (argc, argv, "-irls-iter"))) {
     21    remove_argument (N, &argc, argv);
     22    Niter = atoi (argv[N]);
     23    remove_argument (N, &argc, argv);
     24  }
     25
     26  Buffer *model = NULL;
     27  if ((N = get_argument (argc, argv, "-model"))) {
     28    remove_argument (N, &argc, argv);
     29    if ((model = SelectBuffer (argv[N], ANYBUFFER, TRUE)) == NULL) return (FALSE);
     30    remove_argument (N, &argc, argv);
     31  }
     32
     33  if (argc != 6) {
     34    gprint (GP_ERR, "USAGE: deimos fitslit (window) (variance) (slit) (flux) (sky)\n");
     35    gprint (GP_ERR, "  inputs:  window (observed 2D flux), variance (on 2D flux), slit (model 2D flux)\n");
    1836    gprint (GP_ERR, "  outputs: flux (best-fit 1D flux), sky (best-fit 1D background)\n");
    1937    return FALSE;
    … …  
    2240  // XXX I probably should rename FindSpline as SelectSpline and give it the same behavior
    2341  if ((wind = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    24   if ((slit = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    25   if ((flux = SelectVector (argv[3], ANYVECTOR, TRUE)) == NULL) return (FALSE);
    26   if ((sky  = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
    27 
    28   // XXX confirm slit and wind have same dimensions
     42  if ((vari = SelectBuffer (argv[2], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     43  if ((slit = SelectBuffer (argv[3], OLDBUFFER, TRUE)) == NULL) return (FALSE);
     44  if ((flux = SelectVector (argv[4], ANYVECTOR, TRUE)) == NULL) return (FALSE);
     45  if ((sky  = SelectVector (argv[5], ANYVECTOR, TRUE)) == NULL) return (FALSE);
    2946
    3047  // define the output window
    … …  
    3249  int Ny = wind[0].matrix.Naxis[1];
    3350 
     51  // confirm wind, vari, slit have same dimensions
    3452  if (Nx != slit[0].matrix.Naxis[0]) { gprint (GP_ERR, "size mismatch in fitslit for window and slit\n"); return FALSE; }
    3553  if (Ny != slit[0].matrix.Naxis[1]) { gprint (GP_ERR, "size mismatch in fitslit for window and slit\n"); return FALSE; }
     54  if (Nx != vari[0].matrix.Naxis[0]) { gprint (GP_ERR, "size mismatch in fitslit for window and variance\n"); return FALSE; }
     55  if (Ny != vari[0].matrix.Naxis[1]) { gprint (GP_ERR, "size mismatch in fitslit for window and variance\n"); return FALSE; }
    3656
    3757  ResetVector (flux, OPIHI_FLT, Ny);
    3858  ResetVector (sky,  OPIHI_FLT, Ny);
    3959
     60  if (model) {
     61    if (!CreateBuffer (model, Nx, Ny, -32, 1.0, 0.0)) { gprint (GP_ERR, "error generating output model buffer\n"); return FALSE; }
     62  }
     63
    4064  float *Fwind = (float *) wind[0].matrix.buffer;
     65  float *Fvari = (float *) vari[0].matrix.buffer;
    4166  float *Fprof = (float *) slit[0].matrix.buffer;
     67  float *FmodOut = model ? (float *) model[0].matrix.buffer : NULL;
     68
     69  ALLOCATE_PTR (weight, float, Nx);
     70  ALLOCATE_PTR (rawwgt, float, Nx);
    4271
    4372  // loop over the rows
    4473  for (int iy = 0; iy < Ny; iy++) {
    4574
    46     // calculate elements of the chi-square for this row:
     75    double F = 0.0, S = 0.0;
    4776
    48     double R = 0.0, P = 0.0, P2 = 0.0, f = 0.0, fp = 0.0;
    49     double wt = 1.0;
    50 
     77    // set weight based on the variance
    5178    for (int ix = 0; ix < Nx; ix++) {
    52      
    53       float fxy = Fwind[ix + iy*Nx];
    54       float pxy = Fprof[ix + iy*Nx];
    55 
    56       R  += wt;
    57       P  += pxy*wt;
    58       P2 += pxy*pxy*wt;
    59       f  += fxy*wt;
    60       fp += fxy*pxy*wt;
     79      float vari = Fvari[ix + iy*Nx];
     80      weight[ix] = !isfinite(vari) || (fabs(vari) < 1e-6) ? 1.0 : 1.0 / vari;
     81      rawwgt[ix] = weight[ix];
    6182    }
    6283
    63     double det = R*P2 - P*P;
    64     double S = (P2*f - P*fp) / det;
    65     double F = (R*fp - P*f) / det;
     84    // calculate elements of the chi-square for this row:
     85    fitflux (&Fwind[iy*Nx], &Fprof[iy*Nx], weight, Nx, &F, &S);
     86
     87    // do an IRLS loop here to reject outliers
     88    int converged = FALSE;
     89    for (int iter = 0; (iter < Niter) && !converged; iter++) {
     90
     91      double Flast = F;
     92      double Slast = S;
     93
     94      // calculate weight modification based on distances (squared).
     95      // use modifier to calculate new weighted mean
     96      for (int ix = 0; ix < Nx; ix++) {
     97
     98        // F & S are model parameters; compare the model and observed values for each pixel
     99        float Fmodel = S + F*Fprof[ix + iy*Nx];
     100        float dV = (Fwind[ix + iy*Nx] - Fmodel);
     101        float d2 = SQ(dV) * rawwgt[ix];
     102       
     103        float Mod = weight_cauchy_square_flt (d2);
     104        weight[ix] = Mod * rawwgt[ix];
     105      }
     106      fitflux (&Fwind[iy*Nx], &Fprof[iy*Nx], weight, Nx, &F, &S);
     107
     108      float dF = fabs(F - Flast);
     109      float dS = fabs(S - Slast);
     110
     111      if ((dF < F * IRLS_TOLERANCE) && (dS < S * IRLS_TOLERANCE)) converged = TRUE;
     112    }
     113
     114    if (FmodOut) {
     115      for (int ix = 0; ix < Nx; ix++) {
     116        // F & S are model parameters; compare the model and observed values for each pixel
     117        float Fmodel = S + F*Fprof[ix + iy*Nx];
     118        float dV = (Fwind[ix + iy*Nx] - Fmodel);
     119        float d2 = SQ(dV) * rawwgt[ix];
     120        float Mod = weight_cauchy_square_flt (d2);
     121        if (Mod < 0.1) {
     122          FmodOut[ix + iy*Nx] = NAN;
     123        } else {
     124          FmodOut[ix + iy*Nx] = Fmodel;
     125        }
     126      }
     127    }
    66128
    67129    flux->elements.Flt[iy] = F;
    68130    sky->elements.Flt[iy] = S;
    69131  }
    70 
    71132  return TRUE;
    72133}
    73134
     135// weight is 1 / variance
     136// Fwind, Fprof are pointers to the start of a single row
     137// Fwind is the observed flux in the window
     138// Fprof is the slit profile
     139void fitflux (float *Fwind, float *Fprof, float *weight, int Nx, double *flux, double *sky) {
     140
     141  // calculate elements of the chi-square fit for this row:
     142
     143  double R = 0.0, P = 0.0, P2 = 0.0, f = 0.0, fp = 0.0;
     144
     145  for (int ix = 0; ix < Nx; ix++) {
     146     
     147    double fxy = Fwind[ix];
     148    double pxy = Fprof[ix];
     149
     150    // skip NAN / Inf pixels in either window or slit
     151    if (!isfinite(fxy)) continue;
     152    if (!isfinite(pxy)) continue;
     153
     154    double wt = (weight == NULL) ? 1.0 : weight[ix];
     155    R  += wt;
     156    P  += pxy*wt;
     157    P2 += pxy*pxy*wt;
     158    f  += fxy*wt;
     159    fp += fxy*pxy*wt;
     160  }
     161
     162  double det = R*P2 - P*P;
     163  *sky = (P2*f - P*fp) / det;
     164  *flux = (R*fp - P*f) / det;
     165}

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/mana/deimos_mkalt.c

    r40642 r41170  
    11# include "data.h"
    22# include "deimos.h"
     3
     4XXX: deprecated
    35
    46// Use this to test the deimos_make_object function used by deimos_fitobj
    … …  
    9294
    9395  // generate the model based on the test values
    94   float *model = deimos_make_model (objV, skyV, bckV, psf, profile, trace, Nx, Ny, 0);
     96  float *model = deimos_make_model (objV, skyV, bckV, psf, profile, trace, NULL, Nx, Ny, 0);
    9597
    9698  ResetBuffer (output, Nx, Ny, -32, 0.0, 1.0);

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/mana/deimos_mkmodel.c

    r40659 r41170  
    8686// generate a subimage of the full model for the range row - row + Ny
    8787// obj,sky,bck are subset vectors corresponding to the range row to row+Ny
    88 float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *trace, int Nx, int Ny, int row) {
    89 
    90   float *raw = deimos_make_straight_image (obj, sky, psf, Nx, Ny, row);
     88// the slit_trace_* and psf_trace refer to the full image starting at row
     89float *deimos_make_model (opihi_flt *obj, opihi_flt *sky, opihi_flt *bck, Vector *psf, Vector *profile, Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, Spline *psf_trace, int Nx, int Ny, int row) {
     90
     91  float *raw = deimos_make_straight_image (obj, sky, psf, psf_trace, Nx, Ny, row);
    9192  float *new = deimos_apply_tilt (raw, Nx, Ny);
    9293  deimos_apply_profile (profile, new, Nx, Ny);
    9394  deimos_add_background (bck, new, Nx, Ny, row);
    94   float *out = deimos_apply_trace (trace, new, Nx, Ny, row);
     95  float *out = deimos_apply_trace (slit_trace_red, slit_trace_blu, redlimit, new, Nx, Ny, row);
    9596
    9697  free (raw);
    … …  
    104105// *** generate a raw image with no slit tilt and no trace offsets.  apply the PSF to the
    105106// object flux, add in the sky flux
    106 float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, int Nx, int Ny, int row) {
     107float *deimos_make_straight_image (opihi_flt *obj, opihi_flt *sky, Vector *psf, Spline *psf_trace, int Nx, int Ny, int row) {
    107108  OHANA_UNUSED_PARAM(row);
    108109 
    … …  
    113114
    114115  // Xref is a global supplied by the user (defaults to Nx/2 if < 0)
    115   int Xoff = (int)(Xref - Npsf);
     116  int XoffRef = (int)(Xref - Npsf);
    116117
    117118  // Nx : width of output window
    … …  
    123124   
    124125    // we are generating the image for just the row range row to row+Ny
     126    // note that obj & sky are subset vectors for just this range of pixels
    125127    opihi_flt objVy = obj[iy]; // if we shift to using the row offset: obj[iy + row]
    126128    opihi_flt skyVy = sky[iy]; // if we shift to using the row offset: sky[iy + row]
    127129     
     130    // integral and fractional pixel offsets of PSF (none, by default)
     131    int dxi = 0;
     132    float dxf = 0.0, dxr = 1.0;
     133    int Xoff = XoffRef;
     134   
     135    if (psf_trace) {
     136      // the psf_trace is referenced against the full image, so we need to add row to iy:
     137      float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy + row);
     138      dxi = floor(dx);
     139      dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
     140      dxr = 1 - dxf;
     141      Xoff += dxi;
     142      // fprintf (stderr, "%d: %f : %d\n", iy + row, dx, Xoff);
     143    }
     144
     145    // if fractional offset is small, do not interpolate
     146    int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
     147
    128148    // flux = obj * PSF + sky
    129149    for (int ix = 0; ix < Nx; ix++) {
    … …  
    136156      // only add in the flux if we are in range of the PSF
    137157      if ((n >= 0) && (n < NpsfFull)) {
    138         value += objVy * psfV[n];
     158        if (doInterp) {
     159          if ((n > 0) && (n < NpsfFull)) {
     160            if (isfinite(psfV[n]) && isfinite(psfV[n-1])) {
     161              value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
     162            }
     163          }
     164          if (n == 0) {
     165            if (isfinite(psfV[n])) {
     166              value += objVy*psfV[n]*dxr;
     167            }
     168          }
     169        } else {
     170          if (isfinite(psfV[n])) {
     171            value += objVy * psfV[n];
     172          }
     173        }
    139174      }
    140175      out[ix + iy*Nx] = value;
    … …  
    219254
    220255// shift pixels in x based on the trace
    221 float *deimos_apply_trace (Spline *trace, float *input, int Nx, int Ny, int row) {
     256float *deimos_apply_trace (Spline *slit_trace_red, Spline *slit_trace_blu, float redlimit, float *input, int Nx, int Ny, int row) {
    222257
    223258  ALLOCATE_PTR (output, float, Nx*Ny);
    … …  
    227262    // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
    228263    // NOTE: spline is evaluated at full dispersion coordinate
    229     float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy + row);
     264    int iy_full = iy + row;
     265    Spline *slit_trace = (iy_full < redlimit) ? slit_trace_red : slit_trace_blu;
     266    float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy_full);
    230267     
    231268    // extract the integer pixel offset and the fractional offset

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/mana/deimos_mkobj.c

    r40659 r41170  
    66
    77  // input parameters:
    8   // * trace       : spline fit of slit central x pos vs y-coord
    9   // * profile     : slit window profile (vector)
    10   // * object      : vector of object flux vs y-coord
    11   // * sky         : vector of local sky signal vs y-coord
    12   // * background  : vector of extra-slit background flux vs y-coord
    13   // * PSF         : point-spread function vector (flux normalized, x-dir)
    14   // * LSF         : sline-spread function vector (flux normalized, y-dir)
    15   // * stilt       : slit tilt response : 2D kernel?
     8  // * slit_trace_red : spline fit of slit central x pos vs y-coord for red chip
     9  // * slit_trace_blu : spline fit of slit central x pos vs y-coord for blu chip
     10  // * psf_trace      : spline fit of psf central x pos vs y-coord
     11  // * profile        : slit window profile (vector)
     12  // * object         : vector of object flux vs y-coord
     13  // * sky            : vector of local sky signal vs y-coord
     14  // * backgnd        : vector of extra-slit background flux vs y-coord
     15  // * PSF            : point-spread function vector (flux normalized, x-dir)
     16  // * LSF            : sline-spread function vector (flux normalized, y-dir)
     17  // * stilt          : slit tilt response : 2D kernel?
     18
     19  // * redlimit       : pixel break between red and blu chips
     20  // * Nwave          : number of wavelength pixels in output image
    1621
    1722  int N;
    1823
    1924  // if any of these are not defined, they will have assumed identity values
    20   Spline *trace   = NULL;
    21   Vector *profile = NULL;
    22   Vector *object  = NULL;
    23   Vector *sky     = NULL;
    24   Vector *backgnd = NULL;
    25 
    26   Vector *psf     = NULL;
    27   Vector *lsf     = NULL;
     25  Spline *slit_trace_red = NULL;
     26  Spline *slit_trace_blu = NULL;
     27  Spline *psf_trace  = NULL;
     28  Vector *profile    = NULL;
     29  Vector *object     = NULL;
     30  Vector *sky        = NULL;
     31  Vector *backgnd    = NULL;
     32                     
     33  Vector *psf        = NULL;
     34  Vector *lsf        = NULL;
    2835  // add tilt later
    2936
    30   Buffer *output  = NULL;
     37  Buffer *output     = NULL;
    3138
    3239  float stilt = 0.0; // angle of the slit
    3340
    34   if ((N = get_argument (argc, argv, "-trace"))) {
    35     remove_argument (N, &argc, argv);
    36     if ((trace = FindSpline (argv[N])) == NULL) return (FALSE);
     41  // for a red vs blu spline, we need to specify the split point
     42  // XXX this is REALLY ad-hoc for Deimos.  not sure how to make
     43  // this more generic (need to define the ranges somewhere)
     44  int redlimit = 4096;
     45  if ((N = get_argument (argc, argv, "-redlimit"))) {
     46    remove_argument (N, &argc, argv);
     47    redlimit = atoi (argv[N]);
     48    remove_argument (N, &argc, argv);
     49  }
     50
     51  if ((N = get_argument (argc, argv, "-slit-trace-red"))) {
     52    remove_argument (N, &argc, argv);
     53    if ((slit_trace_red = FindSpline (argv[N])) == NULL) return (FALSE);
     54    remove_argument (N, &argc, argv);
     55  }
     56  if ((N = get_argument (argc, argv, "-slit-trace-blu"))) {
     57    remove_argument (N, &argc, argv);
     58    if ((slit_trace_blu = FindSpline (argv[N])) == NULL) return (FALSE);
     59    remove_argument (N, &argc, argv);
     60  }
     61  if ((N = get_argument (argc, argv, "-psf-trace"))) {
     62    remove_argument (N, &argc, argv);
     63    if ((psf_trace = FindSpline (argv[N])) == NULL) return (FALSE);
    3764    remove_argument (N, &argc, argv);
    3865  }
    … …  
    81108  }
    82109
     110  int VERBOSE = FALSE;
     111  if ((N = get_argument (argc, argv, "-v"))) {
     112    remove_argument (N, &argc, argv);
     113    VERBOSE = TRUE;
     114  }
     115
    83116  if (argc != 3) {
    84     gprint (GP_ERR, "USAGE: deimos mkobj (buffer) (Ncross) [-Nwave N] [-object vector] [-sky vector] [-backgnd vector] [-trace spline] [-profile vector]\n");
     117    gprint (GP_ERR, "USAGE: deimos mkobj (buffer) (Ncross) [-object vector] [-sky vector] [-backgnd vector] [-slit-trace-red spline] [-slit-trace-blu spline] [-psf-trace spline] [-profile vector] [-psf vector] [-lsf vector] [-stilt angle] [-Nwave Npixel] [-redlimit pixel]\n");
     118    gprint (GP_ERR, "  slit_trace_red : spline fit of slit central x pos vs y-coord for red chip\n");
     119    gprint (GP_ERR, "  slit_trace_blu : spline fit of slit central x pos vs y-coord for blu chip\n");
     120    gprint (GP_ERR, "  psf_trace      : spline fit of psf central x pos vs y-coord\n");
     121    gprint (GP_ERR, "  profile        : slit window profile (vector)\n");
     122    gprint (GP_ERR, "  object         : vector of object flux vs y-coord \n");
     123    gprint (GP_ERR, "  sky            : vector of local sky signal vs y-coord\n");
     124    gprint (GP_ERR, "  backgnd        : vector of extra-slit background flux vs y-coord\n");
     125    gprint (GP_ERR, "  PSF            : point-spread function vector (flux normalized, x-dir)\n");
     126    gprint (GP_ERR, "  LSF            : sline-spread function vector (flux normalized, y-dir)\n");
     127    gprint (GP_ERR, "  stilt          : slit tilt response : 2D kernel? \n");
     128    gprint (GP_ERR, "  redlimit       : pixel break between red and blu chips\n");
     129    gprint (GP_ERR, "  Nwave          : number of wavelength pixels in output image\n");
    85130    return FALSE;
    86131  }
    … …  
    123168  // if we are appying a trace offset spline, we need to generate an output window which
    124169  // is Nx + the full swing of the trace, then window back down
    125   if (trace) {
     170  if (slit_trace_red && slit_trace_blu) {
    126171    float dXmin = +1000;
    127172    float dXmax = -1000;
    128173    for (int iy = 0; iy < Ny; iy++) {
    129174      // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
    130       float dx = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
     175      Spline *slit_trace = (iy < redlimit) ? slit_trace_red : slit_trace_blu;
     176      float dx = spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy);
    131177      dXmin = MIN (dx, dXmin);
    132178      dXmax = MAX (dx, dXmax);
    … …  
    181227    if (psf) Npof -= Npsf;
    182228
     229    // integral and fractional pixel offsets of PSF (none, by default)
     230    int dxi = 0;
     231    float dxf = 0.0, dxr = 1.0;
     232   
     233    if (psf_trace) {
     234      float dx = spline_apply_dbl (psf_trace->xk, psf_trace->yk, psf_trace->y2, psf_trace->Nknots, iy);
     235      dxi = floor(dx);
     236      dxf = dx - dxi;  // -1.7 -> +0.3, -0.5 -> +0.5, +0.5 ->+0.5, +1.7 -> +0.7
     237      dxr = 1 - dxf;
     238      Npof += dxi;
     239      // fprintf (stderr, "%d: %f : %d\n", iy, dx, Npof);
     240      if (VERBOSE && (iy % 200 == 0)) { gprint (GP_ERR, "%d : %f : %d : %f %f\n", iy, dx, dxi, dxf, dxr); }
     241    }
     242
     243    // if fractional offset is small, do not interpolate
     244    int doInterp = fabs(dxf) < 1e-5 ? FALSE : TRUE;
     245
    183246    // flux = obj * PSF + sky
    184247    for (int ix = 0; ix < Nx; ix++) {
    … …  
    191254        // only add in the flux if we are in range of the PSF
    192255        if ((n >= 0) && (n < NpsfFull)) {
    193           value += objVy * psfV[n];
     256          if (doInterp) {
     257            if ((n > 0) && (n < NpsfFull)) {
     258              if (isfinite(psfV[n]) && isfinite(psfV[n-1])) {
     259                value += objVy*psfV[n]*dxr + objVy*psfV[n-1]*dxf;
     260              }
     261            }
     262            if (n == 0) {
     263              if (isfinite(psfV[n])) {
     264                value += objVy*psfV[n]*dxr;
     265              }
     266            }
     267          } else {
     268            if (isfinite(psfV[n])) {
     269              value += objVy*psfV[n];
     270            }
     271          }
    194272        }
    195273      } else {
    … …  
    329407
    330408  // shift pixels in x based on the trace
    331   if (trace) {
     409  if (slit_trace_red && slit_trace_blu) {
    332410    ALLOCATE_PTR (outTraceBuffer, char, NxBase*Ny*sizeof(float));
    333411    float *outTrace = (float *) outTraceBuffer;
    … …  
    337415
    338416      // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
    339       float dx = -spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
     417      Spline *slit_trace = (iy < redlimit) ? slit_trace_red : slit_trace_blu;
     418      float dx = -spline_apply_dbl (slit_trace->xk, slit_trace->yk, slit_trace->y2, slit_trace->Nknots, iy);
    340419     
    341420      // extract the integer pixel offset and the fractional offset

  • branches/eam_branches/ipp-20191011/Ohana/src/opihi/mana/deimos_mkslit.c

    r40659 r41170  
    77  // input parameters:
    88  //   slit profile response : vector of fractional flux vs x-coord
    9   //   trace                 : spline fit of slit central x pos vs y-coord
    10   //   flux                  : vector of input signal vs y-coord
     9  //   trace_ref,trace_blu   : spline fit of slit central x pos vs y-coord
     10  //   flux (optional)       : vector of input signal vs y-coord
     11  //   sky  (optional)       : vector of input signal vs y-coord
     12  //   fluxbins (optional)   : size of output image in y-direction
     13  //   NOTE: one of flux, sky, or fluxbins must be defined
     14
     15  // output : buff (an image with Nx defined by profile and Ny defined
     16  //                by flux, sky, or fluxbins)
    1117
    1218  int N;
    1319
    14   Vector *profile = NULL;
    15   Vector *flux    = NULL;
    16   Vector *sky     = NULL;
    17   Spline *trace   = NULL;
    18   Buffer *buff    = NULL;
     20  Vector *profile   = NULL;
     21  Vector *flux      = NULL;
     22  Vector *sky       = NULL;
     23  Spline *trace_red = NULL;
     24  Spline *trace_blu = NULL;
     25  Buffer *buff      = NULL;
    1926
    2027  // user-specified flux value (otherwise assumed to be 1)
    … …  
    4047  }
    4148
    42   if (argc != 4) {
    43     gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace) (buffer) [-flux vector] [-fluxbins N] [-sky vector]\n");
     49  // either flux or fluxbins must be specified to define output size
     50  int profilebins = 0;
     51  if ((N = get_argument (argc, argv, "-profilebins"))) {
     52    remove_argument (N, &argc, argv);
     53    profilebins = atoi (argv[N]);
     54    remove_argument (N, &argc, argv);
     55  }
     56
     57  // for a red vs blu spline, we need to specify the split point
     58  // XXX this is REALLY ad-hoc for Deimos.  not sure how to make
     59  // this more generic (need to define the ranges somewhere)
     60  int redlimit = 4096;
     61  if ((N = get_argument (argc, argv, "-redlimit"))) {
     62    remove_argument (N, &argc, argv);
     63    redlimit = atoi (argv[N]);
     64    remove_argument (N, &argc, argv);
     65  }
     66
     67  if (argc != 5) {
     68    gprint (GP_ERR, "USAGE: deimos mkslit (profile) (trace_red) (trace_blu) (buffer) [-flux vector] [-fluxbins N] [-sky vector]\n");
    4469    return FALSE;
    4570  }
    … …  
    5176
    5277  // XXX I probably should rename FindSpline as SelectSpline and give it the same behavior
    53   if ((profile = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
    54   if ((trace   = FindSpline (argv[2])) == NULL) return (FALSE);
    55   if ((buff    = SelectBuffer (argv[3], ANYBUFFER, TRUE)) == NULL) return (FALSE);
     78  if ((profile   = SelectVector (argv[1], OLDVECTOR, TRUE)) == NULL) return (FALSE);
     79  if ((trace_red = FindSpline (argv[2])) == NULL) return (FALSE);
     80  if ((trace_blu = FindSpline (argv[3])) == NULL) return (FALSE);
     81  if ((buff      = SelectBuffer (argv[4], ANYBUFFER, TRUE)) == NULL) return (FALSE);
    5682
    5783  // define the output window
    58   int Nx = 2*profile[0].Nelements + 1;
     84  // if no profile half-width is specified, use supplied profile
     85  int Nx = profilebins ? 2*profilebins + 1 : 2*profile[0].Nelements + 1;
    5986  int Ny = flux ? flux[0].Nelements : fluxbins;
    6087 
    … …  
    73100
    74101    // evaluate the trace spline at this y-coord to find the x-coord offset of the profile center
     102    Spline *trace = (iy < redlimit) ? trace_red : trace_blu;
    75103    float dx = spline_apply_dbl (trace->xk, trace->yk, trace->y2, trace->Nknots, iy);
    76104   
    … …  
    90118
    91119      // set the sky in the entire window
    92       out[ix + iy*Nx] = sky ? sky->elements.Flt[iy] : 0.0;
     120      out[ix + iy*Nx] = NAN;
    93121
    94122      // equivalent coord in the profile:
    … …  
    98126
    99127      // a default value:
     128      float vsky = sky ? sky->elements.Flt[iy] : 0.0;
    100129      float vout = NAN;
    101130
    … …  
    121150      if (flux) vout *= flux->elements.Flt[iy];
    122151
    123       out[ix + iy*Nx] += vout;
     152      out[ix + iy*Nx] = vsky + vout;
    124153    }
    125154  }

  • branches/eam_branches/ipp-20191011/Ohana/src/relphot/src/reload_catalogs.c

    r40386 r41170  
    6262
    6363    catalog.catflags    = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_MISSING | DVO_LOAD_SECFILT;
     64    // if we are going to update the format, we should update all tables
     65    if (UPDATE_CATFORMAT) {
     66      catalog.catflags |= DVO_LOAD_LENSING | DVO_LOAD_LENSOBJ | DVO_LOAD_STARPAR | DVO_LOAD_GALPHOT;
     67    }
    6468    catalog.Nsecfilt    = GetPhotcodeNsecfilt ();               // set the desired number in case we need to create the catalog
    6569

  • branches/eam_branches/ipp-20191011/Ohana/src/relphot/src/relphot_objects.c

    r40291 r41170  
    7373
    7474    catalog.catflags    = DVO_LOAD_AVERAGE | DVO_LOAD_MEASURE | DVO_LOAD_MISSING | DVO_LOAD_SECFILT;
     75    // if we are going to update the format, we should update all tables
     76    if (UPDATE_CATFORMAT) {
     77      catalog.catflags |= DVO_LOAD_LENSING | DVO_LOAD_LENSOBJ | DVO_LOAD_STARPAR | DVO_LOAD_GALPHOT;
     78    }
    7579    catalog.Nsecfilt    = GetPhotcodeNsecfilt ();
    7680
    … …  
    162166    // we can optionally convert output format here
    163167    // but it would be better to define a dvo crawler program to do this
    164     // catalog.catformat = DVO_FORMAT_PS1_V1; 
     168    // catalog.catformat = DVO_FORMAT_PS1_V1;
     169
     170    // modify the output format as desired (ignore current format on disk)
     171    if (UPDATE_CATFORMAT) {
     172      catalog.catformat = dvo_catalog_catformat (UPDATE_CATFORMAT);
     173    }
     174   
    165175    SetProtect (TRUE);
    166176    if (!dvo_catalog_save (&catalog, VERBOSE)) { fprintf (stderr, "ERROR: failed to save %s\n", catalog.filename); exit (1); }
Note: See TracChangeset for help on using the changeset viewer.