Binary processing

Moments of Regions

Moments

The moment of a region is computed thanks to next equation:
$m_{p q} = \sum_{x = 0}^{M} \sum_{y = 0}^{N} I (x, y) \cdot x^{p} y^{q}$

Aera of regions

The aera of a regions can be compute with the moment of order p=0 et q=0:
$A (R) = \sum_{x = 0}^{M} \sum_{y + 0}^{N} I (x, y) \cdot x^{0} y^{0} = m_{00} (R)$
where A(R) is the aera of the region.
Here is the program:

clear all; close all; %x = imread('objetinv.PNG'); x = imread('objetinv2.PNG'); %x = imread('forme.PNG'); x=x(:,:,1); [M N] = size(x); %binarisation level = graythresh(x); img = im2bw(x,0.1); %labellisation CC = bwconncomp(img); L = labelmatrix(CC); %affihcage figure subplot(221) imshow(img) subplot(222) imshow(label2rgb(L)); %nombre d'objet, number of objects nbobj = max(max(L)); %parametres des objets, object's parameters param=zeros(nbobj,6); for i = 1 : nbobj param(i,1)=i; end %param = [ numero , nb pixel , boite_gauche , boite_droite , boite_haut , %boite_bas] for i = 1 : M for j = 1 : N if L(i,j)~=0 numero = L(i,j); param(numero,2)= param(numero,2) + 1 ; if param(numero,3)==0 || param(numero,3)>j param(numero,3)=j; end if param(numero,4)==0 || param(numero,4)<j param(numero,4)=j; end if param(numero,5)==0 || param(numero,5)>i param(numero,5)=i; end if param(numero,6)==0 || param(numero,6)<i param(numero,6)=i; end end end end % supression petit objet var =1; for i = 1 : nbobj if param(i,2)>10 parametre(var,:)=param(i,:); var=var+1; end end nbobj = size(parametre,1); %########### %########### %%%%%%moment moment_0_0=zeros(nbobj,1); p=0; q=0; for i = 1 : nbobj boite = L(parametre(i,5):parametre(i,6),parametre(i,3):parametre(i,4)); label = parametre(i,1); for m = 1 : size(boite,1) for n = 1 : size(boite,2) if boite(m,n)==label moment_0_0(i,1)=moment_0_0(i,1)+ (m^p*n^q); end end end end % affichage centre de gravité %subplot(223) figure('color',[0.3137 0.3137 0.5098]) imshow(img) hold on for i = 1 : nbobj subplot(3,4,i) imshow(img(parametre(i,5):parametre(i,6),parametre(i,3):parametre(i,4))) title(strcat('A(R)=',num2str(moment_0_0(i)))) end

An example

Gravity center

The gravity center is computed with moments of order p=1 and q=0 for the coordinate on x and p=0 and q=1 for y. The moments found are then normalized by the moment of order 0 of the region (A(R)). we have:
$x_{c e n t r e} = \frac{1}{A (R)} \cdot \sum_{x = 0}^{M} \sum_{y + 0}^{N} I (x, y) \cdot x^{1} y^{0} = \frac{m_{10} (R)}{m_{00} (R)}$
$y_{c e n t r e} = \frac{1}{A (R)} \cdot \sum_{x = 0}^{M} \sum_{y + 0}^{N} I (x, y) \cdot x^{0} y^{1} = \frac{m_{01} (R)}{m_{00} (R)}$
Here is the program:

%centroide moment_1_0=zeros(nbobj,1); p=1; q=0; for i = 1 : nbobj boite = L(parametre(i,5):parametre(i,6),parametre(i,3):parametre(i,4)); label = parametre(i,1); for m = 1 : size(boite,1) for n = 1 : size(boite,2) if boite(m,n)==label moment_1_0(i,1)=moment_1_0(i,1)+ (m^p*n^q); end end end end moment_1_0=round(moment_1_0./moment_0_0); moment_1_0= moment_1_0+parametre(:,5); % moment_0_1=zeros(nbobj,1); p=0; q=1; for i = 1 : nbobj boite = L(parametre(i,5):parametre(i,6),parametre(i,3):parametre(i,4)); label = parametre(i,1); for m = 1 : size(boite,1) for n = 1 : size(boite,2) if boite(m,n)==label moment_0_1(i,1)=moment_0_1(i,1)+ (m^p*n^q); end end end end moment_0_1=round(moment_0_1./moment_0_0); moment_0_1= moment_0_1+parametre(:,3); % affichage centre de gravité %subplot(223) figure('color',[0.3137 0.3137 0.5098]) imshow(img) hold on for i = 1 : nbobj plot(moment_0_1(i),moment_1_0(i),'r+','LineWidth',2) end title('center of gravity')

An example

Central Moments

They are computed for coordinates taken with respect to the coordinates of the center of gravity:
$μ_{p q} = \sum_{x = 0}^{M} \sum_{y + 0}^{N} I (x, y) \cdot {(x - x_{c e n t r e})}^{p} {(y - y_{c e n t r e})}^{q}$

Normalized Central Moments

To normalize a moment, we multiply the inverse of the moment of order 0 (the aera) raised to the power (p+q+2)/2. We find:
${\overline{μ}}_{p q} (R) = μ_{p q} \cdot {(\frac{1}{μ_{00} (R)})}^{(\frac{p + q + 2}{2})}$

Orientation of regions

The orientation of a region shows the axe passing trough its centroïde and by the major part of the region. The angle $θ$ represent the angle of rotation of the axe with respect to the horizontal. We have:
$θ (R) = \frac{1}{2} \cdot t a n^{- 1} (\frac{2 \cdot μ_{11} (R)}{μ_{20} (R) - μ_{02} (R)})$
Thanks to this angle , we can plot the line segment passing through the center of gravity of the region. The factor $λ$ allows to defined the length of the line segment.
$(\binom{x'}{y'}) = (\binom{x_{c} e n t r e}{y_{c} e n t r e}) + λ (\binom{\cos (θ (R))}{\sin (θ (R))})$
Here is the program:

% orientation theta=zeros(nbobj,1); moment_0_0=zeros(nbobj,1); moment_1_0=zeros(nbobj,1); moment_0_1=zeros(nbobj,1); moment_1_1=zeros(nbobj,1); moment_2_0=zeros(nbobj,1); moment_0_2=zeros(nbobj,1); for i = 1 : nbobj boite = L(parametre(i,5):parametre(i,6)+1,parametre(i,3):parametre(i,4)+1); label = parametre(i,1); for m = 1 : size(boite,1) for n = 1 : size(boite,2) if boite(m,n)==label p=0; q=0; moment_0_0(i,1)=moment_0_0(i,1)+ (m^p)*(n^q); p=1; q=0; moment_1_0(i,1)=moment_1_0(i,1)+ (m^p)*(n^q); p=0; q=1; moment_0_1(i,1)=moment_0_1(i,1)+ (m^p)*(n^q); end end end moment_1_0(i,1)=moment_1_0(i,1)/moment_0_0(i,1); moment_0_1(i,1)=moment_0_1(i,1)/moment_0_0(i,1); for m = 1 : size(boite,1) for n = 1 : size(boite,2) if boite(m,n)==label p=1; q=1; moment_1_1(i,1)=moment_1_1(i,1)+ ((m-moment_1_0(i,1))^p)*((n-moment_0_1(i,1))^q); p=2; q=0; moment_2_0(i,1)=moment_2_0(i,1)+ ((m-moment_1_0(i,1))^p)*((n-moment_0_1(i,1))^q); p=0; q=2; moment_0_2(i,1)=moment_0_2(i,1)+ ((m-moment_1_0(i,1))^p)*((n-moment_0_1(i,1))^q); end end end moment_1_1(i,1)=moment_1_1(i,1)*(1/moment_0_0(i,1))^2; moment_2_0(i,1)=moment_2_0(i,1)*(1/moment_0_0(i,1))^2; moment_0_2(i,1)=moment_0_2(i,1)*(1/moment_0_0(i,1))^2; a=2*moment_1_1(i,1); b=moment_2_0(i,1)-moment_0_2(i,1); theta(i)= atan(a/b)/2; if moment_2_0(i,1)<moment_0_2(i,1) theta(i)=theta(i)+pi/2; end end % affichage orientation des momentsd %subplot(223) figure('color',[0.3137 0.3137 0.5098]) a=axes x=0; y=0; lambda=50; imshow(img) hold on for i = 1 : nbobj x(1)=moment_1_0(i)+parametre(i,5); x(2)=(x(1) + cos(theta(i))*lambda); y(1)=moment_0_1(i)+parametre(i,3); y(2)=(y(1) + sin(theta(i))*lambda); x=round(x); y=round(y); plot(y(1),x(1),'r+','LineWidth',2) plot([y(1) y(2)],[x(1) x(2)] ,'r') end title('moment orientation')

An example

Eccentricity of regions

$E c c (R) = \frac{μ_{20} + μ_{02} + \sqrt{{(μ_{20} - μ_{02})}^{2} + 4 \cdot μ_{11}^{2}}}{μ_{20} + μ_{02} - \sqrt{{(μ_{20} - μ_{02})}^{2} + 4 \cdot μ_{11}^{2}}} = \frac{a_{1}}{a_{2}},$
$r a = {(\frac{2 a_{1}}{A (R)})}^{\frac{1}{2}}$
$r b = {(\frac{2 a_{2}}{A (R)})}^{\frac{1}{2}}$
where ra and rb are the ellipse parameters of the region. We can plot the ellipse thanks to ellipse equations:
$(\binom{x (t)}{y (t)}) = (\binom{x_{c} e n t r e}{y_{c} e n t r e}) + (\begin{matrix} \cos (θ) & - \sin (θ) \\ \sin (θ) & \cos (θ) \end{matrix}) \cdot (\binom{r a \cdot \cos (t)}{r b \cdot \sin (t)})$
where $0 \leq t \leq 2 π$
the program:

% excentricité for i = 1 : nbobj a1(i)=moment_2_0(i)+moment_0_2(i)+sqrt( (moment_2_0(i)-moment_0_2(i))^2 + 4*moment_1_1(i)^2); a2(i)=moment_2_0(i)+moment_0_2(i)-sqrt( (moment_2_0(i)-moment_0_2(i))^2 + 4*moment_1_1(i)^2); ra(i)=(2*a1(i)/moment_0_0(i))^0.5; rb(i)=(2*a2(i)/moment_0_0(i))^0.5; end figure('color',[0.3137 0.3137 0.5098]) a=axes x=0; y=0; imshow(img) hold on for i = 1 : nbobj a(1)=moment_1_0(i)+parametre(i,5); b(1)=moment_0_1(i)+parametre(i,3); a=round(a); b=round(b); ra(i)=ra(i)*moment_0_0(i); rb(i)=rb(i)*moment_0_0(i); t=0:1/1000:2*pi; for j = 1 : length(t) [x(1:2,j)]= [a;b]+[cos(theta(i)) -sin(theta(i));sin(theta(i)) cos(theta(i))]*[ra(i)*cos(t(j));rb(i)*sin(t(j))]; end plot(b,a,'r+','LineWidth',2) plot(x(2,:),x(1,:) ,'r') end title('moment orientation')

An example

Complete Program

clear all; close all; x = imread('objetinv2.PNG'); x=x(:,:,1); [M N] = size(x); %binarisation level = graythresh(x); img = im2bw(x,0.1); %labellisation CC = bwconncomp(img); L = labelmatrix(CC); %affihcage figure subplot(221) imshow(img) subplot(222) imshow(label2rgb(L)); %nombre d'objet, number of objects nbobj = max(max(L)); %parametres des objets, object's parameters param=zeros(nbobj,6); for i = 1 : nbobj param(i,1)=i; end %param = [ numero , nb pixel , boite_gauche , boite_droite , boite_haut , %boite_bas] for i = 1 : M for j = 1 : N if L(i,j)~=0 numero = L(i,j); param(numero,2)= param(numero,2) + 1 ; if param(numero,3)==0 || param(numero,3)>j param(numero,3)=j; end if param(numero,4)==0 || param(numero,4)<j param(numero,4)=j; end if param(numero,5)==0 || param(numero,5)>i param(numero,5)=i; end if param(numero,6)==0 || param(numero,6)<i param(numero,6)=i; end end end end % supression petit objet var =1; for i = 1 : nbobj if param(i,2)>10 parametre(var,:)=param(i,:); var=var+1; end end nbobj = size(parametre,1); %%%%%%moment moment_0_0=zeros(nbobj,1); p=0; q=0; for i = 1 : nbobj boite = L(parametre(i,5):parametre(i,6),parametre(i,3):parametre(i,4)); label = parametre(i,1); for m = 1 : size(boite,1) for n = 1 : size(boite,2) if boite(m,n)==label moment_0_0(i,1)=moment_0_0(i,1)+ (m^p*n^q); end end end end % affichage centre de gravité %subplot(223) figure('color',[0.3137 0.3137 0.5098]) imshow(img) hold on for i = 1 : nbobj subplot(3,4,i) imshow(img(parametre(i,5):parametre(i,6),parametre(i,3):parametre(i,4))) title(strcat('A(R)=',num2str(moment_0_0(i)))) end %centroide moment_1_0=zeros(nbobj,1); p=1; q=0; for i = 1 : nbobj boite = L(parametre(i,5):parametre(i,6),parametre(i,3):parametre(i,4)); label = parametre(i,1); for m = 1 : size(boite,1) for n = 1 : size(boite,2) if boite(m,n)==label moment_1_0(i,1)=moment_1_0(i,1)+ (m^p*n^q); end end end end moment_1_0=round(moment_1_0./moment_0_0); moment_1_0= moment_1_0+parametre(:,5); % moment_0_1=zeros(nbobj,1); p=0; q=1; for i = 1 : nbobj boite = L(parametre(i,5):parametre(i,6),parametre(i,3):parametre(i,4)); label = parametre(i,1); for m = 1 : size(boite,1) for n = 1 : size(boite,2) if boite(m,n)==label moment_0_1(i,1)=moment_0_1(i,1)+ (m^p*n^q); end end end end moment_0_1=round(moment_0_1./moment_0_0); moment_0_1= moment_0_1+parametre(:,3); % affichage centre de gravité %subplot(223) figure('color',[0.3137 0.3137 0.5098]) imshow(img) hold on for i = 1 : nbobj plot(moment_0_1(i),moment_1_0(i),'r+','LineWidth',2) end title('center of gravity') % orientation theta=zeros(nbobj,1); moment_0_0=zeros(nbobj,1); moment_1_0=zeros(nbobj,1); moment_0_1=zeros(nbobj,1); moment_1_1=zeros(nbobj,1); moment_2_0=zeros(nbobj,1); moment_0_2=zeros(nbobj,1); for i = 1 : nbobj boite = L(parametre(i,5):parametre(i,6)+1,parametre(i,3):parametre(i,4)+1); label = parametre(i,1); for m = 1 : size(boite,1) for n = 1 : size(boite,2) if boite(m,n)==label p=0; q=0; moment_0_0(i,1)=moment_0_0(i,1)+ (m^p)*(n^q); p=1; q=0; moment_1_0(i,1)=moment_1_0(i,1)+ (m^p)*(n^q); p=0; q=1; moment_0_1(i,1)=moment_0_1(i,1)+ (m^p)*(n^q); end end end moment_1_0(i,1)=moment_1_0(i,1)/moment_0_0(i,1); moment_0_1(i,1)=moment_0_1(i,1)/moment_0_0(i,1); for m = 1 : size(boite,1) for n = 1 : size(boite,2) if boite(m,n)==label p=1; q=1; moment_1_1(i,1)=moment_1_1(i,1)+ ((m-moment_1_0(i,1))^p)*((n-moment_0_1(i,1))^q); p=2; q=0; moment_2_0(i,1)=moment_2_0(i,1)+ ((m-moment_1_0(i,1))^p)*((n-moment_0_1(i,1))^q); p=0; q=2; moment_0_2(i,1)=moment_0_2(i,1)+ ((m-moment_1_0(i,1))^p)*((n-moment_0_1(i,1))^q); end end end moment_1_1(i,1)=moment_1_1(i,1)*(1/moment_0_0(i,1))^2; moment_2_0(i,1)=moment_2_0(i,1)*(1/moment_0_0(i,1))^2; moment_0_2(i,1)=moment_0_2(i,1)*(1/moment_0_0(i,1))^2; a=2*moment_1_1(i,1); b=moment_2_0(i,1)-moment_0_2(i,1); theta(i)= atan(a/b)/2; if moment_2_0(i,1)<moment_0_2(i,1) theta(i)=theta(i)+pi/2; end end % affichage orientation des momentsd %subplot(223) figure('color',[0.3137 0.3137 0.5098]) a=axes x=0; y=0; lambda=50; imshow(img) hold on for i = 1 : nbobj x(1)=moment_1_0(i)+parametre(i,5); x(2)=(x(1) + cos(theta(i))*lambda); y(1)=moment_0_1(i)+parametre(i,3); y(2)=(y(1) + sin(theta(i))*lambda); x=round(x); y=round(y); plot(y(1),x(1),'r+','LineWidth',2) plot([y(1) y(2)],[x(1) x(2)] ,'r') end title('moment orientation') % excentricité for i = 1 : nbobj a1(i)=moment_2_0(i)+moment_0_2(i)+sqrt( (moment_2_0(i)-moment_0_2(i))^2 + 4*moment_1_1(i)^2); a2(i)=moment_2_0(i)+moment_0_2(i)-sqrt( (moment_2_0(i)-moment_0_2(i))^2 + 4*moment_1_1(i)^2); ra(i)=(2*a1(i)/moment_0_0(i))^0.5; rb(i)=(2*a2(i)/moment_0_0(i))^0.5; end figure('color',[0.3137 0.3137 0.5098]) x=0; y=0; imshow(img) hold on for i = 1 : nbobj a(1)=moment_1_0(i)+parametre(i,5); b(1)=moment_0_1(i)+parametre(i,3); a=round(a); b=round(b); ra(i)=ra(i)*moment_0_0(i); rb(i)=rb(i)*moment_0_0(i); t=0:1/1000:2*pi; for j = 1 : length(t) [x(1:2,j)]= [a;b]+[cos(theta(i)) -sin(theta(i));sin(theta(i)) cos(theta(i))]*[ra(i)*cos(t(j));rb(i)*sin(t(j))]; end plot(b,a,'r+','LineWidth',2) plot(x(2,:),x(1,:) ,'r') end title('moment orientation')

Home

My videos

Summary

General information

Acquisition and Save

Gray Level Image Processing

Binary processing

Color Processing

Interest Points Detection

Time-Frequency processing

Multi Scale Processing

Optical Flow

Binary processing

Moments of Regions

Moments

Aera of regions

Gravity center

Central Moments

Normalized Central Moments

Orientation of regions

Eccentricity of regions

Complete Program