C/C++读取SEGY文件(三)

虚幻大学 xuhss 442℃ 0评论

Python微信订餐小程序课程视频

https://edu.csdn.net/course/detail/36074

Python实战量化交易理财系统

https://edu.csdn.net/course/detail/35475

SEGY IO (IBM&PC)

本文档将介绍SEGY的读取与写入过程,其中包括IBMPC两种数据格式的转换。

程序将调用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 //C Language header file 
#include
#include
#include

//调用需要的C++头文件
#include // C++ header file
#include
#include

//调用非标准库
#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 // C Language header file 
#include
#include
#include

//(2)调用需要的C++头文件
#include // C++ header file
#include
#include

//(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 C++文件ReadSeismic.cpp的编写及其规范

2.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)备注、改动信息等

*************************************************************************************************************/

2.2 定义读、写函数

#include "ReadSeismic.h"
bool copySeismicDataIBM(const char *filenameInput, const char *filenameOutput){

    //实现代码 
    ...

}

(1)定义待使用的结构体变量、数值型变量

bhedsegy均为定义在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); // 先写入卷头
  1. fread() 从给定流读取数据到指针所指向的数组中;
  2. fwrite(*ptr, size, nmemb,*stream) 参数与fread()相同,把ptr所指向的数组中的数据写入到给定流stream中;
  3. _fseeki64的用法与fseek相同,表示从文件指定位置偏移一定字节数;前者具有更高的兼容性;
  4. _ftelli64ftell同理,返回给定流的当前文件位置;
  5. 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
  1. calloc(num,size):在内存的动态存储区中分配num个长度为size的连续空间,函数返回一个指向分配起始地址的指针;如果分配不成功,返回NULL;
  2. 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;
}

3 主函数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文件的复制。

转载请注明:xuhss » C/C++读取SEGY文件(三)

喜欢 (8)

您必须 登录 才能发表评论!