Skip to content

Commit

Permalink
Ref ARRUS-113: arbitrary rx apodization (#248)
Browse files Browse the repository at this point in the history 
* iqRaw2Lri - custom apodization (uses texture)

Replaces previous (fixed gaussian) apodization

* iqRaw2Lri - change rxApod position among other parameters

* Us4R & Reconstruction - include rxApod

* Control scripts - include rxApod

* iqRaw2Lri - checkData - validation of int32 arrays added

* Reconstruction - default rxApod (uniform)
  • Loading branch information
@PiotrKarwat
PiotrKarwat authored Dec 13, 2021
1 parent 8076024 commit 9a314f6
Show file tree
Hide file tree
Showing 11 changed files with 95 additions and 54 deletions.
1 change: 1 addition & 0 deletions api/matlab/arrus/Reconstruction.m
View file
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
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
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
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
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
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
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
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
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
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
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

0 comments on commit 9a314f6

Please sign in to comment.