vipick.c

seismic software,very useful

	if(icdppre==0) icdppre = 1;

	/* read in reference velocity curve if any */
        if(getparstring("vgfile",&vgfile)) {
		vgfp = fopen(vgfile,"r");
		ierr = fgetghdr(vgfp,&gh);
		if(!getparfloat("fcdpvg",&fcdpvg)) {
			if(ierr==0) {
				getgval(&gh,"ocdp2",&fcdpvg);
			} else {
				fcdpvg=1;
			}
		}
		if(!getparfloat("dcdpvg",&dcdpvg)) {
			if(ierr==0) {	
				getgval(&gh,"dcdp2",&dcdpvg);
			} else {
				dcdpvg=1;
			}
		}
		if(!getparint("nzvg",&nzvg)) {
			if(ierr==0) {
				nzvg = (int) (float)(gh.n1/gh.scale+0.5);
			} else {
				err(" nzvg must be specified ");
			}
		}
		if(!getparfloat("fzvg",&fzvg)) {
			if(ierr==0) {
				getgval(&gh,"o1",&fzvg);
			} else {
				fzvg=0.;
			}
		}
		if(!getparfloat("dzvg",&dzvg)) {
			if(ierr==0) {
				getgval(&gh,"d1",&dzvg);
			} else {
				err(" dzvg must be specified ");
			}
		}


		if(!getparint("ivgtype",&ivgtype)) ivgtype = 0;
		
		if(ierr==0) {
			nxvg = (int) (float) (gh.n2/gh.scale + 0.5);
		} else {
			fseek(vgfp,0L,2);
			nxvg = ftell(vgfp)/(nzvg*sizeof(float));
		}

	        tmp =  (ppos-fcdpvg)/dcdpvg; 	
		indxvg  = tmp;
		
		if(indxvg<0) {
			indxvg = 0;
		} else if ( indxvg >= nxvg ) {
			indxvg = nxvg-1;
		} 

printf("fcdpvg=%g dcdpvg=%g cdp=%d indxvg=%d \n",fcdpvg,dcdpvg,ppos,indxvg+1);

		zr = (float*) malloc(nzvg*sizeof(float));
		vr = (float*) malloc(nzvg*sizeof(float));

		fseek(vgfp,indxvg*nzvg*sizeof(float),0);
		fread(vr,sizeof(float),nzvg,vgfp);

		/* convert interval velocity to average velocity */
		if(ivgtype==0) {
			zr[0] = vr[0] * fzvg;
			for(iz=1;iz<nzvg;iz++) {
				zr[iz]  = zr[iz-1] + vr[iz] * dzvg;
			}
			for(iz=1;iz<nzvg;iz++) vr[iz] = zr[iz]/(fzvg+iz*dzvg);
		}

		/* compute depth of reference velocity input */
		for(iz=0;iz<nzvg;iz++) zr[iz] = fzvg + iz*dzvg;
		nr = nzvg;
	} else {
		zr = (float*) malloc(1*sizeof(float));
		vr = (float*) malloc(1*sizeof(float));
	        nr=0;	
	}

	/* determine the vifile (automatic picks) from velan */
	if (nzifile>0) {
		if(!getparfloat("avplo",&avplo)) avplo=80.;
		if(!getparfloat("avphi",&avphi)) avphi=120.;

		for(ih=0;ih<nzifile;ih++) {

			if(nr==0) {
				i2min = 0; 
				i2max = n2-1;	
				vifile[ih] = vmin;
			} else {
				
				tmp = (zifile[ih]-fzvg)/dzvg;
				indx = tmp;
				if(indx<0) {
					indx=0;
				} else if (indx>nzvg-1) {
					indx = nzvg-1;
				}
				vifile[ih] = vr[indx]; 

				tmp = vifile[ih] * avplo/100.;
				tmp = (tmp - f2)/d2; 
				i2min = tmp;
				if(i2min<0) {
					i2min=0;
				} else if (i2min>n2-1) {
					i2min = n2-1;
				}
				tmp = vifile[ih] * avphi/100.;
				tmp = (tmp - f2)/d2; 
				i2max = tmp;
				if(i2max<0) {
					i2max=0;
				} else if (i2max>n2-1) {
					i2max = n2-1;
				}
			}

			if(avplo!=avphi) {  
				tmp = (zifile[ih]-f1)/d1;
				indx = tmp;
				if (indx<0) indx=0;
				if (indx>n1) indx=n1-1;
				tmp = 0.;

				for(i2=i2min;i2<=i2max;i2++) {
					if(tmp<z[i2*n1+indx]) {
						tmp = z[i2*n1+indx];
						vifile[ih] = f2 + i2*d2;
					}
				}
			}
			if(vifile[ih]<vmin && vmin>0.) vifile[ih]=vmin;
			if(vifile[ih]>vmax && vmax>0.) vifile[ih]=vmax;
			zp[ih] = zifile[ih];
			vp[ih] = vifile[ih];
		}
		np = nzifile;
	}

	/* adjust x1beg and x1end to fall on sampled values */
	i1beg = NINT((x1beg-f1)/d1);
	i1beg = MAX(0,MIN(n1-1,i1beg));
	x1beg = f1+i1beg*d1;
	i1end = NINT((x1end-f1)/d1);
	i1end = MAX(0,MIN(n1-1,i1end));
	x1end = f1+i1end*d1;

	/* adjust x2beg and x2end to fall on sampled values */
	i2beg = NINT((x2beg-f2)/d2);
	i2beg = MAX(0,MIN(n2-1,i2beg));
	x2beg = f2+i2beg*d2;
	i2end = NINT((x2end-f2)/d2);
	i2end = MAX(0,MIN(n2-1,i2end));
	x2end = f2+i2end*d2;

	/* allocate space for image bytes */
	n1c = 1+abs(i1end-i1beg);
	n2c = 1+abs(i2end-i2beg);
	cz = ealloc1(n1c*n2c,sizeof(unsigned char));

	/* convert data to be imaged into signed characters */
	zscale = (wclip!=bclip)?255.0/(wclip-bclip):1.0e10;
	zoffset = -bclip*zscale;
	i1step = (i1end>i1beg)?1:-1;
	i2step = (i2end>i2beg)?1:-1;
	if (style==NORMAL) {
		for (i2c=0,i2=i2beg; i2c<n2c; i2c++,i2+=i2step) {
			czp = cz+n1c*n2c-(i2c+1)*n1c;
			for (i1c=0,i1=i1beg; i1c<n1c; i1c++,i1+=i1step) {
				zi = zoffset+z[i1+i2*n1]*zscale;
				if (zi<0.0) zi = 0.0;
				if (zi>255.0) zi = 255.0;
				*czp++ = (unsigned char)zi;
			}
		}
	} else {
		czp = cz;
		for (i1c=0,i1=i1beg; i1c<n1c; i1c++,i1+=i1step) {
			for (i2c=0,i2=i2beg; i2c<n2c; i2c++,i2+=i2step) {
				zi = zoffset+z[i1+i2*n1]*zscale;
				if (zi<0.0) zi = 0.0;
				if (zi>255.0) zi = 255.0;
				*czp++ = (unsigned char)zi;
			}
		}
	}
	free1float(z);
	
	/* initialize zoom box parameters */
	nxb = nx = (style==NORMAL ? n1c : n2c);
	nyb = ny = (style==NORMAL ? n2c : n1c);
	ixb = iyb = 0;
	czb = cz;
	x1begb = x1beg;  x1endb = x1end;
	x2begb = x2beg;  x2endb = x2end;

	/* connect to X server */
	if ((dpy=XOpenDisplay(NULL))==NULL)
		err("Cannot connect to display %s!\n",XDisplayName(NULL));
	scr = DefaultScreen(dpy);
	black = BlackPixel(dpy,scr);
	white = WhitePixel(dpy,scr);
	
	/* create window */
	win = xNewWindow(dpy,xbox,ybox,wbox,hbox,black,white,"vpick");

	/* if necessary, create private colormap with gray scale */
	if (STREQ(cmap,"gray")) {
		XSetWindowColormap(dpy,win,xCreateGrayColormap(dpy,win));
	} else if (STREQ(cmap,"hue")) {
		XSetWindowColormap(dpy,win,xCreateHueColormap(dpy,win));
	} else if (STREQ(cmap,"rgb")) {
		XSetWindowColormap(dpy,win,xCreateRGBColormap(dpy,win));
	}
	
	/* determine min and max pixels from standard colormap */
	pmin = xGetFirstPixel(dpy);
	pmax = xGetLastPixel(dpy);
		
	/* make GC for image */
	gci = XCreateGC(dpy,win,0,NULL);
	
	/* set normal event mask */
	XSelectInput(dpy,win,
		StructureNotifyMask |
		ExposureMask |
		KeyPressMask |
		PointerMotionMask |
		ButtonPressMask |
		ButtonReleaseMask |
		Button1MotionMask |
		Button2MotionMask);
	
	/* map window */
	XMapWindow(dpy,win);
					
	/* determine good size for axes box */
	xSizeAxesBox(dpy,win,
		labelfont,titlefont,style,
		&x,&y,&width,&height);
	
	/* clear the window */
	XClearWindow(dpy,win);
	
	/* note that image is out of date */
	imageOutOfDate = 1;
	savebg = 1;
	ia = 0;

	/* main event loop */
	while(True) {
		XNextEvent(dpy,&event);

		/* if window was resized */
		if (event.type==ConfigureNotify &&
			(event.xconfigure.width!=winwidth ||
			 event.xconfigure.height!=winheight)) {
			winwidth = event.xconfigure.width;
			winheight = event.xconfigure.height;
							
			/* determine good size for axes box */
			xSizeAxesBox(dpy,win,
				labelfont,titlefont,style,
				&x,&y,&width,&height);
			
			/* clear background */	
                        if(savebg) XSetWindowBackground(dpy,win,0L);

			/* clear the window */
			XClearWindow(dpy,win);
			
			/* note that image is out of date */
			imageOutOfDate = 1;

		/* else if window exposed */
		} else if (event.type==Expose) {
			
			/* clear all expose events from queue */
			while (XCheckTypedEvent(dpy,Expose,&event));
			
			/* if necessary, make new image */
			if (imageOutOfDate) {
				if (czbi!=NULL) free1(czbi);
				czbi = newInterpBytes(nxb,nyb,czb,
					width,height);
				if (image!=NULL) XDestroyImage(image);
				image = xNewImage(dpy,pmin,pmax,
					width,height,czbi);
				imageOutOfDate = 0;
			}

			/* clear background */	
                        if(savebg) XSetWindowBackground(dpy,win,0L);

			/* draw image (before axes so grid lines visible) */
			XPutImage(dpy,win,gci,image,0,0,x,y,
				image->width,image->height);
			
			/* draw axes on top of image */
			xDrawAxesBox(dpy,win,
				x,y,width,height,
				x1begb,x1endb,0.0,0.0,
				d1num,f1num,n1tic,grid1,label1,
				x2begb,x2endb,0.0,0.0,
				d2num,f2num,n2tic,grid2,label2,
				labelfont,title,titlefont,
				labelcolor,titlecolor,gridcolor,
				style);

		/* else if key down */
		} else if (event.type==KeyPress) {

			XLookupString(&event,keybuf,0,&keysym,
					&keystat);
			/* add, delete or insert picks */
			pkey =99;
			if (keysym==XK_a) {
				pkey = 0; 
			} else if ( keysym==XK_d ) {
				pkey = 1; 
			} else if ( keysym==XK_i ) {
				pkey = 2; 
			} else if (keysym==XK_p) {
				pkey = 3;
			} else if (keysym==XK_r) {
				pkey = 4;
			} else if (keysym==XK_A) {
				pkey = 5;
				/* track pointer and get boundary */
				xxb = xb;
				yyb = yb;
				wwb = wb;
				hhb = hb;
				xRubberBox(dpy,win,event,&xxb,&yyb,&wwb,&hhb);
				zptop[ia] = x1begb+(x1endb-x1begb)*
							(yyb-y)/height;
				zpbot[ia] = x1begb+(x1endb-x1begb)*
							(yyb+hhb-y)/height;
				ia = ia + 1;

			} else if (keysym==XK_D) {
				pkey = 6;
			}
			if (pkey<=6) {
				xMousePicks(dpy,win,event,style,
					vpickcolor,vrefcolor,
					mpicksfp, x,y,width,height,
					x1begb,x1endb,x2begb,x2endb,
					zp,vp,&np,zr,vr,nr,pkey,gci,
					&savebg,ivmode,
					zptop,zpbot,
					ztop,zbot,v0,dvdz,vatop,vabot,&ia,
					zpscrd+(icdppre-1)*maxp,
					vpscrd+(icdppre-1)*maxp,
					npscrd[icdppre-1],vpcrdcolor,
					zifile,vifile,nzifile,ivof);
			} else if (keysym==XK_s) {
				xMousePrint(dpy,win,event,style,
					mpicksfp, x,y,width,height,
					x1begb,x1endb,x2begb,x2endb,
					pcard,ppos,zp,vp,np,fex,
					ztop, zbot, dvdz, ia, ivof);
			} else if (keysym==XK_Q) {
			/* This is the exit from the event loop */
				break;
			} else {
				continue;
			}


		/* else if button down (1 == zoom, 2 == mouse tracking */
		} else if (event.type==ButtonPress) {
			/* if 1st button: zoom */
			if (event.xbutton.button==Button1) {
				savebg = 1;

				/* track pointer and get new box */
				xRubberBox(dpy,win,event,&xb,&yb,&wb,&hb);
			
				/* if new box has tiny width or height */
				if (wb<4 || hb<4) {
				
					/* reset box to initial values */
					x1begb = x1beg;
					x1endb = x1end;
					x2begb = x2beg;
					x2endb = x2end;
					nxb = nx;
					nyb = ny;
					ixb = iyb = 0;
					if (czb!=cz) free1(czb);
					czb = cz;
			
				/* else, if new box has non-zero width
				/* and height */
				} else {
			
					/* calculate new box parameters */
					if (style==NORMAL) {
					    zoomBox(x,y,width,height,
						    xb,yb,wb,hb,
						    nxb,ixb,x1begb,x1endb,
						    nyb,iyb,x2endb,x2begb,
						    &nxb,&ixb,&x1begb,&x1endb,
						    &nyb,&iyb,&x2endb,&x2begb);
					} else {
					    zoomBox(x,y,width,height,
						    xb,yb,wb,hb,
						    nxb,ixb,x2begb,x2endb,
						    nyb,iyb,x1begb,x1endb,
						    &nxb,&ixb,&x2begb,&x2endb,
						    &nyb,&iyb,&x1begb,&x1endb);
					}
			
					/* make new bytes in zoombox */
					if (czb!=cz) free1(czb);
					czb = ealloc1(nxb*nyb,
						sizeof(signed char));
					for (i=0,czbp=czb; i<nyb; i++) {
					    czp = cz+(iyb+i)*nx+ixb;