HaLongBay Rock 3D-Reconstruction

3D-Reconstruction of one of the rocks in HaLongBay, Vietnam, from a series of drive-by images.

Initialization of the Environment


In [1]:
# Magic!
%pylab inline


Populating the interactive namespace from numpy and matplotlib

In [2]:
# Imports
import os

import numpy as np
import scipy as sp

import matplotlib.pylab as plt
import matplotlib.image as mpimg
from matplotlib import cm

from PIL import Image

import cv2

Dataset: The Rock

A Series of 6 images, taken with a Canon SX 280 point-and-shoot camera, full resolution (12.1mpix?), saved as JPG.

Load Images


In [3]:
# setup location and shape of the dataset
img_num = 6

# cwd is at 'base/ipnbs'
base_dir = os.path.normpath("..")
img_dir = os.path.join(base_dir,os.path.normpath("data/halong_rock/"))

img_fnames = [ os.path.join(img_dir,"{}.jpg".format(i)) for i in range(img_num) ]

In [4]:
# finally load the images with matplotlib
imgs = [ mpimg.imread(img_fname) for img_fname in img_fnames ]

This is the first image of the series as example:


In [5]:
plt.imshow(imgs[0]);


Also convert the images to grayscale for later use:


In [6]:
imgs_gray = [ np.asarray(Image.fromarray(img).convert('L')) for img in imgs ]
# OpenCV alternative
#imgs_gray = [ cv2.cvtColor(img,cv2.COLOR_RGB2GRAY) for img in imgs ]

In [7]:
plt.imshow(imgs_gray[0],cmap = cm.Greys_r);


Finding Image Features

Now detect some features to find matches/tracks along the image series. Try different feature detection algorithms and compare reliability, coverage, speed and what ever else may come to mind, say quality of the abbreviation name or googleability.

First, create some containers (dicts) to store the data categorized by feature algorithm. These will be dicts using the algorithms abbreviation as key, e.g. 'FAST' oder 'ORB'. For per-image data, they will contain a list for each key, i.e. keypoints['FAST'] is a list of keypoint-lists.


In [8]:
feature_detectors = {}
keypoint_sets = {}

FAST


In [9]:
feature_detectors['FAST'] = cv2.FastFeatureDetector()
fast = feature_detectors['FAST']
keypoint_sets['FAST'] = [ fast.detect(gray,None) for gray in imgs_gray ]
# do runtime measurement later grouped with other feature detection algorithms
#%timeit keypoints_list = [ fast.detect(gray,None) for gray in imgs_gray ]

In [10]:
print 'Number of features:\t\t', len(keypoint_sets['FAST'][0])
print 'Position of first feature:\t', keypoint_sets['FAST'][0][0].pt
print 'Angle of first feature:\t\t', keypoint_sets['FAST'][0][0].angle
print 'Size of first feature:\t\t', keypoint_sets['FAST'][0][0].size


Number of features:		259365
Position of first feature:	(2371.0, 497.0)
Angle of first feature:		-1.0
Size of first feature:		7.0

ORB


In [11]:
feature_detectors['ORB'] = cv2.ORB()
orb = feature_detectors['ORB']
keypoint_sets['ORB'] = [ orb.detect(gray,None) for gray in imgs_gray ]
# do runtime measurement later grouped with other feature detection algorithms
#%timeit keypoints_list = [ orb.detect(gray,None) for gray in imgs_gray ]

Detected features:


In [12]:
print 'Number of features:\t\t', len(keypoint_sets['ORB'][0])
print 'Position of first feature:\t', keypoint_sets['ORB'][0][0].pt
print 'Angle of first feature:\t\t', keypoint_sets['ORB'][0][0].angle
print 'Size of first feature:\t\t', keypoint_sets['ORB'][0][0].size


Number of features:		500
Position of first feature:	(2350.0, 779.0)
Angle of first feature:		8.09868907928
Size of first feature:		31.0

Display the Features

scale the features - this is for visualization reasons only!


In [57]:
# scale img to change line-thickness, thickness cannot be changed in OpenCV directly (sigh)
img_scale = 0.12
# scale keypoints size to bring them to a reasonable size
# Note: the final size is also affected by 'img_scale'
kp_scales = {}
kp_scales['FAST'] = 0.01
kp_scales['ORB'] = 0.5

feature_names = ['FAST','ORB']
# use first image to display features
disp_img_idx = 0

keypoints_scaled = {}
for feat_name in feature_names:
    keypoints = keypoint_sets[feat_name][disp_img_idx]
    keypoints_scaled[feat_name] = [ cv2.KeyPoint(kp.pt[0],kp.pt[1],kp.size,kp.angle)
                                        for kp in keypoints ]
    for keypoint in keypoints_scaled[feat_name]:
        keypoint.size *= kp_scales[feat_name]
        keypoint.pt = tuple(img_scale * coord for coord in keypoint.pt)

#keypoints_scaled = [ cv2.KeyPoint(kp.pt[0],kp.pt[1],kp.size,kp.angle) for kp in keypoints_list[0] ]
#for keypoint in keypoints_scaled:
#    keypoint.size *= kp_scale
#    keypoint.pt = tuple(img_scale * coord for coord in keypoint.pt)

Prepare (scale) the display image


In [58]:
tmp_img = Image.fromarray(np.array(imgs_gray[0]))
scaled_size = tuple( s * img_scale for s in tmp_img.size)
tmp_img.thumbnail(scaled_size,Image.ANTIALIAS)
img_scaled = np.asarray(tmp_img)

Now draw the scaled features on the scaled display image


In [59]:
#img = imgs_gray[0]
#kps = keypoints_list[0]
display_img = img_scaled
kps = keypoints_scaled

imgs_keypoints = {}
flags = cv2.DRAW_MATCHES_FLAGS_DRAW_RICH_KEYPOINTS
for feat_name in feature_names:
    imgs_keypoints[feat_name] = cv2.drawKeypoints(display_img,kps[feat_name],flags = flags)
    
#plt.imshow(img_keypoints[750:800,2300:2400,:]);
#plt.imshow(img_keypoints[400:1200,1800:2400,:]);
#plt.figure(num=1,figsize=(12,14))
#plt.imshow(imgs_keypoints['FAST']);

fig = plt.figure(figsize=(15, 5))

N = len(feature_names)
for i,feat_name in enumerate(feature_names):
    # max 3 plots per row, number of columns smaller than max if N < max
    plt.subplot(N/3+1,min(N%3+1,N),i+1)
    im = plt.imshow(imgs_keypoints[feat_name])
    im.axes.set_xlabel(feat_name)


Matching the Features


In [ ]:


In [ ]:


In [ ]: