Thursday, December 26, 2013

SURF and BRISK in OpenCV

Today I spent half a day to compare the two feature detector and descriptor SURF and BRISK.

The conclusion is

  • SURF is more accurate but takes much longer
  • BRISK is 10 times faster with comparable accuracy. 
One thing to notice:
when using Brute-force matcher, mind the normType. See HERE
C++: BFMatcher::BFMatcher(int normType=NORM_L2, bool crossCheck=false )
* normType – One of NORM_L1, NORM_L2, NORM_HAMMING,NORM_HAMMING2. L1 and L2 norms are preferable choices for SIFT and SURF descriptors, 
* NORM_HAMMING should be used with ORB, BRISK and BRIEF, NORM_HAMMING2 should be used with ORB when WTA_K==3 or 4.

The other thing to notice is the threshold used to select good matches. Maybe RANSAC should be used.
SURF feature
BRISK feature

Saturday, November 16, 2013

Das ist eine brilliante Idee

Teil 4 Lektion 1

  • Gibt es schon Reaktionen auf den Sprachkurse?
  • Ja, wir haben viele Hörerbriefe bekommen.
  • Und was steht in den Briefen?
  • Die kann ich doch nicht alle vorlesen. Das dauert viel zu lange.
  • Nicht alle, aber einige - bitte!
  • Ja, das interessiert mich auch.
  • Na gut. 
  • Aber bitte, machen Sie es kurz!
  • Hier habe ich einen Brief von Herrn Card aus Amerika - Moment. Mir gefallen die Abenteuer von Andreas - asl Portier im Hotel Europa.
  • Mir auch
  • s
  • Und hier ist ein Brief von Angela aus Kolumbien. Sie schreibt, ' Ich bin so glücklick, weil ich die Grammatik studiert habe. Jetzt verstehe ich den Akkusativ - der war immer ...'
  • Grammatik, Grammatik, Akkusativ - das ist ja langweilig! Schreiben die Hörer den nichts über mich?! Wie finden die Hörer mich, das will ich wissen!
  • Kein Problem, Ex.  Hier ist ein Brief aus England - da steht etwas über dich! ' The intruduction of Ex is a brilliant idea. '
  • Das verstehe ich doch nicht! Was heißt das denn auf deutsch?
  • Du bist eine brillante Idee!
  • Idee? Wieso bin ich eine Idee? Ich bin ich!
  • s
  • Das hier ist noch ganz wichtig: Manche Hörer schreiben, daß sie Ex nicht so gut verstehen. 
  • Wir können ihr ja eine andere Stimme geben.
  • Probieren wir es doch mal! Ex, sprich mal etwas!
  • Ben dem Zauberwort sollte ich das Buch verlassen und ...
  • Okay! Stop! Und noch einmal bitte!
  • Ben dem Zauberwort sollte ich das Buch verlassen und ...
  • Kann ihre Stimme nicht ganz normal bleiben?
  • Nein - Ex ist ja eine besondere Person, ein weiblicher kobold, da braucht sie auch eine besondere Stimme. 
  • Das finde ich auch!
  • Aber das ist ein technisches Problem. Das lösen wir später. 

Saturday, November 2, 2013

Finds an object pose from 3D-2D point correspondences

C++: bool solvePnP(InputArray objectPoints, InputArray imagePoints, InputArray cameraMatrix, InputArray distCoeffs, OutputArray rvec, OutputArray tvec, bool useExtrinsicGuess=false, int flags=ITERATIVE )

  • objectPoints – Array of object points in the object coordinate space, 3xN/Nx3 1-channel or 1xN/Nx1 3-channel, where N is the number of points. vector<Point3f> can be also passed here.
  • imagePoints – Array of corresponding image points, 2xN/Nx2 1-channel or 1xN/Nx1 2-channel, where N is the number of points. vector<Point2f> can be also passed here.
  • cameraMatrix – Input camera matrix A = \vecthreethree{fx}{0}{cx}{0}{fy}{cy}{0}{0}{1} .
  • distCoeffs – Input vector of distortion coefficients (k_1, k_2, p_1, p_2[, k_3[, k_4, k_5, k_6]]) of 4, 5, or 8 elements. If the vector is NULL/empty, the zero distortion coefficients are assumed.
  • rvec – Output rotation vector (see Rodrigues() ) that, together with tvec , brings points from the model coordinate system to the camera coordinate system.
  • tvec – Output translation vector.
  • useExtrinsicGuess – If true (1), the function uses the provided rvec and tvec values as initial approximations of the rotation and translation vectors, respectively, and further optimizes them.
  • flags –
    Method for solving a PnP problem:
    • CV_ITERATIVE Iterative method is based on Levenberg-Marquardt optimization. In this case the function finds such a pose that minimizes reprojection error, that is the sum of squared distances between the observed projections imagePoints and the projected (using projectPoints() ) objectPoints .
    • CV_P3P Method is based on the paper of X.S. Gao, X.-R. Hou, J. Tang, H.-F. Chang “Complete Solution Classification for the Perspective-Three-Point Problem”. In this case the function requires exactly four object and image points.
    • CV_EPNP Method has been introduced by F.Moreno-Noguer, V.Lepetit and P.Fua in the paper “EPnP: Efficient Perspective-n-Point Camera Pose Estimation”.
The function estimates the object pose given a set of object points, their corresponding image projections, as well as the camera matrix and the distortion coefficients.

Friday, October 11, 2013

Tuesday, September 10, 2013

Thursday, August 15, 2013

Image stiching in OpenCV

Reference to

Main steps of the code:
  1. Load two images;
  2. Convert to gray scale;
  3. Using SURF detector to find SURF descriptor in both images;
  4. matching the SURF descriptor using FLANN Matcher;
  5. Postprocessing matches to find good matches
  6. Using RANSAC to estimate the Homography matrix using the matched SURF descriptors;
  7. Warping the images based on the homography matrix

The image shows the definition of Homography which transforms 2d Planar point to the image plane. 

 The following image is shows the initial image and the matched SURF features of the two. Homography is derived from the left image to the right image.

Stiching image result.

Source Code

Tuesday, August 6, 2013

Resource to learn optimization

Theory behind MPC

MPC is based on iterative, finite horizon optimization of a plant model. At time t the current plant state is sampled and a cost minimizing control strategy is computed (via a numerical minimization algorithm) for a relatively short time horizon in the future: [t,t+T]. Specifically, an online or on-the-fly calculation is used to explore state trajectories that emanate from the current state and find (via the solution of Euler-Lagrange equations) a cost-minimizing control strategy until time t+T. Only the first step of the control strategy is implemented, then the plant state is sampled again and the calculations are repeated starting from the current state, yielding a new control and new predicted state path. The prediction horizon keeps being shifted forward and for this reason MPC is also called receding horizon control.

Sunday, May 26, 2013

Hessian Matrix and Jacobian Matrix

Hessian Matrix
In mathematics, the Hessian matrix or Hessian is a square matrix of second-order partial derivatives of a function. It describes the local curvature of a function of many variables. The Hessian matrix was developed in the 19th century by the German mathematician Ludwig Otto Hesse and later named after him. Hesse originally used the term "functional determinants".
Given the real-valued function
f(x_1, x_2, \dots, x_n),\,\!
if all second partial derivatives of f exist and are continuous over the domain of the function, then the Hessian matrix of f is
H(f)_{ij}(\mathbf x) = D_i D_j f(\mathbf x)\,\!
where x = (x1x2, ..., xn) and Di is the differentiation operator with respect to the ith argument. Thus
H(f) = \begin{bmatrix}
\dfrac{\partial^2 f}{\partial x_1^2} & \dfrac{\partial^2 f}{\partial x_1\,\partial x_2} & \cdots & \dfrac{\partial^2 f}{\partial x_1\,\partial x_n} \\[2.2ex]
\dfrac{\partial^2 f}{\partial x_2\,\partial x_1} & \dfrac{\partial^2 f}{\partial x_2^2} & \cdots & \dfrac{\partial^2 f}{\partial x_2\,\partial x_n} \\[2.2ex]
\vdots & \vdots & \ddots & \vdots \\[2.2ex]
\dfrac{\partial^2 f}{\partial x_n\,\partial x_1} & \dfrac{\partial^2 f}{\partial x_n\,\partial x_2} & \cdots & \dfrac{\partial^2 f}{\partial x_n^2}
Because f is often clear from context, H(f)(\mathbf x) is frequently abbreviated to H(\mathbf x).
The Hessian matrix is related to the Jacobian matrix by H(f)(\mathbf x) = J(\nabla \! f)(\mathbf x).
The determinant of the above matrix is also sometimes referred to as the Hessian.[1]
Hessian matrices are used in large-scale optimization problems within Newton-type methods because they are the coefficient of the quadratic term of a local Taylor expansion of a function. That is,
y=f(\mathbf{x}+\Delta\mathbf{x})\approx f(\mathbf{x}) + J(\mathbf{x})\Delta \mathbf{x} +\frac{1}{2} \Delta\mathbf{x}^\mathrm{T} H(\mathbf{x}) \Delta\mathbf{x}
where J is the Jacobian matrix, which is a vector (the gradient) for scalar-valued functions. The full Hessian matrix can be difficult to compute in practice; in such situations, quasi-Newton algorithms have been developed that use approximations to the Hessian. The best-known quasi-Newton algorithm is the BFGS algorithm
Jacobian Matrix

In vector calculus, the Jacobian matrix is the matrix of all first-order partial derivatives of a vector-valued function. Specifically, suppose F : \mathbb{R}^n \rightarrow \mathbb{R}^m is a function (which takes as input real n-tuples and produces as output real m-tuples). Such a function is given by m real-valued component functions, F_1(x_1,\ldots,x_n),\ldots,F_m(x_1,\ldots,x_n). The partial derivatives of all these functions with respect to the variables x_1,\ldots,x_n (if they exist) can be organized in an m-by-n matrix, the Jacobian matrix J of F, as follows:
J=\begin{bmatrix} \dfrac{\partial F_1}{\partial x_1} & \cdots & \dfrac{\partial F_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \dfrac{\partial F_m}{\partial x_1} & \cdots & \dfrac{\partial F_m}{\partial x_n}  \end{bmatrix}.

Friday, May 3, 2013

Implementation of KalmanFilter in Opencv

Two Functions used: gemm , solve
gemm: Performs generalized matrix multiplication.
C++: void gemm(InputArray src1, InputArray src2, double alpha, InputArray src3, double gamma, OutputArray dst, intflags=0 )
The function performs generalized matrix multiplication similar to the gemm functions in BLAS level 3. For example, gemm(src1,src2, alpha, src3, beta, dst, GEMM_1_T + GEMM_3_T) corresponds to
\texttt{dst} =  \texttt{alpha} \cdot \texttt{src1} ^T  \cdot \texttt{src2} +  \texttt{beta} \cdot \texttt{src3} ^T

solve: Solves one or more linear systems or least-squares problems.
C++: bool solve(InputArray src1, InputArray src2, OutputArray dst, int flags=DECOMP_LU)
  • solution (matrix inversion) method.
    • DECOMP_LU Gaussian elimination with optimal pivot element chosen.
    • DECOMP_CHOLESKY Cholesky LL^T factorization; the matrix src1 must be symmetrical and positively defined.
    • DECOMP_EIG eigenvalue decomposition; the matrix src1 must be symmetrical.
    • DECOMP_SVD singular value decomposition (SVD) method; the system can be over-defined and/or the matrix src1 can be singular.
    • DECOMP_QR QR factorization; the system can be over-defined and/or the matrix src1 can be singular.
    • DECOMP_NORMAL while all the previous flags are mutually exclusive, this flag can be used together with any of the previous; it means that the normal equations\texttt{src1}^T\cdot\texttt{src1}\cdot\texttt{dst}=\texttt{src1}^T\texttt{src2} are solved instead of the original system \texttt{src1}\cdot\texttt{dst}=\texttt{src2} .
The function solve solves a linear system or least-squares problem (the latter is possible with SVD or QR methods, or by specifying the flag DECOMP_NORMAL ):
\texttt{dst} =  \arg \min _X \| \texttt{src1} \cdot \texttt{X} -  \texttt{src2} \|

OpenCV: Surf Matching in Video Sequence

OpenCV: SURF Feature matching

  1. Load two images
  2. do SURF feature extraction
  3. Using Flann matching to match the keypoints
  4. Identify good matches
  5. find the object in the scene image

Wednesday, May 1, 2013

OpenCV: SURF Feature extractor

Main steps:
  1. Read Two Images
  2. Resize to half of its size
  3. Detectect Keypoints  using SURF
  4. Calculate Feature Descriptor
  5. Matching Descriptor using Brute Force Matcher
Original Image
 SURF Keypoints
SURF Matches

Sunday, April 28, 2013

Surf Detector

OpenCV Sobel edge detector


Sobel X Gradient

Sobel Y Gradient

Sobel Gradient

OpenCV Tutorial 1: Mat - The basic Image Container

Color Space In OpenCV
There are, however, many other color systems each with their own advantages:
  • RGB is the most common as our eyes use something similar, our display systems also compose colors using these.
  • The HSV and HLS decompose colors into their hue, saturation and value/luminance components, which is a more natural way for us to describe colors. You might, for example, dismiss the last component, making your algorithm less sensible to the light conditions of the input image.
  • YCrCb is used by the popular JPEG image format.
  • CIE L*a*b* is a perceptually uniform color space, which comes handy if you need to measure the distance of a given color to another color.
Sample Code

Tuesday, April 16, 2013

Camera Calibration tool box from Caltech
This is a release of a Camera Calibration Toolbox for Matlab® with a complete documentation. This document may also be used as a tutorial on camera calibration since it includes general information about calibration, references and related links. 

Monday, April 1, 2013

Latex array stretch

Stretch the row size of a table
\renewcommand\arraystretch{MyValue}% (MyValue=1.0 is for standard spacing)

Saturday, March 23, 2013

German 1

wir haben uns kennengelernt

Da begrüßte er mich, wir hatten uns persönlich nie kennengelernt, aber er legte Wert darauf, mich am 4. November kennenzulernen.

He greeted me, we had never met in person, but it was important to him that he meet me on the 4th of November.











【励志语录】1.好多人做不好自己,是因为总想着做别人!2.从不奢求生活能给予我最好的,只是执着于寻求最适合我的!3.宁愿跑起来被拌倒无数次,也不愿规规矩矩走一辈子,就算跌倒也要豪迈的笑 。4.不要生气要争气,不要看破要突破,不要嫉妒要欣赏,不要托延要积极,不要心动要行动。

Saturday, March 2, 2013

Convex Optimization

Reference in Convex optimization


The following problems are all convex minimization problems, or can be transformed into convex minimizations problems via a change of variables:


Convex minimization problems can be solved by the following contemporary methods:[4]
Other methods of interest:
Subgradient methods can be implemented simply and so are widely used.[5] Dual subgradient methods are subgradient methods applied to a dual problem. The drift-plus-penalty method is similar to the dual subgradient method, but takes a time average of the primal variables.

Friday, March 1, 2013

OpenOF Framework for Sparse Non-linear Least Squares Optimization on a GPU

OpenOF Framework for Sparse Non-linear Least Squares Optimization on a GPU
With OpenOF, a framework is presented, which enables developers to design sparse optimizations regarding parameters and measurements and utilize the parallel power of a GPU

This code is written in Python with three major libraries: Thrust, CUSP and SymPy. Code framework is written in Python but can also generate C++ code.

Code website

Process of Nonlinear least squares optimization: 
1. Iterative method
2. Linearize the cost function in each iteration
3. Levenberg-Marquardt(LM) algorithm is standard, combing the Gauss-Newton algorithm with the gradient descent approach. LM guarantees convergence.
4. In each interation , solving linear Ax = b is most intensive.
5. Sparse matrix representation is used: sparseLM (Lourakis, 2010) and g2o (Kummerle et al., 2011), but on CPU
6. Solving Ax =b, many algorithms can achieve, Cholesky docomposition A = LDL'
7. this paper use Conjugate gradient (CG) approach on GPU.

Nonlinear least squares optimization is widely used in SLAM and BA. 
The authors' some comments about three BA libraries:
   1.The SBA library (Lourakis and Argyros,2009) takes advantage of the special structure of the Hessian matrix to apply the Schur complement for solving the linear system. Nevertheless it has several drawbacks. Integrating additional parameters which remain identical for all measurements (e.g. camera calibration) is not possible, as the structure would change such that the Schur complement could not be applied anymore.

2. sparseLM (Lourakis, 2010) is slow.
3. g2o: the Jacobian is evaluated by numerical differentiation which is time consuming and also degrades the convergence rate.
4. ISAM: (Kaess et al., 2011),which address only a subset of problems, have been presented previously for least squares optimization

Overall Comment: this paper is claims to present an open source framework for sparse nonlinear opitmization. The cost functions is described in high level scripting language. It can not be used without GPU yet. It seems for me g2o or iSam would be more useful on CPU. 

Thursday, February 28, 2013


class Compare_int
public :
Compare(int a,int b)
int max( )
return (x>y)?x:y;
int min( )
return (x<y)?x:y;}
private :
int x,y;

class Compare_float
public :
Compare(float a,float b)
float max( )
{return (x>y)?x:y;}
float min( )
{return (x<y)?x:y;}
private :
float x,y;

C++在发展的后期增加了模板(template )的功能,提供了解决这类问题的途径。可以声明一个通用的类模板,它可以有一个或多个虚拟的类型参数,如对以上两个类可以综合写出以下的类模板:
template <class numtype> //声明一个模板,虚拟类型名为numtype
class Compare //类模板名为Compare
public :
Compare(numtype a,numtype b)
numtype max( )
{return (x>y)?x:y;}
numtype min( )
{return (x<y)?x:y;}
private :
numtype x,y;

  1. 声明类模板时要增加一行
    template <class 类型参数名>
  2. 原有的类型名int换成虚拟类型参数名numtype。


Compare_int cmp1(4,7); // Compare_int是已声明的类
Compare cmp(4,7); // Compare是类模板名
Compare <int> cmp(4,7);

例9.14 声明一个类模板,利用它分别实现两个整数、浮点数和字符的比较,求出大数和小数。
#include <iostream>
using namespace std;
template <class numtype>
class Compare
public :
Compare(numtype a,numtype b)
numtype max( )
{return (x>y)?x:y;}
numtype min( )
{return (x<y)?x:y;}
private :
numtype x,y;
int main( )
Compare<int > cmp1(3,7);//定义对象cmp1,用于两个整数的比较
cout<<cmp1.max( )<<″ is the Maximum of two integer numbers.″<<endl;
cout<<cmp1.min( )<<″ is the Minimum of two integer numbers.″<<endl<<endl;
Compare<float > cmp2(45.78,93.6); //定义对象cmp2,用于两个浮点数的比较
cout<<cmp2.max( )<<″ is the Maximum of two float numbers.″<<endl;
cout<<cmp2.min( )<<″ is the Minimum of two float numbers.″<<endl<<endl;
Compare<char> cmp3(′a′,′A′); //定义对象cmp3,用于两个字符的比较
cout<<cmp3.max( )<<″ is the Maximum of two characters.″<<endl;
cout<<cmp3.min( )<<″ is the Minimum of two characters.″<<endl;
return 0;
7 is the Maximum of two integers.
3 is the Minimum of two integers.
93.6 is the Maximum of two float numbers.
45.78 is the Minimum of two float numbers.
a is the Maximum of two characters.
A is the Minimum of two characters.

还有一个问题要说明: 上面列出的类模板中的成员函数是在类模板内定义的。如果改为在类模板外定义,不能用一般定义类成员函数的形式:
numtype Compare::max( ) {…} //不能这样定义类模板中的成员函数
template <class numtype>
numtype Compare<numtype>::max( )
{{return (x>y)?x:y;}

  1. 先写出一个实际的类。由于其语义明确,含义清楚,一般不会出错。
  2. 将此类中准备改变的类型名(如int要改变为float或char)改用一个自己指定的虚拟类型名(如上例中的numtype)。
  3. 在类声明前面加入一行,格式为
    template <class 虚拟类型参数>,如
    template <class numtype> //注意本行末尾无分号
    class Compare
    {…}; //类体
  4. 用类模板定义对象时用以下形式:
    类模板名<实际类型名> 对象名;
    类模板名<实际类型名> 对象名(实参表列);

    Compare<int> cmp;
    Compare<int> cmp(3,7);
  5. 如果在类模板外定义成员函数,应写成类模板形式:
    template <class 虚拟类型参数>
    函数类型 类模板名<虚拟类型参数>::成员函数名(函数形参表列) {…}

  1. 类模板的类型参数可以有一个或多个,每个类型前面都必须加class,如
    template <class T1,class T2>
    class someclass
    someclass<int,double> obj;
  2. 和使用类一样,使用类模板时要注意其作用域,只能在其有效作用域内用它定义对象。
  3. 模板可以有层次,一个类模板可以作为基类,派生出派生模板类。