Paparazzi UAS v7.0_unstable
Paparazzi is a free software Unmanned Aircraft System.
Loading...
Searching...
No Matches
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#include "math/RANSAC.h"
46#include "std.h"
47
48// Is this still necessary?
49#define MAX_COUNT_PT 50
50
51#define MIN_SAMPLES_FIT 3
52
53#define N_PAR_TR_FIT 6
54
69{
70 // Are there enough flow vectors to perform a fit?
71 if (count < MIN_SAMPLES_FIT) {
72 // return that no fit was made:
73 return false;
74 }
75
76 // fit linear flow field:
79 &info->fit_error, &min_error_u, &min_error_v, &info->n_inliers_u, &info->n_inliers_v);
80
81 // extract information from the parameters:
83
84 // surface roughness is equal to fit error:
85 info->surface_roughness = info->fit_error;
86 info->divergence = info->relative_velocity_z;
87 float diverg = (fabsf(info->divergence) < 1E-5) ? 1E-5 : info->divergence;
88 info->time_to_contact = 1.0f / diverg;
89
90 // return successful fit:
91 return true;
92}
93
110 float *parameters_u, float *parameters_v, float *fit_error, float *min_error_u, float *min_error_v, int *n_inliers_u,
111 int *n_inliers_v)
112{
113
114 // We will solve systems of the form A x = b,
115 // where A = [nx3] matrix with entries [x, y, 1] for each optic flow location
116 // and b = [nx1] vector with either the horizontal (bu) or vertical (bv) flow.
117 // x in the system are the parameters for the horizontal (pu) or vertical (pv) flow field.
118
119 // local vars for iterating, random numbers:
120 int sam, p, i_rand, si, add_si;
121
122 // ensure that n_samples is high enough to ensure a result for a single fit:
124 // n_samples should not be higher than count:
125 n_samples = (n_samples < count) ? n_samples : count;
126
127 // initialize matrices and vectors for the full point set problem:
128 // this is used for determining inliers
129 float _AA[count][3];
130 MAKE_MATRIX_PTR(AA, _AA, count);
131 float _bu_all[count][1];
133 float _bv_all[count][1];
135 for (sam = 0; sam < count; sam++) {
136 AA[sam][0] = (float) vectors[sam].pos.x;
137 AA[sam][1] = (float) vectors[sam].pos.y;
138 AA[sam][2] = 1.0f;
139 bu_all[sam][0] = (float) vectors[sam].flow_x;
140 bv_all[sam][0] = (float) vectors[sam].flow_y;
141 }
142
143 // later used to determine the error of a set of parameters on the whole set:
144 float _bb[count][1];
145 MAKE_MATRIX_PTR(bb, _bb, count);
146 float _C[count][1];
147 MAKE_MATRIX_PTR(C, _C, count);
148
149 // ***************
150 // perform RANSAC:
151 // ***************
152
153 // set up variables for small linear system solved repeatedly inside RANSAC:
154 float _A[n_samples][3];
156 float _bu[n_samples][1];
158 float _bv[n_samples][1];
160 float w[n_samples], _v[3][3];
161 MAKE_MATRIX_PTR(v, _v, 3);
162 float _pu[3][1];
164 float _pv[3][1];
166
167 // iterate and store parameters, errors, inliers per fit:
168 float PU[n_iterations * 3];
169 float PV[n_iterations * 3];
170 float errors_pu[n_iterations];
171 errors_pu[0] = 0.0;
172 float errors_pv[n_iterations];
173 errors_pv[0] = 0.0;
176 int it, ii;
177 for (it = 0; it < n_iterations; it++) {
178 // select a random sample of n_sample points:
180 i_rand = 0;
181
182 // sampling without replacement:
183 while (i_rand < n_samples) {
184 si = rand() % count;
185 add_si = 1;
186 for (ii = 0; ii < i_rand; ii++) {
187 if (sample_indices[ii] == si) { add_si = 0; }
188 }
189 if (add_si) {
191 i_rand ++;
192 }
193 }
194
195 // Setup the system:
196 for (sam = 0; sam < n_samples; sam++) {
197 A[sam][0] = (float) vectors[sample_indices[sam]].pos.x;
198 A[sam][1] = (float) vectors[sample_indices[sam]].pos.y;
199 A[sam][2] = 1.0f;
200 bu[sam][0] = (float) vectors[sample_indices[sam]].flow_x;
201 bv[sam][0] = (float) vectors[sample_indices[sam]].flow_y;
202 //printf("%d,%d,%d,%d,%d\n",A[sam][0],A[sam][1],A[sam][2],bu[sam][0],bv[sam][0]);
203 }
204
205 // Solve the small system:
206
207 // for horizontal flow:
208 // decompose A in u, w, v with singular value decomposition A = u * w * vT.
209 // u replaces A as output:
210 pprz_svd_float(A, w, v, n_samples, 3);
211 pprz_svd_solve_float(pu, A, w, v, bu, n_samples, 3, 1);
212 PU[it * 3] = pu[0][0];
213 PU[it * 3 + 1] = pu[1][0];
214 PU[it * 3 + 2] = pu[2][0];
215
216 // for vertical flow:
217 pprz_svd_solve_float(pv, A, w, v, bv, n_samples, 3, 1);
218 PV[it * 3] = pv[0][0];
219 PV[it * 3 + 1] = pv[1][0];
220 PV[it * 3 + 2] = pv[2][0];
221
222 // count inliers and determine their error on all points:
223 errors_pu[it] = 0;
224 errors_pv[it] = 0;
225 n_inliers_pu[it] = 0;
226 n_inliers_pv[it] = 0;
227
228 // for horizontal flow:
229 // bb = AA * pu:
230 MAT_MUL(count, 3, 1, bb, AA, pu);
231 // subtract bu_all: C = 0 in case of perfect fit:
232 MAT_SUB(count, 1, C, bb, bu_all);
233
234 for (p = 0; p < count; p++) {
235 C[p][0] = fabsf(C[p][0]);
236 if (C[p][0] < error_threshold) {
237 errors_pu[it] += C[p][0];
238 n_inliers_pu[it]++;
239 } else {
241 }
242 }
243
244 // for vertical flow:
245 // bb = AA * pv:
246 MAT_MUL(count, 3, 1, bb, AA, pv);
247 // subtract bv_all: C = 0 in case of perfect fit:
248 MAT_SUB(count, 1, C, bb, bv_all);
249
250 for (p = 0; p < count; p++) {
251 C[p][0] = fabsf(C[p][0]);
252 if (C[p][0] < error_threshold) {
253 errors_pv[it] += C[p][0];
254 n_inliers_pv[it]++;
255 } else {
257 }
258 }
259 }
260
261 // After all iterations:
262 // select the parameters with lowest error:
263 // for horizontal flow:
264 int param;
265 int min_ind = 0;
267 for (it = 1; it < n_iterations; it++) {
268 if (errors_pu[it] < *min_error_u) {
270 min_ind = it;
271 }
272 }
273 for (param = 0; param < 3; param++) {
274 parameters_u[param] = PU[min_ind * 3 + param];
275 }
276 *n_inliers_u = n_inliers_pu[min_ind];
277
278 // for vertical flow:
279 min_ind = 0;
281 for (it = 0; it < n_iterations; it++) {
282 if (errors_pv[it] < *min_error_v) {
284 min_ind = it;
285 }
286 }
287 for (param = 0; param < 3; param++) {
288 parameters_v[param] = PV[min_ind * 3 + param];
289 }
290 *n_inliers_v = n_inliers_pv[min_ind];
291
292 // error has to be determined on the entire set without threshold:
293 // bb = AA * pu:
294 MAT_MUL(count, 3, 1, bb, AA, pu);
295 // subtract bu_all: C = 0 in case of perfect fit:
296 MAT_SUB(count, 1, C, bb, bu_all);
297 *min_error_u = 0;
298 for (p = 0; p < count; p++) {
299 *min_error_u += fabsf(C[p][0]);
300 }
301 // bb = AA * pv:
302 MAT_MUL(count, 3, 1, bb, AA, pv);
303 // subtract bv_all: C = 0 in case of perfect fit:
304 MAT_SUB(count, 1, C, bb, bv_all);
305 *min_error_v = 0;
306 for (p = 0; p < count; p++) {
307 *min_error_v += fabsf(C[p][0]);
308 }
309 *fit_error = (*min_error_u + *min_error_v) / (2 * count);
310
311}
323{
324 // This method assumes a linear flow field in x- and y- direction according to the formulas:
325 // u = parameters_u[0] * x + parameters_u[1] * y + parameters_u[2]
326 // v = parameters_v[0] * x + parameters_v[1] * y + parameters_v[2]
327 // where u is the horizontal flow at image coordinate (x,y)
328 // and v is the vertical flow at image coordinate (x,y)
329
330 // relative velocities:
331 info->relative_velocity_z = (parameters_u[0] + parameters_v[1]) / 2.0f; // divergence / 2
332
333 // translation orthogonal to the camera axis:
334 // flow in the center of the image:
335 info->relative_velocity_x = -(parameters_u[2] + (im_width / 2.0f) * parameters_u[0] +
336 (im_height / 2.0f) * parameters_u[1]);
337 info->relative_velocity_y = -(parameters_v[2] + (im_width / 2.0f) * parameters_v[0] +
338 (im_height / 2.0f) * parameters_v[1]);
339
340 float arv_x = fabsf(info->relative_velocity_x);
341 float arv_y = fabsf(info->relative_velocity_y);
342
343 // extract inclination from flow field:
344 float threshold_slope = 1.0;
345 float eta = 0.002;
346
348 // there is no forward motion and not enough vertical motion, but enough horizontal motion:
349 info->slope_x = parameters_u[0] / info->relative_velocity_x;
350 } else if (arv_y >= 2 * threshold_slope) {
351 // there is sufficient vertical motion:
352 info->slope_x = parameters_v[0] / info->relative_velocity_y;
353 } else {
354 // there may be forward motion, then we can do a quadratic fit:
355 // a linear fit provides no information though
356 info->slope_x = 0.0f;
357 }
358
360 // there is no forward motion, little horizontal movement, but sufficient vertical motion:
361 info->slope_y = parameters_v[1] / info->relative_velocity_y;
362 } else if (arv_x >= 2 * threshold_slope) {
363 // there is sufficient horizontal motion:
364 info->slope_y = parameters_u[1] / info->relative_velocity_x;
365 } else {
366 // there could be forward motion, then we can do a quadratic fit:
367 // a linear fit provides no information though
368 info->slope_y = 0.0f;
369 }
370
371 // Focus of Expansion:
372 // the flow planes intersect the flow=0 plane in a line
373 // the FoE is the point where these 2 lines intersect (flow = (0,0))
374 // x:
376 if (fabsf(denominator) > 1E-5) {
377 info->focus_of_expansion_x = ((parameters_u[2] * parameters_v[1] - parameters_v[2] * parameters_u[1]) / denominator);
378 } else { info->focus_of_expansion_x = 0.0f; }
379 // y:
381 if (fabsf(denominator) > 1E-5) {
382 info->focus_of_expansion_y = (-(parameters_u[0] * (info->focus_of_expansion_x) + parameters_u[2]) / denominator);
383 } else { info->focus_of_expansion_y = 0.0f; }
384}
385
399bool analyze_flow_field(struct flow_t *vectors __attribute__((unused)), int count, float error_threshold __attribute__((unused)), int n_iterations __attribute__((unused)), int n_samples __attribute__((unused)),
400 int im_width __attribute__((unused)), int im_height __attribute__((unused)), float focal_length __attribute__((unused)), struct linear_flow_fit_info *info __attribute__((unused)))
401{
402 // Are there enough flow vectors to perform a fit?
403 if (count < N_PAR_TR_FIT) {
404 // return that no fit was made:
405 return false;
406 }
407 return true;
408
409 /*
410 // fit linear flow field:
411 //float parameters_u[N_PAR_TR_FIT], parameters_v[N_PAR_TR_FIT], min_error_u, min_error_v;
412 float parameters_comb[N_PAR_TR_FIT], min_error_comb;
413
414 // normalize the flow vectors. Is supposed to be better for the fit and the parameters will be easier to interpret:
415 struct flow_t *normalized_vectors = (struct flow_t *) malloc(count *sizeof(struct flow_t));
416
417
418 //float targets_x[count];
419 //float targets_y[count];
420
421 // separate horizontal / vertical fit:
422 //float D_x[count][N_PAR_TR_FIT-1]; // bias is added in the RANSAC routine
423 //float D_y[count][N_PAR_TR_FIT-1];
424
425 // combined fit:
426 int total_count = 2*count;
427 float targets[total_count];
428 float D_comb[total_count][N_PAR_TR_FIT];
429 int index;
430 for (int i = 0; i < count; i++) {
431
432 // normalize vectors:
433 normalized_vectors[i].pos.x = (vectors[i].pos.x - (im_width / 2.0f)) / focal_length;
434 normalized_vectors[i].pos.y = (vectors[i].pos.y - (im_height / 2.0f)) / focal_length;
435 normalized_vectors[i].flow_x = vectors[i].flow_x / focal_length;
436 normalized_vectors[i].flow_y = vectors[i].flow_y / focal_length;
437
438
439 // Combine vertical and horizontal flow in one system:
440 // As in "Estimation of Visual Motion Parameters Used for Bio-inspired Navigation", by Alkowatly et al.
441 index = i*2;
442 targets[index] = normalized_vectors[i].flow_x;
443 D_comb[index][0] = 1.0f;
444 D_comb[index][1] = normalized_vectors[i].pos.x;
445 D_comb[index][2] = normalized_vectors[i].pos.y;
446 D_comb[index][3] = normalized_vectors[i].pos.x * normalized_vectors[i].pos.y;
447 D_comb[index][4] = normalized_vectors[i].pos.x * normalized_vectors[i].pos.x;
448 D_comb[index][5] = 0.0f;
449
450 index++;
451 targets[index] = normalized_vectors[i].flow_y;
452 D_comb[index][0] = 0.0f;
453 D_comb[index][1] = normalized_vectors[i].pos.y;
454 D_comb[index][1] = -normalized_vectors[i].pos.x;
455 D_comb[index][2] = normalized_vectors[i].pos.y * normalized_vectors[i].pos.y;
456 D_comb[index][3] = normalized_vectors[i].pos.x * normalized_vectors[i].pos.y;
457 D_comb[index][4] = 1.0f;
458
459 }
460
461 // Perform RANSAC on the horizontal flow:
462 // RANSAC_linear_model(n_samples, n_iterations, error_threshold, targets_x, N_PAR_TR_FIT, D_x, count, parameters_u, &min_error_u);
463
464 // Perform RANSAC on the combined system:
465 bool use_bias = false;
466 RANSAC_linear_model(n_samples, n_iterations, error_threshold, targets, N_PAR_TR_FIT+1, D_comb, total_count, parameters_comb, &min_error_comb, use_bias);
467
468 // extract information from the parameters:
469 info->rotation_X = parameters_comb[3];
470 info->rotation_Y = -parameters_comb[4];
471 info->rotation_Z = parameters_comb[2];
472 info->divergence = parameters_comb[1];
473 info->relative_velocity_x = -parameters_comb[0] - info->rotation_Y;
474 info->relative_velocity_y = -parameters_comb[5] + info->rotation_X;
475 if(fabs(parameters_comb[1]) > 1E-3) {
476 info->focus_of_expansion_x = info->relative_velocity_x / parameters_comb[1];
477 info->focus_of_expansion_y = info->relative_velocity_y / parameters_comb[1];
478 }
479 else {
480 // FoE in the center of the image:
481 info->focus_of_expansion_x = 0.0f;
482 info->focus_of_expansion_y = 0.0f;
483 }
484
485 // free up allocated memory:
486 free(normalized_vectors);
487
488 // return successful fit:
489 return true;
490 */
491}
Perform Random Sample Consensus (RANSAC), a robust fitting method.
int n_samples
Definition detect_gate.c:85
#define MAKE_MATRIX_PTR(_ptr, _mat, _rows)
Make a pointer to a matrix of _rows lines.
#define E
#define A
Definition image.h:78
static float p[2][2]
#define N_PAR_TR_FIT
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,...
bool analyze_flow_field(struct flow_t *vectors, int count, float error_threshold, int n_iterations, int n_samples, int im_width, int im_height, float focal_length, struct linear_flow_fit_info *info)
Analyze a flow field, retrieving information on relative velocities and rotations,...
#define MIN_SAMPLES_FIT
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,...
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.
uint16_t foo
Definition main_demo5.c:58
Paparazzi floating point algebra.
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.
int pprz_svd_float(float **a, float *w, float **v, int m, int n)
SVD decomposition.
Matrix decompositions in floating point.
Simple matrix helper macros.
#define MAT_MUL(_i, _k, _j, C, A, B)
#define MAT_SUB(_i, _j, C, A, B)