Category Archives: Matlab

RANSAC Algorithm parameter explained

In this tutorial I explain the RANSAC algorithm, their corresponding parameters and how to choose the number of samples:

N = number of samples
e = probability that a point is an outlier
s = number of points in a sample
p = desired probability that we get a good sample
N =log(1-p) /log(1- (1- e) s )

ref: 1

Finding Homography Matrix using Singular-value Decomposition and RANSAC in OpenCV and Matlab

Leave a reply

\(\)

Solving a Homography problem leads to solving a set of homogeneous linear equations such below:

\begin{equation}
\left(
\begin{array}{ccccccccc}
-x1 & -y1 & -1 & 0 & 0 & 0 & x1*xp1 & y1*xp1 & xp1\\
0 & 0 & 0 & -x1 & -y1 & -1 & x1*yp1 & y1*yp1 & yp1\\
-x2 & -y2 & -1 & 0 & 0 & 0 & x2*xp2 & y2*xp2 & xp2\\
0 & 0 & 0 & -x2 & -y2 & -1 & x2*yp2 & y2*yp2 & yp2\\
-x3 & -y3 & -1 & 0 & 0 & 0 & x3*xp3 & y3*xp3 & xp3\\
0 & 0 & 0 & -x3 & -y3 & -1 & x3*yp3 & y3*yp3 & yp3\\
-x4 & -y4 & -1 & 0 & 0 & 0 & x4*xp4 & y4*xp4 & xp4\\
0 & 0 & 0 & -x4 & -y4 & -1 & x4*yp4 & y4*yp4 & yp4\\
\end{array}
\right) *H=0
\end{equation}

\( H^{*} \underset{H}{\mathrm{argmin}}= \|AH\|^{2} \) subject to \( \|H\|=1 \)

We can’t use least square since it’s a homogeneous linear equations (the other side of equation is 0 therfore we can’t just multyly it by the psudo inverse). To solve this problem we use Singular-value Decomposition (SVD).

For any given matrix \( A_{M{\times}N} \)

\begin{equation}
\underbrace{\mathbf{A}}_{M \times N} = \underbrace{\mathbf{U}}_{M \times M} \times \underbrace{\mathbf{\Sigma}}_{M\times N} \times \underbrace{\mathbf{V}^{\text{T}}}_{N \times N}
\end{equation}

\(U\) is an \( {M\times M} \) matrix with orthogonal matrix (columns are eigen vectors of A).
\(\Sigma\) is an \( {M\times N} \) matrix with non-negative entries, termed the singular values (diagonal entries are eigen values of A).
\(V \) is an \( {N\times N} \) orthogonal matrix.

\( H^{*} \) is the last column of \(V \)

Finding homography matrix in Matlab between 4 pairs of points

function H=homography(x1,y1,x2,y2,x3,y3,x4,y4 , xp1,yp1,xp2,yp2,xp3,yp3,xp4,yp4)
%This function will find the homography betweeb 4 points using svd


A=[
-x1  -y1  -1   0    0    0   x1*xp1   y1*xp1   xp1;
 0    0    0 -x1   -y1  -1   x1*yp1   y1*yp1   yp1;
-x2  -y2  -1   0    0    0   x2*xp2   y2*xp2   xp2;
 0    0    0 -x2   -y2  -1   x2*yp2   y2*yp2   yp2;
-x3  -y3  -1   0    0    0   x3*xp3   y3*xp3   xp3;
 0    0    0 -x3   -y3  -1   x3*yp3   y3*yp3   yp3;
-x4  -y4   -1  0    0    0   x4*xp4   y4*xp4   xp4;
 0    0    0  -x4  -y4  -1   x4*yp4   y4*yp4   yp4];

[U,S,V] = svd(A);


H=V(:,end);
H=reshape(H,3,3);

function H=homography(x1,y1,x2,y2,x3,y3,x4,y4 , xp1,yp1,xp2,yp2,xp3,yp3,xp4,yp4)

%This function will find the homography betweeb 4 points using svd

A=[

-x1 -y1 -1 0 0 0 x1*xp1 y1*xp1 xp1;

0 0 0 -x1 -y1 -1 x1*yp1 y1*yp1 yp1;

-x2 -y2 -1 0 0 0 x2*xp2 y2*xp2 xp2;

0 0 0 -x2 -y2 -1 x2*yp2 y2*yp2 yp2;

-x3 -y3 -1 0 0 0 x3*xp3 y3*xp3 xp3;

0 0 0 -x3 -y3 -1 x3*yp3 y3*yp3 yp3;

-x4 -y4 -1 0 0 0 x4*xp4 y4*xp4 xp4;

0 0 0 -x4 -y4 -1 x4*yp4 y4*yp4 yp4];

[U,S,V] = svd(A);

H=V(:,end);

H=reshape(H,3,3);

%This script will find the homography between 4 points and their repsective
%projection 

%These points are for vertex of a square in world cordinate system, the
%lengh of square is 150
x0=0;
y0=0;

x1=150;
y1=0;

x2=150;
y2=150;

x3=0;
y3=150;

%The projected position of points
xp0=100;
yp0=100;

xp1=200;
yp1=80;

xp2=220;
yp2=80;

xp3=100;
yp3=200; 

fprintf('Estimated homography matrix is:')
H=homography(x0,y0, x1,y1, x2,y2, x3,y3, xp0,yp0, xp1,yp1, xp2,yp2,xp3,yp3)

p0=[x0,y0; x1,y1 ; x2,y2 ; x3,y3];
p1=[xp0,yp0; xp1,yp1; xp2,yp2; xp3,yp3];

%We project the point (x2 y2) by H and we should get (xp2 yp2)
projected_point=[x2 y2 1] * H;
projected_point(1,1)/projected_point(1,3)
projected_point(1,2)/projected_point(1,3)

%This script will find the homography between 4 points and their repsective

%projection

%These points are for vertex of a square in world cordinate system, the

%lengh of square is 150

x0=0;

y0=0;

x1=150;

y1=0;

x2=150;

y2=150;

x3=0;

y3=150;

%The projected position of points

xp0=100;

yp0=100;

xp1=200;

yp1=80;

xp2=220;

yp2=80;

xp3=100;

yp3=200;

fprintf('Estimated homography matrix is:')

H=homography(x0,y0, x1,y1, x2,y2, x3,y3, xp0,yp0, xp1,yp1, xp2,yp2,xp3,yp3)

p0=[x0,y0; x1,y1 ; x2,y2 ; x3,y3];

p1=[xp0,yp0; xp1,yp1; xp2,yp2; xp3,yp3];

%We project the point (x2 y2) by H and we should get (xp2 yp2)

projected_point=[x2 y2 1] * H;

projected_point(1,1)/projected_point(1,3)

projected_point(1,2)/projected_point(1,3)

Finding homography matrix in OpenCV between 4 pairs of points

void get_homographyMatrix()
{
    std::vector<cv::Point2f> obj;
    std::vector<cv::Point2f> scene;
    std::vector<cv::Point2f> obj_projection;

    cv::Point2f A,B,C,D,A_P,B_P,C_P,D_P ;
    A.x=0;
    A.y=0;

    B.x=150;
    B.y=0;

    C.x=150;
    C.y=150;

    D.x=0;
    D.y=150;

    obj.push_back(A);
    obj.push_back(B);
    obj.push_back(C);
    obj.push_back(D);


    A_P.x=100;
    A_P.y=100;

    B_P.x=200;
    B_P.y=80;

    C_P.x=220;
    C_P.y=80;

    D_P.x=100;
    D_P.y=200;

    scene.push_back(A_P);
    scene.push_back(B_P);
    scene.push_back(C_P);
    scene.push_back(D_P);

    cv::Mat homography_matrix= cv::getPerspectiveTransform(obj,scene);
    std::cout<<"Estimated Homography Matrix is:" <<std::endl;
    std::cout<< homography_matrix <<std::endl;

    perspectiveTransform( obj, obj_projection, homography_matrix);
    for(std::size_t i=0;i<obj_projection.size();i++)
    {
        std::cout<<obj_projection.at(i).x <<"," <<obj_projection.at(i).y<<std::endl;
    }
}

void get_homographyMatrix()

{

std::vector<cv::Point2f> obj;

std::vector<cv::Point2f> scene;

std::vector<cv::Point2f> obj_projection;

cv::Point2f A,B,C,D,A_P,B_P,C_P,D_P ;

A.x=0;

A.y=0;

B.x=150;

B.y=0;

C.x=150;

C.y=150;

D.x=0;

D.y=150;

obj.push_back(A);

obj.push_back(B);

obj.push_back(C);

obj.push_back(D);

A_P.x=100;

A_P.y=100;

B_P.x=200;

B_P.y=80;

C_P.x=220;

C_P.y=80;

D_P.x=100;

D_P.y=200;

scene.push_back(A_P);

scene.push_back(B_P);

scene.push_back(C_P);

scene.push_back(D_P);

cv::Mat homography_matrix= cv::getPerspectiveTransform(obj,scene);

std::cout<<"Estimated Homography Matrix is:" <<std::endl;

std::cout<< homography_matrix <<std::endl;

perspectiveTransform( obj, obj_projection, homography_matrix);

for(std::size_t i=0;i<obj_projection.size();i++)

{

std::cout<<obj_projection.at(i).x <<"," <<obj_projection.at(i).y<<std::endl;

}

Finding homography matrix in OpenCV between N pairs of points

void findHomography_Test(int argc, char** argv)
{
//	The input arrays src_points and dst_points can be either
//		1)N-by-2(pixel coordinates) matrices
//		2)N-by-3 (homogeneous coordinates)	matrices
//	The final argument, homography, is just a 	3-by-3 matrix to


//	void cvFindHomography(const CvMat* src_points, 	const CvMat* dst_points, 	CvMat* 	homography 	);
    if( argc != 3 )
    {
        printf(" Usage: ./main_Homography <img1> <img2>\n");
        return;
    }

    Mat img_object = imread( argv[1], CV_LOAD_IMAGE_GRAYSCALE );
    Mat img_scene = imread( argv[2], CV_LOAD_IMAGE_GRAYSCALE );

    if( !img_object.data || !img_scene.data )
    {
        printf(" --(!) Error reading images \n"); return ;
    }

    //-- Step 1: Detect the keypoints using SURF Detector
    int minHessian = 400;

    SurfFeatureDetector detector( minHessian );

    std::vector<KeyPoint> keypoints_object, keypoints_scene;

    detector.detect( img_object, keypoints_object );
    detector.detect( img_scene, keypoints_scene );

    //-- Step 2: Calculate descriptors (feature vectors)
    SurfDescriptorExtractor extractor;

    Mat descriptors_object, descriptors_scene;

    extractor.compute( img_object, keypoints_object, descriptors_object );
    extractor.compute( img_scene, keypoints_scene, descriptors_scene );

    //-- Step 3: Matching descriptor vectors using FLANN matcher
    FlannBasedMatcher matcher;
    std::vector< DMatch > matches;
    matcher.match( descriptors_object, descriptors_scene, matches );

    double max_dist = 0; double min_dist = 100;

    //-- Quick calculation of max and min distances between keypoints
    for( int i = 0; i < descriptors_object.rows; i++ )
    {   double dist = matches[i].distance;
        if( dist < min_dist ) min_dist = dist;
        if( dist > max_dist ) max_dist = dist;
    }

    printf("-- Max dist : %f \n", max_dist );
    printf("-- Min dist : %f \n", min_dist );

    //-- Draw only "good" matches (i.e. whose distance is less than 3*min_dist )
    std::vector< DMatch > good_matches;

    for( int i = 0; i < descriptors_object.rows; i++ )
    { if( matches[i].distance < 1.5*min_dist )
        {
            good_matches.push_back( matches[i]);
        }
    }

    Mat img_matches;
    drawMatches( img_object, keypoints_object, img_scene, keypoints_scene,
    good_matches, img_matches, Scalar::all(-1), Scalar::all(-1),
    vector<char>(), DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS );


    //-- Localize the object from img_1 in img_2
    std::vector<Point2f> obj;
    std::vector<Point2f> scene;

    for( size_t i = 0; i < good_matches.size(); i++ )
    {
    //-- Get the keypoints from the good matches
        obj.push_back( keypoints_object[ good_matches[i].queryIdx ].pt );
        scene.push_back( keypoints_scene[ good_matches[i].trainIdx ].pt );
    }

    Mat H = findHomography( obj, scene, CV_RANSAC );

    //-- Get the corners from the image_1 ( the object to be "detected" )
    std::vector<Point2f> obj_corners(4);
    obj_corners[0] = Point(0,0); obj_corners[1] = Point( img_object.cols, 0 );
    obj_corners[2] = Point( img_object.cols, img_object.rows ); obj_corners[3] = Point( 0, img_object.rows );
    std::vector<Point2f> scene_corners(4);

    perspectiveTransform( obj_corners, scene_corners, H);


    //-- Draw lines between the corners (the mapped object in the scene - image_2 )
    Point2f offset( (float)img_object.cols, 0);
    line( img_matches, scene_corners[0] + offset, scene_corners[1] + offset, Scalar(0, 255, 0), 4 );
    line( img_matches, scene_corners[1] + offset, scene_corners[2] + offset, Scalar( 0, 255, 0), 4 );
    line( img_matches, scene_corners[2] + offset, scene_corners[3] + offset, Scalar( 0, 255, 0), 4 );
    line( img_matches, scene_corners[3] + offset, scene_corners[0] + offset, Scalar( 0, 255, 0), 4 );

    //-- Show detected matches
    imshow( "Good Matches & Object detection", img_matches );

    waitKey(0);

    return ;
}

100

101

102

103

104

105

106

107

108

109

110

void findHomography_Test(int argc, char** argv)

{

// The input arrays src_points and dst_points can be either

// 1)N-by-2(pixel coordinates) matrices

// 2)N-by-3 (homogeneous coordinates) matrices

// The final argument, homography, is just a 3-by-3 matrix to

// void cvFindHomography(const CvMat* src_points, const CvMat* dst_points, CvMat* homography );

if( argc != 3 )

{

printf(" Usage: ./main_Homography <img1> <img2>\n");

return;

}

Mat img_object = imread( argv[1], CV_LOAD_IMAGE_GRAYSCALE );

Mat img_scene = imread( argv[2], CV_LOAD_IMAGE_GRAYSCALE );

if( !img_object.data || !img_scene.data )

{

printf(" --(!) Error reading images \n"); return ;

}

//-- Step 1: Detect the keypoints using SURF Detector

int minHessian = 400;

SurfFeatureDetector detector( minHessian );

std::vector<KeyPoint> keypoints_object, keypoints_scene;

detector.detect( img_object, keypoints_object );

detector.detect( img_scene, keypoints_scene );

//-- Step 2: Calculate descriptors (feature vectors)

SurfDescriptorExtractor extractor;

Mat descriptors_object, descriptors_scene;

extractor.compute( img_object, keypoints_object, descriptors_object );

extractor.compute( img_scene, keypoints_scene, descriptors_scene );

//-- Step 3: Matching descriptor vectors using FLANN matcher

FlannBasedMatcher matcher;

std::vector< DMatch > matches;

matcher.match( descriptors_object, descriptors_scene, matches );

double max_dist = 0; double min_dist = 100;

//-- Quick calculation of max and min distances between keypoints

for( int i = 0; i < descriptors_object.rows; i++ )

{ double dist = matches[i].distance;

if( dist < min_dist ) min_dist = dist;

if( dist > max_dist ) max_dist = dist;

}

printf("-- Max dist : %f \n", max_dist );

printf("-- Min dist : %f \n", min_dist );

//-- Draw only "good" matches (i.e. whose distance is less than 3*min_dist )

std::vector< DMatch > good_matches;

for( int i = 0; i < descriptors_object.rows; i++ )

{ if( matches[i].distance < 1.5*min_dist )

{

good_matches.push_back( matches[i]);

}

Mat img_matches;

drawMatches( img_object, keypoints_object, img_scene, keypoints_scene,

good_matches, img_matches, Scalar::all(-1), Scalar::all(-1),

vector<char>(), DrawMatchesFlags::NOT_DRAW_SINGLE_POINTS );

//-- Localize the object from img_1 in img_2

std::vector<Point2f> obj;

std::vector<Point2f> scene;

for( size_t i = 0; i < good_matches.size(); i++ )

{

//-- Get the keypoints from the good matches

obj.push_back( keypoints_object[ good_matches[i].queryIdx ].pt );

scene.push_back( keypoints_scene[ good_matches[i].trainIdx ].pt );

}

Mat H = findHomography( obj, scene, CV_RANSAC );

//-- Get the corners from the image_1 ( the object to be "detected" )

std::vector<Point2f> obj_corners(4);

obj_corners[0] = Point(0,0); obj_corners[1] = Point( img_object.cols, 0 );

obj_corners[2] = Point( img_object.cols, img_object.rows ); obj_corners[3] = Point( 0, img_object.rows );

std::vector<Point2f> scene_corners(4);

perspectiveTransform( obj_corners, scene_corners, H);

//-- Draw lines between the corners (the mapped object in the scene - image_2 )

Point2f offset( (float)img_object.cols, 0);

line( img_matches, scene_corners[0] + offset, scene_corners[1] + offset, Scalar(0, 255, 0), 4 );

line( img_matches, scene_corners[1] + offset, scene_corners[2] + offset, Scalar( 0, 255, 0), 4 );

line( img_matches, scene_corners[2] + offset, scene_corners[3] + offset, Scalar( 0, 255, 0), 4 );

line( img_matches, scene_corners[3] + offset, scene_corners[0] + offset, Scalar( 0, 255, 0), 4 );

//-- Show detected matches

imshow( "Good Matches & Object detection", img_matches );

waitKey(0);

return ;

}

Applying homography matrix in OpenCV

// We need 4 corresponding 2D points(x,y) to calculate homography.
vector<Point2f> left_image;      // Stores 4 points(x,y) of the logo image. Here the four points are 4 corners of image.
vector<Point2f> right_image;    // stores 4 points that the user clicks(mouse left click) in the main image.

// Image containers for main and logo image
Mat imageMain;
Mat imageLogo;

// Function to add main image and transformed logo image and show final output.
// Icon image replaces the pixels of main image in this implementation.
void showFinal(Mat src1,Mat src2)
{

    Mat gray,gray_inv,src1final,src2final;
    cvtColor(src2,gray,CV_BGR2GRAY);
    threshold(gray,gray,0,255,CV_THRESH_BINARY);
    //adaptiveThreshold(gray,gray,255,ADAPTIVE_THRESH_MEAN_C,THRESH_BINARY,5,4);
    bitwise_not ( gray, gray_inv );
    src1.copyTo(src1final,gray_inv);
    src2.copyTo(src2final,gray);
    Mat finalImage = src1final+src2final;
    namedWindow( "output", WINDOW_AUTOSIZE );
    imshow("output",finalImage);
    cvWaitKey(0);

}

// Here we get four points from the user with left mouse clicks.
// On 5th click we output the overlayed image.
void on_mouse( int e, int x, int y, int d, void *ptr )
{
    if (e == EVENT_LBUTTONDOWN )
    {
        if(right_image.size() < 4 )
        {

            right_image.push_back(Point2f(float(x),float(y)));
            cout << x << " "<< y <<endl;
        }
        else
        {
            cout << " Calculating Homography " <<endl;
            // Deactivate callback
            cv::setMouseCallback("Display window", NULL, NULL);
            // once we get 4 corresponding points in both images calculate homography matrix
            Mat H = findHomography(  left_image,right_image,0 );
            Mat logoWarped;
            // Warp the logo image to change its perspective
            warpPerspective(imageLogo,logoWarped,H,imageMain.size() );
            showFinal(imageMain,logoWarped);

        }

    }
}

/*
How to run:
./main  ../images/homography/main.jpg ../images/homography/logo.jpg
*/
int homographyPerspective_Test( int argc, char** argv )
{
//  We need tow argumemts. "Main image" and "logo image"
    if( argc != 3)
    {
        cout <<" Usage: error" << endl;
        return -1;
    }


// Load images from arguments passed.
    imageMain = imread(argv[1], CV_LOAD_IMAGE_COLOR);
    imageLogo = imread(argv[2], CV_LOAD_IMAGE_COLOR);
// Push the 4 corners of the logo image as the 4 points for correspondence to calculate homography.
    left_image.push_back(Point2f(float(0),float(0)));
    left_image.push_back(Point2f(float(0),float(imageLogo.rows)));
    left_image.push_back(Point2f(float(imageLogo.cols),float(imageLogo.rows)));
    left_image.push_back(Point2f(float(imageLogo.cols),float(0)));



    namedWindow( "Display window", WINDOW_AUTOSIZE );// Create a window for display.
    imshow( "Display window", imageMain );


    setMouseCallback("Display window",on_mouse, NULL );


//  Press "Escape button" to exit
    while(1)
    {
        int key=cvWaitKey(10);
        if((char)key==(char)27) break;
    }


    return 0;
}

// We need 4 corresponding 2D points(x,y) to calculate homography.

vector<Point2f> left_image; // Stores 4 points(x,y) of the logo image. Here the four points are 4 corners of image.

vector<Point2f> right_image; // stores 4 points that the user clicks(mouse left click) in the main image.

// Image containers for main and logo image

Mat imageMain;

Mat imageLogo;

// Function to add main image and transformed logo image and show final output.

// Icon image replaces the pixels of main image in this implementation.

void showFinal(Mat src1,Mat src2)

{

Mat gray,gray_inv,src1final,src2final;

cvtColor(src2,gray,CV_BGR2GRAY);

threshold(gray,gray,0,255,CV_THRESH_BINARY);

//adaptiveThreshold(gray,gray,255,ADAPTIVE_THRESH_MEAN_C,THRESH_BINARY,5,4);

bitwise_not ( gray, gray_inv );

src1.copyTo(src1final,gray_inv);

src2.copyTo(src2final,gray);

Mat finalImage = src1final+src2final;

namedWindow( "output", WINDOW_AUTOSIZE );

imshow("output",finalImage);

cvWaitKey(0);

}

// Here we get four points from the user with left mouse clicks.

// On 5th click we output the overlayed image.

void on_mouse( int e, int x, int y, int d, void *ptr )

{

if (e == EVENT_LBUTTONDOWN )

{

if(right_image.size() < 4 )

{

right_image.push_back(Point2f(float(x),float(y)));

cout << x << " "<< y <<endl;

}

else

{

cout << " Calculating Homography " <<endl;

// Deactivate callback

cv::setMouseCallback("Display window", NULL, NULL);

// once we get 4 corresponding points in both images calculate homography matrix

Mat H = findHomography( left_image,right_image,0 );

Mat logoWarped;

// Warp the logo image to change its perspective

warpPerspective(imageLogo,logoWarped,H,imageMain.size() );

showFinal(imageMain,logoWarped);

}

How to run:

./main ../images/homography/main.jpg ../images/homography/logo.jpg

int homographyPerspective_Test( int argc, char** argv )

{

// We need tow argumemts. "Main image" and "logo image"

if( argc != 3)

{

cout <<" Usage: error" << endl;

return -1;

}

// Load images from arguments passed.

imageMain = imread(argv[1], CV_LOAD_IMAGE_COLOR);

imageLogo = imread(argv[2], CV_LOAD_IMAGE_COLOR);

// Push the 4 corners of the logo image as the 4 points for correspondence to calculate homography.

left_image.push_back(Point2f(float(0),float(0)));

left_image.push_back(Point2f(float(0),float(imageLogo.rows)));

left_image.push_back(Point2f(float(imageLogo.cols),float(imageLogo.rows)));

left_image.push_back(Point2f(float(imageLogo.cols),float(0)));

namedWindow( "Display window", WINDOW_AUTOSIZE );// Create a window for display.

imshow( "Display window", imageMain );

setMouseCallback("Display window",on_mouse, NULL );

// Press "Escape button" to exit

while(1)

{

int key=cvWaitKey(10);

if((char)key==(char)27) break;

}

return 0;

}

Examples of Dynamic Programming with C++ and Matlab

Leave a reply

In this tutorial, I will give you examples of using dynamic programming for solving the following problems:

1)Minimum number of coins for summing X.

int c[]={1,3,4}; // value of coins
int k=3; //number of cpins that we have

/*
The d[20] array is for memorization, for instance if we computed f_coin_problem_memo(3),
we put this value in d[3] and we don't compute it again, we just use d[3] instead of f_coin_problem_memo(3)
*/
int d[20];
int f_coin_problem_memo(int x)
{
    if (x < 0) return 1e9;
    if (x == 0) return 0;
    if (d[x]) return d[x];
    int u = 1e9;
    int u_old = 1e9;
    for (int i = 0; i < k; i++)
    {
        u_old=u;
        u = std::min(u, f_coin_problem_memo(x-c[i]) +1);
    }
    d[x] = u;
    return d[x];
}

int c[]={1,3,4}; // value of coins

int k=3; //number of cpins that we have

The d[20] array is for memorization, for instance if we computed f_coin_problem_memo(3),

we put this value in d[3] and we don't compute it again, we just use d[3] instead of f_coin_problem_memo(3)

int d[20];

int f_coin_problem_memo(int x)

{

if (x < 0) return 1e9;

if (x == 0) return 0;

if (d[x]) return d[x];

int u = 1e9;

int u_old = 1e9;

for (int i = 0; i < k; i++)

{

u_old=u;

u = std::min(u, f_coin_problem_memo(x-c[i]) +1);

}

d[x] = u;

return d[x];

}

The minimum number of coins from the set {1,3,4} to sum up value of 12 is:
3

1 2	The minimum number of coins from the set {1,3,4} to sum up value of 12 is: 3

2)The most (least) costly path on a grid (dynamic time warping).

finding a path in an n × n grid from the upper-left corner to the lower-right corner. We are only allowed
 to move right or down.

1,1  1,2 .... 1,X
2,1  2,2 .... 2,X
.
.
.
Y,1  Y,2 .... Y,X


3 7 9 2 7
9 8 3 5 5
1 7 9 8 5
3 8 6 4 10
6 3 9 7 8

*/

int grid[5][5]={{3,7, 9, 2, 7},{9, 8, 3 ,5 ,5},{1, 7 ,9 ,8 ,5},{3, 8, 6 ,4 ,10},{6, 3 ,9 ,7 ,8}};
int memo[5][5]{{-1,-1,-1,-1,-1},{-1,-1,-1,-1,-1},{-1,-1,-1,-1,-1},{-1,-1,-1,-1,-1},{-1,-1,-1,-1,-1}};

std::vector<int>x_solution,y_solution;

int f_paths_in_grid(int y,int x)
{
    if((x==0)&& (y==0))
    {
        memo[y][x]=grid[0][0];
    }
    else if (x==0)
    {
        if(memo[y][x]<0)
        {
            memo[y][x] =f_paths_in_grid(y-1,x) +grid[y][x];
        }
    }
    else if(y==0)
    {
       if(memo[y][x]<0)
       {
            memo[y][x]=f_paths_in_grid(y,x-1) +grid[y][x];
       }
    }
   else if(memo[y][x]<0)
   {
        memo[y][x]=std::max(f_paths_in_grid(y,x-1)  ,f_paths_in_grid(y-1,x)  )+grid[y][x];
    }
    return  memo[y][x];
}

finding a path in an n × n grid from the upper-left corner to the lower-right corner. We are only allowed

to move right or down.

1,1 1,2 .... 1,X

2,1 2,2 .... 2,X

Y,1 Y,2 .... Y,X

3 7 9 2 7

9 8 3 5 5

1 7 9 8 5

3 8 6 4 10

6 3 9 7 8

int grid[5][5]={{3,7, 9, 2, 7},{9, 8, 3 ,5 ,5},{1, 7 ,9 ,8 ,5},{3, 8, 6 ,4 ,10},{6, 3 ,9 ,7 ,8}};

int memo[5][5]{{-1,-1,-1,-1,-1},{-1,-1,-1,-1,-1},{-1,-1,-1,-1,-1},{-1,-1,-1,-1,-1},{-1,-1,-1,-1,-1}};

std::vector<int>x_solution,y_solution;

int f_paths_in_grid(int y,int x)

{

if((x==0)&& (y==0))

{

memo[y][x]=grid[0][0];

}

else if (x==0)

{

if(memo[y][x]<0)

{

memo[y][x] =f_paths_in_grid(y-1,x) +grid[y][x];

}

else if(y==0)

{

if(memo[y][x]<0)

{

memo[y][x]=f_paths_in_grid(y,x-1) +grid[y][x];

}

else if(memo[y][x]<0)

{

memo[y][x]=std::max(f_paths_in_grid(y,x-1) ,f_paths_in_grid(y-1,x) )+grid[y][x];

}

return memo[y][x];

}

maximum sum on a path from the upper-left corner to the lower-right corner is:
67
3|10|19|21|28|
---------------
12|20|23|28|33|
---------------
13|27|36|44|49|
---------------
16|35|42|48|59|
---------------
22|38|51|58|67|
---------------

maximum sum on a path from the upper-left corner to the lower-right corner is:

3|10|19|21|28|

---------------

12|20|23|28|33|

---------------

13|27|36|44|49|

---------------

16|35|42|48|59|

---------------

22|38|51|58|67|

---------------

3)Levenshtein edit distance.

int LevenshteinDistance(std::string s, std::string t )
{
  int cost;

  /* base case: empty strings */
  if (s.length() == 0) return t.length();
  if (t.length() == 0) return s.length();

  /* test if last characters of the strings match */
  if (s.at(s.length()-1)  == t.at(t.length()-1) )
      cost = 0;
  else
      cost = 1;

  /* return minimum of delete char from s, delete char from t, and delete char from both */

  return std::min(   std::min(LevenshteinDistance(s.substr(0, s.size()-1 ),t ) + 1,LevenshteinDistance(s  , t.substr(0, t.size()-1 )) + 1 )    ,  LevenshteinDistance(s.substr(0, s.size()-1 ), t.substr(0, t.size()-1 )) + cost );
}

int LevenshteinDistance(std::string s, std::string t )

{

int cost;

/* base case: empty strings */

if (s.length() == 0) return t.length();

if (t.length() == 0) return s.length();

/* test if last characters of the strings match */

if (s.at(s.length()-1) == t.at(t.length()-1) )

cost = 0;

else

cost = 1;

/* return minimum of delete char from s, delete char from t, and delete char from both */

return std::min( std::min(LevenshteinDistance(s.substr(0, s.size()-1 ),t ) + 1,LevenshteinDistance(s , t.substr(0, t.size()-1 )) + 1 ) , LevenshteinDistance(s.substr(0, s.size()-1 ), t.substr(0, t.size()-1 )) + cost );

}

The Levenshtein distance between saturday and sunday is:
3

1 2	The Levenshtein distance between saturday and sunday is: 3

4)Seam Carving. I have written a tutorial on that here and the recursive part is in the following lines:

for x=2:rows
    for y=1:cols  
        M(x+1,y+1)=score_matrix(x,y)+ min( [M(x-1+1,y-1 +1) , M(x-1+1,y+1), M(x-1+1,y+1+1) ] );
    end
end

for x=2:rows

for y=1:cols

M(x+1,y+1)=score_matrix(x,y)+ min( [M(x-1+1,y-1 +1) , M(x-1+1,y+1), M(x-1+1,y+1+1) ] );

end

The examples are taken from “Competitive Programmer’s Handbook” written by Antti Laaksonen.

Peak Signal-to-Noise Ratio (PSNR) in Image using OpenCV and Matlab

Leave a reply

\(\)\(\)

Peak signal-to-noise ratio (PSNR) shows the ratio between the maximum possible power of a signal and the power of the same image with noise. PSNR is usually expressed in logarithmic decibel scale.

\( MSE =1/m*n \sum_{i=0}^{m-1} \sum_{j=0}^{n-1} [ Image( i,j) -NoisyImage( i,j) ] ^2 \)
\( PSNR =10* \log_{10} \Bigg(MAXI^2/MSE\Bigg) \)

MSE is Mean Square Error and MAXI is the maximum possible pixel value of the image. For instance if the image is uint8, it will be 255.

ref = imread('pout.tif');
A = imnoise(ref,'salt & pepper', 0.02);
%%%%%%%%%%%%%%%%or we can apply gaussian noise%%%%%%%%%%%%%%%%
%A = imnoise(ref, 'gaussian', 0, 0.003);

[M,N]=size(ref);

SE = (double(A) - double(ref)) .^ 2; % Squared Error
MSE=sum( sum( SE  ) ) /(M*N); %Mean Square Error
MAXI=255;
fprintf('\n PSNR is:')
10*log10((MAXI^2)/MSE)

figure,imshow(A, []);
figure,imshow(ref, []);
figure,imshow(SE, []);


%%%%%%%%%%%%%%%%% using Matlab API %%%%%%%%%%%%%%%%
%MSE = immse(A, ref)
[peaksnr, snr] = psnr(A, ref);
fprintf('\n The Peak-SNR value is %0.4f', peaksnr);
fprintf('\n The SNR value is %0.4f \n', snr);

ref = imread('pout.tif');

A = imnoise(ref,'salt & pepper', 0.02);

%%%%%%%%%%%%%%%%or we can apply gaussian noise%%%%%%%%%%%%%%%%

%A = imnoise(ref, 'gaussian', 0, 0.003);

[M,N]=size(ref);

SE = (double(A) - double(ref)) .^ 2; % Squared Error

MSE=sum( sum( SE ) ) /(M*N); %Mean Square Error

MAXI=255;

fprintf('\n PSNR is:')

10*log10((MAXI^2)/MSE)

figure,imshow(A, []);

figure,imshow(ref, []);

figure,imshow(SE, []);

%%%%%%%%%%%%%%%%% using Matlab API %%%%%%%%%%%%%%%%

%MSE = immse(A, ref)

[peaksnr, snr] = psnr(A, ref);

fprintf('\n The Peak-SNR value is %0.4f', peaksnr);

fprintf('\n The SNR value is %0.4f \n', snr);

double getPSNR(const Mat& I1, const Mat& I2)
{
 Mat s1;
 absdiff(I1, I2, s1);       // |I1 - I2|
 s1.convertTo(s1, CV_32F);  // cannot make a square on 8 bits
 s1 = s1.mul(s1);           // |I1 - I2|^2

 Scalar s = sum(s1);         // sum elements per channel

 double sse = s.val[0] + s.val[1] + s.val[2]; // sum channels

 if( sse <= 1e-10) // for small values return zero
     return 0;
 else
 {
     double  mse =sse /(double)(I1.channels() * I1.total());
     double psnr = 10.0*log10((255*255)/mse);
     return psnr;
 }
}

double getPSNR(const Mat& I1, const Mat& I2)

{

Mat s1;

absdiff(I1, I2, s1); // |I1 - I2|

s1.convertTo(s1, CV_32F); // cannot make a square on 8 bits

s1 = s1.mul(s1); // |I1 - I2|^2

Scalar s = sum(s1); // sum elements per channel

double sse = s.val[0] + s.val[1] + s.val[2]; // sum channels

if( sse <= 1e-10) // for small values return zero

return 0;

else

{

double mse =sse /(double)(I1.channels() * I1.total());

double psnr = 10.0*log10((255*255)/mse);

return psnr;

}

Ref [1],[2], [3], 4

Creating Hybrid Image with Matlab

Leave a reply

In this tutorial, I’m gonna show you how to create hybrid images . Basically, you need to filter your first image with high pass filter and the second image with a low pass filter and then add theses two images. The following code shows how to do that. The full code is available at my GitHub.

first_image_path='../images/obama.jpg';
second_image_path='../images/michelle_obama.jpg';

first_image=im2double( rgb2gray( imread(first_image_path) ) );
second_image=im2double( rgb2gray( imread(second_image_path) ) );

sigma=6;
gaussian_dimension=3*sigma*2+1;


first_image_high_pass_filtered=highPassFilter(first_image,gaussian_dimension);


figure('Name','first_image_high_pass_filtered'), imshow(first_image_high_pass_filtered,[]);

sigma=1;
gaussian_dimension=3*sigma*2+1;
second_image_low_pass_filtered=lowPassFilter(second_image,gaussian_dimension);

figure('Name','second_image_low_pass_filtered'), imshow(second_image_low_pass_filtered,[]);

final_image=second_image_low_pass_filtered+first_image_high_pass_filtered;
figure('Name','second_image_low_pass_filtered'), imshow(final_image, []);

first_image_path='../images/obama.jpg';

second_image_path='../images/michelle_obama.jpg';

first_image=im2double( rgb2gray( imread(first_image_path) ) );

second_image=im2double( rgb2gray( imread(second_image_path) ) );

sigma=6;

gaussian_dimension=3*sigma*2+1;

first_image_high_pass_filtered=highPassFilter(first_image,gaussian_dimension);

figure('Name','first_image_high_pass_filtered'), imshow(first_image_high_pass_filtered,[]);

sigma=1;

gaussian_dimension=3*sigma*2+1;

second_image_low_pass_filtered=lowPassFilter(second_image,gaussian_dimension);

figure('Name','second_image_low_pass_filtered'), imshow(second_image_low_pass_filtered,[]);

final_image=second_image_low_pass_filtered+first_image_high_pass_filtered;

figure('Name','second_image_low_pass_filtered'), imshow(final_image, []);

Stitching image using SIFT and Homography

Leave a reply

This Matlab tutorial I use SIFT, RANSAC, and homography to find corresponding points between two images. Here I have used vlfeat to find SIFT features.

Full code is available at my GitHub repository

Major steps are:

0.Adding vlfeat to your Matlab workspace:

run('<path_to_vlfeat>/toolbox/vl_setup')

1	run('<path_to_vlfeat>/toolbox/vl_setup')

1.Detect key points and extract descriptors. In the image below you can see some SIFT key points and their descriptor.

image_left=imread('images/image-left.jpg');
scale=0.20;
image_left=imresize(image_left,scale);
image_left=single(rgb2gray(image_left)) ;
[f_image_left,d_image_left] = vl_sift(image_left) ;

image_left=imread('images/image-left.jpg');

scale=0.20;

image_left=imresize(image_left,scale);

image_left=single(rgb2gray(image_left)) ;

[f_image_left,d_image_left] = vl_sift(image_left) ;

2.Matching features:

[dims, image_right_features]=size(d_image_right);
[dims, image_left_features]=size(d_image_left);

M=zeros(image_right_features,image_left_features);

for i=1:image_right_features
  for j=1:image_left_features


    diff = double(d_image_right(:,i)) - double(d_image_left(:,j)) ;
    distance = sqrt(diff' * diff);
    M(i,j)=distance ;
  end
end

[dims, image_right_features]=size(d_image_right);

[dims, image_left_features]=size(d_image_left);

M=zeros(image_right_features,image_left_features);

for i=1:image_right_features

for j=1:image_left_features

diff = double(d_image_right(:,i)) - double(d_image_left(:,j)) ;

distance = sqrt(diff' * diff);

M(i,j)=distance ;

end

3.Pruning features

In this step, I only took first k top matches.

number_of_top_matches_to_select=90;

M_Vec=sort(M(:));
M=(M < M_Vec( number_of_top_matches_to_select ) );

1 2	M_Vec=sort(M(:)); M=(M < M_Vec( number_of_top_matches_to_select ) );

4.Estimating transformation using RANSAC

I used RANSAC to estimate the transformation.

for i=1:M_Rows
  for j=1:M_Cols
    if M(i,j)>0
      x_right=round(f_image_right(1,i));
      y_right=round(f_image_right(2,i));
      x_left=round(f_image_left(1,j));
      y_left=round(f_image_left(2,j));
      left_matches=vertcat(left_matches, [x_left,y_left]);
      right_matches=vertcat(right_matches, [x_right,y_right]);
    end
  end
end

p=99/100;
e=50/100;
[s,dim]=size(left_matches);
s=3
N= (log (1-p))/(log(1-(1-e)^s))

iter=N;%we choose this so with probability of 99% at least one random sample is free from ouliers
num=3;%minmum number of samples to make the model i.e for line is 2 and for affine is 3
threshDist=5; %5 pixel
[best_number_of_inliers,best_M,index_of_matches]=affine_ransac_est(left_matches,right_matches,num,iter,threshDist);

for i=1:M_Rows

for j=1:M_Cols

if M(i,j)>0

x_right=round(f_image_right(1,i));

y_right=round(f_image_right(2,i));

x_left=round(f_image_left(1,j));

y_left=round(f_image_left(2,j));

left_matches=vertcat(left_matches, [x_left,y_left]);

right_matches=vertcat(right_matches, [x_right,y_right]);

end

p=99/100;

e=50/100;

[s,dim]=size(left_matches);

s=3

N= (log (1-p))/(log(1-(1-e)^s))

iter=N;%we choose this so with probability of 99% at least one random sample is free from ouliers

num=3;%minmum number of samples to make the model i.e for line is 2 and for affine is 3

threshDist=5; %5 pixel

[best_number_of_inliers,best_M,index_of_matches]=affine_ransac_est(left_matches,right_matches,num,iter,threshDist);

5.Compute optimal transformation

Finally, I used the least square method to find the affine transformation between two images.

best_M=affine_least_square_generalized(best_left_matches, best_right_matches);

1	best_M=affine_least_square_generalized(best_left_matches, best_right_matches);

Analytic distance metric for Gaussian mixture models

Leave a reply

In many applications, you need to compare two or more data sets with each other to see how much they are similar or different. For instance, you have measured the height of men and women in Japan and Netherlands and now you like to know how much they are different.

\( \)

Two commonly used method for measuring distances are the Kullback-Liebler divergence and the Bhattacharyya distance.

KL divergence (Kullback-Liebler divergence) measures the difference between two probability distributions p and q.

\begin{equation} \label{Kullback_Liebler_divergence}
D_{KL}(p||q)=\int_{-\infty}^\infty p(x)\log\frac{p(x)}{q(x)} \,\mathrm{d}x
\end{equation}

But it only works if your data is made of a single Gaussian and it is not applicable If your data is made of a mixture of Gaussians.

Sfikas et al [1] have extended the Kullback Liebler divergence for GMM and proposed a distance metric using the values \(\mu ,\Sigma,\pi \) for each one of the two distributions in the following form:

\begin{equation} \label{analytical_Kullback_Liebler_divergence}
C2(p||q)=-\log \large[ \frac{2\sum_{i,j}\pi_{i}\pi_{j}^{\prime} \sqrt{ \frac{|V_{ij}|}{e^{k_{ij}}|\sum_{i}| |\sum_{j}^{\prime}|} } }
{
\sum_{i,j}\pi_{i}\pi_{j} \sqrt{ \frac{|V_{ij}|}{e^{k_{ij}}|\sum_{i}| |\sum_{j}|} }+
\sum_{i,j}\pi_{i}^{\prime}\pi_{j}^{\prime} \sqrt{ \frac{|V_{ij}|}{e^{k_{ij}}|\sum_{i}^{\prime}| |\sum_{j}^{\prime}|} }
}
\large]
\end{equation}
Where:
\begin{equation}\label{Kullback_Liebler_divergence_Details1}
V_{ij}=(\Sigma_{i}^{-1} +\Sigma_{j}^{-1})^{-1}
\end{equation}
and
\begin{equation}\label{Kullback_Liebler_divergence_Details2}
K_{ij}=\mu_{i}^{T}\Sigma_{i}^{-1}(\mu_{i}-\mu_{j}^{\prime})+\mu_{j}^{\prime T}\Sigma_{j}^{\prime -1}(\mu_{j}^{\prime}-\mu_{i})
\end{equation}

Code in matlab:

function distance = GetDistanceBetweenTwoDists(means1, pis1, covs1, means2, pis2, covs2)
    %V11 and K11 
    for i=1:size(means1,1)
      for j=1:size(means1,1)
         cov1= covs1(i,:,:);
         cov2= covs1(j,:,:);
         cov1 = reshape(cov1,size(covs1,2), size(covs1,2));
         cov2 = reshape(cov2,size(covs1,2), size(covs1,2)); 
         V11(i,j) = det (inv(inv(cov1)+inv(cov2)));
         K11(i,j) = means1(i,:)*inv(cov1)*(means1(i,:)-means1(j,:))'+means1(j,:)*inv(cov2)*(means1(j,:)-means1(i,:))';
      end
    end

    %V12 and K12
    for i=1:size(means1,1)
      for j=1:size(means2,1)
         cov1= covs1(i,:,:);
         cov2= covs2(j,:,:);
         cov1 = reshape(cov1,size(covs1,2), size(covs1,2));
         cov2 = reshape(cov2,size(covs1,2), size(covs1,2));
         V12(i,j) = det (inv(inv(cov1)+inv(cov2)));
         K12(i,j) = means1(i,:)*inv(cov1)*(means1(i,:)-means2(j,:))'+means2(j,:)*inv(cov2)*(means2(j,:)-means1(i,:))';
      end
    end
    %end

    %V22 and K22
    for i=1:size(means2,1)
      for j=1:size(means2,1)
         cov1= covs2(i,:,:);
         cov2= covs2(j,:,:);
         cov1 = reshape(cov1,size(covs1,2), size(covs1,2));
         cov2 = reshape(cov2,size(covs1,2), size(covs1,2));
         V22(i,j) = det (inv(inv(cov1)+inv(cov2)));
         K22(i,j) = means2(i,:)*inv(cov1)*(means2(i,:)-means2(j,:))'+means2(j,:)*inv(cov2)*(means2(j,:)-means2(i,:))';
      end
    end
    %end

    %Sum11
    Sum11 = 0;

    for i=1:size(means1,1)
      for j=1:size(means1,1)
         cov1= covs1(i,:,:);
         cov2= covs1(j,:,:);
         cov1 = reshape(cov1,size(covs1,2), size(covs1,2));
         cov2 = reshape(cov2,size(covs1,2), size(covs1,2));
         Sum11 = Sum11 + pis1(i)*pis1(j)*sqrt(V11(i,j)/(exp(K11(i,j))*det(cov1)*det(cov2)));
      end
    end

    %Sum12
    Sum12 = 0;

    for i=1:size(means1,1)
      for j=1:size(means2,1)
         cov1= covs1(i,:,:);
         cov2= covs2(j,:,:);
         cov1 = reshape(cov1,size(covs1,2), size(covs1,2));
         cov2 = reshape(cov2,size(covs1,2), size(covs1,2));
         Sum12 = Sum12 + pis1(i)*pis2(j)*sqrt(V12(i,j)/(exp(K12(i,j))*det(cov1)*det(cov2)));
      end
    end

    %Sum22
    Sum22 = 0;

    for i=1:size(means2,1)
      for j=1:size(means2,1)
         cov1= covs2(i,:,:);
         cov2= covs2(j,:,:);
         cov1 = reshape(cov1,size(covs1,2), size(covs1,2));
         cov2 = reshape(cov2,size(covs1,2), size(covs1,2));
         Sum22 = Sum22 + pis2(i)*pis2(j)*sqrt(V22(i,j)/(exp(K22(i,j))*det(cov1)*det(cov2)));
      end
    end

distance = -log(2*Sum12/(Sum11+Sum22));
%distance = 2*Sum12/(Sum11+Sum22);

function distance = GetDistanceBetweenTwoDists(means1, pis1, covs1, means2, pis2, covs2)

%V11 and K11

for i=1:size(means1,1)

for j=1:size(means1,1)

cov1= covs1(i,:,:);

cov2= covs1(j,:,:);

cov1 = reshape(cov1,size(covs1,2), size(covs1,2));

cov2 = reshape(cov2,size(covs1,2), size(covs1,2));

V11(i,j) = det (inv(inv(cov1)+inv(cov2)));

K11(i,j) = means1(i,:)*inv(cov1)*(means1(i,:)-means1(j,:))'+means1(j,:)*inv(cov2)*(means1(j,:)-means1(i,:))';

end

%V12 and K12

for i=1:size(means1,1)

for j=1:size(means2,1)

cov1= covs1(i,:,:);

cov2= covs2(j,:,:);

cov1 = reshape(cov1,size(covs1,2), size(covs1,2));

cov2 = reshape(cov2,size(covs1,2), size(covs1,2));

V12(i,j) = det (inv(inv(cov1)+inv(cov2)));

K12(i,j) = means1(i,:)*inv(cov1)*(means1(i,:)-means2(j,:))'+means2(j,:)*inv(cov2)*(means2(j,:)-means1(i,:))';

end

%end

%V22 and K22

for i=1:size(means2,1)

for j=1:size(means2,1)

cov1= covs2(i,:,:);

cov2= covs2(j,:,:);

cov1 = reshape(cov1,size(covs1,2), size(covs1,2));

cov2 = reshape(cov2,size(covs1,2), size(covs1,2));

V22(i,j) = det (inv(inv(cov1)+inv(cov2)));

K22(i,j) = means2(i,:)*inv(cov1)*(means2(i,:)-means2(j,:))'+means2(j,:)*inv(cov2)*(means2(j,:)-means2(i,:))';

end

%end

%Sum11

Sum11 = 0;

for i=1:size(means1,1)

for j=1:size(means1,1)

cov1= covs1(i,:,:);

cov2= covs1(j,:,:);

cov1 = reshape(cov1,size(covs1,2), size(covs1,2));

cov2 = reshape(cov2,size(covs1,2), size(covs1,2));

Sum11 = Sum11 + pis1(i)*pis1(j)*sqrt(V11(i,j)/(exp(K11(i,j))*det(cov1)*det(cov2)));

end

%Sum12

Sum12 = 0;

for i=1:size(means1,1)

for j=1:size(means2,1)

cov1= covs1(i,:,:);

cov2= covs2(j,:,:);

cov1 = reshape(cov1,size(covs1,2), size(covs1,2));

cov2 = reshape(cov2,size(covs1,2), size(covs1,2));

Sum12 = Sum12 + pis1(i)*pis2(j)*sqrt(V12(i,j)/(exp(K12(i,j))*det(cov1)*det(cov2)));

end

%Sum22

Sum22 = 0;

for i=1:size(means2,1)

for j=1:size(means2,1)

cov1= covs2(i,:,:);

cov2= covs2(j,:,:);

cov1 = reshape(cov1,size(covs1,2), size(covs1,2));

cov2 = reshape(cov2,size(covs1,2), size(covs1,2));

Sum22 = Sum22 + pis2(i)*pis2(j)*sqrt(V22(i,j)/(exp(K22(i,j))*det(cov1)*det(cov2)));

end

distance = -log(2*Sum12/(Sum11+Sum22));

%distance = 2*Sum12/(Sum11+Sum22);

Lucas–Kanade method optical flow with MATLAB

Leave a reply

In this tutorial, I will show you how to estimate optical flow based on Lucas–Kanade method. This project has the following scripts: Optical_flow_estimation, myFlow, myWarp, computeColor, flowToColor.

The myFlow does the main job, it gets two images and a window length (patch length) and a threshold for accepting the optical flow. In the following, you see the myFlow. You can uncomment figure function calls to see output result of each step. The other scripts are just for visualization. you can access the code and image in my Github repository.

first image

second image

estimated optical flow

color map of optical flow

Displaying Adjusted Frequencies of signal with Fast Fourier transform

Leave a reply

In this tutorial, I will show you how to get Fast Fourier transform of a signal and then correctly display the signal. Link to the code in my Github repository.

%% Creating signal
Fs=1000; %Sampling fresquncy
T_inc=1/Fs; %Time increament
T_measure=1.5; %Duration of measurment
time=0:T_inc:T_measure - T_inc; %Vector contains sampling time
L=length(time); %Lenght of signal


%% Sum of a 60Hz and 120 Hz sinusoidal
f1=120;
f2=60;
phi1=1.5;
phi2=2.4;
amplitude1=1.5;
amplitude2=2.2;
y=amplitude1*sin(2.*pi.*f1*time+phi1) + amplitude2*sin(2.*pi.*f2*time+phi2);

%% Creating signal

Fs=1000; %Sampling fresquncy

T_inc=1/Fs; %Time increament

T_measure=1.5; %Duration of measurment

time=0:T_inc:T_measure - T_inc; %Vector contains sampling time

L=length(time); %Lenght of signal

%% Sum of a 60Hz and 120 Hz sinusoidal

f1=120;

f2=60;

phi1=1.5;

phi2=2.4;

amplitude1=1.5;

amplitude2=2.2;

y=amplitude1*sin(2.*pi.*f1*time+phi1) + amplitude2*sin(2.*pi.*f2*time+phi2);

Signal in time domain

Getting Fourier transform of the signal:

signal_FFT=fft(y,L);

1	signal_FFT=fft(y,L);

The signal in frequency domain (Bins)

signal_FFT=fft(y,L)/L;
amplitude=2*abs( signal_FFT(1: floor(L/2) +1) );
frequency=Fs/2*linspace(0,1, floor(L/2) +1);

signal_FFT=fft(y,L)/L;

amplitude=2*abs( signal_FFT(1: floor(L/2) +1) );

frequency=Fs/2*linspace(0,1, floor(L/2) +1);

The signal in frequency domain consists of 60Hz and 120Hz.