使用Matlab測量影象目標尺寸
在傳統的數字影象處理當中，
邊緣檢測
與
形態學
為兩門非常重要的技術，在筆者的第一篇文章中已經重點介紹了各種邊緣檢測運算元，因此這次筆者將結合一些較為簡單的形態學演算法，使用Matlab為大家介紹一個很有意思的測量目標尺寸的小專案，效果如下
圖1 效果圖
1。測距原理
在數字影象處理當中我們可知，在計算機眼中，每一張圖片都實際上表現為一個龐大的矩陣，若在不知道測距物體距離的情況下，是不可能對影象中物體進行大小（size）的測量計算的，因此我們需要引入一個和比例尺類似的概念：
pixels per metric ratio
意為給定度量單位的畫素比率，在本篇文章中我們將給點度量單位設定為英寸（inch），可以理解為參考物的作用`，給定影象中一參考物大小，便可測得其它目標物體的尺寸大小
其中，參考物需要有兩個重要性質：
性質1：
參考物尺寸
我們應該知道物體的尺寸（就是寬或高）包括測量的單位（如mm、英寸等）
性質2：
易於識別
我們應該能夠很容易地在圖片中找到參照物體，無論是基於物體的位置（例如，參考物體總是放在圖片的左上角）還是透過外觀（例如，獨特的顏色或形狀，不同與圖片中的其他物體）。無論是哪種情況，我們的參照物都應該以某種方式具有唯一的可識別性。
在本篇文章中，我們將硬幣作為我們的參考物，已知其尺寸大小為1 inch*1 inch，並且為滿足性質2，我們確保其始終置於
影象最左側
圖2 測距原理
因此我們得到
給定度量單位畫素比例
計算公式：
pixels per metric ratio = 硬幣畫素數/物體實際尺寸
已知硬幣長寬均為
1英寸
，假設其在影象中畫素寬為157px（基於其關聯的檢測框），得：
pixels per metric ratio = 157px/1.000in = 157px/in
透過使用這一比例，我們便可計算影象中其它物體的尺寸大小資訊了
2。
利用計算機視覺測量物體的大小
現在我們理解了
pixels per metric ratio
比率的含義，但我們還需要對目標進行檢測進行檢測框長寬的提取，這一步我們將用到諸如
灰度圖變換、邊緣檢測、形態學等演算法
首先我們定義硬幣長寬，並且讀取原始影象
coin_width
=
1。000
；
coin_height
=
1。000
；
I
=
imread
（
‘path\to\your\image’
）；
figure
（
1
），
imshow
（
I
）；
title
（
‘原始影象’
）；
圖3 原始影象
之後我們使用rgb2gray（image）函式進行
灰度圖轉換
，並且透過imfilter（f， w，boundary_options）函式對影象進行
高斯濾波
，其中w為濾波器，由fspecial（type，parameters，sigma）生成，其中將type設為“gaussian”，sigma設為1，程式碼及效果如下
%轉換為灰度影象
I1
=
rgb2gray
（
I
）；
%figure（2）；imshow（I1）；title（‘灰度影象’）；
sigma
=
1
；
gausFilter
=
fspecial
（
‘gaussian’
，［
5
5
］，
sigma
）；
I2
=
imfilter
（
I1
，
gausFilter
，
‘replicate’
）；
figure
（
2
）；
imshow
（
I2
）；
title
（
‘高斯濾波後圖像’
）；
圖4 高斯影象
對高斯濾波後的影象進行
Canny邊緣檢測
，使用edge（I，method，threshold），，其中I為輸入影象，method為指定演算法，文章中使用“canny”進行邊緣檢測，而threshold為閾值，通常設為0。1，具體想了解各邊緣檢測演算法詳解請看這篇文章
邊緣檢測程式碼及效果如下所示
I3
=
edge
（
I2
，
‘Canny’
，
0。1
）；
figure
（
3
），
imshow
（
I3
）；
title
（
‘ Canny邊緣檢測影象’
）；
圖5 Canny邊緣檢測影象
透過觀察邊緣影象可以發現各目標內部還具有較細的紋理，對後續八連通區域的檢測造成較大幹擾，因此我們可以透過
孔洞填充
操作去除內部紋理，其次提取孔洞填充影象外圍邊緣，最後使用形態學演算法去除較小物體，具體程式碼如下：
% 孔洞填充
I41
=
imfill
（
I3
，
‘holes’
）；
figure
（
4
），
imshow
（
I41
）；
title
（
‘孔洞填充影象’
）；
% 提取最外圍邊緣
I4
=
bwperim
（
I41
）；
figure
（
5
），
imshow
（
I4
）；
title
（
‘邊緣影象’
）；
% 去除面積小於150px物體
I5
=
bwareaopen
（
I4
，
150
）；
figure
（
6
），
imshow
（
I5
）；
其中bwperim函式與傳統的邊緣檢測概念不同，邊緣檢測為提取影象邊緣特徵影象，而bwperim只保留目標最外層邊緣，bwareaopen（I4，150）意為去除影象中小於150px的八連通區域（至於為什麼選擇150，emmmm筆者是把它當成引數直接人為調的哈）
最後效果如下：
圖6 孔洞填充影象
圖7 提取外圍邊緣
圖8 去除小物體
至此，我們便完成了對影象的處理階段，接下來我們需要對物體進行標記，並且繪製相關的檢測框圖，其中我們使用 ［labelpic，num］ = bwlabel（I5，8）來對影象
從左至右對目標進行標記
，並配合find函式進行特定目標的操作，如接下來我們需要計算
pixels per metric ratio
，我們需要將最左邊硬幣目標提取出來，我們可以使用以下程式碼進行提取
［
labelpic
，
num
］
=
bwlabel
（
I5
，
8
）；
% 標記
［
r
c
］=
find
（
labelpic
==
1
）；
其中r，c分別為目標物件（coin）的八連通行列座標，然後我們使用函式minboundrect（c，r，‘p’）將得到的r，c進行最小外接矩形的計算，其中引數p表示按邊長最小原則進行計算，最後使用minboxing（）函式計算最小外接矩形的長寬，具體程式碼如下：
% 獲取標記物體最小外接矩形座標點
［
rectx
，
recty
，
area
，
perimeter
］
=
minboundrect
（
c
，
r
，
‘p’
）；
% 計算最小外接矩形長寬
［
length
width
］
=
minboxing
（
rectx
（
1
：
end
-
1
），
recty
（
1
：
end
-
1
））；
計算
pixels per metric ratio:
% 單位英寸畫素點比例計算
pixels_length_rate
=
length
/
coin_height
；
pixels_width_rate
=
width
/
coin_width
；
其中minboundrect與minboxing函式並不是Matlab官方函式，具體程式碼會在文末給出
我們使用同樣的原理進行各目標最小外接矩形的計算，並使用line函式在影象中繪製出矩形框與中點連線
for
v
=
2
：
num
［
r
c
］=
find
（
labelpic
==
v
）；
［
rectx
，
recty
，
area
，
perimeter
］=
minboundrect
（
c
，
r
，
‘p’
）；
［
length
width
］
=
minboxing
（
rectx
（
1
：
end
-
1
），
recty
（
1
：
end
-
1
））；
% 繪製目標檢測框
line
（
rectx
，
recty
，
‘color’
，
‘y’
，
‘linewidth’
，
2
）；
midpointx
（
1
）=（
rectx
（
1
）
+
rectx
（
2
））
/
2
；
midpointx
（
2
）=（
rectx
（
3
）
+
rectx
（
4
））
/
2
；
midpointx
（
3
）=（
rectx
（
2
）
+
rectx
（
3
））
/
2
；
midpointx
（
4
）=（
rectx
（
4
）
+
rectx
（
1
））
/
2
；
midpointy
（
1
）=（
recty
（
1
）
+
recty
（
2
））
/
2
；
midpointy
（
2
）=（
recty
（
3
）
+
recty
（
4
））
/
2
；
midpointy
（
3
）=（
recty
（
2
）
+
recty
（
3
））
/
2
；
midpointy
（
4
）=（
recty
（
4
）
+
recty
（
1
））
/
2
；
line
（
midpointx
，
midpointy
，
‘color’
，
‘m’
，
‘linewidth’
，
2
）；
最後，我們根據
pixels per metric ratio
與最小外接矩形長寬計算目標實際尺寸，並透過text函式在影象中進行顯示，程式碼如下：
line
（
midpointx
，
midpointy
，
‘color’
，
‘m’
，
‘linewidth’
，
2
）；
target_float_length
=
length
/
pixels_length_rate
；
target_length
=
num2str
（
target_float_length
）；
target_float_width
=
width
/
pixels_width_rate
；
target_width
=
num2str
（
target_float_width
）；
if
（（
rectx
（
2
）
-
rectx
（
1
））
<
=（
recty
（
2
）
-
recty
（
1
）））
text
（
midpointx
（
1
），
midpointy
（
1
）
-
10
，
target_length
，
‘Color’
，
‘white’
）；
text
（
midpointx
（
3
）
+
10
，
midpointy
（
3
），
target_width
，
‘Color’
，
‘white’
）；
else
text
（
midpointx
（
1
），
midpointy
（
1
）
-
10
，
target_width
，
‘Color’
，
‘white’
）；
text
（
midpointx
（
3
）
+
10
，
midpointy
（
3
），
target_length
，
‘Color’
，
‘white’
）；
end
好了現在我們完成了所有工作，最終的效果圖便是這樣啦
圖9 效果圖
最後在放兩個動圖讓大家看看效果吧
圖10 效果演示
圖11 效果演示
如上圖所示，我們已經成功的計算出圖片中每個物體的尺寸
然而，並非所有的結果都是完美的，可能的原因：
拍攝的角度並非是一個完美的90°俯視。如果不是90°拍攝，物體的尺寸可能會出現扭曲
沒有使用相機內在和外在引數來校準。當無法確定這些引數時，照片很容易發生徑向和切向的透鏡變形。執行額外的校準步驟來找到這些引數可以消除圖片中的失真並得到更好的物體大小的近似值
參考
1。Measuring size of objects in an image with OpenCV
程式碼如下：
coin_width
=
1。000
；
coin_height
=
1。000
；
I
=
imread
（
‘path\to\your\image’
）；
figure
（
1
），
imshow
（
I
）；
title
（
‘原影象’
）；
%轉換為灰度影象
I1
=
rgb2gray
（
I
）；
%figure（2）；imshow（I1）；title（‘灰度影象’）；
sigma
=
1
；
gausFilter
=
fspecial
（
‘gaussian’
，［
5
5
］，
sigma
）；
I2
=
imfilter
（
I1
，
gausFilter
，
‘replicate’
）；
%figure（2）；imshow（I2）；title（‘高斯濾波後圖像’）；
I3
=
edge
（
I2
，
‘Canny’
，
0。1
）；
%figure（3），imshow（I3）；title（‘ Canny邊緣檢測影象’）；
% 孔洞填充
I41
=
imfill
（
I3
，
‘holes’
）；
%figure（4），imshow（I41）；title（‘孔洞填充影象’）；
% 提取最外圍邊緣
I4
=
bwperim
（
I41
）；
%figure（5），imshow（I4）； title（‘邊緣影象’）；
% 去除面積小於150px物體
I5
=
bwareaopen
（
I4
，
150
）；
%figure（6），imshow（I5）；
［
labelpic
，
num
］
=
bwlabel
（
I5
，
8
）；
% 標記
［
r
c
］=
find
（
labelpic
==
1
）；
% 獲取標記物體最小外接矩形座標點
［
rectx
，
recty
，
area
，
perimeter
］
=
minboundrect
（
c
，
r
，
‘p’
）；
% 計算最小外接矩形長寬
［
length
width
］
=
minboxing
（
rectx
（
1
：
end
-
1
），
recty
（
1
：
end
-
1
））；
%figure（7），imshow（I5）；%hold on
% 單位英寸畫素點比例計算
pixels_length_rate
=
length
/
coin_height
；
pixels_width_rate
=
width
/
coin_width
；
midpointx
=
zeros
（
2
，
2
）；
midpointy
=
zeros
（
2
，
2
）；
target_mat
=
zeros
（
1
，
num
）；
for
v
=
2
：
num
［
r
c
］=
find
（
labelpic
==
v
）；
［
rectx
，
recty
，
area
，
perimeter
］=
minboundrect
（
c
，
r
，
‘p’
）；
［
length
width
］
=
minboxing
（
rectx
（
1
：
end
-
1
），
recty
（
1
：
end
-
1
））；
% 繪製目標檢測框
line
（
rectx
，
recty
，
‘color’
，
‘y’
，
‘linewidth’
，
2
）；
midpointx
（
1
）=（
rectx
（
1
）
+
rectx
（
2
））
/
2
；
midpointx
（
2
）=（
rectx
（
3
）
+
rectx
（
4
））
/
2
；
midpointx
（
3
）=（
rectx
（
2
）
+
rectx
（
3
））
/
2
；
midpointx
（
4
）=（
rectx
（
4
）
+
rectx
（
1
））
/
2
；
midpointy
（
1
）=（
recty
（
1
）
+
recty
（
2
））
/
2
；
midpointy
（
2
）=（
recty
（
3
）
+
recty
（
4
））
/
2
；
midpointy
（
3
）=（
recty
（
2
）
+
recty
（
3
））
/
2
；
midpointy
（
4
）=（
recty
（
4
）
+
recty
（
1
））
/
2
；
% 繪製目標長寬中點間連線
line
（
midpointx
，
midpointy
，
‘color’
，
‘m’
，
‘linewidth’
，
2
）；
target_float_length
=
length
/
pixels_length_rate
；
target_length
=
num2str
（
target_float_length
）；
target_float_width
=
width
/
pixels_width_rate
；
target_width
=
num2str
（
target_float_width
）；
% 顯示目標物體長寬資訊
if
（（
rectx
（
2
）
-
rectx
（
1
））
<
=（
recty
（
2
）
-
recty
（
1
）））
text
（
midpointx
（
1
），
midpointy
（
1
）
-
10
，
target_length
，
‘Color’
，
‘white’
）；
text
（
midpointx
（
3
）
+
10
，
midpointy
（
3
），
target_width
，
‘Color’
，
‘white’
）；
else
text
（
midpointx
（
1
），
midpointy
（
1
）
-
10
，
target_width
，
‘Color’
，
‘white’
）；
text
（
midpointx
（
3
）
+
10
，
midpointy
（
3
），
target_length
，
‘Color’
，
‘white’
）；
end
end
minboundrect函式
function
［rectx，recty，area，perimeter］
=
minboundrect
（
x，y，metric
）
% minboundrect： Compute the minimal bounding rectangle of points in the plane
% usage： ［rectx，recty，area，perimeter］ = minboundrect（x，y，metric）
%
% arguments： （input）
% x，y - vectors of points， describing points in the plane as
% （x，y） pairs。 x and y must be the same lengths。
%
% metric - （OPTIONAL） - single letter character flag which
% denotes the use of minimal area or perimeter as the
% metric to be minimized。 metric may be either ‘a’ or ‘p’，
% capitalization is ignored。 Any other contraction of ‘area’
% or ‘perimeter’ is also accepted。
%
% DEFAULT： ‘a’ （‘area’）
%
% arguments： （output）
% rectx，recty - 5x1 vectors of points that define the minimal
% bounding rectangle。
%
% area - （scalar） area of the minimal rect itself。
%
% perimeter - （scalar） perimeter of the minimal rect as found
%
%
% Note： For those individuals who would prefer the rect with minimum
% perimeter or area， careful testing convinces me that the minimum area
% rect was generally also the minimum perimeter rect on most problems
% （with one class of exceptions）。 This same testing appeared to verify my
% assumption that the minimum area rect must always contain at least
% one edge of the convex hull。 The exception I refer to above is for
% problems when the convex hull is composed of only a few points，
% most likely exactly 3。 Here one may see differences between the
% two metrics。 My thanks to Roger Stafford for pointing out this
% class of counter-examples。
%
% Thanks are also due to Roger for pointing out a proof that the
% bounding rect must always contain an edge of the convex hull， in
% both the minimal perimeter and area cases。
%
%
% See also： minboundcircle， minboundtri， minboundsphere
%
%
% default for metric
if
（
nargin
<
3
）
||
isempty
（
metric
）
metric
=
‘a’
；
elseif
~
ischar
（
metric
）
error
‘metric must be a character flag if it is supplied。’
else
% check for ‘a’ or ‘p’
metric
=
lower
（
metric
（：）
‘
）；
ind
=
strmatch
（
metric
，{
’area‘
，
’perimeter‘
}）；
if
isempty
（
ind
）
error
’metric does not match either ‘’area‘’ or ‘’perimeter‘’‘
end
% just keep the first letter。
metric
=
metric
（
1
）；
end
% preprocess data
x
=
x
（：）；
y
=
y
（：）；
% not many error checks to worry about
n
=
length
（
x
）；
if
n
~=
length
（
y
）
error
’x and y must be the same sizes‘
end
% if var（x）==0
% start out with the convex hull of the points to
% reduce the problem dramatically。 Note that any
% points in the interior of the convex hull are
% never needed， so we drop them。
if
n
>
3
%%%%%%%%%%%%%%%%%%%%%%%%%
if
（
var
（
x
）
==
0
||
var
（
y
）
==
0
）
if
var
（
x
）
==
0
x
=
［
x
-
1
；
x
（
1
）；
x
+
1
］；
y
=
［
y
；
y
（
1
）；
y
］；
flag
=
1
；
else
y
=
［
y
-
1
；
y
（
1
）；
y
+
1
］；
x
=
［
x
；
x
（
1
）；
x
］；
flag
=
1
；
end
else
flag
=
0
；
%%%%%%%%%%%%%%%%%%%%%%
edges
=
convhull
（
x
，
y
）；
% ’Pp‘ will silence the warnings
end
% exclude those points inside the hull as not relevant
% also sorts the points into their convex hull as a
% closed polygon
%%%%%%%%%%%%%%%%%%%%
if
flag
==
0
%%%%%%%%%%%%%%%%%%%%
x
=
x
（
edges
）；
y
=
y
（
edges
）；
%%%%%%%%%%%%%%%%%%
end
%%%%%%%%%%%%%
% probably fewer points now， unless the points are fully convex
nedges
=
length
（
x
）
-
1
；
elseif
n
>
1
% n must be 2 or 3
nedges
=
n
；
x
（
end
+
1
）
=
x
（
1
）；
y
（
end
+
1
）
=
y
（
1
）；
else
% n must be 0 or 1
nedges
=
n
；
end
% now we must find the bounding rectangle of those
% that remain。
% special case small numbers of points。 If we trip any
% of these cases， then we are done， so return。
switch
nedges
case
0
% empty begets empty
rectx
=
［］；
recty
=
［］；
area
=
［］；
perimeter
=
［］；
return
case
1
% with one point， the rect is simple。
rectx
=
repmat
（
x
，
1
，
5
）；
recty
=
repmat
（
y
，
1
，
5
）；
area
=
0
；
perimeter
=
0
；
return
case
2
% only two points。 also simple。
rectx
=
x
（［
1
2
2
1
1
］）；
recty
=
y
（［
1
2
2
1
1
］）；
area
=
0
；
perimeter
=
2
*
sqrt
（
diff
（
x
）
。^
2
+
diff
（
y
）
。^
2
）；
return
end
% 3 or more points。
% will need a 2x2 rotation matrix through an angle theta
Rmat
=
@（
theta
）
［
cos
（
theta
）
sin
（
theta
）；
-
sin
（
theta
）
cos
（
theta
）］；
% get the angle of each edge of the hull polygon。
ind
=
1
：（
length
（
x
）
-
1
）；
edgeangles
=
atan2
（
y
（
ind
+
1
）
-
y
（
ind
），
x
（
ind
+
1
）
-
x
（
ind
））；
% move the angle into the first quadrant。
edgeangles
=
unique
（
mod
（
edgeangles
，
pi
/
2
））；
% now just check each edge of the hull
nang
=
length
（
edgeangles
）；
area
=
inf
；
perimeter
=
inf
；
met
=
inf
；
xy
=
［
x
，
y
］；
for
i
=
1
：
nang
% rotate the data through -theta
rot
=
Rmat
（
-
edgeangles
（
i
））；
xyr
=
xy
*
rot
；
xymin
=
min
（
xyr
，［］，
1
）；
xymax
=
max
（
xyr
，［］，
1
）；
% The area is simple， as is the perimeter
A_i
=
prod
（
xymax
-
xymin
）；
P_i
=
2
*
sum
（
xymax
-
xymin
）；
if
metric
==
’a‘
M_i
=
A_i
；
else
M_i
=
P_i
；
end
% new metric value for the current interval。 Is it better？
if
M_i
<
met
% keep this one
met
=
M_i
；
area
=
A_i
；
perimeter
=
P_i
；
rect
=
［
xymin
；［
xymax
（
1
），
xymin
（
2
）］；
xymax
；［
xymin
（
1
），
xymax
（
2
）］；
xymin
］；
rect
=
rect
*
rot
’
；
rectx
=
rect
（：，
1
）；
recty
=
rect
（：，
2
）；
end
end
% get the final rect
% all done
end
% mainline end
minboxing函式
function
［ wid hei ］
=
minboxing
（
d_x ， d_y
）
%minboxing Summary of this function goes here
% Detailed explanation goes here
dd
=
［
d_x
，
d_y
］；
dd1
=
dd
（［
4
1
2
3
］，：）；
ds
=
sqrt
（
sum
（（
dd
-
dd1
）
。^
2
，
2
））；
wid
=
min
（
ds
（
1
：
2
））；
hei
=
max
（
ds
（
1
：
2
））；
end
碼字不易，各位看官點個贊再走呀嘻嘻嘻
上一篇:《放鶴亭記》_東坡文選
下一篇：臺兒莊｜古城裡的戰爭與和平
