sunwen/UR
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

feat: new recontest logic

Browse Source
This commit is contained in:
kradchen
2025-03-27 09:06:50 +08:00
parent ed10e822ec
commit 4b878976bb
1 changed files with 140 additions and 23 deletions
Download Patch File Download Diff File Expand all files Collapse all files

163
test/Reconstruction_Test.cpp
View File

@@ -2,6 +2,7 @@
#include <gtest/gtest.h>
#include <limits>
#include "CudaMatrix.h"
#include "Function1D.h"
#include "Function3D.h"
#include "MatlabReader.h"
@@ -14,6 +15,9 @@
#include "transmissionReconstruction/reconstruction/buildMatrix/DGradient.h"
#include "transmissionReconstruction/reconstruction/buildMatrix/FMM.h"
#include "transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.h"
#include "transmissionReconstruction/reconstruction/buildMatrix/buildMatrix.cuh"
#include "transmissionReconstruction/reconstruction/reconstruction.h"
#include "transmissionReconstruction/reconstruction/solvingEquationSystem/TVAL/TVAL.h"
@@ -23,6 +27,9 @@
#include "reflectionReconstruction/preprocessData/determineOptimalPulse.h"
#include <cusparse.h>
#include <mkl_spblas.h>
@@ -52,6 +59,14 @@ protected:
}
};
TEST_F(Reconstruction_Test, GPUMajor) {
int cudaDevice;
cudaGetDevice(&cudaDevice);
cudaDeviceProp prop;
cudaGetDeviceProperties(&prop, cudaDevice);
EXPECT_GT(prop.major,2)<<"GPU major:"<<prop.major;
}
TEST_F(Reconstruction_Test, determineOptimalPulse) {
Recon::reflectParams::imageResolution = 8.925316499483377e-04;
Recon::reflectParams::optPulseFactor = 48;
@@ -171,30 +186,85 @@ TEST_F(Reconstruction_Test, checkAndScale) {
}
TEST_F(Reconstruction_Test, buildMatrix) {
MatlabReader m("/home/sun/testData/buildMatrix.mat");
auto j1 = m.read("j1");
auto i1 = m.read("i1");
auto s1 = m.read("s1");
auto dims = m.read("dims");
Recon::transParams::bresenham = 1;
MatlabReader m("/home/krad/TestData/buildmatrixFloat.mat");
auto j1 = m.read("jCopy");
auto i1 = m.read("iCopy");
auto s1 = m.read("sCopy");
auto dims = Aurora::Matrix::fromRawData(new float[3]{128, 128, 73}, 1, 3);;
auto receiverList = m.read("receiverList");
auto potentialMap = m.read("potentialMap");
auto senderList = m.read("senderList");
auto res = m.read("res");
auto result = Recon::buildMatrix(senderList, receiverList, res, dims, false, potentialMap);
size_t count = 0;
for(size_t i=0; i<i1.getDataSize(); ++i)
{
EXPECT_FLOAT_AE(i1[i] - 1, result.M.getColVector()[i]);
// EXPECT_FLOAT_AE(i1[i] - 1, result.M.getRowVector()[i])<<"ROW error,"<<" index:"<<i;
if ((int)(i1[i] - 1)!= (int)result.M.getRowVector()[i]){
count++;
}
}
EXPECT_EQ(count, 0)<<"ROW error percent,"<<((double)count/(double)i1.getDataSize())*100.0;
count=0;
for(size_t i=0; i<j1.getDataSize(); ++i)
{
EXPECT_FLOAT_AE(j1[i]- 1, result.M.getRowVector()[i]);
EXPECT_FLOAT_AE(j1[i]- 1, result.M.getColVector()[i])<<"COL error,"<<" index:"<<i;
if ((int)(j1[i] - 1)!= (int)result.M.getColVector()[i]){
count++;
}
}
EXPECT_EQ(count, 0)<<"COL error percent,"<<((double)count/(double)j1.getDataSize())*100.0;
count=0;
for(size_t i=0; i<s1.getDataSize(); ++i)
{
EXPECT_FLOAT_AE(s1[i], result.M.getValVector()[i]);
EXPECT_FLOAT_AE(s1[i], result.M.getValVector()[i])<<"Value error,"<<" index:"<<i;
if ((int)(s1[i] - 1)!= (int)result.M.getValVector()[i]){
count++;
}
}
EXPECT_EQ(count, 0)<<"Value error percent,"<<((double)count/(double)j1.getDataSize())*100.0;
}
TEST_F(Reconstruction_Test, buildMatrixCuda) {
Recon::transParams::bresenham = 1;
MatlabReader m("/home/krad/TestData/buildmatrixFloat.mat");
auto j1 = m.read("jCopy");
auto i1 = m.read("iCopy");
auto s1 = m.read("sCopy");
auto dims = Aurora::Matrix::fromRawData(new float[3]{128, 128, 73}, 1, 3).toDeviceMatrix();
auto receiverList = m.read("receiverList").toDeviceMatrix();
auto potentialMap = m.read("potentialMap").toDeviceMatrix();
auto senderList = m.read("senderList").toDeviceMatrix();
auto res = m.read("res").toDeviceMatrix();
auto result = Recon::buildMatrix(senderList, receiverList, res, dims, false, potentialMap);
size_t count = 0;
for(size_t i=0; i<i1.getDataSize(); ++i)
{
EXPECT_FLOAT_AE(i1[i] - 1, result.M.getRowVector()[i])<<"ROW error,"<<" index:"<<i;
if ((int)(i1[i] - 1)!= (int)result.M.getRowVector()[i]){
count++;
}
}
EXPECT_EQ(count, 0)<<"ROW error percent,"<<((double)count/(double)i1.getDataSize())*100.0;
count=0;
for(size_t i=0; i<j1.getDataSize(); ++i)
{
EXPECT_FLOAT_AE(j1[i]- 1, result.M.getColVector()[i])<<"COL error,"<<" index:"<<i;
if ((int)(j1[i] - 1)!= (int)result.M.getColVector()[i]){
count++;
}
}
EXPECT_EQ(count, 0)<<"COL error percent,"<<((double)count/(double)j1.getDataSize())*100.0;
for(size_t i=0; i<s1.getDataSize(); ++i)
{
EXPECT_FLOAT_EQ(s1[i], result.M.getValVector()[i])<<"Value error,"<<" index:"<<i;
}
EXPECT_EQ(count, 0)<<"Value error percent,"<<((double)count/(double)j1.getDataSize())*100.0;
}
@@ -371,17 +441,64 @@ TEST_F(Reconstruction_Test,callTval3){
opt.mu = 100;
opt.beta = 1;
opt.beta0 = 1;
MatlabReader m("/home/sun/testData/buildMatrixTest.mat");
auto j1 = m.read("j1");
auto i1 = m.read("i1");
auto s1 = m.read("s1");
Aurora::Sparse M(i1-1,j1-1,s1,734989,1196032);
MatlabReader m2("/home/sun/testData/tval3gpu3d.mat");
auto b = m2.read("b");
auto dims = Aurora::Matrix::fromRawData(new float[3]{128,128,73}, 1, 3);
auto result = Recon::callTval3(M, b, dims, 0,opt);
auto outSOS = Recon::slownessToSOS(result, 1.498206569328594e+03) ;
MatlabWriter w2("/home/krad/transmissionSOS111.mat");
w2.setMatrix(outSOS, "SOS");
w2.write();
MatlabReader m("/home/krad/TestData/ijs.mat");
auto j = m.read("j");
auto i = m.read("i");
auto s = m.read("s");
Aurora::Sparse M(i-1,j-1,s,1900092,1196032);
MatlabReader m2("/home/sun/TVALGPU.mat");
auto j1 = m2.read("colPtr");
auto i1 = m2.read("rowPtr");
auto s1 = m2.read("data");
int * yIdxs = new int[M.getColVector().getDataSize()];
std::copy(M.getColVector().getData(),M.getColVector().getData()+M.getColVector().getDataSize(),yIdxs);
int * xIdxs = new int[M.getRowVector().getDataSize()];
std::copy(M.getRowVector().getData(),M.getRowVector().getData()+M.getRowVector().getDataSize(),xIdxs);
Aurora::Matrix values = M.getValVector();
size_t rows = M.getM(), cols = M.getN();
int nz = std::max(M.getColVector().getDataSize(),std::max(M.getRowVector().getDataSize(),M.getValVector().getDataSize()));
sparse_matrix_t A;
sparse_matrix_t csrA;
sparse_matrix_t csrA_t;
sparse_status_t status;
status = mkl_sparse_s_create_coo(&A, sparse_index_base_t::SPARSE_INDEX_BASE_ZERO, rows,cols, nz, xIdxs,yIdxs,values.getData());
status = mkl_sparse_convert_csr(A, sparse_operation_t::SPARSE_OPERATION_TRANSPOSE, &csrA);
// status = mkl_sparse_order(csrA);
int n_rows = 0;
int n_cols = 0;
int *col_start,*col_end,*row_idx;
float * csrValues;
sparse_index_base_t index;
status = mkl_sparse_s_export_csr(csrA, &index, &n_cols, &n_rows, &col_start, &col_end, &row_idx, &csrValues);
int *col_idx = new int[n_cols+1];
std::copy(col_start,col_start+ n_cols,col_idx);
col_idx[n_cols] = nz;
printf("j1 count:%zu,i1 count:%zu,s1 count:%zu\r\n",j1.getDataSize(),i1.getDataSize(),s1.getDataSize());
printf("n_rows:%d,n_cols:%d\r\n", n_rows, n_cols);
for (size_t k = 0; k < j1.getDataSize(); k++)
{
ASSERT_FLOAT_EQ(j1.getData()[k],col_idx[k])<<"col error, index:"<<k;
}
for (size_t k = 0; k < i1.getDataSize(); k++)
{
ASSERT_FLOAT_EQ(i1.getData()[k],row_idx[k])<<"row error, index:"<<k;
}
for (size_t k = 0; k < s1.getDataSize(); k++)
{
ASSERT_FLOAT_EQ(s1.getData()[k],csrValues[k])<<"value error, index:"<<k;
}
mkl_sparse_destroy(A);
mkl_sparse_destroy(csrA);
delete [] xIdxs;
delete [] yIdxs;
}