C/C++ 地震資料SEGY檔案 IBM格式與PC格式的轉換
阿新 • • 發佈：2022-03-16
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檔案的複製。