psimage.c

su 的源代码库

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

/* PSIMAGE: $Revision: 1.26 $ ; $Date: 2006/01/27 18:14:14 $	*/

#include "par.h"
#include "psplot.h"
/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" PSIMAGE - PostScript IMAGE plot of a uniformly-sampled function f(x1,x2)",
"									",
" psimage n1= [optional parameters] <binaryfile >postscriptfile		",
"									",
" Required Parameters:							",
" n1			 number of samples in 1st (fast) dimension	",
"									",
" Optional Parameters:							",
" d1=1.0		 sampling interval in 1st dimension		",
" f1=0.0		 first sample in 1st dimension			",
" n2=all		 number of samples in 2nd (slow) dimension	",
" d2=1.0		 sampling interval in 2nd dimension		",
" f2=0.0		 first sample in 2nd 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",
"                        bclip and wclip will be set to be inside       ",
"                        [lbeg,lend] if lbeg and/or lend are supplied   ",
" threecolor=1		 supply 3 color values instead of only two,	",
"                        using not only black and white, but f.e. red,	",
"                        green and blue					",
" brgb=0.0,0.0,0.0	 red, green, blue values corresponding to black	",
" grgb=1.0,1.0,1.0	 red, green, blue values corresponding to grey	",
" 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",
" ghls=0.0,1.0,0.0	 hue, lightness, saturation corresponding to grey",
" 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	",
" 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	",
" width=6.0		 width in inches of axes box			",
" height=8.0		 height in inches of axes box			",
" x1beg=x1min		 value at which axis 1 begins			",
" 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			",
" x2end=x2max		 value at which axis 2 ends			",
" 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				",
" 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					",
" axeswidth=1            width (in points) of axes                      ",
" ticwidth=axeswidth     width (in points) of tic marks			",
" gridwidth=axeswidth    width (in points) of grid lines		",
" style=seismic		 normal (axis 1 horizontal, axis 2 vertical) or	",
"			 seismic (axis 1 vertical, axis 2 horizontal)	",
" legend=0	         =1 display the color scale			",
" lnice=0                =1 nice legend arrangement                     ",
"                        (overrides ybox,lx,width,height parameters)    ",
" 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",
"									",
" curve=curve1,curve2,...  file(s) containing points to draw curve(s)   ",
" npair=n1,n2,n2,...            number(s) of pairs in each file         ",
" curvecolor=black,..	 color of curve(s)				",
" curvewidth=axeswidth	 width (in points) of curve(s)			",
"									",
" NOTES:								",
" The curve file is an ascii file with the points specified as x1 x2 	",
" pairs, one pair to a line.  A \"vector\" of curve files and curve	",
" colors may be specified as curvefile=file1,file2,etc.			",
" and curvecolor=color1,color2,etc, and the number of pairs of values   ",
" in each file as npair=npair1,npair2,... .				",
"									",
" You may eliminate the blocky appearance of psimages by adjusting the  ",
" d1s= and d2s= parameters, so that psimages appear similar to ximages.	",
"									",
" All color specifications may also be made in X Window style Hex format",
" example:   axescolor=#255						",
" 									",
" Some example colormap settings:					",
" red white blue: wrgb=1.0,0,0 grgb=1.0,1.0,1.0 brgb=0,0,1.0 		",
" white red blue: wrgb=1.0,1.0,1.0 grgb=1.0,0.0,0.0 brgb=0,0,1.0 	",
" blue red white: wrgb=0.0,0.0,1.0 grgb=1.0,0.0,0.0 brgb=1.0,1.0,1.0 	",
" red green blue: wrgb=1.0,0,0 grgb=0,1.0,0 brgb=0,0,1.0		",
" orange light-blue green: wrgb=1.0,.5,0 grgb=0,.7,1.0 brgb=0,1.0,0	",
" red light-blue dark blue: wrgb=0.0,0,1.0 grgb=0,1.0,1.0 brgb=0,0,1.0 	",
"									",
NULL};

/*
 * AUTHOR:  Dave Hale, Colorado School of Mines, 05/29/90
 * MODIFIED:  Craig Artley, Colorado School of Mines, 08/30/91
 *	    BoundingBox moved to top of PostScript output
 * MODIFIED:  Craig Artley, Colorado School of Mines, 12/16/93
 *	    Added color options (Courtesy of Dave Hale, Advance Geophysical).
 * Modified: Morten Wendell Pedersen, Aarhus University, 23/3-97
 *           Added ticwidth,axeswidth, gridwidth parameters 
 * MODIFIED: Torsten Schoenfelder, Koeln, Germany 006/07/97
 *          colorbar (legend) (as in ximage (by Berend Scheffers, Delft))
 * MODIFIED: Brian K. Macy, Phillips Petroleum, 01/14/99
 *	    Added curve plotting option
 * MODIFIED: Torsten Schoenfelder, Koeln, Germany 02/10/99
 *          color scale with interpolation of three colors
 */
/**************** end self doc ********************************/

/* 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 RIGHT 0
#define LEFT 1
#define FMAX(x,y) (float) (x) > (y) ? (x) : (y) 
#define FMIN(x,y) (float) (x) > (y) ? (y) : (x) 
#define EPS_PSPLOT .000001	/* epsilon value to stabilize legend */


/* functions defined and used internally */
static void drawimage(int hls, float colors[3][3], 
	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,n1tic,n2tic,nfloats,bbox[4],
	  i1,i2,grid1,grid2,style,
	  n1c,n2c,n1s,n2s,i1beg,i1end,i2beg,i2end,i1c,i2c,
	  nz,iz,i1step,i2step,verbose,hls,bps,
	  legend,ugrid=SOLID,lstyle=VERTLEFT,lz,lbegsup=0,lendsup=0,ln=256,
	  lbbox[4], threecolor=0; /* BEREND, Schoenfelder */
        int lnice; /* c liner */
	float labelsize,titlesize,perc,clip,bperc,wperc,bclip,wclip,
		d1,f1,d2,f2,*z,*temp,zscale,zoffset,zi,
		xbox,ybox,width,height,
		x1beg,x1end,x2beg,x2end,
		x1min,x1max,x2min,x2max,
		d1num,f1num,d2num,f2num,
		p1beg,p1end,p2beg,p2end,matrix[6],colors[3][3], /* for 3 color mode */
		d1s,d2s,
	  lwidth,lheight,lx,ly,lbeg,lend,lmin=(float) FLT_MAX,lmax=(float) -FLT_MAX,
	  ldnum,lfnum,ld,lf=0,labmatrix[6]; /* BEREND, Schoenfelder */
	float axeswidth, ticwidth, gridwidth;
	unsigned char *cz,*czp,*sz,*data_legend=NULL;
	char *label1="",*label2="",*title="",*units="",
	  *legendfont="times_roman10",
	  *labelfont="Helvetica",*titlefont="Helvetica-Bold",
	  *styles="seismic",*grid1s="none",*grid2s="none",
	  *titlecolor="black",*axescolor="black",*gridcolor="black",
	  *lstyles="vertleft",*lgrids="none";
	FILE *infp=stdin;

	float **x1curve=NULL,**x2curve=NULL,*curvewidth=NULL;
	int i,j,curve=0,*npair=NULL,ncurvecolor=0,ncurvewidth=0;
	char **curvecolor=NULL,**curvefile=NULL;
	FILE *curvefp=NULL;
	cwp_Bool is_curve = cwp_false;

	/* initialize getpar */
	initargs(argc,argv);
	requestdoc(1);

	/* get parameters describing 1st dimension sampling */
	if (!getparint("n1",&n1)) err("must specify n1!\n");
	d1 = 1.0;  getparfloat("d1",&d1);
	f1 = 0.0;  getparfloat("f1",&f1);
	x1min = (d1>0.0)?f1:f1+(n1-1)*d1;
	x1max = (d1<0.0)?f1:f1+(n1-1)*d1;

	/* get parameters describing 2nd dimension sampling */
	if (!getparint("n2",&n2)) {
		if (efseeko(infp,(off_t) 0,SEEK_END)!=0)
			err("must specify n2 if in a pipe!\n");
		nfloats = (int) (eftello(infp)/((off_t) sizeof(float)));
		efseeko(infp,(off_t) 0,SEEK_SET);
		n2 = nfloats/n1;
	}
	d2 = 1.0;  getparfloat("d2",&d2);
	f2 = 0.0;  getparfloat("f2",&f2);
	x2min = (d2>0.0)?f2:f2+(n2-1)*d2;
	x2max = (d2<0.0)?f2:f2+(n2-1)*d2;

	/* read color parameters */
	if (!getparint("threecolor",&threecolor)) threecolor=1;
	bps = 8;
	hls = 0;
	/* color[][0] is black, color[][2] is white in 2 color mode */
	colors[R][0] = colors[G][0] = colors[B][0] = 0.0;
	colors[R][1] = colors[G][1] = colors[B][1] = 0.5;
	colors[R][2] = colors[G][2] = colors[B][2] = 1.0;
	if (countparval("brgb") || countparval("wrgb")) {
		float brgb[3],grgb[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]);
		grgb[R] = (brgb[R] + wrgb[R])/2.;
		grgb[G] = (brgb[G] + wrgb[G])/2.;
		grgb[B] = (brgb[B] + wrgb[B])/2.;
		if (threecolor==1)
		  getparfloat("grgb",&grgb[0]);
		brgb[R] = MAX(0.0,MIN(1.0,brgb[R]));
		grgb[R] = MAX(0.0,MIN(1.0,grgb[R]));
		wrgb[R] = MAX(0.0,MIN(1.0,wrgb[R]));
		brgb[G] = MAX(0.0,MIN(1.0,brgb[G]));
		grgb[G] = MAX(0.0,MIN(1.0,grgb[G]));
		wrgb[G] = MAX(0.0,MIN(1.0,wrgb[G]));
		brgb[B] = MAX(0.0,MIN(1.0,brgb[B]));
		grgb[B] = MAX(0.0,MIN(1.0,grgb[B]));
		wrgb[B] = MAX(0.0,MIN(1.0,wrgb[B]));
		colors[R][0] = brgb[R];	 colors[R][1] = grgb[R];  colors[R][2] = wrgb[R];
		colors[G][0] = brgb[G];	 colors[G][1] = grgb[G];  colors[G][2] = wrgb[G];
		colors[B][0] = brgb[B];	 colors[B][1] = grgb[B];  colors[B][2] = 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],ghls[3],whls[3];
		hls = 1;
		bhls[H] = ghls[H] = whls[H] = 0.0;
		bhls[L] = 0.0;	ghls[L] = 0.5;	whls[L] = 1.0;
		bhls[S] = ghls[S] = whls[S] = 0.0;
		getparfloat("bhls",&bhls[0]);
		getparfloat("whls",&whls[0]);
		ghls[H] = (bhls[H] + whls[H])/2.;
		ghls[L] = (bhls[L] + whls[L])/2.;
		ghls[S] = (bhls[S] + whls[S])/2.;
		if (threecolor==1)
		  getparfloat("ghls",&ghls[0]);
		bhls[L] = MAX(0.0,MIN(1.0,bhls[L]));
		ghls[L] = MAX(0.0,MIN(1.0,ghls[L]));
		whls[L] = MAX(0.0,MIN(1.0,whls[L]));
		bhls[S] = MAX(0.0,MIN(1.0,bhls[S]));
		ghls[S] = MAX(0.0,MIN(1.0,ghls[S]));
		whls[S] = MAX(0.0,MIN(1.0,whls[S]));
		colors[H][0] = bhls[0];	 colors[H][1] = ghls[0];  colors[H][2] = whls[0];
		colors[L][0] = bhls[1];	 colors[L][1] = ghls[1];  colors[L][2] = whls[1];
		colors[S][0] = bhls[2];	 colors[S][1] = ghls[2];  colors[S][2] = whls[2];
		if (!getparint("bps",&bps)) bps = 12;
		if (bps!=12 && bps!=24)
			err("bps must equal 12 or 24 for color plots!\n");
	}

	/* get legend specs BEREND, Schoenfelder */
	legend = 0; getparint("legend", &legend); /* BEREND, Schoenfelder */
	getparstring("units", &units); /* BEREND, Schoenfelder */
	getparstring("legendfont", &legendfont);     /* BEREND, Schoenfelder */

	/* set up curve plotting */
	if ((curve=countparval("curve"))!=0) {
		curvefile=(char**)ealloc1(curve,sizeof(void*));
		getparstringarray("curve",curvefile);
		if ((x1curve=(float**)malloc(curve*sizeof(void*)))==NULL)
			err("Could not allocate x1curve pointers\n");
		if ((x2curve=(float**)malloc(curve*sizeof(void*)))==NULL)
			err("Could not allocate x2curve pointers\n");
		npair=ealloc1int(curve);
		getparint("npair",npair);
		is_curve = cwp_true;
	} else {
		npair=(int *)NULL;
		curvefile=(char **)NULL;
		x1curve=(float **)NULL;
		x2curve=(float **)NULL;
		is_curve = cwp_false;
	}
	if (is_curve) {
	 if ((ncurvecolor=countparval("curvecolor"))<curve) {
		curvecolor=(char**)ealloc1(curve,sizeof(void*));
		if (!getparstringarray("curvecolor",curvecolor)) {
			curvecolor[0]=(char *)cwp_strdup("black\0");
			ncurvecolor=1;
		}
		for (i=ncurvecolor; i<curve; i++)
			curvecolor[i]=(char *)cwp_strdup(curvecolor[ncurvecolor-1]);
	 } else if (ncurvecolor) {
		curvecolor=(char**)ealloc1(ncurvecolor,sizeof(void*));
		getparstringarray("curvecolor",curvecolor);
	 }
	 for (j=0; j<curve; j++) {
		curvefp=efopen(curvefile[j],"r");
		x1curve[j]=ealloc1float(npair[j]);
		x2curve[j]=ealloc1float(npair[j]);
		for (i=0; i<npair[j]; i++) {
			fscanf(curvefp,"%f",&x1curve[j][i]);
			fscanf(curvefp,"%f",&x2curve[j][i]);
		}
		efclose(curvefp);
	 }
	}

	/* read binary data to be plotted */
	nz = n1*n2;
	z = ealloc1float(nz);
	if (fread(z,sizeof(float),nz,infp)!=nz)
		err("error reading input file!\n");

	/* if necessary, determine clips from percentiles */
	if (getparfloat("clip",&clip)) {
		bclip = clip;
		wclip = -clip;
	}
	if ((!getparfloat("bclip",&bclip) || !getparfloat("wclip",&wclip)) &&
		!getparfloat("clip",&clip)) {
		perc = 100.0;  getparfloat("perc",&perc);
		temp = ealloc1float(nz);
		for (iz=0; iz<nz; iz++)
			temp[iz] = z[iz];
		if (!getparfloat("bclip",&bclip)) {
			bperc = perc;	getparfloat("bperc",&bperc);
			iz = (nz*bperc/100.0);
			if (iz<0) iz = 0;
			if (iz>nz-1) iz = nz-1;
			qkfind(iz,nz,temp);
			bclip = temp[iz];
		}
		if (!getparfloat("wclip",&wclip)) {
			wperc = 100.0-perc;  getparfloat("wperc",&wperc);
			iz = (nz*wperc/100.0);
			if (iz<0) iz = 0;
			if (iz>nz-1) iz = nz-1;
			qkfind(iz,nz,temp);
			wclip = temp[iz];
		}
		free1float(temp);
	}
	verbose = 1;  getparint("verbose",&verbose);
	if (verbose) warn("bclip=%g wclip=%g",bclip,wclip);

	/* get scaled sampling intervals */
	d1s = 1.0;  getparfloat("d1s",&d1s);
	d2s = 1.0;  getparfloat("d2s",&d2s);
	d1s = fabs(d1s);  d1s *= d1;
	d2s = fabs(d2s);  d2s *= d2;

	/* get axes parameters */
	xbox = 1.5; getparfloat("xbox",&xbox); /* if psimage is called by ximage, it */
	ybox = 1.5; getparfloat("ybox",&ybox); /* will xbox=1.166 and ybox=1.167 */
	width = 6.0; getparfloat("wbox",&width); getparfloat("width",&width);
	height = 8.0;getparfloat("hbox",&height);getparfloat("height",&height);
         /* begin c liner */
	lnice = 0;  getparint("lnice",&lnice); 
        if (lnice==1) {
            ybox = 2.2;
            /* lx=8 is set below, after getpar on lx ... c liner */
            width = 5.4;
            height = 7.2;
        }
         /* end c liner */
	x1beg = x1min; getparfloat("x1beg",&x1beg);
	x1end = x1max; getparfloat("x1end",&x1end);
	d1num = 0.0; getparfloat("d1num",&d1num);
	f1num = x1min; getparfloat("f1num",&f1num);
	n1tic = 1; getparint("n1tic",&n1tic);
	getparstring("grid1",&grid1s);
	if (STREQ("dot",grid1s))
		grid1 = DOT;
	else if (STREQ("dash",grid1s))
		grid1 = DASH;
	else if (STREQ("solid",grid1s))
		grid1 = SOLID;
	else
		grid1 = NONE;
	getparstring("label1",&label1);
	x2beg = x2min; getparfloat("x2beg",&x2beg);
	x2end = x2max; getparfloat("x2end",&x2end);
	d2num = 0.0; getparfloat("d2num",&d2num);
	f2num = 0.0; getparfloat("f2num",&f2num);
	n2tic = 1; getparint("n2tic",&n2tic);
	getparstring("grid2",&grid2s);
	if (STREQ("dot",grid2s))
		grid2 = DOT;
	else if (STREQ("dash",grid2s))
		grid2 = DASH;
	else if (STREQ("solid",grid2s))
		grid2 = SOLID;
	else
		grid2 = NONE;
	getparstring("label2",&label2);
	getparstring("labelfont",&labelfont);
	labelsize = 18.0; getparfloat("labelsize",&labelsize);
	getparstring("title",&title);
	getparstring("titlefont",&titlefont);
	titlesize = 24.0; getparfloat("titlesize",&titlesize);
	getparstring("titlecolor",&titlecolor);
	getparstring("axescolor",&axescolor);
	getparstring("gridcolor",&gridcolor);

	/* axes and tic width */
        if(!getparfloat("axeswidth",&axeswidth)) axeswidth=1;
        if (!getparfloat("ticwidth",&ticwidth)) ticwidth=axeswidth;
        if(!getparfloat("gridwidth",&gridwidth)) gridwidth =axeswidth;

	if (is_curve) {
	 if ((ncurvewidth=countparval("curvewidth"))<curve) {
		curvewidth=ealloc1float(curve);
		if (!getparfloat("curvewidth",curvewidth)) {
			curvewidth[0]=axeswidth;
			ncurvewidth=1;
		}
		for (i=ncurvewidth; i<curve; i++)
			curvewidth[i]=curvewidth[ncurvewidth-1];
	 } else {
		curvewidth=ealloc1float(ncurvewidth);
		getparfloat("curvewidth",curvewidth);
	 }
	}

	getparstring("style",&styles);

	if (STREQ("normal",styles))
		style = NORMAL;
	else
		style = SEISMIC;

	/* Get or calc legend parameters */
	/* Legend min and max: Calc from data read in */
	if (legend) {