Since the 1950s Synthetic aperture radar (SAR) systems have gained extreme popularity in both civilian and military domains due to their all-weather, day-or-night capabilities as well as the ability to render different views of a “target”. However, the raw SAR data (phase-history data) must be preprocessed since all point targets at each pulse instance are superimposed and create a complex interference that is not very useful for target location. SAR image formation algorithms compress this target information in range (frequency) and along-track (azimuth) directions to obtain interpretable images.
In the paper titled “SAR image formation toolbox for MATLAB®“, Gorham L.A. and Moore L.J. of the Air Force Research Lab discuss the implementation of the matched filter and backprojection image formation algorithms in MATLAB®. We demonstrate how easily these algorithms can be accelerated on the GPU using Jacket. For a data set, we have taken November 96 Collection X-band SAR imagery (MSTAR data) available for download at SDMS Public Web Site. Each data set consists of complex phase-history data for 128 azimuth points and 128 frequency samples. Range and cross-range resolutions have been set to 0.304 m. The data is HH-polarimetric with a frequency bandwidth of 0.591 GHz at a center frequency of 9.6 GHz. Two reconstructed SAR images (using both algorithms) are shown below:
Matched filter (left) and backprojection (right) images using MSTAR dataset
Matched filter (left) and backprojection (right) images using MSTAR dataset
We have conducted experiments on a parallel cluster with the following specs:
- 4x Intel(R) Xeon(R) CPU E7-4860 (40 cores / 80 threads via Hyper-Threading)
- 256 GB RAM (32 * 8GB)
- 4x NVIDIA Tesla S2050 (GF100)
running MATLAB® R2012a under 64-bit Linux platform. We have run both reconstruction algorithms on the same data set for increasing resolutions:
Note that, for the matched filter algorithm (left diagram), the default MATLAB® implementation can benefit from the large number of CPU cores available (that is why the running time does not grow linearly with increasing resolution). But, even with this, Jacket achieves about 20-50x speed-up over the MATLAB® CPU implementation, while ArrayFire achieves even larger speed-up.
The code that was used to produce these results is shared below.
Here is the implementation for both of these SAR algorithms with ArrayFire:
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
};
void matched_filter_af(SAR_data& data) {
double c = 299792458.0; // speed of light
// Determine the size of the phase history data
int K = data.phdata.dims(0); //! The number of frequency bins per pulse
int Np = data.phdata.dims(1); //! The number of pulses
//! Initialize the image with all zero values (complex)
data.im_final = zeros(data.x_mat.dims(), c64);
array im_slices = zeros(K, data.x_mat.dims(0),
data.x_mat.dims(1), c64);
array fspan = array(seq(0.0, K-1)) * data.deltaF;
for(int ii = 0; ii < Np; ii++) {
//! Calculate differential range for each pixel in the image (m)
array dx = data.AntX(ii)-data.x_mat,
dy = data.AntY(ii)-data.y_mat,
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) { // parallel gfor loop
im_slices(jj,span,span) =
tt(jj)*exp(i*(4.0*Pi/c*freq(jj)*dR));
}
tt = sum(im_slices,0);
data.im_final = data.im_final + reshape(tt, data.x_mat.dims());
}
}
void backprojection_af(SAR_data& data) {
double c = 299792458.0; // speed of light
// Determine the size of the phase history data
int K = data.phdata.dims(0); //! The number of frequency bins per pulse
int Np = data.phdata.dims(1); //! The number of pulses
int Nfft = K * 10; // set FFT size roughly 10 times larger than K
//! Initialize the image with all zero values (complex)
data.im_final = zeros(data.x_mat.dims(), c64);
//! Calculate the range to every bin in the range profile (m)
array r_vec = array(seq(-Nfft/2, Nfft/2-1)) *data.maxWr/Nfft;
for(int ii = 0; ii < Np; ii++) {
//! Form the range profile with zero padding added
array rc = shift(ifft(data.phdata(span,ii), Nfft), Nfft/2);
//! Calculate differential range for each pixel in the image (m)
array dx = data.AntX(ii)-data.x_mat,
dy = data.AntY(ii)-data.y_mat,
dz = data.AntZ(ii)-data.z_mat;
array dR = sqrt(dx*dx + dy*dy + dz*dz) - data.R0(ii);
//! Calculate phase correction for image
array phCorr = exp(i*(4.0*Pi*data.minF(ii)/c*dR));
//! Determine which pixels fall within the range swath
array I = where((dR > min(r_vec)) && (dR < max(r_vec)));
//! Update the image using linear interpolation
array slice = dR(I).copy();
array res = interp(r_vec, rc, slice, af_interp_linear, af::NaN);
data.im_final(I) = data.im_final(I) + res * phCorr(I);
}
}
Here is the implementation for both of these SAR algorithms with Jacket:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% For both algorithms, the following fields need to be populated:
%
% data.deltaF (scalar): Step size of frequency data (Hz)
% data.minF (Np x 1): Vector containing the start frequency of each pulse (Hz) %
% data.x_mat (Sx x Sy): The x-position of each pixel (m)
% data.y_mat (Sx x Sy): The y-position of each pixel (m)
% data.z_mat (Sx x Sy): The z-position of each pixel (m)
% data.AntX (Np x 1): The x-position of the sensor at each pulse (m)
% data.AntY (Np x 1): The y-position of the sensor at each pulse (m)
% data.AntZ (Np x 1): The z-position of the sensor at each pulse (m)
% data.R0 (Np x 1): The range to scene center (m)
% data.phdata (K x Np): Phase history data (frequency domain)
% Fast time in rows, slow time in columns
% The output is:
% data.im_final (Sx x Sy): The complex image value at each pixel
% data.maxWr (scalar): maximum scene size in range direction (m)
% data.Nfft (scalar): size of the FFT to form the range profile
% (backprojection only)
%
% MATLAB," Algorithms for Synthetic Aperture Radar Imagery XVII
% 7669, SPIE (2010).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 'data' fields must be cast to GPU (gdouble) before calling this function
function data = matched_filter_gpu(data)
c = 299792458; % Define speed of light (m/s)
% Determine the size of the phase history data
data.K = size(data.phdata,1); % The number of frequency bins per pulse
data.Np = size(data.phdata,2); % The number of pulses
% Initialize the image with all zero values
data.im_final = gzeros(size(data.x_mat), 'gdouble');
data.im_slices = gzeros([data.K,size(data.x_mat)], 'gdouble');
% Loop through every pulse
for ii = 1:data.Np
% Calculate differential range for each pixel in the image (m)
dR = sqrt((data.AntX(ii)-data.x_mat).^2 + ...
(data.AntY(ii)-data.y_mat).^2 + ...
(data.AntZ(ii)-data.z_mat).^2) - data.R0(ii);
% Calculate the frequency of each sample in the pulse (Hz)
freq = data.minF(ii) + (0:(data.K-1)) * data.deltaF;
% Perform the Matched Filter operation
gfor jj = 1:data.K
data.im_slices(jj,:,:) = data.phdata(jj,ii) * exp(1i*4*pi*freq(jj)/c*dR);
gend
data.im_final = data.im_final + squeeze(sum(data.im_slices,1));
end
end
% 'data' fields must be cast to GPU (gdouble) before calling this function
function data = backprojection_gpu(data)
c = 299792458; % Define speed of light (m/s)
% Determine the size of the phase history data
data.K = size(data.phdata,1); % The number of frequency bins per pulse
data.Np = size(data.phdata,2); % The number of pulses
% Calculate the range to every bin in the range profile (m)
data.r_vec = linspace(-data.Nfft/2,data.Nfft/2-1,data.Nfft)*data.maxWr/data.Nfft;
% Initialize the image with all zero values
data.im_final = gzeros(size(data.x_mat), 'gdouble');
% Loop through every pulse
for ii = 1:data.Np
% Form the range profile with zero padding added
rc = fftshift(ifft(data.phdata(:,ii), data.Nfft));
% Calculate differential range for each pixel in the image (m)
dR = sqrt((data.AntX(ii)-data.x_mat).^2 + ...
(data.AntY(ii)-data.y_mat).^2 + ...
(data.AntZ(ii)-data.z_mat).^2) - data.R0(ii);
% Calculate phase correction for image
phCorr = exp(1i*4*pi*data.minF(ii)/c*dR);
% Determine which pixels fall within the range swath
I = find(and(dR > min(data.r_vec), dR < max(data.r_vec)));
% Update the image using linear interpolation
data.im_final(I) = data.im_final(I) + interp1(data.r_vec,rc,dR(I),'linear') .* phCorr(I);
end
end
Comments 1
Pingback: Walking Randomly » A Month of Math Software – August 2012