// Sample particle filter construction, input, and output using ArrayFire
//
// Note: array data is stored on the device and is marked as free when it is out of scope
// for this reason you do not need to manually deallocate device data
#define _USE_MATH_DEFINES
#define FILENAME_BUFFER_SIZE 25
#include <arrayfire.h>
#include <cstdio>
#include <cstdlib>
#include <math.h>
#include <time.h>
using namespace af;
int main(int argc, char *argv[]) {
try {
// Select a device and display arrayfire info
int device = argc > 1 ? atoi(argv[1]) : 0;
setDevice(device);
info();
// Set seed for af::randn
setSeed(time(NULL));
// Image and filter variables
int tempSize = 39;
int patchSize = (tempSize - 1) / 2;
int resY = 720;
int resX = 1280;
int padResY = resY + tempSize;
int padResX = resX + tempSize;
int startFrame = 100;
int numFrames = 50;
int numObjects = 1;
int numParticles = 500;
// Distribution
float variance = 4;
float sigma = 10000;
float mu = 0;
// Array of images in set to be tracked
array imArray = constant(0, padResY, padResX, numFrames);
char fileName[FILENAME_BUFFER_SIZE];
for (int f = startFrame; f < startFrame + numFrames; f++) {
array pad = constant(255, padResY, padResX);
snprintf(fileName, FILENAME_BUFFER_SIZE, "Input/frame_%i.png", f);
pad(seq(patchSize, resY + (patchSize - 1)), seq(patchSize, resX + (patchSize - 1))) = rgb2gray(loadImageNative(fileName));
imArray(span, span, f - startFrame) = pad;
}
// Ground Truth
int hostGT[] = { 608 + patchSize, 532 + patchSize };
array deviceGT = constant(0, numObjects, 2, s32);
deviceGT(0, 1) = hostGT[0];
deviceGT(0, 0) = hostGT[1];
// Create vector of template image
array templateArray = constant(255, tempSize, tempSize, numObjects);
templateArray(span, span, 0) = imArray(seq(hostGT[1] - patchSize, hostGT[1] + patchSize), seq(hostGT[0] - patchSize, hostGT[0] + patchSize), 0);
array templateVector = tile(templateArray(span, span, 0), 1, 1, numParticles);
// Best particle
array best = constant(0, numFrames, 2, numObjects);
// Linear model of image patch
array model = (iota(dim4(tempSize, 1), dim4(1, tempSize)) - patchSize);
model += tile((range(dim4(1, tempSize), 1) - patchSize) * padResY, tempSize);
// Init particles randomly from ground truth
array particles = constant(0, numParticles, 3, numObjects);
array noise = randn(numParticles, 2) * variance;
particles(span, seq(2), 0) = round(tile(deviceGT(0, seq(2)), numParticles) + noise);
particles(span, 2, 0) = 0;
// Frame loop
for (int f = 0; f < numFrames; f++) {
array im = imArray(span, span, f);
// Create translate matrix for each particle, declare vector for each particles patch
array translate = tile(moddims(particles(span, 0, 0) + (particles(span, 1, 0) * padResY), 1, 1, numParticles), tempSize, tempSize);
array patchVector = constant(0, tempSize, tempSize, numParticles);
// Collect patches for each particle
gfor(seq n, numParticles) {
patchVector = moddims(im(model + translate(span, span, n)), tempSize, tempSize, numParticles);
}
// Sum of Absolute Differences
array sad = moddims(sum(sum(abs(patchVector - templateVector)), 1), numParticles);
// Probability that the object is at that particle
array probs = (1 / (sqrt(2 * M_PI)*sigma) * exp(-(sad(span, 0) * sad(span, 0)) / (2 * (sigma*sigma))));
particles(span, 2, 0) = probs(span, 0) / tile(sum(probs(span, 0), 0), numParticles);
// Find maximum of the particle probabilities
// val here is need for the max function to return, as I don't
array val = constant(0, 1);
array idx = constant(0, 1);
max(val, idx, particles(span, 2, 0));
best(f, span, 0) = particles(idx, seq(2), 0);
// Redistribute particles randomly
particles(span, seq(2), 0) = round(tile(particles(idx, seq(2), 0), numParticles) + variance * randn(numParticles, 2));
// Video output for example only
// As hard disk R/W is a bottleneck, ideally this program would only return a matrix of object locations
// If video is needed it is better suited on its own
snprintf(fileName, FILENAME_BUFFER_SIZE, "Input/frame_%i.png", f+100);
array imOut = loadImageNative(fileName);
for (int n = 0; n < numParticles; n++) {
imOut(*particles(n, 0, 0).host<float>() - patchSize, *particles(n, 1, 0).host<float>() - patchSize, 0) = 255;
imOut(*particles(n, 0, 0).host<float>() - patchSize, *particles(n, 1, 0).host<float>() - patchSize, 1) = 255;
imOut(*particles(n, 0, 0).host<float>() - patchSize, *particles(n, 1, 0).host<float>() - patchSize, 2) = 0;
}
snprintf(fileName, FILENAME_BUFFER_SIZE, "Output/frame_%i.png", f);
saveImageNative(fileName, imOut);
} // Frame Loop
} catch (exception& e) {
fprintf(stderr, "%s\n", e.what());
throw;
}
return 0;
