iipick.c

seismic software,very useful

#include "comva.h"
#include "par.h"
#include "xplot.h"
#include <X11/Xatom.h>
#include <X11/keysym.h>

/* self-documentation */
char *sdoc = 
"IIPICK - X IMAGE plot and interface pick \n"
"\n"
"iipick n1= [optional parameters] <binaryfile\n"
"\n"
"X Functionality:\n"
"Button 1	Zoom with rubberband box\n"
"Button 2	Show mouse (x1,x2) coordinates while pressed\n"
"Button 3	Same as key 'a' 				\n"
"Q key		Quit (can also use Motif Action button)\n"
"a key		Append current mouse (x1,x2) location to end of picks \n"
"i key		Insert current mouse (x1,x2) location into middle of picks \n"
"d key		Delete current mouse (x1,x2) location from picks \n"
"p key		Display picks \n"
"n key		Do not display picks (seismic section display only) \n"
"s key		Append saved mouse (x1,x2) picks to pick file \n"
"c key		Clear all saved mouse (x1,x2) picks of the current horizon\n"
"f key          go to next horizon	\n" 
"b key          go to previous horizon	\n" 
"1-9 key        go to horizon 1-9	\n"
"0 key          go to horizon 10	\n"
"F1 key         go to horizon 11	\n"
"F2 key         go to horizon 12	\n"
"Fi key         go to horizon i+10	(i=1,12) \n"
"F12 key        go to horizon 22        \n"
"T key          edit top velocity of current horizon \n"
"               (Dashed line --- previous; Solid line --- edited.	\n"
"                to save edited velocity, use S key; to exit, use Q key; \n"
"		 to edit, use right button of mouse;			\n" 
"                to show mouse location, use button 2 of mouse) 	\n" 
"B key          edit bottom velocity of current horizon \n"
"               (Dashed line --- previous; Solid line --- edited.	\n"
"                to save edited velocity, use S key; to exit, use Q key; \n"
"		 to edit, use right button of mouse; 			\n" 
"                to show mouse location, use button 2 of mouse) 	\n" 
"\n"
"Required Parameters:\n"
"n1                     number of samples in 1st (fast) dimension\n"
"\n"
"Optional Parameters:\n"
"d1=1.0                 sampling interval in 1st dimension\n"
"f1=d1                  first sample in 1st dimension\n"
"n2=all                 number of samples in 2nd (slow) dimension\n"
"d2=1.0                 sampling interval in 2nd dimension\n"
"f2=d2                  first sample in 2nd dimension\n"
"perc=100.0             percentile used to determine clip\n"
"clip=(perc percentile) clip used to determine bclip and wclip\n"
"bperc=perc             percentile for determining black clip value\n"
"wperc=100.0-perc       percentile for determining white clip value\n"
"bclip=clip             data values outside of [bclip,wclip] are clipped\n"
"wclip=-clip            data values outside of [bclip,wclip] are clipped\n"
"cmap=gray              gray, hue, or rgb colormaps may be specified\n"
"verbose=1              =1 for info printed on stderr (0 for no info)\n"
"xbox=0                 x in pixels of upper left corner of window\n"
"ybox=50                y in pixels of upper left corner of window\n"
"wbox=550               width in pixels of window\n"
"hbox=700               height in pixels of window\n"
"x1beg=x1min            value at which axis 1 begins\n"
"x1end=x1max            value at which axis 1 ends\n"
"d1num=0.0              numbered tic interval on axis 1 (0.0 for automatic)\n"
"f1num=x1min            first numbered tic on axis 1 (used if d1num not 0.0)\n"
"n1tic=1                number of tics per numbered tic on axis 1\n"
"grid1=none             grid lines on axis 1 - none, dot, dash, or solid\n"
"label1=                label on axis 1\n"
"x2beg=x2min            value at which axis 2 begins\n"
"x2end=x2max            value at which axis 2 ends\n"
"d2num=0.0              numbered tic interval on axis 2 (0.0 for automatic)\n"
"f2num=x2min            first numbered tic on axis 2 (used if d2num not 0.0)\n"
"n2tic=1                number of tics per numbered tic on axis 2\n"
"grid2=none             grid lines on axis 2 - none, dot, dash, or solid\n"
"label2=                label on axis 2\n"
"labelfont=Erg14        font name for axes labels\n"
"title=                 title of plot\n"
"titlefont=Rom22        font name for title\n"
"labelcolor=blue        color for axes labels\n"
"titlecolor=red         color for title\n"
"gridcolor=SteelBlue    color for grid lines\n"
"style=seismic          normal (axis 1 horizontal, axis 2 vertical) or\n"
"labelfont=Erg14        font name for axes labels\n"
"                       seismic (axis 1 vertical, axis 2 horizontal)\n"
"mpicks=stderr          file name to save mouse picks (ifile or hfile)\n"
"                       (where stderr is standard error output --- screen)\n"
"ifile=mpicks           old ifile or hfile name 			\n"
"xtol=d2                the horizons will extend, at both ends, a distance \n" 
"                         of xtol to avoid discontinuity		\n" 
"ism=1                  smoothing option; to laterally smooth interval  \n"
"                         velocities at top and bottom, and the interval \n"
"                         velocity gradient, along interface, before output\n" 
"                         0=average all values at picks along interface and \n"
"                           this single value is assigned to all output \n"
"                           points along the interface \n"
"                         1=no smoothing \n"
"                         n=smooting values at picks along interface using \n"
"                           n data points \n"
"VELO=NULL              average velocity VELO card and interval velocity \n"
"                         gradient DVDZ card file name              \n"
"                       (if not given, i.e., NULL, the picking will output \n"
"                        the positions of the horizons with velocity \n"
"                        fields being interpolated from those 		\n"
"                        in input if present) 				\n"
"lhnvelo=0              last interface/horizon number to use VELO for 	\n"
"                       velocity updating				\n" 
"                       (for i-th horizons,				\n" 
"                              i < lhnvelo, velocity will not be updated \n"   
"                              i >= lhnvelo, velocity updated if needed) \n"   
"dif=1                  DIF values between two ends of an interface in	\n" 
"                       the output ifile (either 1 or 0)		\n"
"itopbot=0              mode of editing top and bottom velocity of same layer\n"
"                       0=top and bottom change independently 		\n"
"                       1=top and bottom change togather 		\n"
"                         (when editing top velocity of one layer, 	\n" 
"                          bottom velocity of the same layer gets 	\n"
"			   updated automatically)			\n"
"\n"
"NOTE: \n"
" 1. picks output format will follow the hfile format with velocity at top  \n"
"    and at bottom redefined (IFILE format). Average velocity and velocity \n"
"    vertical gradient are also added in IFILE format.			\n"
" 2. if itopbot=1, depths of interfaces (horizons) in input ifile must 	\n"
"    incease in order, i.e., shallow interface before deep interfaces.	\n"
"\n"
"	AUTHOR:	Zhiming Li	      	9/91	\n"
"\n";


/* define maximum number of horizons and maximum number of picks per horizon */
#define maxhs 128
#define maxpks 256

/* functions defined and used internally */
static void zoomBox (int x, int y, int w, int h, 
	int xb, int yb, int wb, int hb,
	int nx, int ix, float x1, float x2,
	int ny, int iy, float y1, float y2,
	int *nxb, int *ixb, float *x1b, float *x2b,
	int *nyb, int *iby, float *y1b, float *y2b);
static unsigned char *newInterpBytes (int n1in, int n2in, unsigned char *bin,
	int n1out, int n2out);
void xMouseLoc(Display *dpy, Window win, XEvent event, int style, Bool show,
	int x, int y, int width, int height,
	float x1begb, float x1endb, float x2begb, float x2endb);
void xMousePrint(Display *dpy, Window win, XEvent event, int style,
	FILE *mpicksfp, int x, int y, int width, int height,
	float x1begb, float x1endb, float x2begb, float x2endb,
	float *x1picks, float *x2picks, float *veltop, float *velbot, int *dif,
	float *velavg, float *dvdz,
	int *npicks, int nhs, int ism, 
	float xmin, float xmax, float zmin, float zmax,
	int cdpxmin, int cdpxmax, float dcdp, 
	int *hnums, char hnames[maxhs][80],
	char *dunits, char *zattrib, int idif, int ihfile);
void xMousePicks(Display *dpy, Window win, XEvent event, int style,
	char hcolors[maxhs][10],
	FILE *mpicksfp, int x, int y, int width, int height,
	float x1begb, float x1endb, float x2begb, float x2endb,
	float *x1picks, float *x2picks, float *veltop, float *velbot, int *dif,
	int *npicks, int nhs, int ih, GC gc, int pkey, int *savebg);
void vtopbot(int maxp,int maxcdp, 
     	int ncdp_velo,float *x_velo,float *z_velo,float *v_velo,int *nps_velo,
	int ncdp_dvdz,float *x_dvdz,float *zt_dvdz,float *zb_dvdz, 
	float *dvdz_dvdz,int *nps_dvdz,
	int nhs, float *xpicks, float *zpicks, int *npicks, float xtol, 
	float *veltops, float *velbots, float *velavgs, float *dvdzs,
	int lhnvelo);

void ivedit(float *xpick, float *vpick, char *hcolor, int nvi, char *title,
		float xmin, float dcdp, int cdpxmin, int *isave);

void interp(float *xold, float *xnew, int nold, int nnew, float *vtop, 
	float *vbot, float *vavg, float *dvdz);

static Display *dpy;
static Window win;

main (argc,argv)
int argc; char **argv;
{
	int n1,n2,n1tic,n2tic,nfloats,
		i1,i2,grid1,grid2,style,
		n1c,n2c,i1beg,i1end,i2beg,i2end,i1c,i2c,
		nz,iz,i1step,i2step,verbose,
		xbox,ybox,wbox,hbox,
		xb,yb,wb,hb,
		x,y,width,height,
		i,j,nx,ny,nxb,nyb,ixb,iyb,
		imageOutOfDate,winwidth=-1,winheight=-1,
		showloc=0,*npicks,savebg,pkey,fex,
		cdpxmin,cdpxmax,*hnums,ih,nhs,*difs,*indx,
		cdpxmini,cdpxmaxi,ihsave,itopbot,
		ism,ncdps_velo,ncdps_dvdz,maxp,maxcdp,ihfile=0,
		*nps_velo,*nps_dvdz,*sortindex,*isort,isave,
		*cdps_velo,*cdps_dvdz,cdp, ivtopbot=0,nvi,jh,nold,idif,lhnvelo;
	float labelsize,titlesize,perc,clip,bperc,wperc,bclip,wclip,
		d1,f1,d2,f2,*z,*temp,zscale,zoffset,zi,
		x1beg,x1end,x2beg,x2end,
		x1min,x1max,x2min,x2max,
		d1num,f1num,d2num,f2num,
		x1begb,x1endb,x2begb,x2endb,
		*x1picks,*x2picks,*veltops,*velbots,
		*dens,*qval,*pois, 
		*velavgs,*dvdzs, 
		xmin,xmax,zmin,zmax,dcdp,
		xmini,xmaxi,zmini,zmaxi,dcdpi,vmini,vmaxi, xtol,
		*z_velo,*x_velo,*v_velo,*sorts,*cdpsort,
		*zt_dvdz,*zb_dvdz,*dvdz_dvdz,*x_dvdz,xcdp,*xold,
		*vtopold, *vbotold, *vavgold, *dvdzold;
	unsigned char *cz,*czp,*czb,*czbp,*czbi=NULL;
	char *label1="depth",*label2="cdp",
		*title="Seismic Interface Picking",
		*labelfont="Erg14",*titlefont="Rom22",
		*styles="seismic",*grid1s="none",*grid2s="none",
		*labelcolor="blue",*titlecolor="red",
		*gridcolor="SteelBlue",*cmap="gray",keybuf[256],*mpicks,
		hnames[maxhs][80],strhnum[3],hcolors[maxhs][10],
		*hname="NAME",*VELO="NULL",
		*dunits="METER",*zattrib="DEPTH SECTION",
		dunitsi[80], zattribi[80], *vtype, vtypei[80],
		ivtitle[80],veloyes,*ifile;
	FILE *infp=stdin, *mpicksfp, *velofp, *tty, *ifilefp;
/*
	Display *dpy;
	Window win;
*/
	XEvent event;
	KeySym keysym;
	XComposeStatus keystat;
	XImage *image=NULL;
	GC gci;
	int scr;
	unsigned long black,white,pmin,pmax;


	float cdpbegb, cdpendb;

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

	/* get parameters describing 1st dimension sampling */
	if (!getparint("n1",&n1))
		err("Must specify number of samples in 1st dimension!\n");
	d1 = 1.0;  getparfloat("d1",&d1);
	f1 = d1;  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 (fseek(infp,0L,2)!=0)
			err("must specify n2 if in a pipe!\n");
		nfloats = eftell(infp)/sizeof(float);
		efseek(infp,0L,0);
		n2 = nfloats/n1;
	}
	d2 = 1.0;  getparfloat("d2",&d2);
	f2 = d2;  getparfloat("f2",&f2);
	x2min = (d2>0.0)?f2:f2+(n2-1)*d2;
	x2max = (d2<0.0)?f2:f2+(n2-1)*d2;


	/* set up file to save mouse picks */
	if (!getparstring("mpicks", &mpicks)) {
		mpicksfp = stderr;
	} else {
		/* if file exist */
                if((mpicksfp = fopen(mpicks,"r"))!=NULL) {
			fex = 1;
                        fclose(mpicksfp);
                        mpicksfp = efopen(mpicks,"r");
                } else {
                        mpicksfp = efopen(mpicks,"w");
                        fex = 0;
                }
	}

	if (!getparstring("ifile", &ifile)) {
		ifilefp = mpicksfp;
	} else {
                ifilefp = efopen(ifile,"r");
		fex = 1;
	}

	/* memory allocations */
	x1picks = (float *) malloc(maxhs*maxpks*sizeof(float));
	x2picks = (float *) malloc(maxhs*maxpks*sizeof(float));
	veltops = (float *) malloc(maxhs*maxpks*sizeof(float));
	velbots = (float *) malloc(maxhs*maxpks*sizeof(float));
	npicks = (int *) malloc(maxhs*sizeof(int));
	hnums = (int *) malloc(maxhs*sizeof(int));
	difs = (int *) malloc(maxhs*maxpks*sizeof(int));
	dens = (float *) malloc(maxhs*maxpks*sizeof(float));
	qval = (float *) malloc(maxhs*maxpks*sizeof(float));
	pois = (float *) malloc(maxhs*maxpks*sizeof(float));
	velavgs = (float *) malloc(maxhs*maxpks*sizeof(float));
	dvdzs = (float *) malloc(maxhs*maxpks*sizeof(float));
	xold = (float *) malloc(maxpks*sizeof(float));
	indx = (int *) malloc(maxpks*sizeof(int));
	vtopold = (float *) malloc(maxpks*sizeof(float));
	vbotold = (float *) malloc(maxpks*sizeof(float));
	vavgold = (float *) malloc(maxpks*sizeof(float));
	dvdzold = (float *) malloc(maxpks*sizeof(float));

	for(ih=0;ih<maxhs;ih++) npicks[ih]=0;
	for(ih=0;ih<maxhs*maxpks;ih++) {
		veltops[ih]=0.;
		velbots[ih]=0.;
		difs[ih] = 0;
		velavgs[ih] = 0.;
		dvdzs[ih] = 0.;
	}

	/* horizon color assignment */
	strncpy(hcolors[0],"blue",4);
	strncpy(hcolors[1],"cyan",4);
	strncpy(hcolors[2],"green",5);
	strncpy(hcolors[3],"brown",5);
	strncpy(hcolors[4],"orange",6);
	strncpy(hcolors[5],"red",3);
	strncpy(hcolors[6],"violet",6);
	strncpy(hcolors[7],"purple",6);
	for(ih=8;ih<maxhs;ih++) strcpy(hcolors[ih],hcolors[ih-ih/8*8]);

	

	/* 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");

	/* 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 colormap specification */
	getparstring("cmap",&cmap);

	/* get axes parameters */
	xbox = 0; getparint("xbox",&xbox);
	ybox = 50; getparint("ybox",&ybox);
	wbox = 550; getparint("wbox",&wbox);
	hbox = 700; getparint("hbox",&hbox);
	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("style",&styles);
	if (STREQ("normal",styles)) style = NORMAL;
	else style = SEISMIC;
	getparstring("titlecolor",&titlecolor);
	getparstring("labelcolor",&labelcolor);
	getparstring("gridcolor",&gridcolor);

	if(!getparint("cdpxmin",&cdpxmin)) cdpxmin=1;
	if(!getparint("cdpxmax",&cdpxmax)) cdpxmax=n2;

	if(!getparfloat("xmin",&xmin)) xmin=f2;
	if(!getparfloat("xmax",&xmax)) xmax=f2+(n2-1)*d2;
	if(!getparfloat("zmin",&zmin)) zmin=f1;
	if(!getparfloat("zmax",&zmax)) zmax=f1+(n1-1)*d1;

	dcdp = (xmax-xmin)/(cdpxmax-cdpxmin);

	getparstring("dunits",&dunits);
	getparstring("zattrib",&zattrib);

	if(!getparint("ism",&ism)) ism=1;
	if(!getparfloat("xtol",&xtol)) xtol=d2;
	if(!getparint("dif",&idif)) idif=1;
	if(!getparint("lhnvelo",&lhnvelo)) {
		lhnvelo=0;
	} else {
		lhnvelo=lhnvelo-1;
	}
	if(idif!=1) idif=0;
	if(!getparint("itopbot",&itopbot)) itopbot=0;

	/* read in velocity card dataset, if any */
	ncdps_velo = 0;
	ncdps_dvdz = 0;
	maxp = 256;