/*
* Copyright (C) 2011-2015 Karlsruhe Institute of Technology
*
* This file is part of Ufo.
*
* This library is free software: you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation, either
* version 3 of the License, or (at your option) any later version.
*
* This library is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library. If not, see .
*
* Authored by: Alexandre Lewkowicz (lewkow_a@epita.fr)
*/
#include "config.h"
#include
#include
#include
#include "ufo-filter-particle-task.h"
#include "ufo-ring-coordinates.h"
struct _UfoFilterParticleTaskPrivate {
unsigned *img;
size_t num_pixels;
float threshold;
float min;
};
/*Maximal number of particles that can be tagged*/
#define nVectMax 100000
static void ufo_task_interface_init (UfoTaskIface *iface);
G_DEFINE_TYPE_WITH_CODE (UfoFilterParticleTask, ufo_filter_particle_task, UFO_TYPE_TASK_NODE,
G_IMPLEMENT_INTERFACE (UFO_TYPE_TASK,
ufo_task_interface_init))
#define UFO_FILTER_PARTICLE_TASK_GET_PRIVATE(obj) (G_TYPE_INSTANCE_GET_PRIVATE((obj), UFO_TYPE_FILTER_PARTICLE_TASK, UfoFilterParticleTaskPrivate))
enum {
PROP_0,
PROP_MIN,
PROP_THRESHOLD,
N_PROPERTIES
};
static GParamSpec *properties[N_PROPERTIES] = { NULL, };
UfoNode *
ufo_filter_particle_task_new (void)
{
return UFO_NODE (g_object_new (UFO_TYPE_FILTER_PARTICLE_TASK, NULL));
}
static void
ufo_filter_particle_task_setup (UfoTask *task,
UfoResources *resources,
GError **error)
{
}
static void
ufo_filter_particle_task_get_requisition (UfoTask *task,
UfoBuffer **inputs,
UfoRequisition *requisition,
GError **error)
{
UfoFilterParticleTaskPrivate *priv = UFO_FILTER_PARTICLE_TASK_GET_PRIVATE (task);
requisition->n_dims = 1;
ufo_buffer_get_requisition(inputs[0], requisition);
size_t num_pixels = requisition->dims[0] * requisition->dims[1];
/* Allocate internal buffer to hold copy of image */
/* Reuse buffer when size is the smae */
if (priv->num_pixels != num_pixels) {
priv->num_pixels = num_pixels;
priv->img = realloc (priv->img, sizeof (unsigned) * num_pixels);
}
/* Output size varies according to data whithin the input */
/* I don't know the output size, I'll allocate it myself once I have
* processed the data */
requisition->dims[0] = 0;
requisition->n_dims = 1;
}
static guint
ufo_filter_particle_task_get_num_inputs (UfoTask *task)
{
return 1;
}
static guint
ufo_filter_particle_task_get_num_dimensions (UfoTask *task,
guint input)
{
return 2;
}
static UfoTaskMode
ufo_filter_particle_task_get_mode (UfoTask *task)
{
return UFO_TASK_MODE_PROCESSOR | UFO_TASK_MODE_CPU;
}
static void
part_cen_radius(float *img, unsigned *bin, unsigned nc, unsigned nr,
unsigned pc, double *cx, double *cy, double *rr)
{
/* nr == number of rows */
/* nc == number of columns */
/* pc == pixel count, (i-e nb of detected rings) */
/* cx, cy, rr, (x, y, r) coordinates of each ring */
unsigned i,j,k, itop, ibot, jleft, jright;
unsigned b;
unsigned *tops,*bots,*lefts,*rights;
float cenx, radx, ceny, rady, wei;
tops = malloc((pc+1)*sizeof(unsigned) * 100);
bots = malloc((pc+1)*sizeof(unsigned) * 100);
lefts = malloc((pc+1)*sizeof(unsigned) * 100);
rights = malloc((pc+1)*sizeof(unsigned) * 100);
for (k=0;k<=pc;k++) {
tops[k]= (unsigned)nr; bots[k]=0; lefts[k]=(unsigned)nc; rights[k]=0;
}
/* Find extremities of each clustered pixels */
for (i=0;ii ) tops[b]=i;
if ( bots[b]j ) lefts[b]=j;
if ( rights[b] threshold) {
/* Bottom left pixel */
if (i < (nr - 1) && j > 0 && input[(i + 1) * nc + (j - 1)] > threshold)
{ ncase += 1; iodd = 1; }
/* Bottom right pixel */
if (i < (nr - 1) && j < (nc - 1) && input[(i + 1) * nc + (j + 1)] > threshold)
ncase += 2;
/* Top right pixel */
if (i > 0 && j < (nc - 1) && input[(i - 1) * nc + (j + 1)] > threshold)
ncase += 4;
/* Top left pixel */
if (i > 0 && j > 0 && input[(i - 1) * nc + (j - 1)] > threshold)
ncase += 8;
/* Bottom pixel */
if (i < (nr - 1) && input[(i + 1) * nc + j] > threshold)
ncase += 16;
/* Right pixel */
if (j < (nc - 1) && input[i * nc + (j + 1)] > threshold)
ncase += 32;
/* Top pixel */
if (i > 0 && input[(i - 1) * nc + j] > threshold)
ncase += 64;
/* Left pixel */
if (j > 0 && input[i * nc + (j - 1)] > threshold)
ncase += 128;
/* Pixel has no : Bottom, Right, Top, Left */
if (ncase <= 15) {
// if (NO1 == 1)
// bin[i][j] = ++pc;
// else
out[i * nc + j] = 0;
}
/* Pixel has no : Right, Top, Left */
/* and Pixel has : Bottom */
else if (ncase <= 31) {
if(iodd == 1)
out[i * nc + j] = out[(i + 1) * nc + (j - 1)]; /* Clone id */
else
out[i * nc + j] = ++pc;
}
/* Pixel has no : Top, Left, Bottom */
/* Pixel has : Right */
else if (ncase <= 47)
out[i * nc + j] = ++pc;
/* Pixel has no : Top, Left */
/* Pixel has : Right, Bottom */
else if (ncase <= 63) {
if (iodd == 1)
out[i * nc + j] = out[(i + 1) * nc + (j - 1)];
else
out[i * nc + j] = ++pc;
}
/* Pixel has no : Left */
/* Pixel has : Top */
else if (ncase <= 127)
out[i * nc + j] = out[(i - 1) * nc + j];
/* Pixel has : Left */
else
out[i * nc + j] = out[i * nc + (j - 1)];
/* Bottom Left && Bottom && Top && not(Left) */
if ((ncase & 1) && (ncase & 16) && (ncase & 64) && !(ncase & 128)) {
/* Find duplicates, and locate them by x, and y list */
/* These duplicates shall be joined later */
k1 = out[i * nc + j];
/* 11
* 1
* 22
* in this case, bottom left was processed later, which
* explains the comparison */
k2 = out[(i + 1) * nc + (j - 1)];
dupx[ndup] = (k1 > k2 ? k2 : k1);
dupy[ndup] = (k1 > k2 ? k1 : k2);
ndup++;
}
}
else
out[i * nc + j] = 0;
}
}
dupindex = malloc ((unsigned) (pc + 1) * sizeof (unsigned));
/* binloc stores the 1D position (on row basis, like matlab) of each
* candidate block of pixel (Most left-top of block)
* or MAX_POS_VAL otherwise */
binloc = malloc ((unsigned) (pc + 1) * sizeof(unsigned));
for(k=0;k<=pc;k++) {
binloc[k]= nr * nc;
dupindex[k] = k;
}
for(j=0;j 0 && binloc[b] > (j*nr+i) )
binloc[b]=j*nr+i;
}
}
for(k=0; k < ndup; k++) {
/* 11
* 1
* 22 in this case*/
if( binloc[dupx[k]] > binloc[dupy[k]] ) {
/* Dupx has always been afected the smaller id
* ==> with row wise basis, locations dupx > dupy is always false */
temp = dupx[k]; dupx[k]=dupy[k]; dupy[k]=temp;
}
}
/* Sort dupx to be in crescent order */
/* i-e
* 22
* 2
* 33
* 3
*111
* Here dupx = { 2, 1 } ==> Not in crescent order */
if (ndup != 0) {
for (k=0;kbinloc[dupx[j]]) {
temp=dupx[k];
dupx[k]=dupx[j];
dupx[j]=temp;
temp=dupy[k];
dupy[k]=dupy[j];
dupy[j]=temp;
}
else if( dupx[k]== dupx[j] && binloc[dupy[k]] > binloc[dupy[j]]) {
temp=dupy[k];
dupy[k]=dupy[j];
dupy[j]=temp;
}
}
}
}
/* dupindex gives the new index for each index, when two index are duplicates
* dupindex will affect those two indexes a same value */
for(k=0;kdupx[j]) {
temp=dupx[k];
dupx[k]=dupx[j];
dupx[j]=temp;
}
}
}
ncc=0;
dupy[0]=0;
/* map dupindex values to values from 1 to number_of_detected_ring, this
* allows to have a constant step of one, eg: converting {1,3,6} -> {1,2,3} */
for(k=1;k<=pc;k++) {
if(dupx[k]>dupx[k-1])
dupy[ dupx[k] ] = ++ncc;
}
*pcount=ncc;
/* set dupindex to its newly mapped value, (making dupindex have a step of
* one */
for (k=0;k<=pc;k++)
dupindex[k]= dupy[ dupindex[k] ];
/* Make all duplicates have the same ID */
for(i=0;i b ? b : a;
}
static gboolean
ufo_filter_particle_task_process (UfoTask *task,
UfoBuffer **inputs,
UfoBuffer *output,
UfoRequisition *requisition)
{
(void) requisition;
UfoFilterParticleTaskPrivate *priv = UFO_FILTER_PARTICLE_TASK_GET_PRIVATE (task);
UfoRequisition req;
ufo_buffer_get_requisition(inputs[0], &req);
unsigned pcount = 0;
float *input = ufo_buffer_get_host_array(inputs[0], NULL);
unsigned radius;
get_radius_metadata(inputs[0], &radius);
/* Give unique ID to each cluster of pixels */
double max_value = get_max_value(input, req.dims[0], req.dims[1]);
double threshold = priv->threshold * max_value;
// When max_value is too low, it is unlikely for there to be any rings
if (max_value < priv->min) {
g_print("Filter_particle: min : Ignoring radius %u. Inputs max value is"
" %f, at least %f expected\n", radius, max_value, priv->min);
UfoRequisition out_req = {
.n_dims = 1,
.dims[0] = 1
};
ufo_buffer_resize (output, &out_req);
float *out = ufo_buffer_get_host_array(output, NULL);
*out = 0;
return TRUE;
}
pfind (input, priv->img, (unsigned) req.dims[0], (unsigned) req.dims[1], &pcount, threshold);
/* Compute center of each cluster and store coordinates in cx, cy */
double *cx = malloc (sizeof (double) * pcount);
double *cy = malloc (sizeof (double) * pcount);
double *rr = malloc (sizeof (double) * pcount);
part_cen_radius(input, priv->img, (unsigned) req.dims[0],
(unsigned) req.dims[1], pcount, cx, cy, rr);
/* a list of pcount UfoRingCoordinate are generated in output */
/* first element is the number of rings detected */
UfoRequisition out_req = {
.n_dims = 1,
.dims[0] = 1 + pcount * sizeof (UfoRingCoordinate) / sizeof (float),
};
ufo_buffer_resize (output, &out_req);
float *out = ufo_buffer_get_host_array(output, NULL);
UfoRingCoordinate *res = (UfoRingCoordinate *) (out + 1);
*out = (float) pcount;
for (unsigned i = 0; i < pcount; ++i) {
res[i].x = (float) cx[i];
res[i].y = (float) cy[i];
res[i].r = (float) radius;
unsigned x = min((unsigned)round(cx[i]), (unsigned) (req.dims[0] - 1));
unsigned y = min((unsigned)round(cy[i]), (unsigned) (req.dims[1] - 1));
res[i].intensity = input[y * req.dims[0] + x];
}
free (cx);
free (cy);
free (rr);
return TRUE;
}
static void
ufo_filter_particle_task_set_property (GObject *object,
guint property_id,
const GValue *value,
GParamSpec *pspec)
{
UfoFilterParticleTaskPrivate *priv = UFO_FILTER_PARTICLE_TASK_GET_PRIVATE (object);
switch (property_id) {
case PROP_MIN:
priv->min = g_value_get_float(value);
break;
case PROP_THRESHOLD:
priv->threshold = g_value_get_float(value);
break;
default:
G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
break;
}
}
static void
ufo_filter_particle_task_get_property (GObject *object,
guint property_id,
GValue *value,
GParamSpec *pspec)
{
UfoFilterParticleTaskPrivate *priv = UFO_FILTER_PARTICLE_TASK_GET_PRIVATE (object);
switch (property_id) {
case PROP_MIN:
g_value_set_float(value, priv->min);
break;
case PROP_THRESHOLD:
g_value_set_float(value, priv->threshold);
break;
default:
G_OBJECT_WARN_INVALID_PROPERTY_ID (object, property_id, pspec);
break;
}
}
static void
ufo_filter_particle_task_finalize (GObject *object)
{
G_OBJECT_CLASS (ufo_filter_particle_task_parent_class)->finalize (object);
UfoFilterParticleTaskPrivate *priv = UFO_FILTER_PARTICLE_TASK_GET_PRIVATE (object);
free (priv->img);
}
static void
ufo_task_interface_init (UfoTaskIface *iface)
{
iface->setup = ufo_filter_particle_task_setup;
iface->get_num_inputs = ufo_filter_particle_task_get_num_inputs;
iface->get_num_dimensions = ufo_filter_particle_task_get_num_dimensions;
iface->get_mode = ufo_filter_particle_task_get_mode;
iface->get_requisition = ufo_filter_particle_task_get_requisition;
iface->process = ufo_filter_particle_task_process;
}
static void
ufo_filter_particle_task_class_init (UfoFilterParticleTaskClass *klass)
{
GObjectClass *gobject_class = G_OBJECT_CLASS (klass);
gobject_class->set_property = ufo_filter_particle_task_set_property;
gobject_class->get_property = ufo_filter_particle_task_get_property;
gobject_class->finalize = ufo_filter_particle_task_finalize;
properties[PROP_MIN] =
g_param_spec_float ("min",
"Ignore image when it's maximal value is less than min",
"Ignore image when it's maximal value is less than min",
0, 1, 0.125f,
G_PARAM_READWRITE);
properties[PROP_THRESHOLD] =
g_param_spec_float ("threshold",
"Set threshold value which is relative to the image maximum value",
"Set threshold value which is relative to the image maximum value",
0, 1, 0.8f,
G_PARAM_READWRITE);
for (guint i = PROP_0 + 1; i < N_PROPERTIES; i++)
g_object_class_install_property (gobject_class, i, properties[i]);
g_type_class_add_private (gobject_class, sizeof(UfoFilterParticleTaskPrivate));
}
static void
ufo_filter_particle_task_init(UfoFilterParticleTask *self)
{
self->priv = UFO_FILTER_PARTICLE_TASK_GET_PRIVATE(self);
self->priv->img = NULL;
self->priv->num_pixels = 0;
self->priv->min = 0.125f;
self->priv->threshold = 0.8f;
}