Detecting Circular Objects using Template Matching

Motivation: The rationale behind the selection of this project was to explore the opportunity to learn new areas of skills in Detecting Circular Objects. This concept can be integrated to the ‘Automotive’ industry to detect circular traffic boards and perform task intelligently.

Problem Statement: The challenge of this project is to Detect Circular Objects from a range of objects. The detection of circular object can be done through several mechanisms. Some of them are already available as inbuilt function of MATLAB®. But some of them have their own disadvantages. As an example Hough transform works worst when there are large numbers of edges.

Solution at a glance: The proposed solution contains mainly two parts. First is to detect the edges which are appeared on the image. Second part which works on template matching contains two sub parts which are generated concurrently as firstly creating the template and secondly detecting the circles.

Method: In this project the edges are detected using a customized code where the Sobel Kernels are used in a moderated level. See [1]. The kernels are sequenced throughout the image and results are stored in matrix ‘result’.

 Image 

Why we ignore some terms in the equation is, our aim is to detect circular objects and not to detect all the edges in the image as in [1]. The produced output is then processed using template matching. In template matching every possible circles have to be matched with the image. For that possible circles are created using ‘meshgrid’ and updated as template image. (Circles also can be created using edge detection)

The created template images are then matched and calculated their dot products as in Cross correlation unlike in [3]. In this method when the template is placed at a particular position then each of the counterparts are multiplied and get the summation. The results are stored in ‘final’ matrix. When the template is highly correlated with the image at a particular pixel then the pixel will be brighter than others as shown in [2].

The next main stream was to define the threshold value to extract possible pixel points which have been highly correlated. The results are stored in ‘g’ matrix. One thing we have to be more precise is the deviation in the intensity levels when highly correlated and hence differ the threshold. For that an equation is generated as follows.

g = final > 250* (diameter+var) <–(1) 

Variable ‘var’ is to change the parameters of the threshold and the diameter is of the circle we are trying to match with. Other notations are as described above.

Then as to the highly correlated values on ‘g’, another walk through on the image will create possible circles around the objected circles.

Results:

 Image

            Image 1: Detecting circles of same size

The results are largely depended on the threshold value we defined to extract Cross Correlated results. As you can see when the brightness level, object’s size and clear stand out are differed across the image then the threshold levels are varied.

Image 1

(17 circles)

Image 2

(3 circles)

Image 3

(1 traffic board)

Image 4

(8 coins)

‘var’

# Detected circles

‘var’

# Detected circles

‘var’

Detected? (yes/no)

‘var’

# Detected coins

80

17

12

3

45

Yes

100

8 (2 non circle)

82

13

20

2

85

Yes

120

6 (1 non circle)

85

2

30

1

100

Yes(Not clear)

130

5 (1 non circle)

90

0

50

1

110

No

145

4 (1 non circle)

 Image

Image 2: Detecting circles of different sizes

When the ‘var’ value in (1) equation is incremented, circles that can be detected by the program is reduced. In image 2, the circle which is at the middle is not disappeared as ‘var’ increases as in image 1. This is because the contrast level at which this circle is differentiated is large.

 Image

 Image 3: Detecting traffic signs

In image 3, the detection of traffic sign is illustrated. In image 4 even there are only 6 out of 8 are circular (C) coins, since the non circular (NC) coins are round shape those are susceptible to be recognized as circular coins when ‘var’ is 100. But when ‘var’ is 120 the 6th coin (which is an NC) will not be detected. Since 3rd coin (which is an NC) is highly contrasted it is detected. But 7th C coin is not detected due to low contrast. To overcome this, the threshold level set at edge detection should be reduced. Another reason for this is the objects are close to each other and hence edges are not clear.

 Image

Image 4: Detecting coins of different sizes

Conclusion: The extraction of the circular object on the image has extensively used the edge detection technique which is on the concept of Sobel operation. But the Sobel operation has been used in this project in a more moderated level. Even the moderated level does not significantly affect the edge detection of circular objects, the enhancement in the speed of the operation has been enormously reduced.

The point we have to keep in mind here is Generally Circular object detection through Template matching is noisier than Circular Hough Transform. See [4]. An extraction on comparing these two methods has been provided in [4].

The main criticism in this program is the Threshold level we define on each image. This can be overcome by normalizing the edge detected image’s brightness and suitably defining a z value for the normal distribution. But in this project it was not expected and hence it was not examined. Therefore the Threshold level varies with picture to picture.

A clear example on this project can be visualized as follows.

 Image

Image 5: Detecting a traffic sign

References:

[1] http://angeljohnsy.blogspot.com/2011/12/sobel-edge-detection.html

[2] http://docs.opencv.org/doc/tutorials/imgproc/histograms/template_matching/template_matching.html

[3] http://en.wikipedia.org/wiki/Template_matching

[4]http://www.intechopen.com/books/international_journal_of_advanced_robotic_systems/bolting-cabin-assistance-system-using-a-sensor-network

Designing a Custom Constellation Modulation Scheme

We were asked to design a custom constellation modulation module in 5th semester as a partial fulfillment of the module Communication II. The implementation of the module should be in MATLAB. All the groups including our groups were given with a different constellation diagram than other groups. This is to individually evaluate each of the groups to ensure the students in each group have gained knowledge in modulation properly. The Constellation diagram we were given with was displayed below.

First of all the main step we had to take was giving appropriate value (bit pattern) to each given point in constellation. For that we used gray coding. In the given constellation diagram there were only 8 points. It means 3 symbols have to compare with a single point in constellation diagram. Gray code makes one symbol error to one bit error. One thing we have to keep in mind is the gray code won’t work as worked in PSKs and QAMs. Therefore we tried to implement the gray code starting from minimum distant points and go to higher distant points.

25

When googled the requirement of the project, the first useful post we encountered on this topic was: http://www.mathworks.com/matlabcentral/fileexchange/17263-developing-custom-modulation-schemes/content/publishedCustomModulationAgilentLogo/custommodulationagilentlogo.html

In this post the designing of a modulation scheme has been illustrated in a descriptive way using MATLAB. The thing we had to do from this point was just manipulating the code and modifying it. But the modifying part was not as easy as described above. We screwed up in some parts in the modulation design. This post is about how the designing of the modulation was carried out by using the above post. (Throughout the discussion the 8-PSK is referenced to the custom modulation scheme for better comparability)

1. Defining the parameters

%% Setup
% Define parameters.
M = 8; % Size of signal constellation
k = log2(M); % Number of bits per symbol
n = 3e4; % Number of bits to process
nSamp = 1; % Oversampling rate

A = 1;
N = 0.02;

%% Customizing constellation.
inphase = [A/2 0 A A*cos(pi/4) 0 -A/2 -A*cos(pi/4) -A];
quadr = [0 A/2 0 A*sin(pi/4) -A/2 0 -A*sin(pi/4) 0];
const = inphase + 1i*quadr;

% Create a scatter plot
scatterPlot = commscope.ScatterPlot(‘SamplesPerSymbol’,1,…
‘Constellation’,const);
% Show constellation
scatterPlot.PlotSettings.Constellation = ‘on’;
scatterPlot.PlotSettings.ConstellationStyle = ‘*’;
title(‘Customized Constellation for My-QAM’);

2. Design the modulator and demodulator

%% Create Modulator and Demodulator
% modulator and demodulator of custom constellation
hMod = modem.genqammod(‘Constellation’,const); % Create a My-QAM modulator
hDemod = modem.genqamdemod(hMod); % Create a My-QAM demodulator based on the modulator

% modulator and demodulator of 8-PSK constellation
hMod_ref = modem.pskmod(8); % Create a 8-PSK modulator
hDemod_ref = modem.pskdemod(8); % Create a 8-PSK demodulator based on the modulator

3. Generating bits using random function

%% Signal Source
% Create a binary data stream as a column vector.
x = randi([0 1],n,1); % Random binary data stream

4. Convert the bit patterns in to symbol patterns

%% Bit-to-Symbol Mapping
% Convert the bits in x into k-bit symbols.
xsym = bi2de(reshape(x,k,length(x)/k).’,’left-msb’);

5. Modulation using defined modulator (match the symbols onto the points in constellation)

%% Modulation
% Modulate using custom constellation.
y = modulate(hMod,xsym);

% Modulate using 8-PSK.
y_ref = modulate(hMod_ref,xsym);

6Direct the modulated signals into an AWGN noised channel

%% Transmitted Signals
yTx = y;
yTx_ref = y_ref;

%% Channel
% Send signal over an AWGN channel.

% EsNo = EbNo + 10*log10(k);In dB
% EsNo = 10*log(T_sym/T_samp) + SNR <–(A)
% nSamp = T_sym/T_samp
% Therefore, SNR + 10*log10(nSamp) = EbNo + 10*log10(k)
% Or from (A), SNR = EsNo – 10*log10(nSamp)
Es = 1/M*(4*(A/2)^2 + 4*A^2);
EsNo = Es/N;
SNR = 10*log10(EsNo) – 10*log10(nSamp);

% Send signal of custom constellation
yNoisy = awgn(yTx,SNR,’measured’);

% Send signal of 8-PSK
yNoisy_ref = awgn(yTx_ref,SNR,’measured’);

7. Demodulate the signals using defined demodulator and Convert the detected symbols into bits

%% Demodulation
% Demodulate signal using custom constellation.
zsym = demodulate(hDemod,yRx);
% convert integers to bits.
z = de2bi(zsym, ‘left-msb’);
% convert z from a matrix to a vector.
z = reshape(z.’, numel(z), 1);

% Demodulate signal using 8-PSK.
zsym_ref = demodulate(hDemod_ref,yRx_ref);
% convert integers to bits.
z_ref = de2bi(zsym_ref, ‘left-msb’);
% convert z from a matrix to a vector.
z_ref = reshape(z_ref.’, numel(z_ref), 1);

8.  Compare with the generated random bit stream (from step 3) and calculate the number of errors and so as the probability of errors

%% SER computation
% Symbol error rate of custom constellation
[number_of_errors_of_symbols_of_custom,symbol_error_rate_of_custom] = symerr(xsym,zsym)

% Symbol error rate of 8-PSK
[number_of_errors_of_symbols_of_8PSK,symbol_error_rate_of_8PSK] = symerr(xsym,zsym_ref)

9.  Plot Probability of errors vs Eb/No­  where each symbol implies usual meaning

%% Probability of symbol error for different values of EsNo
%No = 0.015:0.002:0.025;
%EsNo = Es./No;
%SNR = 10*log10(EsNo) – 10*log10(nSamp)
SNR = -5:1:15;
for i=1:length(SNR)
% Send signal of custom constellation
yNoisy(:,i) = awgn(yTx,SNR(i),’measured’);
end

yRx = yNoisy;
zsym = demodulate(hDemod,yRx);
[number_of_errors_of_sym,sym_error_rate] = symerr(xsym,zsym);

% Reference’s SER curve
SNR = -5 : 1 :15;
[prob,ser] = berawgn(SNR,’psk’,8,’nondiff’);

The end result of the Custom Constellation’s Sybol Error Rate (SER) with respect to 8-PSK SER is shown below

endresult

Hope you find this post useful. Let me know if you have any question on this post. Good day!

Setting Recursion Limit in XST

XST(Xilinx®) is the tool used in ISE® software to synthesize VHDL, Verilog or mixed languages and to create Xilinx® specific netlist files known as NGC files. To have an overview on XST® synthesis please refer the following link.

http://www.xilinx.com/support/documentation/sw_manuals/xilinx11/ise_c_using_xst_for_synthesis.htm

There was a misconception that hardware design is incapable of implementing recursion functions. Of course Verilog ’95 did not support recursion but later Verilog 2001 introduced Recursion to be generated in designs. But one thing that a real hardware programmer must keep in mind is hardware programming is more concerned on parallel execution rather than being sequential as in typical programming languages. But for hardware designs with tree structure it is better to have recursion in the design. One useful place related to recursion can be found in the following link.

http://www.beyond-circuits.com/wordpress/2009/01/recursive-modules/

I have a question to ask from you. Have you ever tried a recursive function (having a considerable depth) ? Then you would have noticed following error message (as i myself encountered).

Function ‘<your_function>’ has iterated 64 times. Use “set -recursion_iteration_limit XX” to iterate more.

And then probably you would have googled it and found the following link.

http://www.xilinx.com/support/answers/18429.htm

As a novice I felt bit hard on the implementation of what have been given in the above link. Yes again I googled it. Again I googled it and again and again. Even the error message itself and the link above are more than enough for a ‘veteran’ in ISE® environment, I had to. After I am done (not as fast as I am writing) with this issue I thought to blog this out (nothing miracle happened, believe me).

Please note that my version of ISE design suite is 14.6 and my OS is windows 7. The steps I recommend you to resolve the issue under this topic is listed down as follows.

Step I: Go to the working directory of your design (project).

Step II: Find and open the (say) “main_module.xst” file in your directory. (this file’s icon looks like an Microsoft Office Outlook file)

file

project’s .xst extension file – first few lines

Step III:  Copy and paste this file to a new file and save under a different name (say) “main_module_new.xst” in the same directory. (This is because to avoid the file being overwritten by the Project Navigator.)

Step IV: Insert the following line at third line in the file you created.

                                                                          set -recursion_iteration_limit 256

file2

Project’s newly modified .xst extension file – first few lines

Now what all you have to do is, generate all synthesized files. To do this,

Step V: Go to Start -> All Programs -> Xilinx Design Tools -> ISE Design Suite 14.6 -> Accessories and open ISE Design Suite Command Prompt.

Step VI: First of all navigate to the directory where your project lies.

Step VII: Run XST® on newly created .xst file (main_module_new.xst) using following command.

                                                                             xst -ifn main_module_new.xst

It will take some time to compile.

Step VIII: After that ‘Implement Design’ and other further activities can be carried out based on modified project files.

Thats all. Feel free to raise any question under this topic.

Determining The Resultant of Einthoven triangle

Einthoven triangle is a conceptual idea which comes in Bio medical Science. This concept is used for ECG interpretation as it makes an underline platform to determine resultant electrical axis of the heart. If you have a basic knowledge in ECG interpretation this short wordings would not be a matter of wasting time, i believe.

To have a nice explanation on axis determination in ECG have a look at below article.

http://ems12lead.com/tag/einthovens-triangle/

Though the article is about ECG interpretation, I won’t go through every aspects under this topic. There are lots of confusions has been built around interpreting the axis using Einthovens triangle as many a students have misinterpreted that this triangle is exactly an equilateral triangle as the name suggests. Even the triangle is called “Einthovens equilateral triangle”, the ‘equilateral’ is used to indicate the electrical equality which comes under Einthovens law (I + III = II).

V1 and V2 are components of ROne of the other misinterpretation (which i also came through) is interpreting the resultant electrical axis using the above triangle. Actually a fascinating feature of this triangle is that two leads’ electrode magnitudes are more than enough to get done the work under this topic. But the resultant won’t come as we think in basic physics. This is a point where the most of the confusion comes alive. Look at the example provided herewith to grasp the idea.

V1 = R*cos(a) and V2  = R*cos(b) if we have to determine the components of R along a and b angle apart as depicted above. Now how are you going to take the R from V1 and V2 ? You must be brave enough to understand the concept now. You won’t put Resultant equation to derive the R from V1 and V2. Why ? V1 and V2 are components of R itself and clearly a+b is not exactly 90 degrees (rather provided as such). To determine the R from V1 and V2 you should work backward.

Working backward would not be possible if we could not construct the picture in mind. (provided there are other methods as well). Only tricky part here is to determine the angles ‘a’ and ‘b’.

More precisely the Einthovens concept comes in to play under above conception. To determine the resultant electrical axis of the heart using Einthovens triangle is not just a calculation of resultant vector of lead vectors. One main thing that we have to keep in mind is that the Einthovens law is an scalar addition and not a vector addition. But in the calculation of the resultant electric axis the characteristics of vector should be considered. Therefore the above two calculations are deferred in their very nature.

Resultant of the Einthoven triangle is calculated as follows.

An example to overview the resultant calculation

The resultant deflections are plotted on the Einthoven Triangle as vectors, appropriate units are used to represent the millimeter or millivolt height of each deflection. Each vector is drawn along the side of the triangle corresponding to the lead it represents.  The vector is drawn from the centre of each side and toward the (+) pole or electrode.  If Einthoven’s law is valid, perpendicular lines dropped from the point inside the triangle.  A line drawn from this intersection to the centre of the triangle represents the electrical axis of the heart.  The position of this axis is measured in degrees from the horizontal line drawn through the centre of the triangle (see below).  The length of the axis vector is a measure of the overall electrical potential of the heart in the axis direction.

Image

(source:http://www.acad.polyu.edu.hk/~bcbio/doc/CVS2b.doc)

As depicted above the resultant electric axis of the heart can be found using the three Bi-polar limb electrodes. One last thing that should be noted here is as I suggested at the beginning of this article, for the determination of resultant only two BI-polar vectors are enough. And calculating the resultant axis using only two vectors can be used to derive the remaining lead vector from characteristics of vectors and hence can prove the law of Einthovens.

Hope things went fine. I am pleased to accept any views of you upon this article. Nice day!

“If we look inside the atom, any atom, we will see a sun in its core.” – Ali ibn Abi Thalib