segy.cpp

本代码提供读取Segy格式地震数据的SegyReader类（C++） 在segy.h文件中有segy卷头/道头的每一个字段的中文详细说明,方便大家查询使用. 另外, 使用QT3.x制作的G

#include "segy.h"
#include <numeric>
#include <functional>
#include <vector>
#include <iterator>
#include <math.h>
#include <ctype.h>
#include <iostream>
using namespace std;
//下面两行用于在生成索引文件时,在"炮号"与"当前炮下排列数"两个参数的前缀
//shotNo=100,arraySums=12
#define INDEX_FILE_HEAD "segyIndexFile:component shotSums shotMin shotMax\n"
#define SHOT_NO_PREFIX "shotNo="
#define ARRAY_SUMS_PREFIX "arraySums="
//打印调试信息
#define DEBUG_SEGY 1
//使用测试数据(非炮集数据)
#define DEBUG_EVAL 1

void * SegyReader::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
}

void SegyReader::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;
	}
}

void SegyReader::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;
	}
}

int SegyReader::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;
}

SegyReader::SegyReader(const char * fileNameData, const char * fileNameIndex, short componentFlag, int arraySum)
{
	//判断数据文件的存在性
	FILE *infile;
	if (!strlen(fileNameData)) {
		this->fileNameData[0] = '\0';
	}
	if ((infile = fopen(fileNameData, "rb")) == NULL) {
		fprintf(stdout,"ReadSegy: Cannot open segy file - %s\n",fileNameData);
		this->fileNameData[0] = '\0';			//文件名
	} else {
		fclose(infile);
		strcpy(this->fileNameData , fileNameData);			//文件名
	}

	//判断索引文件的存在性
	if (!strlen(fileNameIndex)) {
		this->fileNameIndex[0] = '\0';
	}
	if ((infile = fopen(fileNameIndex, "rb")) == NULL) {
		fprintf(stdout,"ReadSegy: Cannot open segy index file - %s\n",fileNameIndex);
		if(DEBUG_EVAL) { 
			strcpy(this->fileNameIndex , fileNameIndex) ;
		} else {
			this->fileNameIndex[0] = '\0';
		}
	} else {
		fclose(infile);
		strcpy(this->fileNameIndex , fileNameIndex) ;			
	}

	if (componentFlag>=0 && componentFlag<=2) {
		this->componentFlag = componentFlag; 
	} else {
		this->componentFlag = -1 ; 
	}
	this->arraySum = arraySum; 
	shotMin = -1;		
	shotMax = -1;		
	shotSum = shotMax - shotMin;		
	if (DEBUG_EVAL) {
		traceSumEachArray = 3;	
	} else {
		traceSumEachArray = 300;  
	}

	int trace_count = 0;

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

	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);	

	fseek(infile,3600, SEEK_SET);	
	if (fread(&trace_hdr, SEGY_TRACE_HDR_SIZE, 1, infile)) {
		ConvertBinary(&trace_hdr.original_field_record_number, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.samples_in_this_trace, sizeof(int16_t), 1);	
		shotMin = trace_hdr.original_field_record_number;
	}

	fseek(infile, -(240 + trace_hdr.samples_in_this_trace * 4L) , SEEK_END);
	if (fread(&trace_hdr, SEGY_TRACE_HDR_SIZE, 1, infile)) {
		ConvertBinary(&trace_hdr.original_field_record_number, sizeof(int32_T), 1);
		ConvertBinary(&trace_hdr.samples_in_this_trace, sizeof(int16_t), 1);	
		shotMax = trace_hdr.original_field_record_number;
	}

	fseek(infile,0,SEEK_END); 
	fpos_t filepos;
	fgetpos(infile, &filepos);      
	shotSum =  (filepos - 3600L ) / (240 + trace_hdr.samples_in_this_trace * 4L );  
	shotSum = shotSum / arraySum / traceSumEachArray; 	

	fclose(infile);
}

//生成索引文件
bool  SegyReader::buildIndexFile()
{
	FILE * dataFile;		//数据文件,用于读
	FILE * indexFile;		//索引文件,用于写
	long traceNo=0;			
	long curOffect=0;		
	if ((dataFile = fopen(fileNameData,"rb")) == NULL) {
		fprintf(stdout,"ReadSegy: Cannot open segy file - %s\n",fileNameData);
		return false;
	}

	if ((indexFile = fopen(fileNameIndex,"w")) == NULL ) {
		fprintf(stdout,"ReadSegy: Cannot open segy-index file - %s\n", fileNameIndex);
		return false;
	}
	fprintf(indexFile,"%d %ld %ld %ld\n",componentFlag, shotSum, shotMin, shotMax);
	curOffect = 3600L;
	fseek(dataFile,curOffect,SEEK_SET);

	while (!feof(dataFile)&& traceNo< shotSum*arraySum) {
		//读道头
		if (!fread(&trace_hdr, SEGY_TRACE_HDR_SIZE, 1, dataFile)) 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.mute_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 (DEBUG_SEGY) {
			cout<<trace_hdr.trace_sequence_number_within_line<<"\t\t\t";
			cout<<trace_hdr.trace_sequence_number_within_reel<<"\t\t";
			cout<<trace_hdr.original_field_record_number<<"\t\t";
			cout<<trace_hdr.trace_sequence_number_within_original_field_record<<"\t\t";
			cout<<trace_hdr.energy_source_point_number<<"\t\t\t";
			cout<<trace_hdr.cdp_ensemble_number<<"\t\t";
			cout<<trace_hdr.distance_from_source_point_to_receiver_group<<"\t\t\t";
			cout<<trace_hdr.samples_in_this_trace<<"\t\t\t";
			cout<<trace_hdr.sample_intervall;
			cout<<endl;
		}
		if (traceNo%12 ==0) {
			fprintf(indexFile,"%ld %ld %d\n", traceNo/12,trace_hdr.original_field_record_number,arraySum);	 //无文字提示版本,如0 125 12
		}
		long iTraceNoInArray = ( traceNo % arraySum ) +1 ;
		fprintf(indexFile,"%ld %ld\n",iTraceNoInArray,curOffect );
		
		curOffect += traceSumEachArray * (240 + trace_hdr.samples_in_this_trace * 4L) ;	
		fseek(dataFile,curOffect,SEEK_SET); 
		traceNo += 1 ;
	}

	fclose(dataFile);
	fclose(indexFile);
}

bool  SegyReader::phaseIndexFile()
{
	FILE *infile;
	if ((infile = fopen(fileNameIndex, "r")) == NULL) {
		fprintf(stdout,"ReadSegy: Cannot open segy-index file - %s\n",fileNameIndex);
		return false;
	}

	fscanf(infile, "%ld %ld %ld %ld\n",&componentFlag, &shotSum, &shotMin, &shotMax);
	if (DEBUG_SEGY) {
		printf("compont=%ld,shotSum=%ld,shotMin=%ld,shotMax=%ld\n",componentFlag, shotSum, shotMin, shotMax);
	}
    cHash = new CHash(shotSum); 	

	
	if (DEBUG_SEGY) {
		printf("====动态生成二维float指针, shotSum=%d,arraySum=%d===\n",shotSum, arraySum);
	}
	indexBuf =  new long * [shotSum];
	for (int ii=0; ii< shotSum ; ii++) {
		indexBuf[ii] = new long[arraySum];
	}

	long indexFileShotIndex, indexFileShotNo, indexFileArraySum;
	long  arrayNo, arrayOffset;
	for ( int i=0; i < shotSum;  i++) {
		fscanf(infile, "%ld %ld %ld\n",&indexFileShotIndex, &indexFileShotNo, &indexFileArraySum);
		if (DEBUG_SEGY) {
			printf("indexFileShotIndex=%ld,indexFileShotNo=%ld,indexFileArraySum=%ld\n",indexFileShotIndex, indexFileShotNo, indexFileArraySum);
		}
	    cHash->Insert(indexFileShotNo, indexFileShotIndex);

		for (int j=0; j<indexFileArraySum; j++) {
			fscanf(infile, "%ld %ld\n",&arrayNo, &arrayOffset);
			if (DEBUG_SEGY) {
				printf("arrayNo=%ld,arrayOffset=%ld\n", arrayNo, arrayOffset);
			}
			indexBuf[indexFileShotIndex][arrayNo - 1] = arrayOffset;
		}
	}
	fclose(infile);

}

void SegyReader::printIndexbuf()
{
	printf("\n================printIndexbuf=====================\n");
	for(int s=0; s< shotSum; s++) {
		printf("---shotIndex=%ld---\n", s);
		for (int a=0; a< arraySum; a++) {
			printf("array= %d, offset= %ld\n", a,indexBuf[s][a]);
		}
	}
}

float * SegyReader::getSegyData(int shotNo, int arrayNo)
{
	FILE *infile;
	float *trace_vals;		
	float *vals ,*vals1;    
	int shotIndex;
	long  indexOffset;
	cHash->QueryValue(shotNo,shotIndex);
	if (shotIndex>=0 && shotIndex<shotSum) {
		indexOffset = indexBuf[shotIndex][arrayNo]; 
	}

	if ((infile = fopen(fileNameData, "rb")) == NULL) {
		fprintf(stdout,"ReadSegy: Cannot open segy file - %s\n",fileNameData);
		return NULL;
	}
	int maxVals;
	int trace_count = 0;
	int YDIM = reel_hdr.traces_per_record;  
	int ZDIM = reel_hdr.samples_per_trace;	
	if (ZDIM<=0)
	{
		fprintf(stdout,"ReadSegy: Cannot read segy file - %s\n",fileNameData);
		return NULL;
	};


	trace_vals = (float*) malloc (ZDIM*sizeof(float));	

	maxVals = 300;                                     
	vals = (float*) malloc (maxVals*ZDIM*sizeof(float)); 

	while (!feof(infile)&& trace_count < traceSumEachArray) {
		fseek(infile, indexOffset + SEGY_TRACE_HDR_SIZE , SEEK_SET);
	
		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",fileNameData);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));
		++trace_count;
	};
	fclose(infile);

	if (trace_vals)
	{
		free(trace_vals);
		trace_vals=NULL;
	}
	return vals;
}