I am trying to adapt the matched filter algorithm (given here) for arrayfire version 3.6.4. Here is what I ended up with:
#include <arrayfire.h>
using namespace af;
struct SAR_data { //! SAR data structure format
double deltaF; //! step size of frequency data (Hz)
double maxWr; //! maximum scene size in range direction (m)
double maxWx; //! maximum scene size in cross-range direction (m)
array minF; //! (Np x 1): Vector containing the start frequency of
//! each pulse (Hz)
array x_mat; //! (Sx x Sy): the x-positions of each pixel (m)
array y_mat; //! (Sx x Sy): the y-positions of each pixel (m)
array z_mat; //! (Sx x Sy): the z-positions of each pixel (m)
array AntX; //! (Np x 1): the x-position of the sensor at each pulse (m)
array AntY; //! (Np x 1): the y-position of the sensor at each pulse (m)
array AntZ; //! (Np x 1): the z-position of the sensor at each pulse (m)
array R0; //! (Np x 1): the range to scene center (m)
array phdata; //! (K x Np): phase history data (frequency domain),
//! fast time in rows, slow time in columns
array im_final; //! (Sx x Sy): the complex image value at each pixel
};
// matched filter realization using ArrayFire
void matched_filter_ArrayFire(SAR_data& data) {
double c = 299792458.0; // speed of light
// Determine the size of the phase history data
int K = data.phdata.dims(0); // # of frequency bins per pulse
int Np = data.phdata.dims(1); // # of pulses
// Initialize the image with all zero values (complex)
data.im_final = constant<cdouble>(0, data.x_mat.dims(), c64);
array im_slices = constant<cdouble>(0, K, data.x_mat.dims(0),
data.x_mat.dims(1), c64);
array fspan = array(seq(0.0, K-1)) * data.deltaF;
cdouble unit = {0, 1};
for (int ii = 0; ii < Np; ii++) {
// compute differential range for each image pixel (m)
array dx = data.AntX(ii) - data.x_mat;
array dy = data.AntY(ii) - data.y_mat;
array dz = data.AntZ(ii) - data.z_mat;
array dR = sqrt(dx*dx + dy*dy + dz*dz) - data.R0(ii);
// calculate the frequency of each sample in the pulse (Hz)
array freq = data.minF(ii) + fspan;
array tt = data.phdata(span,ii);
// perform the Matched Filter operation
gfor (array jj, K) {
im_slices(jj,span,span) =
tt(jj)*exp(unit*((double)(4.0*Pi/c)*freq(jj)*dR));
}
tt = sum(im_slices,0);
data.im_final = data.im_final + moddims(tt, data.x_mat.dims());
}
}
The code compiles well but when I try to run it, it throws an exception:
What() is:ArrayFire Exception (Invalid input size:203):
In function class af::dim4 __cdecl getOutDims(const class af::dim4 &,const class af::dim4 &,bool)
In file src\backend\common\ArrayInfo.cpp:130
Invalid dimension for argument 1
Expected: ldims == rdims
In function class af::array __cdecl af::operator -(const class af::array &,const class af::array &)
In file src\api\cpp\array.cpp:838
As I found out, it complains about the following lines:
array dx = data.AntX(ii) - data.x_mat;
array dy = data.AntY(ii) - data.y_mat;
array dz = data.AntZ(ii) - data.z_mat;
since the dimensions of data.AntX(ii) (which is 1x1) and data.x_mat (which is Sx x Sy) do not match.. Though, in older versions of ArrayFire it used to work fine. Apparently, I need some way to tell ArrayFire to "expand" the value of data.AntX(ii) for the whole matrix dimension of data.x_mat. But how to do it ?
Working version of the algorithm using the hint of @pradeep:
// matched filter realization using ArrayFire
void matched_filter_ArrayFire(SAR_data& data) {
double c = 299792458.0; // speed of light
// Determine the size of the phase history data
int K = data.phdata.dims(0); // # of frequency bins per pulse
int Np = data.phdata.dims(1); // # of pulses
// Initialize the image with all zero values (complex)
data.im_final = constant<cdouble>({0,0}, data.x_mat.dims(), c64);
array im_slices = constant<cdouble>({0,0}, data.x_mat.dims(0),
data.x_mat.dims(1), K, c64);
array fspan = array(seq(0.0, K-1)) * data.deltaF;
cdouble unit = {0, 1};
for (int ii = 0; ii < Np; ii++) {
// compute differential range for each image pixel (m)
array dx = tile(data.AntX(ii), data.x_mat.dims()) - data.x_mat;
array dy = tile(data.AntY(ii), data.y_mat.dims()) - data.y_mat;
array dz = tile(data.AntZ(ii), data.x_mat.dims()) - data.z_mat;
array dR = sqrt(dx*dx + dy*dy + dz*dz) - tile(data.R0(ii), dx.dims());
// calculate the frequency of each sample in the pulse (Hz)
array freq = tile(data.minF(ii), fspan.dims()) + fspan;
// perform the Matched Filter operation
gfor(seq jj, K) {
im_slices(span,span,jj) = tile(data.phdata(jj,ii), dR.dims())*
exp(unit*(double)(4.0*Pi/c)*tile(freq(jj), dR.dims())*dR);
}
// sum over the last dimension
data.im_final += moddims(sum(im_slices,2),data.x_mat.dims());
}
}