C/C++讀取SEGY檔案（三）

GeoFXR發表於2022-03-16
C++

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 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）定義待使用的結構體變數、數值型變數

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); // 先寫入卷頭
  1. fread() 從給定流讀取資料到指標所指向的陣列中；
  2. fwrite(*ptr, size, nmemb,*stream) 引數與fread()相同，把ptr所指向的陣列中的資料寫入到給定流stream中；
  3. _fseeki64的用法與fseek相同，表示從檔案指定位置偏移一定位元組數；前者具有更高的相容性；
  4. _ftelli64與ftell同理，返回給定流的當前檔案位置；
  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檔案的複製。

相關文章