viewer3.c

su 的源代码库

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

/* VIEWER3: $Revision: 1.5 $ ; $Date: 1997/07/31 17:10:19 $	*/

#include "par.h" /*this line must appear first*/
#include "GL/glut.h"
#include "MGL/trackball.h"
#include "MGL/viewer3.h"

/*********************** self documentation *****************************/
char *sdoc[] = {
"									",
" VIEWER3 - VIEWER for 3-dimensional model				", 
"									",
" viewer3 [parameters] < hzfile		 				",
" 									",
" Optional Parameters:							",	
" 									",
" hue=1		=1 for hue =0 for black/white				",
" q=-0.6,0.06,-0.06,0.8 define the quaternion			   	",
" tbs=0.8 	the lager the trackball size is, the slower it rotates  ",
" verbose=0     =1 print some useful information			",
" rayfile=NULL	ray path file to read in 				",
" wffile=NULL	wavefront file to read in				",
" sttfile=NULL  surface traveltime file to read in                      ",
"									",
" Menu: The rightmost mouse button, with cursor placed on the graphics  ",
"       window activates the following menu:				",
" Quit:			quits the application				",
" Full Screen:		expands graphics window to full screen		",
" White/Color Rays:		selects colored (default) or white rays	", 
" Plot Rays:			show the rays (default: off)		",
" Surface Traveltimes:	show surface traveltimes (default: off)		",
" Wired or Solid WFs:	style of wavefront display wireframe (default)	",
"			or solid					",	
" Plot Wavefronts:	turn on wavefronts default: off			",
" TRI or TETRA or LAYER or HORZ:	toggles display from wireframe	",
" 			triangulated horizons, to wireframe tetrahedra, ",
"			to solid layers, to solid horizons		",
" Layer1:		toggle layer on or off (default is on)		",
" ....									",
" 									",
" Notes:								",
" In general, the tetrahedra model should be generated by tetramod	",
" (tetrahedra modle builder) which outputs a file called hzfile, so	",
" that viewer3 can read input from stdin.				",
"									",
" If, in addition, you have 3D raypath information (rayfile), wavefronts",
" (wffile) and surface traveltimes (sttfile), typically generated by    ",
" sutetraray (tetrahedral wavefront construction ray tracing),	  	",
" then rays, wavefronts and surface traveltimes can be displayed as well",
" as the model.							 	",
"									",
" The plot may be rotated by depressing the leftmost mouse button	",
" and dragging the cursor in the desired direction of rotation. The speed",
" of rotation is controlled by the speed of the cursor.			",
" 									",
" The plot may be rescaled by depressing the shift key, while also	",
" depressing and dragging the cursor up or down on the display window.	",
" 									",
NULL};

/*
 * Credits:
 *  	CWP: Zhaobo Meng, 1996
 */
/**************** end self doc *******************************************/

int 
main(int argc, char **argv)
{
      int ix;	/*index for nx*/
      int iy;	/*index for ny*/
      int nx;	/*number of sampels in horizon*/
      int ny;	/*number of sampels in horizon*/

      float xmin,xmax;
      float ymin,ymax;
      float zmin,zmax;

      float ***databot;	/*data for plotting*/
      float ***datatop;
      float ***emisbot; /*color on top horizon*/
      float ***emistop; /*color right above base horizon*/
      float v0;

      int verbose;    /*if =1 print some useful information*/

      float eyez;

      int ihz;	/*index for interfaces*/

      int *ntris;	/*number of triangles*/
      int nt;		/*number of samples in each ray*/
      int iray; 	/*index for nrays*/
      int it;		/*index for nt*/
      int iw,iwf,nwf;
      int is,ns;        /*number of sources*/

      float q0[4];

      char *rayfile=""; /*ray file*/
      char *wffile="";
      char *sttfile="";

      FILE *rayfp=NULL;
      FILE *wffp=NULL;
      FILE *sttfp=NULL;

      Layer *horz;
      Layer *ray;
      Layer *wf;
      float vmin=99999.0;
      float vmax=0.0;

      float tt; /*debugging information in the ray file*/

      int itri;

      char names[10];

      int iflag;	/*flag: =1 means ray effective*/
      float emission[4];
      float tmax=0.0,tmin=FLT_MAX;
      float ***stt=NULL;
      float **ttt=NULL;
      int ntr;

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

      /* get parameters */
      if (!getparint("verbose",&verbose)) verbose=0;

      /******************************************
      Read model parameters from hzfile
      ******************************************/
      fread(&nhz,sizeof(int),1,stdin);
      fread(&nx,sizeof(int),1,stdin);
      fread(&ny,sizeof(int),1,stdin);
      fread(&xmin,sizeof(float),1,stdin);
      fread(&xmax,sizeof(float),1,stdin);
      fread(&ymin,sizeof(float),1,stdin);
      fread(&ymax,sizeof(float),1,stdin);
      fread(&zmin,sizeof(float),1,stdin);
      fread(&zmax,sizeof(float),1,stdin);

      if (verbose)
            fprintf(stderr,"xmin=%f\nxmax=%f\nymin=%f\nymax=%f\nzmin=%f\nzmax=%f\n",
	          xmin,xmax,ymin,ymax,zmin,zmax);

      if (getparstring("rayfile",&rayfile))  
	    if ((rayfp=fopen(rayfile,"r"))==NULL)
		  err("Can not open rayfile %s",rayfile);

      if (getparstring("wffile",&wffile))
	    if ((wffp=fopen(wffile,"r"))==NULL)
		  err("Can not open wffile %s",wffile);

      if (getparstring("sttfile",&sttfile))
	    if ((sttfp=fopen(sttfile,"r"))==NULL)
		  err("Can not open sttfile %s",sttfile);

      if (!getparfloat("tbs",&tbs))   tbs=0.8;
      if (!getparint("hue",&glb_hue)) glb_hue=1; /*1 for glb_hue*/

      if (verbose) 
	    warn("nhz=%d, nx=%d, ny=%d\n",nhz,nx,ny);

      glb_on_or_off=(enum On_or_Off *)ealloc1int(3*nhz+6);
      for (ihz=0;ihz<nhz;ihz++) glb_on_or_off[ihz]=ON;

      horz=(Layer *)alloc1float(sizeof(Layer)*(nhz+1));

      /*********************************************************
      Do not use GLUT_INDEX, which gives no image;
      GLUT_SINGLE will cause redrawing every time you retate it;
      *********************************************************/
      glutInit(&argc, argv);
      glutInitWindowSize(768,768);
      glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH);
      glutCreateWindow("viewer3");
      glutDisplayFunc(redraw);
      glutIdleFunc(NULL);

      if (!getparfloat("q",q0)){
	    q0[0]=-0.6; 
	    q0[1]=0.05;
	    q0[2]=-0.06;
	    q0[3]=0.8;
      }

      normalize_quat(q0);

      curquat[0]=q0[0];
      curquat[1]=q0[1];
      curquat[2]=q0[2];
      curquat[3]=q0[3];

      glutReshapeFunc(myReshape);
      glutVisibilityFunc(vis);
      glutMouseFunc(mouse);
      glutMotionFunc(motion);
      glutCreateMenu(controlLights);
      glutAddMenuEntry("Quit",-1);
      glutAddMenuEntry("Full Screen",0);      
      glutAddMenuEntry("White/Color Rays",1);
      glutAddMenuEntry("Plot Rays",2);
      glutAddMenuEntry("Surface Traveltimes",3);
      glutAddMenuEntry("Wired or Solid WFs",4);
      glutAddMenuEntry("Plot Wavefronts",5);
      glutAddMenuEntry("TRI or TETRA or LAYER or HORZ",6);
      for (ihz=0;ihz<nhz;ihz++) {
	    sprintf(names,"Layer %d",ihz+1);
	    glutAddMenuEntry(names,ihz+7);
      }
      glutAttachMenu(GLUT_RIGHT_BUTTON);

      glShadeModel(GL_SMOOTH);
      glEnable(GL_DEPTH_TEST);
      glEnable(GL_LIGHTING);

      eyez=25;
      glMatrixMode(GL_PROJECTION);
      gluPerspective(
            40.0,   /*fovy: view angle in y direction*/
            1.0,    /*aspect: ratio of width (x) to y (height)*/
            eyez-DIAMETER,  /*near clipping plane*/
            eyez+DIAMETER); /*far clipping plane*/

      glMatrixMode(GL_MODELVIEW);
      gluLookAt(
            0.0, 0.0, eyez, /*(eyex,eyey,eyez): the eye position*/
            0.0, 0.0, 0.0,  /*(centerx,centery,centerz): the center*/
            0.0, 1.0, 0.0); /*(upx,upy,upz): the up direction*/
      glPushMatrix();

      /*the order that tetramod uses is like this*/
      for (ihz=0;ihz<nhz;ihz++) {
            fprintf(stderr,"reading horizon information %d\n",ihz);
	    /**********************************************************
	    input the horizon information from file hzfile:
	    **********************************************************/
	    horz[ihz].x=ealloc2float(nx,ny);
	    horz[ihz].y=ealloc2float(nx,ny);
	    horz[ihz].z=ealloc2float(nx,ny);
	    horz[ihz].v0=ealloc2float(nx,ny);
	    horz[ihz].v1=ealloc2float(nx,ny);

            fprintf(stderr,"read horz[%d].x...\n",ihz);
	    if (fread(horz[ihz].x[0],sizeof(float),nx*ny,
		  stdin)!=nx*ny)
		  err("Can not read x to stdin");
           
            fprintf(stderr,"read horz[%d].y...\n",ihz);
	    if (fread(horz[ihz].y[0],sizeof(float),nx*ny,
		  stdin)!=nx*ny)
		  err("Can not read y to stdin");

	    fprintf(stderr,"read horz[%d].z...\n",ihz);
	    if (fread(horz[ihz].z[0],sizeof(float),nx*ny,
		  stdin)!=nx*ny)
		  err("Can not read z to stdin");

	    fprintf(stderr,"read horz[%d].v0...\n",ihz);
	    if (fread(horz[ihz].v0[0],sizeof(float),nx*ny,stdin)!=
		  nx*ny)
		  err("Can not read v0 to stdin");
 
	    fprintf(stderr,"read horz[%d].v1...\n",ihz);
	    if (fread(horz[ihz].v1[0],sizeof(float),nx*ny,stdin)!=
	      	  nx*ny)
		  err("Can not read v1 to stdin");

	    for (iy=0;iy<ny;iy++) {
	          for (ix=0;ix<nx;ix++) {
		        vmin=MIN(vmin,horz[ihz].v0[iy][ix]);
			vmax=MAX(vmax,horz[ihz].v0[iy][ix]);
		        vmin=MIN(vmin,horz[ihz].v1[iy][ix]);
			vmax=MAX(vmax,horz[ihz].v1[iy][ix]);
		  }
	    }
      }

      if (verbose)
	    fprintf(stderr,"vmin=%f, vmax=%f\n",vmin,vmax);

      horz[nhz].x=ealloc2float(nx,ny);
      horz[nhz].y=ealloc2float(nx,ny);
      horz[nhz].z=ealloc2float(nx,ny);

      fprintf(stderr,"assign horz[%d].x,y,z\n",nhz);
      for (ix=0;ix<nx;ix++) {
	    for (iy=0;iy<ny;iy++) {
	       	  horz[nhz].x[iy][ix]=horz[nhz-1].x[iy][ix];
		  horz[nhz].y[iy][ix]=horz[nhz-1].y[iy][ix];
		  horz[nhz].z[iy][ix]=zmax;
	    }
      }

      databot=ealloc3float(3,nx,ny);
      emisbot=ealloc3float(4,nx,ny);
      datatop=ealloc3float(3,nx,ny);
      emistop=ealloc3float(4,nx,ny);

      for (ihz=0;ihz<nhz;ihz++) {
	    fprintf(stderr,"assigning datatop for ihz=%d\n",ihz);
	    for (ix=0;ix<nx;ix++) {
		  for (iy=0;iy<ny;iy++) {

			datatop[iy][ix][0]=(
			      (horz[ihz].x[iy][ix]-xmin)/
			      (xmax-xmin)-0.5)*DIAMETER;
			datatop[iy][ix][1]=-(
			      (horz[ihz].y[iy][ix]-ymin)/
			      (ymax-ymin)-0.5)*DIAMETER;
			datatop[iy][ix][2]=(
			      (horz[ihz].z[iy][ix]-zmin)/
			      (zmax-zmin)-0.5)*DIAMETER;
			v0=horz[ihz].v0[iy][ix];

			vEmission(v0,vmin,vmax,emistop[iy][ix]);
		  }
	    }

	    fprintf(stderr,"assigning databot for ihz=%d\n",ihz);
	    for (ix=0;ix<nx;ix++) {
		  for (iy=0;iy<ny;iy++) {
			databot[iy][ix][0]=(
			      (horz[ihz+1].x[iy][ix]-xmin)
			      /(xmax-xmin)-0.5)*DIAMETER;
			databot[iy][ix][1]=-(
			      (horz[ihz+1].y[iy][ix]-ymin)
			      /(ymax-ymin)-0.5)*DIAMETER;	
			databot[iy][ix][2]=(
			      (horz[ihz+1].z[iy][ix]-zmin)
			      /(zmax-zmin)-0.5)*DIAMETER;
                                
			v0=horz[ihz].v1[iy][ix];
			vEmission(v0,vmin,vmax,emisbot[iy][ix]);
		  }
	    }

	    showLayer(ihz,databot,datatop,nx,ny,emisbot,emistop);
	    showHorz(ihz,datatop,nx,ny,emistop);
	    showTetra(ihz,databot,datatop,nx,ny,emisbot,emistop);
	    showTri(ihz,datatop,nx,ny,emistop);
      }

      free3float(databot);
      free3float(datatop);
      free3float(emisbot);
      free3float(emistop);

      /*******************************************************************
      The ray positions are generated by sutetraray, named by rayfile. 
      This part will be ignored if rayfile not specified.
      ********************************************************************/
      if (rayfp!=NULL) {

            fscanf(rayfp,
                 "%d =Number of shots\n",&ns);
            fprintf(stderr,"ns=%d\n",ns);

            if (ns<=0 || ns>100) {
                 ns=0;
                 rayfp=NULL;
            }

            ray=(Layer *)alloc1float(sizeof(Layer)*ns);

	    tmax=0.0;
            for (is=0;is<ns;is++) {
                  fscanf(rayfp,
                        "%d =Maximum number of segments\n",&nt);

                  fprintf(stderr,"%d =Maximum number of segments\n",nt);
                  fscanf(rayfp,
                        "%d =Number of rays\n",&ray[is].nrays);
		  fprintf(stderr,"%d =Number of rays\n",ray[is].nrays);

	          ray[is].x=ealloc2float(ray[is].nrays,nt);
	          ray[is].y=ealloc2float(ray[is].nrays,nt);
	          ray[is].z=ealloc2float(ray[is].nrays,nt);
                  ray[is].v0=ealloc2float(ray[is].nrays,nt);
		  ray[is].nseg=ealloc1int(ray[is].nrays);

                  for (iray=0;iray<ray[is].nrays;iray++) {
                        fscanf(rayfp,"%d=nseg %f=ttotal\n",&ray[is].nseg[iray],&tt);

                        if (nt<ray[is].nseg[iray]) err("nt should >=ray[is].nseg[iray]");
                        for (it=0;it<ray[is].nseg[iray];it++) {
                              fscanf(rayfp,"%f %f %f %f %f\n",
                                    &ray[is].x[it][iray],
                                    &ray[is].y[it][iray],
                                    &ray[is].z[it][iray],
				    &ray[is].v0[it][iray],&tt);
			      tmax=MAX(tmax,ray[is].v0[it][iray]);
                        }
                        
                        ray[is].z[ray[is].nseg[iray]-1][iray]=
			      MAX(0.001,ray[is].z[ray[is].nseg[iray]-1][iray]);
           
	                for (it=0;it<ray[is].nseg[iray];it++) {
		      	      ray[is].x[it][iray]=((ray[is].x[it][iray]-xmin)/
		       	            (xmax-xmin)-0.5)*DIAMETER;
		    	      ray[is].y[it][iray]=-((ray[is].y[it][iray]-ymin)/
		       	            (ymax-ymin)-0.5)*DIAMETER;
		       	      ray[is].z[it][iray]=((ray[is].z[it][iray]-zmin)/
		       	            (zmax-zmin)-0.5)*DIAMETER;
	                }
                  }
            }
	    fclose(rayfp);

            /*white rays*/
            glNewList(nhz*4+3,GL_COMPILE);
            emission[0]=emission[1]=emission[2]=emission[3]=1.0;
            glMaterialfv(GL_FRONT_AND_BACK,GL_EMISSION,emission);
            for (is=0;is<ns;is++) {
                  for (iray=0;iray<ray[is].nrays;iray++) {
                        iflag=0;
	                glBegin(GL_LINE_STRIP);
	                for (it=0;it<ray[is].nseg[iray];it++) {
		              if (fabs(ray[is].x[it][iray])<RADIUS &&
		                  fabs(ray[is].y[it][iray])<RADIUS &&
		                  fabs(ray[is].z[it][iray])<RADIUS) {
			            glVertex3f(ray[is].x[it][iray],ray[is].y[it][iray],
			                  ray[is].z[it][iray]);

                                    iflag=1;
		              } else if (iflag) break; /*once good, now bad*/    
	                }
	                glEnd();
                  }