segy.cpp

segy读取原码

#include "stdafx.h"
#include "segy.h"
#include <numeric>
#include <functional>
#include <vector>
#include <iterator>
#include <math.h>

#include <ctype.h>
using namespace std;

void * SEGY_realloc (void* ptr, size_t size, size_t oldsize)
{
#ifdef MS_DOS
	/* fix for missing realloc in MS VC6.0 */
	void *new;
	if (new=(void*)malloc(size))
	{
		memcpy(new,ptr,oldsize);
		free(ptr);
		ptr=new;
	}
	return new;
#else
	return realloc(ptr,size);
#endif
}

class ShotPoint
{
public:
    int X;
    int Y;
	int Trace;
    ShotPoint() : X(0), Y(0), Trace(0) {}
    ShotPoint(int NewX, int NewY, int NewTrace) : X(NewX), Y(NewY), Trace(NewTrace) {}
};

bool operator==(const ShotPoint& a, const ShotPoint& b)
{
    return (a.X == b.X) && (a.Y == b.Y);
}

bool operator<(const ShotPoint& a, const ShotPoint& b)
{
	if (a.X == b.X) return (a.Y < b.Y);
    return a.X < b.X;
}


static void SwapWordByteOrder (char * data_array, int word_count)
{
  int  i, j;
  char tmpchar;
  for (i=0; i<word_count; i++)  {
      j = (i*4);
      tmpchar         = data_array[j+0];
      data_array[j+0] = data_array[j+3];
      data_array[j+3] = tmpchar;
      tmpchar         = data_array[j+2];
      data_array[j+2] = data_array[j+1];
      data_array[j+1] = tmpchar;
  }
}


static void SwapNybbleByteOrder (char * data_array, int nybble_count)
{
  int  i, j;
  char tmpchar;
  for (i=0; i<nybble_count; i++)  {
      j = (i*2);
      tmpchar         = data_array[j+0];
      data_array[j+0] = data_array[j+1];
      data_array[j+1] = tmpchar;
  }
}

static int ConvertBinary (void * data, int size, int count)
{
    if (size == sizeof (short))
	SwapNybbleByteOrder ((char *) data, count);
    else if (size != sizeof (char))
	SwapWordByteOrder ((char *) data, count);
	
    return 0;
}
float * SegyReader::getSegyDataAndChar(SegyDim3s &dim,char *filename_buf)
{
	int maxvals;
	
    // Declare a dynamically allocated vector of IDs.
	vector<ShotPoint> spArray;

    // Iterator is used to loop through the vector.	
    vector<ShotPoint>::iterator spIterator;

	/***********************/
	/* Function's Body     */
	/***********************/

	// Use AVS/Express function to map enviroment variables

	FILE *infile;
	struct segy_ebcdic_hdr ebcdic_hdr;
	struct segy_reel_hdr reel_hdr;
	struct segy_trace_hdr trace_hdr;
	float *trace_vals;
	float *vals ,*vals1;

	int trace_count = 0;
	int index;
	float minx, miny, maxx, maxy;
	float intervall=0.0;

	minx = miny = maxx = maxy = 0.0;

	if (!strlen(filename_buf)) return NULL;

	if ((infile = fopen(filename_buf, "rb")) == NULL) {
		fprintf(stdout,"ReadSegy: Cannot open segy file - %s\n",filename_buf);
		return NULL;
	}

	fread(&ebcdic_hdr, SEGY_EBCDIC_HDR_SIZE, 1, infile);
	fread(&reel_hdr, SEGY_REEL_HDR_SIZE, 1, infile);
	ConvertBinary(&reel_hdr.job_id_number, sizeof(int32_T), 1);
	ConvertBinary(&reel_hdr.line_number, sizeof(int32_T), 1);
	ConvertBinary(&reel_hdr.reel_number, sizeof(int32_T), 1);
	ConvertBinary(&reel_hdr.traces_per_record, sizeof(int16_t), 1);
	ConvertBinary(&reel_hdr.aux_traces_per_record, sizeof(int16_t), 1);
	ConvertBinary(&reel_hdr.sample_data_interval_ms, sizeof(int16_t), 1);
	ConvertBinary(&reel_hdr.original_data_interval_ms, sizeof(int16_t), 1);
	ConvertBinary(&reel_hdr.samples_per_trace, sizeof(int16_t), 1);
	ConvertBinary(&reel_hdr.original_samples_per_trace, sizeof(int16_t), 1);
	ConvertBinary(&reel_hdr.data_sample_format_code, sizeof(int16_t), 1);

	int XDIM = 1;
	int YDIM = reel_hdr.traces_per_record;
	int ZDIM = reel_hdr.samples_per_trace;

	// we might have to juggle the coordiates a bit 

	if (ZDIM<=0) 
	{	
		fprintf(stdout,"ReadSegy: Cannot read segy file - %s\n",filename_buf);
		return false;
	};
	
	trace_vals = (float*) malloc (ZDIM*sizeof(float));

	maxvals = 200;
	vals = (float*) malloc (maxvals*ZDIM*sizeof(float));
	

	int first=1,YIndex=0,YLong1=0,YLong2=0,YFlag=0;
	while (!feof(infile)) {
		if (!fread(&trace_hdr, SEGY_TRACE_HDR_SIZE, 1, infile)) break;

		ConvertBinary(&trace_hdr.trace_sequence_number_within_line, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.trace_sequence_number_within_reel, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.original_field_record_number, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.trace_sequence_number_within_original_field_record, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.energy_source_point_number, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.cdp_ensemble_number, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.trace_sequence_number_within_cdp_ensemble, sizeof(int32_T), 1);
		//_4_*_2_=_8_bytes
		ConvertBinary(&trace_hdr.trace_identification_code, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.number_of_vertically_summed_traces_yielding_this_trace, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.number_of_horizontally_stacked_traced_yielding_this_trace, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.data_use, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.distance_from_source_point_to_receiver_group, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.receiver_group_elevation, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.surface_elevation_at_source, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.source_depth_below_surface, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.datum_elevation_at_receiver_group, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.datum_elevation_at_source, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.water_depth_at_source, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.water_depth_at_receiver_group, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr. scalar_for_elevations_and_depths, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr. scalar_for_coordinates, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.x_source_coordinate, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.y_source_coordinate, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.x_receiver_group_coordinate, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.y_receiver_group_coordinate, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.coordinate_units, sizeof(int16_t), 1);

		ConvertBinary(&trace_hdr.weathering_velocity, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.subweathering_velocity, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.uphole_time_at_source, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.uphole_time_at_group, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.source_static_correction, sizeof(int16_t), 1);

		ConvertBinary(&trace_hdr.group_static_correction, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.total_static_applied, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.lag_time_a, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.lag_time_b, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.delay_according_time, sizeof(int16_t), 1);

		ConvertBinary(&trace_hdr.brute_time_start, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.mute_time_end, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.samples_in_this_trace, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.sample_intervall, sizeof(int16_t), 1);
		ConvertBinary(&trace_hdr.gain_type_instruments, sizeof(int16_t), 1);


		if (first) {
			first=0;
			minx = trace_hdr.x_source_coordinate;
			miny = trace_hdr.y_source_coordinate;
			maxx = trace_hdr.x_source_coordinate;
			maxy = trace_hdr.y_source_coordinate;
		}
		if (minx > trace_hdr.x_source_coordinate) {
			minx = trace_hdr.x_source_coordinate;
		}
		if (miny > trace_hdr.y_source_coordinate) {
			miny = trace_hdr.y_source_coordinate;
		}
		if (maxx < trace_hdr.x_source_coordinate) {
			maxx = trace_hdr.x_source_coordinate;
		}
		if (maxy < trace_hdr.y_source_coordinate) {
			maxy = trace_hdr.y_source_coordinate;
		}
		
		intervall=trace_hdr.sample_intervall;

		if (!fread(trace_vals, sizeof(float), ZDIM, infile))
			break;

		ConvertBinary(trace_vals, sizeof(float), ZDIM);

		if (trace_count>=maxvals)
		{
			vals1 = (float*) SEGY_realloc (vals, (maxvals+200)*ZDIM*sizeof(float), (maxvals+200)*ZDIM*sizeof(float)); 
			if (!vals1) {
				// not enough memory
				fprintf(stderr,"ReadSegy: Not Enough Memory to read %s\n",filename_buf);fflush(stderr);			
				if (vals)
				{
					free(vals);
					vals=NULL;
				}
				if (trace_vals)
				{
					free(trace_vals);
					trace_vals=NULL;
				}
				return NULL;
			};
			vals=vals1;
			maxvals+=200;
		};
		memcpy(&vals[trace_count*ZDIM],trace_vals,ZDIM*sizeof(float));
		spArray.push_back(ShotPoint(trace_hdr.x_source_coordinate, trace_hdr.y_source_coordinate,trace_count));
		if(YIndex==0)
		{
            YLong1=trace_hdr.original_field_record_number;
		}
		if(YFlag==0)
		{
			YLong2=YLong1;
			YLong1=trace_hdr.original_field_record_number;
			if(YLong1==YLong2)
			{
                YIndex=YIndex+1;
			}
			else
			{
				YFlag=1;
			}
		}

		++trace_count;
	};
    YDIM=YIndex;
	XDIM=trace_count/YDIM;
	fclose(infile);	
	dim.XDIM=ZDIM;
	dim.YDIM=YDIM;
	dim.ZDIM=XDIM;
	
	if (trace_vals)
	{
		free(trace_vals);
		trace_vals=NULL;
	}
	
	return vals;
}