<jupyter_start><jupyter_text>Augmented Reality This notebook aims to give a brief explaination of the theory and math behind the present *augmented reality* project.  The content is based on the articles ([part1](https://bitesofcode.wordpress.com/2017/09/12/augmented-reality-with-python-and-opencv-part-1/) and [part2](https://bitesofcode.wordpress.com/2018/09/16/augmented-reality-with-python-and-opencv-part-2/)) written by *<NAME>*, while the code that you can find in this project is a fork from his [project](https://github.com/juangallostra/augmented-reality) too.<jupyter_code># import all the required module for the notebook
import cv2
import numpy as np<jupyter_output><empty_output><jupyter_text>Other interesting articles:  [glyph recognition](https://rdmilligan.wordpress.com/2015/07/19/glyph-recognition-using-opencv-and-python/)  [AR with python OpenCV and OpenGL](https://rdmilligan.wordpress.com/2015/10/15/augmented-reality-using-opencv-opengl-and-blender/)  [pygame OBJ loader](http://www.pygame.org/wiki/OBJFileLoader)  [3D object search](https://clara.io/library) ---  IntroductionThe aim of the project is to display in a screen a 3D model of a figure whose position and orientation matches the position and orientation of some predefined flat surface. Furthermore, we want to do it in real time, so that if the surface changes its position or orientation the projected model does so accordingly.  To achieve this we first have to be able to identify the flat surface of reference in an image or video frame. Once identified, we can easily determine the transformation from the reference surface image (2D) to the target image (2D). This transformation is called [homography](https://en.wikipedia.org/wiki/Homography_(computer_vision)). However, if what we want is to project a 3D model placed on top of the reference surface to the target image we need to extend the previous transformation to handle cases were the height of the point to project in the reference surface coordinate system is different than zero. This can be achieved with a bit of algebra. Finally, we should apply this transformation to our 3D model and draw it on the screen. Bearing the previous points in mind our project can be divided into:1.  [Recognize the reference flat surface.](1.-Recognition-of-the-reference-target)2.  [Estimate the homography.](2.-Homography-estimation)3.  [Derive from the homography the transformation (pose estimation) from the reference surface coordinate system to the target image coordinate system.](3.-Pose-estimation-from-a-plane)4.  [Project our 3D model in the image (pixel space) and draw it.](4.-Model-projection)          The main tools we will use are *Python* and *OpenCV* because they are both open source, easy to set up and use and it is fast to build prototypes with them. For the needed algebra bit I will be using *numpy*. **insert image of 1.2.3.4. scheme here**  1. Recognition of the reference targetFrom the many possible techniques that exist to perform object recognition I decided to tackle the problem with a feature based recognition method. This kind of methods, without going into much detail, consist in three main steps: 1. [feature detection or extraction.](1.1.-Feature-extraction)2. [feature description.](1.2.-Feature-description)3. [feature matching.](1.3.-Feature-matching)  1.1. Feature extractionRoughly speaking, this step consists in first looking in both the reference and target images for features that stand out and, in some way, describe part the object to be recognized. This features can be later used to find the reference object in the target image. We will assume we have found the object when a certain number of positive feature matches are found between the target and reference images. For this to work it is important to have a reference image where the only thing seen is the object (or surface, in this case) to be found.  We don’t want to detect features that are not part of the surface. And, although we will deal with this later, we will use the dimensions of the reference image when estimating the pose of the surface in a scene.  For a region or point of an image to be labeled as feature it should fulfill two important properties: first of all, it should present some uniqueness at least locally. Good examples of this could be corners or edges. Secondly, since we don’t know beforehand which will be, for example, the orientation, scale or brightness conditions of this same object in the image where we want to recognize it a feature should, ideally, be invariant to transformations; i.e, invariant against scale, rotation or brightness changes. As a rule of thumb, the more invariant the better.  *Example of features extracted (right) from a reference target (left).*  1.2. [Feature description](https://docs.opencv.org/3.0-beta/doc/py_tutorials/py_feature2d/py_orb/py_orb.html)Once features have been found we should find a suitable representation of the information they provide. This will allow us to look for them in other images and also to obtain a measure of how similar two detected features are when being compared. This is were descriptors roll in. A descriptor provides a representation of the information given by a feature and its surroundings. Once the descriptors have been computed the object to be recognized can then be abstracted to a [feature vector](https://en.wikipedia.org/wiki/Feature_(machine_learning)),  which is a vector that contains the descriptors of the keypoints found in the image with the reference object.  This is for sure a nice idea, but how can it actually be done? There are many algorithms that extract image features and compute its descriptors and, since I won’t go into much more detail (a whole post could be devoted only to this) if you are interested in knowing more take a look at [SIFT](http://aishack.in/tutorials/sift-scale-invariant-feature-transform-introduction/), [SURF](http://www.vision.ee.ethz.ch/~surf/eccv06.pdf), or [Harris](http://aishack.in/tutorials/harris-corner-detector/). The one we will be using was developed at the OpenCV Lab and it is called [ORB](https://docs.opencv.org/3.0-beta/doc/py_tutorials/py_feature2d/py_orb/py_orb.html) (Oriented FAST and Rotated BRIEF). The shape and values of the descriptor depend on the algorithm used and, in our case,  the descriptors obtained will be binary strings.  With *OpenCV*, extracting features and its descriptors via the ORB detector is as easy as:<jupyter_code>img = cv2.imread('../reference/glyph_01.png', 0)

# Initiate ORB detector
orb = cv2.ORB_create()

# find the keypoints with ORB
kp = orb.detect(img, None)

# compute the descriptors with ORB
kp, des = orb.compute(img, kp)

# draw keypoints location, size and orientation
img2 = cv2.drawKeypoints(img, kp, img, color=(0,255,0), flags=4)

# resize image
width  = int(img2.shape[1] * 0.8)
height = int(img2.shape[0] * 0.8)
img3 = cv2.resize(img2, (width, height))

# display image
cv2.imshow('keypoints', img3)
cv2.waitKey(0)
cv2.destroyAllWindows()<jupyter_output><empty_output><jupyter_text>1.3. [Feature matching](https://docs.opencv.org/3.0-beta/doc/py_tutorials/py_feature2d/py_matcher/py_matcher.html)Once we have found the features of both the object and the scene were the object is to be found and computed its descriptors it is time to look for matches between them. The simplest way of doing this is to take the descriptor of each feature in the first set, compute the distance to all the descriptors in the second set and return the closest one as the best match (I should state here that it is important to choose a way of measuring distances suitable with the descriptors being used. Since our descriptors will be binary strings we will use [Hamming distance](https://en.wikipedia.org/wiki/Hamming_distance)). This is a brute force approach, and more sophisticated methods exist.For example, and this is what we will be also using, we could check that the match found as explained before is also the best match when computing matches the other way around, from features in the second set to features in the first set. This means that both features match each other. Once the matching has finished in both directions we will take as valid matches only the ones that fulfilled the previous condition.Another option to reduce the number of false positives would be to check if the distance to the second to best match is below a certain threshold.  If it is, then the match is considered valid.**image**Finally, after matches have been found, we should define some criteria to decide if the object has been found or not. For this I defined a threshold on the minimum number of matches that should be found. If the number of matches is above the threshold, then we assume the object has been found. Otherwise we consider that there isn’t enough evidence to say that the recognition was successful.With OpenCV all this recognition process can be done in a few lines of code:<jupyter_code>MIN_MATCHES = 30

# read images in grayscale mode
cap   = cv2.imread('images/scene_ticket.jpg', 0)    
model = cv2.imread('../reference/ticket.jpg', 0)

# ORB keypoint detector
orb = cv2.ORB_create()              
# create brute force  matcher object
bf = cv2.BFMatcher(cv2.NORM_HAMMING, crossCheck=True)  

# Compute model keypoints and its descriptors
kp_model, des_model = orb.detectAndCompute(model, None)  
# Compute scene keypoints and its descriptors
kp_frame, des_frame = orb.detectAndCompute(cap, None)

# Match frame descriptors with model descriptors
matches = bf.match(des_model, des_frame)
# Sort them in the order of their distance
matches = sorted(matches, key=lambda x: x.distance)

if len(matches) > MIN_MATCHES:
    # draw first 15 matches.
    cap = cv2.drawMatches(model, kp_model, cap, kp_frame,
                          matches[:MIN_MATCHES], 0, flags=2)
    # resize
    width  = int(cap.shape[1] * 0.3)
    height = int(cap.shape[0] * 0.3)
    cap2 = cv2.resize(cap, (width, height))
    # show result
    cv2.imshow('frame', cap2)
    cv2.waitKey(0)
    cv2.destroyAllWindows()
else:
    print("Not enough matches found - %d/%d" % (len(matches), MIN_MATCHES))<jupyter_output><empty_output><jupyter_text>On a final note and before stepping into the next step of the process I must point out that, since we want a real time application, it would have been better to implement a tracking technique and not just plain recognition. This is due to the fact that object recognition will be performed in each frame independently without taking into account previous frames that could add valuable information about the location of the reference object. Another thing to take into account is that, the easier to found the reference surface the more robust detection will be. In this particular sense, the reference surface I’m using might not be the best option, but it helps to understand the process.  2. [Homography estimation](https://docs.opencv.org/3.0-beta/doc/py_tutorials/py_feature2d/py_feature_homography/py_feature_homography.html)Once we have identified the reference surface in the current frame and have a set of valid matches we can proceed to estimate the homography between both images. As explained before, we want to find the transformation that maps points from the target surface plane to the image plane. This transformation will have to be updated each new frame we process.  Since we have already found a set of matches between both images we can certainly find directly by any of the existing methods (I advance we will be using [RANSAC](https://en.wikipedia.org/wiki/Random_sample_consensus)) an homogeneous transformation that performs the mapping.What we have is an object (a plane in this case) with known coordinates in the *World* coordinate system and we take a picture of it with a camera located at a certain position and orientation with respect to the World coordinate system. We will assume the camera works following the [pinhole model](https://en.wikipedia.org/wiki/Pinhole_camera_model), which roughly means that the rays passing through a 3D point *p* and the corresponding 2D point *u* intersect at *c*, the camera center. A good resource if you are interested in knowing more about the pinhole model can be found [here](http://alumni.media.mit.edu/~maov/classes/comp_photo_vision08f/lect/09_image_formation.pdf).  Although not entirely true, the pinhole model assumption eases our calculations and works well enough for our purposes. The *u*, *v* coordinates (coordinates in the image plane) of a point *p* expressed in the Camera coordinate system if we assume a pinhole camera can be computed as (the derivation of the equation is left as an exercise to the reader):$$  \begin{equation}    \begin{bmatrix}        u \cdot k \\        v \cdot k \\        k    \end{bmatrix}    =    A    \begin{bmatrix}        x_{cam} \\        y_{cam} \\        z_{cam}    \end{bmatrix}    \label{eq:pin_hole_model}    \tag{1}\end{equation}$$ where the matrix $$\begin{equation}    A =    \begin{bmatrix}        f_u & 0   & u_0\\        0   & f_v & v_0\\        0   & 0   & 1    \end{bmatrix}\end{equation}$$is know as *calibration matrix*, $(u_0, v_0)$ are the *projection of the optical center*, which is the position of the optical center in the image plane, $k$ is a scaling factor and $f_u, f_v$ are the *focal lengths*, the distance from the pinhole to the image plane.Equation ($\ref{eq:pin_hole_model}$) tells us how the image is formed. However, as stated before, we know the coordinates of the point *p* in the World coordinate system and not in the Camera coordinate system, so we have to add another transformation that maps points from the World coordinate system to the Camera coordinate system. The transformation that tells us the coordinates in the image plane of a point *p* in the World coordinate system is then:  $$  \begin{equation}    \begin{bmatrix}        u \cdot k \\        v \cdot k \\        k    \end{bmatrix}    =    A     \begin{bmatrix}        R & t    \end{bmatrix}    \begin{bmatrix}        x_{world} \\        y_{world} \\        z_{world} \\        1    \end{bmatrix}    \label{eq:pin_hole_model_extended}    \tag{2}\end{equation}$$ Where $R$ is a $3 \times 3$ rotation matrix and $t$ is the translation vector that specifies the position of the reference target.Luckily for us, since the points in the reference surface plane do always have its *z* coordinate equal to 0 (it is a surface!) we can simplify the transformation that we found above. It can be easily seen that the product of the *z* coordinate and the third column of the projection matrix will always be $0$ so we can drop this column and the *z* coordinate from the previous equation. $$  \begin{equation}    \begin{bmatrix}        u \cdot k \\        v \cdot k \\        k    \end{bmatrix}    =    A     \begin{bmatrix}        R_1 & R_2 & R_3 & t    \end{bmatrix}    \begin{bmatrix}        x_{world} \\        y_{world} \\        0 \\        1    \end{bmatrix}    =    A    \begin{bmatrix}        R_1 & R_2 & t    \end{bmatrix}    \begin{bmatrix}        x_{world} \\        y_{world} \\        1    \end{bmatrix}    =     H    \begin{bmatrix}        x_{world} \\        y_{world} \\        1    \end{bmatrix}    \label{eq:pin_hole_model_simplified}    \tag{3}\end{equation}$$ There are several methods that allow us to estimate the values of the homography matrix $H$, and you might be familiar with some of them. The one we will be using is **RANdom SAmple Consensus** ([RANSAC](https://en.wikipedia.org/wiki/Random_sample_consensus)).  RANSAC is an iterative algorithm used for model fitting in the presence of a large number of outliers and can be schematize as follows.1. Choose a small subset of points uniformly at random.2. Fit a model to that subset.3. Find all remaining points that are "close" to the model and reject the rest as outliers.4. Do this many times and choose the best model.Since we cannot guarantee that all the matches we have found are actually valid matches we have to consider that there might be some false matches (which will be our outliers) and, hence, we have to use an estimation method that is robust against outliers.For homography estimation the algorithm can be outlined with the following steps:1. Randomly sample 4 matches.2. Estimate Homography H.3. Verify homography: search for other matches consistent with H.4. Iterate until convergence.Before seeing how OpenCV can handle this for us we should  discuss one final aspect of the algorithm, which is what does it mean that a match is consistent with $H$. What this mainly means is that if after estimating an homography we project into the target image the matches that were not used to estimate it then the projected points from the reference surface should be close to its matches in the target image. How close they should be to be considered consistent is up to you.In OpenCV estimating the homography with RANSAC is as easy as:<jupyter_code># assuming matches stores the matches found and 
# returned by bf.match(des_model, des_frame)
# differenciate between source points and destination points
src_pts = np.float32([kp_model[m.queryIdx].pt for m in matches]).reshape(-1, 1, 2)
dst_pts = np.float32([kp_frame[m.trainIdx].pt for m in matches]).reshape(-1, 1, 2)
# compute Homography
M, mask = cv2.findHomography(src_pts, dst_pts, cv2.RANSAC, 5.0)<jupyter_output><empty_output><jupyter_text>Where 5.0 is the threshold distance to determine if a match is consistent with the estimated homography. If after estimating the homography we project the four corners of the reference surface on the target image and connect them with a line we should expect the resulting lines to enclose the reference surface in the target image. We can do this by:<jupyter_code>model = cv2.imread('../reference/ticket.jpg', 0)
frame = cv2.imread('images/scene_ticket.jpg', 0)    
# Draw a rectangle that marks the found model in the frame
h, w = model.shape
pts = np.float32([[0, 0], [0, h - 1], 
                  [w - 1, h - 1], [w - 1, 0]]).reshape(-1, 1, 2)
# project corners into frame
dst = cv2.perspectiveTransform(pts, M)  
# connect them with lines
img = cv2.polylines(frame, [np.int32(dst)], True, 255, 3, cv2.LINE_AA)


img2 = cv2.drawMatches(model, kp_model, img, kp_frame, matches[:MIN_MATCHES], 0, flags=2)
# resize
width  = int(img2.shape[1] * 0.4)
height = int(img2.shape[0] * 0.4)
img2 = cv2.resize(img2, (width, height))
# show
cv2.imshow('frame', img2)
cv2.waitKey(0)
cv2.destroyAllWindows()<jupyter_output><empty_output><jupyter_text>3. Pose estimation from a planeTo locate the 3D model in the camera frame we need its pose estimation with respect the camera frame, that is a rotation matrix $R$ and a position vector $t$, that is a homogeneus transformation from points in the *World* frame to the *Camera* frame.  This values are the same of equation ($\ref{eq:pin_hole_model_simplified}$), and since at this point the homography matrix $H$ has been already estimated we could invert equation ($\ref{eq:pin_hole_model_simplified}$) and obtain: $$  \begin{equation}    \begin{bmatrix}        R_1 & R_2 & t    \end{bmatrix}    =     A^{-1}H    \label{eq:homogeneus_transformation}    \tag{4}\end{equation}$$ Equation ($\ref{eq:homogeneus_transformation}$) has still some problem.  First of all matrix $A$ is unknown, there are methods to evaluate an estimation (see [this](https://opencv-python-tutroals.readthedocs.io/en/latest/py_tutorials/py_calib3d/py_calibration/py_calibration.html), [this](https://www.learnopencv.com/approximate-focal-length-for-webcams-and-cell-phone-cameras/) and [this](http://ksimek.github.io/2012/08/13/introduction/)), but for now the code assume that the matrix is known *a priori* with some "rational" values.  Another problem is the absence of $R_3$, which is the estmate of the *z-axis* of the *World* frame. So it is necessary to add a step to evalluate $R_3$. It can be evaluated as $R_1 \times R_2$, however the values of $R_1$ and $R_2$ are affected by some estimation error and they can not result to be *orthonormal*, so its needed to correct them too.  *Remark: the correction of $R_1$ and $R_2$ is not important to the estimate of $R_3$ but for the representation of the entire matrix $R$! if $R_1$, $R_2$ and $R_3$ are not orthonormal the transformation $R$ rotates and **deforms** the points transformated resulting in the rotation and **deformation** of the 3D object rendered!*The estimation of the basis considers the assumption that, since $R_1$ and $R_2$ are estimates of the real ${R_1}'$ and ${R_2}'$ (which are orthonormal), the angle between $R_1$ and $R_2$ will be not 90 degrees. Furthermore, the modulus of each of this vectors will be close to 1. From $R_1$ and $R_2$ we can easily compute an orthogonal basis -this meaning that the angle between the basis vectors will exactly be 90 degrees- that will be rotated approximately 45° degrees clockwise with respect to the basis formed by $R_1$ and $R_2$. This basis is the one formed by $c=R_1+R_2$ and  $d=c \times p=(R_1+R_2)\times (R_1\times R_2)$ in figure below. If the vectors that form this new basis $(c,d)$ are made unit vectors and rotated 45° degrees counterclockwise (note that once the vectors have been transformed into unit vectors rotating the basis is as easy as d’ = c / ||c|| + d / ||d|| and  c’ = c / ||c|| – d / ||d||), guess what? We will have an orthogonal basis which is pretty close to our original basis $(R_1,R_2)$. If we normalize this rotated basis we will finally get the pair of vectors we were looking for.  Once this basis $({R_1}',{R_2}')$ has been obtained it is trivial to get the value of $R_3$ as $R_1' \times R_2'$.  Below you can find the function which compute the estimate.<jupyter_code>def projection_matrix(camera_parameters, homography):
"""
 From the camera calibration matrix and the estimated homography
 compute the 3D projection matrix
 """
# Compute rotation along the x and y axis as well as the translation
homography = homography * (-1)
rot_and_transl = np.dot(np.linalg.inv(camera_parameters), homography)
col_1 = rot_and_transl[:, 0]
col_2 = rot_and_transl[:, 1]
col_3 = rot_and_transl[:, 2]
# normalise vectors
l = math.sqrt(np.linalg.norm(col_1, 2) * np.linalg.norm(col_2, 2))
rot_1 = col_1 / l
rot_2 = col_2 / l
translation = col_3 / l
# compute the orthonormal basis
c = rot_1 + rot_2
p = np.cross(rot_1, rot_2)
d = np.cross(c, p)
rot_1 = np.dot(c / np.linalg.norm(c, 2) + d / np.linalg.norm(d, 2), 1 / math.sqrt(2))
rot_2 = np.dot(c / np.linalg.norm(c, 2) - d / np.linalg.norm(d, 2), 1 / math.sqrt(2))
rot_3 = np.cross(rot_1, rot_2)
# finally, compute the 3D projection matrix from the model to the current frame
projection = np.stack((rot_1, rot_2, rot_3, translation)).T
return np.dot(camera_parameters, projection)<jupyter_output><empty_output><jupyter_text>4. Model projectionThe program currently only uses simple models in `.obj` format, because they are easy to process and render directly with bare Python without having to make use of other libraries such as OpenGL. The problem with complex models is that the amount of processing they require is way more than what my computer can handle. Since I want my application to be real-time, this limits the complexity of the models I am able to render.  The code I used to load the models is based on this [OBJFileLoader](http://www.pygame.org/wiki/OBJFileLoader) script that I found on Pygame’s website. I stripped out any references to OpenGL and left only the code that loads the geometry of the model. Once the model is loaded we just have to implement a function that reads this data and projects it on top of the video frame with the projection matrix we obtained in the previous section. To do so we take every point used to define the model and multiply it by the projection matrix. One this has been done, we only have to fill with color the faces of the model. The following function can be used to do so:<jupyter_code>def render(img, obj, projection, model, color=False):
    vertices = obj.vertices
    scale_matrix = np.eye(3) * 3
    h, w = model.shape

    for face in obj.faces:
        face_vertices = face[0]
        points = np.array([vertices[vertex - 1] for vertex in face_vertices])
        points = np.dot(points, scale_matrix)
        # render model in the middle of the reference surface. To do so,
        # model points must be displaced
        points = np.array([[p[0] + w / 2, p[1] + h / 2, p[2]] for p in points])
        dst = cv2.perspectiveTransform(points.reshape(-1, 1, 3), projection)
        imgpts = np.int32(dst)
        if color is False:
            cv2.fillConvexPoly(img, imgpts, (137, 27, 211))
        else:
            color = hex_to_rgb(face[-1])
            color = color[::-1] # reverse
            cv2.fillConvexPoly(img, imgpts, color)

    return img<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.07892409134770023, 'fraction_numerical': 0.014956577677709875, 'mean_word_length': 4.0637214983713354, 'pattern_counts': {'":': 0, '<': 30, '<?xml version=': 0, '>': 31, 'https://': 18, 'lorem ipsum': 0, 'www.': 4, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '7', 'dir': 'jupyter-structured-clean-dedup', 'id': '217772', 'license': 'MIT', 'n_tokens_mistral': 7130, 'n_tokens_neox': 6794, 'n_words': 3337, 'path': 'notebook/Augmented_Reality.ipynb', 'repo_name': 'IIRAI/augmented-reality', 'size': '24880'}