427 lines
13 KiB
C++
427 lines
13 KiB
C++
#include "stdafx.h"
|
||
#include "TBPImageAlgCls.h"
|
||
#include <QDateTime>
|
||
#include <QDebug>
|
||
#include <QString>
|
||
#include <cmath>
|
||
#include <QProgressDialog>
|
||
|
||
|
||
|
||
void TBPImageProcess(QString echofile, QString outImageFolder, QString imagePlanePath,long num_thread)
|
||
{
|
||
std::shared_ptr<EchoL0Dataset> echoL0ds(new EchoL0Dataset);
|
||
echoL0ds->Open(echofile);
|
||
|
||
std::shared_ptr< SARSimulationImageL1Dataset> imagL1(new SARSimulationImageL1Dataset);
|
||
imagL1->setCenterAngle(echoL0ds->getCenterAngle());
|
||
imagL1->setCenterFreq(echoL0ds->getCenterFreq());
|
||
imagL1->setNearRange(echoL0ds->getNearRange());
|
||
imagL1->setRefRange((echoL0ds->getNearRange() + echoL0ds->getFarRange()) / 2);
|
||
imagL1->setFarRange(echoL0ds->getFarRange());
|
||
imagL1->setFs(echoL0ds->getFs());
|
||
imagL1->setLookSide(echoL0ds->getLookSide());
|
||
imagL1->OpenOrNew(outImageFolder, echoL0ds->getSimulationTaskName(), echoL0ds->getPluseCount(), echoL0ds->getPlusePoints());
|
||
|
||
|
||
TBPImageAlgCls TBPimag;
|
||
TBPimag.setEchoL0(echoL0ds);
|
||
TBPimag.setImageL1(imagL1);
|
||
TBPimag.setImagePlanePath(imagePlanePath);
|
||
|
||
TBPimag.Process(num_thread);
|
||
}
|
||
|
||
void TBPImageAlgCls::setImagePlanePath(QString INimagePlanePath)
|
||
{
|
||
this->imagePlanePath = INimagePlanePath;
|
||
}
|
||
|
||
QString TBPImageAlgCls::getImagePlanePath()
|
||
{
|
||
return this->imagePlanePath;
|
||
}
|
||
|
||
void TBPImageAlgCls::setEchoL0(std::shared_ptr<EchoL0Dataset> inL0ds)
|
||
{
|
||
this->L0ds = inL0ds;
|
||
}
|
||
|
||
void TBPImageAlgCls::setImageL1(std::shared_ptr<SARSimulationImageL1Dataset> inL1ds)
|
||
{
|
||
this->L1ds = inL1ds;
|
||
}
|
||
|
||
std::shared_ptr<EchoL0Dataset> TBPImageAlgCls::getEchoL1()
|
||
{
|
||
return this->L0ds;
|
||
}
|
||
|
||
std::shared_ptr<SARSimulationImageL1Dataset> TBPImageAlgCls::getImageL0()
|
||
{
|
||
return this->L1ds;
|
||
}
|
||
|
||
ErrorCode TBPImageAlgCls::Process(long num_thread)
|
||
{
|
||
if (GPURUN) {
|
||
return this->ProcessGPU();
|
||
}
|
||
else {
|
||
return this->ProcessCPU(num_thread);
|
||
}
|
||
}
|
||
|
||
|
||
ErrorCode TBPImageAlgCls::ProcessCPU(long num_thread)
|
||
{
|
||
omp_set_num_threads(num_thread);
|
||
// ³£ÓòÎÊý
|
||
long rowCount = this->L1ds->getrowCount();
|
||
long colCount = this->L1ds->getcolCount();
|
||
long pixelCount = rowCount * colCount;
|
||
long PRFCount = this->L0ds->getPluseCount();
|
||
long PlusePoints = this->L0ds->getPlusePoints();
|
||
|
||
|
||
double Rnear = this->L1ds->getNearRange();
|
||
double Rfar = this->L1ds->getFarRange();
|
||
double fs = this->L1ds->getFs();
|
||
double dx = LIGHTSPEED / 2 / fs;
|
||
double factorj = this->L1ds->getCenterFreq() * 4 * M_PI / LIGHTSPEED * 1e9;
|
||
|
||
|
||
Eigen::MatrixXcd echo = Eigen::MatrixXcd::Zero(PRFCount, PlusePoints);
|
||
{
|
||
std::shared_ptr<std::complex<double>> echodata = this->L0ds->getEchoArr();
|
||
for (long i = 0; i < PRFCount; i++) {
|
||
for (long j = 0; j < PlusePoints; j++) {
|
||
echo(i, j) = echodata.get()[i * PlusePoints + j];
|
||
}
|
||
}
|
||
echodata.reset();
|
||
}
|
||
|
||
Eigen::MatrixXd pixelX = Eigen::MatrixXd::Zero(rowCount, colCount);
|
||
Eigen::MatrixXd pixelY = Eigen::MatrixXd::Zero(rowCount, colCount);
|
||
Eigen::MatrixXd pixelZ = Eigen::MatrixXd::Zero(rowCount, colCount);
|
||
|
||
Eigen::MatrixXd Pxs = Eigen::MatrixXd::Zero(this->L0ds->getPluseCount(), 1);
|
||
Eigen::MatrixXd Pys = Eigen::MatrixXd::Zero(this->L0ds->getPluseCount(), 1);
|
||
Eigen::MatrixXd Pzs = Eigen::MatrixXd::Zero(this->L0ds->getPluseCount(), 1);
|
||
|
||
// ͼÏñÍø¸ñ×ø±ê
|
||
{
|
||
std::shared_ptr<double> antpos = this->L0ds->getAntPos();
|
||
double time = 0;
|
||
double Px = 0;
|
||
double Py = 0;
|
||
double Pz = 0;
|
||
double Vx = 0;
|
||
double Vy = 0;
|
||
double Vz = 0;
|
||
double AntDirectX = 0;
|
||
double AntDirectY = 0;
|
||
double AntDirectZ = 0;
|
||
double AVx = 0;
|
||
double AVy = 0;
|
||
double AVz = 0;
|
||
|
||
double R = 0;
|
||
double NormAnt = 0;
|
||
|
||
for (long i = 0; i < rowCount; i++) {
|
||
time = antpos.get()[i *19 + 0];
|
||
Px = antpos.get()[i *19 + 1];
|
||
Py = antpos.get()[i *19 + 2];
|
||
Pz = antpos.get()[i *19 + 3];
|
||
Vx = antpos.get()[i *19 + 4];
|
||
Vy = antpos.get()[i *19 + 5];
|
||
Vz = antpos.get()[i *19 + 6];
|
||
AntDirectX = antpos.get()[i *19 + 13]; // Zero doppler
|
||
AntDirectY = antpos.get()[i *19 + 14];
|
||
AntDirectZ = antpos.get()[i *19 + 15];
|
||
AVx = antpos.get()[i *19 + 10];
|
||
AVy = antpos.get()[i *19 + 11];
|
||
AVz = antpos.get()[i *19 + 12];
|
||
|
||
NormAnt = std::sqrt(AntDirectX * AntDirectX + AntDirectY * AntDirectY + AntDirectZ * AntDirectZ);
|
||
AntDirectX = AntDirectX / NormAnt;
|
||
AntDirectY = AntDirectY / NormAnt;
|
||
AntDirectZ = AntDirectZ / NormAnt;// ¹éÒ»»¯
|
||
|
||
antpos.get()[i *19 + 13] = AntDirectX;
|
||
antpos.get()[i *19 + 14] = AntDirectY;
|
||
antpos.get()[i *19 + 15] = AntDirectZ;
|
||
Pxs(i, 0) = Px;
|
||
Pys(i, 0) = Py;
|
||
Pzs(i, 0) = Pz;
|
||
|
||
for (long j = 0; j < colCount; j++) {
|
||
R = j * dx + Rnear;
|
||
pixelX(i, j) = Px + AntDirectX * R;
|
||
pixelY(i, j) = Py + AntDirectY * R;
|
||
pixelZ(i, j) = Pz + AntDirectZ * R;
|
||
}
|
||
}
|
||
|
||
this->L1ds->saveAntPos(antpos);
|
||
antpos.reset();
|
||
}
|
||
|
||
// BP³ÉÏñ
|
||
long BlockLine = Memory1MB * 10 / 16 / rowCount;
|
||
if (rowCount / BlockLine / num_thread < 3) {
|
||
BlockLine = rowCount / num_thread / 3;
|
||
}
|
||
BlockLine = BlockLine > 10 ? BlockLine : 10;
|
||
|
||
std::shared_ptr<std::complex<double>> imagarr = this->L1ds->getImageRaster();
|
||
{
|
||
for (long i = 0; i < pixelCount; i++) {
|
||
imagarr.get()[i] = imagarr.get()[i];
|
||
}
|
||
}
|
||
omp_lock_t lock; // ¶¨ÒåËø
|
||
omp_init_lock(&lock); // ³õʼ»¯Ëø
|
||
|
||
long writeImageCount = 0;
|
||
qDebug() << "block line:\t" << BlockLine;
|
||
|
||
|
||
long startLine = 0;
|
||
long processValue = 0;
|
||
long processNumber = 0;
|
||
QProgressDialog progressDialog(u8"RTPC»Ø²¨Éú³ÉÖÐ", u8"ÖÕÖ¹", 0, rowCount);
|
||
progressDialog.setWindowTitle(u8"RTPC»Ø²¨Éú³ÉÖÐ");
|
||
progressDialog.setWindowModality(Qt::WindowModal);
|
||
progressDialog.setAutoClose(true);
|
||
progressDialog.setValue(0);
|
||
progressDialog.setMaximum(rowCount);
|
||
progressDialog.setMinimum(0);
|
||
progressDialog.show();
|
||
|
||
#pragma omp parallel for
|
||
for (startLine = 0; startLine < rowCount; startLine = startLine + BlockLine) { // ͼÏñ´óС
|
||
long stepLine = startLine + BlockLine < rowCount ? BlockLine : rowCount - startLine;
|
||
long imageRowID = startLine; //
|
||
|
||
//Eigen::MatrixXd R = Eigen::MatrixXd::Zero(rowCount, 1);
|
||
long pluseId = 0;
|
||
std::complex<double> factPhas(0, 0);
|
||
std::complex < double> sign(0, 0);
|
||
Eigen::MatrixXd R = Eigen::MatrixXd::Zero(rowCount, 1);
|
||
Eigen::MatrixXcd Rphi = Eigen::MatrixXd::Zero(rowCount, 1);
|
||
|
||
long PluseIDs = 0;
|
||
double mask = 0;
|
||
|
||
for (long i = 0; i < stepLine; i++) { // ͼÏñÐÐ
|
||
imageRowID = startLine + i;
|
||
for (long j = 0; j < colCount; j++) { //ͼÏñÁÐ
|
||
R = ((pixelX(i, j) - Pxs.array()).array().pow(2) + (pixelY(i, j) - Pys.array()).array().pow(2) + (pixelZ(i, j) - Pzs.array()).array().pow(2)).array().sqrt();
|
||
Rphi = Rphi.array() * 0;
|
||
Rphi.imag() = R.array() * factorj;
|
||
Rphi = Rphi.array().exp();
|
||
|
||
for (long prfid = 0; prfid < rowCount; prfid++) { // Âö³åÐÐ
|
||
PluseIDs = std::floor((R(prfid, 0) - Rnear) / dx);
|
||
mask = (PluseIDs < 0 || PluseIDs >= PlusePoints) ? 0 : 1;
|
||
PluseIDs = (PluseIDs < 0 || PluseIDs >= PlusePoints) ? 0 : PluseIDs;
|
||
imagarr.get()[imageRowID * colCount + j] =
|
||
imagarr.get()[imageRowID * colCount + j] +
|
||
mask * echo(prfid, PluseIDs) * Rphi(prfid, 0);// ÐźÅ* ÏàλУÕý
|
||
}
|
||
}
|
||
|
||
}
|
||
omp_set_lock(&lock); // ±£´æÎļþ
|
||
processValue = processValue + BlockLine;
|
||
this->L1ds->saveImageRaster(imagarr, 0, rowCount);
|
||
qDebug() << QDateTime::currentDateTime().toString("yyyy-MM-dd HH:mm:ss.zzz").toUtf8().constData() << "\t" << processValue * 100.0 / rowCount << "%\t" << startLine << "\t-\t" << startLine + BlockLine << "\tend\t\t";
|
||
processNumber = processNumber + BlockLine;
|
||
processNumber = processNumber < progressDialog.maximum() ? processNumber : progressDialog.maximum();
|
||
progressDialog.setValue(processNumber);
|
||
omp_unset_lock(&lock); // ½âËø
|
||
}
|
||
omp_destroy_lock(&lock); // Ïú»ÙËø
|
||
|
||
this->L1ds->saveImageRaster(imagarr, 0, rowCount);
|
||
this->L1ds->saveToXml();
|
||
progressDialog.close();
|
||
return ErrorCode::SUCCESS;
|
||
}
|
||
|
||
ErrorCode TBPImageAlgCls::ProcessGPU()
|
||
{
|
||
// ³£ÓòÎÊý
|
||
long rowCount = this->L1ds->getrowCount();
|
||
long colCount = this->L1ds->getcolCount();
|
||
long pixelCount = rowCount * colCount;
|
||
long PRFCount = this->L0ds->getPluseCount();
|
||
long PlusePoints = this->L0ds->getPlusePoints();
|
||
|
||
|
||
float Rnear = this->L1ds->getNearRange();
|
||
float Rfar = this->L1ds->getFarRange();
|
||
float fs = this->L1ds->getFs();
|
||
double dx = LIGHTSPEED / 2 / fs;
|
||
float freq = this->L1ds->getCenterFreq()*1.0*1e9;
|
||
double factorj = freq * 4 * M_PI / LIGHTSPEED ;
|
||
|
||
std::shared_ptr<float> pixelX(new float[rowCount*colCount],delArrPtr); // ͼÏñ³ÉÏñÍø¸ñ
|
||
std::shared_ptr<float> pixelY(new float[rowCount*colCount],delArrPtr);
|
||
std::shared_ptr<float> pixelZ(new float[rowCount*colCount],delArrPtr);
|
||
|
||
std::shared_ptr<float> Pxs (new float[this->L0ds->getPluseCount()]);
|
||
std::shared_ptr<float> Pys (new float[this->L0ds->getPluseCount()]);
|
||
std::shared_ptr<float> Pzs (new float[this->L0ds->getPluseCount()]);
|
||
std::shared_ptr<float> Vxs (new float[this->L0ds->getPluseCount()]);
|
||
std::shared_ptr<float> Vys (new float[this->L0ds->getPluseCount()]);
|
||
std::shared_ptr<float> Vzs (new float[this->L0ds->getPluseCount()]);
|
||
|
||
|
||
|
||
// ͼÏñÍø¸ñ×ø±ê
|
||
{
|
||
std::shared_ptr<double> antpos = this->L0ds->getAntPos();
|
||
double time = 0;
|
||
double Px = 0;
|
||
double Py = 0;
|
||
double Pz = 0;
|
||
double Vx = 0;
|
||
double Vy = 0;
|
||
double Vz = 0;
|
||
double AntDirectX = 0;
|
||
double AntDirectY = 0;
|
||
double AntDirectZ = 0;
|
||
double AVx = 0;
|
||
double AVy = 0;
|
||
double AVz = 0;
|
||
|
||
double R = 0;
|
||
double NormAnt = 0;
|
||
|
||
for (long i = 0; i < rowCount; i++) {
|
||
time = antpos.get()[i *19 + 0];
|
||
Px = antpos.get()[i *19 + 1];
|
||
Py = antpos.get()[i *19 + 2];
|
||
Pz = antpos.get()[i *19 + 3];
|
||
Vx = antpos.get()[i *19 + 4];
|
||
Vy = antpos.get()[i *19 + 5];
|
||
Vz = antpos.get()[i *19 + 6];
|
||
AntDirectX = antpos.get()[i *19 + 13];// zero doppler
|
||
AntDirectY = antpos.get()[i *19 + 14];
|
||
AntDirectZ = antpos.get()[i *19 + 15];
|
||
AVx = antpos.get()[i *19 + 10];
|
||
AVy = antpos.get()[i *19 + 11];
|
||
AVz = antpos.get()[i *19 + 12];
|
||
|
||
NormAnt = std::sqrt(AntDirectX * AntDirectX + AntDirectY * AntDirectY + AntDirectZ * AntDirectZ);
|
||
AntDirectX = AntDirectX / NormAnt;
|
||
AntDirectY = AntDirectY / NormAnt;
|
||
AntDirectZ = AntDirectZ / NormAnt;// ¹éÒ»»¯
|
||
|
||
antpos.get()[i *19 + 13] = AntDirectX;
|
||
antpos.get()[i *19 + 14] = AntDirectY;
|
||
antpos.get()[i *19 + 15] = AntDirectZ;
|
||
|
||
Pxs.get()[i] = Px;
|
||
Pys.get()[i] = Py;
|
||
Pzs.get()[i] = Pz;
|
||
Vxs.get()[i] = Vx;
|
||
Vys.get()[i] = Vy;
|
||
Vzs.get()[i] = Vz;
|
||
for (long j = 0; j < colCount; j++) {
|
||
R = j * dx + Rnear;
|
||
pixelX.get()[i*colCount+ j] = Px + AntDirectX * R;
|
||
pixelY.get()[i*colCount+ j] = Py + AntDirectY * R;
|
||
pixelZ.get()[i*colCount+ j] = Pz + AntDirectZ * R;
|
||
}
|
||
}
|
||
qDebug()<<"R: " << R;
|
||
this->L1ds->saveAntPos(antpos);
|
||
antpos.reset();
|
||
}
|
||
|
||
|
||
|
||
std::shared_ptr<std::complex<double>> Rasterarr = this->L1ds->getImageRaster();
|
||
std::shared_ptr<std::complex<double>> echodataPtr = this->L0ds->getEchoArr();
|
||
|
||
for (long i = 0; i < rowCount; i++) {
|
||
for (long j = 0; j < colCount; j++) {
|
||
Rasterarr.get()[i * colCount + j] = Rasterarr.get()[i * colCount + j] * 0.0;
|
||
}
|
||
}
|
||
|
||
|
||
TBPImageGPUAlg(Pxs, Pys, Pzs, // ÌìÏß×ø±ê
|
||
Vxs, Vys, Vzs,
|
||
pixelX, pixelY, pixelZ, // ͼÏñ×ø±ê
|
||
echodataPtr, Rasterarr,
|
||
freq, fs, Rnear, Rfar,
|
||
rowCount, colCount, this->L1ds);
|
||
|
||
qDebug() << "image writing:\t" << this->L1ds->getxmlFilePath();
|
||
|
||
return ErrorCode::SUCCESS;
|
||
}
|
||
|
||
void TBPImageAlgCls::setGPU(bool flag)
|
||
{
|
||
this->GPURUN = flag;
|
||
}
|
||
|
||
bool TBPImageAlgCls::getGPU( )
|
||
{
|
||
return this->GPURUN;
|
||
}
|
||
|
||
/// <summary>
|
||
/// TBP GPU´úÂë
|
||
/// </summary>
|
||
/// <param name="antpos_ptr">ÎÀÐǹìµÀ×ø±ê</param>
|
||
/// <param name="echoArr">»Ø²¨¾ØÕó</param>
|
||
/// <param name="img_arr">ͼÏñ¾ØÕó</param>
|
||
void TBPImageGPUAlg(std::shared_ptr<float> antPx, std::shared_ptr<float> antPy, std::shared_ptr<float> antPz, // ÌìÏß×ø±ê
|
||
std::shared_ptr<float> antVx, std::shared_ptr<float> antVy, std::shared_ptr<float> antVz,
|
||
std::shared_ptr<float> img_x, std::shared_ptr<float> img_y, std::shared_ptr<float> img_z, // ͼÏñ×ø±ê
|
||
std::shared_ptr<std::complex<double>> echoArr, std::shared_ptr<std::complex<double>> img_arr,
|
||
float freq, float fs, float Rnear, float Rfar,
|
||
long rowcount, long colcount, std::shared_ptr<SARSimulationImageL1Dataset> L1ds) {
|
||
|
||
float factorj = freq * 4 * PI / LIGHTSPEED;
|
||
|
||
|
||
qDebug() << "factorj:\t" << factorj;
|
||
qDebug() << "freq:\t" << freq;
|
||
qDebug() << "fs:\t" << fs;
|
||
qDebug() << "Rnear:\t" << Rnear;
|
||
qDebug() << "Rfar:\t" << Rfar;
|
||
|
||
qDebug() << "img_x:\t" << img_x.get()[0];
|
||
qDebug() << "img_y:\t" << img_y.get()[0];
|
||
qDebug() << "img_z:\t" << img_z.get()[0];
|
||
|
||
|
||
long blockline = Memory1MB * 1000 / sizeof(float) / colcount;
|
||
blockline = blockline < 10 ? 10 : blockline;
|
||
for (long startline = 0; startline < rowcount; startline = startline + blockline) {
|
||
long stepline = startline + blockline < rowcount ? blockline : rowcount - startline;
|
||
std::cout << startline << " \ " << rowcount << " "<< stepline << " start " << std::endl;
|
||
//TBPImageGPUBlock(antPx.get(), antPy.get(), antPz.get(), img_x.get(), img_y.get(), img_z.get(),
|
||
// echoArr, rowcount, colcount,
|
||
// img_arr,
|
||
// freq, fs, Rnear, Rfar, factorj, startline, stepline,
|
||
// stepline, colcount);
|
||
//std::cout << startline << " \ " << rowcount << " " << stepline << " end " << std::endl;
|
||
//L1ds->saveImageRaster(img_arr, 0, rowcount);
|
||
}
|
||
L1ds->saveImageRaster(img_arr, 0, rowcount);
|
||
|
||
L1ds->saveToXml();
|
||
}
|
||
|
||
|