Compare commits
2 Commits
| Author | SHA1 | Date | |
|---|---|---|---|
|
330fedfefc | ||
|
|
99251268af |
@@ -1,78 +0,0 @@
|
||||
#include "Data3DCluster.h"
|
||||
#include "SegImgFunction.h"
|
||||
#include "Function2D.h"
|
||||
#include "Function3D.h"
|
||||
#include "maxConnectedComponent3D.cuh"
|
||||
|
||||
#include <algorithm>
|
||||
#include <float.h>
|
||||
|
||||
#include <MatlabWriter.h>
|
||||
|
||||
using namespace Aurora;
|
||||
using namespace Recon;
|
||||
|
||||
Data3DClusterResult Recon::Data3DCluster(const Aurora::Matrix& aMatrix, float aTopZ)
|
||||
{
|
||||
int clusNum = 2;
|
||||
Aurora::Matrix imgSamp = downSample(aMatrix);
|
||||
int ite = 1;
|
||||
int LStartNum = 101;
|
||||
float nipThre = 0;
|
||||
Aurora::Matrix clusImg;
|
||||
Aurora::Matrix nipRe;
|
||||
while(LStartNum > 100)
|
||||
{
|
||||
if(ite < 5)
|
||||
{
|
||||
size_t size = imgSamp.getDataSize();
|
||||
imgSamp.forceReshape(size,1,1);
|
||||
KmeansResult kResult = kmeans(imgSamp, 2);
|
||||
clusImg = kResult.Label;
|
||||
nipThre = max(kResult.CenterPoint, FunctionDirection::All)[0];
|
||||
ite++;
|
||||
}
|
||||
else
|
||||
{
|
||||
nipThre+=0.02;
|
||||
Aurora::Matrix oneMatrix = zeros(clusImg.getDataSize(),1) + 1;
|
||||
for(int i=0; i<clusImg.getDataSize(); ++i)
|
||||
{
|
||||
if(imgSamp[i] > nipThre)
|
||||
{
|
||||
clusImg[i] += 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
nipRe = zeros(aMatrix.getDimSize(0), aMatrix.getDimSize(1),aMatrix.getDimSize(2));
|
||||
size_t effectiveDataSize = aMatrix.getDataSize() - ((aTopZ + 1) * aMatrix.getDimSize(0) * aMatrix.getDimSize(1));
|
||||
for(size_t i=0; i<effectiveDataSize; ++i)
|
||||
{
|
||||
if(aMatrix[i] > nipThre)
|
||||
{
|
||||
nipRe[i] = 1;
|
||||
}
|
||||
}
|
||||
size_t nipRe20Size = nipRe.getDimSize(0) * nipRe.getDimSize(1) * 20;
|
||||
float* nipRe20Data = new float[nipRe20Size];
|
||||
std::copy(nipRe.getData(), nipRe.getData() + nipRe20Size, nipRe20Data);
|
||||
Aurora::Matrix nipRe20 = Aurora::Matrix::fromRawData(nipRe20Data, nipRe.getDimSize(0), nipRe.getDimSize(1), 20);
|
||||
LStartNum = bwconncomp3D(nipRe20);
|
||||
nipRe = maxConnectedComponent3D(nipRe);
|
||||
}
|
||||
std::vector<float> C3D(clusNum, FLT_MAX);
|
||||
for(size_t i=0; i<imgSamp.getDataSize(); ++i)
|
||||
{
|
||||
if( imgSamp[i] < C3D[clusImg[i] - 1])
|
||||
{
|
||||
C3D[clusImg[i] - 1] = imgSamp[i];
|
||||
}
|
||||
}
|
||||
auto it = std::max_element(C3D.begin(), C3D.end());
|
||||
|
||||
Data3DClusterResult result;
|
||||
result.nipRe = nipRe;
|
||||
result.thre3D = *it;
|
||||
return result;
|
||||
}
|
||||
|
||||
@@ -1,14 +0,0 @@
|
||||
#ifndef DATA3DCLUSTER_H
|
||||
#define DATA3DCLUSTER_H
|
||||
#include <Matrix.h>
|
||||
namespace Recon
|
||||
{
|
||||
struct Data3DClusterResult
|
||||
{
|
||||
Aurora::Matrix nipRe;
|
||||
float thre3D;
|
||||
};
|
||||
Data3DClusterResult Data3DCluster(const Aurora::Matrix& aMatrix, float aTopZ);
|
||||
|
||||
}
|
||||
#endif //DATA3DCLUSTER_H
|
||||
@@ -1,17 +0,0 @@
|
||||
#include "DealBW.h"
|
||||
|
||||
#include "SegImgFunction.h"
|
||||
|
||||
using namespace Recon;
|
||||
using namespace Aurora;
|
||||
|
||||
Aurora::Matrix Recon::DealBW(const Aurora::Matrix& aMatrix)
|
||||
{
|
||||
Aurora::Matrix tmp3 = maxConnRegoin(aMatrix);
|
||||
tmp3 = imdilate(tmp3);
|
||||
tmp3 = imerode(tmp3);
|
||||
tmp3 = imfill(tmp3);
|
||||
tmp3 = maxConnRegoin(tmp3);
|
||||
tmp3 = bwconvhull(tmp3);
|
||||
return tmp3;
|
||||
}
|
||||
@@ -1,8 +0,0 @@
|
||||
#ifndef DEALBW_H
|
||||
#define DEALBW_H
|
||||
#include <Matrix.h>
|
||||
namespace Recon
|
||||
{
|
||||
Aurora::Matrix DealBW(const Aurora::Matrix& aMatrix);
|
||||
}
|
||||
#endif //DEALBW_H
|
||||
@@ -1,38 +0,0 @@
|
||||
#include "FindStartZSliceBaseCluster.h"
|
||||
|
||||
#include "Function2D.h"
|
||||
#include "SegImgFunction.h"
|
||||
|
||||
#include <math.h>
|
||||
|
||||
using namespace Recon;
|
||||
using namespace Aurora;
|
||||
|
||||
unsigned int Recon::FindStartZSliceBaseCluster(const Aurora::Matrix& aMatrix, float aReso, int aTopZ)
|
||||
{
|
||||
unsigned int startZ;
|
||||
float nipThre = M_PI * pow(roundf32(5e-3 / aReso), 2);
|
||||
int i=1;
|
||||
while(i < aMatrix.getDimSize(2) - aTopZ)
|
||||
{
|
||||
Aurora::Matrix label = aMatrix.block(2, i-1, i-1);
|
||||
if(sum(label, FunctionDirection::All)[0] > 0)
|
||||
{
|
||||
label = maxConnRegoin(label);
|
||||
if(sum(label, FunctionDirection::All)[0] < nipThre)
|
||||
{
|
||||
++i;
|
||||
}
|
||||
else
|
||||
{
|
||||
startZ = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
++i;
|
||||
}
|
||||
}
|
||||
return startZ;
|
||||
}
|
||||
@@ -1,8 +0,0 @@
|
||||
#ifndef FINDSTARTZSLICEBASECLUSTER_H
|
||||
#define FINDSTARTZSLICEBASECLUSTER_H
|
||||
#include <Matrix.h>
|
||||
namespace Recon
|
||||
{
|
||||
unsigned int FindStartZSliceBaseCluster(const Aurora::Matrix& aMatrix, float aReso, int aTopZ);
|
||||
}
|
||||
#endif //FINDSTARTZSLICEBASECLUSTER_H
|
||||
@@ -1,45 +0,0 @@
|
||||
#include "PreSegImg.h"
|
||||
#include "SegImgFunction.h"
|
||||
#include "Function3D.h"
|
||||
#include "Function2D.h"
|
||||
#include "DealBW.h"
|
||||
|
||||
using namespace Recon;
|
||||
using namespace Aurora;
|
||||
|
||||
PreSegImgResult Recon::PreSegImg(const Aurora::Matrix& aNipRe, unsigned int aTopZ, unsigned int aStartZ)
|
||||
{
|
||||
int x = aNipRe.getDimSize(0);
|
||||
int y = aNipRe.getDimSize(1);
|
||||
int z = aNipRe.getDimSize(2);
|
||||
Aurora::Matrix SegRes = zeros(aNipRe.getDimSize(0), aNipRe.getDimSize(1), aNipRe.getDimSize(2));
|
||||
unsigned int endZ = aStartZ;
|
||||
for(int i= z - aTopZ; i>=aStartZ; --i)
|
||||
{
|
||||
Aurora::Matrix label = aNipRe.block(2, i-1, i-1);
|
||||
Aurora::Matrix labelH = imfill(label);
|
||||
Aurora::Matrix hole = labelH - label;
|
||||
if(sum(hole, FunctionDirection::All)[0] > 2*sum(label, FunctionDirection::All)[0])
|
||||
{
|
||||
endZ = i;
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
for(int i=aStartZ; i<=endZ; ++i)
|
||||
{
|
||||
Aurora::Matrix label = aNipRe.block(2, i-1, i-1);
|
||||
Aurora::Matrix maskIn = SegRes.block(2, i-2, i-2);
|
||||
Aurora::Matrix maskOut = imerode(maskIn);
|
||||
label = label * (1-maskOut);
|
||||
Aurora::Matrix SegI = (label + maskIn) > 0;
|
||||
SegI = maxConnRegoin(SegI);
|
||||
SegI = DealBW(SegI);
|
||||
SegRes.setBlock(2, i-1, i-1, SegI);
|
||||
}
|
||||
|
||||
PreSegImgResult result;
|
||||
result.endZ = endZ;
|
||||
result.SegRes = SegRes;
|
||||
return result;
|
||||
}
|
||||
@@ -1,14 +0,0 @@
|
||||
#ifndef PRESEGIMG_H
|
||||
#define PRESEGIMG_H
|
||||
#include <Matrix.h>
|
||||
namespace Recon
|
||||
{
|
||||
struct PreSegImgResult
|
||||
{
|
||||
Aurora::Matrix SegRes;
|
||||
unsigned int endZ;
|
||||
};
|
||||
|
||||
PreSegImgResult PreSegImg(const Aurora::Matrix& aNipRe, unsigned int aTopZ, unsigned int aStartZ);
|
||||
}
|
||||
#endif //PRESEGIMG_H
|
||||
@@ -1,170 +0,0 @@
|
||||
#include "SegImg.h"
|
||||
#include "config/config.h"
|
||||
#include "Function2D.h"
|
||||
#include "Function3D.h"
|
||||
#include "SegImgFunction.h"
|
||||
#include "DealBW.h"
|
||||
#include "Data3DCluster.h"
|
||||
#include "FindStartZSliceBaseCluster.h"
|
||||
#include "PreSegImg.h"
|
||||
|
||||
#include <float.h>
|
||||
|
||||
using namespace Aurora;
|
||||
|
||||
namespace
|
||||
{
|
||||
Aurora::Matrix create1ToXMatrix(int aEndValue)
|
||||
{
|
||||
float* data = new float[aEndValue];
|
||||
for(int i=1; i<=aEndValue; ++i)
|
||||
{
|
||||
data[i-1] = i;
|
||||
}
|
||||
return Aurora::Matrix::fromRawData(data, 1, aEndValue);
|
||||
}
|
||||
}
|
||||
|
||||
Aurora::Matrix Recon::SegImgMain(const Aurora::Matrix& aMatrix, float aImageResolution, const Aurora::Matrix& aImageStartPoint, const Aurora::Matrix& aImageEndPoint)
|
||||
{
|
||||
int topZ = roundf32((aImageEndPoint.getData()[2]-0) / aImageResolution);
|
||||
Aurora::Matrix temp = aMatrix.deepCopy();
|
||||
temp.forceReshape(temp.getDataSize(), 1, 1);
|
||||
Aurora::Matrix img = (aMatrix - min(temp,FunctionDirection::All)) / (max(temp,FunctionDirection::All) - min(temp,FunctionDirection::All));
|
||||
Data3DClusterResult data3DClusterResult = Data3DCluster(img, topZ);
|
||||
unsigned int startZ = FindStartZSliceBaseCluster(data3DClusterResult.nipRe, aImageResolution, topZ);
|
||||
PreSegImgResult preSegImgResult = PreSegImg(data3DClusterResult.nipRe, topZ, startZ);
|
||||
Aurora::Matrix segRes = SegImg(preSegImgResult.SegRes, img, preSegImgResult.endZ, data3DClusterResult.thre3D, aImageResolution, aImageStartPoint, aImageEndPoint);
|
||||
return aMatrix * segRes;
|
||||
}
|
||||
|
||||
Aurora::Matrix Recon::SegImg(const Aurora::Matrix& aSegRes, const Aurora::Matrix& aImg, int aEndZ, float aThre3D,
|
||||
float aImageResolution, const Aurora::Matrix& aImageStartPoint, const Aurora::Matrix& aImageEndPoint)
|
||||
{
|
||||
Aurora::Matrix result = aSegRes.deepCopy();
|
||||
int clustNumZ = 2;
|
||||
int x = result.getDimSize(0);
|
||||
int y = result.getDimSize(1);
|
||||
int z = result.getDimSize(2);
|
||||
Aurora::Matrix xIdx = transpose(repmat(create1ToXMatrix(x),y,1));
|
||||
Aurora::Matrix yIdx = repmat(create1ToXMatrix(y),x,1);
|
||||
float CenterP1 = roundf32(-aImageStartPoint[0] / aImageResolution);
|
||||
float CenterP2 = roundf32(-aImageStartPoint[1] / aImageResolution);
|
||||
Aurora::Matrix xDistance = (xIdx-CenterP1)^2;
|
||||
Aurora::Matrix yDistance = (yIdx-CenterP2)^2;
|
||||
Aurora::Matrix centerDis = Aurora::sqrt(xDistance + yDistance);
|
||||
for(int i=aEndZ+1; i<=z; ++i)
|
||||
{
|
||||
int updateFlag = 1;
|
||||
Aurora::Matrix tmp = result.block(2, i-2, i-2);
|
||||
float zDis = aImageEndPoint[2] - (z-i)*aImageResolution;
|
||||
Aurora::Matrix centerDisZ = centerDis*tmp;
|
||||
float centerDis_Z = Aurora::max(centerDisZ)[0]*aImageResolution;
|
||||
float ThSlice;
|
||||
if(centerDis_Z > 0.08 || zDis > -0.03)
|
||||
{
|
||||
tmp = aImg.block(2, i-1, i-1);
|
||||
Aurora::Matrix imgData = downSample(tmp);
|
||||
Aurora::Matrix imgDataCopy = imgData.deepCopy();
|
||||
imgDataCopy.forceReshape(imgData.getDataSize(),1,1);
|
||||
KmeansResult kResult = kmeans(imgDataCopy, clustNumZ);
|
||||
Aurora::Matrix clusImg = kResult.Label;
|
||||
Aurora::Matrix clus = zeros(clustNumZ,2);
|
||||
for(int k=1; k<=clustNumZ;++k)
|
||||
{
|
||||
float cIDx = 0;
|
||||
float minValue = MAXFLOAT;
|
||||
for(unsigned int j=0; j<imgData.getDataSize(); ++j)
|
||||
{
|
||||
if(clusImg[j] == k)
|
||||
{
|
||||
cIDx++;
|
||||
if(minValue > imgData[j])
|
||||
{
|
||||
minValue = imgData[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
clus[k-1] = minValue;
|
||||
clus[k+1] = cIDx;
|
||||
}
|
||||
Aurora::Matrix idx;
|
||||
clus = sortrows(clus, &idx);
|
||||
if(clus[0] < clus[1] && clus[2] > clus[3])
|
||||
{
|
||||
ThSlice = clus[1];
|
||||
}
|
||||
else
|
||||
{
|
||||
tmp = aImg.block(2, i-1, i-1);
|
||||
Aurora::Matrix mask = result.block(2, i-2, i-2);
|
||||
mask = imerode(mask);
|
||||
tmp = tmp*(1-mask) + mask * mean(tmp,FunctionDirection::All)[0];
|
||||
imgData = downSample(tmp);
|
||||
Aurora::Matrix imgDataCopy = imgData.deepCopy();
|
||||
imgDataCopy.forceReshape(imgData.getDataSize(),1,1);
|
||||
kResult = kmeans(imgDataCopy, 2);
|
||||
clusImg = kResult.Label;
|
||||
clus = zeros(clustNumZ,2);
|
||||
for(int k=1; k<=clustNumZ;++k)
|
||||
{
|
||||
float cIDx = 0;
|
||||
float minValue = MAXFLOAT;
|
||||
for(unsigned int j=0; j<imgData.getDataSize(); ++j)
|
||||
{
|
||||
if(clusImg[j] == k)
|
||||
{
|
||||
cIDx++;
|
||||
if(minValue > imgData[j])
|
||||
{
|
||||
minValue = imgData[j];
|
||||
}
|
||||
}
|
||||
}
|
||||
clus[k-1] = minValue;
|
||||
clus[k+1] = cIDx;
|
||||
}
|
||||
Aurora::Matrix idx;
|
||||
clus = sortrows(clus, &idx);
|
||||
if(clus[0] < clus[1] && clus[2] > clus[3])
|
||||
{
|
||||
ThSlice = clus[1];
|
||||
}
|
||||
else
|
||||
{
|
||||
updateFlag = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
else
|
||||
{
|
||||
ThSlice = aThre3D;
|
||||
}
|
||||
|
||||
if(zDis > -0.03)
|
||||
{
|
||||
Aurora::Matrix tmp1 = result.block(2, i-2, i-2);
|
||||
Aurora::Matrix tmp2 = result.block(2, i-1, i-1);
|
||||
float idx = ((sum(tmp2,FunctionDirection::All) - sum(tmp1,FunctionDirection::All)) / sum(tmp1,FunctionDirection::All)).getData()[0];
|
||||
if(idx>0.15)
|
||||
{
|
||||
updateFlag = 0;
|
||||
}
|
||||
}
|
||||
|
||||
Aurora::Matrix segI;
|
||||
if(updateFlag == 1)
|
||||
{
|
||||
Aurora::Matrix ref = aImg > ThSlice;
|
||||
Aurora::Matrix refI = ref.block(2, i-1, i-1);
|
||||
segI = (refI + result.block(2, i-2, i-2)) > 0;
|
||||
segI = DealBW(segI);
|
||||
}
|
||||
else
|
||||
{
|
||||
segI = result.block(2, i-2, i-2);
|
||||
}
|
||||
result.setBlock(2, i-1, i-1, segI);
|
||||
}
|
||||
return result;
|
||||
}
|
||||
@@ -1,10 +0,0 @@
|
||||
#ifndef SEGIMG_H
|
||||
#define SEGIMG_H
|
||||
#include <Matrix.h>
|
||||
namespace Recon
|
||||
{
|
||||
Aurora::Matrix SegImgMain(const Aurora::Matrix& aMatrix, float aImageResolution, const Aurora::Matrix& aImageStartPoint, const Aurora::Matrix& aImageEndPoint);
|
||||
Aurora::Matrix SegImg(const Aurora::Matrix& aSegRes, const Aurora::Matrix& aImg, int aEndZ, float aThre3D,
|
||||
float aImageResolution, const Aurora::Matrix& aImageStartPoint, const Aurora::Matrix& aImageEndPoint);
|
||||
}
|
||||
#endif //SEGIMG_H
|
||||
@@ -1,218 +0,0 @@
|
||||
#include "SegImgFunction.h"
|
||||
|
||||
#include <opencv2/opencv.hpp>
|
||||
|
||||
#include "Function3D.h"
|
||||
|
||||
|
||||
using namespace Recon;
|
||||
|
||||
Aurora::Matrix Recon::maxConnRegoin(const Aurora::Matrix& aMatrix)
|
||||
{
|
||||
int column = aMatrix.getDimSize(0);
|
||||
int row = aMatrix.getDimSize(1);
|
||||
size_t size = column * row;
|
||||
uint8_t* data = new uint8_t[size];
|
||||
|
||||
std::copy(aMatrix.getData(), aMatrix.getData() + size, data);
|
||||
cv::Mat src = cv::Mat(row, column, CV_8U, data);
|
||||
cv::Mat labels, stats, centroids;;
|
||||
int nLabels = connectedComponentsWithStats(src, labels,stats,centroids, 4, CV_32S);
|
||||
if(nLabels<=1)
|
||||
{
|
||||
return Aurora::Matrix();
|
||||
}
|
||||
int maxArea = 0;
|
||||
int maxAreaIndex = 1;
|
||||
for (int i = 1; i < nLabels; ++i)
|
||||
{
|
||||
int area = stats.at<int>(i, cv::CC_STAT_AREA);
|
||||
if (area > maxArea) {
|
||||
maxArea = area;
|
||||
maxAreaIndex = i;
|
||||
}
|
||||
}
|
||||
|
||||
Aurora::Matrix result = Aurora::zeros(column, row);
|
||||
for(int i=0; i<row; ++i)
|
||||
{
|
||||
for(int j=0; j<column; ++j)
|
||||
{
|
||||
if (labels.at<int>(i, j) == maxAreaIndex)
|
||||
{
|
||||
result[i*column + j] = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
return result;
|
||||
}
|
||||
|
||||
Aurora::Matrix Recon::imfill(const Aurora::Matrix& aMatrix)
|
||||
{
|
||||
|
||||
int width = aMatrix.getDimSize(0);
|
||||
int height = aMatrix.getDimSize(1);
|
||||
size_t dataSize = aMatrix.getDataSize();
|
||||
uint8_t* data = new uint8_t[aMatrix.getDataSize()];
|
||||
std::copy(aMatrix.getData(), aMatrix.getData() + dataSize, data);
|
||||
cv::Mat src = cv::Mat(height, width, CV_8U, data);
|
||||
cv::Mat inv = 1 - src;
|
||||
cv::Mat labels;
|
||||
int nLabels = connectedComponents(inv, labels, 8, CV_32S);
|
||||
cv::Mat filled = src.clone();
|
||||
std::vector<bool> isHole(nLabels, true);
|
||||
|
||||
int rows = src.rows;
|
||||
int cols = src.cols;
|
||||
|
||||
for (int x = 0; x < cols; x++) {
|
||||
isHole[ labels.at<int>(0, x) ] = false;
|
||||
isHole[ labels.at<int>(rows-1, x) ] = false;
|
||||
}
|
||||
for (int y = 0; y < rows; y++) {
|
||||
isHole[ labels.at<int>(y, 0) ] = false;
|
||||
isHole[ labels.at<int>(y, cols-1) ] = false;
|
||||
}
|
||||
|
||||
for (int y = 0; y < rows; y++) {
|
||||
for (int x = 0; x < cols; x++) {
|
||||
int label = labels.at<int>(y, x);
|
||||
if (isHole[label]) {
|
||||
filled.at<uchar>(y, x) = 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
float* resultData = new float[aMatrix.getDataSize()];
|
||||
|
||||
std::copy(filled.data, filled.data + dataSize, resultData);
|
||||
Aurora::Matrix result = Aurora::Matrix::fromRawData(resultData, width, height);
|
||||
delete[] data;
|
||||
return result;
|
||||
}
|
||||
|
||||
Aurora::Matrix Recon::bwconvhull(const Aurora::Matrix& aMatrix)
|
||||
{
|
||||
int width = aMatrix.getDimSize(0);
|
||||
int height = aMatrix.getDimSize(1);
|
||||
size_t dataSize = aMatrix.getDataSize();
|
||||
uint8_t* data = new uint8_t[aMatrix.getDataSize()];
|
||||
std::copy(aMatrix.getData(), aMatrix.getData() + dataSize, data);
|
||||
cv::Mat input = cv::Mat(height, width, CV_8U, data);
|
||||
std::vector<cv::Point> points;
|
||||
cv::findNonZero(input, points);
|
||||
|
||||
if (points.empty())
|
||||
{
|
||||
delete[] data;
|
||||
return Aurora::zeros(width, height);;
|
||||
}
|
||||
|
||||
std::vector<cv::Point> hull;
|
||||
cv::convexHull(points, hull);
|
||||
|
||||
|
||||
cv::Mat hullMask = cv::Mat::zeros(input.size(), CV_8UC1);
|
||||
cv::fillConvexPoly(hullMask, hull, cv::Scalar(1));
|
||||
float* resultData = new float[aMatrix.getDataSize()];
|
||||
|
||||
std::copy(hullMask.data, hullMask.data + dataSize, resultData);
|
||||
Aurora::Matrix result = Aurora::Matrix::fromRawData(resultData, width, height);
|
||||
|
||||
delete[] data;
|
||||
return result;
|
||||
}
|
||||
|
||||
KmeansResult Recon::kmeans(const Aurora::Matrix& aMatrix, int aClusNum)
|
||||
{
|
||||
KmeansResult result;
|
||||
auto input = aMatrix.deepCopy();
|
||||
size_t dataSize = input.getDataSize();
|
||||
cv::Mat data(dataSize, 1, CV_32F, input.getData());
|
||||
|
||||
int k = aClusNum;
|
||||
cv::Mat labels;
|
||||
cv::Mat centers;
|
||||
cv::TermCriteria criteria(cv::TermCriteria::EPS + cv::TermCriteria::MAX_ITER, 100, 1e-6);
|
||||
int attempts = 5;
|
||||
int flags = cv::KMEANS_PP_CENTERS;
|
||||
double compactness = cv::kmeans(data, k, labels, criteria, attempts, flags, centers);
|
||||
|
||||
// std::cout << "聚类中心: " << std::endl << centers << std::endl;
|
||||
// std::cout << "紧密度: " << compactness << std::endl;
|
||||
|
||||
float* labelDataResult = new float[dataSize];
|
||||
float* centerPointDataResult = new float[k];
|
||||
int* labelData = reinterpret_cast<int*>(labels.data);
|
||||
float* centerPointData = reinterpret_cast<float*>(centers.data);
|
||||
|
||||
std::copy(labelData, labelData + dataSize, labelDataResult);
|
||||
std::copy(centerPointData, centerPointData + k, centerPointDataResult);
|
||||
Aurora::Matrix label = Aurora::Matrix::fromRawData(labelDataResult, dataSize, 1) + 1;
|
||||
Aurora::Matrix center = Aurora::Matrix::fromRawData(centerPointDataResult, k, 1);
|
||||
result.Label = label;
|
||||
result.CenterPoint = center;
|
||||
return result;
|
||||
}
|
||||
|
||||
Aurora::Matrix Recon::imerode(const Aurora::Matrix& aMatrix)
|
||||
{
|
||||
float* data = aMatrix.getData();
|
||||
unsigned int width = aMatrix.getDimSize(0);
|
||||
unsigned int heigh = aMatrix.getDimSize(1);
|
||||
unsigned int dataSize = width * heigh;
|
||||
cv::Mat kernel = cv::getStructuringElement(cv::MORPH_RECT, cv::Size(3,3));
|
||||
cv::Mat mat(heigh, width, CV_32F, data );
|
||||
cv::Mat out;
|
||||
cv::erode(mat, out, kernel);
|
||||
float* outData = reinterpret_cast<float*>(out.data);
|
||||
float* resultData = new float[dataSize];
|
||||
std::copy(outData, outData + dataSize, resultData);
|
||||
return Aurora::Matrix::fromRawData(resultData, width, heigh);
|
||||
}
|
||||
|
||||
Aurora::Matrix Recon::imdilate(const Aurora::Matrix& aMatrix)
|
||||
{
|
||||
float* data = aMatrix.getData();
|
||||
unsigned int width = aMatrix.getDimSize(0);
|
||||
unsigned int heigh = aMatrix.getDimSize(1);
|
||||
unsigned int dataSize = width * heigh;
|
||||
cv::Mat kernel = cv::getStructuringElement(cv::MORPH_RECT, cv::Size(3,3));
|
||||
cv::Mat mat(heigh, width, CV_32F, data );
|
||||
cv::Mat out;
|
||||
cv::dilate(mat, out, kernel);
|
||||
float* outData = reinterpret_cast<float*>(out.data);
|
||||
float* resultData = new float[dataSize];
|
||||
std::copy(outData, outData + dataSize, resultData);
|
||||
return Aurora::Matrix::fromRawData(resultData, width, heigh);
|
||||
}
|
||||
|
||||
|
||||
Aurora::Matrix Recon::downSample(const Aurora::Matrix& aMatrix)
|
||||
{
|
||||
int row = aMatrix.getDimSize(0);
|
||||
int column = aMatrix.getDimSize(1);
|
||||
int slice = aMatrix.getDimSize(2);
|
||||
int newRow = (row + 3) / 4;
|
||||
int newColumn = (column +3) / 4;
|
||||
int newSlice = (slice + 3) / 4;
|
||||
float* resultData = new float[newColumn * newRow * newSlice];
|
||||
for (int s = 0; s < newSlice; ++s)
|
||||
{
|
||||
int oldSlice = s * 4;
|
||||
if (oldSlice >= slice) oldSlice = slice - 1;
|
||||
for (int c = 0; c < newColumn; ++c)
|
||||
{
|
||||
int oldColumn = c * 4;
|
||||
if (oldColumn >= column) oldColumn = column - 1;
|
||||
for (int r = 0; r < newRow; ++r)
|
||||
{
|
||||
int oldRow = r * 4;
|
||||
if (oldRow >= row) oldRow = row - 1;
|
||||
resultData[s * newRow * newColumn + c * newRow + r] = aMatrix.getData()[oldSlice * row * column + oldColumn * row + oldRow];
|
||||
}
|
||||
}
|
||||
}
|
||||
return Aurora::Matrix::fromRawData(resultData, newRow, newColumn, newSlice);
|
||||
}
|
||||
@@ -1,20 +0,0 @@
|
||||
#ifndef SEGIMGFUNCTION_H
|
||||
#define SEGIMGFUNCTION_H
|
||||
#include <Matrix.h>
|
||||
namespace Recon
|
||||
{
|
||||
struct KmeansResult
|
||||
{
|
||||
Aurora::Matrix Label;
|
||||
Aurora::Matrix CenterPoint;
|
||||
};
|
||||
|
||||
Aurora::Matrix maxConnRegoin(const Aurora::Matrix& aMatrix);
|
||||
Aurora::Matrix imfill(const Aurora::Matrix& aMatrix);
|
||||
Aurora::Matrix bwconvhull(const Aurora::Matrix& aMatrix);
|
||||
Aurora::Matrix imdilate(const Aurora::Matrix& aMatrix);
|
||||
KmeansResult kmeans(const Aurora::Matrix& aMatrix, int aClusNum);
|
||||
Aurora::Matrix imerode(const Aurora::Matrix& aMatrix);
|
||||
Aurora::Matrix downSample(const Aurora::Matrix& aMatrix);
|
||||
}
|
||||
#endif //SEGIMGFUNCTION_H
|
||||
@@ -1,263 +0,0 @@
|
||||
#include "maxConnectedComponent3D.cuh"
|
||||
#include <cuda_runtime.h>
|
||||
#include <cstdio>
|
||||
#include <cstdint>
|
||||
#include <unordered_set>
|
||||
#include <unordered_map>
|
||||
#include <vector>
|
||||
#include <algorithm>
|
||||
|
||||
#define INDEX(x,y,z,W,C,S) ((x)*(C) + (y) + (z)*(W)*(C))
|
||||
|
||||
__constant__ int NEIGH_OFFS[26*3];
|
||||
|
||||
static int HOST_OFFSET[26*3];
|
||||
|
||||
__global__ void floatToMask(const float* __restrict__ in,
|
||||
uint8_t* __restrict__ mask,
|
||||
int W, int C, int S)
|
||||
{
|
||||
int x = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
int y = blockIdx.y * blockDim.y + threadIdx.y;
|
||||
int z = blockIdx.z * blockDim.z + threadIdx.z;
|
||||
if (x>=W || y>=C || z>=S) return;
|
||||
|
||||
size_t idx = INDEX(x,y,z,W,C,S);
|
||||
float v = in[idx];
|
||||
mask[idx] = (v > 0.5f) ? 1 : 0;
|
||||
}
|
||||
|
||||
__global__ void initLabels(const uint8_t* __restrict__ mask,
|
||||
uint32_t* __restrict__ labels,
|
||||
int W, int C, int S)
|
||||
{
|
||||
int x = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
int y = blockIdx.y * blockDim.y + threadIdx.y;
|
||||
int z = blockIdx.z * blockDim.z + threadIdx.z;
|
||||
if (x>=W || y>=C || z>=S) return;
|
||||
|
||||
size_t idx = INDEX(x,y,z,W,C,S);
|
||||
labels[idx] = mask[idx] ? (uint32_t)(idx + 1) : 0u;
|
||||
}
|
||||
|
||||
__global__ void propagateLabels(const uint8_t* __restrict__ mask,
|
||||
uint32_t* __restrict__ labels,
|
||||
int W, int C, int S,
|
||||
int* __restrict__ d_changed)
|
||||
{
|
||||
int x = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
int y = blockIdx.y * blockDim.y + threadIdx.y;
|
||||
int z = blockIdx.z * blockDim.z + threadIdx.z;
|
||||
if (x>=W || y>=C || z>=S) return;
|
||||
|
||||
size_t idx = INDEX(x,y,z,W,C,S);
|
||||
if (!mask[idx]) return;
|
||||
|
||||
uint32_t cur = labels[idx];
|
||||
uint32_t minLabel = cur;
|
||||
|
||||
// 遍历 26 邻域
|
||||
#pragma unroll
|
||||
for (int k = 0; k < 26; ++k)
|
||||
{
|
||||
int nx = x + NEIGH_OFFS[3*k + 0];
|
||||
int ny = y + NEIGH_OFFS[3*k + 1];
|
||||
int nz = z + NEIGH_OFFS[3*k + 2];
|
||||
if ((unsigned)nx < (unsigned)W &&
|
||||
(unsigned)ny < (unsigned)C &&
|
||||
(unsigned)nz < (unsigned)S)
|
||||
{
|
||||
size_t nidx = INDEX(nx,ny,nz,W,C,S);
|
||||
if (mask[nidx]) {
|
||||
uint32_t nl = labels[nidx];
|
||||
if (nl != 0u && nl < minLabel) minLabel = nl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
if (minLabel < cur)
|
||||
{
|
||||
labels[idx] = minLabel;
|
||||
*d_changed = 1;
|
||||
}
|
||||
}
|
||||
|
||||
__global__ void writeLargestComponent(const uint32_t* __restrict__ labels,
|
||||
float* __restrict__ out,
|
||||
uint32_t keep_label,
|
||||
int W, int C, int S)
|
||||
{
|
||||
int x = blockIdx.x * blockDim.x + threadIdx.x;
|
||||
int y = blockIdx.y * blockDim.y + threadIdx.y;
|
||||
int z = blockIdx.z * blockDim.z + threadIdx.z;
|
||||
if (x>=W || y>=C || z>=S) return;
|
||||
|
||||
size_t idx = INDEX(x,y,z,W,C,S);
|
||||
out[idx] = (labels[idx] == keep_label) ? 1.0f : 0.0f;
|
||||
}
|
||||
|
||||
static uint32_t findLargestLabel(const std::vector<uint32_t>& h_labels)
|
||||
{
|
||||
std::unordered_map<uint32_t, uint32_t> cnt;
|
||||
cnt.reserve(h_labels.size()/16 + 1);
|
||||
for (auto L : h_labels) {
|
||||
if (L) ++cnt[L];
|
||||
}
|
||||
uint32_t best_label = 0;
|
||||
uint32_t best_size = 0;
|
||||
for (auto& kv : cnt) {
|
||||
if (kv.second > best_size) {
|
||||
best_size = kv.second;
|
||||
best_label = kv.first;
|
||||
}
|
||||
}
|
||||
return best_label;
|
||||
}
|
||||
|
||||
unsigned int bwconncomp3D(const Aurora::Matrix aInput, int max_iters)
|
||||
{
|
||||
int W = aInput.getDimSize(0);
|
||||
int C = aInput.getDimSize(1);
|
||||
int S = aInput.getDimSize(2);
|
||||
float* inVol = aInput.getData();
|
||||
size_t N = W * C * S;
|
||||
if (S <= 1) return 0;
|
||||
|
||||
|
||||
{
|
||||
int t = 0;
|
||||
for (int dz = -1; dz <= 1; ++dz)
|
||||
for (int dx = -1; dx <= 1; ++dx)
|
||||
for (int dy = -1; dy <= 1; ++dy) {
|
||||
if (dx==0 && dy==0 && dz==0) continue;
|
||||
HOST_OFFSET[3*t+0] = dx;
|
||||
HOST_OFFSET[3*t+1] = dy;
|
||||
HOST_OFFSET[3*t+2] = dz;
|
||||
++t;
|
||||
}
|
||||
cudaMemcpyToSymbol(NEIGH_OFFS, HOST_OFFSET, sizeof(HOST_OFFSET));
|
||||
}
|
||||
|
||||
float *d_in = nullptr, *d_out = nullptr;
|
||||
uint8_t *d_mask = nullptr;
|
||||
uint32_t *d_labels = nullptr;
|
||||
int *d_changed = nullptr;
|
||||
|
||||
cudaMalloc(&d_in, N * sizeof(float));
|
||||
cudaMalloc(&d_out, N * sizeof(float));
|
||||
cudaMalloc(&d_mask, N * sizeof(uint8_t));
|
||||
cudaMalloc(&d_labels, N * sizeof(uint32_t));
|
||||
cudaMalloc(&d_changed,sizeof(int));
|
||||
|
||||
cudaMemcpy(d_in, inVol, N*sizeof(float), cudaMemcpyHostToDevice);
|
||||
|
||||
dim3 block(8,8,8);
|
||||
dim3 grid( (W+block.x-1)/block.x,
|
||||
(C+block.y-1)/block.y,
|
||||
(S+block.z-1)/block.z );
|
||||
|
||||
floatToMask<<<grid, block>>>(d_in, d_mask, W, C, S);
|
||||
|
||||
initLabels<<<grid, block>>>(d_mask, d_labels, W, C, S);
|
||||
|
||||
int h_changed = 0;
|
||||
int iter = 0;
|
||||
do
|
||||
{
|
||||
h_changed = 0;
|
||||
cudaMemcpy(d_changed, &h_changed, sizeof(int), cudaMemcpyHostToDevice);
|
||||
propagateLabels<<<grid, block>>>(d_mask, d_labels, W, C, S, d_changed);
|
||||
cudaMemcpy(&h_changed, d_changed, sizeof(int), cudaMemcpyDeviceToHost);
|
||||
++iter;
|
||||
if (iter >= max_iters)
|
||||
{
|
||||
fprintf(stderr, "[Warn] Reached max_iters=%d before convergence.\n", max_iters);
|
||||
break;
|
||||
}
|
||||
} while (h_changed != 0);
|
||||
std::vector<uint32_t> h_labels(N);
|
||||
cudaMemcpy(h_labels.data(), d_labels, N*sizeof(uint32_t), cudaMemcpyDeviceToHost);
|
||||
std::unordered_set<uint32_t> uniqueSet(h_labels.begin(), h_labels.end());
|
||||
return uniqueSet.size() - 1;
|
||||
}
|
||||
|
||||
Aurora::Matrix maxConnectedComponent3D(const Aurora::Matrix aInput, int max_iters)
|
||||
{
|
||||
int W = aInput.getDimSize(0);
|
||||
int C = aInput.getDimSize(1);
|
||||
int S = aInput.getDimSize(2);
|
||||
float* inVol = aInput.getData();
|
||||
size_t N = W * C * S;
|
||||
if (S <= 1) return Aurora::Matrix();
|
||||
|
||||
|
||||
{
|
||||
int t = 0;
|
||||
for (int dz = -1; dz <= 1; ++dz)
|
||||
for (int dx = -1; dx <= 1; ++dx)
|
||||
for (int dy = -1; dy <= 1; ++dy) {
|
||||
if (dx==0 && dy==0 && dz==0) continue;
|
||||
HOST_OFFSET[3*t+0] = dx;
|
||||
HOST_OFFSET[3*t+1] = dy;
|
||||
HOST_OFFSET[3*t+2] = dz;
|
||||
++t;
|
||||
}
|
||||
cudaMemcpyToSymbol(NEIGH_OFFS, HOST_OFFSET, sizeof(HOST_OFFSET));
|
||||
}
|
||||
|
||||
float *d_in = nullptr, *d_out = nullptr;
|
||||
uint8_t *d_mask = nullptr;
|
||||
uint32_t *d_labels = nullptr;
|
||||
int *d_changed = nullptr;
|
||||
|
||||
cudaMalloc(&d_in, N * sizeof(float));
|
||||
cudaMalloc(&d_out, N * sizeof(float));
|
||||
cudaMalloc(&d_mask, N * sizeof(uint8_t));
|
||||
cudaMalloc(&d_labels, N * sizeof(uint32_t));
|
||||
cudaMalloc(&d_changed,sizeof(int));
|
||||
|
||||
cudaMemcpy(d_in, inVol, N*sizeof(float), cudaMemcpyHostToDevice);
|
||||
|
||||
dim3 block(8,8,8);
|
||||
dim3 grid( (W+block.x-1)/block.x,
|
||||
(C+block.y-1)/block.y,
|
||||
(S+block.z-1)/block.z );
|
||||
|
||||
floatToMask<<<grid, block>>>(d_in, d_mask, W, C, S);
|
||||
|
||||
initLabels<<<grid, block>>>(d_mask, d_labels, W, C, S);
|
||||
|
||||
int h_changed = 0;
|
||||
int iter = 0;
|
||||
do
|
||||
{
|
||||
h_changed = 0;
|
||||
cudaMemcpy(d_changed, &h_changed, sizeof(int), cudaMemcpyHostToDevice);
|
||||
propagateLabels<<<grid, block>>>(d_mask, d_labels, W, C, S, d_changed);
|
||||
cudaMemcpy(&h_changed, d_changed, sizeof(int), cudaMemcpyDeviceToHost);
|
||||
++iter;
|
||||
if (iter >= max_iters)
|
||||
{
|
||||
fprintf(stderr, "[Warn] Reached max_iters=%d before convergence.\n", max_iters);
|
||||
break;
|
||||
}
|
||||
} while (h_changed != 0);
|
||||
std::vector<uint32_t> h_labels(N);
|
||||
cudaMemcpy(h_labels.data(), d_labels, N*sizeof(uint32_t), cudaMemcpyDeviceToHost);
|
||||
|
||||
uint32_t keep_label = findLargestLabel(h_labels);
|
||||
|
||||
writeLargestComponent<<<grid, block>>>(d_labels, d_out, keep_label, W, C, S);
|
||||
|
||||
float* outputData = new float[N];
|
||||
cudaMemcpy(outputData, d_out, N*sizeof(float), cudaMemcpyDeviceToHost);
|
||||
|
||||
|
||||
cudaFree(d_in);
|
||||
cudaFree(d_out);
|
||||
cudaFree(d_mask);
|
||||
cudaFree(d_labels);
|
||||
cudaFree(d_changed);
|
||||
|
||||
return Aurora::Matrix::fromRawData(outputData, W, C, S);
|
||||
}
|
||||
@@ -1,5 +0,0 @@
|
||||
#include <Matrix.h>
|
||||
|
||||
Aurora::Matrix maxConnectedComponent3D(const Aurora::Matrix aInput, int max_iters = 5000);
|
||||
|
||||
unsigned int bwconncomp3D(const Aurora::Matrix aInput, int max_iters = 5000);
|
||||
143
src/common/SoundSpeedWaterCalculation/estimateSOSWater.cpp
Normal file
143
src/common/SoundSpeedWaterCalculation/estimateSOSWater.cpp
Normal file
@@ -0,0 +1,143 @@
|
||||
#include "estimateSOSWater.h"
|
||||
#include "Function1D.cuh"
|
||||
#include "Function1D.h"
|
||||
#include "Function3D.h"
|
||||
#include "Matrix.h"
|
||||
#include "CudaMatrix.h"
|
||||
#include "common/dataBlockCreation/getAscanBlock.h"
|
||||
#include "common/getMeasurementMetaData.h"
|
||||
#include "transmissionReconstruction/detection/getTransmissionData.cuh"
|
||||
#include "transmissionReconstruction/detection/detection.h"
|
||||
#include "config/config.h"
|
||||
#include "config/geometryConfig.h"
|
||||
|
||||
#include "Function2D.cuh"
|
||||
#include "Function2D.h"
|
||||
#include "common/dataBlockCreation/getAscanBlock.h"
|
||||
#include "common/dataBlockCreation/removeDataFromArrays.h"
|
||||
#include "log/log.h"
|
||||
#include <cstring>
|
||||
|
||||
using namespace Recon;
|
||||
|
||||
namespace
|
||||
{
|
||||
const Aurora::Matrix EMIT_TAS = Aurora::Matrix::fromRawData(new float[17]{94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110}, 17);
|
||||
const Aurora::Matrix RECEIVE_TAS = Aurora::Matrix::fromRawData(new float[17]{103, 104, 105, 106, 107, 108, 109, 110, 94, 95, 96, 97, 98, 99, 100, 101, 102}, 17);
|
||||
const Aurora::Matrix ALL_RECEIVER_TAS_LIST = Aurora::Matrix::fromRawData(new float[128] {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128},128); //1~128
|
||||
const Aurora::Matrix ALL_RECEIVER_ELEMENT_LIST = Aurora::Matrix::fromRawData(new float[18] {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18},18);
|
||||
const float THRE = 0.5;
|
||||
inline int findIndexFromEmitterAndReceiver(const Aurora::Matrix& aEmitter, float aEmitterValue, const Aurora::Matrix& aReceiver, float aReceiverValue)
|
||||
{
|
||||
for (int i = 0; i < aEmitter.getDataSize(); ++i)
|
||||
{
|
||||
if (aEmitter[i] == aEmitterValue && aReceiver[i] == aReceiverValue)
|
||||
{
|
||||
return i;
|
||||
}
|
||||
}
|
||||
return -1;
|
||||
}
|
||||
|
||||
Aurora::Matrix matchFilt(const Aurora::CudaMatrix& aSig, const Aurora::CudaMatrix& aMod)
|
||||
{
|
||||
Aurora::CudaMatrix fftMod = fft(aMod);
|
||||
Aurora::CudaMatrix fftSig = fft(aSig);
|
||||
Aurora::CudaMatrix mfFftSig = getTransmissionDataSubFunction(fftSig, fftMod);
|
||||
Aurora::Matrix mfSig = Aurora::real(ifft(mfFftSig)).toHostMatrix();
|
||||
std::memmove(mfSig.getData(), mfSig.getData() + 1, (mfSig.getDataSize() - 1) * sizeof(float));
|
||||
return mfSig;
|
||||
}
|
||||
|
||||
Aurora::Matrix addWin(const Aurora::Matrix& aSig, float aOffset, float aMinSos, float aMaxSos, float aDist, float aFq)
|
||||
{
|
||||
float offset = aOffset * aFq;
|
||||
int start_ = round(aDist / aMaxSos * aFq);
|
||||
int end_ = round(aDist / aMinSos * aFq);
|
||||
float start = start_ + offset;
|
||||
float end = end_ + offset;
|
||||
Aurora::Matrix result = Aurora::zeros(aSig.getDimSize(0), aSig.getDimSize(1));
|
||||
result.setBlock(0, start-1, end-1, aSig.block(0, start-1, end-1));
|
||||
return result;
|
||||
}
|
||||
|
||||
inline int findFirstvalueGreaterThanGivenValue(const Aurora::Matrix& aMatrix, float aValue)
|
||||
{
|
||||
for(int i=0; i<aMatrix.getDataSize(); ++i)
|
||||
{
|
||||
if(aMatrix[i] > aValue)
|
||||
{
|
||||
return i;
|
||||
}
|
||||
}
|
||||
return -1;
|
||||
}
|
||||
}
|
||||
|
||||
Aurora::Matrix Recon::estimateSOSWater(Parser *aParser, const CEInfo &aCE, const PreComputes& aPreComputes)
|
||||
{
|
||||
float offset = -aPreComputes.offset;
|
||||
float thre = 0.99;
|
||||
Aurora::Matrix tof = Aurora::zeros(1, EMIT_TAS.getDataSize() * 18 * 18);
|
||||
Aurora::Matrix distance = tof.deepCopy();
|
||||
Aurora::Matrix mp = Aurora::Matrix::fromRawData(new float[1]{1}, 1);
|
||||
Aurora::Matrix snMatrix = Aurora::Matrix::fromRawData(new float[18]{1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18}, 18);
|
||||
|
||||
for(int idx=0; idx<EMIT_TAS.getDataSize(); ++idx)
|
||||
{
|
||||
Aurora::Matrix slMatrix = Aurora::Matrix::fromRawData(new float[1]{EMIT_TAS[idx]}, 1);
|
||||
AscanBlock ascanTotal = getAscanBlock(aParser, mp, slMatrix, snMatrix, ALL_RECEIVER_TAS_LIST, ALL_RECEIVER_ELEMENT_LIST);
|
||||
#pragma omp parallel for
|
||||
for(int sn=1; sn<=18; ++sn)
|
||||
{
|
||||
// Aurora::Matrix snMatrix = Aurora::Matrix::fromRawData(new float[1]{sn}, 1);
|
||||
// Aurora::Matrix slMatrix = Aurora::Matrix::fromRawData(new float[1]{EMIT_TAS[idx]}, 1);
|
||||
// AscanBlock ascan = getAscanBlock(aParser, mp, slMatrix, snMatrix, ALL_RECEIVER_TAS_LIST, ALL_RECEIVER_ELEMENT_LIST);
|
||||
// size_t dataSize = ascan.ascanBlock.getDataSize();
|
||||
// short* data = new short[dataSize];
|
||||
// std::copy(ascan.ascanBlock.getData(), ascan.ascanBlock.getData() + dataSize, data);
|
||||
// Aurora::Matrix rfSet = Aurora::convertfp16tofloat(data, ascan.ascanBlock.getDimSize(0), ascan.ascanBlock.getDimSize(1));
|
||||
// delete[] data;
|
||||
// Aurora::CudaMatrix rfSetDevice = rfSet.toDeviceMatrix() / ascan.gainBlock.toDeviceMatrix();
|
||||
//AscanBlock ascan = getAscanBlock(aParser, mp, slMatrix, snMatrix, ALL_RECEIVER_TAS_LIST, ALL_RECEIVER_ELEMENT_LIST);
|
||||
size_t dataSize = ascanTotal.ascanBlock.getDataSize() / 18;
|
||||
short* data = new short[dataSize];
|
||||
std::copy(ascanTotal.ascanBlock.getData() + (sn - 1) * dataSize, ascanTotal.ascanBlock.getData() + sn * dataSize, data);
|
||||
Aurora::Matrix rfSet = Aurora::convertfp16tofloat(data, ascanTotal.ascanBlock.getDimSize(0), ascanTotal.ascanBlock.getDimSize(1) / 18);
|
||||
delete[] data;
|
||||
Aurora::CudaMatrix rfSetDevice = rfSet.toDeviceMatrix() / ascanTotal.gainBlock.block(1, sn-1, sn-1).toDeviceMatrix();
|
||||
for(int rn=1; rn<=18; ++rn)
|
||||
{
|
||||
int cnt = idx * 18 * 18 + 18 * (sn - 1) + rn - 1;
|
||||
Aurora::Matrix positionOfEmitter = emitterPositions[EMIT_TAS[idx]-1](sn-1).toMatrix();
|
||||
Aurora::Matrix positionOfReceiver = emitterPositions[RECEIVE_TAS[idx]-1](rn-1).toMatrix();
|
||||
distance[cnt] = Aurora::norm(positionOfEmitter - positionOfReceiver, Aurora::Norm2);
|
||||
int rfIndex = findIndexFromEmitterAndReceiver(aCE.tasIndices, RECEIVE_TAS[idx], aCE.receiverIndices, rn);
|
||||
Aurora::CudaMatrix rf = rfSetDevice.block(1, rfIndex, rfIndex);
|
||||
Aurora::Matrix rfMf;
|
||||
if(aPreComputes.matchedFilter.getDimSize(1) == 1) //Use CE
|
||||
{
|
||||
rfMf = matchFilt(rf, Aurora::real(ifft(aPreComputes.matchedFilter.toDeviceMatrix())));
|
||||
}
|
||||
else // Use CEM
|
||||
{
|
||||
rfMf = matchFilt(rf, Aurora::real(ifft(aPreComputes.matchedFilter.block(1, rfIndex, rfIndex).toDeviceMatrix())));
|
||||
}
|
||||
Aurora::Matrix rfWin = addWin(rfMf, transParams::offsetElectronic, transParams::minSpeedOfSound, transParams::maxSpeedOfSound, distance[cnt], transParams::aScanReconstructionFrequency);
|
||||
Aurora::Matrix rfEnv = Aurora::abs(Aurora::hilbert(rfWin));
|
||||
Aurora::Matrix rfNorm = rfEnv / Aurora::max(rfEnv);
|
||||
tof[cnt] = findFirstvalueGreaterThanGivenValue(rfNorm, thre) + 1;
|
||||
}
|
||||
}
|
||||
}
|
||||
tof = offset + tof/transParams::aScanReconstructionFrequency;
|
||||
Aurora::Matrix sosRough = distance / tof;
|
||||
Aurora::Matrix mu = Aurora::mean(sosRough);
|
||||
Aurora::Matrix sigma = Aurora::std(sosRough);
|
||||
Aurora::Matrix idx_ = Aurora::abs(sosRough - mu) <= sigma;
|
||||
Aurora::Matrix sosClean = removeDataFromArrays(sosRough, idx_);
|
||||
Aurora::Matrix sos = mean(sosClean);
|
||||
RECON_INFO("SOS Value");
|
||||
RECON_INFO(std::to_string(sos[0]));
|
||||
return sos;
|
||||
}
|
||||
@@ -8,7 +8,7 @@ class Parser;
|
||||
|
||||
namespace Recon
|
||||
{
|
||||
Aurora::Matrix estimateSOSWater(Parser* aParser, const CEInfo& aCE);
|
||||
Aurora::Matrix estimateSOSWater(Parser* aParser, const CEInfo& aCE, const PreComputes& aPreComputes);
|
||||
}
|
||||
|
||||
#endif
|
||||
@@ -0,0 +1,11 @@
|
||||
#include "getSoundSpeedWaterValue.h"
|
||||
#include "estimateSOSWater.h"
|
||||
|
||||
using namespace Recon;
|
||||
|
||||
SOSInfo Recon::getSoundSpeedWaterValue(Parser* aParser, const CEInfo& aCE, const PreComputes& aPreComputes)
|
||||
{
|
||||
SOSInfo result;
|
||||
result.expectedSOSWater = Recon::estimateSOSWater(aParser, aCE, aPreComputes);
|
||||
return result;
|
||||
}
|
||||
@@ -0,0 +1,19 @@
|
||||
#ifndef GET_SOUND_SPEED_WATER_VALUE_H
|
||||
#define GET_SOUND_SPEED_WATER_VALUE_H
|
||||
|
||||
#include "Matrix.h"
|
||||
#include "common/getMeasurementMetaData.h"
|
||||
|
||||
class Parser;
|
||||
|
||||
struct SOSInfo
|
||||
{
|
||||
Aurora::Matrix expectedSOSWater;
|
||||
};
|
||||
|
||||
namespace Recon
|
||||
{
|
||||
SOSInfo getSoundSpeedWaterValue(Parser* aParser, const CEInfo& aCE, const PreComputes& aPreComputes);
|
||||
}
|
||||
|
||||
#endif
|
||||
@@ -0,0 +1,29 @@
|
||||
#include "temperatureToSoundSpeed.h"
|
||||
#include "estimateSOSWater.h"
|
||||
|
||||
#include "Function1D.h"
|
||||
|
||||
using namespace Recon;
|
||||
|
||||
Aurora::Matrix Recon::temperatureToSoundSpeed(const Aurora::Matrix& aTemperature, const std::string& aMethod)
|
||||
{
|
||||
Aurora::Matrix result;
|
||||
if (aMethod == "marczak")
|
||||
{
|
||||
float* kData = new float[6] {2.78786e-9, -1.398845e-6, 3.287156e-4, -5.799136e-2, 5.038813, 1.402385e3};
|
||||
Aurora::Matrix k = Aurora::Matrix::fromRawData(kData, 6);
|
||||
result = polyval(k, aTemperature);
|
||||
}
|
||||
else if(aMethod == "mader")
|
||||
{
|
||||
float* kData = new float[6] {3.14643e-9, -1.478e-6, 0.000334199, -0.0580852, 5.03711, 1402.39};
|
||||
Aurora::Matrix k = Aurora::Matrix::fromRawData(kData, 6);
|
||||
result = polyval(k, aTemperature);
|
||||
}
|
||||
else if(aMethod == "jan")
|
||||
{
|
||||
float speedData = 1557 - 0.0245 * pow((74 - aTemperature.getData()[0]), 2);
|
||||
result = Aurora::Matrix::fromRawData(new float[1] {speedData}, 1);
|
||||
}
|
||||
return result;
|
||||
}
|
||||
@@ -0,0 +1,11 @@
|
||||
#ifndef TEMPERATURE_TO_SOUND_SPEED_H
|
||||
#define TEMPERATURE_TO_SOUND_SPEED_H
|
||||
|
||||
#include "Matrix.h"
|
||||
|
||||
namespace Recon
|
||||
{
|
||||
Aurora::Matrix temperatureToSoundSpeed(const Aurora::Matrix& aTemperature, const std::string& aMethod);
|
||||
}
|
||||
|
||||
#endif
|
||||
@@ -17,7 +17,7 @@
|
||||
#include "ceMatchedFilterHandling.h"
|
||||
#include "ShotList/ShotList.h"
|
||||
#include "config/config.h"
|
||||
#include "temperatureCalculation/estimateSOSWater.h"
|
||||
#include "common/SoundSpeedWaterCalculation/getSoundSpeedWaterValue.h"
|
||||
|
||||
|
||||
using namespace Recon;
|
||||
@@ -127,36 +127,13 @@ TransFormInfo Recon::getTransformationMatrix(Parser* aParser, const Matrix& aMot
|
||||
return result;
|
||||
}
|
||||
|
||||
Matrix Recon::temperatureToSoundSpeed(const Matrix& aTemperature, const std::string& aMethod)
|
||||
{
|
||||
Matrix result;
|
||||
if (aMethod == "marczak")
|
||||
{
|
||||
float* kData = new float[6] {2.78786e-9, -1.398845e-6, 3.287156e-4, -5.799136e-2, 5.038813, 1.402385e3};
|
||||
Matrix k = Matrix::fromRawData(kData, 6);
|
||||
result = polyval(k, aTemperature);
|
||||
}
|
||||
else if(aMethod == "mader")
|
||||
{
|
||||
float* kData = new float[6] {3.14643e-9, -1.478e-6, 0.000334199, -0.0580852, 5.03711, 1402.39};
|
||||
Matrix k = Matrix::fromRawData(kData, 6);
|
||||
result = polyval(k, aTemperature);
|
||||
}
|
||||
else if(aMethod == "jan")
|
||||
{
|
||||
float speedData = 1557 - 0.0245 * pow((74 - aTemperature.getData()[0]), 2);
|
||||
result = Matrix::fromRawData(new float[1] {speedData}, 1);
|
||||
}
|
||||
return result;
|
||||
}
|
||||
|
||||
//已验证，完全正确
|
||||
TempInfo Recon::getTemperatureInfo(Parser* aParser, const CEInfo& aCEInfo)
|
||||
{
|
||||
TempInfo result;
|
||||
result.expectedSOSWater = estimateSOSWater(aParser,aCEInfo);
|
||||
return result;
|
||||
}
|
||||
//该版本删除该函数
|
||||
// TempInfo Recon::getTemperatureInfo(Parser* aParser, const CEInfo& aCEInfo, const PreComputes& aPreComputes)
|
||||
// {
|
||||
// TempInfo result;
|
||||
// result.expectedSOSWater = getSoundSpeedWaterValue(aParser,aCEInfo, aPreComputes);
|
||||
// return result;
|
||||
// }
|
||||
|
||||
//已验证 完全正确
|
||||
CEInfo Recon::getCEInfo(Parser* aParser, const MeasurementInfo aInfo)
|
||||
|
||||
@@ -71,9 +71,8 @@ namespace Recon
|
||||
Aurora::Matrix getAvailableMotorPositions(Parser* aParser);
|
||||
MeasurementInfo loadMeasurementInfos(Parser* aParser);
|
||||
TransFormInfo getTransformationMatrix(Parser* aParser, const Aurora::Matrix& aMotorPosList);
|
||||
TempInfo getTemperatureInfo(Parser* aParser, const CEInfo& aCEInfo);
|
||||
//TempInfo getTemperatureInfo(Parser* aParser, const CEInfo& aCEInfo, const PreComputes& aPreComputes);
|
||||
CEInfo getCEInfo(Parser* aParser, const MeasurementInfo aInfo);
|
||||
Aurora::Matrix temperatureToSoundSpeed(const Aurora::Matrix& aTemperature, const std::string& aMethod);
|
||||
|
||||
}
|
||||
#endif // !GET_MEASUREMENT_METADATA_H
|
||||
@@ -1,102 +0,0 @@
|
||||
#include "estimateSOSWater.h"
|
||||
#include "Function1D.cuh"
|
||||
#include "Function1D.h"
|
||||
#include "Matrix.h"
|
||||
#include "CudaMatrix.h"
|
||||
#include "common/dataBlockCreation/getAscanBlock.h"
|
||||
#include "transmissionReconstruction/detection/getTransmissionData.cuh"
|
||||
#include "transmissionReconstruction/detection/detection.h"
|
||||
#include "config/config.h"
|
||||
#include "config/geometryConfig.h"
|
||||
|
||||
#include "Function2D.cuh"
|
||||
#include "Function2D.h"
|
||||
#include "common/dataBlockCreation/getAscanBlock.h"
|
||||
|
||||
#include <cstring>
|
||||
|
||||
using namespace Recon;
|
||||
|
||||
namespace
|
||||
{
|
||||
const Aurora::Matrix EMIT_TAS = Aurora::Matrix::fromRawData(new float[17]{94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110}, 17);
|
||||
const Aurora::Matrix RECEIVE_TAS = Aurora::Matrix::fromRawData(new float[17]{103, 104, 105, 106, 107, 108, 109, 110, 94, 95, 96, 97, 98, 99, 100, 101, 102}, 17);
|
||||
const Aurora::Matrix ALL_RECEIVER_TAS_LIST = Aurora::Matrix::fromRawData(new float[128] {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128},128); //1~128
|
||||
const Aurora::Matrix ALL_RECEIVER_ELEMENT_LIST = Aurora::Matrix::fromRawData(new float[18] {1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18},18);
|
||||
const float SOS_INITIAL = 1509.3;
|
||||
const float THRE = 0.5;
|
||||
inline int findIndexFromEmitterAndReceiver(const Aurora::Matrix& aEmitter, float aEmitterValue, const Aurora::Matrix& aReceiver, float aReceiverValue)
|
||||
{
|
||||
for (int i = 0; i < aEmitter.getDataSize(); ++i)
|
||||
{
|
||||
if (aEmitter[i] == aEmitterValue && aReceiver[i] == aReceiverValue)
|
||||
{
|
||||
return i;
|
||||
}
|
||||
}
|
||||
return -1;
|
||||
}
|
||||
|
||||
Aurora::Matrix matchFilt(const Aurora::CudaMatrix& aSig, const Aurora::CudaMatrix& aMod)
|
||||
{
|
||||
Aurora::CudaMatrix fftMod = fft(aMod);
|
||||
Aurora::CudaMatrix fftSig = fft(aSig);
|
||||
Aurora::CudaMatrix mfFftSig = getTransmissionDataSubFunction(fftSig, fftMod);
|
||||
Aurora::Matrix mfSig = Aurora::real(ifft(mfFftSig)).toHostMatrix();
|
||||
std::memmove(mfSig.getData(), mfSig.getData() + 1, (mfSig.getDataSize() - 1) * sizeof(float));
|
||||
return mfSig;
|
||||
}
|
||||
|
||||
inline int findFirstvalueGreaterThanGivenValue(const Aurora::Matrix& aMatrix, float aValue)
|
||||
{
|
||||
for(int i=0; i<aMatrix.getDataSize(); ++i)
|
||||
{
|
||||
if(aMatrix[i] > aValue)
|
||||
{
|
||||
return i;
|
||||
}
|
||||
}
|
||||
return -1;
|
||||
}
|
||||
}
|
||||
|
||||
Aurora::Matrix Recon::estimateSOSWater(Parser *aParser, const CEInfo &aCE)
|
||||
{
|
||||
int offset = transParams::priorWvOffset;
|
||||
int rn = 18;
|
||||
int emitterSize = EMIT_TAS.getDataSize();
|
||||
Aurora::Matrix distance = Aurora::Matrix::fromRawData(new float[emitterSize]{0}, emitterSize);
|
||||
Aurora::Matrix mp = Aurora::Matrix::fromRawData(new float[1]{1}, 1);
|
||||
Aurora::Matrix sn = Aurora::Matrix::fromRawData(new float[1]{18}, 1);
|
||||
Aurora::Matrix tof = Aurora::Matrix::fromRawData(new float[emitterSize], 1, emitterSize);
|
||||
float snValue = 18;
|
||||
#pragma omp parallel for
|
||||
for(int i=0; i<emitterSize; ++i)
|
||||
{
|
||||
Aurora::Matrix sl = Aurora::Matrix::fromRawData(new float[1]{EMIT_TAS[i]}, 1);
|
||||
float slValue = EMIT_TAS[i];
|
||||
float rl = RECEIVE_TAS[i];
|
||||
Aurora::Matrix positionOfEmitter = emitterPositions[slValue-1](snValue-1).toMatrix();
|
||||
Aurora::Matrix positionOfReceiver = emitterPositions[rl-1](rn-1).toMatrix();
|
||||
distance[i] = Aurora::norm(positionOfEmitter - positionOfReceiver, Aurora::Norm2);
|
||||
AscanBlock ascan = getAscanBlock(aParser, mp, sl, sn, ALL_RECEIVER_TAS_LIST, ALL_RECEIVER_ELEMENT_LIST);
|
||||
size_t dataSize = ascan.ascanBlock.getDataSize();
|
||||
short* data = new short[dataSize];
|
||||
std::copy(ascan.ascanBlock.getData(), ascan.ascanBlock.getData() + dataSize, data);
|
||||
Aurora::Matrix rf = Aurora::convertfp16tofloat(data, ascan.ascanBlock.getDimSize(0), ascan.ascanBlock.getDimSize(1));
|
||||
delete[] data;
|
||||
Aurora::CudaMatrix rfDevice = rf.toDeviceMatrix() / ascan.gainBlock.toDeviceMatrix();
|
||||
int rfIndex = findIndexFromEmitterAndReceiver(aCE.tasIndices, rl, aCE.receiverIndices, rn);
|
||||
rfDevice = rfDevice.block(1, rfIndex, rfIndex);
|
||||
Aurora::Matrix rfMf = matchFilt(rfDevice, aCE.ceRef.toDeviceMatrix());
|
||||
Aurora::Matrix rfWin = applyTimeWindowing(rfMf, transParams::aScanReconstructionFrequency, distance(i).toMatrix(), SOS_INITIAL,
|
||||
transParams::offsetElectronic * transParams::aScanReconstructionFrequency, transParams::detectionWindowSOS,
|
||||
transParams::minSpeedOfSound, transParams::maxSpeedOfSound, transParams::gaussWindow).AscanBlockProcessed;
|
||||
Aurora::Matrix rfEnv = abs(Aurora::hilbert(rfWin));
|
||||
Aurora::Matrix rfNorm = rfEnv / max(rfEnv);
|
||||
tof[i] = findFirstvalueGreaterThanGivenValue(rfNorm, THRE) + 1;
|
||||
}
|
||||
tof = (tof + offset) / transParams::aScanReconstructionFrequency;
|
||||
auto sos = mean(distance / tof);
|
||||
return sos;
|
||||
}
|
||||
@@ -18,7 +18,7 @@ using namespace Recon;
|
||||
|
||||
PreprocessReflectionData Recon::preprocessTransmissionReconstructionForReflection(const Aurora::Matrix& aSOSMap, const Aurora::Matrix& aATTMap,
|
||||
const Aurora::Matrix& aDdims, const GeometryInfo& aGeometryInfo,
|
||||
const TempInfo& aTemp)
|
||||
const SOSInfo& aSOSInfo)
|
||||
{
|
||||
PreprocessReflectionData result;
|
||||
int soundSpeedCorrection = reflectParams::soundSpeedCorrection;
|
||||
@@ -57,7 +57,7 @@ PreprocessReflectionData Recon::preprocessTransmissionReconstructionForReflectio
|
||||
{
|
||||
if(std::abs(sppedMap3d[i]) < 1)
|
||||
{
|
||||
sppedMap3d[i] = aTemp.expectedSOSWater[0];
|
||||
sppedMap3d[i] = aSOSInfo.expectedSOSWater[0];
|
||||
}
|
||||
}
|
||||
result.transRecos.speedMap3D = sppedMap3d;
|
||||
@@ -116,7 +116,7 @@ PreprocessReflectionData Recon::preprocessTransmissionReconstructionForReflectio
|
||||
// transRecos.SpeedMap3D = reshape(temperatureToSoundSpeed(polyvaln(temp.TemperatureModel4D.p,[XNEW(:),YNEW(:), ZNEW(:), zeros(size(XNEW(:)))])),size(XNEW));
|
||||
// catch
|
||||
|
||||
result.transRecos.speedMap3D = zeros(x.getDimSize(1), y.getDimSize(1), z.getDimSize(1)) + aTemp.expectedSOSWater;
|
||||
result.transRecos.speedMap3D = zeros(x.getDimSize(1), y.getDimSize(1), z.getDimSize(1)) + aSOSInfo.expectedSOSWater;
|
||||
|
||||
|
||||
result.transRecos.beginTransMap = begin_TransMap;
|
||||
|
||||
@@ -4,6 +4,7 @@
|
||||
#include "Matrix.h"
|
||||
#include "common/getGeometryInfo.h"
|
||||
#include "common/getMeasurementMetaData.h"
|
||||
#include "common/SoundSpeedWaterCalculation/getSoundSpeedWaterValue.h"
|
||||
|
||||
namespace Recon
|
||||
{
|
||||
@@ -24,7 +25,7 @@ namespace Recon
|
||||
};
|
||||
PreprocessReflectionData preprocessTransmissionReconstructionForReflection(const Aurora::Matrix& aSOSMap, const Aurora::Matrix& aATTMap,
|
||||
const Aurora::Matrix& aDdims, const GeometryInfo& aGeometryInfo,
|
||||
const TempInfo& aTemp);
|
||||
const SOSInfo& aSOSInfo);
|
||||
}
|
||||
|
||||
#endif
|
||||
@@ -20,9 +20,10 @@
|
||||
#include "common/estimatePulseLength.h"
|
||||
#include "common/ceMatchedFilterHandling.h"
|
||||
#include "common/createPositioningImage.h"
|
||||
#include "common/estimateOffset.h"
|
||||
#include "common/SoundSpeedWaterCalculation/getSoundSpeedWaterValue.h"
|
||||
#include "transmissionReconstruction/startTransmissionReconstruction.h"
|
||||
#include "transmissionReconstruction/saveTransmissionImagesInReflCoords.h"
|
||||
#include "reflectionReconstruction/preprocessData/estimateOffset.h"
|
||||
#include "reflectionReconstruction/preprocessData/precalcImageParameters.h"
|
||||
#include "reflectionReconstruction/preprocessData/preprocessTransmissionReconstructionForReflection.h"
|
||||
#include "reflectionReconstruction/startReflectionReconstruction.h"
|
||||
@@ -86,21 +87,38 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
|
||||
float maxNumTAS = max(auroraUnion(slList, rlList)).getData()[0];
|
||||
MeasurementInfo expInfo = loadMeasurementInfos(&dataParser);
|
||||
CEInfo ce = getCEInfo(&dataParser, expInfo);
|
||||
TempInfo temp = getTemperatureInfo(&dataParser, ce);
|
||||
|
||||
PreComputes preComputes;
|
||||
if(ce.measuredCEUsed)
|
||||
{
|
||||
preComputes.matchedFilter = createMatchedFilter(ce.ce, ce.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
|
||||
}
|
||||
else
|
||||
{
|
||||
preComputes.matchedFilter = createMatchedFilter(ce.ceRef, ce.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
|
||||
}
|
||||
preComputes.timeInterval = (float)1 / reflectParams::aScanReconstructionFrequency;
|
||||
preComputes.measuredCEUsed = ce.measuredCEUsed;
|
||||
preComputes.measuredCE_TASIndices = ce.tasIndices;
|
||||
preComputes.measuredCE_receiverIndices = ce.receiverIndices;
|
||||
preComputes.offset = estimateOffset(expInfo, ce);
|
||||
|
||||
SOSInfo sosValue = getSoundSpeedWaterValue(&dataParser, ce, preComputes);
|
||||
SOSInfo sosValueRef;
|
||||
|
||||
TransFormInfo transformationInfo = getTransformationMatrix(&dataParser, motorPosTotal);
|
||||
Matrix transformationMatrices = transformationInfo.rotationMatrix;
|
||||
Matrix motorPosAvailable = transformationInfo.motorPos;
|
||||
|
||||
Matrix motorPosAvailableRef;
|
||||
MeasurementInfo expInfoRef;
|
||||
TempInfo tempRef;
|
||||
CEInfo ceRef;
|
||||
Matrix transformationMatricesRef;
|
||||
if(transParams::runTransmissionReco)
|
||||
{
|
||||
expInfoRef = loadMeasurementInfos(&refParser);
|
||||
ceRef = getCEInfo(&refParser, expInfoRef);
|
||||
tempRef = getTemperatureInfo(&refParser, ceRef);
|
||||
sosValueRef = getSoundSpeedWaterValue(&refParser, ceRef, preComputes);
|
||||
transformationInfo = getTransformationMatrix(&refParser, motorPosTotal);
|
||||
transformationMatricesRef = transformationInfo.rotationMatrix;
|
||||
motorPosAvailableRef = transformationInfo.motorPos;
|
||||
@@ -142,18 +160,6 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
|
||||
}
|
||||
|
||||
GeometryInfo geom = getGeometryInfo(motorPosAvailable, transformationMatrices, rlList, rnList, slList, snList);
|
||||
PreComputes preComputes;
|
||||
if((reflectParams::matchedFilterCeAScan && reflectParams::runReflectionReco) || transParams::runTransmissionReco)
|
||||
{
|
||||
if(ce.measuredCEUsed)
|
||||
{
|
||||
preComputes.matchedFilter = createMatchedFilter(ce.ce, ce.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
|
||||
}
|
||||
else
|
||||
{
|
||||
preComputes.matchedFilter = createMatchedFilter(ce.ceRef, ce.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
|
||||
}
|
||||
}
|
||||
if(expInfo.sampleRate != reflectParams::aScanReconstructionFrequency)
|
||||
{
|
||||
reflectParams::expectedAScanDataLength = ceil(expInfo.numberSamples * ((float)reflectParams::aScanReconstructionFrequency / expInfo.sampleRate));
|
||||
@@ -172,15 +178,8 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
|
||||
{
|
||||
preComputes.matchedFilterRef = createMatchedFilter(ceRef.ceRef, ceRef.measuredCEUsed, reflectParams::findDefects, reconParams::removeOutliersFromCEMeasured, expInfo.Hardware);
|
||||
}
|
||||
preComputes.timeInterval = (float)1 / reflectParams::aScanReconstructionFrequency;
|
||||
preComputes.measuredCEUsed = ce.measuredCEUsed;
|
||||
preComputes.measuredCE_TASIndices = ce.tasIndices;
|
||||
preComputes.measuredCE_receiverIndices = ce.receiverIndices;
|
||||
preComputes.offset = estimateOffset(expInfo, ce);
|
||||
|
||||
expInfo.matchedFilter = preComputes.matchedFilter;
|
||||
expInfoRef.matchedFilter = preComputes.matchedFilterRef;
|
||||
|
||||
CeEstimatePulseLength ceEstimatePulseLength;
|
||||
CeEstimatePulseLength ceEstimatePulseLengthRef;
|
||||
if(ce.measuredCEUsed)
|
||||
@@ -217,7 +216,7 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
|
||||
|
||||
GeometryInfo geomRef = getGeometryInfo(motorPosAvailableRef, transformationMatricesRef, rlList, rnList, slList, snList);
|
||||
RECON_INFO("Start transmissionRecostruction.");
|
||||
transmissionResult = startTransmissionReconstruction(mp_inter, mpRef_inter, slList_inter, snList_inter, rlList_inter, rnList_inter, temp, tempRef, geom, geomRef, expInfo, expInfoRef, preComputes, &dataParser, &refParser);
|
||||
transmissionResult = startTransmissionReconstruction(mp_inter, mpRef_inter, slList_inter, snList_inter, rlList_inter, rnList_inter, sosValue, sosValueRef, geom, geomRef, expInfo, expInfoRef, preComputes, &dataParser, &refParser);
|
||||
attAvailable = true;
|
||||
sosAvailable = true;
|
||||
Matrix minPos = Matrix::fromRawData(new float[3],3,1,1);
|
||||
@@ -238,17 +237,12 @@ int Recon::startReconstructions( const std::string& aDataPath, const std::string
|
||||
Matrix recoATT = transmissionResult.recoATT;
|
||||
//检测可使用内存是否足够,输出警报用,todo
|
||||
//checkEnvAndMemory(reflectParams.imageInfos.IMAGE_XYZ);
|
||||
auto preProcessData = preprocessTransmissionReconstructionForReflection(recoSOS, recoATT, transmissionResult.ddmis, geom, temp);
|
||||
auto preProcessData = preprocessTransmissionReconstructionForReflection(recoSOS, recoATT, transmissionResult.ddmis, geom, sosValue);
|
||||
Matrix mp_inter = intersect(motorPosAvailable, reflectParams::motorPos);
|
||||
Matrix slList_inter = intersect(slList, reflectParams::senderTasList);
|
||||
Matrix snList_inter = intersect(snList, reflectParams::senderElementList);
|
||||
Matrix rlList_inter = intersect(rlList, reflectParams::receiverTasList);
|
||||
Matrix rnList_inter = intersect(rnList, reflectParams::receiverElementList);
|
||||
preComputes.timeInterval = (float)1 / reflectParams::aScanReconstructionFrequency;
|
||||
preComputes.measuredCEUsed = ce.measuredCEUsed;
|
||||
preComputes.measuredCE_TASIndices = ce.tasIndices;
|
||||
preComputes.measuredCE_receiverIndices = ce.receiverIndices;
|
||||
preComputes.offset = estimateOffset(expInfo, ce);
|
||||
|
||||
reflectParams::gpuSelectionList = reconParams::gpuSelectionList;
|
||||
RECON_INFO("Start reflectionRecostruction.");
|
||||
|
||||
@@ -82,11 +82,12 @@ namespace Recon {
|
||||
}
|
||||
|
||||
TimeWindowResult applyTimeWindowing(const Aurora::Matrix &AscanBlock, float sampleRate,
|
||||
const Aurora::Matrix &distBlock, float sosWater,
|
||||
const Aurora::Matrix &distBlock, const Aurora::Matrix& sosBlock,float sosWater,
|
||||
float startOffset, float segmentLenOffset,
|
||||
float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow)
|
||||
{
|
||||
Aurora::Matrix sosOffset = Aurora::Matrix::fromRawData(new float[1]{0}, 1);
|
||||
// Aurora::Matrix sosOffset = Aurora::Matrix::fromRawData(new float[1]{0}, 1);
|
||||
Aurora::Matrix sosOffset = calculateSOSOffset(sosBlock, sosWater, distBlock, sampleRate);
|
||||
auto calcResult = calculateStarEndSearchPosition(distBlock, minSpeedOfSound, maxSpeedOfSound, sampleRate,AscanBlock.getDimSize(0), sosOffset, startOffset, segmentLenOffset);
|
||||
|
||||
auto AscanBlockProcessed = zeros(AscanBlock.getDimSize(0),AscanBlock.getDimSize(1));
|
||||
@@ -129,7 +130,7 @@ namespace Recon {
|
||||
result.AscanBlockProcessed = AscanBlockProcessed;
|
||||
return result;
|
||||
}
|
||||
|
||||
|
||||
Aurora::Matrix detectAttVectorized(const Aurora::Matrix &Ascan, const Aurora::Matrix &AscanRef,
|
||||
const Aurora::Matrix &distRef,
|
||||
float sosWaterRef,
|
||||
@@ -173,6 +174,7 @@ namespace Recon {
|
||||
DetectResult detectTofVectorized(
|
||||
const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
|
||||
const Aurora::Matrix &distBlock, const Aurora::Matrix &distBlockRef,
|
||||
const Aurora::Matrix &sosWaterBlock,const Aurora::Matrix &sosWaterRefBlock,
|
||||
float aSOSWater, float aSOSWaterRef,
|
||||
int useTimeWindowing, int aScanReconstructionFrequency,
|
||||
float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound,
|
||||
@@ -187,11 +189,11 @@ namespace Recon {
|
||||
timeResult2.AscanBlockProcessed = AscanRefBlock;
|
||||
if (useTimeWindowing == 1) {
|
||||
timeResult1 = applyTimeWindowing(
|
||||
AscanBlock, sampleRate, distBlock,
|
||||
AscanBlock, sampleRate, distBlock, sosWaterBlock,
|
||||
aSOSWater, offsetElectronicSamples, detectionWindowSOS,
|
||||
minSpeedOfSound, maxSpeedOfSound, gaussWindow);
|
||||
timeResult2 = applyTimeWindowing(
|
||||
AscanRefBlock, sampleRate, distBlockRef,
|
||||
AscanRefBlock, sampleRate, distBlockRef, sosWaterRefBlock,
|
||||
aSOSWaterRef, offsetElectronicSamples, detectionWindowSOS,
|
||||
minSpeedOfSound, maxSpeedOfSound, gaussWindow);
|
||||
|
||||
@@ -244,6 +246,7 @@ namespace Recon {
|
||||
DetectResult detectTofAndAtt(
|
||||
const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
|
||||
const Aurora::Matrix &distBlock, const Aurora::Matrix &distRefBlock,
|
||||
const Aurora::Matrix &sosWaterBlock,const Aurora::Matrix &sosWaterRefBlock,
|
||||
int resampleFactor,int nthreads, float aSOSWater, float aSOSWaterRef,
|
||||
int useTimeWindowing, int aScanReconstructionFrequency,int detectionWindowATT,
|
||||
float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound,
|
||||
@@ -258,11 +261,11 @@ namespace Recon {
|
||||
timeResult2.AscanBlockProcessed = AscanRefBlock;
|
||||
if (useTimeWindowing == 1) {
|
||||
timeResult1 = applyTimeWindowing(
|
||||
AscanBlock, sampleRate, distBlock,
|
||||
AscanBlock, sampleRate, distBlock, sosWaterBlock,
|
||||
aSOSWater, offsetElectronicSamples, detectionWindowSOS,
|
||||
minSpeedOfSound, maxSpeedOfSound, gaussWindow);
|
||||
timeResult2 = applyTimeWindowing(
|
||||
AscanRefBlock, sampleRate, distRefBlock,
|
||||
AscanRefBlock, sampleRate, distRefBlock, sosWaterRefBlock,
|
||||
aSOSWaterRef, offsetElectronicSamples, detectionWindowSOS,
|
||||
minSpeedOfSound, maxSpeedOfSound, gaussWindow);
|
||||
|
||||
@@ -420,7 +423,7 @@ namespace Recon {
|
||||
const Aurora::CudaMatrix &AscanRefBlock,
|
||||
const Aurora::CudaMatrix &distBlock,
|
||||
const Aurora::CudaMatrix &distRefBlock,
|
||||
float aSOSWater, float aSOSWaterRef)
|
||||
float aSOSWater, float aSOSWaterRef, float aOffset)
|
||||
{
|
||||
switch (Recon::transParams::version) {
|
||||
// case 1: {
|
||||
@@ -433,19 +436,40 @@ namespace Recon {
|
||||
// Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow);
|
||||
// }
|
||||
// case 2:
|
||||
// {
|
||||
// auto SOSWaterBlock = zerosCuda(1,AscanBlock.getDimSize(1));
|
||||
// SOSWaterBlock.setBlockValue(1,0,AscanBlock.getDimSize(1)-1,aSOSWater);
|
||||
// auto SOSWaterRefBlock = zerosCuda(1,AscanBlock.getDimSize(1));
|
||||
// SOSWaterRefBlock.setBlockValue(1,0,AscanBlock.getDimSize(1)-1,aSOSWaterRef);
|
||||
// auto r = detectTofAndAtt(
|
||||
// AscanBlock, AscanRefBlock, distBlock, distRefBlock,SOSWaterBlock,SOSWaterRefBlock,
|
||||
// Recon::transParams::resampleFactor, Recon::transParams::nThreads,
|
||||
// aSOSWater, aSOSWaterRef, Recon::transParams::useTimeWindowing,
|
||||
// Recon::transParams::aScanReconstructionFrequency, Recon::transParams::detectionWindowATT,
|
||||
// Recon::transParams::offsetElectronic, Recon::transParams::detectionWindowSOS, Recon::transParams::minSpeedOfSound,
|
||||
// Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow);
|
||||
// DetectResult ret;
|
||||
// ret.att = r.att.toHostMatrix();
|
||||
// ret.sosValue = r.sosValue.toHostMatrix();
|
||||
// ret.tof = r.tof.toHostMatrix();
|
||||
// return ret;
|
||||
// }
|
||||
// case 3:
|
||||
default:
|
||||
auto r = detectTofAndAtt(
|
||||
AscanBlock, AscanRefBlock, distBlock, distRefBlock,
|
||||
Recon::transParams::resampleFactor, Recon::transParams::nThreads,
|
||||
aSOSWater, aSOSWaterRef, Recon::transParams::useTimeWindowing,
|
||||
Recon::transParams::aScanReconstructionFrequency, Recon::transParams::detectionWindowATT,
|
||||
Recon::transParams::offsetElectronic, Recon::transParams::detectionWindowSOS, Recon::transParams::minSpeedOfSound,
|
||||
Recon::transParams::maxSpeedOfSound, Recon::transParams::gaussWindow);
|
||||
{
|
||||
float sampleRate = Recon::transParams::aScanReconstructionFrequency;
|
||||
auto tof = detectTofThreshold(AscanBlock, distBlock, aOffset, aSOSWater, sampleRate,
|
||||
Recon::transParams::offsetElectronic*sampleRate,
|
||||
Recon::transParams::gaussWindow, Recon::transParams::minSpeedOfSound,
|
||||
Recon::transParams::maxSpeedOfSound);
|
||||
auto att = detectAttVectorizedCuda(AscanBlock, AscanRefBlock,distRefBlock, aSOSWaterRef, tof,
|
||||
sampleRate, Recon::transParams::offsetElectronic,
|
||||
Recon::transParams::detectionWindowATT);
|
||||
DetectResult ret;
|
||||
ret.att = r.att.toHostMatrix();
|
||||
ret.sosValue = r.sosValue.toHostMatrix();
|
||||
ret.tof = r.tof.toHostMatrix();
|
||||
ret.tof = (tof - distBlock/aSOSWater).toHostMatrix();
|
||||
ret.att = att.toHostMatrix();
|
||||
return ret;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
@@ -99,6 +99,115 @@ CudaMatrix Recon::calculateSOSOffset(const CudaMatrix &sosBlock, float reference
|
||||
return offset;
|
||||
}
|
||||
|
||||
__global__ void addWinKernel(float* aData, int aRowCount, int aColCount, float aOffset, float aMinSOS, float aMaxSOS, float* aDist, float aSamepleRate)
|
||||
{
|
||||
__shared__ int startIndex;
|
||||
__shared__ int endIndex;
|
||||
unsigned int i = blockIdx.x;
|
||||
unsigned int base_offset = i * aRowCount;
|
||||
float * dataPointer = aData + base_offset;
|
||||
float dist = aDist[i];
|
||||
//确定每个thread的index
|
||||
if (threadIdx.x == 0)
|
||||
{
|
||||
float offset = aOffset * aSamepleRate;
|
||||
startIndex = (int)(round(dist / aMaxSOS * aSamepleRate) + offset)-1;
|
||||
endIndex = (int)(round(dist / aMinSOS * aSamepleRate) + offset)-1;
|
||||
}
|
||||
for (size_t j = threadIdx.x; j < aRowCount; j+=blockDim.x)
|
||||
{
|
||||
if(j<startIndex) dataPointer[j] = .0;
|
||||
if(j>endIndex) dataPointer[j] = .0;
|
||||
}
|
||||
}
|
||||
|
||||
__global__ void genTofKernel(float* aData, int aRowCount, int aColCount, float aOffset, float* aTOFRef, float aSamepleRate, int* index, int aNpeaks)
|
||||
{
|
||||
unsigned int i = blockIdx.x;
|
||||
unsigned int base_offset = i * aRowCount ;
|
||||
float * dataPointer = aData + base_offset;
|
||||
int startIndex = index[i*aNpeaks]-15;
|
||||
startIndex = startIndex < 0 ? 0 : startIndex;
|
||||
int endIndex = index[i*aNpeaks]+15;
|
||||
endIndex = endIndex >= aRowCount? aRowCount-1 : endIndex;
|
||||
int maxIndex = 0;
|
||||
float maxValue = FLT_MIN;
|
||||
//find max in index[i] + (-15,15)
|
||||
for (size_t j = startIndex; j <= endIndex; j++)
|
||||
{
|
||||
if (dataPointer[j] > maxValue){
|
||||
maxIndex = j;
|
||||
maxValue = dataPointer[j];
|
||||
}
|
||||
}
|
||||
index[i] = maxIndex;
|
||||
//attention: this remove redudent -aTOFRef[i] ,because it will be add back in upper called level
|
||||
aTOFRef[i] = ((float)(maxIndex+1))/aSamepleRate + aOffset ;
|
||||
// aTOFRef[i] = ((float)(maxIndex+1))/aSamepleRate + aOffset - aTOFRef[i];
|
||||
|
||||
}
|
||||
|
||||
CudaMatrix Recon::detectTofThreshold(const CudaMatrix& aAscan, const CudaMatrix &aDistRef, float aOffset, float aSOSWater,
|
||||
float aAScanReconstructionFrequency, float aOffsetElectronic, int aUseTimeWindowing, float aMinSpeedOfSound,
|
||||
float aMaxSpeedOfSound)
|
||||
{
|
||||
const int THREADS_PER_BLOCK_ = 32;
|
||||
float offsetWV = -aOffset;
|
||||
float sampleRate = aAScanReconstructionFrequency;
|
||||
float offsetElectronSamples = aOffsetElectronic * sampleRate;
|
||||
// if (inits.detection.useTimeWindowing == 1)
|
||||
// parfor idx__ = 1 : size(AscanBlock, 2)
|
||||
// AscanBlock(:, idx__) = add_win(AscanBlock(:, idx__) ,...
|
||||
// offsetElectronicSamples,...
|
||||
// inits.detection.minSpeedOfSound,...
|
||||
// inits.detection.maxSpeedOfSound,...
|
||||
// distBlock(idx__), sampleRate);
|
||||
// end
|
||||
// end
|
||||
int * index = nullptr;
|
||||
if (aUseTimeWindowing){
|
||||
addWinKernel<<<aAscan.getDimSize(1), THREADS_PER_BLOCK_>>>(aAscan.getData(), aAscan.getDimSize(0),
|
||||
aAscan.getDimSize(1), offsetElectronSamples,
|
||||
aMinSpeedOfSound, aMaxSpeedOfSound, aDistRef.getData(), sampleRate);
|
||||
cudaDeviceSynchronize();
|
||||
|
||||
}
|
||||
// AScan_env = envelope(AscanBlock);
|
||||
// pkValue_mea = max(AScan_env, [], 1);
|
||||
// AScan_env_norm = AScan_env ./ pkValue_mea;
|
||||
// parfor idx__ = 1 : size(AScan_env_norm, 2)
|
||||
// %-find 1st wvfrom based on envelope
|
||||
// [~, locs] = findpeaks(AScan_env_norm(:, idx__),...
|
||||
// 'MinPeakHeight', 0.2,...
|
||||
// 'Npeaks', 10,...
|
||||
// 'MInPeakProminence', 0.05);
|
||||
// tof_env = locs(1);
|
||||
// %-then select maximum pt on match-filted signal
|
||||
// Ascan_valid_part = zeros(size(AscanBlock, 1), 1);
|
||||
// Ascan_valid_part((-15:15) + tof_env) = AscanBlock((-15:15) + tof_env, idx__);
|
||||
// loc = find(Ascan_valid_part == max(Ascan_valid_part), 1);
|
||||
// try
|
||||
// tof_mea(idx__) = loc;
|
||||
// catch ME
|
||||
// error(ME)
|
||||
// end
|
||||
// end
|
||||
{
|
||||
auto envAScanTemp = abs(hilbert(aAscan));//envelope(Ascanblock);
|
||||
auto pkValueMea = max(envAScanTemp);
|
||||
envAScanTemp = envAScanTemp / pkValueMea;
|
||||
auto pearks = Aurora::findPeaks(envAScanTemp, 10, 0.2, 0.05, &index);
|
||||
}
|
||||
|
||||
auto tofRef = aDistRef.deepCopy();
|
||||
genTofKernel<<<aAscan.getDimSize(1),1>>>(aAscan.getData(), aAscan.getDimSize(0),
|
||||
aAscan.getDimSize(1), offsetWV, tofRef.getData(),
|
||||
sampleRate,index, 10);
|
||||
cudaDeviceSynchronize();
|
||||
cudaFree(index);
|
||||
return tofRef;
|
||||
}
|
||||
|
||||
Recon::SearchPositionC Recon::calculateStarEndSearchPosition(const CudaMatrix &aVDistBlock,
|
||||
float minSpeedOfSound, float maxSpeedOfSound,
|
||||
float sampleRate, float maxSample,
|
||||
@@ -170,11 +279,13 @@ __global__ void guassWindowKernel(float* aStartSearch,float* aEndSearch,
|
||||
}
|
||||
|
||||
Recon::TimeWindowResultC Recon::applyTimeWindowing(const Aurora::CudaMatrix &AscanBlock, float sampleRate,
|
||||
const Aurora::CudaMatrix &distBlock,
|
||||
const Aurora::CudaMatrix &distBlock,const Aurora::CudaMatrix &sosWaterBlock,
|
||||
float aSOSWater, float startOffset, float segmentLenOffset,
|
||||
float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow)
|
||||
{
|
||||
Aurora::CudaMatrix sosOffset = Aurora::zerosCuda(1,1);
|
||||
// Aurora::CudaMatrix sosOffset = Aurora::zerosCuda(1,1);
|
||||
Aurora::CudaMatrix sosOffset = calculateSOSOffset(sosWaterBlock,aSOSWater, distBlock, sampleRate);
|
||||
|
||||
auto calcResult = calculateStarEndSearchPosition(distBlock, minSpeedOfSound, maxSpeedOfSound, sampleRate,AscanBlock.getDimSize(0), sosOffset, startOffset, segmentLenOffset);
|
||||
|
||||
auto AscanBlockProcessed = zerosCuda(AscanBlock.getDimSize(0),AscanBlock.getDimSize(1));
|
||||
@@ -253,6 +364,7 @@ int nextpow2(unsigned int value){
|
||||
Recon::DetectResultC Recon::detectTofAndAtt(
|
||||
const Aurora::CudaMatrix &AscanBlock, const Aurora::CudaMatrix &AscanRefBlock,
|
||||
const Aurora::CudaMatrix &distBlock, const Aurora::CudaMatrix &distRefBlock,
|
||||
const Aurora::CudaMatrix &soswaterBlock,const Aurora::CudaMatrix &soswaterRefBlock,
|
||||
int resampleFactor,int nthreads, float aSOSWater, float aSOSWaterRef,
|
||||
int useTimeWindowing, int aScanReconstructionFrequency,int detectionWindowATT,
|
||||
float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound,
|
||||
@@ -267,11 +379,11 @@ Recon::DetectResultC Recon::detectTofAndAtt(
|
||||
timeResult2.AscanBlockProcessed = AscanRefBlock;
|
||||
if (useTimeWindowing == 1) {
|
||||
timeResult1 = applyTimeWindowing(
|
||||
AscanBlock, sampleRate, distBlock,
|
||||
AscanBlock, sampleRate, distBlock, soswaterBlock,
|
||||
aSOSWater, offsetElectronicSamples, detectionWindowSOS,
|
||||
minSpeedOfSound, maxSpeedOfSound, gaussWindow);
|
||||
timeResult2 = applyTimeWindowing(
|
||||
AscanRefBlock, sampleRate, distRefBlock,
|
||||
AscanRefBlock, sampleRate, distRefBlock, soswaterRefBlock,
|
||||
aSOSWaterRef, offsetElectronicSamples, detectionWindowSOS,
|
||||
minSpeedOfSound, maxSpeedOfSound, gaussWindow);
|
||||
|
||||
|
||||
@@ -30,7 +30,7 @@ struct TimeWindowResultC {
|
||||
float offsetElectronic, int detectionWindowATT);
|
||||
|
||||
TimeWindowResultC applyTimeWindowing(const Aurora::CudaMatrix &AscanBlock, float sampleRate,
|
||||
const Aurora::CudaMatrix &distBlock,
|
||||
const Aurora::CudaMatrix &distBlock,const Aurora::CudaMatrix &sosWaterBlock,
|
||||
float aSOSWater, float startOffset, float segmentLenOffset,
|
||||
float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow);
|
||||
SearchPositionC calculateStarEndSearchPosition(const CudaMatrix &aVDistBlock,
|
||||
@@ -39,9 +39,14 @@ struct TimeWindowResultC {
|
||||
const CudaMatrix &aVSosOffsetBlock,
|
||||
float startOffset, float segmentLenOffset);
|
||||
CudaMatrix calculateSOSOffset(const CudaMatrix &sosBlock, float referenceSOS, const CudaMatrix &distBlock, float sampleRate);
|
||||
|
||||
CudaMatrix detectTofThreshold(const CudaMatrix& aAscan, const CudaMatrix &aDistRef,float aOffset, float aSOSWater,
|
||||
float aAScanReconstructionFrequency, float aOffsetElectronic, int aUseTimeWindowing, float aMinSpeedOfSound, float aMaxSpeedOfSound);
|
||||
|
||||
DetectResultC detectTofAndAtt(
|
||||
const Aurora::CudaMatrix &AscanBlock, const Aurora::CudaMatrix &AscanRefBlock,
|
||||
const Aurora::CudaMatrix &distBlock, const Aurora::CudaMatrix &distRefBlock,
|
||||
const Aurora::CudaMatrix &soswaterBlock, const Aurora::CudaMatrix &soswaterRefBlock,
|
||||
int resampleFactor,int nthreads, float aSOSWater, float aSOSWaterRef,
|
||||
int useTimeWindowing, int aScanReconstructionFrequency,int detectionWindowATT,
|
||||
float offsetElectronic, int detectionWindowSOS, float minSpeedOfSound,
|
||||
|
||||
@@ -32,8 +32,8 @@ calculateStarEndSearchPosition(const Aurora::Matrix &aVDistBlock,
|
||||
|
||||
TimeWindowResult applyTimeWindowing(
|
||||
const Aurora::Matrix &AscanBlock, float sampleRate,
|
||||
const Aurora::Matrix &distBlock, float sosWater,
|
||||
float startOffset, float segmentLenOffset,
|
||||
const Aurora::Matrix &distBlock, const Aurora::Matrix &sosWaterBlock,
|
||||
float sosWater, float startOffset, float segmentLenOffset,
|
||||
float minSpeedOfSound, float maxSpeedOfSound, bool gaussWindow);
|
||||
|
||||
Aurora::Matrix
|
||||
@@ -56,6 +56,7 @@ DetectResult
|
||||
detectTofAndAtt(
|
||||
const Aurora::Matrix &AscanBlock, const Aurora::Matrix &AscanRefBlock,
|
||||
const Aurora::Matrix &distBlock, const Aurora::Matrix &distRefBlock,
|
||||
const Aurora::Matrix &sosWaterBlock,const Aurora::Matrix &sosWaterRefBlock,
|
||||
int resampleFactor, int nthreads, float aSOSWater, float aSOSWaterRef,
|
||||
int useTimeWindowing, int aScanReconstructionFrequency,
|
||||
int detectionWindowATT, float offsetElectronic, int detectionWindowSOS,
|
||||
@@ -77,7 +78,7 @@ DetectResult
|
||||
transmissionDetection(
|
||||
const Aurora::CudaMatrix &AscanBlock, const Aurora::CudaMatrix &AscanRefBlock,
|
||||
const Aurora::CudaMatrix &distBlock, const Aurora::CudaMatrix &distRefBlock,
|
||||
float aSOSWater, float aSOSWaterRef);
|
||||
float aSOSWater, float aSOSWaterRef, float aOffset);
|
||||
|
||||
} // namespace Recon
|
||||
|
||||
|
||||
@@ -78,7 +78,7 @@ namespace
|
||||
}
|
||||
|
||||
BlockOfTransmissionData getBlockOfTransmissionData(const Matrix& aMp, const Matrix& aMpRef, const Matrix& aSl, const Matrix& aSn, const Matrix& aRlList, const Matrix& aRnList,
|
||||
const TempInfo& aTemp, const TempInfo& aTempRef, GeometryInfo aGeom, GeometryInfo& aGeomRef,
|
||||
const SOSInfo& aSOS, const SOSInfo& aSOSRef, GeometryInfo aGeom, GeometryInfo& aGeomRef,
|
||||
const Matrix& aSnrRmsNoise, const Matrix& aSnrRmsNoiseRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
|
||||
const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
|
||||
{
|
||||
@@ -162,13 +162,13 @@ namespace
|
||||
}
|
||||
|
||||
void getBlockOfTransmissionDataInThread(size_t aIndex, const Matrix& aMp, const Matrix& aMpRef, const Matrix& aSl, const Matrix& aSn, const Matrix& aRlList, const Matrix& aRnList,
|
||||
const TempInfo& aTemp, const TempInfo& aTempRef, GeometryInfo aGeom, GeometryInfo aGeomRef,
|
||||
const SOSInfo& aSOS, const SOSInfo& aSOSRef, GeometryInfo aGeom, GeometryInfo aGeomRef,
|
||||
const Matrix& aSnrRmsNoise, const Matrix& aSnrRmsNoiseRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
|
||||
const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef, unsigned int aGPUId)
|
||||
{
|
||||
cudaSetDevice(aGPUId);
|
||||
auto buffer = getBlockOfTransmissionData(aMp, aMpRef, aSl, aSn, aRlList, aRnList, aTemp,
|
||||
aTempRef, aGeom, aGeomRef, aSnrRmsNoise, aSnrRmsNoiseRef,
|
||||
auto buffer = getBlockOfTransmissionData(aMp, aMpRef, aSl, aSn, aRlList, aRnList, aSOS,
|
||||
aSOSRef, aGeom, aGeomRef, aSnrRmsNoise, aSnrRmsNoiseRef,
|
||||
aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef);
|
||||
std::unique_lock<std::mutex> lock(CREATE_BUFFER_MUTEX);
|
||||
BLOCK_OF_TRANSIMISSIONDARA_BUFFER[std::to_string(aIndex)] = buffer;
|
||||
@@ -178,7 +178,7 @@ void getBlockOfTransmissionDataInThread(size_t aIndex, const Matrix& aMp, const
|
||||
}
|
||||
|
||||
void createThreadForGetBlockOfTransmissionData(const Matrix& aMotorPos, const Matrix& aMotoPosRef, const Matrix& aSlList, const Matrix& aSnList, const Matrix& aRlList, const Matrix& aRnList,
|
||||
const TempInfo& aTemp, const TempInfo& aTempRef, GeometryInfo aGeom, GeometryInfo aGeomRef,
|
||||
const SOSInfo& aSOS, const SOSInfo& aSOSRef, GeometryInfo aGeom, GeometryInfo aGeomRef,
|
||||
const Matrix& aSnrRmsNoise, const Matrix& aSnrRmsNoiseRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
|
||||
const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
|
||||
{
|
||||
@@ -202,7 +202,7 @@ void createThreadForGetBlockOfTransmissionData(const Matrix& aMotorPos, const M
|
||||
CREATE_BUFFER_CONDITION.wait(lock, []{return BUFFER_COUNT<BUFFER_SIZE;});
|
||||
++BUFFER_COUNT;
|
||||
lock.unlock();
|
||||
speedUpThread[index] = std::thread(getBlockOfTransmissionDataInThread,index,mp,mpRef,sl,sn,aRlList,aRnList,aTemp,aTempRef,aGeom,aGeomRef,aSnrRmsNoise,aSnrRmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef, index % BUFFER_SIZE);
|
||||
speedUpThread[index] = std::thread(getBlockOfTransmissionDataInThread,index,mp,mpRef,sl,sn,aRlList,aRnList,aSOS,aSOSRef,aGeom,aGeomRef,aSnrRmsNoise,aSnrRmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef, index % BUFFER_SIZE);
|
||||
}
|
||||
}
|
||||
}
|
||||
@@ -216,7 +216,7 @@ void createThreadForGetBlockOfTransmissionData(const Matrix& aMotorPos, const M
|
||||
|
||||
TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, const Aurora::Matrix& aMotoPosRef, const Aurora::Matrix& aSlList,
|
||||
const Aurora::Matrix& aSnList, const Aurora::Matrix& aRlList, const Aurora::Matrix& aRnList,
|
||||
const TempInfo& aTemp, const TempInfo& aTempRef, GeometryInfo& aGeom,
|
||||
const SOSInfo& aSOS, const SOSInfo& aSOSRef, GeometryInfo& aGeom,
|
||||
GeometryInfo& aGeomRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
|
||||
const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
|
||||
{
|
||||
@@ -259,7 +259,7 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
|
||||
rnBlockTotal = zeros(1,numScans,1);
|
||||
}
|
||||
|
||||
std::thread speedUpThread = std::thread(createThreadForGetBlockOfTransmissionData,aMotorPos,aMotoPosRef,aSlList,aSnList,aRlList,aRnList,aTemp,aTempRef,aGeom,aGeomRef,rmsNoise,rmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef);
|
||||
std::thread speedUpThread = std::thread(createThreadForGetBlockOfTransmissionData,aMotorPos,aMotoPosRef,aSlList,aSnList,aRlList,aRnList,aSOS,aSOSRef,aGeom,aGeomRef,rmsNoise,rmsNoiseRef,aExpInfo,aExpInfoRef,aPreComputes,aParser, aParserRef);
|
||||
int numData = 0;
|
||||
int numPossibleScans = 0;
|
||||
sched_param sch;
|
||||
@@ -288,7 +288,7 @@ TransmissionData Recon::getTransmissionData(const Aurora::Matrix& aMotorPos, con
|
||||
cudaSetDevice(index % BUFFER_SIZE);
|
||||
DetectResult detect = transmissionDetection( blockData.ascanBlock, blockData.ascanBlockRef,
|
||||
blockData.dists.toDeviceMatrix(), blockData.distRefBlock.toDeviceMatrix(),
|
||||
aTemp.expectedSOSWater[0], aTempRef.expectedSOSWater[0]);
|
||||
aSOS.expectedSOSWater[0], aSOSRef.expectedSOSWater[0], aPreComputes.offset);
|
||||
blockData.attData = detect.att;
|
||||
blockData.tofData = detect.tof;
|
||||
BlockOfTransmissionData transmissionBlock=blockData;
|
||||
|
||||
@@ -4,6 +4,7 @@
|
||||
#include "Matrix.h"
|
||||
#include "common/getMeasurementMetaData.h"
|
||||
#include "common/getGeometryInfo.h"
|
||||
#include "common/SoundSpeedWaterCalculation/getSoundSpeedWaterValue.h"
|
||||
|
||||
class Parser;
|
||||
namespace Recon
|
||||
@@ -42,7 +43,7 @@ namespace Recon
|
||||
|
||||
TransmissionData getTransmissionData(const Aurora::Matrix& aMotorPos, const Aurora::Matrix& aMotoPosRef, const Aurora::Matrix& aSlList,
|
||||
const Aurora::Matrix& aSnList, const Aurora::Matrix& aRlList, const Aurora::Matrix& aRnList,
|
||||
const TempInfo& aTemp, const TempInfo& aTempRef, GeometryInfo& aGeom,
|
||||
const SOSInfo& aSOS, const SOSInfo& aSOSRef, GeometryInfo& aGeom,
|
||||
GeometryInfo& aGeomRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
|
||||
const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef);
|
||||
}
|
||||
|
||||
@@ -19,16 +19,16 @@ using namespace Recon;
|
||||
|
||||
TransmissionReconstructionResult Recon::startTransmissionReconstruction(const Aurora::Matrix& aMotorPos, const Aurora::Matrix& aMotoPosRef, const Aurora::Matrix& aSlList,
|
||||
const Aurora::Matrix& aSnList, const Aurora::Matrix& aRlList, const Aurora::Matrix& aRnList,
|
||||
const TempInfo& aTemp, const TempInfo& aTempRef, Recon::GeometryInfo& aGeom,
|
||||
const SOSInfo& aSOS, const SOSInfo& aSOSRef, Recon::GeometryInfo& aGeom,
|
||||
GeometryInfo& aGeomRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
|
||||
const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef)
|
||||
{
|
||||
RECON_INFO("Start getTransmissionData.");
|
||||
auto transmissionData = getTransmissionData(aMotorPos, aMotoPosRef, aSlList, aSnList, aRlList, aRnList, aTemp, aTempRef,
|
||||
auto transmissionData = getTransmissionData(aMotorPos, aMotoPosRef, aSlList, aSnList, aRlList, aRnList, aSOS, aSOSRef,
|
||||
aGeom, aGeomRef, aExpInfo, aExpInfoRef, aPreComputes, aParser, aParserRef);
|
||||
Matrix dists = Recon::distanceBetweenTwoPoints(transmissionData.senderList, transmissionData.receiverList);
|
||||
|
||||
Matrix valid = Recon::checkTofDetections(transmissionData.tofDataTotal, dists, aTemp.expectedSOSWater[0],
|
||||
Matrix valid = Recon::checkTofDetections(transmissionData.tofDataTotal, dists, aSOS.expectedSOSWater[0],
|
||||
Recon::transParams::minSpeedOfSound,Recon::transParams::maxSpeedOfSound).valid;
|
||||
if(transParams::qualityCheck)
|
||||
{
|
||||
@@ -40,7 +40,7 @@ TransmissionReconstructionResult Recon::startTransmissionReconstruction(const Au
|
||||
Matrix senderList = removeDataFromArrays(positionValues.senderCoordList, valid);
|
||||
Matrix reveiverList = removeDataFromArrays(positionValues.receiverCoordList, valid);
|
||||
RECON_INFO("Start reconstructArt.");
|
||||
auto transmissionReon = reconstructArt(tofData, attData, positionValues.dims, senderList, reveiverList, positionValues.res, aTemp.expectedSOSWater[0]);
|
||||
auto transmissionReon = reconstructArt(tofData, attData, positionValues.dims, senderList, reveiverList, positionValues.res, aSOS.expectedSOSWater[0]);
|
||||
|
||||
TransmissionReconstructionResult result;
|
||||
result.recoATT = transmissionReon.outATT;
|
||||
|
||||
@@ -5,6 +5,7 @@
|
||||
#include "common/getMeasurementMetaData.h"
|
||||
#include "common/getGeometryInfo.h"
|
||||
#include "transmissionReconstruction/detection/getTransmissionData.h"
|
||||
#include "common/SoundSpeedWaterCalculation/getSoundSpeedWaterValue.h"
|
||||
|
||||
class Parser;
|
||||
|
||||
@@ -19,7 +20,7 @@ namespace Recon
|
||||
|
||||
TransmissionReconstructionResult startTransmissionReconstruction(const Aurora::Matrix& aMotorPos, const Aurora::Matrix& aMotoPosRef, const Aurora::Matrix& aSlList,
|
||||
const Aurora::Matrix& aSnList, const Aurora::Matrix& aRlList, const Aurora::Matrix& aRnList,
|
||||
const TempInfo& aTemp, const TempInfo& aTempRef, Recon::GeometryInfo& aGeom,
|
||||
const SOSInfo& aTemp, const SOSInfo& aTempRef, Recon::GeometryInfo& aGeom,
|
||||
GeometryInfo& aGeomRef, const MeasurementInfo& aExpInfo, const MeasurementInfo& aExpInfoRef,
|
||||
const PreComputes& aPreComputes, Parser* aParser, Parser* aParserRef);
|
||||
}
|
||||
|
||||
Reference in New Issue
Block a user