Paparazzi UAS  v5.15_devel-113-g1b57ff1
Paparazzi is a free software Unmanned Aircraft System.
 All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Modules Pages
linear_flow_fit.c
Go to the documentation of this file.
1 /*
2  * Copyright (C) 2014 Hann Woei Ho
3  * Guido de Croon
4  *
5  *
6  * This file is part of Paparazzi.
7  *
8  * Paparazzi is free software; you can redistribute it and/or modify
9  * it under the terms of the GNU General Public License as published by
10  * the Free Software Foundation; either version 2, or (at your option)
11  * any later version.
12  *
13  * Paparazzi is distributed in the hope that it will be useful,
14  * but WITHOUT ANY WARRANTY; without even the implied warranty of
15  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  * GNU General Public License for more details.
17  *
18  * You should have received a copy of the GNU General Public License
19  * along with Paparazzi; see the file COPYING. If not, write to
20  * the Free Software Foundation, 59 Temple Place - Suite 330,
21  * Boston, MA 02111-1307, USA.
22  */
23 
24 /*
25  * @file modules/computer_vision/opticflow/linear_flow_fit.c
26  * @brief Takes a set of optical flow vectors and extracts information such as relative velocities and surface roughness.
27  *
28  * A horizontal and vertical linear fit is made with the optical flow vectors, and from the fit parameters information is extracted such as relative velocities (useful for time-to-contact determination), slope angle, and surface roughness.
29  *
30  * Code based on the article:
31  * "Optic-flow based slope estimation for autonomous landing",
32  * de Croon, G.C.H.E., and Ho, H.W., and De Wagter, C., and van Kampen, E., and Remes B., and Chu, Q.P.,
33  * in the International Journal of Micro Air Vehicles, Volume 5, Number 4, pages 287 – 297, (2013)
34  */
35 
36 #include <stdlib.h>
37 #include <stdio.h>
38 #include <math.h>
39 #include <string.h>
40 //#include "defs_and_types.h"
41 #include "linear_flow_fit.h"
45 
46 // Is this still necessary?
47 #define MAX_COUNT_PT 50
48 
49 #define MIN_SAMPLES_FIT 3
50 
63 bool analyze_linear_flow_field(struct flow_t *vectors, int count, float error_threshold, int n_iterations, int n_samples, int im_width, int im_height, struct linear_flow_fit_info *info)
64 {
65  // Are there enough flow vectors to perform a fit?
66  if (count < MIN_SAMPLES_FIT) {
67  // return that no fit was made:
68  return false;
69  }
70 
71  // fit linear flow field:
72  float parameters_u[3], parameters_v[3], min_error_u, min_error_v;
73  fit_linear_flow_field(vectors, count, n_iterations, error_threshold, n_samples, parameters_u, parameters_v, &info->fit_error, &min_error_u, &min_error_v, &info->n_inliers_u, &info->n_inliers_v);
74 
75  // extract information from the parameters:
76  extract_information_from_parameters(parameters_u, parameters_v, im_width, im_height, info);
77 
78  // surface roughness is equal to fit error:
79  info->surface_roughness = info->fit_error;
80  info->divergence = info->relative_velocity_z;
81  float diverg = (fabsf(info->divergence) < 1E-5) ? 1E-5 : info->divergence;
82  info->time_to_contact = 1.0f / diverg;
83 
84  // return successful fit:
85  return true;
86 }
87 
103 void fit_linear_flow_field(struct flow_t *vectors, int count, float error_threshold, int n_iterations, int n_samples, float *parameters_u, float *parameters_v, float *fit_error, float *min_error_u, float *min_error_v, int *n_inliers_u, int *n_inliers_v)
104 {
105 
106  // We will solve systems of the form A x = b,
107  // where A = [nx3] matrix with entries [x, y, 1] for each optic flow location
108  // and b = [nx1] vector with either the horizontal (bu) or vertical (bv) flow.
109  // x in the system are the parameters for the horizontal (pu) or vertical (pv) flow field.
110 
111  // local vars for iterating, random numbers:
112  int sam, p, i_rand, si, add_si;
113 
114  // ensure that n_samples is high enough to ensure a result for a single fit:
115  n_samples = (n_samples < MIN_SAMPLES_FIT) ? MIN_SAMPLES_FIT : n_samples;
116  // n_samples should not be higher than count:
117  n_samples = (n_samples < count) ? n_samples : count;
118 
119  // initialize matrices and vectors for the full point set problem:
120  // this is used for determining inliers
121  float _AA[count][3];
122  MAKE_MATRIX_PTR(AA, _AA, count);
123  float _bu_all[count][1];
124  MAKE_MATRIX_PTR(bu_all, _bu_all, count);
125  float _bv_all[count][1];
126  MAKE_MATRIX_PTR(bv_all, _bv_all, count);
127  for (sam = 0; sam < count; sam++) {
128  AA[sam][0] = (float) vectors[sam].pos.x;
129  AA[sam][1] = (float) vectors[sam].pos.y;
130  AA[sam][2] = 1.0f;
131  bu_all[sam][0] = (float) vectors[sam].flow_x;
132  bv_all[sam][0] = (float) vectors[sam].flow_y;
133  }
134 
135  // later used to determine the error of a set of parameters on the whole set:
136  float _bb[count][1];
137  MAKE_MATRIX_PTR(bb, _bb, count);
138  float _C[count][1];
139  MAKE_MATRIX_PTR(C, _C, count);
140 
141  // ***************
142  // perform RANSAC:
143  // ***************
144 
145  // set up variables for small linear system solved repeatedly inside RANSAC:
146  float _A[n_samples][3];
147  MAKE_MATRIX_PTR(A, _A, n_samples);
148  float _bu[n_samples][1];
149  MAKE_MATRIX_PTR(bu, _bu, n_samples);
150  float _bv[n_samples][1];
151  MAKE_MATRIX_PTR(bv, _bv, n_samples);
152  float w[n_samples], _v[3][3];
153  MAKE_MATRIX_PTR(v, _v, 3);
154  float _pu[3][1];
155  MAKE_MATRIX_PTR(pu, _pu, 3);
156  float _pv[3][1];
157  MAKE_MATRIX_PTR(pv, _pv, 3);
158 
159  // iterate and store parameters, errors, inliers per fit:
160  float PU[n_iterations * 3];
161  float PV[n_iterations * 3];
162  float errors_pu[n_iterations];
163  errors_pu[0] = 0.0;
164  float errors_pv[n_iterations];
165  errors_pv[0] = 0.0;
166  int n_inliers_pu[n_iterations];
167  int n_inliers_pv[n_iterations];
168  int it, ii;
169  for (it = 0; it < n_iterations; it++) {
170  // select a random sample of n_sample points:
171  int sample_indices[n_samples];
172  i_rand = 0;
173 
174  // sampling without replacement:
175  while (i_rand < n_samples) {
176  si = rand() % count;
177  add_si = 1;
178  for (ii = 0; ii < i_rand; ii++) {
179  if (sample_indices[ii] == si) { add_si = 0; }
180  }
181  if (add_si) {
182  sample_indices[i_rand] = si;
183  i_rand ++;
184  }
185  }
186 
187  // Setup the system:
188  for (sam = 0; sam < n_samples; sam++) {
189  A[sam][0] = (float) vectors[sample_indices[sam]].pos.x;
190  A[sam][1] = (float) vectors[sample_indices[sam]].pos.y;
191  A[sam][2] = 1.0f;
192  bu[sam][0] = (float) vectors[sample_indices[sam]].flow_x;
193  bv[sam][0] = (float) vectors[sample_indices[sam]].flow_y;
194  //printf("%d,%d,%d,%d,%d\n",A[sam][0],A[sam][1],A[sam][2],bu[sam][0],bv[sam][0]);
195  }
196 
197  // Solve the small system:
198 
199  // for horizontal flow:
200  // decompose A in u, w, v with singular value decomposition A = u * w * vT.
201  // u replaces A as output:
202  pprz_svd_float(A, w, v, n_samples, 3);
203  pprz_svd_solve_float(pu, A, w, v, bu, n_samples, 3, 1);
204  PU[it * 3] = pu[0][0];
205  PU[it * 3 + 1] = pu[1][0];
206  PU[it * 3 + 2] = pu[2][0];
207 
208  // for vertical flow:
209  pprz_svd_solve_float(pv, A, w, v, bv, n_samples, 3, 1);
210  PV[it * 3] = pv[0][0];
211  PV[it * 3 + 1] = pv[1][0];
212  PV[it * 3 + 2] = pv[2][0];
213 
214  // count inliers and determine their error on all points:
215  errors_pu[it] = 0;
216  errors_pv[it] = 0;
217  n_inliers_pu[it] = 0;
218  n_inliers_pv[it] = 0;
219 
220  // for horizontal flow:
221  // bb = AA * pu:
222  MAT_MUL(count, 3, 1, bb, AA, pu);
223  // subtract bu_all: C = 0 in case of perfect fit:
224  MAT_SUB(count, 1, C, bb, bu_all);
225 
226  for (p = 0; p < count; p++) {
227  C[p][0] = abs(C[p][0]);
228  if (C[p][0] < error_threshold) {
229  errors_pu[it] += C[p][0];
230  n_inliers_pu[it]++;
231  } else {
232  errors_pu[it] += error_threshold;
233  }
234  }
235 
236  // for vertical flow:
237  // bb = AA * pv:
238  MAT_MUL(count, 3, 1, bb, AA, pv);
239  // subtract bv_all: C = 0 in case of perfect fit:
240  MAT_SUB(count, 1, C, bb, bv_all);
241 
242  for (p = 0; p < count; p++) {
243  C[p][0] = abs(C[p][0]);
244  if (C[p][0] < error_threshold) {
245  errors_pv[it] += C[p][0];
246  n_inliers_pv[it]++;
247  } else {
248  errors_pv[it] += error_threshold;
249  }
250  }
251  }
252 
253  // After all iterations:
254  // select the parameters with lowest error:
255  // for horizontal flow:
256  int param;
257  int min_ind = 0;
258  *min_error_u = (float)errors_pu[0];
259  for (it = 1; it < n_iterations; it++) {
260  if (errors_pu[it] < *min_error_u) {
261  *min_error_u = (float)errors_pu[it];
262  min_ind = it;
263  }
264  }
265  for (param = 0; param < 3; param++) {
266  parameters_u[param] = PU[min_ind * 3 + param];
267  }
268  *n_inliers_u = n_inliers_pu[min_ind];
269 
270  // for vertical flow:
271  min_ind = 0;
272  *min_error_v = (float)errors_pv[0];
273  for (it = 0; it < n_iterations; it++) {
274  if (errors_pv[it] < *min_error_v) {
275  *min_error_v = (float)errors_pv[it];
276  min_ind = it;
277  }
278  }
279  for (param = 0; param < 3; param++) {
280  parameters_v[param] = PV[min_ind * 3 + param];
281  }
282  *n_inliers_v = n_inliers_pv[min_ind];
283 
284  // error has to be determined on the entire set without threshold:
285  // bb = AA * pu:
286  MAT_MUL(count, 3, 1, bb, AA, pu);
287  // subtract bu_all: C = 0 in case of perfect fit:
288  MAT_SUB(count, 1, C, bb, bu_all);
289  *min_error_u = 0;
290  for (p = 0; p < count; p++) {
291  *min_error_u += abs(C[p][0]);
292  }
293  // bb = AA * pv:
294  MAT_MUL(count, 3, 1, bb, AA, pv);
295  // subtract bv_all: C = 0 in case of perfect fit:
296  MAT_SUB(count, 1, C, bb, bv_all);
297  *min_error_v = 0;
298  for (p = 0; p < count; p++) {
299  *min_error_v += abs(C[p][0]);
300  }
301  *fit_error = (*min_error_u + *min_error_v) / (2 * count);
302 
303 }
313 void extract_information_from_parameters(float *parameters_u, float *parameters_v, int im_width, int im_height, struct linear_flow_fit_info *info)
314 {
315  // This method assumes a linear flow field in x- and y- direction according to the formulas:
316  // u = parameters_u[0] * x + parameters_u[1] * y + parameters_u[2]
317  // v = parameters_v[0] * x + parameters_v[1] * y + parameters_v[2]
318  // where u is the horizontal flow at image coordinate (x,y)
319  // and v is the vertical flow at image coordinate (x,y)
320 
321  // relative velocities:
322  info->relative_velocity_z = (parameters_u[0] + parameters_v[1]) / 2.0f; // divergence / 2
323 
324  // translation orthogonal to the camera axis:
325  // flow in the center of the image:
326  info->relative_velocity_x = -(parameters_u[2] + (im_width / 2.0f) * parameters_u[0] + (im_height / 2.0f) * parameters_u[1]);
327  info->relative_velocity_y = -(parameters_v[2] + (im_width / 2.0f) * parameters_v[0] + (im_height / 2.0f) * parameters_v[1]);
328 
329  float arv_x = abs(info->relative_velocity_x);
330  float arv_y = abs(info->relative_velocity_y);
331 
332  // extract inclination from flow field:
333  float threshold_slope = 1.0;
334  float eta = 0.002;
335 
336  if (abs(parameters_v[1]) < eta && arv_y < threshold_slope && arv_x >= 2 * threshold_slope) {
337  // there is no forward motion and not enough vertical motion, but enough horizontal motion:
338  info->slope_x = parameters_u[0] / info->relative_velocity_x;
339  } else if (arv_y >= 2 * threshold_slope) {
340  // there is sufficient vertical motion:
341  info->slope_x = parameters_v[0] / info->relative_velocity_y;
342  } else {
343  // there may be forward motion, then we can do a quadratic fit:
344  // a linear fit provides no information though
345  info->slope_x = 0.0f;
346  }
347 
348  if (abs(parameters_u[0]) < eta && arv_x < threshold_slope && arv_y >= 2 * threshold_slope) {
349  // there is no forward motion, little horizontal movement, but sufficient vertical motion:
350  info->slope_y = parameters_v[1] / info->relative_velocity_y;
351  } else if (arv_x >= 2 * threshold_slope) {
352  // there is sufficient horizontal motion:
353  info->slope_y = parameters_u[1] / info->relative_velocity_x;
354  } else {
355  // there could be forward motion, then we can do a quadratic fit:
356  // a linear fit provides no information though
357  info->slope_y = 0.0f;
358  }
359 
360  // Focus of Expansion:
361  // the flow planes intersect the flow=0 plane in a line
362  // the FoE is the point where these 2 lines intersect (flow = (0,0))
363  // x:
364  float denominator = parameters_v[0] * parameters_u[1] - parameters_u[0] * parameters_v[1];
365  if (abs(denominator) > 1E-5) {
366  info->focus_of_expansion_x = ((parameters_u[2] * parameters_v[1] - parameters_v[2] * parameters_u[1]) / denominator);
367  } else { info->focus_of_expansion_x = 0.0f; }
368  // y:
369  denominator = parameters_u[1];
370  if (abs(denominator) > 1E-5) {
371  info->focus_of_expansion_y = (-(parameters_u[0] * (info->focus_of_expansion_x) + parameters_u[2]) / denominator);
372  } else { info->focus_of_expansion_y = 0.0f; }
373 }
374 
float relative_velocity_y
Relative velocity in y-direction, i.e., vy / z, where z is the depth in direction of the camera's pri...
void fit_linear_flow_field(struct flow_t *vectors, int count, float error_threshold, int n_iterations, int n_samples, float *parameters_u, float *parameters_v, float *fit_error, float *min_error_u, float *min_error_v, int *n_inliers_u, int *n_inliers_v)
Analyze a linear flow field, retrieving information such as divergence, surface roughness, focus of expansion, etc.
Simple matrix helper macros.
#define MAKE_MATRIX_PTR(_ptr, _mat, _rows)
Make a pointer to a matrix of _rows lines.
float fit_error
Error of the fit (same as surface roughness)
int n_inliers_u
Number of inliers in the horizontal flow fit.
Definition: image.h:66
float relative_velocity_z
Relative velocity in z-direction, i.e., vz / z, where z is the depth in direction of the camera's pri...
float time_to_contact
Basically, 1 / relative_velocity_z.
int n_samples
Definition: detect_gate.c:85
float focus_of_expansion_x
Image x-coordinate of the focus of expansion (contraction)
float focus_of_expansion_y
Image y-coordinate of the focus of expansion (contraction)
Paparazzi floating point algebra.
bool analyze_linear_flow_field(struct flow_t *vectors, int count, float error_threshold, int n_iterations, int n_samples, int im_width, int im_height, struct linear_flow_fit_info *info)
Analyze a linear flow field, retrieving information such as divergence, surface roughness, focus of expansion, etc.
float slope_y
Slope of the surface in y-direction - given sufficient lateral motion.
void pprz_svd_solve_float(float **x, float **u, float *w, float **v, float **b, int m, int n, int l)
SVD based linear solver.
float relative_velocity_x
Relative velocity in x-direction, i.e., vx / z, where z is the depth in direction of the camera's pri...
#define AA
Definition: LPC21xx.h:179
struct point_t pos
The original position the flow comes from in subpixels.
Definition: image.h:67
float slope_x
Slope of the surface in x-direction - given sufficient lateral motion.
uint32_t y
The y coordinate of the point.
Definition: image.h:59
void extract_information_from_parameters(float *parameters_u, float *parameters_v, int im_width, int im_height, struct linear_flow_fit_info *info)
Extract information from the parameters that were fit to the optical flow field.
Matrix decompositions in floating point.
int n_inliers_v
Number of inliers in the vertical flow fit.
float divergence
Basically, relative_velocity_z. Actual divergence of a 2D flow field is 2 * relative_velocity_z.
int pprz_svd_float(float **a, float *w, float **v, int m, int n)
SVD decomposition.
#define MAT_SUB(_i, _j, C, A, B)
float surface_roughness
The error of the linear fit is a measure of surface roughness.
#define E
static float p[2][2]
#define A
#define MAT_MUL(_i, _k, _j, C, A, B)
#define MIN_SAMPLES_FIT