2016年6月2日 星期四

camera calibration


看了Learning OpenCV的第11章Camera Models and Calibration,知道每個鏡頭都有鏡面曲率,這會影響到攝影的結果,讓原本是方形的物體成像後變成桶狀,這對影像處理影響深遠,所以有必要做相機矯正 (Camera Calibration)
這裡使用一個鏡頭和一張棋盤,如圖所示,擷取多張不同角度的棋盤影像。
接著程式會找尋影像中的「角」,判斷是否為我們所指定的棋盤,並且在影像上做標記。
最後使用相機內部參數來矯正影像。





Camera Calibration

Goal

In this section,
  • We will learn about distortions in camera, intrinsic and extrinsic parameters of camera etc.
  • We will learn to find these parameters, undistort images etc.

Basics

Today’s cheap pinhole cameras introduces a lot of distortion to images. Two major distortions are radial distortion and tangential distortion.
Due to radial distortion, straight lines will appear curved. Its effect is more as we move away from the center of image. For example, one image is shown below, where two edges of a chess board are marked with red lines. But you can see that border is not a straight line and doesn’t match with the red line. All the expected straight lines are bulged out. Visit Distortion (optics) for more details.
Radial Distortion
This distortion is solved as follows:
x_{corrected} = x( 1 + k_1 r^2 + k_2 r^4 + k_3 r^6) \\
y_{corrected} = y( 1 + k_1 r^2 + k_2 r^4 + k_3 r^6)
Similarly, another distortion is the tangential distortion which occurs because image taking lense is not aligned perfectly parallel to the imaging plane. So some areas in image may look nearer than expected. It is solved as below:
x_{corrected} = x + [ 2p_1xy + p_2(r^2+2x^2)] \\
y_{corrected} = y + [ p_1(r^2+ 2y^2)+ 2p_2xy]
In short, we need to find five parameters, known as distortion coefficients given by:
Distortion \; coefficients=(k_1 \hspace{10pt} k_2 \hspace{10pt} p_1 \hspace{10pt} p_2 \hspace{10pt} k_3)
In addition to this, we need to find a few more information, like intrinsic and extrinsic parameters of a camera. Intrinsic parameters are specific to a camera. It includes information like focal length (f_x,f_y), optical centers (c_x, c_y) etc. It is also called camera matrix. It depends on the camera only, so once calculated, it can be stored for future purposes. It is expressed as a 3x3 matrix:
camera \; matrix = \left [ \begin{matrix}   f_x & 0 & c_x \\  0 & f_y & c_y \\   0 & 0 & 1 \end{matrix} \right ]
Extrinsic parameters corresponds to rotation and translation vectors which translates a coordinates of a 3D point to a coordinate system.
For stereo applications, these distortions need to be corrected first. To find all these parameters, what we have to do is to provide some sample images of a well defined pattern (eg, chess board). We find some specific points in it ( square corners in chess board). We know its coordinates in real world space and we know its coordinates in image. With these data, some mathematical problem is solved in background to get the distortion coefficients. That is the summary of the whole story. For better results, we need at least 10 test patterns.

Code

As mentioned above, we need atleast 10 test patterns for camera calibration. OpenCV comes with some images of chess board (see samples/cpp/left01.jpg -- left14.jpg), so we will utilize it. For sake of understanding, consider just one image of a chess board. Important input datas needed for camera calibration is a set of 3D real world points and its corresponding 2D image points. 2D image points are OK which we can easily find from the image. (These image points are locations where two black squares touch each other in chess boards)
What about the 3D points from real world space? Those images are taken from a static camera and chess boards are placed at different locations and orientations. So we need to know (X,Y,Z)values. But for simplicity, we can say chess board was kept stationary at XY plane, (so Z=0 always) and camera was moved accordingly. This consideration helps us to find only X,Y values. Now for X,Y values, we can simply pass the points as (0,0), (1,0), (2,0), ... which denotes the location of points. In this case, the results we get will be in the scale of size of chess board square. But if we know the square size, (say 30 mm), and we can pass the values as (0,0),(30,0),(60,0),..., we get the results in mm. (In this case, we don’t know square size since we didn’t take those images, so we pass in terms of square size).
3D points are called object points and 2D image points are called image points.

Setup

So to find pattern in chess board, we use the function, cv2.findChessboardCorners(). We also need to pass what kind of pattern we are looking, like 8x8 grid, 5x5 grid etc. In this example, we use 7x6 grid. (Normally a chess board has 8x8 squares and 7x7 internal corners). It returns the corner points and retval which will be True if pattern is obtained. These corners will be placed in an order (from left-to-right, top-to-bottom)
See also
This function may not be able to find the required pattern in all the images. So one good option is to write the code such that, it starts the camera and check each frame for required pattern. Once pattern is obtained, find the corners and store it in a list. Also provides some interval before reading next frame so that we can adjust our chess board in different direction. Continue this process until required number of good patterns are obtained. Even in the example provided here, we are not sure out of 14 images given, how many are good. So we read all the images and take the good ones.
See also
Instead of chess board, we can use some circular grid, but then use the functioncv2.findCirclesGrid() to find the pattern. It is said that less number of images are enough when using circular grid.
Once we find the corners, we can increase their accuracy using cv2.cornerSubPix(). We can also draw the pattern using cv2.drawChessboardCorners(). All these steps are included in below code:
import numpy as np
import cv2
import glob

# termination criteria
criteria = (cv2.TERM_CRITERIA_EPS + cv2.TERM_CRITERIA_MAX_ITER, 30, 0.001)

# prepare object points, like (0,0,0), (1,0,0), (2,0,0) ....,(6,5,0)
objp = np.zeros((6*7,3), np.float32)
objp[:,:2] = np.mgrid[0:7,0:6].T.reshape(-1,2)

# Arrays to store object points and image points from all the images.
objpoints = [] # 3d point in real world space
imgpoints = [] # 2d points in image plane.

images = glob.glob('*.jpg')

for fname in images:
    img = cv2.imread(fname)
    gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)

    # Find the chess board corners
    ret, corners = cv2.findChessboardCorners(gray, (7,6),None)

    # If found, add object points, image points (after refining them)
    if ret == True:
        objpoints.append(objp)

        corners2 = cv2.cornerSubPix(gray,corners,(11,11),(-1,-1),criteria)
        imgpoints.append(corners2)

        # Draw and display the corners
        img = cv2.drawChessboardCorners(img, (7,6), corners2,ret)
        cv2.imshow('img',img)
        cv2.waitKey(500)

cv2.destroyAllWindows()
One image with pattern drawn on it is shown below:
Calibration Pattern

Calibration

So now we have our object points and image points we are ready to go for calibration. For that we use the function, cv2.calibrateCamera(). It returns the camera matrix, distortion coefficients, rotation and translation vectors etc.
ret, mtx, dist, rvecs, tvecs = cv2.calibrateCamera(objpoints, imgpoints, gray.shape[::-1],None,None)

Undistortion

We have got what we were trying. Now we can take an image and undistort it. OpenCV comes with two methods, we will see both. But before that, we can refine the camera matrix based on a free scaling parameter using cv2.getOptimalNewCameraMatrix(). If the scaling parameter alpha=0, it returns undistorted image with minimum unwanted pixels. So it may even remove some pixels at image corners. If alpha=1, all pixels are retained with some extra black images. It also returns an image ROI which can be used to crop the result.
So we take a new image (left12.jpg in this case. That is the first image in this chapter)
img = cv2.imread('left12.jpg')
h,  w = img.shape[:2]
newcameramtx, roi=cv2.getOptimalNewCameraMatrix(mtx,dist,(w,h),1,(w,h))

1. Using cv2.undistort()

This is the shortest path. Just call the function and use ROI obtained above to crop the result.
# undistort
dst = cv2.undistort(img, mtx, dist, None, newcameramtx)

# crop the image
x,y,w,h = roi
dst = dst[y:y+h, x:x+w]
cv2.imwrite('calibresult.png',dst)

2. Using remapping

This is curved path. First find a mapping function from distorted image to undistorted image. Then use the remap function.
# undistort
mapx,mapy = cv2.initUndistortRectifyMap(mtx,dist,None,newcameramtx,(w,h),5)
dst = cv2.remap(img,mapx,mapy,cv2.INTER_LINEAR)

# crop the image
x,y,w,h = roi
dst = dst[y:y+h, x:x+w]
cv2.imwrite('calibresult.png',dst)
Both the methods give the same result. See the result below:
Calibration Result
You can see in the result that all the edges are straight.
Now you can store the camera matrix and distortion coefficients using write functions in Numpy (np.savez, np.savetxt etc) for future uses.

Re-projection Error

Re-projection error gives a good estimation of just how exact is the found parameters. This should be as close to zero as possible. Given the intrinsic, distortion, rotation and translation matrices, we first transform the object point to image point using cv2.projectPoints(). Then we calculate the absolute norm between what we got with our transformation and the corner finding algorithm. To find the average error we calculate the arithmetical mean of the errors calculate for all the calibration images.
mean_error = 0
for i in xrange(len(objpoints)):
    imgpoints2, _ = cv2.projectPoints(objpoints[i], rvecs[i], tvecs[i], mtx, dist)
    error = cv2.norm(imgpoints[i],imgpoints2, cv2.NORM_L2)/len(imgpoints2)
    tot_error += error

print "total error: ", mean_error/len(objpoints)




Ref : http://docs.opencv.org/2.4/modules/calib3d/doc/camera_calibration_and_3d_reconstruction.html

https://cg2010studio.com/2012/01/04/opencv-%E7%9B%B8%E6%A9%9F%E7%9F%AF%E6%AD%A3-camera-calibration/

2016年6月1日 星期三

fastlab on 3D scan.... open source code

code ref : https://github.com/mariolukas/FabScanPi-Server/
R
The fastest way to use the FabScanPi is to use the FabScanPi Raspbian ISO image. 
Connect all needed hardware parts, download and flash the latest ISO image: 

Latest: fabscanpi_basic_jessie_v2016-05-15-14-05.zip

(655 MB) 
All list of all released images can be found here. 
Scince release v2016-05-15-14-05 the image files are zipped. You have to unzip the file to get the .img file. 
After flashing the image to an SD card point your browser to 
http://[ ip-address-of-your-raspberry-pi ]


More information about getting started with the software and hardware setup can be found here.






Ref : http://www.fabscan.org/

kinect on rpi




Ref : http://www.mariolukas.de/2015/04/proof-of-concept-3d-scanner-with-kinect-and-raspberry-pi2/

kinect linux device driver



https://github.com/xxorde/librekinect

神魔之塔的演算法研究...


Ref : http://35around.blogspot.tw/2014/07/ui-test-script-for-android-puzzle-and_14.html


接下來是最重要的演算法的部分, 如何設計一個好的演算法來求得傷害最高的路徑, 說真的, 我覺得要想到一個很好的演算法頗困難的, 還好目前的電腦cpu都還蠻強的, 加上其實盤面並不會很大, 只有6*5的大小, 所以在這邊我直接使用BFS方式來求解。實作的方式為: 針對盤面上的每一個珠子, 計算它可以往外前進的八個方向, 每走一步之後就計算目前盤面可以消除區域, 並計算落下之後是否造成下次的消除, 再代入上述的公式, 並且設定一個珠子最多可以走幾步, 最後就我們就可以得到每一個起始珠子走過的路徑的Combo數及攻擊力的權重表格, 透過這個表格我們就可以取到一個攻擊力最高的路徑了, 下面為pseudocode:
for row=1, 5 
    for col=1, 6
        solutions.add(new Solution(cursor(row, col));
while not solutions.isEmpty()
    solution = solutions.pop()
    if solution.depth >= MAX_DEPTH
        continue
    for move=1, 8
        newSolution = moveOrbs(solution, move)
        calculateWeight(newSolution) 
        solutions.add(newSolution) 


這個就是BFS 的演算法.......   這邊有提到是  Ref

https://github.com/kennytm/pndopt

接下來就是  code 的 study



row ->  y
col  ->  x


function evolve_solutions(solutions, weights, dir_step) {
    var new_solutions = [];
    solutions.forEach(function(s) {
        if (s.is_done) {
            return;
        }
        for (var dir = 0; dir < 8; dir += dir_step) {
            if (!can_move_orb_in_solution(s, dir)) {
                continue;
            }
            var solution = copy_solution(s);
            in_place_swap_orb_in_solution(solution, dir);
            in_place_evaluate_solution(solution, weights);
            new_solutions.push(solution);
        }
        s.is_done = true;
    });
    solutions = solutions.concat(new_solutions);
    solutions.sort(function(a, b) { return b.weight - a.weight; });
    return solutions.slice(0, MAX_SOLUTIONS_COUNT);
}

然後   can_move_orb_in_solution 

function can_move_orb_in_solution(solution, dir) {
    // Don't allow going back directly. It's pointless.  把直接退後的solution remove
    if (solution.path[solution.path.length-1] == (dir + 4) % 8) {
        return false;
    }
    return can_move_orb(solution.cursor, dir);
}

然後  can_move_orb() 的定義如下 :

function can_move_orb(rc, dir) {
    switch (dir) {
        case 0: return                    rc.col < COLS-1;
        case 1: return rc.row < ROWS-1 && rc.col < COLS-1;
        case 2: return rc.row < ROWS-1;
        case 3: return rc.row < ROWS-1 && rc.col > 0;
        case 4: return                    rc.col > 0;
        case 5: return rc.row > 0      && rc.col > 0;
        case 6: return rc.row > 0;
        case 7: return rc.row > 0      && rc.col < COLS-1;
    }
    return false;
}

col ->  y
row -> x

0  -> top ,  1 -> right and top ,  2 -> right,  3 -> right and down,  4 -> down
5  ->  left and down,   6 -> left ,     7  ->
           
 7    0    1
 6    p    2
 5    4    3

再來看   in_place_swap_orb_in_solution()

function in_place_swap_orb_in_solution(solution, dir) {
    var res = in_place_swap_orb(solution.board, solution.cursor, dir);
    solution.cursor = res.rc;
    solution.path.push(dir);
}

再來看     in_place_swap_orb()

function in_place_swap_orb(board, rc, dir) {
    var old_rc = copy_rc(rc);
    in_place_move_rc(rc, dir);
    var orig_type = board[old_rc.row][old_rc.col];
    board[old_rc.row][old_rc.col] = board[rc.row][rc.col];
    board[rc.row][rc.col] = orig_type;
    return {board: board, rc: rc};
}
 

function find_matches(board) {
    var match_board = create_empty_board();

    // 1. filter all 3+ consecutives.
    //  (a) horizontals
    for (var i = 0; i < ROWS; ++ i) {
        var prev_1_orb = 'X';
        var prev_2_orb = 'X';
        for (var j = 0; j < COLS; ++ j) {
            var cur_orb = board[i][j];
            if (prev_1_orb == prev_2_orb && prev_2_orb == cur_orb && cur_orb != 'X') {
                match_board[i][j] = cur_orb;
                match_board[i][j-1] = cur_orb;
                match_board[i][j-2] = cur_orb;
            }
            prev_1_orb = prev_2_orb;
            prev_2_orb = cur_orb;
        }
    }
    //  (b) verticals
    for (var j = 0; j < COLS; ++ j) {
        var prev_1_orb = 'X';
        var prev_2_orb = 'X';
        for (var i = 0; i < ROWS; ++ i) {
            var cur_orb = board[i][j];
            if (prev_1_orb == prev_2_orb && prev_2_orb == cur_orb && cur_orb != 'X') {
                match_board[i][j] = cur_orb;
                match_board[i-1][j] = cur_orb;
                match_board[i-2][j] = cur_orb;
            }
            prev_1_orb = prev_2_orb;
            prev_2_orb = cur_orb;
        }
    }

    var scratch_board = copy_board(match_board);

    // 2. enumerate the matches by flood-fill.
    var matches = [];
    for (var i = 0; i < ROWS; ++ i) {
        for (var j = 0; j < COLS; ++ j) {
            var cur_orb = scratch_board[i][j];
            if (typeof(cur_orb) == 'undefined') { continue; }
            var stack = [make_rc(i, j)];
            var count = 0;
            while (stack.length) {
                var n = stack.pop();
                if (scratch_board[n.row][n.col] != cur_orb) { continue; }
                ++ count;
                scratch_board[n.row][n.col] = undefined;
                if (n.row > 0) { stack.push(make_rc(n.row-1, n.col)); }
                if (n.row < ROWS-1) { stack.push(make_rc(n.row+1, n.col)); }
                if (n.col > 0) { stack.push(make_rc(n.row, n.col-1)); }
                if (n.col < COLS-1) { stack.push(make_rc(n.row, n.col+1)); }
            }
            matches.push(make_match(cur_orb, count));
        }
    }

    return {matches: matches, board: match_board};
}

function make_rc(row, col) {
    return {row: row, col: col};
}


function in_place_evaluate_solution(solution, weights) {
    var current_board = copy_board(solution.board);
    var all_matches = [];
    while (true) {
        var matches = find_matches(current_board);
        if (matches.matches.length == 0) {
            break;
        }
        in_place_remove_matches(current_board, matches.board);
        in_place_drop_empty_spaces(current_board);
        all_matches = all_matches.concat(matches.matches);
    }
    solution.weight = compute_weight(all_matches, weights);
    solution.matches = all_matches;
    return current_board;
}


//    rc ->  代表目前的   (x,y)   ,  dir -> 要移動的方向....  這個function 計算出下個點的col,row
function in_place_move_rc(rc, dir) {
    switch (dir) {
        case 0:              rc.col += 1; break;
        case 1: rc.row += 1; rc.col += 1; break;
        case 2: rc.row += 1;              break;
        case 3: rc.row += 1; rc.col -= 1; break;
        case 4:              rc.col -= 1; break;
        case 5: rc.row -= 1; rc.col -= 1; break;
        case 6: rc.row -= 1;              break;
        case 7: rc.row -= 1; rc.col += 1; break;
    }
}


//   update  掉下來的 code
function in_place_drop_empty_spaces(board) {
    for (var j = 0; j < COLS; ++ j) {
        var dest_i = ROWS-1;
        for (var src_i = ROWS-1; src_i >= 0; -- src_i) {
            if (board[src_i][j] != 'X') {
                board[dest_i][j] = board[src_i][j];
                -- dest_i;
            }
        }
        for (; dest_i >= 0; -- dest_i) {
            board[dest_i][j] = 'X';
        }
    }
    return board;
}


舉例  

src       -dest
0 O2
1 O1
2 O0      O2
3 X        O1
4 X        O0

src 從 4 開始  (row-1)
一開始  board[3][0]  為  X,      if (board[src_i][j] != 'X') 不成立

src 變 3
            board[3][0]  為  X,        if (board[src_i][j] != 'X') 不成立

src 變 2
            board[2][0]  為  O0,        if (board[src_i][j] != 'X') 成立

  dest 為 row            
board[dest-i][0] = O0  

依序下去....  就會造成消去的現象


function in_place_remove_matches(board, match_board) {
    for (var i = 0; i < ROWS; ++ i) {
        for (var j = 0; j < COLS; ++ j) {
            if (typeof(match_board[i][j]) != 'undefined') {
                board[i][j] = 'X';  ///  代表可以消去的
            }
        }
    }
    return board;
}


function compute_weight(matches, weights) {
    var total_weight = 0;
    matches.forEach(function(m) {
        var base_weight = weights[m.type][m.count >= 5 ? 'mass' : 'normal'];
        var multi_orb_bonus = (m.count - 3) * MULTI_ORB_BONUS + 1;
        total_weight += multi_orb_bonus * base_weight;
    });
    var combo_bonus = (matches.length - 1) * COMBO_BONUS + 1;
    return total_weight * combo_bonus;
}

 Normal (3+)
 Mass (5+)


當按下  rand 的 button

$('#randomize').click(function() {
        var types = $('#randomization-type').val().split(/,/);
        $('#grid > div').each(function() {
            var index = Math.floor(Math.random() * types.length); //
            show_element_type($(this), types[index]);
        });
        clear_canvas();
    });

 舉例來說....  當你選擇    5 + Heal  ,  randomization-type ->  "0,1,2,3,4,5"  ....

types.length ->  6

<select id="randomization-type">
        <option value="0,1,2">3-color</option>
        <option value="0,1,2,5">3 + Heal</option>
        <option value="0,1,2,3,4">5-color</option>
        <option value="0,1,2,3,4,5" selected="selected">5 + Heal</option>
        <option value="0,1,2,3,4,5,6">All</option>
    </select>

當按下  clear 的按鈕.... 把每個方格填入  'X'  ->  也就是問號....
 
 $('#clear').click(function() {
        $('#grid > div').each(function() { show_element_type($(this), 'X'); });
        clear_canvas();
    });

當按下   drop  按鈕 ....
$('#drop').click(function() {
        var solution = global_solutions[global_index];
        if (!solution) {
            return;
        }
        var board = in_place_evaluate_solution(solution, get_weights());
        show_board(board);
        clear_canvas();
    });


從   global_solution  抓出 solution......  


當按下   solve ...

 $('#solve').click(function() {
        var solver_button = this;
        var board = get_board();
        global_board = board;
        solver_button.disabled = true;
        solve_board(board, function(p, max_p) {
            $('#status').text('Solving (' + p + '/' + max_p + ')...');
        }, function(solutions) {
            var html_array = [];
            solutions = simplify_solutions(solutions);
            global_solutions = solutions;
            solutions.forEach(function(solution) {
                add_solution_as_li(html_array, solution, board);
            });
            $('#solutions > ol').html(html_array.join(''));
            solver_button.disabled = false;
        });
    });


//   給定新的dir...  然後和原本的互換....
function in_place_swap_orb(board, rc, dir) {
    var old_rc = copy_rc(rc);
    in_place_move_rc(rc, dir);
    var orig_type = board[old_rc.row][old_rc.col];
    board[old_rc.row][old_rc.col] = board[rc.row][rc.col];
    board[rc.row][rc.col] = orig_type;
    return {board: board, rc: rc};
}

//    呼叫剛剛   in_place_swap_orb   ...  然後把  dir push 到  path queue裡面....
function in_place_swap_orb_in_solution(solution, dir) {
    var res = in_place_swap_orb(solution.board, solution.cursor, dir);
    solution.cursor = res.rc;
    solution.path.push(dir);
}

function in_place_evaluate_solution(solution, weights) {
    var current_board = copy_board(solution.board);
    var all_matches = [];
    while (true) {
        var matches = find_matches(current_board);
        if (matches.matches.length == 0) {
            break;
        }
        in_place_remove_matches(current_board, matches.board);
        in_place_drop_empty_spaces(current_board);
        all_matches = all_matches.concat(matches.matches);
    }
    solution.weight = compute_weight(all_matches, weights);
    solution.matches = all_matches;
    return current_board;
}


//   遞迴....自己 call 自己
function solve_board_step(solve_state) {
    if (solve_state.p >= solve_state.max_length) {   //  結束條件
        solve_state.finish_callback(solve_state.solutions);
        return;
    }

    ++ solve_state.p;
    solve_state.solutions = evolve_solutions(solve_state.solutions,
                                             solve_state.weights,
                                             solve_state.dir_step);
    solve_state.step_callback(solve_state.p, solve_state.max_length);

    setTimeout(function() { solve_board_step(solve_state); }, 0);
}


//   感覺就是8各方向都去計算出來 score.....然後利用排序找出最佳的方向...
//   然後只留下那個方向....

function evolve_solutions(solutions, weights, dir_step) {
    var new_solutions = [];
    solutions.forEach(function(s) {
        if (s.is_done) {
            return;
        }
        for (var dir = 0; dir < 8; dir += dir_step) {
            if (!can_move_orb_in_solution(s, dir)) {
                continue;
            }
            var solution = copy_solution(s);
            in_place_swap_orb_in_solution(solution, dir);
            in_place_evaluate_solution(solution, weights);
            new_solutions.push(solution);
        }
        s.is_done = true;
    });
    solutions = solutions.concat(new_solutions);
    solutions.sort(function(a, b) { return b.weight - a.weight; });
    return solutions.slice(0, MAX_SOLUTIONS_COUNT);
}