Programming

Simple and effective coin segmentation using Python and OpenCV

The new generation of OpenCV bindings for Python are getting better and better with the hard work of the community. The new bindings, called “cv2″ are the replacement of the old “cv” bindings; in this new generation of bindings, almost all operations returns now native Python objects or Numpy objects, which is pretty nice since it simplified a lot and also improved performance on some areas due to the fact that you can now also use the optimized operations from Numpy and also enabled the integration with other frameworks like the scikit-image which also uses Numpy arrays for image representation.

In this example, I’ll show how to segment coins present in images or even real-time video capture with a simple approach using thresholding, morphological operators and contour approximation. This approach is a lot simpler than the approach using Otsu’s thresholding and Watershed segmentation here in OpenCV Python tutorials, which I highly recommend you to read due to its robustness. Unfortunatelly, the approach using Otsu’s thresholding is highly dependent on a illumination normalization. One could extract small patches of the image to implement something similar to an adaptive Otsu’s binarization (like the one implemented in Letptonica – the framework used by Tesseract OCR) to overcome this problem, but let’s see another approach. For reference, see the output of the Otsu’s thresholding using an image taken with my webcam with a non-normalized illumination:

 

Original image vs Otsu binarization

Original image vs Otsu binarization

1. Setting the Video Capture configuration

The first step to create a real-time Video Capture using the Python bindings is to instantiate the VideoCapture class, set the properties and then start reading frames from the camera:

1
2
3
4
5
6
import numpy as np
import cv2

cap = cv2.VideoCapture(0)
cap.set(cv2.cv.CV_CAP_PROP_FRAME_WIDTH, 1280)
cap.set(cv2.cv.CV_CAP_PROP_FRAME_HEIGHT, 720)

In newer versions (unreleased yet), the constants for CV_CAP_PROP_FRAME_WIDTH are now in the cv2 module, for now let’s just use the cv2.cv module.

2. Reading image frames

The next step is to use the VideoCapture object to read the frames and then convert them to gray color (we are not going to use color information to segment the coins):

1
2
3
4
while True:
    ret, frame = cap.read()
    roi = frame[0:500, 0:500]
    gray = cv2.cvtColor(roi, cv2.COLOR_BGR2GRAY)

Note that here I’m extracting a smal portion of the complete image (where the coins are located), but you don’t have to do that if you have only coins on your image. At this moment, we have the following gray image:

The original Gray image captured.

The original Gray image captured.

3. Applying adaptive thresholding

In this step we will apply the Adaptive Thresholding after applying a Gaussian Blur kernel to eliminate the noise that we have in the image:

1
2
3
gray_blur = cv2.GaussianBlur(gray, (15, 15), 0)
thresh = cv2.adaptiveThreshold(gray_blur, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C,
    cv2.THRESH_BINARY_INV, 11, 1)

See the effect of the Gaussian Kernel in the image:

The original gray image and the image after applying the Gaussian Kernel.

The original gray image and the image after applying the Gaussian Kernel.

And now the effect of the Adaptive Thresholding with the blurry image:

The effect of the adaptive thresholding into the blurry image

Note that at that moment we already have the coins segmented except for the small noisy inside the center of the coins and also in some places around them.

4. Morphology

The Morphological Operators are used to dilate, erode and other operations on the pixels of the image. Here, due to the fact that sometimes the camera can present some artifacts we will use the Morphological Operation of Closing to make sure that the borders of the coins are always close, otherwise we may found a coin with a semi circle or something like that. To understand the effect of the Closing operation (which is the operation of erosion of the pixels already dilated) see the image below:

Morphological Closing

You can see that after some iterations of the operation, the circles starts to become filled. To use the Closing operation, we’ll use the morphologyEx function from the OpenCV Python bindings:

1
2
kernel = np.ones((3, 3), np.uint8)
closing = cv2.morphologyEx(thresh, cv2.MORPH_CLOSE, kernel, iterations=4)

See now the effect of the Closing operation on our coins:

Closing Operation in the coins

The operations of Morphological Operators are very simple, the main principle is the application of a element (in our case we have a block element of 3×3) into the pixels of the image. If you want to understand it, please see this animation explaning the operation of Erosion.

5. Contour detection and filtering

After applying the morphological operators, all we have to do is to find the contour of each coin and then filter the contours having an area smaller or larger than a coin area. You can imagine the procedure of finding contours in OpenCV as the operation of finding connected components and their boundaries. To do that, we’ll use the OpenCV findContours function.

1
2
3
cont_img = closing.copy()
contours, hierarchy = cv2.findContours(cont_img, cv2.RETR_EXTERNAL,
    cv2.CHAIN_APPROX_SIMPLE)

Note that we made a copy of the closing image because the function findContours will change the image passed as the first parameter, we’re also using the RETR_EXTERNAL flag, which means that the contours returned are only the extreme outer contours. The parameter CHAIN_APPROX_SIMPLE will also return a compact representation of the contour, for more information see here.

After finding the contours, we need to iterate into each one and check the area of them to filter the contours containing an area greater or smaller than the area of a coin. We also need to fit an ellipse to the contour found. We could have done this using the minimum enclosing circle, but since my cameara isn’t perfectly above the coins, the coins appear with a small inclination describing an ellipse.

1
2
3
4
5
6
7
8
9
10
for cnt in contours:
    area = cv2.contourArea(cnt)
    if area < 2000 or area > 4000:
        continue

    if len(cnt) < 5:
        continue

    ellipse = cv2.fitEllipse(cnt)
    cv2.ellipse(roi, ellipse, (0,255,0), 2)

Note that in the code above we are iterating on each contour, filtering coins with area smaller than 2000 or greater than 4000 (these are hardcoded values I found for the brazilian coins at this distance from the camera), later we check for the number of points of the contour because the function fitEllipse needs a number of points greater or equal than 5 and finally we use the ellipse function to draw the ellipse in green over the original image.

To show the final image with the contours we just use the imshow function to show a new window with the image:

1
cv2.imshow('final result', roi)

And finally, this is the result in the end of all steps described above:

The final image with the contours detected

The final image with the contours detected

The complete source-code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
import numpy as np
import cv2

def run_main():
    cap = cv2.VideoCapture(0)
    cap.set(cv2.cv.CV_CAP_PROP_FRAME_WIDTH, 1280)
    cap.set(cv2.cv.CV_CAP_PROP_FRAME_HEIGHT, 720)

    while(True):
        ret, frame = cap.read()
        roi = frame[0:500, 0:500]
        gray = cv2.cvtColor(roi, cv2.COLOR_BGR2GRAY)

        gray_blur = cv2.GaussianBlur(gray, (15, 15), 0)
        thresh = cv2.adaptiveThreshold(gray_blur, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C,
            cv2.THRESH_BINARY_INV, 11, 1)

        kernel = np.ones((3, 3), np.uint8)
        closing = cv2.morphologyEx(thresh, cv2.MORPH_CLOSE,
            kernel, iterations=4)

        cont_img = closing.copy()
        contours, hierarchy = cv2.findContours(cont_img, cv2.RETR_EXTERNAL,
            cv2.CHAIN_APPROX_SIMPLE)
       
        for cnt in contours:
            area = cv2.contourArea(cnt)
            if area < 2000 or area > 4000:
                continue

            if len(cnt) < 5:
                continue

            ellipse = cv2.fitEllipse(cnt)
            cv2.ellipse(roi, ellipse, (0,255,0), 2)

        cv2.imshow("Morphological Closing", closing)
        cv2.imshow("Adaptive Thresholding", thresh)
        cv2.imshow('Contours', roi)

        if cv2.waitKey(1) & 0xFF == ord('q'):
            break

    cap.release()
    cv2.destroyAllWindows()

if __name__ == "__main__":
    run_main()

The beauty of Bitcoin P2P network

So, in the last days I just released Protocoin, a framework in pure Python with a Bitcoin P2P network implementation. While I’m in process of development of the v.0.2 of the framework (with new and nice features like Bitcoin keys management – you can see some preview here) I would like to show a real-time visualization I’ve made with Protocoin and Ubigraph of a node connecting to a seed node and then issuing GetAddr message for each node and connecting on the received nodes in a breadth-first search fashion. I’ll release the code used to create this visualization in the next release of Protocoin as soon as possible. I hope you enjoy it !

Color legend

Yellow = Connecting
Green = Connected
Blue = Disconnected after connection

Video

Protocoin – a pure Python Bitcoin protocol implementation

Just release the first version of Protocoin, a pure Python Bitcoin protocol implementation, for more information see the documentation or the project in Github.

If you want to donate for the project, my Bitcoin address is: 1Q6JZEE5turJXaTn1fburWkqjhC4oMJ4yV.

I hope you like it !

- Christian S. Perone

Mapa de calor dos dados de acidentes de transito do DataPoa

Esta semana será disponibilzada a nova versão do Django GIS Brasil, segue abaixo um exemplo de mapa criado usando os novos dados do Django GIS Brasil importados do DataPoa.

O exemplo abaixo é um mapa de calor utilizando os dados de acidentes de trânsito em Porto Alegre /RS durante os anos de 2000 até 2012. Os eixos (ruas, avenidas, etc.) também estarão presentes no Django GIS Brasil.

 

Heatmap POA

Mapa de Acidentes de Trânsito em PoA/RS. Clique para ampliar.

Machine Learning :: Cosine Similarity for Vector Space Models (Part III)

* It has been a long time since I wrote the TF-IDF tutorial (Part I and Part II) and as I promissed, here is the continuation of the tutorial. Unfortunately I had no time to fix the previous tutorials for the newer versions of the scikit-learn (sklearn) package nor to answer all the questions, but I hope to do that in a close future.

So, on the previous tutorials we learned how a document can be modeled in the Vector Space, how the TF-IDF transformation works and how the TF-IDF is calculated, now what we are going to learn is how to use a well-known similarity measure (Cosine Similarity) to calculate the similarity between different documents.

The Dot Product

Let’s begin with the definition of the dot product for two vectors: \vec{a} = (a_1, a_2, a_3, \ldots) and \vec{b} = (b_1, b_2, b_3, \ldots), where a_n and b_n are the components of the vector (features of the document, or TF-IDF values for each word of the document in our example) and the \mathit{n} is the dimension of the vectors:

  \vec{a} \cdot \vec{b} = \sum_{i=1}^n a_ib_i = a_1b_1 + a_2b_2 + \cdots + a_nb_n

As you can see, the definition of the dot product is a simple multiplication of each component from the both vectors added together. See an example of a dot product for two vectors with 2 dimensions each (2D):

  \vec{a} = (0, 3) \\   \vec{b} = (4, 0) \\   \vec{a} \cdot \vec{b} = 0*4 + 3*0 = 0

The first thing you probably noticed is that the result of a dot product between two vectors isn’t another vector but a single value, a scalar.

This is all very simple and easy to understand, but what is a dot product ? What is the intuitive idea behind it ? What does it mean to have a dot product of zero ? To understand it, we need to understand what is the geometric definition of the dot product:

  \vec{a} \cdot \vec{b} = \|\vec{a}\|\|\vec{b}\|\cos{\theta}

Rearranging the equation to understand it better using the commutative property, we have:

   \vec{a} \cdot \vec{b} = \|\vec{b}\|\|\vec{a}\|\cos{\theta}

So, what is the term \displaystyle \|\vec{a}\|\cos{\theta} ? This term is the projection of the vector \vec{a} into the vector \vec{b} as shown on the image below:

The projection of the vector A into the vector B. By Wikipedia.

Now, what happens when the vector \vec{a} is orthogonal (with an angle of 90 degrees) to the vector \vec{b} like on the image below ?

Two orthogonal vectors (with 90 degrees angle).

There will be no adjacent side on the triangle, it will be equivalent to zero, the term \displaystyle \|\vec{a}\|\cos{\theta} will be zero and the resulting multiplication with the magnitude of the vector \vec{b} will also be zero. Now you know that, when the dot product between two different vectors is zero, they are orthogonal to each other (they have an angle of 90 degrees), this is a very neat way to check the orthogonality of different vectors. It is also important to note that we are using 2D examples, but the most amazing fact about it is that we can also calculate angles and similarity between vectors in higher dimensional spaces, and that is why math let us see far than the obvious even when we can’t visualize or imagine what is the angle between two vectors with twelve dimensions for instance.

The Cosine Similarity

The cosine similarity between two vectors (or two documents on the Vector Space) is a measure that calculates the cosine of the angle between them. This metric is a measurement of orientation and not magnitude, it can be seen as a comparison between documents on a normalized space because we’re not taking into the consideration only the magnitude of each word count (tf-idf) of each document, but the angle between the documents. What we have to do to build the cosine similarity equation is to solve the equation of the dot product for the \cos{\theta}:

   \displaystyle   \vec{a} \cdot \vec{b} = \|\vec{a}\|\|\vec{b}\|\cos{\theta} \\   \cos{\theta} = \frac{\vec{a} \cdot \vec{b}}{\|\vec{a}\|\|\vec{b}\|}

And that is it, this is the cosine similarity formula. Cosine Similarity will generate a metric that says how related are two documents by looking at the angle instead of magnitude, like in the examples below:

The Cosine Similarity values for different documents, 1 (same direction), 0 (90 deg.), -1 (opposite directions).

Note that even if we had a vector pointing to a point far from another vector, they still could have an small angle and that is the central point on the use of Cosine Similarity, the measurement tends to ignore the higher term count on documents. Suppose we have a document with the word “sky” appearing 200 times and another document with the word “sky” appearing 50, the Euclidean distance between them will be higher but the angle will still be small because they are pointing to the same direction, which is what matters when we are comparing documents.

Now that we have a Vector Space Model of documents (like on the image below) modeled as vectors (with TF-IDF counts) and also have a formula to calculate the similarity between different documents in this space, let’s see now how we do it in practice using scikit-learn (sklearn).

Vector Space Model

Practice Using Scikit-learn (sklearn)

* In this tutorial I’m using the Python 2.7.5 and Scikit-learn 0.14.1.

The first thing we need to do is to define our set of example documents:

1
2
3
4
5
6
documents = (
"The sky is blue",
"The sun is bright",
"The sun in the sky is bright",
"We can see the shining sun, the bright sun"
)

And then we instantiate the Sklearn TF-IDF Vectorizer and transform our documents into the TF-IDF matrix:

1
2
3
4
5
from sklearn.feature_extraction.text import TfidfVectorizer
tfidf_vectorizer = TfidfVectorizer()
tfidf_matrix = tfidf_vectorizer.fit_transform(documents)
print tfidf_matrix.shape
(4, 11)

Now we have the TF-IDF matrix (tfidf_matrix) for each document (the number of rows of the matrix) with 11 tf-idf terms (the number of columns from the matrix), we can calculate the Cosine Similarity between the first document (“The sky is blue”) with each of the other documents of the set:

1
2
3
from sklearn.metrics.pairwise import cosine_similarity
cosine_similarity(tfidf_matrix[0:1], tfidf_matrix)
array([[ 1.        ,  0.36651513,  0.52305744,  0.13448867]])

The tfidf_matrix[0:1] is the Scipy operation to get the first row of the sparse matrix and the resulting array is the Cosine Similarity between the first document with all documents in the set. Note that the first value of the array is 1.0 because it is the Cosine Similarity between the first document with itself. Also note that due to the presence of similar words on the third document (“The sun in the sky is bright”), it achieved a better score.

If you want, you can also solve the Cosine Similarity for the angle between vectors:

   \cos{\theta} = \frac{\vec{a} \cdot \vec{b}}{\|\vec{a}\|\|\vec{b}\|}

We only need to isolate the angle (\theta) and move the \cos to the right hand of the equation:

   \theta = \arccos{\frac{\vec{a} \cdot \vec{b}}{\|\vec{a}\|\|\vec{b}\|}}

The \arccos is the same as the inverse of the cosine (\cos^-1).

 Lets for instance, check the angle between the first and third documents:
1
2
3
4
5
6
import math
# This was already calculated on the previous step, so we just use the value
cos_sim = 0.52305744
angle_in_radians = math.acos(cos_sim)
print math.degrees(angle_in_radians)
58.462437107432784
 And that angle of ~58.5 is the angle between the first and the third document of our document set.
That is it, I hope you liked this third tutorial !

Related Material

A video about Dot Product on The Khan Academy

Wikipedia: Dot Product

Wikipedia: Cosine Similarity

Scikit-learn (sklearn) – The de facto Machine Learning package for Python

Accessing HP Cloud OpenStack Nova using Python and Requests

So, my request to enter on the free and private beta season of the new HP Cloud Services was gently accepted by the HP Cloud team, and today I finally got some time to play with the OpenStack API at HP Cloud. I’ll start with the first impressions I had with the service:

The user interface of the management is very user-friendly, the design is much like of the Twitter Bootstrap, see the screenshot below of the “Compute” page from the “Manage” section:

As you can see, they have a set of 4 Ubuntu images and a CentOS, I think that since they are still in the beta period, soon we’ll have more default images to use.

Here is a screenshot of the instance size set:

Since they are using OpenStack, I really think that they should have imported the vocabulary of the OpenStack into the user interface, and instead of calling it “Size”, it would be more sensible to use “Flavour“.

The user interface still doesn’t have many features, something that I would really like to have is a “Stop” or something like that for the instances, only the “Terminate” function is present on the Manage interface, but those are details that they should be still working on since they’re only in beta.

Another important info to cite is that the access to the instances are done through SSH using a generated RSA key that they provide to you.

Let’s dig into the OpenStack API now.

OpenStack API

To access the OpenStack API you’ll need the credentials for the authentication, HP Cloud services provide these keys on the Manage interface for each zone/service you have, see the screenshot below (with keys anonymized of course):

Now, OpenStack authentication could be done in different schemes, the scheme that I know that HP supports is the token authentication. I know that there is a lot of clients already supporting the OpenStack API (some have no documentation, some have weird API design, etc.), but the aim of this post is to show how easy would be to create a simple interface to access the OpenStack API using Python and Requests (HTTP for Humans !).

Let’s start defining our authentication scheme by sub-classing Requests AuthBase:

1
2
3
4
5
6
7
8
9
class OpenStackAuth(AuthBase):
    def __init__(self, auth_user, auth_key):
        self.auth_key = auth_key
        self.auth_user = auth_user

def __call__(self, r):
    r.headers['X-Auth-User'] = self.auth_user
    r.headers['X-Auth-Key'] = self.auth_key
    return r

As you can see, we’re defining the X-Auth-User and the X-Auth-Key in the header of the request with the parameters. These parameters are respectively your Account ID and  Access Key we cited earlier. Now, all you have to do is to make the request itself using the authentication scheme, which is pretty easy using Requests:

1
2
3
4
ENDPOINT_URL = 'https://az-1.region-a.geo-1.compute.hpcloudsvc.com/v1.1/'
ACCESS_KEY = 'Your Access Key'
ACCOUNT_ID = 'Your Account ID'
response = requests.get(ENDPOINT_URL, auth=OpenStackAuth(ACCOUNT_ID, ACCESS_KEY))

And that is it, you’re done with the authentication mechanism using just a few lines of code, and this is how the request is going to be sent to the HP Cloud service server:

 This request is sent to the HP Cloud Endpoint URL (https://az-1.region-a.geo-1.compute.hpcloudsvc.com/v1.1/). Let’s see now how the server answered this authentication request:

You can show this authentication response using Requests by printing the header attribute of the request Response object. You can see that the server answered our request with two important header items: X-Server-Management-URL and the X-Auth-Token. The management URL is now our new endpoint, is the URL we should use to do further requests to the HP Cloud services and the X-Auth-Token is the authentication Token that the server generated based on our credentials, these tokens are usually valid for 24 hours, although I haven’t tested it.

What we need to do now is to sub-class the Requests AuthBase class again but this time defining only the authentication token that we need to use on each new request we’re going to make to the management URL:

1
2
3
4
5
6
7
class OpenStackAuthToken(AuthBase):
    def __init__(self, request):
        self.auth_token = request.headers['x-auth-token']

def __call__(self, r):
    r.headers['X-Auth-Token'] = self.auth_token
    return r

Note that the OpenStackAuthToken is receiving now a response request as parameter, copying the X-Auth-Token and setting it on the request.

Let’s consume a service from the OpenStack API v.1.1, I’m going to call the List Servers API function, parse the results using JSON and then show the results on the screen:

1
2
3
4
5
6
7
8
9
10
# Get the management URL from the response header
mgmt_url = response.headers['x-server-management-url']

# Create a new request to the management URL using the /servers path
# and the OpenStackAuthToken scheme we created
r_server = requests.get(mgmt_url + '/servers', auth=OpenStackAuthToken(response))

# Parse the response and show it to the screen
json_parse = json.loads(r_server.text)
print json.dumps(json_parse, indent=4)

And this is what we get in response to this request:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
{
    "servers": [
        {
            "id": 22378,
            "uuid": "e2964d51-fe98-48f3-9428-f3083aa0318e",
            "links": [
                {
                    "href": "https://az-1.region-a.geo-1.compute.hpcloudsvc.com/v1.1/20817201684751/servers/22378",
                    "rel": "self"
                },
                {
                    "href": "https://az-1.region-a.geo-1.compute.hpcloudsvc.com/20817201684751/servers/22378",
                    "rel": "bookmark"
                }
            ],
            "name": "Server 22378"
        },
        {
            "id": 11921,
            "uuid": "312ff473-3d5d-433e-b7ee-e46e4efa0e5e",
            "links": [
                {
                    "href": "https://az-1.region-a.geo-1.compute.hpcloudsvc.com/v1.1/20817201684751/servers/11921",
                    "rel": "self"
                },
                {
                    "href": "https://az-1.region-a.geo-1.compute.hpcloudsvc.com/20817201684751/servers/11921",
                    "rel": "bookmark"
                }
            ],
            "name": "Server 11921"
        }
    ]
}

And that is it, now you know how to use Requests and Python to consume OpenStack API. If you wish to read more information about the API and how does it works, you can read the documentation here.

- Christian S. Perone

C++0x :: Introduction to some amazing features

I’ve made those slides in March of this year to a training session, I was expecting to get time to update it to cover more features but I wasn’t able to do that yet, so I’m publishing it here for those who are interested in some of the new language features.