segyread.c

这是一个读SEGY文件的源码

	itr = 0;
	while (itr < trmax) {
		int nread;

		if (buff) {
			if (-1 == 
			   (nread = (int) read(tapefd, (char *) &tapetr, nsegy))){
				if (verbose)
				      warn("tape read error on trace %d", itr);
				if (++errcount > errmax)
				      err("exceeded maximum io errors");
			} else { /* Reset counter on successful tape IO */
				errcount = 0;
			}
		} else {
			nread = (int) fread((char *) &tapetr, 1, nsegy, tapefp);
			if (ferror(tapefp)) {
				if (verbose)
				      warn("tape read error on trace %d", itr);
				if (++errcount > errmax)
				      err("exceeded maximum io errors");
				clearerr(tapefp);
			} else { /* Reset counter on successful tape IO */
				errcount = 0;
			}
		}
			
	if (!nread) break; /* middle exit loop instead of mile-long while */
	
		/* Convert from bytes to ints/shorts */
		tapesegy_to_segy(&tapetr, &tr);
	
		/* If little endian machine, then swap bytes in trace header */
		if (endian==0)
			for (i = 0; i < SEGY_NKEYS; ++i) swaphval(&tr,i);
	
		/* Check tr.ns field */
		if (!nsflag && ns != tr.ns) {
			warn("discrepant tr.ns = %d with tape/user ns = %d\n"
				"\t... first noted on trace %d",
				tr.ns, ns, itr + 1);
			nsflag = cwp_true;
		}

		/* loop over key fields and remap */
		for (ikey=0; ikey<nkeys; ++ikey) {
			
			/* get header values */
			
			ugethval(type1[ikey], &val1, 
				 type2[ikey], ubyte[ikey]-1, 
				 (char*) &tapetr, endian, conv);
			puthval(&tr, index1[ikey], &val1);
		}

		/* Convert and write desired traces */
		if (++itr >= trmin) {

			/* Convert IBM floats to native floats */
			  if (conv) {
				switch (bh.format) {
				case 1:
			  /* Convert IBM floats to native floats */
					ibm_to_float((int *) tr.data,
						(int *) tr.data, ns, endian);
					break;
				case 2:
			  /* Convert 4 byte integers to native floats */
					long_to_float((long *) tr.data,
						(float *) tr.data, ns, endian);
					break;
				case 3:
			  /* Convert 2 byte integers to native floats */
					short_to_float((short *) tr.data,
						(float *) tr.data, ns, endian);
					break;
				case 5:
			  /* IEEE floats.  Byte swap if necessary. */
                                        if (endian==0)
                                           for (i = 0; i < ns ; ++i)
                                                swap_float_4(&tr.data[i]);
					break;
				case 8:
			  /* Convert 1 byte integers to native floats */
					integer1_to_float((signed char *)tr.data,
					(float *) tr.data, ns);
                                        break;
				}

				/* Apply trace weighting. */
				if (trcwt && tr.trwf!=0) {
			  		float scale = pow(2.0,-tr.trwf);
				  	int i;
				  	for (i=0; i<ns; ++i) {
				    		tr.data[i] *= scale;
				  	}
				}
			} else if (conv==0) {
				/* don't convert, if not appropriate */

                                switch (bh.format) {
                                case 1: /* endian=0 byte swapping */
				case 5:
                                        if (endian==0)
                                           for (i = 0; i < ns ; ++i)
                                                swap_float_4(&tr.data[i]);
                                        break;
                                case 2: /* convert longs to floats */
                                        /* SU has no provision for reading */
                                        /* data as longs */
                                        long_to_float((long *) tr.data,
                                                (float *) tr.data, ns, endian);
                                        break;
                                case 3: /* shorts are the SHORTPAC format */
					/* used by supack2 and suunpack2 */
                                        if (endian==0)/* endian=0 byte swap */
                                           for (i = 0; i < ns ; ++i)
                                                swap_short_2((short *) &tr.data[i]);
                                        /* Set trace ID to SHORTPACK format */
                                        tr.trid = SHORTPACK;
                                        break;
                                case 8: /* convert bytes to floats */
                                        /* SU has no provision for reading */
                                        /* data as bytes */
					integer1_to_float((signed char *)tr.data,
					(float *) tr.data, ns);
                                        break;
                                }
                        }

			/* Write the trace to disk */
			tr.ns = ns;
			puttr(&tr);

			/* Echo under verbose option */
			if (verbose && itr % vblock == 0)
				warn(" %d traces from tape", itr);
		}
	}



	/* Re-iterate error in case not seen during run */
	if (nsflag) warn("discrepancy found in header and trace ns values\n"
		"the value (%d) was used to extract traces", ns);


	/* Clean up (binary & header files already closed above) */
	(buff) ? eclose(tapefd):
		 efclose(tapefp);
	if (verbose) warn("tape closed successfully");


	return(CWP_Exit());
}

#ifdef _HPUX_SOURCE
static void ibm_to_float(int from[], int to[], int n, int endian)
{
	/* HP version of ibm_to_float */

    register int fconv, fmant, i, t, dummy;

	dummy = endian;

    for (i=0;i<n;++i) {
        fconv = from[i];

        /* next lines modified (M.J.Rutty 20/9/92) */
        /* if (fconv) { */
        /* fmant = 0x00ffffff & fconv; */

        fmant = 0x00ffffff & fconv;
        if (!fmant)
          fconv=0;
        else {
          /* end modifications */
            t = (int) ((0x7f000000 & fconv) >> 22) - 130;
            while (!(fmant & 0x00800000)) { --t; fmant <<= 1; }
            if (t > 254) fconv = (0x80000000 & fconv) | 0x7f7fffff;
            else if (t <= 0) fconv = 0;
            else fconv = (0x80000000 & fconv) |(t << 23)|(0x007fffff & fmant);
        }
        to[i] = fconv;
    }
    return;
}

#else /* use the regular ibm_to_float routine */

static void ibm_to_float(int from[], int to[], int n, int endian)
/***********************************************************************
ibm_to_float - convert between 32 bit IBM and IEEE floating numbers
************************************************************************
Input::
from		input vector
to		output vector, can be same as input vector
endian		byte order =0 little endian (DEC, PC's)
			    =1 other systems 
************************************************************************* 
Notes:
Up to 3 bits lost on IEEE -> IBM

Assumes sizeof(int) == 4

IBM -> IEEE may overflow or underflow, taken care of by 
substituting large number or zero

Only integer shifting and masking are used.
************************************************************************* 
Credits: CWP: Brian Sumner,  c.1985
*************************************************************************/
{
    register int fconv, fmant, i, t;

    for (i=0;i<n;++i) {

	fconv = from[i];

	/* if little endian, i.e. endian=0 do this */
	if (endian==0) fconv = (fconv<<24) | ((fconv>>24)&0xff) |
		((fconv&0xff00)<<8) | ((fconv&0xff0000)>>8);

	if (fconv) {
	    fmant = 0x00ffffff & fconv;
	    /* The next two lines were added by Toralf Foerster */
	    /* to trap non-IBM format data i.e. conv=0 data  */
	    if (fmant == 0)
		warn("mantissa is zero data may not be in IBM FLOAT Format !");	
	    t = (int) ((0x7f000000 & fconv) >> 22) - 130;
	    while (!(fmant & 0x00800000)) { --t; fmant <<= 1; }
	    if (t > 254) fconv = (0x80000000 & fconv) | 0x7f7fffff;
	    else if (t <= 0) fconv = 0;
	    else fconv = (0x80000000 & fconv) |(t << 23)|(0x007fffff & fmant);
	}
	to[i] = fconv;
    }
    return;
}
#endif

static void tapebhed_to_bhed(const tapebhed *tapebhptr, bhed *bhptr)
/****************************************************************************
tapebhed_to_bhed -- converts the seg-y standard 2 byte and 4 byte
	integer header fields to, respectively, the
	machine's short and int types. 
*****************************************************************************
Input:
tapbhed		pointer to array of 
*****************************************************************************
Notes:
The present implementation assumes that these types are actually the "right"
size (respectively 2 and 4 bytes), so this routine is only a placeholder for
the conversions that would be needed on a machine not using this convention.
*****************************************************************************
Author: CWP: Jack  K. Cohen, August 1994
****************************************************************************/

{
	register int i;
	Value val;
	
	/* convert binary header, field by field */
	for (i = 0; i < BHED_NKEYS; ++i) {
		gettapebhval(tapebhptr, i, &val);
		putbhval(bhptr, i, &val);
	}
}

static void tapesegy_to_segy(const tapesegy *tapetrptr, segy *trptr)
/****************************************************************************
tapesegy_to_segy -- converts the seg-y standard 2 byte and 4 byte
		    integer header fields to, respectively, the machine's
		    short and int types. 
*****************************************************************************
Input:
tapetrptr	pointer to trace in "tapesegy" (SEG-Y on tape) format

Output:
trptr		pointer to trace in "segy" (SEG-Y as in	 SU) format
*****************************************************************************
Notes:
Also copies float data byte by byte.  The present implementation assumes that
the integer types are actually the "right" size (respectively 2 and 4 bytes),
so this routine is only a placeholder for the conversions that would be needed
on a machine not using this convention.	 The float data is preserved as
four byte fields and is later converted to internal floats by ibm_to_float
(which, in turn, makes additonal assumptions).
*****************************************************************************
Author: CWP:Jack K. Cohen,  August 1994
****************************************************************************/
{
	register int i;
	Value val;
	
	/* convert header trace header fields */
	for (i = 0; i < SEGY_NKEYS; ++i) {
		gettapehval(tapetrptr, i, &val);
		puthval(trptr, i, &val);
	}

	/* copy the optional portion */
	memcpy((char *)&(trptr->otrav)+2, tapetrptr->unass, 60);

	/* copy data portion */
	memcpy(trptr->data, tapetrptr->data, 4*SU_NFLTS);
}

static void long_to_float(long from[], float to[], int n, int endian)
/****************************************************************************
Author:	J.W. de Bruijn, May 1995
****************************************************************************/
{
	register int i;
  
	if (endian == 0) {
		for (i = 0; i < n; ++i) {
			swap_long_4(&from[i]);
			to[i] = (float) from[i];
		}
	} else {
		for (i = 0; i < n; ++i) {
			to[i] = (float) from[i];
		}
	}
}

static void short_to_float(short from[], float to[], int n, int endian)
/****************************************************************************
short_to_float - type conversion for additional SEG-Y formats
*****************************************************************************
Author: Delft: J.W. de Bruijn, May 1995
Modified by: Baltic Sea Reasearch Institute: Toralf Foerster, March 1997 
****************************************************************************/
{
	register int i;
  
	if (endian == 0) {
		for (i = n-1; i >= 0 ; --i) {
			swap_short_2(&from[i]);
			to[i] = (float) from[i];
		}
	} else { 
		for (i = n-1; i >= 0 ; --i)
			to[i] = (float) from[i];
	}
}

static void integer1_to_float(signed char from[], float to[], int n)
/****************************************************************************
integer1_to_float - type conversion for additional SEG-Y formats
*****************************************************************************
Author: John Stockwell,  2005
****************************************************************************/
{
  	while (n--) {
		to[n] =  from[n];
	}
}

void ugethval(cwp_String type1, Value *valp1, 
	      char type2, int ubyte, 
	      char *ptr2, int endian, int conv)
{	double dval1 = 0;
	char   c = 0;
	short  s = 0;
	int    l = 0;
	float  f = 0.0;
	char	*ptr1;
	
	#if 0
	fprintf(stderr, "start ugethval %d %c\n", ubyte, type2);
	#endif
	
	switch (type2) {
	case 'b':
		ptr1 = (char*) &c;
		ptr1[0] = ptr2[ubyte];
		dval1 = (double) c;
	break;
	case 's':
		ptr1 = (char*) &s;
		ptr1[0] = ptr2[ubyte];
		ptr1[1] = ptr2[ubyte+1];
		if (endian==0)
		   swap_short_2(&s);
		dval1 = (double) s;
	break;
	case 'l':
		ptr1 = (char*) &l;
		ptr1[0] = ptr2[ubyte];
		ptr1[1] = ptr2[ubyte+1];
		ptr1[2] = ptr2[ubyte+2];
		ptr1[3] = ptr2[ubyte+3];
		if (endian==0)
		   swap_long_4((long *)&l);
		dval1 = (double) l;
	break;
	case 'f':
		ptr1 = (char*) &f;
		ptr1[0] = ptr2[ubyte];
		ptr1[1] = ptr2[ubyte+1];
		ptr1[2] = ptr2[ubyte+2];
		ptr1[3] = ptr2[ubyte+3];
		if (conv)
		   ibm_to_float((int*) &f, (int*) &f, 1, endian);
		else if (conv==0 && endian==0)
		   swap_float_4(&f);
		dval1 = (double) f;
	break;
	default:
		err("unknown type %s", type2);
	break;
	}

	#if 0
	fprintf(stderr, "value %lf\n", dval1);
	#endif	
	
	switch (*type1) {
	case 's':
		err("can't change char header word");
	break;
	case 'h':
		valp1->h = (short) dval1;
	break;
	case 'u':
		valp1->u = (unsigned short) dval1;
	break;
	case 'l':
		valp1->l = (long) dval1;
	break;
	case 'v':
		valp1->v = (unsigned long) dval1;
	break;
	case 'i':
		valp1->i = (int) dval1;
	break;
	case 'p':
		valp1->p = (unsigned int) dval1;
	break;
	case 'f':
		valp1->f = (float) dval1;
	break;
	case 'd':
		valp1->d = (double) dval1;
	break;
	default:
		err("unknown type %s", type1);
	break;
	}
}