Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Ref ARRUS-113: arbitrary rx apodization #248

Merged
merged 6 commits into from
Dec 13, 2021
Merged

Ref ARRUS-113: arbitrary rx apodization #248

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
231f825
iqRaw2Lri - custom apodization (uses texture)
PiotrKarwat Nov 22, 2021
eff0bcf
iqRaw2Lri - change rxApod position among other parameters
PiotrKarwat Nov 22, 2021
a52383c
Us4R & Reconstruction - include rxApod
PiotrKarwat Nov 22, 2021
cbd1bf5
Control scripts - include rxApod
PiotrKarwat Nov 22, 2021
824c1e8
iqRaw2Lri - checkData - validation of int32 arrays added
PiotrKarwat Nov 23, 2021
05e4351
Reconstruction - default rxApod (uniform)
PiotrKarwat Nov 23, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions api/matlab/arrus/Reconstruction.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@
decimation
xGrid
zGrid
rxApod = [1 1]
bmodeEnable = true
colorEnable = false
vectorEnable = false
Expand Down
4 changes: 4 additions & 0 deletions api/matlab/arrus/Us4R.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -203,6 +203,7 @@ function upload(obj, sequenceOperation, reconstructOperation)
'decimation', reconstructOperation.decimation, ...
'xGrid', reconstructOperation.xGrid, ...
'zGrid', reconstructOperation.zGrid, ...
'rxApod', reconstructOperation.rxApod, ...
'bmodeEnable', reconstructOperation.bmodeEnable, ...
'colorEnable', reconstructOperation.colorEnable, ...
'vectorEnable', reconstructOperation.vectorEnable, ...
Expand Down Expand Up @@ -548,6 +549,7 @@ function setRecParams(obj,varargin)
'decimation', 'dec'; ...
'xGrid', 'xGrid'; ...
'zGrid', 'zGrid'; ...
'rxApod', 'rxApod'; ...
'bmodeEnable', 'bmodeEnable'; ...
'colorEnable', 'colorEnable'; ...
'vectorEnable', 'vectorEnable'; ...
Expand Down Expand Up @@ -599,6 +601,7 @@ function setRecParams(obj,varargin)
obj.sys.tangElem = gpuArray(single(obj.sys.tangElem));
obj.rec.zGrid = gpuArray(single(obj.rec.zGrid));
obj.rec.xGrid = gpuArray(single(obj.rec.xGrid));
obj.rec.rxApod = gpuArray(single(obj.rec.rxApod));
obj.seq.txFoc = gpuArray(single(obj.seq.txFoc));
obj.seq.txAngZX = gpuArray(single(obj.seq.txAngZX));
obj.seq.txApCentZ = gpuArray(single(obj.seq.txApCentZ));
Expand Down Expand Up @@ -1165,6 +1168,7 @@ function closeSequence(obj)
obj.sys.tangElem, ...
obj.rec.zGrid, ...
obj.rec.xGrid, ...
obj.rec.rxApod, ...
obj.seq.txFoc(selFrames), ...
obj.seq.txAngZX(selFrames), ...
obj.seq.txApCentZ(selFrames), ...
Expand Down
122 changes: 75 additions & 47 deletions api/matlab/arrus/mexcuda/iqRaw2Lri.cu
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ __constant__ float xElemConst[256];
__constant__ float tangElemConst[256];

texture <float2, cudaTextureType1DLayered, cudaReadModeElementType> iqRawTex;

texture <float, cudaTextureType1D, cudaReadModeElementType> rxApodTex;

__forceinline__ __device__ float ownHypotf(float x, float y)
{
Expand Down Expand Up @@ -43,10 +43,7 @@ __global__ void iqRaw2Lri( float2 * iqLri, float const * zPix, float const * xP
float const omega = 2 * M_PI * fn;
float const sosInv = 1 / sos;
// float const zDistInv = 1 / zPix[z];
float const nSigma = 3; // number of sigmas in half of the apodization Gaussian curve
float const twoSigSqrInv = nSigma * nSigma * 0.5f;
float const rngRxTangInv = 2 / (maxRxTang - minRxTang); // inverted half range
float const centRxTang = (maxRxTang + minRxTang) * 0.5f;
float const rngRxTangInv = 1 / (maxRxTang - minRxTang); // inverted tangent range

for (int iTx=0; iTx<nTx; iTx++) {

Expand Down Expand Up @@ -109,8 +106,8 @@ __global__ void iqRaw2Lri( float2 * iqLri, float const * zPix, float const * xP
rxTang = __fdividef(xPix[x] - xElemConst[iElem], zPix[z] - zElemConst[iElem]);
rxTang = __fdividef(rxTang-tangElemConst[iElem], 1.f+rxTang*tangElemConst[iElem]);
if (rxTang < minRxTang || rxTang > maxRxTang) continue;
rxApod = (rxTang-centRxTang)*rngRxTangInv;
rxApod = __expf(-rxApod*rxApod*twoSigSqrInv);
rxApod = (rxTang-minRxTang)*rngRxTangInv; // <0,1>, needs normalized texture fetching, errors at aperture sided
rxApod = tex1D(rxApodTex, rxApod);

time = (txDist + rxDist) * sosInv + initDel;
iSamp = time * fs;
Expand All @@ -133,27 +130,36 @@ __global__ void iqRaw2Lri( float2 * iqLri, float const * zPix, float const * xP
}
}

__host__ void checkData(mxGPUArray const * const data, char const * const name, bool const isComplex, int const nDims, char const * const invalidInputMsgId)
__host__ void checkData(mxGPUArray const * const data,
char const * const name,
bool const mustBeInt,
bool const mustBeComplex,
int const mustBeNDim,
char const * const invalidInputMsgId)
{
std::string invalidInputMsgTxt(name);

if (mxGPUGetClassID(data) != mxSINGLE_CLASS)
if (mustBeInt && mxGPUGetClassID(data) != mxINT32_CLASS)
invalidInputMsgTxt += std::string(" must be int32.");

else if (!mustBeInt && mxGPUGetClassID(data) != mxSINGLE_CLASS)
invalidInputMsgTxt += std::string(" must be single.");

else if (!isComplex && mxGPUGetComplexity(data))
else if (!mustBeComplex && mxGPUGetComplexity(data))
invalidInputMsgTxt += std::string(" must be real.");

else if (isComplex && !mxGPUGetComplexity(data))
else if (mustBeComplex && !mxGPUGetComplexity(data))
invalidInputMsgTxt += std::string(" must be complex.");

else if (nDims==1 && !( mxGPUGetNumberOfDimensions(data) == 1 ||
(mxGPUGetNumberOfDimensions(data) == 2 && mxGPUGetDimensions(data)[0] == 1)))
else if (mustBeNDim==1 && !( mxGPUGetNumberOfDimensions(data) == 1 ||
(mxGPUGetNumberOfDimensions(data) == 2 &&
mxGPUGetDimensions(data)[0] == 1)))
invalidInputMsgTxt += std::string(" must be at most 1D vector.");

else if (nDims==2 && !(mxGPUGetNumberOfDimensions(data) <= 2))
else if (mustBeNDim==2 && !(mxGPUGetNumberOfDimensions(data) <= 2))
invalidInputMsgTxt += std::string(" must be at most 2D array.");

else if (nDims==3 && !(mxGPUGetNumberOfDimensions(data) <= 3))
else if (mustBeNDim==3 && !(mxGPUGetNumberOfDimensions(data) <= 3))
invalidInputMsgTxt += std::string(" must be at most 3D array.");

else
Expand All @@ -177,6 +183,7 @@ void mexFunction(int nlhs, mxArray * plhs[],
mxGPUArray const * tangElem;
mxGPUArray const * zPix;
mxGPUArray const * xPix;
mxGPUArray const * rxApod;
mxGPUArray const * foc;
mxGPUArray const * ang;
mxGPUArray const * centZ;
Expand All @@ -193,6 +200,7 @@ void mexFunction(int nlhs, mxArray * plhs[],
float const * dev_tangElem;
float const * dev_zPix;
float const * dev_xPix;
float const * dev_rxApod;
float const * dev_foc;
float const * dev_ang;
float const * dev_centZ;
Expand All @@ -215,6 +223,7 @@ void mexFunction(int nlhs, mxArray * plhs[],
int nXPix;
int nRx;
int nTx;
int nRxApodSamp;

dim3 const threadsPerBlock = {16, 16, 1};
dim3 blocksPerGrid;
Expand All @@ -224,15 +233,15 @@ void mexFunction(int nlhs, mxArray * plhs[],
char const * const invalidOutputMsgId = "iqRaw2Lri:InvalidOutput";

/* Validate mex inputs/outputs */
if (nrhs!=20) {
mexErrMsgIdAndTxt( invalidInputMsgId, "20 inputs required");
if (nrhs!=21) {
mexErrMsgIdAndTxt( invalidInputMsgId, "21 inputs required");
}

if (nlhs>1) {
mexErrMsgIdAndTxt( invalidOutputMsgId, "One output allowed");
}

// for (int i=13; i<19; i++) {
// for (int i=15; i<20; i++) {
// if (!mxIsSingle(prhs[i]) || mxIsComplex(prhs[i]) || mxGetNumberOfElements(prhs[i]) != 1) {
// mexErrMsgIdAndTxt( invalidInputMsgId, "Last 6 inputs must be single, real scalars");
// }
Expand All @@ -246,44 +255,47 @@ void mexFunction(int nlhs, mxArray * plhs[],
tangElem = mxGPUCreateFromMxArray(prhs[3]);
zPix = mxGPUCreateFromMxArray(prhs[4]);
xPix = mxGPUCreateFromMxArray(prhs[5]);
foc = mxGPUCreateFromMxArray(prhs[6]);
ang = mxGPUCreateFromMxArray(prhs[7]);
centZ = mxGPUCreateFromMxArray(prhs[8]);
centX = mxGPUCreateFromMxArray(prhs[9]);
elemFst = mxGPUCreateFromMxArray(prhs[10]);
elemLst = mxGPUCreateFromMxArray(prhs[11]);
rxElemOrig = mxGPUCreateFromMxArray(prhs[12]);
nSampOmit = mxGPUCreateFromMxArray(prhs[13]);

minRxTang = mxGetScalar(prhs[14]);
maxRxTang = mxGetScalar(prhs[15]);
fs = mxGetScalar(prhs[16]);
fn = mxGetScalar(prhs[17]);
sos = mxGetScalar(prhs[18]);
initDel = mxGetScalar(prhs[19]);
rxApod = mxGPUCreateFromMxArray(prhs[6]);
foc = mxGPUCreateFromMxArray(prhs[7]);
ang = mxGPUCreateFromMxArray(prhs[8]);
centZ = mxGPUCreateFromMxArray(prhs[9]);
centX = mxGPUCreateFromMxArray(prhs[10]);
elemFst = mxGPUCreateFromMxArray(prhs[11]);
elemLst = mxGPUCreateFromMxArray(prhs[12]);
rxElemOrig= mxGPUCreateFromMxArray(prhs[13]);
nSampOmit = mxGPUCreateFromMxArray(prhs[14]);

minRxTang = mxGetScalar(prhs[15]);
maxRxTang = mxGetScalar(prhs[16]);
fs = mxGetScalar(prhs[17]);
fn = mxGetScalar(prhs[18]);
sos = mxGetScalar(prhs[19]);
initDel = mxGetScalar(prhs[20]);

/* Validate inputs */
checkData(iqRaw, "iqRaw", true, 3, invalidInputMsgId);
checkData(zElem, "zElem", false, 1, invalidInputMsgId);
checkData(xElem, "xElem", false, 1, invalidInputMsgId);
checkData(tangElem, "tangElem", false, 1, invalidInputMsgId);
checkData(zPix, "zPix", false, 1, invalidInputMsgId);
checkData(xPix, "xPix", false, 1, invalidInputMsgId);
checkData(foc, "foc", false, 1, invalidInputMsgId);
checkData(ang, "ang", false, 1, invalidInputMsgId);
checkData(centZ, "centZ", false, 1, invalidInputMsgId);
checkData(centX, "centX", false, 1, invalidInputMsgId);
// checkData(elemFst, "elemFst", false, 1, invalidInputMsgId); //int
// checkData(elemLst, "elemLst", false, 1, invalidInputMsgId); //int
// checkData(rxElemOrig,"rxElemOrig",false, 1, invalidInputMsgId); //int
// checkData(nSampOmit,"nSampOmit",false, 1, invalidInputMsgId); //int
checkData(iqRaw, "iqRaw", false, true, 3, invalidInputMsgId);
checkData(zElem, "zElem", false, false, 1, invalidInputMsgId);
checkData(xElem, "xElem", false, false, 1, invalidInputMsgId);
checkData(tangElem, "tangElem", false, false, 1, invalidInputMsgId);
checkData(zPix, "zPix", false, false, 1, invalidInputMsgId);
checkData(xPix, "xPix", false, false, 1, invalidInputMsgId);
checkData(rxApod, "rxApod", false, false, 1, invalidInputMsgId);
checkData(foc, "foc", false, false, 1, invalidInputMsgId);
checkData(ang, "ang", false, false, 1, invalidInputMsgId);
checkData(centZ, "centZ", false, false, 1, invalidInputMsgId);
checkData(centX, "centX", false, false, 1, invalidInputMsgId);
checkData(elemFst, "elemFst", true, false, 1, invalidInputMsgId);
checkData(elemLst, "elemLst", true, false, 1, invalidInputMsgId);
checkData(rxElemOrig,"rxElemOrig",true, false, 1, invalidInputMsgId);
checkData(nSampOmit, "nSampOmit", true, false, 1, invalidInputMsgId);

/* Get some additional information */
nSamp = mxGPUGetDimensions(iqRaw)[0];
nRx = mxGPUGetDimensions(iqRaw)[1];
nElem = mxGPUGetNumberOfElements(xElem);
nZPix = mxGPUGetNumberOfElements(zPix);
nXPix = mxGPUGetNumberOfElements(xPix);
nRxApodSamp = mxGPUGetNumberOfElements(rxApod);
if (mxGPUGetNumberOfDimensions(iqRaw)<3) {
nTx = 1;
}
Expand Down Expand Up @@ -313,6 +325,7 @@ void mexFunction(int nlhs, mxArray * plhs[],
dev_tangElem = static_cast<float const *>(mxGPUGetDataReadOnly(tangElem));
dev_zPix = static_cast<float const *>(mxGPUGetDataReadOnly(zPix));
dev_xPix = static_cast<float const *>(mxGPUGetDataReadOnly(xPix));
dev_rxApod = static_cast<float const *>(mxGPUGetDataReadOnly(rxApod));
dev_foc = static_cast<float const *>(mxGPUGetDataReadOnly(foc));
dev_ang = static_cast<float const *>(mxGPUGetDataReadOnly(ang));
dev_centZ = static_cast<float const *>(mxGPUGetDataReadOnly(centZ));
Expand All @@ -330,6 +343,17 @@ void mexFunction(int nlhs, mxArray * plhs[],
cudaMemcpyToSymbol( xElemConst, dev_xElem, nElem*sizeof(float), 0, cudaMemcpyDeviceToDevice);
cudaMemcpyToSymbol(tangElemConst, dev_tangElem, nElem*sizeof(float), 0, cudaMemcpyDeviceToDevice);

/* configure texture reference (apodization) */
rxApodTex.normalized = true;
rxApodTex.addressMode[0] = cudaAddressModeBorder;
rxApodTex.filterMode = cudaFilterModeLinear;

cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc(32, 0, 0, 0, cudaChannelFormatKindFloat);
cudaArray* cuArrayApod;
cudaMallocArray(&cuArrayApod, &channelDesc, nRxApodSamp, 0);
cudaMemcpyToArray(cuArrayApod, 0, 0, dev_rxApod, nRxApodSamp*sizeof(float), cudaMemcpyDeviceToDevice);
cudaBindTextureToArray(rxApodTex, cuArrayApod, channelDesc);

/* configure texture reference */
iqRawTex.normalized = false;
iqRawTex.addressMode[0] = cudaAddressModeBorder;
Expand Down Expand Up @@ -386,13 +410,17 @@ void mexFunction(int nlhs, mxArray * plhs[],
cudaUnbindTexture(iqRawTex);
cudaFreeArray(cuArray);

cudaUnbindTexture(rxApodTex);
cudaFreeArray(cuArrayApod);

mxGPUDestroyGPUArray(iqLri);
mxGPUDestroyGPUArray(iqRaw);
mxGPUDestroyGPUArray(zElem);
mxGPUDestroyGPUArray(xElem);
mxGPUDestroyGPUArray(tangElem);
mxGPUDestroyGPUArray(zPix);
mxGPUDestroyGPUArray(xPix);
mxGPUDestroyGPUArray(rxApod);
mxGPUDestroyGPUArray(foc);
mxGPUDestroyGPUArray(ang);
mxGPUDestroyGPUArray(centZ);
Expand Down
3 changes: 2 additions & 1 deletion api/matlab/examples/Us4R_ATL_control.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,8 @@
'cicOrder', 2, ...
'decimation', 4, ...
'xGrid', (-20:0.10:20)*1e-3, ...
'zGrid', ( 0:0.10:100)*1e-3);
'zGrid', ( 0:0.10:100)*1e-3, ...
'rxApod', hamming(20).');

us.upload(seqPWI,rec);
% us.upload(seqSTA,rec);
Expand Down
3 changes: 2 additions & 1 deletion api/matlab/examples/Us4R_ATL_control_batches.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -92,7 +92,8 @@
'cicOrder', 1, ...
'decimation', 1, ...
'xGrid', xGrid, ...
'zGrid', zGrid);
'zGrid', zGrid, ...
'rxApod', hamming(20).');
end

%% Prepare the acquisition (& reconstruction)
Expand Down
3 changes: 2 additions & 1 deletion api/matlab/examples/Us4R_Olympus_control.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,8 @@
'cicOrder', 2, ...
'decimation', 4, ...
'xGrid', (-40:0.20:40)*1e-3, ...
'zGrid', ( 0:0.20:110)*1e-3);
'zGrid', ( 0:0.20:110)*1e-3, ...
'rxApod', hamming(20).');

% us.upload(seqSTA, rec);
us.upload(seqPWI, rec);
Expand Down
3 changes: 2 additions & 1 deletion api/matlab/examples/Us4R_Ultrasonix_control.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,7 +59,8 @@
'cicOrder', 2, ...
'decimation', 4, ...
'xGrid', (-20:0.10:20)*1e-3, ...
'zGrid', ( 0:0.10:50)*1e-3);
'zGrid', ( 0:0.10:50)*1e-3, ...
'rxApod', hamming(20).');

us.upload(seqPWI,rec);
% us.upload(seqSTA,rec);
Expand Down
3 changes: 2 additions & 1 deletion api/matlab/examples/Us4R_Vermon_control.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -71,7 +71,8 @@
'cicOrder', 2, ...
'decimation', 1, ...
'xGrid', (-5:0.05:5)*1e-3, ...
'zGrid', ( 0:0.05:10)*1e-3);
'zGrid', ( 0:0.05:10)*1e-3, ...
'rxApod', hamming(20).');

% us.upload(seqSTA, rec);
us.upload(seqPWI, rec);
Expand Down
3 changes: 2 additions & 1 deletion api/matlab/examples/Us4R_control.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -72,7 +72,8 @@
'cicOrder', 2, ...
'decimation', 4, ...
'xGrid', (-20:0.10:20)*1e-3, ...
'zGrid', ( 0:0.10:50)*1e-3);
'zGrid', ( 0:0.10:50)*1e-3, ...
'rxApod', hamming(20).');

% us.upload(seqSTA, rec);
us.upload(seqPWI, rec);
Expand Down
1 change: 1 addition & 0 deletions api/matlab/examples/Us4R_duplex.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -74,6 +74,7 @@
'decimation', 2, ...
'xGrid', (-7:0.025:7)*1e-3, ...
'zGrid', ( 0:0.025:7)*1e-3, ...
'rxApod', hamming(20).', ...
'colorEnable', true, ...
'vectorEnable', false, ...
'bmodeFrames', 1:9, ...
Expand Down
3 changes: 2 additions & 1 deletion api/matlab/examples/Us4R_maxSequence.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -70,7 +70,8 @@
'cicOrder', 2, ...
'decimation', 4, ...
'xGrid', (-20:0.10:20)*1e-3, ...
'zGrid', ( 0:0.10:50)*1e-3);
'zGrid', ( 0:0.10:50)*1e-3, ...
'rxApod', hamming(20).');

us.upload(seqSTA,rec);
% us.upload(seqPWI,rec);
Expand Down