#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <iomanip>
#include <math.h>
#include <armadillo>
#include "boost/multi_array.hpp"
typedef long int int32;
typedef boost::multi_array<float, 4> multiarray_type;
struct Position
{
std::vector<float> theta;
std::vector<float> sdepth;
std::vector<float> rdepth;
std::vector<float> rrange;
};
//函数接口
bool read_shd(const std::string &filePath, float depth, float &freq, std::vector<float> &rkm, std::vector<float> &tslice)
{
freq = 0.0f;
rkm.clear();
tslice.clear();
//打开文件
std::ifstream in(filePath, std::ios::in | std::ios::binary);
if (!in)
{
in.close();
return false;
}
//读取变量
int32 recl;
in.read((char *)&recl, sizeof(int32));
char title[81] = "\0";
in.read(title, 80 * sizeof(char));
//设置读取位置
in.seekg(4 * recl * sizeof(char), std::ios_base::beg);
char plottype[11] = "\0";
in.read(plottype, 10 * sizeof(char));
float xs = 0.0f;
in.read((char *)&xs, sizeof(float));
float ys = 0.0f;
in.read((char *)&ys, sizeof(float));
//设置读取位置,读取freq
in.seekg(2 * 4 * recl * sizeof(char), std::ios_base::beg);
in.read((char *)&freq, sizeof(float));
int32 Ntheta;
in.read((char *)&Ntheta, sizeof(int32));
int32 Nsd;
in.read((char *)&Nsd, sizeof(int32));
int32 Nrd;
in.read((char *)&Nrd, sizeof(int32));
int32 Nrr;
in.read((char *)&Nrr, sizeof(int32));
float atten = 0.0f;
in.read((char *)&atten, sizeof(float));
//设置读取位置
in.seekg(3 * 4 * recl * sizeof(char), std::ios_base::beg);
Position pos;
pos.theta.resize(Ntheta, 0.0f);
for (int i = 0; i < Ntheta; ++i)
{
in.read((char *)&pos.theta[i], sizeof(float));
}
//设置读取位置
in.seekg(4 * 4 * recl * sizeof(char), std::ios_base::beg);
pos.sdepth.resize(Nsd, 0.0f);
for (int i = 0; i < Nsd; ++i)
{
in.read((char *)&pos.sdepth[i], sizeof(float));
}
//设置读取位置,读取r.depth
in.seekg(5 * 4 * recl * sizeof(char), std::ios_base::beg);
pos.rdepth.resize(Nrd, 0.0f);
for (int i = 0; i < Nrd; ++i)
{
in.read((char *)&pos.rdepth[i], sizeof(float));
}
//设置读取位置,读取,r.range
in.seekg(6 * 4 * recl * sizeof(char), std::ios_base::beg);
pos.rrange.resize(Nrr, 0.0f);
for (int i = 0; i < Nrr; ++i)
{
in.read((char *)&pos.rrange[i], sizeof(float));
}
int Nrcvrs_per_range = 0;
multiarray_type abs_pressure;
//定义压力场维度
if (std::string(plottype) == "irregular ")
{
abs_pressure.resize(boost::extents[Ntheta][Nsd][1][Nrr]);
Nrcvrs_per_range = 1;
}
else
{
abs_pressure.resize(boost::extents[Ntheta][Nsd][Nrd][Nrr]);
Nrcvrs_per_range = Nrd;
}
//初始化压力最大值
float abs_max = 0.0f;
std::vector<float> tempArr(2 * Nrr, 0);
//读取文件中的压力数据
for (int itheta = 0; itheta < Ntheta; ++itheta)
{
for (int isd = 0; isd < Nsd; ++isd)
{
for (int ird = 0; ird < Nrcvrs_per_range; ++ird)
{
int recnum = 7 + itheta*Nsd*Nrcvrs_per_range + isd*Nrcvrs_per_range + ird;
in.seekg(recnum * 4 * recl * sizeof(char), std::ios_base::beg);
for (int i = 0; i < 2 * Nrr; ++i)
{
in.read((char *)&tempArr[i], sizeof(float));
}
for (int irr = 0; irr < Nrr; ++irr)
{
abs_pressure[itheta][isd][ird][irr] = sqrtf(tempArr[2 * irr] * tempArr[2 * irr] + tempArr[2 * irr + 1] * tempArr[2 * irr + 1]);
abs_max = abs_max > abs_pressure[itheta][isd][ird][irr] ? abs_max : abs_pressure[itheta][isd][ird][irr];
}
}
}
}
//把0值改成一个较小的浮点数,并计算dB值
for (multiarray_type::iterator it_d1 = abs_pressure.begin(); it_d1 != abs_pressure.end(); ++it_d1)
{
for (auto it_d2 = it_d1->begin(); it_d2 != it_d1->end(); ++it_d2)
{
for (auto it_d3 = it_d2->begin(); it_d3 != it_d2->end(); ++it_d3)
{
for (auto it_d4 = it_d3->begin(); it_d4 != it_d3->end(); ++it_d4)
{
*it_d4 = (*it_d4 == 0.0f) ? (abs_max / 1.0e+10) : (*it_d4);
*it_d4 = -20.0f * log10((*it_d4));
}
}
}
}
//计算rkm
std::transform(pos.rrange.begin(), pos.rrange.end(), std::back_inserter(rkm), [](float v)->float {return v / 1000.0f; });
//插值计算
if (pos.rdepth.size() == 1)
{
//只有一个元素时不插值
for (int i = 0; i < abs_pressure[0][0][0].size(); ++i)
tslice.push_back(abs_pressure[0][0][0][i]);
}
else
{
//x vector
arma::fvec x;
x.resize(pos.rdepth.size());
for (int i = 0; i < pos.rdepth.size(); ++i)
x(i) = pos.rdepth[i];
//y vectors
arma::fmat yVectors;
int rows = abs_pressure[0][0].size();
int cols = abs_pressure[0][0][0].size();
yVectors.resize(rows, cols);
for (int i = 0; i < rows; ++i)
{
for (int j = 0; j < cols; ++j)
{
yVectors(i, j) = abs_pressure[0][0][i][j];
}
}
//插值结果变量定义
arma::fvec x1(1);
x1(0) = 80.0f;
arma::fvec y1(1);
//计算插值结果
yVectors.each_col([&](const arma::fvec &col) {
arma::interp1(x, col, x1, y1);
tslice.push_back(y1(0));
});
}
in.close();
return true;
}
void output_result(const std::string & filePath, float freq, const std::vector<float> &rkm, const std::vector<float> &tslice)
{
std::ofstream out(filePath, std::ios_base::out);
if (!out)
return;
out << std::setprecision(8) << "Freq:" << freq << std::endl;
out << "rkm vector" << std::endl;
for (auto e : rkm)
out << std::setprecision(8) << e << std::endl;
out << "tslice vector" << std::endl;
for(auto e:tslice)
out << std::setprecision(8) << e << std::endl;
}
int main()
{
std::string shdFile = "D:/ReadShdFile/Isosonic3_0.shd";
float freq;
float depth;
std::vector<float> rkm;
std::vector<float> tslice;
std::cout << "Input shd file path (include file name)" << std::endl;
std::cin >> shdFile;
std::cout << "Input depth value" << std::endl;
std::cin >> depth;
if (read_shd(shdFile, depth, freq, rkm, tslice))
{
std::string resultFile = "result.txt";
std::cout << "file read successfully" << std::endl;
output_result(resultFile, freq, rkm, tslice);
}
else
{
std::cout << "file read Fail" << std::endl;
}
system("pause");
return 0;
}