pscube.c

su 的源代码库

/* Copyright (c) Colorado School of Mines, 2006.*/
/* All rights reserved.                       */

/* PSCUBE: $Revision: 1.12 $ ; $Date: 2001/06/12 17:41:50 $	*/


#include "par.h"
#include "psplot.h"

/*********************** self documentation **********************/
char *sdoc[] = {
" 									",
" PSCUBE - PostScript image plot with Legend of a data CUBE       ",
" 									",
" pscube n1= n2= n3= [optional parameters] <binaryfile >postscriptfile	",
"    or									",
" pscube n1= n2= n3= front= side= top= [optional parameters] >postscriptfile",
" 									",
" Data formats supported:						",
"	1. Entire cube read from stdin (n1*n2*n3 floats) [default format]",
"	2. Faces read from stdin (n1*n2 floats for front, followed by n1*n3",
"	   floats for side, and n2*n3 floats for top) [specify faces=1]	",
"	3. Faces read from separate data files [specify filenames]	",
" 									",
" Required Parameters:							",
" n1                     number of samples in 1st (fastest) dimension	",
" n2                     number of samples in 2nd dimension		",
" n3                     number of samples in 3rd (slowest) dimension	",
" 									",
" Optional Parameters:							",
" front                  name of file containing front panel		",
" side                   name of file containing side panel		",
" top                    name of file containing top panel		",
" faces=0                =1 to read faces from stdin (data format 2)	",
" d1=1.0                 sampling interval in 1st dimension		",
" f1=0.0                 first sample in 1st dimension			",
" d2=1.0                 sampling interval in 2nd dimension		",
" f2=0.0                 first sample in 2nd dimension			",
" d3=1.0                 sampling interval in 3rd dimension		",
" f3=0.0                 first sample in 3rd dimension			",
" perc=100.0             percentile used to determine clip		",
" clip=(perc percentile) clip used to determine bclip and wclip		",
" bperc=perc             percentile for determining black clip value	",
" wperc=100.0-perc       percentile for determining white clip value	",
" bclip=clip             data values outside of [bclip,wclip] are clipped",
" wclip=-clip            data values outside of [bclip,wclip] are clipped",
" brgb=0.0,0.0,0.0       red, green, blue values corresponding to black	",
" wrgb=1.0,1.0,1.0       red, green, blue values corresponding to white	",
" bhls=0.0,0.0,0.0       hue, lightness, saturation corresponding to black",
" whls=0.0,1.0,0.0       hue, lightness, saturation corresponding to white",
" bps=12                 bits per sample for color plots, either 12 or 24",
" d1s=1.0                factor by which to scale d1 before imaging	",
" d2s=1.0                factor by which to scale d2 before imaging	",
" d3s=1.0                factor by which to scale d3 before imaging	",
" verbose=1              =1 for info printed on stderr (0 for no info)	",
" xbox=1.5               offset in inches of left side of axes box	",
" ybox=1.5               offset in inches of bottom side of axes box	",
" size1=4.0              size in inches of 1st axes (vertical)		",
" size2=4.0              size in inches of 2nd axes (horizontal)	",
" size3=3.0              size in inches of 3rd axes (projected)		",
" angle=45               projection angle of cube in degrees (0<angle<90)",
"                        (angle between 2nd axis and projected 3rd axis)",
" x1end=x1max            value at which axis 1 ends			",
" d1num=0.0              numbered tic interval on axis 1 (0.0 for automatic)",
" f1num=x1min            first numbered tic on axis 1 (used if d1num not 0.0)",
" n1tic=1                number of tics per numbered tic on axis 1	",
" grid1=none             grid lines on axis 1 - none, dot, dash, or solid",
" label1=                label on axis 1				",
" x2beg=x2min            value at which axis 2 begins			",
" d2num=0.0              numbered tic interval on axis 2 (0.0 for automatic)",
" f2num=x2min            first numbered tic on axis 2 (used if d2num not 0.0)",
" n2tic=1                number of tics per numbered tic on axis 2	",
" grid2=none             grid lines on axis 2 - none, dot, dash, or solid",
" label2=                label on axis 2				",
" x3end=x3max            value at which axis 3 ends			",
" d3num=0.0              numbered tic interval on axis 3 (0.0 for automatic)",
" f3num=x3min            first numbered tic on axis 3 (used if d3num not 0.0)",
" n3tic=1                number of tics per numbered tic on axis 3	",
" grid3=none             grid lines on axis 3 - none, dot, dash, or solid",
" label3=                label on axis 3				",
" labelfont=Helvetica    font name for axes labels			",
" labelsize=18           font size for axes labels			",
" title=                 title of plot					",
" titlefont=Helvetica-Bold font name for title				",
" titlesize=24           font size for title				",
" titlecolor=black       color of title					",
" axescolor=black        color of axes					",
" gridcolor=black        color of grid					",
" legend=0               =1 display the color scale                     ",
"                        if ==1, resize xbox,ybox,width,height          ",
" lstyle=vertleft       Vertical, axis label on left side               ",
"                        =vertright (Vertical, axis label on right side)",
"                        =horibottom (Horizontal, axis label on bottom) ",
" units=                 unit label for legend                          ",
" legendfont=times_roman10    font name for title                       ",
" following are defaults for lstyle=0. They are changed for other lstyles",
" lwidth=1.2             colorscale (legend) width in inches            ",
" lheight=height/3       colorscale (legend) height in inches           ",
" lx=1.0                 colorscale (legend) x-position in inches       ",
" ly=(height-lheight)/2+xybox    colorscale (legend) y-position in pixels",
" lbeg= lmin or wclip-5*perc    value at which legend axis begins       ",
" lend= lmax or bclip+5*perc    value at which legend axis ends         ",
" ldnum=0.0      numbered tic interval on legend axis (0.0 for automatic)",
" lfnum=lmin     first numbered tic on legend axis (used if d1num not 0.0)",
" lntic=1        number of tics per numbered tic on legend axis ",
" lgrid=none     grid lines on legend axis - none, dot, dash, or solid",
" 									",
" All color specifications may also be made in X Window style Hex format",
" example:   axescolor=#255						",
"									",
" Note: The values of x1beg=x1min, x2end=x2max and x3beg=x3min cannot   ",
" be changed.								",
" 									",
NULL};
/**************** end self doc ********************************/

/*
 * AUTHOR:  Craig Artley, Colorado School of Mines, 03/12/93
 * NOTE:  Original written by Zhiming Li & Dave Hale, CSM, 07/01/90
 *	  Completely rewritten, the code now bears more similarity to
 *	  psimage than the previous pscube.  Faces of cube now rendered
 *	  as three separate images, rather than as a single image.  The
 *	  output no longer suffers from stretching artifacts, and the
 *	  code is simpler.  -Craig
 * MODIFIED:  Craig Artley, Colorado School of Mines, 12/17/93
 * 	  Added color options.
 * MODIFIED: made up from pscube and psimage to generate 3-D plots 
 *  with legend  by Claudia Vanelle, Institute of Geophysics,
 *  University of Hamburg, Germany somewhen in 2000
 */

/* legend stuff added from psimage by Claudia Vanelle, 
 * Applied Geophysics Group Hamburg somewhen in 2000 */

/* color specification array indices for RGB and HLS */
#define R 0
#define G 1
#define B 2
#define H 0
#define L 1
#define S 2
#define FMAX(x,y) (float) (x) > (y) ? (x) : (y) 
#define FMIN(x,y) (float) (x) > (y) ? (y) : (x) 

/* functions defined and used internally */
static void drawimage(int hls, float colors[3][2], 
	int width, int height, int bps, float matrix[], unsigned char* z);
static float rgbvalue (float n1, float n2, float hue);
static void hlsrgb (float h, float l, float s, float *r, float *g, float *b);

int main (int argc, char **argv)
{
  int n1,n2,n3,n1s,n2s,n3s,n1c,n2c,n3c,i1,i2,i3,i1c,i2c,i3c,
    i1beg,i1end,i2beg,i2end,i3beg,i3end,i1step,i2step,i3step,
    n1tic,n2tic,n3tic,grid1,grid2,grid3,nz,iz,
    verbose,faces,hls,bps,style=SEISMIC,bbox[4],
    legend,ugrid=SOLID,lstyle=VERTLEFT,lz,lbegsup=0,lendsup=0,ln=256,
    lbbox[4];
  float d1,d2,d3,d1s,d2s,d3s,f1,f2,f3,size1,size2,size3,xbox,ybox,angle,
    x1min,x1max,x2min,x2max,x3min,x3max,
    x1beg,x1end,x2beg,x2end,x3beg,x3end,
    d1num,f1num,d2num,f2num,d3num,f3num,
    p1beg,p1end,p2beg,p2end,p3beg,p3end,
    clip,bclip,wclip,perc,bperc,wperc,
    zscale,zoffset,zi,labelsize,titlesize,
    *z,*zfront,*zside,*ztop,*temp,matrix[6],colors[3][2],
    lwidth,lheight,lx,ly,lbeg,lend,
    lmin=(float)FLT_MAX,lmax=(float)-FLT_MAX,
    ldnum,lfnum,ld,lf=0,labmatrix[6];
  unsigned char *czfront,*czside,*cztop,
    *szfront,*szside,*sztop,*czp,
    *data_legend=NULL;
  char *label1="",*label2="",*label3="",*title="",
    *labelfont="Helvetica",*titlefont="Helvetica-Bold",
    *grid1s="none",*grid2s="none",*grid3s="none",
    *titlecolor="black",*axescolor="black",*gridcolor="black",
    *frontf,*sidef,*topf,
    *units="", *legendfont="times_roman10", *lstyles="vertleft",*lgrids="none";
  FILE *infp=stdin,*frontfp,*sidefp,*topfp;
  
  /* initialize getpar */
  initargs(argc,argv);
  requestdoc(1);
  
  /* get parameters describing 1st dimension sampling */
  if (!getparint("n1",&n1)) err("must specify n1!\n");
  if (!getparfloat("d1",&d1)) d1 = 1.0;
  if (!getparfloat("f1",&f1)) f1 = 0.0;
  
  /* get parameters describing 2nd dimension sampling */
  if (!getparint("n2",&n2)) err("must specify n2!\n");
  if (!getparfloat("d2",&d2)) d2 = 1.0;
  if (!getparfloat("f2",&f2)) f2 = 0.0;
  
  /* get parameters describing 3rd dimension sampling */
  if (!getparint("n3",&n3)) err("must specify n3!\n");
  if (!getparfloat("d3",&d3)) d3 = 1.0;
  if (!getparfloat("f3",&f3)) f3 = 0.0;
  
  /* determine input type */
  if (!getparint("faces",&faces)) faces = 0;
  
  /* read color parameters */
  bps = 8;
  hls = 0;
  colors[R][0] = colors[G][0] = colors[B][0] = 0.0;
  colors[R][1] = colors[G][1] = colors[B][1] = 1.0;
  if (countparval("brgb") || countparval("wrgb")) {
    float brgb[3],wrgb[3];
    brgb[R] = brgb[G] = brgb[B] = 0.0;
    wrgb[R] = wrgb[G] = wrgb[B] = 1.0;
    getparfloat("brgb",&brgb[0]);
    getparfloat("wrgb",&wrgb[0]);
    brgb[R] = MAX(0.0,MIN(1.0,brgb[R]));
    wrgb[R] = MAX(0.0,MIN(1.0,wrgb[R]));
    brgb[G] = MAX(0.0,MIN(1.0,brgb[G]));
    wrgb[G] = MAX(0.0,MIN(1.0,wrgb[G]));
    brgb[B] = MAX(0.0,MIN(1.0,brgb[B]));
    wrgb[B] = MAX(0.0,MIN(1.0,wrgb[B]));
    colors[R][0] = brgb[R];  colors[R][1] = wrgb[R];
    colors[G][0] = brgb[G];  colors[G][1] = wrgb[G];
    colors[B][0] = brgb[B];  colors[B][1] = wrgb[B];
		if (!getparint("bps",&bps)) bps = 12;
		if (bps!=12 && bps!=24)
		  err("bps must equal 12 or 24 for color plots!\n");
  } else if (countparval("bhls") || countparval("whls")) {
    float bhls[3],whls[3];
    hls = 1;
    bhls[H] = whls[H] = 0.0;
    bhls[L] = 0.0;  whls[L] = 1.0;
    bhls[S] = whls[S] = 0.0;
    getparfloat("bhls",&bhls[0]);
    getparfloat("whls",&whls[0]);
    bhls[L] = MAX(0.0,MIN(1.0,bhls[L]));
    whls[L] = MAX(0.0,MIN(1.0,whls[L]));
    bhls[S] = MAX(0.0,MIN(1.0,bhls[S]));
    whls[S] = MAX(0.0,MIN(1.0,whls[S]));
    colors[H][0] = bhls[0];  colors[H][1] = whls[0];
    colors[L][0] = bhls[1];  colors[L][1] = whls[1];
    colors[S][0] = bhls[2];  colors[S][1] = whls[2];
    if (!getparint("bps",&bps)) bps = 12;
    if (bps!=12 && bps!=24)
      err("bps must equal 12 or 24 for color plots!\n");
  }
  
  /* allocate space */
  nz = n1*n2+n1*n3+n2*n3;
  z = ealloc1float(nz);
  zfront = z;
  zside = zfront+n1*n2;
  ztop = zside+n1*n3;
  
  /* read data */
  if (getparstring("front",&frontf)
      && getparstring("side",&sidef)
      && getparstring("top",&topf)) {
    
    /* read face files */
    if ((frontfp = fopen(frontf,"r")) == NULL)
      err("error opening front file!\n");
    if (fread(zfront,sizeof(float),n1*n2,frontfp)!=n1*n2)
      err("error reading front file!\n");
    if ((sidefp = fopen(sidef,"r")) == NULL)
      err("error opening side file!\n");
    if (fread(zside,sizeof(float),n1*n3,sidefp)!=n1*n3)
      err("error reading side file!\n");
    if ((topfp = fopen(topf,"r")) == NULL)
      err("error opening top file!\n");
    if (fread(ztop,sizeof(float),n2*n3,topfp)!=n2*n3)
      err("error reading top file!\n");
    
  } else if (getparstring("front",&frontf)
	     || getparstring("side",&sidef)
	     || getparstring("top",&topf)) {
    
    err("must specify all or none of face, side, and top!\n");
    
  } else if (faces) {
    /* read faces from stdin */
    if (fread(zfront,sizeof(float),n1*n2,infp)!=n1*n2)
      err("error reading front from input!\n");
		if (fread(zside,sizeof(float),n1*n3, infp)!=n1*n3)
		  err("error reading side from input!\n");
		if (fread(ztop,sizeof(float),n2*n3, infp)!=n2*n3)
		  err("error reading top from input!\n");
  } else {
    /* read cube from stdin, pick off faces */
    temp = ealloc1float(n1);
    for (i3=0; i3<n3; i3++) {
      for (i2=0; i2<n2; i2++) {
	if (fread(temp,sizeof(float),n1,infp)!=n1)
	  err("error reading cube from input!\n");
	if (i3==0) 
	  for (i1=0; i1<n1; i1++)
	    zfront[i1+i2*n1] = temp[i1];
	if (i2==n2-1)
	  for (i1=0; i1<n1; i1++)
	    zside[i1+i3*n1] = temp[i1];
	ztop[i2+i3*n2] = temp[0];
      }
    }
    free1float(temp);
  }
  
  /* if necessary, determine clips from percentiles */
  if (getparfloat("clip",&clip)) {