Skip to content

Commit

Permalink
Optical flow filter and derotation fix (#1855)
Browse files Browse the repository at this point in the history
Here is some improvements to the optical flow code.

-  Kalman filter: I made a generic kalman filter just for the use of vision. The problem of using the HHF filter is that it also effect the reaction of the controls (for some reason there is too much dependency on the process model but I can't change the noise factors). I think that it is wise to filter vision results separate from the state filters. Right now it is implemented for optical flow, but later it can be used for tracking as well. I hope you will agree
 -  derotation correction factor: there seems to be an error between the calculated derotation flow by the gyro and the calculated flow. Now it has been fitted together and for the bebop a correction term has to be added (credits to @mmoiozo for figuring this out).

This has been test flown with two bebops.
  • Loading branch information
knmcguire authored and flixr committed Sep 29, 2016
1 parent 10eb920 commit ca14142
Show file tree
Hide file tree
Showing 9 changed files with 488 additions and 73 deletions.
4 changes: 2 additions & 2 deletions conf/airframes/TUDELFT/tudelft_ardrone2_opticflow_stereo.xml
Expand Up @@ -23,8 +23,8 @@

<modules main_freq="512">
<module name="bat_voltage_ardrone2"/>
<module name="cv_opticflow"/>
<module name="opticflow_hover"/>
<!--module name="cv_opticflow"/>
<module name="opticflow_hover"/-->
<module name="stereocam">
<configure name="STEREO_UART" value="UART1"/>
<configure name="STEREO_BAUD" value="B1000000"/>
Expand Down
24 changes: 8 additions & 16 deletions conf/airframes/TUDELFT/tudelft_bebop_opticflow.xml
Expand Up @@ -35,13 +35,13 @@

<module name="cv_opticflow">
<define name="OPTICFLOW_CAMERA" value="bottom_camera"/>
<define name="MAX_HORIZON" value="2"/>
<!--define name="CAMERA_ROTATED_180" value="1"/-->
<!--define name="OPTICFLOW_DEROTATION" value="0"/-->
<define name="MAX_HORIZON" value="10"/>
<define name="OPTICFLOW_FX" value="347.22222222"/> <!-- 2.5 / (3.6 * 2.0) * 1000 -->
<define name="OPTICFLOW_FY" value="347.22222222"/> <!-- 2.5 / (3.6 * 2.0) * 1000 -->
<define name="OPTICFLOW_FOV_W" value="0.665499265"/> <!-- 2 * arctan(240 / (2*347.22222222)) -->
<define name="OPTICFLOW_FOV_H" value="0.665499265"/> <!-- 2 * arctan(240 / (2*347.22222222)) -->
<define name="OPTICFLOW_DEROTATION_CORRECTION_FACTOR_X" value="0.8"/> <!--Obtained from a linefit-->
<define name="OPTICFLOW_DEROTATION_CORRECTION_FACTOR_Y" value="0.85"/> <!--Obtained from a linefit-->
</module>

<module name="video_capture">
Expand Down Expand Up @@ -144,19 +144,11 @@
</section>

<section name="STABILIZATION_ATTITUDE_INDI" prefix="STABILIZATION_INDI_">
<!-- control effectiveness >
<define name="G1_P" value="0.0639"/>
<define name="G1_Q" value="0.0361"/>
<define name="G1_R" value="0.0022"/>
<define name="G2_R" value="0.1450"/>
<define name="G1_P" value="0.038703"/>
<define name="G1_Q" value="0.033092"/>
<define name="G1_R" value="0.000892"/>
<define name="G2_R" value="0.110938"/-->
<define name="G1_P" value="0.047332"/>
<define name="G1_Q" value="0.032168"/>
<define name="G1_R" value="0.002326"/>
<define name="G2_R" value="0.122961"/>
<!-- control effectiveness -->
<define name="G1_P" value="0.052669"/>
<define name="G1_Q" value="0.03326"/>
<define name="G1_R" value="0.002377"/>
<define name="G2_R" value="0.136936"/>

<!-- reference acceleration for attitude control -->
<define name="REF_ERR_P" value="600.0"/>
Expand Down
13 changes: 12 additions & 1 deletion conf/modules/cv_opticflow.xml
Expand Up @@ -10,7 +10,10 @@

<!-- Optical flow calculation parameters -->
<section name="OPTICFLOW" prefix="OPTICFLOW_">
<define name="AGL_ID" value="ABI_SENDER_ID" description="ABI sender id for AGL message (sonar measurement) (default: ABI_BROADCAST)"/>
<define name="AGL_ID" value="ABI_BROADCAST" description="ABI receive id for ABI message (sonar measurement) (default: ABI_BROADCAST)"/>
<define name="IMU_ID" value="ABI_BROADCAST" description="ABI receive id for ABI message (accelerometers) (default: ABI_BROADCAST)"/>
<define name="BODY_TO_IMU_ID" value="ABI_BROADCAST" description="ABI receive id for ABI message (body to imu quat) (default: ABI_BROADCAST)"/>
<define name="SEND_ABI_ID" value="1" description="ABI sender id for opticflow messages (default: 1)"/>

<!-- Video device parameters -->
<define name="CAMERA" value="bottom_camera|front_camera" description="The V4L2 camera device that is used for the calculations"/>
Expand All @@ -20,13 +23,18 @@
<define name="FOV_H" value="0.67020643276" description="The field of view height of the bottom camera (Defaults are from an ARDrone 2)"/>
<define name="FX" value="343.1211" description="Field in the x direction of the camera (Defaults are from an ARDrone 2)"/>
<define name="FY" value="348.5053" description="Field in the y direction of the camera (Defaults are from an ARDrone 2)"/>
<define name="DEROTATION_CORRECTION_FACTOR_X" value="1.0" description="Correction factor for derotation (in x direction), estimated from a fit between the gyro's rates and the resulting flow (caused by the camera not being exactly in the middle (Defaults are from an ARDrone 2)"/>
<define name="DEROTATION_CORRECTION_FACTOR_Y" value="1.0" description="Correction factor for derotation (in y direction), estimated from a fit between the gyro's rates and the resulting flow (caused by the camera not being exactly in the middle (Defaults are from an ARDrone 2)"/>

<!-- General optical flow calculation parameters -->
<define name="METHOD" value="0" description="Method used to calculate optical flow"/>
<define name="WINDOW_SIZE" value="10" description="Window size used for block matching (pixels)"/>
<define name="SEARCH_DISTANCE" value="10" description="Maximum search distance for blockmatching (pixels)"/>
<define name="SUBPIXEL_FACTOR" value="10" description="Amount of subpixels per pixel, used for more precise (subpixel) calculations of the flow"/>
<define name="DEROTATION" value="1" description="Derotation either turned on or off (depended on gyroscope measurements)"/>
<define name="MEDIAN_FILTER" value="0" description="A median filter on the resulting velocities to be turned on or off (last 5 measurements)"/>
<define name="KALMAN_FILTER" value="1" description="A kalman filter on the resulting velocities to be turned on or off (fused with accelerometers)"/>
<define name="KALMAN_FILTER_PROCESS_NOISE" value="0.01" description="The expected variance of the error of the model's prediction in the kalman filter"/>

<!-- Lucas Kanade optical flow calculation parameters -->
<define name="MAX_TRACK_CORNERS" value="25" description="The maximum amount of corners the Lucas Kanade algorithm is tracking between two frames"/>
Expand All @@ -52,6 +60,8 @@
<dl_setting var="opticflow.subpixel_factor" module="computer_vision/opticflow_module" min="0" step="10" max="1000" shortname="subpixel_factor" param="OPTICFLOW_SUBPIXEL_FACTOR"/>
<dl_setting var="opticflow.derotation" min="0" step="1" max="1" module="computer_vision/opticflow_module" values="OFF|ON" shortname="derotation" param="OPTICFLOW_DEROTATION"/>
<dl_setting var="opticflow.median_filter" module="computer_vision/opticflow_module" min="0" step="1" max="1" values="OFF|ON" shortname="median_filter" param="OPTICFLOW_MEDIAN_FILTER"/>
<dl_setting var="opticflow.kalman_filter" module="computer_vision/opticflow_module" min="0" step="1" max="1" values="OFF|ON" shortname="kalman_filter" param="OPTICFLOW_KALMAN_FILTER"/>
<dl_setting var="opticflow.kalman_filter_process_noise" module="computer_vision/opticflow_module" min="0.0001" step="0.0001" max="0.1" shortname="KF_process_noise" param="OPTICFLOW_KALMAN_FILTER_PROCESS_NOISE"/>

<!-- Specifically for Lucas Kanade and FAST9 -->
<dl_setting var="opticflow.max_track_corners" module="computer_vision/opticflow_module" min="0" step="1" max="500" shortname="max_trck_corners" param="OPTICFLOW_MAX_TRACK_CORNERS"/>
Expand Down Expand Up @@ -100,6 +110,7 @@
<file name="fast_rosten.c" dir="modules/computer_vision/lib/vision"/>
<file name="lucas_kanade.c" dir="modules/computer_vision/lib/vision"/>
<file name="edge_flow.c" dir="modules/computer_vision/lib/vision"/>
<file name="kalman_filter_vision.c" dir="modules/computer_vision/lib/filters"/>
</makefile>


Expand Down
205 changes: 205 additions & 0 deletions sw/airborne/modules/computer_vision/lib/filters/kalman_filter_vision.c
@@ -0,0 +1,205 @@
/*
* kalman_filter_vision.c
*
* Created on: Sep 12, 2016
* Author: knmcguire
*/
#include "kalman_filter_vision.h"
#include "math/pprz_algebra_float.h"

//TODO: implement kalman filter for larger states
//TODO: implement extended kalman filter



/**
* A simple linear 2D kalman filter, computed using floats and matrices. To be used for vision task like optical flow, tracking markers etc.
* @param[in] *model The process model for the prediction of the next state (array size of [1 x 4])
* @param[in] *measurements Values of the measurements (array size of [1 x 2])
* @param[in] *covariance Covariance matrix of previous iteration (array size of [1 x 4])
* @param[out] *covariance Updated Covariance matrix (array size of [1 x 4])
* @param[in] *state Previous state vector with estimates(array size of [1 x 2])
* @param[out] *state Updated state vector with estimates
* @param[in] *measurement_noise Expected variance of the noise of the measurements
* @param[in] *process_noise Expected variance of the noise of the process model
*/
void kalman_filter_linear_2D_float(float *model, float *measurements, float *covariance, float *state
, float *process_noise, float *measurement_noise)
{
//......................Preparation kalman filter ..................... //

// process model (linear)
float _G[2][2];
MAKE_MATRIX_PTR(G, _G, 2);
G[0][0] = model[0];
G[0][1] = model[1];
G[1][0] = model[2];
G[1][1] = model[3];

// transpose of G
float _Gtrans[2][2];
MAKE_MATRIX_PTR(Gtrans, _Gtrans, 2);
float_mat_copy(Gtrans, G, 2, 2);
float_mat_transpose(Gtrans, 2);

// Observation model (linear)
// note: right now both velocity and acceleration are observed
float _H[2][2];
MAKE_MATRIX_PTR(H, _H, 2);
H[0][0] = 1.0f;
H[0][1] = 0.0f;
H[1][0] = 0.0f;
H[1][1] = 1.0f;

//transpose of H
float _Htrans[2][2];
MAKE_MATRIX_PTR(Htrans, _Htrans, 2);
float_mat_copy(Htrans, H, 2, 2);
float_mat_transpose(Htrans, 2);

//Previous state
float _Xprev[2][1];
MAKE_MATRIX_PTR(Xprev, _Xprev, 2);
Xprev[0][0] = state[0];
Xprev[1][0] = state[1]; //state[1];

//Previous covariance
float _Pprevious[2][2];
MAKE_MATRIX_PTR(Pprevious, _Pprevious, 2);
Pprevious[0][0] = covariance[0];
Pprevious[0][1] = covariance[1];
Pprevious[1][0] = covariance[2];
Pprevious[1][1] = covariance[3];

//measurements;
float _Z[2][1];
MAKE_MATRIX_PTR(Z, _Z, 2);
Z[0][0] = measurements[0];
Z[1][0] = measurements[1];

//Process noise model
float _Q[2][2];
MAKE_MATRIX_PTR(Q, _Q, 2);
Q[0][0] = process_noise[0];
Q[0][1] = 0.0f;
Q[1][0] = 0.0f;
Q[1][1] = process_noise[1];

//measurement nosie model
float _R[2][2];
MAKE_MATRIX_PTR(R, _R, 2);
R[0][0] = measurement_noise[0];
R[0][1] = 0.0f;
R[1][0] = 0.0f;
R[1][1] = measurement_noise[1];

//Variables during kalman computation:
float _Xpredict[2][1];
MAKE_MATRIX_PTR(Xpredict, _Xpredict, 2);
float _Xnext[2][1];
MAKE_MATRIX_PTR(Xnext, _Xnext, 2);

float _Ppredict[2][2];
MAKE_MATRIX_PTR(Ppredict, _Ppredict, 2);
float _Pnext[2][2];
MAKE_MATRIX_PTR(Pnext, _Pnext, 2);

float _K[2][2];
MAKE_MATRIX_PTR(K, _K, 2);

float _eye[2][2];
MAKE_MATRIX_PTR(eye, _eye, 2);
eye[0][0] = 1.0f;
eye[0][1] = 0.0f;
eye[1][0] = 0.0f;
eye[1][1] = 1.0f;

float _temp_mat[2][2];
MAKE_MATRIX_PTR(temp_mat, _temp_mat, 2);
float _temp_mat2[2][2];
MAKE_MATRIX_PTR(temp_mat2, _temp_mat2, 2)
float _temp_mat3[2][2];
MAKE_MATRIX_PTR(temp_mat3, _temp_mat3, 2)

float _temp_vec[2][1];
MAKE_MATRIX_PTR(temp_vec, _temp_vec, 2);
float _temp_vec2[2][1];
MAKE_MATRIX_PTR(temp_vec2, _temp_vec2, 2);


//......................KALMAN FILTER ..................... //


// 1. calculate state predict

//Xpredict = G* Xprev;
float_mat_mul(Xpredict, G, Xprev, 2, 2, 1);
//......................KALMAN FILTER ..................... //

// 2. calculate covariance predict

// Ppredict = G*Pprevious*Gtrans + Q
//...Pprevious*Gtrans...
float_mat_mul(temp_mat, Pprevious, Gtrans, 2, 2, 2);
//G*Pprevious*Gtrans...
float_mat_mul(temp_mat2, G, temp_mat, 2, 2, 2);
//G*Pprevious*Gtrans+Q
float_mat_sum(Ppredict, temp_mat2, Q, 2, 2);

// 3. Calculate Kalman gain

// K = Ppredict * Htrans /( H * Ppredict * Htrans + R)
// ... Ppredict * Htrans ...
float_mat_mul(temp_mat, Ppredict, Htrans, 2, 2, 2);
//... H * Predict * Htrans
float_mat_mul(temp_mat2, H, temp_mat, 2, 2, 2);
//..( H * Ppredict * Htrans + R)
float_mat_sum(temp_mat3, temp_mat2, R, 2, 2);
//...inv( H * Ppredict * Htrans + R)
//TODO: Make a matrix inverse function for more than 2x2 matrix!

float det_temp2 = 1 / (temp_mat3[0][0] * temp_mat3[1][1] - temp_mat3[0][1] * temp_mat3[1][0]);
temp_mat2[0][0] = det_temp2 * (temp_mat3[1][1]);
temp_mat2[0][1] = det_temp2 * (-1 * temp_mat3[1][0]);
temp_mat2[1][0] = det_temp2 * (-1 * temp_mat3[0][1]);
temp_mat2[1][1] = det_temp2 * (temp_mat3[0][0]);
// K = Ppredict * Htrans / *inv( H * Ppredict * Htrans + R)
float_mat_mul(K, temp_mat, temp_mat2, 2, 2, 2);

// 4. Update state estimate

//Xnext = Xpredict + K *(Z - Htrans * Xpredict)
// ... Htrans * Xpredict)
float_mat_mul(temp_vec, Htrans, Xpredict, 2, 2, 1);

//... (Z - Htrans * Xpredict)
float_mat_diff(temp_vec2, Z, temp_vec, 2, 1);

// ... K *(Z - Htrans * Xpredict)
float_mat_mul(temp_vec, K, temp_vec2, 2, 2, 1);


//Xnext = Xpredict + K *(Z - Htrans * Xpredict)
float_mat_sum(Xnext, Xpredict, temp_vec, 2, 1);

// 5. Update covariance matrix

// Pnext = (eye(2) - K*H)*P_predict
// ...K*H...
float_mat_mul(temp_mat, K, H, 2, 2, 2);
//(eye(2) - K*H)
float_mat_diff(temp_mat2, eye, temp_mat, 2, 2);
// Pnext = (eye(2) - K*H)*P_predict
float_mat_mul(Pnext, temp_mat2, Ppredict, 2, 2, 2);


//save values for next state
covariance[0] = Pnext[0][0];
covariance[1] = Pnext[0][1];;
covariance[2] = Pnext[1][0];;
covariance[3] = Pnext[1][1];;


state[0] = Xnext[0][0];
state[1] = Xnext[1][0];
}
@@ -0,0 +1,15 @@
/*
* kalman_filter_vision.h
*
* Created on: Sep 12, 2016
* Author: knmcguire
*/

#ifndef KALMAN_FILTER_VISION_H_
#define KALMAN_FILTER_VISION_H_

void kalman_filter_linear_2D_float(float *model, float *measurements, float *covariance, float *state,
float *process_noise, float *measurement_noise);


#endif /* KALMAN_FILTER_VISION_H_ */
Expand Up @@ -30,6 +30,8 @@
#define _INTER_THREAD_DATA_H

#include "math/pprz_algebra_float.h"
#include "math/pprz_algebra_int.h"


/* The result calculated from the opticflow */
struct opticflow_result_t {
Expand Down Expand Up @@ -60,6 +62,8 @@ struct opticflow_result_t {
struct opticflow_state_t {
struct FloatRates rates; ///< Body rates
float agl; ///< height above ground [m]
struct FloatQuat imu_to_body_quat; ///< imu to body quaternion
struct Int32Vect3 accel_imu_meas; ///< imu acceleration in imu's coordinates
};

#endif
Expand Down

0 comments on commit ca14142

Please sign in to comment.