From: Greg Burri <greg.burri@gmail.com>
Date: Mon, 23 Nov 2015 10:33:25 +0000 (+0100)
Subject: Add some tests on different variations of Hough algorithm.
X-Git-Tag: 1.0.11~87
X-Git-Url: http://git.euphorik.ch/?p=master-thesis.git;a=commitdiff_plain;h=0ff8fb82457bd5a858b2218ab07f69c81323537e

Add some tests on different variations of Hough algorithm.
---

diff --git a/Documentation/ma_autocount/plasmodium-autocount-0.1.py b/Documentation/ma_autocount/plasmodium-autocount-0.1.py
index b2fa3d2..d97e059 100644
--- a/Documentation/ma_autocount/plasmodium-autocount-0.1.py
+++ b/Documentation/ma_autocount/plasmodium-autocount-0.1.py
@@ -1,3 +1,4 @@
+# -*- coding: utf-8 -*-
 #!/usr/bin/env python
 
 #    Copyright (C) 2010 Monash University
@@ -23,7 +24,6 @@ Identify red blood cells and red blood cells infected with plasmodium in a set o
 """
 
 from PIL import Image
-
 import sys, numpy, os, random
 from scipy import ndimage
 
@@ -58,20 +58,32 @@ N_processors = 1
 #
 # What to show in diagnostic output images
 
-Show_edges = False
+Show_edges = True
 Show_stain = False
 Show_infection = False
 
-Show_hand_marks = True
+Show_hand_marks = False
 Show_computer_marks = True
-Show_cell_shapes = False
+Show_cell_shapes = True
 
 
 # =================
 # Image scaling
 
-Resize_width = 680
-Resize_height = 512
+Resize = False
+
+# Tailles pour les images dans '08.10.2015'
+Resize_width = 1296 #680
+Resize_height = 972 #512
+
+# Tailles pour les images dans '22.05.2015' (50%)
+#Resize_width = 1296
+#Resize_height = 972
+
+
+# [GB] Hemato
+# Resize_width = 2864
+# Resize_height = 2885
 
 
 # ==========================================
@@ -109,7 +121,7 @@ Infection_pixels_required = 1
 
 # Size of cells, for Hough transform
 Hough_min_radius = 5  
-Hough_max_radius = 20 
+Hough_max_radius = 32 # 20
 
 # Portion of cell used for classification as infected/weird
 # (may be abutting an infected cell, so be cautious)
@@ -117,7 +129,7 @@ Hough_max_radius = 20
 Cell_inside_radius = 0.9
 
 # Do not allow another cell center to be detected within this radius
-Cell_suppression_radius = 1.25
+Cell_suppression_radius = 1.05 # 1.25
 
 
 # =======================
@@ -127,17 +139,13 @@ Cell_suppression_radius = 1.25
 Max_dark_stain_proportion = 0.728
 
 # Minimum/maximum radius (as determined by Hough transform)
-Min_cell_radius = 9
-Max_cell_radius = 20
+Min_cell_radius = 11 # 9
+Max_cell_radius = 32 # 20
 
 # Proportion of radius that center of mass can be away from Hough determined center
 Max_offcenter = 0.5
 
 
-
-
-
-
 _offset_cache = { }
 def offsets(inner_radius, outer_radius):
     """
@@ -244,14 +252,14 @@ def process(filename_in, filename_out, filename_coords, name):
     print 'Processing', filename_in, '->', filename_out
 
     # Resize (mostly done for speed increase)
-    i = i.resize((Resize_width,Resize_height), Image.ANTIALIAS)
+    i = i.resize((Resize_width,Resize_height), Image.ANTIALIAS) if Resize else i
 
     # Convert to numpy array
-    a = numpy.asarray(i).astype('float64')
+    a = numpy.asarray(i).astype('float64') # [GB] Shape: (512, 680, 3)
     height = a.shape[0]
     width = a.shape[1]
 
-    a_gray = numpy.sum(a, 2)
+    a_gray = numpy.sum(a, 2) # [GB] Somme de chaque composante RGB (non-normalisé)
     
     # Make a copy of the array to doodle various diagnostic markings on
     a_output = a.copy()    
@@ -262,33 +270,39 @@ def process(filename_in, filename_out, filename_coords, name):
 
 
     # Denoise/unsharp mask
-    a_enhanced = gaussian_blur(a, Enhance_a) - Enhance_b*gaussian_blur(a,Enhance_c)
-    
+    # [GB] (Concerne les trois channels RGB).
+    a_enhanced = gaussian_blur(a, Enhance_a) - Enhance_b * gaussian_blur(a, Enhance_c)
     
     # Split into foreground and background using k-medians
     status('fg/bg')
     
+    # [GB] Chaque canal de couleur est linéarisé dans un tableau.
     a_raveled = numpy.reshape(a_enhanced,(width*height,a_enhanced.shape[2]))
-    average = numpy.average(a_raveled, axis=0)
-    # Initial guess
-    color_bg = average*1.5
-    color_fg = average*0.5
     
+    # [GB] Moyenne de chaque canal.
+    average = numpy.average(a_raveled, axis=0) 
+    
+    # Initial guess of the medians
+    color_bg = average * 1.5 # [GB] Plus clair
+    color_fg = average * 0.5 # [GB] Plus sombre
+    
+    # Fifty times: (0 to 4).
     for i in xrange(5):
-        d_bg = a_raveled-color_bg[None,:]
+        d_bg = a_raveled - color_bg[None,:]
         d_bg *= d_bg
-        d_bg = numpy.sum(d_bg,1)
-        d_fg = a_raveled-color_fg[None,:]
+        d_bg = numpy.sum(d_bg, 1)  # [GB] Somme des carrés des couleurs
+        d_fg = a_raveled - color_fg[None,:]
         d_fg *= d_fg
-        d_fg = numpy.sum(d_fg,1)
-        fg = d_fg*Background_weight < d_bg
-        color_fg = numpy.median(a_raveled[fg,:],axis=0)
+        d_fg = numpy.sum(d_fg, 1)
+        fg = d_fg * Background_weight < d_bg
+        color_fg = numpy.median(a_raveled[fg,:],axis=0) # [GB] axis = 0 -> par colonne (ordonnées). Le résultat est alors 3 valeurs, une pour chaque couleur
         color_bg = numpy.median(a_raveled[~fg,:],axis=0)
-        
-    fg = numpy.reshape(fg, (height,width))
+
+    fg = numpy.reshape(fg, (height,width)) # [GB] Image binaire du premier plan.
     
-    edge = ndimage.maximum_filter(fg,size=3) != ndimage.minimum_filter(fg,size=3)
+    edge = ndimage.maximum_filter(fg,size=3) != ndimage.minimum_filter(fg,size=3) # [GB] image binaire des contours du premier plan.
 
+    # [GB] délinéarisation des tableaux résultant du k-median (niveau de gris).
     d_fg = numpy.sqrt(numpy.reshape(d_fg, (height,width)))
     d_bg = numpy.sqrt(numpy.reshape(d_bg, (height,width)))
     
@@ -298,10 +312,15 @@ def process(filename_in, filename_out, filename_coords, name):
     # (should be invariant under changes in brightness) 
     mag_bg = numpy.average(color_bg)
     
+    # [GB] Les "dark stain" :
+    # * sont d'une intensité plus élevée que la moyenne du bg multiplié par un facteur (Dark_stain_level)
+    # * les valeurs rouges sont plus petites que la mediane du premier-plan pour les rouges
+    # * les valeurs vertes sont plus petites que la médiane du premier-plan pour les verts
+    # * appartiennent au foreground
     dark_stain = (
         (d_fg > mag_bg*Dark_stain_level) &
-        (a_enhanced[:,:,0] < color_fg[0]) & 
-        (a_enhanced[:,:,1] < color_fg[1]) &
+        (a_enhanced[:,:,0] < color_fg[0]) & # [GB] Red channel
+        (a_enhanced[:,:,1] < color_fg[1]) & # [GB] Green channel
         fg
     )
     
@@ -332,6 +351,8 @@ def process(filename_in, filename_out, filename_coords, name):
     if Show_infection:
         a_output[infection,0] += 128
         
+    #a_output[stain, 1] += 128
+        
         
     # Hough transform
     status('hough')
@@ -409,9 +430,11 @@ def process(filename_in, filename_out, filename_coords, name):
     candidates = [ ]
     for y in xrange(height):
         for x in xrange(width):
-            if hough[y,x] > thresh:
-                candidates.append((hough[y,x],y,x))
-    candidates.sort(reverse=True)
+            if hough[y, x] > thresh:
+                candidates.append((hough[y,x], y, x))
+    candidates.sort(reverse = True)
+    hough
+    
     
     status('classify')
     
@@ -420,21 +443,20 @@ def process(filename_in, filename_out, filename_coords, name):
     
     suppress = numpy.zeros((height,width), 'bool')
     claim_strength = numpy.zeros((height,width), 'float64') + 1e30
-    cells = [ ] # [(y,x)]
+    cells = [] # [(y,x)]
     n = 0
     n_infected = 0
-    n_late = 0
     for _,y,x in candidates:
         radius = hough_radius[y,x]        
-
+        
         # Make sure candidate is not near a previously detected cell,
         # and does not abut the edge of the image        
         if not suppress[y,x] and \
            y >= radius and x >= radius and y < height-radius and x < width-radius:
-            coords = clipped_coordinates(0,radius*Cell_inside_radius, y,x,height,width)
+            coords = clipped_coordinates(0, radius * Cell_inside_radius, y, x, height, width)
             
             #Only include foreground
-            this_fg = fg[coords[:,0],coords[:,1]]
+            this_fg = fg[coords[:,0], coords[:,1]] # [GB] Booleans
             coords = coords[this_fg]
             
             this_strength = ((coords[:,0]-y)**2 + (coords[:,1]-x)**2)
@@ -443,13 +465,12 @@ def process(filename_in, filename_out, filename_coords, name):
                 this_strength
             )
             
-            cells.append((y,x, coords,this_strength))
+            cells.append((y, x, coords, this_strength))
 
-           
+       
         # Suppress detection of cells too close to here
-        # (note: we do this irrespective of whether we treated this location as a hit)
-        
-        coords = clipped_coordinates(0,radius*Cell_suppression_radius, y,x,height,width)
+        # (note: we do this irrespective of whether we treated this location as a hit)        
+        coords = clipped_coordinates(0,radius * Cell_suppression_radius, y,x,height,width)
         suppress[coords[:,0],coords[:,1]] = True
         
         #for oy, ox in offsets(0, radius*Cell_suppression_radius): #suppression_offsets:
@@ -525,17 +546,18 @@ def process(filename_in, filename_out, filename_coords, name):
             color = (0,0,255)
             file_coords.write('%d,%d,Uninfected\n' % (x,y))
         
+        if Show_cell_shapes:            
+            a_output[coords[:,0],coords[:,1]] += [[ random.random()*64+32 for i in (0,1,2) ]]
+            
         # Mark detected cell on output image
         if Show_computer_marks:
-            set(a_output,y,x,color)
+            set(a_output, y, x, color)
             for i in xrange(1,3):
                 set(a_output,y-i,x,color)  
                 set(a_output,y+i,x,color)  
                 set(a_output,y,x-i,color)  
                 set(a_output,y,x+i,color)
 
-        if Show_cell_shapes:            
-            a_output[coords[:,0],coords[:,1]] += [[ random.random()*64+32 for i in (0,1,2) ]]
             
         
     
diff --git a/src/Tests_hough/HTCircle.m b/src/Tests_hough/HTCircle.m
new file mode 100644
index 0000000..3088c24
--- /dev/null
+++ b/src/Tests_hough/HTCircle.m
@@ -0,0 +1,42 @@
+% Hough Transform for circles.
+% img : logical image, zeros are the circles edges
+% radius_range is a vector : [min, max]
+function [acc_votes, acc_radius] = HTCircle(img, radius_range)
+
+    [rows, columns] = size(img);
+    
+    acc_votes = zeros(rows, columns);
+    acc_radius = zeros(rows, columns);
+    
+    phi_range = linspace(0, 2 * pi, 100);
+    phi_range(end) = [];
+    
+    indexes_non_zero = find(img)';
+    s = size(img);
+    
+    for r = radius_range(1):radius_range(2)
+        acc = zeros(rows, columns);
+        for i = indexes_non_zero
+            [y, x] = ind2sub(s, i);            
+            for phi = phi_range
+                x0 = round(x - r * cos(phi));
+                y0 = round(y - r * sin(phi));
+                if x0 < columns && x0 > 0 && y0 < rows && y0 > 0
+                    acc(y0, x0) = acc(y0, x0) + 1;
+                end
+            end
+        end
+        
+        for x = 1:columns
+            for y = 1:rows
+                if acc(y, x) > acc_votes(y, x)
+                    acc_votes(y, x) = acc(y, x);
+                    acc_radius(y, x) = r;
+                end
+            end
+        end        
+    end
+    
+    acc_votes = imgaussfilt(acc_votes, 1.0);
+end
+
diff --git a/src/Tests_hough/HTCircleGradient.m b/src/Tests_hough/HTCircleGradient.m
new file mode 100644
index 0000000..fe20e20
--- /dev/null
+++ b/src/Tests_hough/HTCircleGradient.m
@@ -0,0 +1,42 @@
+% Hough Transform for circles with gradient.
+% img : gradient image
+% radius_range is a vector : [min, max]
+function [acc_votes, acc_radius] = HTCircleGradient(img, radius_range)
+
+    [rows, columns] = size(img);
+    
+    acc_votes = zeros(rows, columns);
+    acc_radius = zeros(rows, columns);
+    
+    phi_range = linspace(0, 2 * pi, 100);
+    phi_range(end) = [];
+    
+    indexes_non_zero = find(img)';
+    s = size(img);
+    
+    for r = radius_range(1):radius_range(2)
+        acc = zeros(rows, columns);
+        for i = indexes_non_zero            
+            [y, x] = ind2sub(s, i);  
+            for phi = phi_range
+                x0 = round(x - r * cos(phi));
+                y0 = round(y - r * sin(phi));
+                if x0 < columns && x0 > 0 && y0 < rows && y0 > 0
+                    acc(y0, x0) = acc(y0, x0) + img(y, x);
+                end
+            end
+        end
+        
+        for x = 1:columns
+            for y = 1:rows
+                if acc(y, x) > acc_votes(y, x)
+                    acc_votes(y, x) = acc(y, x);
+                    acc_radius(y, x) = r;
+                end
+            end
+        end        
+    end
+    
+    acc_votes = imgaussfilt(acc_votes, 1.0);
+end
+
diff --git a/src/Tests_hough/HTCircleMa.m b/src/Tests_hough/HTCircleMa.m
new file mode 100644
index 0000000..3afa97f
--- /dev/null
+++ b/src/Tests_hough/HTCircleMa.m
@@ -0,0 +1,41 @@
+% Hough Transform for circles.
+% x_dir : the x component of the normalized gradient
+% y_dir : the y component of the normalized gradient
+% edges : a binary image with edge = 1
+% radius_range is a vector : [min, max]
+function [acc_votes, acc_radius] = HTCircleMa(x_dir, y_dir, edges, radius_range)
+
+    [rows, columns] = size(edges);
+    
+    acc_votes = zeros(rows, columns);
+    acc_radius = zeros(rows, columns);    
+    
+    indexes_edges = find(edges)';
+    s = size(edges);
+    
+    for r = radius_range(1):radius_range(2)
+        acc = zeros(rows, columns);
+        for i = indexes_edges
+            [y, x] = ind2sub(s, i);      
+            
+            x0 = round(x + x_dir(i) * r);
+            y0 = round(y + y_dir(i) * r);
+            
+            if x0 < columns && x0 > 0 && y0 < rows && y0 > 0
+                acc(y0, x0) = acc(y0, x0) + 1;
+            end
+        end
+        
+        for x = 1:columns
+            for y = 1:rows
+                if acc(y, x) > acc_votes(y, x)
+                    acc_votes(y, x) = acc(y, x);
+                    acc_radius(y, x) = r;
+                end
+            end
+        end        
+    end
+    
+    acc_votes = imgaussfilt(acc_votes, 1.0);
+end
+
diff --git a/src/Tests_hough/HTEllipse.m b/src/Tests_hough/HTEllipse.m
new file mode 100644
index 0000000..908f230
--- /dev/null
+++ b/src/Tests_hough/HTEllipse.m
@@ -0,0 +1,56 @@
+% Hough Transform for circles.
+% img : logical image, zeros are the circles edges
+% radius_range is a vector : [min, max]
+function [acc_votes, acc_radius1, acc_radius2, acc_alpha] = HTEllipse(img, radius_range)
+
+    [rows, columns] = size(img);
+    
+    acc_votes = zeros(rows, columns);
+    acc_radius1 = zeros(rows, columns);
+    acc_radius2 = zeros(rows, columns);
+    acc_alpha = zeros(rows, columns);
+    
+    alpha_range = linspace(0, pi / 2, 8);
+    alpha_range(end) = [];
+    
+    phi_range = linspace(0, 2 * pi, 32);
+    phi_range(end) = [];
+    
+    indexes_non_zero = find(img)';
+    s = size(img);
+    
+    for alpha = alpha_range
+        transform_mat = [cos(alpha) -sin(alpha); sin(alpha) cos(alpha)];        
+        for r1 = radius_range(1):radius_range(2)
+            for r2 = radius_range(1):radius_range(2)
+                acc = zeros(rows, columns);
+                for i = indexes_non_zero
+                    [y, x] = ind2sub(s, i);
+                    for phi = phi_range
+                        p0 = [x; y] - transform_mat * [r1 * cos(phi); r2 * sin(phi)];
+                        x0 = round(p0(1));
+                        y0 = round(p0(2));
+
+                        if x0 < columns && x0 > 0 && y0 < rows && y0 > 0
+                            acc(y0, x0) = acc(y0, x0) + 1;
+                        end
+                    end
+                end
+
+                for x = 1:columns
+                    for y = 1:rows
+                        if acc(y, x) > acc_votes(y, x)
+                            acc_votes(y, x) = acc(y, x);
+                            acc_radius1(y, x) = r1;
+                            acc_radius2(y, x) = r2;
+                            acc_alpha(y, x) = alpha;
+                        end
+                    end
+                end      
+            end
+        end
+    end
+    
+    acc_votes = imgaussfilt(acc_votes, 1.0);
+end
+
diff --git a/src/Tests_hough/HTEllipse2.m b/src/Tests_hough/HTEllipse2.m
new file mode 100644
index 0000000..bea1243
--- /dev/null
+++ b/src/Tests_hough/HTEllipse2.m
@@ -0,0 +1,75 @@
+% Hough Transform for circles.
+% img : logical image, zeros are the circles edges
+% radius_range is a vector : [min, max]
+function [acc_votes, acc_radius1, acc_radius2, acc_alpha] = HTEllipse2(x_dir, y_dir, edges, radius_range)
+
+    [rows, columns] = size(edges);
+    
+    acc_votes = zeros(rows, columns);
+    acc_radius1 = zeros(rows, columns);
+    acc_radius2 = zeros(rows, columns);
+    acc_alpha = zeros(rows, columns);
+    
+    alpha_range = linspace(0, pi / 2, 8);
+    alpha_range(end) = [];
+%     alpha_range = pi/8;
+        
+    indexes_edges = find(edges)';
+    s = size(edges);
+    
+    for alpha = alpha_range
+        rot_mat = [cos(alpha) -sin(alpha); sin(alpha) cos(alpha)];
+        rot_mat_inv = [cos(-alpha) -sin(-alpha); sin(-alpha) cos(-alpha)];      
+        for r1 = radius_range(1):radius_range(2)
+            for r2 = radius_range(1):radius_range(2)
+                acc = zeros(rows, columns);
+                for i = indexes_edges
+                    [y, x] = ind2sub(s, i);
+                    v = [-x_dir(i); -y_dir(i)];
+                    u = rot_mat_inv * v;
+                    if u(1) == 0
+                        if u(2) > 0
+                            theta = pi / 2.;
+                        else
+                            theta = 3. * pi / 2.;
+                        end
+                    elseif u(2) == 0
+                        if u(1) > 0
+                            theta = 0;
+                        else
+                            theta = pi;
+                        end
+                    elseif u(2) < 0
+                        t = u(2) / u(1);
+                        theta = 2*pi + 2 * atan(...
+                            (-r1 - r2*t*sqrt((r1^2+r2^2*t^2)/(r2^2*t^2)))...
+                            /(r2*t));
+                    else
+                        t = u(2) / u(1);
+                        theta = 2 * atan(...
+                            (-r1 + r2*t*sqrt((r1^2+r2^2*t^2)/(r2^2*t^2)))...
+                            /(r2*t));
+                    end
+                    C = round([x; y] - rot_mat * [r1*cos(theta); r2*sin(theta)]);
+                    if C(1) < columns && C(1) > 0 && C(2) < rows && C(2) > 0
+                        acc(C(2),C(1)) = acc(C(2),C(1)) + 1;
+                    end
+                end
+
+                for x = 1:columns
+                    for y = 1:rows
+                        if acc(y, x) > acc_votes(y, x)
+                            acc_votes(y, x) = acc(y, x);
+                            acc_radius1(y, x) = r1;
+                            acc_radius2(y, x) = r2;
+                            acc_alpha(y, x) = alpha;
+                        end
+                    end
+                end      
+            end
+        end
+    end
+    
+    acc_votes = imgaussfilt(acc_votes, 1.0);
+end
+
diff --git a/src/Tests_hough/draw_ellipses.m b/src/Tests_hough/draw_ellipses.m
new file mode 100644
index 0000000..d793f4d
--- /dev/null
+++ b/src/Tests_hough/draw_ellipses.m
@@ -0,0 +1,44 @@
+% center : [x y]
+function draw_ellipses(axe, ellipses, plot_3d, plot_center)
+    if ~exist('plot_3d', 'var')
+        plot_3d = false;
+    end
+    
+    if ~exist('plot_center', 'var')
+        plot_center = false;
+    end
+
+    hold on;
+    
+    if plot_3d
+        z_space = [0 1];
+    else
+        z_space = [0];
+    end
+    
+    for i = 1:length(ellipses)
+        e = ellipses{i};
+        transform_mat = [cos(e.alpha) -sin(e.alpha); sin(e.alpha) cos(e.alpha)]; 
+
+        p_previous = [];
+        for phi = linspace(0, 2 * pi, 36);
+           for z = z_space
+                p = [e.x0; e.y0] + transform_mat * [e.r1 * cos(phi); e.r2 * sin(phi)]; % p : [x y]
+                if ~isempty(p_previous)
+                    line = plot3(axe, [p(1) p_previous(1)], [p(2) p_previous(2)], [z z], 'LineWidth', 1, 'Color', 'yellow');
+                    line.Color(4) = 0.5;
+                end
+                if z == 1
+                   line = plot3(axe, [p(1) p(1)], [p(2) p(2)], [0 1], 'LineWidth', 0.5, 'Color', 'yellow'); 
+                   line.Color(4) = 0.5;
+                end
+           end
+           p_previous = p;
+        end
+        
+        if plot_center
+            plot(e.x0, e.y0, 'r.', 'MarkerSize', 20);
+        end
+    end
+end
+
diff --git a/src/Tests_hough/expand_matrix.m b/src/Tests_hough/expand_matrix.m
new file mode 100644
index 0000000..b165902
--- /dev/null
+++ b/src/Tests_hough/expand_matrix.m
@@ -0,0 +1,4 @@
+function [m] = expand_matrix(m)
+    m = [m(:,1) m m(:,end)];
+    m = [m(1,:); m; m(end,:)];
+end
\ No newline at end of file
diff --git a/src/Tests_hough/find_ellipses.m b/src/Tests_hough/find_ellipses.m
new file mode 100644
index 0000000..4c5e7a0
--- /dev/null
+++ b/src/Tests_hough/find_ellipses.m
@@ -0,0 +1,60 @@
+function [ellipses] = find_ellipses(acc_votes, acc_radius_1, acc_radius_2, acc_alpha)
+    if ~exist('acc_radius_2', 'var')
+        acc_radius_2 = acc_radius_1;    
+    end
+
+    if ~exist('acc_alpha', 'var')
+        acc_alpha = zeros(size(acc_votes));
+    end
+
+    ellipses = {}; % struct('x0', 'y0', 'r1', 'r2', 'alpha')
+    max_votes = max(acc_votes(:));
+    min_votes = min(acc_votes(:));
+    threshold_votes = (max_votes - min_votes) * 0.7 + min_votes;
+    
+    acc_votes_suppressed = acc_votes;
+    while true
+        [votes, index_max] = max(acc_votes_suppressed(:));
+        if votes <= threshold_votes
+            break;
+        end
+        
+        acc_votes_suppressed(index_max) = 0; % Suppress the vote.
+        
+        [y, x] = ind2sub(size(acc_votes_suppressed), index_max);        
+        
+        % If the center ellipse isn't in a known ellipse then add it to 'ellipses'.
+        accept_ellipse = true;
+        for i = 1:length(ellipses)
+            ellipse = ellipses{i};
+%             module = sqrt((x - ellipse.x)^2 + (y - ellipse.y)^2);
+%             phi = asin((y - ellipse.y) / module);
+%             
+%             p = [cos(a) -sin(a); sin(a) cos(a)] * [ellipse.r1 * cos(phi); ellipse.r2 * sin(phi)];
+%             module_ellipse = sqrt(p(1)^2 + p(2)^2);                
+            
+            a = ellipse.alpha;
+            x0 = ellipse.x0;
+            y0 = ellipse.y0;
+            r1 = ellipse.r1;
+            r2 =  ellipse.r2;
+            n = ((x - x0) * cos(a) + (y - y0) * sin(a))^2 / r1^2 + ((x - x0) * sin(a) - (y - y0) * cos(a))^2 / r2^2;
+            if n <= 1
+                accept_ellipse = false;
+                break;
+            end
+        end
+        
+        if accept_ellipse            
+            e = struct(...
+                'x0', x,...
+                'y0', y,...
+                'r1', acc_radius_1(index_max),...
+                'r2', acc_radius_2(index_max),...
+                'alpha', acc_alpha(index_max));            
+            ellipses{end+1} = e;
+            fprintf('Ellipse : (x0: %.1f, y0: %.1f, r1: %.1f, r2: %.1f, alpha: %.2f)\n', e.x0, e.y0, e.r1, e.r2, e.alpha);
+        end
+    end
+end
+
diff --git a/src/Tests_hough/images/heap_1.png b/src/Tests_hough/images/heap_1.png
new file mode 100644
index 0000000..beaea3c
Binary files /dev/null and b/src/Tests_hough/images/heap_1.png differ
diff --git a/src/Tests_hough/images/rbc_single.png b/src/Tests_hough/images/rbc_single.png
new file mode 100644
index 0000000..dc9424a
Binary files /dev/null and b/src/Tests_hough/images/rbc_single.png differ
diff --git a/src/Tests_hough/images/rbc_single_blurred.png b/src/Tests_hough/images/rbc_single_blurred.png
new file mode 100644
index 0000000..77afb77
Binary files /dev/null and b/src/Tests_hough/images/rbc_single_blurred.png differ
diff --git a/src/Tests_hough/images/rbc_single_oblong_1.png b/src/Tests_hough/images/rbc_single_oblong_1.png
new file mode 100644
index 0000000..549488a
Binary files /dev/null and b/src/Tests_hough/images/rbc_single_oblong_1.png differ
diff --git a/src/Tests_hough/images/rbc_single_oblong_2.png b/src/Tests_hough/images/rbc_single_oblong_2.png
new file mode 100644
index 0000000..4fb2fb6
Binary files /dev/null and b/src/Tests_hough/images/rbc_single_oblong_2.png differ
diff --git a/src/Tests_hough/images/rbc_single_oblong_3.png b/src/Tests_hough/images/rbc_single_oblong_3.png
new file mode 100644
index 0000000..6e8abdb
Binary files /dev/null and b/src/Tests_hough/images/rbc_single_oblong_3.png differ
diff --git a/src/Tests_hough/images/two_rbc_1.png b/src/Tests_hough/images/two_rbc_1.png
new file mode 100644
index 0000000..bb67d0d
Binary files /dev/null and b/src/Tests_hough/images/two_rbc_1.png differ
diff --git a/src/Tests_hough/images/two_rbc_2.png b/src/Tests_hough/images/two_rbc_2.png
new file mode 100644
index 0000000..23f68b1
Binary files /dev/null and b/src/Tests_hough/images/two_rbc_2.png differ
diff --git a/src/Tests_hough/images/two_rbc_3.png b/src/Tests_hough/images/two_rbc_3.png
new file mode 100644
index 0000000..1148f53
Binary files /dev/null and b/src/Tests_hough/images/two_rbc_3.png differ
diff --git a/src/Tests_hough/images/wbc_1.png b/src/Tests_hough/images/wbc_1.png
new file mode 100644
index 0000000..947ad39
Binary files /dev/null and b/src/Tests_hough/images/wbc_1.png differ
diff --git a/src/Tests_hough/images/wbc_2.png b/src/Tests_hough/images/wbc_2.png
new file mode 100644
index 0000000..611bd3a
Binary files /dev/null and b/src/Tests_hough/images/wbc_2.png differ
diff --git a/src/Tests_hough/results/notes.txt b/src/Tests_hough/results/notes.txt
new file mode 100644
index 0000000..e69de29
diff --git a/src/Tests_hough/show_figures.m b/src/Tests_hough/show_figures.m
new file mode 100644
index 0000000..ed6de4e
--- /dev/null
+++ b/src/Tests_hough/show_figures.m
@@ -0,0 +1,9 @@
+function [hlink] = show_figures(name, images, original, ellipses)    
+    [hlink, axes] = show_votes_figure(name, images);    
+    draw_ellipses(axes(4), ellipses, true);
+    savefig(['results/' name '.fig']);
+
+    show_image_ellipses(name, original, ellipses);
+    savefig(['results/' name ' - output.fig']);
+end
+
diff --git a/src/Tests_hough/show_image_ellipses.m b/src/Tests_hough/show_image_ellipses.m
new file mode 100644
index 0000000..229d524
--- /dev/null
+++ b/src/Tests_hough/show_image_ellipses.m
@@ -0,0 +1,11 @@
+function  show_image_ellipses(name, img, ellipses)
+    f = figure;
+    f.Position = [10 10 700 700];
+    f.Name = name;
+    grayscale = linspace(0, 1, 255);
+    colormap([grayscale' grayscale' grayscale']);
+    image(img);
+    axis image;
+    draw_ellipses(f.CurrentAxes, ellipses, false, true);
+end
+
diff --git a/src/Tests_hough/show_votes_figure.m b/src/Tests_hough/show_votes_figure.m
new file mode 100644
index 0000000..f6856d4
--- /dev/null
+++ b/src/Tests_hough/show_votes_figure.m
@@ -0,0 +1,26 @@
+function [hlink, axes] = show_votes_figure(name, images)
+    f = figure;
+    f.Name = name;
+    f.Position = [100 100 1300 900];
+    
+    norm = @(A) (A - min(A(:)))/(max(A(:)) - min(A(:)));
+        
+    for i = 1:size(images, 2)            
+        axes(i) = subplot(2, 2, i);        
+        %colormap(hsv);
+        surf(norm(double(images{i}{1})), 'FaceColor','interp', 'FaceLighting', 'gouraud', 'EdgeColor', 'none'); 
+        camlight left;
+        shading interp;
+        colormap(jet(10))
+        axis tight;
+        title(images{i}{2});
+        view(axes(i), [-39.5 70]);        
+    end
+    
+    hlink = linkprop(axes, {'CameraPosition', 'CameraUpVector'});
+    
+    rotate3d on;    
+end
+
+
+
diff --git a/src/Tests_hough/test_all.m b/src/Tests_hough/test_all.m
new file mode 100644
index 0000000..5348923
--- /dev/null
+++ b/src/Tests_hough/test_all.m
@@ -0,0 +1,38 @@
+clear;
+
+images_path = 'D:\Master Thesis\src\Tests_hough\images\';
+image_names = {
+%     'heap_1.png'
+%     
+%     'rbc_single.png'
+%     'rbc_single_blurred.png'
+    
+%     'rbc_single_oblong_1.png'
+    'rbc_single_oblong_2.png'
+%     'rbc_single_oblong_3.png'
+%     
+%     'two_rbc_1.png'
+%     'two_rbc_2.png'
+%     'two_rbc_3.png'
+%     'wbc_1.png'
+%     'wbc_2.png'
+    };
+
+for k = 1:length(image_names)
+    image_name = image_names{k};
+    
+%     test_hough_simple(images_path, image_name);
+%     close all;
+
+%     test_hough_gradient(images_path, image_name, false)
+%     close all;
+
+%    test_hough_ma(images_path, image_name)
+% 	close all;
+    
+%     test_hough_ellipses(images_path, image_name);
+%     close all;
+    
+     test_hough_ellipses2(images_path, image_name);
+%     close all;
+end
\ No newline at end of file
diff --git a/src/Tests_hough/test_hough_ellipses.m b/src/Tests_hough/test_hough_ellipses.m
new file mode 100644
index 0000000..490b67d
--- /dev/null
+++ b/src/Tests_hough/test_hough_ellipses.m
@@ -0,0 +1,22 @@
+function test_hough_ellipses(images_path, image_name)
+    i = imread([images_path image_name]);
+    green = i(:,:,2);
+    green_filtered = imgaussfilt(double(green), 1.0);
+    edges = edge(green_filtered);
+
+    disp('Applying Houg...');
+    tic
+    [acc_votes, acc_radius1, acc_radius2, acc_alpha] = HTEllipse(edges, [20, 40]);
+    toc
+
+    disp('Finding ellipses...');
+    tic
+    ellipses = find_ellipses(acc_votes, acc_radius1, acc_radius2, acc_alpha);
+    toc
+    
+    show_figures(['Hough ellipses - ' image_name], {...
+        {green, 'original'},...
+        {green_filtered, 'original filtered with a gaussian, sigma=1'},...
+        {edges, 'edges (sobel)'},...
+        {expand_matrix(acc_votes), 'votes'}}, green, ellipses);
+end
\ No newline at end of file
diff --git a/src/Tests_hough/test_hough_ellipses2.m b/src/Tests_hough/test_hough_ellipses2.m
new file mode 100644
index 0000000..41be4a8
--- /dev/null
+++ b/src/Tests_hough/test_hough_ellipses2.m
@@ -0,0 +1,49 @@
+function test_hough_ellipses2(images_path, image_name)
+    i = imread([images_path image_name]);
+
+    green = i(:,:,2);
+    green_filtered = imgaussfilt(green, 1.0);
+    otsu = ~im2bw(green_filtered, graythresh(green_filtered));
+
+    se = strel('rectangle', [3 3]);
+    edges = imdilate(otsu, se) - imerode(otsu, se);
+
+    s = [-1 0 1; -2 0 2; -1 0 1];
+    x_edge = conv2(double(green_filtered), s, 'valid');
+    y_edge = conv2(double(green_filtered), s', 'valid');
+
+    magnitudes = sqrt(x_edge .^ 2 + y_edge .^ 2);
+
+    % Unit gradient.
+    x_dir = x_edge ./ magnitudes; 
+    y_dir = y_edge ./ magnitudes;
+
+    disp('Applying Houg...');
+    tic
+    [acc_votes, acc_radius1, acc_radius2, acc_alpha] = HTEllipse2(x_dir, y_dir, edges(2:end-1, 2:end-1), [20, 40]);
+    toc
+
+    disp('Finding ellipses...');
+    tic
+    ellipses = find_ellipses(acc_votes, acc_radius1, acc_radius2, acc_alpha);
+    toc
+    
+    name = ['Hough ellipses 2 - ' image_name];
+    show_figures(name, {...
+        {green_filtered, 'original filtered with a gaussian, sigma=1'},...
+        {edges, 'edges (otsu -> dilate - erode)'},...
+        {expand_matrix(x_edge), 'composante x du gradient'},...
+        {expand_matrix(acc_votes), 'votes'}}, green, ellipses);
+    
+    f = figure;
+    img = green_filtered;
+    grayscale = linspace(0, 1, 255);
+    colormap([grayscale' grayscale' grayscale']);
+    img(~edges) = 0;
+    image(img);
+    axis image;
+    hold on;
+    [X, Y] = meshgrid(1:size(green, 2), 1:size(green, 1));
+    quiver(X, Y, expand_matrix(x_dir), expand_matrix(y_dir), 'Color', 'yellow');
+    savefig(['results/' name ' - gradient.fig']);
+end
\ No newline at end of file
diff --git a/src/Tests_hough/test_hough_gradient.m b/src/Tests_hough/test_hough_gradient.m
new file mode 100644
index 0000000..15ca92b
--- /dev/null
+++ b/src/Tests_hough/test_hough_gradient.m
@@ -0,0 +1,45 @@
+function test_hough_gradient(images_path, image_name, use_edges_mask)
+    i = imread([images_path image_name]);
+
+    green = i(:,:,2);
+    green_filtered = imgaussfilt(green, 1.0);
+
+    edges = zeros(size(green));
+    if use_edges_mask
+        otsu = ~im2bw(green_filtered, graythresh(green_filtered));  
+        se = strel('rectangle', [3 3]);
+        edges = imdilate(otsu, se) - imerode(otsu, se);
+    end
+
+    s = [-1 0 1; -2 0 2; -1 0 1];
+    x_edge = conv2(cast(green_filtered, 'double'), s , 'valid');
+    y_edge = conv2(cast(green_filtered, 'double'), s', 'valid');
+
+    magnitudes = sqrt(x_edge .^ 2 + y_edge .^ 2);
+
+    if use_edges_mask
+        magnitudes(~edges(2:end-1, 2:end-1)) = 0;
+    end
+
+    disp('Applying Houg...');
+    tic
+    [acc_votes, acc_radius] = HTCircleGradient(magnitudes, [25, 35]);
+    toc
+
+    disp('Finding ellipses...');
+    tic
+    ellipses = find_ellipses(acc_votes, acc_radius);
+    toc
+
+    name = ['Hough gradient - ' image_name];
+    if use_edges_mask
+        name = [name ' (with edges)'];
+    end
+
+    show_figures(name, {...
+        {green_filtered, 'original filtered with a gaussian, sigma=1'},...
+        {edges, 'edges (otsu -> dilate - erode)'},...
+        {magnitudes, 'magnitudes'},...
+        {expand_matrix(acc_votes), 'votes'}}, green, ellipses);
+end
+
diff --git a/src/Tests_hough/test_hough_ma.m b/src/Tests_hough/test_hough_ma.m
new file mode 100644
index 0000000..690d3b2
--- /dev/null
+++ b/src/Tests_hough/test_hough_ma.m
@@ -0,0 +1,49 @@
+function test_hough_ma(images_path, image_name)
+    i = imread([images_path image_name]);
+
+    green = i(:,:,2);
+    green_filtered = imgaussfilt(green, 1.0);
+    otsu = ~im2bw(green_filtered, graythresh(green_filtered));
+
+    se = strel('rectangle', [3 3]);
+    edges = imdilate(otsu, se) - imerode(otsu, se);
+
+    s = [-1 0 1; -2 0 2; -1 0 1];
+    x_edge = conv2(double(green_filtered), s, 'valid');
+    y_edge = conv2(double(green_filtered), s', 'valid');
+
+    magnitudes = sqrt(x_edge .^ 2 + y_edge .^ 2);
+
+    % Unit gradient.
+    x_dir = x_edge ./ magnitudes; 
+    y_dir = y_edge ./ magnitudes;
+
+    disp('Applying Houg...');
+    tic
+    [acc_votes, acc_radius] = HTCircleMa(x_dir, y_dir, edges(2:end-1, 2:end-1), [25, 35]);
+    toc
+
+    disp('Finding ellipses...');
+    tic
+    ellipses = find_ellipses(acc_votes, acc_radius);
+    toc
+
+    name = ['Hough Ma - ' image_name];
+    show_figures(name,...
+        {{green_filtered, 'original filtered with a gaussian, sigma=1'},...
+        {edges, 'edges (otsu -> dilate - erode)'},...
+        {expand_matrix(x_edge), 'composante x du gradient'},...
+        {expand_matrix(acc_votes), 'votes'}}, green, ellipses);
+
+    f = figure;
+    img = green_filtered;
+    grayscale = linspace(0, 1, 255);
+    colormap([grayscale' grayscale' grayscale']);
+    img(~edges) = 0;
+    image(img);
+    axis image;
+    hold on;
+    [X, Y] = meshgrid(1:size(green, 2), 1:size(green, 1));
+    quiver(X, Y, expand_matrix(x_dir), expand_matrix(y_dir), 'Color', 'yellow');
+    savefig(['results/' name ' - gradient.fig']);
+end
diff --git a/src/Tests_hough/test_hough_simple.m b/src/Tests_hough/test_hough_simple.m
new file mode 100644
index 0000000..2725fc9
--- /dev/null
+++ b/src/Tests_hough/test_hough_simple.m
@@ -0,0 +1,22 @@
+function test_hough_simple(images_path, image_name)
+    i = imread([images_path image_name]);
+    green = i(:, :, 2);
+    green_filtered = imgaussfilt(cast(i(:, :, 2), 'double'), 1.0);
+    edges = edge(green_filtered);
+
+    disp('Applying Houg...');
+    tic
+    [acc_votes, acc_radius] = HTCircle(edges, [20, 40]);
+    toc
+
+    disp('Finding ellipses...');
+    tic
+    ellipses = find_ellipses(acc_votes, acc_radius);
+    toc
+
+    show_figures(['Hough simple - ' image_name], {...
+        {green, 'original'},...
+        {green_filtered, 'original filtered with a gaussian, sigma=1'},...
+        {edges, 'edges (sobel)'},...
+        {acc_votes, 'votes'}}, green, ellipses);
+end
\ No newline at end of file