Example of extracting grid nodes in MATLAB

For last few posts, I am writing about morphological image analysis with python with open source package scikit-image. In this post I want to show benefit of morphological image analysis with some real application. I been working with road extraction from satellite images for past few weeks and  one of the common task is to extract the grid corner nodes of the road network. So let’s assume the road grid has been already extracted and we want to extract intersection of grid network. Here I am showing three approaches for that based on synthetic data. Obviously with real data, the task is more complex and might require many slight modifications.

Here is the grid. You can generate such grid with checkerboard function in MATLAB with minor processing (edge detection and hole filling). Intersections of of grid lines have a property that it has 4 white pixels in 4-N neighborhood.  We can exploit this property to extract it. I am going to show to ways by which it can be done and in addition, another method based on extraction of corners will also be shown.

%make a checkerboard
I=checkerboard(20, 10,10)>0.5;
imshow(I)
% detect edge with canny
BW = edge(I,'canny');
figure, imshow(BW)
%make square structural element to fill holes with closing
SE = strel('square', 3);
BW1 = imclose(BW, SE);
figure, imshow(BW1)
% make a skeleton of the to lines 1 pixel thick
BW2 = bwmorph(BW1,'skel',Inf);
figure, imshow(BW2)

 


In mathematics and, in particular, functional analysis, convolution is a mathematical operation on two functions f and g, producing a third function that is typically viewed as a modified version of one of the original function .In simple word, convolution of image is a process where you scan your image from left to right and top to bottom within local neighborhood defined by user, and do mathematical calculation within the local neighborhood. Here we are going to use it for simply sum the number of pixels which are which are white within 3x3 neighborhood. As shown by above figure, if the point is an intersection grid point then at that point sum should be 5. Any points with less than 5 white pixels are not intersection grid points.

% with block processing and inline function with convolution
tic;
kernel = [0 1 0; ... %# Convolution kernel
1 1 1; ...
0 1 0];
sumX = conv2(double(BW2),kernel,'same');
result=sumX;
% only consider pixels which are in grid image
result (BW2==0)= 0
result(sumX<5)=0;
[r,c] = find(result>0);
toc;
figure,imshow((BW2)), hold on, plot(r,c,'r*');title ('with Convolution')
hold off

Way Two: With Mathematical Morphology (MM)

If you are working in the field of Geo-information and you still saying what the heck this ‘MM’ then seriously you should be acquainted with this field of image processing. You can use MM for classification, edge detection, filtering, segmentation, building detection and many other things. Grab a cup of coffee and Google MM in remote sensing; there is tons of stuff to read.  There are courses in MM taught in different universities which are one semester long, so you got an idea how vast the subject is. If you want to seriously know MM, then there is no better book than ‘ Morphological Image Analysis’ by P Soillie. Here we use ‘Erosion’ method of MM.  Time: 0.011184 seconds.

% method two with morphological analysis
tic;
SE = strel ('disk', 1);
result2 = imerode(BW, SE);
[r,c] = find(result>0);
toc;
figure,imshow((BW2)), hold on, plot(r,c,'b*'), title ('With MM Erosion')
hold off

Way Two:  Mathematical Morphology (MM)

There are many corner detection techniques such as Moravec, Harris, SIFT etc. Here, I am going to use Harris corner detector for intersection grid points. Time: 0.094848 seconds.

% detect corner by harris corner detector
tic;
C = corner(BW2,500 );
timet2 = toc;
toc;
figure,imshow(BW2);
hold on
plot(C(:,1), C(:,2), 'g*'); title('with harris corner detector')
hold off

Way Two:  With Harris Corner Detector

So you have seen there are number of ways to solve a particular problem. All above methods detected intersection points successfully. Among three methods, convolution required longest time for processing and MM required shortest time. The gain in computation cost between MM and Corner detection is 9 folds.  I did all this in MATLAB but it can easily coded in python with scikit-image.

Example of Opening by reconstruction in scikit-image

This blog is a continuation of the last blog that I have written. The aim is to give you potential of the mathematical morphology (MM) using ski-image and their application.  MM is branch of image processing which has wide application in many diverse field. Erosion and dilation are basic operators of MM which can be combined in different way to form many different powerful MM operators. MM was devised for binary images but nowadays it has been extended to grayscale images as well.

Here I would show you application of erosion and opening by reconstruction to do some image processing. For definition of erosion and opening by reconstruction please follow any image processing books or look in the web there are many online tutorials.

So we have an image. Let us assume that out the different blobs that we have we want to have a blob which should have a minimum length of 200 pixels in horizontal direction. Such expert knowledge is handy in many images processing domain to improve the image analysis. E.g building cannot be less than certain width, road should be of certain width, bridge should be of certain width etc. in this case lets assume, the shown blobs are result of some building detection algorithm. Now we want to improve the result by incorporating expert knowledge using MM.

So first we want to erode the image with linear structural element of length 200 pixels and whatever is remaining after that, we want to recover original shape of the blobs which was distorted by using erosion. For that purpose we use opening by reconstruction. Here are the codes, enjoy and dont forget to play around with python.




# -*- coding: utf-8 -*-
Created on Mon Nov 18 15:42:21 2013
@author: shailesh
"""
#%%
% call necessary libraries
import math
import matplotlib.pyplot as plt
import numpy as np
from skimage.draw import ellipse
from skimage.draw import polygon
from skimage.draw import circle
from skimage.morphology import label
from skimage.measure import regionprops
from skimage.transform import rotate
import matplotlib.patches as mpatches
import skimage.morphology as MM
% draw some arbitary shapes
image = np.zeros((1000, 1000),dtype=uint8)
&nbsp;
% create an ellipse with centre (350,350) with minor axis = 100 and major = 220
rr, cc = ellipse(300, 750, 100, 220)
image[rr,cc] = 1
&nbsp;
&nbsp;
% create a polygon
x = np.array([1, 7, 4, 1])
x = np.array([1, 70, 40, 1])
y = np.array([1, 20, 80, 1])
rr, cc = polygon(y, x)
image[rr, cc] = 1
&nbsp;
% create a polygon
rr, cc = circle(200, 200, 50)
image[rr, cc] = 1
fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(6, 6))
&nbsp;
&nbsp;
&nbsp;
rr, cc = ellipse(500, 800, 30, 40)
image[rr, cc] = 1
ax.imshow(image)
&nbsp;
#lable connected regions
label_img = label(image)
#find properties of connected regions
regions = regionprops(label_img,['Centroid','Area','BoundingBox'])
&nbsp;
%loop through connected regions and find
% Centroid', 'BoundingBox', 'Area' for each blobs
i= 0
&nbsp;
NOS= len(np.unique(label_img))-1
X = np.zeros((1,NOS))
Y = np.zeros((1,NOS))
&nbsp;
for props in regions:
y0, x0 = props['Centroid']
print x0, y0
%show figure
plt.gray()
plt.axis((0, 1000, 1000, 0))
plt.title('original image')
plt.show()
&nbsp;
%perform erosion with linear structuring element
SE= np.ones((1,200))
figure()
imageE = MM.binary_erosion(image, SE)
plt.imshow(imageE)
plt.title('Eroded image')
plt.show()
&nbsp;
figure()
%perform opening by reconstruction to recover blob
imageRECON = MM.reconstruction(imageE, image)
plt.imshow(imageRECON)
plt.title('Reconstructed image')
plt.show()

Morphological reconstruction is a very powerful method. People have use it for many different applications such as watershed delieanation, filtering, change detection of building after earthquakes, building detection, bridge detection etc. For my work also, i have used the technique for unsupervised change detection using high resolution images and change detection of high rise buildings.

A little teaser for using scikit learn for classification of remotesensing images

I work with classification of remote sensing images a lot both with supervised as well as unsupervised classification. Unsupervised classifications don’t need any external input where as supervised classifications need samples or training areas for an algorithm to learn. For processing remote sensing images, there are many proprietary software like ENVI, ERDAS, PCI Geomatica, Global Mapper, eCognition and many more. Even after paying thousands of Euros, classifications algorithm available in Costs-off-the-self (COTS) software is far from satisfactory.

Software	Available Algorithms
eCognition 8.7 (cost >1000 Euros/year)	KNN Decision Tree DT SVM (no Y and C parameters available for tunning, no way to perform grid search for optimal determination of Y and C, Only linear and rbf kernel available) Random Forest  RF (only available in 8.8
ENVI (cost >1000 Euros/year)	Maximum Likelihood (ML) SVM (no way to perform grid search for optima, unless you get your hands dirty with IDL programming Neural Network NN ( Numbers of hidden nodes cannot be assigned) Some otheralgorithms more suitable for Hyperspectral imagesy SAM, SID, Spectral Unmixing.
Scikit learn Free Opensoure	Neighrest Neigbour (NN) Decision Tree (DT) SVM ( Grid search, cross validation flexibility) Random Forest RF AdaBoost Naives Bayes Linear Discriminant Analysis (LDA) Quadratic Discriminant Analysis (QDA)

Here is the little teaser of classification accuracy with many algorithms that are available in scikit-learn for a remote sensing imagery. In near future, I will blog with more illustration and with code. Till then go and make your hands dirty with Python and Scikit-Learn. Make that your new year resolution and trust me, you will thank me for that.



Here, algorithms hyperparameters were not optimally tuned hence superior machine learning algorithms like SVM has very low accuracy for test samples which are not seen by trained model.

Image processing in python with scikit-learn and pymorph

I have been using MATLAB quite a bit for some time. But being a proprietary language it has some drawbacks. One has to fork out lot of money to buy it and renew it yearly license. On top of that one hast to busy many toolboxes.Check this blog for finding various advantages of switching over from MATLAB to Python. You will find very helpful insights into key difference  between MATLAB and python numpy here.

Recently for couples of weeks now, I have to experimenting with python which is an open source software for some of my tasks that I do mainly related to classification, segmentation of remote sensing image etc during my weekend free time. So after sacrificing two of my weekends time for scouring python and different libraries which otherwise I would spend for my leisure activity, I am now pretty comfortable with python language as well for image processing tasks. Sweet :-) .

In this post I will show some morphological image analysis with python. I can see this will be series of post related to python. I plan to cover classification, clustering, edge detection in my subsequent posts.

For installing python, I have used pythonXY which comes with many libraries that might be needed for the scientific data analysis such as scipy, numpy, matplotlib and many others.PythoXY also comes with an Interactive development environment (IDE) for python which is called Spyder. With Spyder, you don’t have run python from command prompt. Other benefit includes integrated help, variable explorer, history console etc. In addition its all FREE.

Here is the interface how it looks.

I will be using library called “scikit-image” which comes in a bundle with pythonXY and in addition I would be using a library called “PyMorph” for morphological image analysis. Scikit-image also has a morphologicl module but there you would only find basic morphological operators like opening, closing, erosion and dialation. In PyMorph there are many advance morphological operator such opening by reconstruction and closing by reconstruction, Auto sequential filter(ASF), ASF by reconstruction etc. So after installing pythonXY, download the PYMoprh library and install it.

So, in this post I would show you how to use scikit-image for creating some basic shapes, calculating i area for each shapes or blobs and writing areas in the figure, finding bounding box of each blobs and plotting it in the figures. After finding areas of each figures, I will show how to remove certain blob  which has less area t than specified area.

% -*- coding: utf-8 -*-
Created on Mon Nov 18 15:42:21 2013
@author: shailesh
"""
%%%
% call necessary libraries
import math
import matplotlib.pyplot as plt
import numpy as np
from skimage.draw import ellipse
from skimage.draw import polygon
from skimage.draw import circle
from skimage.morphology import label
from skimage.measure import regionprops
from skimage.transform import rotate
import matplotlib.patches as mpatches
import pymorph as MM

%%%
% draw some arbitary shapes
image = np.zeros((1000, 1000),dtype=uint8)
% create an ellipse with centre (350,350) with minor axis = 100 and major = 220
rr, cc = ellipse(300, 350, 100, 220)
image[rr,cc] = 1

% create a polygon
x = np.array([1, 7, 4, 1])
x = np.array([1, 70, 40, 1])
y = np.array([1, 20, 80, 1])
rr, cc = polygon(y, x)
image[rr, cc] = 1

% create a polygon
rr, cc = circle(200, 200, 50)
image[rr, cc] = 1
fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(6, 6))
ax.imshow(image)

%lable connected regions
label_img = label(image)
%find properties of connected regions
regions = regionprops(label_img,['Centroid', 'BoundingBox', 'Area'])

%loop through connected regions and find
% Centroid', 'BoundingBox', 'Area' for each blobs
for props in regions:
y0, x0 = props['Centroid']
minr, minc, maxr, maxc = props['BoundingBox']
area =props['Area']
rect = mpatches.Rectangle((minc, minr), maxc - minc, maxr - minr,
fill=False, edgecolor='red', linewidth=2)
ax.add_patch(rect)
ax.text(x0, y0, str(area), fontsize=10, color="blue")
%show figure
plt.gray()
plt.axis((0, 1000, 1000, 0))
plt.show()

%%% remove small blobs
image = image.astype(np.uint8)
%remove blobs with area less than 10000 '
% MM.areaopen is coming from PyMorph
b=MM.areaopen(image,10000)
fig, ax = plt.subplots(ncols=1, nrows=1, figsize=(6, 6))
ax.imshow(b)
ax.title('blobs with areas greater than 10000 ')
plt.gray()
plt.axis((0, 1000, 1000, 0))
plt.show()
% to close the pixels
plt.close("all")

Here are the figures produced by above code.

Initial figure after some blobs were made. The area of each blobs and bounding box were plotted as well


Figure where blobs with area less than 10000 pixels are remove with pymorph.


So this is a hypothetical simple illustration. But the real application of this kind of analysis could be many in remote sensing field. For example, after detection buildings you can easily remove buildings with size less than minimum mapping unit (say 100 m2). Similarly it can be used for removing some artifacts which are small in area that you might not be interested in. And the list can go on.

In next few blogs, i will be writing more about other morphological image processing operators as well as classfiers such Support Vector Machine (SVM), Random Forest (RF) etc.

MATLAB linking figure for simultaneous zooming and panning

Sometime back I have written a Matlab GUI, which allows you select variables from your work space and plot them in linked figures so that when you pan and zoom, other figures zoom automatically. Lately I have received a lot of request to share that GUI with other but unfortunately I am unable to do so as I developed it for some research works and I am not allowed to distribute it. But, I am going to share you a function which does that linking.

To use that function, download it from here. There is a test.m for testing purpose as well.

%first image
Im1 = imread('cameraman.tif');
% show first figure
imshow(Im1)
% threshold image
Im2 = Im1>50;
%show second image in another figure
figure,
imshow(Im2);
% then call linkaxes_shailesh ()
% remember that matriz size must be same</code>

Put it in a safe place in a folder, set path in Matlab to that folder so that the function is available in every Matlab session ( file> set path). Open two or more figure you want to link in but bear in mind that the variables used in the figures must of same size i.e. it should have same number of rows and columns. After that, call that function by typing linkaxes_shailesh in your command line or you can even call it within one of your m file or function files by linkaxes_shailesh ()  and viola, your figures will be linked.  Its simple as that, don’t believe? Try out yourself. Someone said rightly ‘ Seeing is believing’. I hope you find it useful.

<iframe width="560" height="315" src="//www.youtube.com/embed/QhqnGdk2yQU?list=UUHQMpdmbYQPJxqsnKuz9DUw" frameborder="0" allowfullscreen></iframe>

Neighborhood analysis in QGIS.

QGIS is an open source software which anyone can download for free. The community is growing each day and it is completely capable of doing common basic GIS tasks as well as more advance stuff. In this tutorial, I will show you how to do some neighbourhood analysis in QGIS with the help of sextante graphical modeler. If you guys are more ArcGIS person then you might be familiar with ArcGIS modeler. QGIS provides similar graphical interface as ArcGIS through which all the available algorithms in QGIS are accessible. On top of that you can combine algorithms from different other providers such as GDAL, Grass, OrfeoToolbox, R and SagaGIS. So it can be really powerful tool to automate your daily GIS task. Go check it out. Scour help file of QGIS, you will find plenty of readings.

Different programs that are available with QGIS

So the task is: There is a polygon shape file of my country Nepal which is shown in green. Each polygon represents a district. There is a point shp file which represent point of interest (assume it to be major infrastructures such as big hospitals or airports or universities etc). Each point of interest is represented by unique ID. Now the task is to assign each district to closest point of interest based on how far the district is from point of interest

Input shpfiles: 

The graphical model in QGIS.

The model is basically doing following main task.

• ·         Convert shp file districts to its centroid
• ·         Calculate nearest point of interest from each district centroid.

After creating the model one could just hit the run button and you are good to go. The model can be saved and used it for later purpose. As you can see in the modeller, the first algorithm is coming from SAGA GIS whereas second algorithms is coming from core QGIS. This is also a good illustration to show how algorithms can be combined from different source.

After running the model, you will get a point file. Let’s examine its attribute table. As you can see, each district is associated with @HubName which in our case is point of interest. It also has distance information. I forgot to converts projection EPSG: 4632 to a projected one so the distance is in terms of latitude and longitude. If you want distance in m or km then use a projected EPSG and you will be fine.

Now we will join this information to the original district shp file based on the column @Name. To do that do the following:

Right click the polygon shp file and click properties. Click join tab and following figure pops up. Click green + button to define a join as shown in next figure. The join column is @Name which is common is both polygon shp files as well point shp files that we calculated from our graphical modeler

Join attribute interface

.

Defining join layer and columns

After the join, the attribute of districts polygon has additional information. Now it’s the matter of applying colour symbology to differentiate which districts are closer to which one of the point of interest based on fourth column.

Attribute table after join

The final result is :

I have not shown each and every step in this tutorial. My main purpose is to show how you can use powerful QGIS and sextante modeler to perform complex GIS task such as neighborhood analysis without using ArcGIS which costs a lot and stimulate you to try QGIS. If any questions, just ask and I will reply. Of course, if I am not busy with my shits J. 

Maximum reflectance in a spectra of Multispectral or Hyperspectral image in MATLAB

This is a very short post, something related to processing of remote sensing images. The image can be either multi-spectral or hyper-spectral images. One of my colleagues asked me this simple help in MATLAB  He was trying that in ENVI with band-math and IDL, with no avail.

I wanted to produce 2 images from my hyper-spectral image: 1) showing the maximum reflectance value across the bands and 2) the band number where the maximum came from.

First step is to read the image in MATLAB, of course :). That's no-brainer. For that you have to use multibandread function. I have a written post about it earlier. Search in the site. Now you have your image read in the MATLAB.

Here is the code, that i have written to perform what my friend asked.Hope this is useful for someone else.

%find the number of rows and columns

[rows,cols,bands]= size(imt1);

%reshape it 

imt1= reshape(imt1,[],4)';

%find max_value and band which has maximum value

[max_value,idx]= max(imt1,[],1);

%reconstruct the image

band_image= reshape(idx,rows,cols);

max_image= reshape(max_value,rows,cols);