Add blockingGeometryInfos.

This commit is contained in:
sunwen
2023-05-24 14:25:23 +08:00
parent db86d973a6
commit 5bd7a15238
5 changed files with 145 additions and 8 deletions

View File

@@ -0,0 +1,74 @@
#include "blockingGeometryInfo.h"
#include "Function.h"
#include "Function2D.h"
#include "Function3D.h"
using namespace Aurora;
using namespace Recon;
GeometryBlock Recon::blockingGeometryInfos(const GeometryInfo& aGeom, const Aurora::Matrix& aRnBlock, const Aurora::Matrix& aRlBlock,
const Aurora::Matrix& aSnBlock, const Aurora::Matrix& aSlBlock, const Aurora::Matrix& aMpBlock)
{
GeometryBlock result;
int numData = aRnBlock.getDimSize(1);
Matrix arrayOfDataSize = zeros(1, numData, 1) + 1;
int dim = aGeom.receiverPositions[0].getDimSize(0);
double* receiverPositionBlockData = Aurora::malloc(dim * numData);
double* receiverNormalBlockData = Aurora::malloc(dim * numData);
result.receiverPositionBlock = Matrix::New(receiverPositionBlockData, dim, numData);
result.receiverNormalBlock = Matrix::New(receiverNormalBlockData, dim, numData);
double* sizeData = Aurora::malloc(4);
sizeData[0] = aGeom.receiverPositions[0].getDimSize(0);
sizeData[1] = aGeom.receiverPositions[0].getDimSize(1);
sizeData[2] = aGeom.receiverPositions[0].getDimSize(2);
sizeData[3] = aGeom.receiverPositions.size();
size_t matrixSize = aGeom.receiverPositions[0].getDataSize();
Matrix receiverPositionsSize = Matrix::New(sizeData, 1, 4);
for(int i=0; i<dim; ++i)
{
Matrix idx = sub2ind(receiverPositionsSize, {arrayOfDataSize * (i+1), aRnBlock, aRlBlock, aMpBlock});
for(int j=0; j<idx.getDataSize(); ++j)
{
unsigned long long index4D = (idx[j] - 1) / matrixSize;
unsigned long long index3D = idx[j] - 1;
if(index4D > 0)
{
index3D -= matrixSize;
}
receiverPositionBlockData[dim*j + i] = aGeom.receiverPositions[index4D][index3D];
receiverNormalBlockData[dim*j + i] = aGeom.receiverNormals[index4D][index3D];
}
}
double* senderPositionBlockData = Aurora::malloc(dim * numData);
double* senderNormalBlockData = Aurora::malloc(dim * numData);
result.senderPositionBlock = Matrix::New(senderPositionBlockData, dim, numData);
result.senderNormalBlock = Matrix::New(senderNormalBlockData, dim, numData);
sizeData = Aurora::malloc(4);
sizeData[0] = aGeom.senderPositions[0].getDimSize(0);
sizeData[1] = aGeom.senderPositions[0].getDimSize(1);
sizeData[2] = aGeom.senderPositions[0].getDimSize(2);
sizeData[3] = aGeom.senderPositions.size();
matrixSize = aGeom.receiverPositions[0].getDataSize();
Matrix senderPositionsSize = Matrix::New(sizeData, 1, 4);
dim = aGeom.senderPositions[0].getDimSize(0);
for(int i=0; i<dim; ++i)
{
Matrix idx = sub2ind(senderPositionsSize, {arrayOfDataSize * (i+1), aSnBlock, aSlBlock, aMpBlock});
for(int j=0; j<idx.getDataSize(); ++j)
{
unsigned long long index4D = (idx[j] - 1) / matrixSize;
unsigned long long index3D = idx[j] - 1;
if(index4D > 0)
{
index3D -= matrixSize;
}
senderPositionBlockData[dim*j + i] = aGeom.senderPositions[index4D][index3D];
senderNormalBlockData[dim*j + i] = aGeom.senderNormals[index4D][index3D];
}
}
return result;
}

View File

@@ -0,0 +1,22 @@
#ifndef BLOCKING_GEOMETRYINFO_H
#define BLOCKING_GEOMETRYINFO_H
#include "Matrix.h"
#include "src/common/getGeometryInfo.h"
namespace Recon
{
struct GeometryBlock
{
Aurora::Matrix senderPositionBlock;
Aurora::Matrix senderNormalBlock;
Aurora::Matrix receiverPositionBlock;
Aurora::Matrix receiverNormalBlock;
};
GeometryBlock blockingGeometryInfos(const GeometryInfo& aGeom, const Aurora::Matrix& aRnBlock, const Aurora::Matrix& aRlBlock,
const Aurora::Matrix& aSnBlock, const Aurora::Matrix& aSlBlock, const Aurora::Matrix& aMpBlock);
}
#endif

View File

@@ -0,0 +1,102 @@
#include "getAscanBlock.h"
#include "Function.h"
#include "Matrix.h"
#include "Parser.h"
#include "Data/OneTasAScanData.h"
#include "Data/AScanData.h"
#include "Data/TasElementIndex.h"
#include "Data/ElementIndex.h"
#include "Data/GeometryIndex.h"
#include "Data/ElectricIndex.h"
#include "Data/MetaData.h"
#include "ShotList/ShotList.h"
#include <vector>
using namespace Recon;
using namespace Aurora;
namespace
{
std::vector<int> findTasIndex(TasIndicesPointer aTasIndices, int aTasNum)
{
std::vector<int> result;
for(int i=0; i<aTasIndices.getLength(); ++i)
{
if (aTasIndices.get()[i] == aTasNum)
{
result.push_back(i);
}
}
return result;
}
}
AscanBlock Recon::getAscanBlock(Parser* aParser, const Aurora::Matrix& aMp, const Aurora::Matrix& aSl, const Aurora::Matrix& aSn,
const Aurora::Matrix& aRl, const Aurora::Matrix& aRn)
{
AscanBlock result;
int numScans = aMp.getDataSize() * aSl.getDataSize() * aSn.getDataSize() * aRl.getDataSize() * aRn.getDataSize();
int numScansIndex = 0;
int ascanBlockSize = numScans * aParser->getMetaData().getSampleNumber();
double* ascanBlockData = new double[ascanBlockSize];
double* mpBlockData = new double[numScans];
double* slBlockData = new double[numScans];
double* snBlockData = new double[numScans];
double* rlBlockData = new double[numScans];
double* rnBlockData = new double[numScans];
double* gainBlockData = new double[numScans];
result.ascanBlock = Matrix::fromRawData(ascanBlockData, aParser->getMetaData().getSampleNumber(), numScans);
result.mpBlock = Matrix::fromRawData(mpBlockData, 1, numScans);
result.slBlock = Matrix::fromRawData(slBlockData, 1, numScans);
result.snBlock = Matrix::fromRawData(snBlockData, 1, numScans);
result.rlBlock = Matrix::fromRawData(rlBlockData, 1, numScans);
result.rnBlock = Matrix::fromRawData(rnBlockData, 1, numScans);
result.gainBlock = Matrix::fromRawData(gainBlockData, 1, numScans);
OneTasAScanData oenTasData = aParser->getOneTasAscanDataOfMotorPosition(1);
auto tasIndices = aParser->getMetaData().getTasIndices();
auto receiverIndices = aParser->getMetaData().getReceiverIndices();
for(int mpIndex=0; mpIndex<aMp.getDataSize();++mpIndex)
{
for(int slIndex=0; slIndex<aSl.getDataSize();++slIndex)
{
OneTasAScanData oenTasData = aParser->getOneTasAscanDataOfMotorPosition(aSl[slIndex]);
for(int snIndex=0; snIndex<aSn.getDataSize();++snIndex)
{
for(int rlIndex=0; rlIndex<aRl.getDataSize();++rlIndex)
{
auto rElementIndex = findTasIndex(tasIndices, aRl[rlIndex]);
for(int rnIndex=0; rnIndex<rElementIndex.size(); ++rnIndex)
{
MotorPosition mp = aMp[mpIndex] == 1 ? MotorPosition::MotorPosition1 : MotorPosition::MotorPosition2;
AScanData ascan = aParser->searchAscanDataFromOneTasAscanDataOfMP(oenTasData, ElementIndex(GeometryIndex(aSn[snIndex])), TasElementIndex(aRl[rlIndex], ElementIndex(GeometryIndex(receiverIndices.get()[rElementIndex[rnIndex]]))), mp);
std::copy(ascan.get() ,ascan.get() + ascan.getAscanDataLength(), ascanBlockData);
ascanBlockData += ascan.getAscanDataLength();
mpBlockData[numScansIndex] = aMp[mpIndex];
slBlockData[numScansIndex] = aSl[slIndex];
snBlockData[numScansIndex] = aSn[snIndex];
rlBlockData[numScansIndex] = aRl[rlIndex];
rnBlockData[numScansIndex] = receiverIndices.get()[rElementIndex[rnIndex]];
gainBlockData[numScansIndex] = ascan.getAmplification()[0];
++numScansIndex;
}
}
}
}
}
// 改分支不会进,暂不实现 todo
// % initialized and calculated numbers agree?, otherwise reduce block sizes
// if num ~= numScans
// warning([mfilename, ':MissingSenderReceiverInfos'],'Not all defined data available.')
// writeReconstructionLog('Not all defined sender/receiver available', 3);
// usedDataIdxs = 1:num;
// [slBlock, snBlock, rlBlock, rnBlock, mpBlock, gainBlock] = removeDataFromArrays(usedDataIdxs, slBlock, snBlock, rlBlock, rnBlock, mpBlock, gainBlock);
// AscanBlock = AscanBlock(:,usedDataIdxs);
// end
return result;
}

View File

@@ -0,0 +1,25 @@
#ifndef GET_ASCAN_BLOCK_H
#define GET_ASCAN_BLOCK_H
#include "Matrix.h"
#include "src/common/getMeasurementMetaData.h"
class Parser;
namespace Recon
{
struct AscanBlock
{
Aurora::Matrix ascanBlock;
Aurora::Matrix mpBlock;
Aurora::Matrix slBlock;
Aurora::Matrix snBlock;
Aurora::Matrix rlBlock;
Aurora::Matrix rnBlock;
Aurora::Matrix gainBlock;
};
AscanBlock getAscanBlock(Parser* aParser, const Aurora::Matrix& mp, const Aurora::Matrix& aSL, const Aurora::Matrix& aSn,
const Aurora::Matrix& aRl, const Aurora::Matrix& aRn);
}
#endif