SEGY IO (IBM&PC)
本文档将介绍SEGY的读取与写入过程，其中包括IBM与PC两种数据格式的转换。
程序将调用IEEE2IBM.cpp文件完成IBM与PC格式的互相转换。
新建头文件ReadSeismic.h与C++文件ReadSeismic.cpp，以及主程序main.cpp。
1 头文件ReadSeismic.h的编写及其规范
1.1 程序描述、调用、声明、定义
/**********************************************************************
* Copyright(C) 2018,Company All Rights Reserved (1)版权说明
*
* @file : ReadSeismic.cpp (2) 文件名
*
* @brief : 实现地震数据的读、写操作 (3) 该文件主要功能简介
*
* @version : 1.0 (4) 版本信息
*
* @author : Fan XinRan (5) 创建作者
*
* @date : 2022/2/8 星期二 (6) 创建时间
*
* Others : (7) 备注、改动信息等
**********************************************************************/
//调用需要的C头文件
#include<stdio.h> //C Language header file
#include<stdlib.h>
#include<string.h>
#include<math.h>
//调用需要的C++头文件
#include<iostream> // C++ header file
#include<vector>
#include<algorithm>
//调用非标准库
#include"alloc.h" // 用于创建多维数组
#include"segy.h" // 包含segy与bhed结构体，用于提取卷头和道头中采集、存储的信息
// 定义全局变量及命名空间
#define PI 3.141592654 //Constant Number Definition
#define EPS 0.0000001
using namespace std;
1.2 声明函数
unsigned short exchangeLowHigh16(unsigned short Data_temp);//16位高低位转换函数 short占2字节，2*8
unsigned int exchangeLowHigh32(unsigned int Data_temp); //32位高低位转换函数 4*8
float ibm2pc(unsigned int Data_temp); //IBM转PC数据
unsigned int pc2ibm(float input); //PC转IBM数据
float ieee2pc(unsigned int Data_temp); //IEEE转为PC
void trace_ibm2pc(float *data_output, int *data_input, int nt); //地震道数据由IBM转换为PC格式
void trace_pc2ibm(float *data_input, int *data_output, int nt); //地震道数据由PC转换为IBM格式
bool copySeismicDataIBM(const char *filenameInput, const char *filenameOutput); //Copy seismic data from Inputfile to Outputfile
完整代码
/**********************************************************************
* Copyright(C) 2018,Company All Rights Reserved
*
* @file : ReadSeismic.cpp
*
* @brief : 实现地震数据的读、写操作
*
* @version : 1.0
*
* @author : Fan XinRan
*
* @date : 2022/2/8 星期二
*
* Others :
**********************************************************************/
//(1)调用需要的C头文件
#include<stdio.h> // C Language header file
#include<stdlib.h>
#include<string.h>
#include<math.h>
//(2)调用需要的C++头文件
#include<iostream> // C++ header file
#include<vector>
#include<algorithm>
//(3)调用需要的非标准库头文件
#include"alloc.h" // project header file
#include"segy.h"
#include
//(4)定义全局常量
#define PI 3.141592654 // Constant Number Definition
#define EPS 0.0000001
//(5)声明命名空间
using namespace std;
//(6)声明函数名、输入、输出及其类型
unsigned short exchangeLowHigh16(unsigned short Data_temp);//16位高低位转换函数 short占2字节，2*8
unsigned int exchangeLowHigh32(unsigned int Data_temp); //32位高低位转换函数 4*8
float ibm2pc(unsigned int Data_temp); //IBM转PC数据
unsigned int pc2ibm(float input); //PC转IBM数据
float ieee2pc(unsigned int Data_temp); //IEEE转为PC
void trace_ibm2pc(float *data_output, int *data_input, int nt); //地震道数据由IBM转换为PC格式
void trace_pc2ibm(float *data_input, int *data_output, int nt); //地震道数据由PC转换为IBM格式
bool copySeismicDataIBM(const char *filenameInput, const char *filenameOutput); //Copy seismic data from Inputfile to Outputfile
2 格式转换工具——IEEE2IBM.cpp
#include"ReadSeismic.h"
unsigned short exchangeLowHigh16(unsigned short Data_temp)//16位高低位转换函数 short占2字节，2*8
{
unsigned short result;
unsigned short temp1=Data_temp>>8;//右移位
unsigned short temp2=Data_temp<<8;
result=temp1^temp2;
return result;//返回高低位转换结果
}
unsigned int exchangeLowHigh32(unsigned int Data_temp)//32位高低位转换函数 4*8
{
signed int *pp=(signed int *)&Data_temp;//读取Data_temp地址
char *b;
b=(char *)pp;
unsigned char temp = b[0]; //8位为单位
b[0] = b[3]; b[3] = temp;
temp = b[1]; b[1] = b[2]; b[2] = temp;
unsigned int result=*pp;
return result;//返回高低位转换结果
}
float ibm2pc(unsigned int Data_temp)//IBM转PC数据
{
signed int *pp=(signed int *)&Data_temp;//无符号整型方式读取Data_temp地址
char *b;
b=(char *)pp;
unsigned char temp = b[0]; //8位为单位
b[0] = b[3]; b[3] = temp;
temp = b[1]; b[1] = b[2]; b[2] = temp;
double sign = (double)( Data_temp >>31) ; //取符号位
double exp =(double )( ( Data_temp &0x7f000000) >>24) ;
exp = exp - 64 ;
double frac1 = (double )( Data_temp &0x00ffffff ) ;
double frac = frac1/ (pow(2.0,24) );
float result;
result = (float)(( 1-2*sign)*( pow( 16 ,exp) ) *frac);
return result; //返回PC格式数据
}
unsigned int pc2ibm(float input)//PC转IBM数据
{
if(input==0)
{
unsigned int result;
result=0;
return result;
}
else
{
int sign;//符号位
sign = ( input<0?1:0 ) ;
int exp;//
float input1 ; // attention : cannot use long input1;
input = input * pow(-1.0, sign);//
exp=0;
input1 = input;
if (input >0 ) //
{
if( (int)input>0)// sgyReadWrite();
{
exp++;
while ((int) input1/16 > 0)
{
exp++;
input1= input1/16;
}
}
else
{
while ( (int)input1*16 ==0)
{
exp--;
input1=input1*16;
}
exp++;//
}
}
int e;
e = ( exp + 64 ) ;
double fm = input * pow(16.0,-exp);////////////////
int fmant=(int) ( fm * pow(2.0,24) ) ;//
unsigned int result ;
result = ( sign<<31) | ( e <<24 ) | fmant ;
char *b;
//unsigned int *b;
b=(char*)&result;
unsigned char temp = b[0];//8位为单位
b[0] = b[3]; b[3] = temp;
temp = b[1]; b[1] = b[2]; b[2] = temp;
return result;//返回IBM格式数据
}
}
float ieee2pc(unsigned int Data_temp) //IEEE转为PC
{
signed int *pp=(signed int *)&Data_temp;//无符号整型方式读取Data_temp地址
char *b;
b=(char *)pp;
unsigned char temp = b[0]; //8位为单位
b[0] = b[3]; b[3] = temp;
temp = b[1]; b[1] = b[2]; b[2] = temp;
double sign = (double)( Data_temp >>31) ; //取符号位
int e;//
e=(int)((Data_temp & 0x7f800000 )>>23)-127;
unsigned int x ; //
x = (unsigned int)((Data_temp & 0x007fffff) -sign); //- sign : for value < 0
float x0 ;
x0 = x* pow(2.0,-23);
float result;
if ( x0 ==0 && e + 127 ==0 ) //
result = 0;
else
result = pow(-1.0,sign)*(1+x0)*pow(2.0,e);
return result;//返回PC格式数据
}
void trace_ibm2pc(float *data_output, int *data_input, int nt) {
for (unsigned int it = 0; it < nt; it++) {
data_output[it] = ibm2pc(data_input[it]);
}
}
void trace_pc2ibm(float *data_input, int *data_output, int nt) {
for (unsigned int it = 0; it < nt; it++) {
data_output[it] = pc2ibm(data_input[it]);
}
}
3 C++文件ReadSeismic.cpp的编写及其规范
3.1 必要的说明
/*************************************************************************************************************
Function: copySeismicDataIBM (1)函数名
Description: copy segy file from input data to output data (2)简要描述其功能
Input:
const char *filenameInput [in] input filename (.segy) (3)输入变量及输入文件类型
Output:
const char *filenameOutput[out] output filename (.segy) (4)输出变量及输出文件类型
Return:
bool true program success
bool false program failed (5)返回值及其说明
Author: Fan XinRan (6)创建作者
Date : 2022/2/8 (7)创建时间
Others: (8)备注、改动信息等
*************************************************************************************************************/
3.2 定义读、写函数
#include "ReadSeismic.h"
bool copySeismicDataIBM(const char *filenameInput, const char *filenameOutput){
//实现代码
...
}
（1）定义待使用的结构体变量、数值型变量
bhed 与 segy均为定义在segy.h中的结构体(structure)，分别包含了二进制卷头信息与道头信息，使用成员访问运算符(.)获取其内容；
使用unsigned int 声明整型变量，使用long long 或__int64声明SEGY文件字节数、地震道字节数，防止数据量超出范围，并且尽可能初始化变量。
bhed fileheader; // file header 卷头
segy traceheader; // trace header 道头
unsigned int nt=0; // number of sample 采样点数
unsigned int sizefileheader=sizeof(fileheader); // size of fileheader;
unsigned int sizetraceheader=sizeof(traceheader); // size of traceheader;
unsigned int ncdp = 0; // number of cdp 道数
long long size_file = 0; //size of input file
long long size_trace = 0; //size of per-trace
（2）新建指针变量
在读、写地震道数据这一任务中，需要用到输入指针、输出指针、地震道数据指针、道头指针以及一个临时指针变量，共五个指针变量。
FILE *fpinput = NULL; // input file pointer
FILE *fpoutput = NULL; // output file pointer
float *dataInput = NULL; // input data pointer
segy *traceheaderArray = NULL; // traceheader pointer
int* temp = NULL; // temp pointer
之前的任务中fread(dataInput[itrace],nt * sizeof(float),1,fpinput)，
整个读、写流程为：fpinput-->读取nt个数据-->dataInput[itrace]-->写入nt个数据-->fpoutput
本次任务中需要对数据进行变换，流程变为：fpinput-->读取nt个数据-->temp-->IBM/PC转换-->dataInput[itrace]-->写入nt个数据-->fpoutput。
（3）打开输入、输出文件指针
fpinput = fopen(filenameInput, "rb"); //open input file pointer
fpoutput = fopen(filenameOutput,"wb"); //open output file pointer
//读写操作
...
//fopen()与fclose()成对出现，在对文件的操作完成后切记关闭文件
fclose(fpinput); //close input file pointer
fclose(fpoutput); //close output file pointer
（4）判断文件打开情况
if(fpinput==NULL){ //如果文件指针为NULL
printf("Cannot open %s file\n", filenameInput); //打印“文件打开失败”
return false; //结束程序
}
if(fpoutput==NULL){
printf("Cannot open %s file\n", filenameOutput);
return false;
}
（5）读取/计算卷、道信息
fread(&fileheader,sizefileheader,1,fpinput); // 从输入流(fpinput)中读取卷头信息到指定地址---->fileheader
nt = exchangeLowHigh16(fileheader.hns) // 高低位转换
_fseeki64(fpinput,0,SEEK_END); // 从文件末尾偏移这个结构体0个长度给文件指针fpinput，即fpinput此时指向文件尾
size_file = _ftelli64(fpinput); // 返回当前文件位置，即文件总字节数
size_trace = nt*sizeof(float)+sizetraceheader; // 每一道的字节数 = 采样点字节数+道头字节数
ncdp = (size_file - (long long)sizefileheader)/size_trace; // 道数 = （文件总字节数 - 卷头字节数）/每一道的字节数
_fseeki64(fpinput,sizefileheader,SEEK_SET); // 从文件开头偏移sizefileheader（卷头字节数）个长度给指针fpinput，即fpinput此时指向第一道的开始
fwrite(&fileheader, sizefileheader, 1, fpoutput); // 先写入卷头
fread()从给定流读取数据到指针所指向的数组中；fwrite(*ptr, size, nmemb,*stream)参数与fread()相同，把ptr所指向的数组中的数据写入到给定流stream中；_fseeki64的用法与fseek相同，表示从文件指定位置偏移一定字节数；前者具有更高的兼容性；_ftelli64与ftell同理，返回给定流的当前文件位置；exchangeLowHigh16()完成short型的高低位转换，int型的高低位转换使用exchangeLowHigh64()。
（6）遍历每一条地震道，读、写数据
dataInput=alloc2float(nt,ncdp); // 分配nt*ncdp（采样点数×道数）所需的内存空间，用来存放二维地震道数据
//其中，alloc2float是卸载alloc.h中的函数，创建一个float型的二维数组
//dataInput为二级指针，可简记为dataInput指向每一行的开头
memset(dataInput[0], 0, nt*ncdp * sizeof(float)); // 从第一行的开头开始，将内存块中nt*ncdp个字符赋值为0
// dataInput指向每行开头，而dataInput[0]则为整个二维数组的起始位置
// 在内存的动态存储区中分配ncdp个长度为sizetraceheader的连续空间
traceheaderArray = (segy*)calloc(ncdp,sizetraceheader);
temp = (int*)calloc(nt,sizeof(int));
//逐道读取道头与地震道数据
for(int itrace = 0; itrace < ncdp; itrace++){
fread(&traceheaderArray[itrace],sizetraceheader,1,fpinput); // &traceheaderArray[itrace]为第itrace道的地址，读取该道头信息
fread(temp,nt * sizeof(float),1,fpinput); // 读取nt个采样点的信息并将结果指向temp
// 使用trace_ibm2pc将temp位置后nt个采样点的数据，进行IBM-->PC的转换，结果指向每一道开头(dataInput[itrace])
trace_ibm2pc(dataInput[itrace], temp, nt);
}//end for(int itrace = 0; itrace < ncdp; itrace++)
//逐道写入道头与地震道数据
for (int itrace = 0; itrace < ncdp; itrace++) {
fwrite(&traceheaderArray[itrace], sizetraceheader, 1, fpoutput); // 写入该道头信息
// 使用trace_pc2ibm将temp位置后nt个采样点的数据，进行PC-->IBM的转换，结果指向每一道开头(dataInput[itrace])
trace_pc2ibm(dataInput[itrace],temp,nt);
fwrite(temp, nt * sizeof(int), 1, fpoutput); // 以IBM的格式存回fpoutput
}//end for(int itrace = 0; itrace < ncdp; itrace++)
// 在每个循环末尾的"}"添加备注，便于寻找和区分
//在写操作完成后释放内存
free(temp); // free temp pointer
free(traceheaderArray); // free traceheader pointer
free2float(dataInput); // free data input pointer
calloc(num,size)：在内存的动态存储区中分配num个长度为size的连续空间，函数返回一个指向分配起始地址的指针；如果分配不成功，返回NULL；malloc(size)：功能与calloc()相似，不同之处是malloc()不会将内存值初始化为0，而calloc()会将新申请的内存填充0。
完整代码
/*****************************************************************************
Function: CopySeismicData
Description: copy segy file from input data to output data
Input:
const char *filenameInput [in] input filename (.segy)
Output:
const char *filenameOutput[out] output filename (.segy)
Return:
bool true program success
bool false program failed
Author: Fan XinRan
Date : 2022/2/8
Others:
*****************************************************************************/
#include "ReadSeismic.h"
bool copySeismicDataIBM(const char *filenameInput, const char *filenameOutput){
segy traceheader; // trace header
bhed fileheader; // file header
unsigned int nt = 0; // number of sample
unsigned int sizetraceheader = sizeof(traceheader); // size of traceheader;
unsigned int sizefileheader = sizeof(fileheader); // size of fileheader;
unsigned int ncdp = 0; // number of cdp
long long size_file = 0; //size of input file
long long size_trace = 0; //size of per-trace
FILE *fpinput = NULL; // input file pointer
FILE *fpoutput = NULL; //output file pointer
float **dataInput = NULL; //input data pointer
segy *traceheaderArray = NULL; //
int* temp = NULL;
fpinput = fopen(filenameInput, "rb"); //open input file pointer
fpoutput = fopen(filenameOutput, "wb"); //open output file pointer
if (fpinput == NULL) {
printf("Cannot open %s file\n", filenameInput);
return false;
}
if (fpoutput == NULL) {
printf("Cannot open %s file\n", filenameOutput);
return false;
}
fread(&fileheader, sizefileheader, 1, fpinput);
nt = fileheader.hns;
nt = exchangeLowHigh16(fileheader.hns);
_fseeki64(fpinput, 0, SEEK_END);
size_file = _ftelli64(fpinput);
size_trace = nt * sizeof(float) + sizetraceheader;
ncdp = (size_file - (long long)sizefileheader) / size_trace;
_fseeki64(fpinput, sizefileheader, SEEK_SET);
dataInput = alloc2float(nt, ncdp);
memset(dataInput[0], 0, nt*ncdp * sizeof(float));
traceheaderArray = (segy*)calloc(ncdp, sizetraceheader);
temp = (int*)calloc(nt,sizeof(int));
fwrite(&fileheader,sizefileheader,1, fpoutput);
for (int itrace = 0; itrace < ncdp; itrace++) {
fread(&traceheaderArray[itrace], sizetraceheader, 1, fpinput);
fread(temp, nt * sizeof(int), 1, fpinput);
trace_ibm2pc(dataInput[itrace], temp, nt);
}//end for(int itrace = 0; itrace < ncdp; itrace++)
for (int itrace = 0; itrace < ncdp; itrace++) {
fwrite(&traceheaderArray[itrace], sizetraceheader, 1, fpoutput);
trace_pc2ibm(dataInput[itrace],temp,nt);
fwrite(temp, nt * sizeof(int), 1, fpoutput);
}//end for(int itrace = 0; itrace < ncdp; itrace++)
free(temp);
free(traceheaderArray);
free2float(dataInput); // free data input pointer
fclose(fpoutput); //close output file pointer
fclose(fpinput); //close input file pointer
return true;
}
4 主函数main.cpp及运行结果
#include"ReadSeismic.h"
void main(){
copySeismicDataIBM("Azi6-Ang35-BZ19-6-1.segy","Outputibm.segy");
}
运行主函数后，程序会读入Azi6-Ang35-BZ19-6-1.segy，这是一个IBM格式的数据。再写入到Outputibm.segy，从而完成对SEGY文件的复制。