خانه / مجلات / مجله شماره 250– مهر 93 / پردازش تصویر در متلب

پردازش تصویر در متلب

این مقاله به آموزش استفاده از متلب در پردازش تصویر اختصاص دارد. آشنایی اولیه با متلب مورد نیاز است (باید نحوه استفاده از ماتریس ها و نحوه نوشتن یک M-file را بدانید).

بهتر است جعبه ابزار پردازش تصویر (Image Processing Toolbox) را داشته باشید ولی خوشبختانه برای بیشتر عملیات نیاز به هیچ جعبه ابزاری نیست. دستوراتی که نیازمند جعبه ابزار پردازش تصویر باشند با [Image Toolbox] علامت گذاری شده اند.

نمایش تصویر
۵ نوع تصویر در متلب وجود دارد.
۱- Grayscale: یک تصویر سیاه و سفید با طول M و عرض N پیکسل بصورت ماتریسی از دو دیتاتایپ دابل با اندازه MxN نمایش داده می شود. ارزش عناصر (برای مثال MyImage(m,n) ) نشان دهنده شدت grayscale پیکسل در [۰,۱] با =۰ سیاه و=۰ سفیدمی باشد.

۲- RGB: یک تصویر قرمز- سبز-آبی (Red-Green-Blue) بصورت یک ماتریس دابل سه بعدی MxNx3 نمایش داده می شود. هر پیکسل دارای اعضای قرمز، سبز و آبی به همراه بعد سوم است که دارای ارزش ۰ و ۱ است برای مثال اجزای رنگی پیکسل (m,n) به این ترتیب هستند:
MyImage(m,n,1) = red
MyImage(m,n,2) = green
MyImage(m,n,3) = blue
۳- Indexed: تصاویر Indexed (paletted) با یک ماتریس شاخص با سایز M×N و یک ماتریس colormap با سایز Kx3 نمایش داده می شود. Colormap همه رنگهای موجود در تصویر را در بر دارد و ماتریس شاخص نمایانگر پیکسل ها با اشاره به رنگ های موجود در colormap می باشد. برای مثال اگر رنگ بیست و دوم قرمز مایل به ارغوانی (magenta) باشد: MyColormap(22,:) = [1,0,1] پس MyImage(m,n) = 22 یک پیکسل با رنگ قرمز است.

۴- Binary: یک تصویر باینری با یک مارتیس منطقی MxN نمایش داده می شود که ارزش پیکسل ها ۱ (True) و ۰ (False) هستند.

۵- uint8: این نوع تصویر حافظه کمتری اشغال می کند و برخی عملیات نسبت به نوع دابل سریع تر انجام می شوند. برای سادگی در این آموزش uint8 بیشتر از این شرح داده نمی شود.

Grayscale معمولا فرمت مورد علاقه در پردازش تصویر است. در مواردی که به رنگ نیاز باشد، می توان یک تصویر رنگی RGB را تجزیه کرده و بصورت سه تصویر سیاه و سفید مجزا با آن رفتار کرد. تصاویر indexed برای بیشتر عملیات باید به سیاه و سفید یا RGB تبدیل شوند.
در ادامه چند تبدیل نمایش داده شده است. چند دستور نیازمند جعبه ابزار پردازش تصویر هستند و با [Image Toolbox] نمایش داده شده اند.

نمایش تصویر grayscale یا باینری
image(MyGray*255);
axis image
colormap(gray(256));

نمایش تصویر RGB (نمایش خطا در صورتی که یک عنصر خارج از [۰,۱] باشد)
image(MyRGB);
axis image
نمایش تصویر RGB
image(min(max(MyRGB,0),1));
axis image

نمایش تصویر indexed
image(MyIndexed);
axis image
colormap(MyColormap);

جداسازی کانال های یک تصویر RGB
MyRed = MyRGB(:,:,1);
MyGreen = MyRGB(:,:,2);
MyBlue = MyRGB(:,:,3);
برگرداندن کانالها به کنار هم
MyRGB = cat(3,MyRed,MyGreen,MyBlue);

تبدیل grayscale به RGB
MyRGB = cat(3,MyGray,MyGray,MyGray);

تبدیل RGB به grayscale با استفاده از میانگین ساده
MyGray = mean(MyRGB,3);
تبدیل RGB به grayscale با استفاده از NTSC weighting [Image Toolbox] MyGray = rgb2gray(MyRGB);

تبدیل RGB به grayscale با استفاده از NTSC weighting
MyGray = 0.299*MyRGB(:,:,1) + 0.587*MyRGB(:,:,2) + 0.114*MyRGB(:,:,3);

تبدیل indexed image بهRGB [Image Toolbox] MyRGB = ind2rgb(MyIndexed,MyColormap);
تبدیل indexed image بهRGB
MyRGB = reshape(cat(3,MyColormap(MyIndexed,1),MyColormap(MyIndexed,2),…
MyColormap(MyIndexed,3)),size(MyIndexed,1),size(MyIndexed,2),3);

تبدیل تصویر RGB به indexed با استفاده از K colors [Image Toolbox] [MyIndexed,MyColormap] = rgb2ind(MyRGB,K);

تبدیل باینری به grayscale
MyGray = double(MyBinary);
تبدیل grayscale به باینری
MyBinary = (MyGray > 0.5);

خواندن و نوشتن فایل های تصاویر
متلب می تواند با دستوات imread و imwrite تصاویر را بخواند و بنویسد. اگرچه تعدادی از فرمت های فایل ها پشتیبانی می شوند، برخی از آنها نیز پشتیبانی نمی شوند. از دستور imformats استفاده کنید تا ببینید که متلب شما از چه فرمت هایی پشتیبانی می کند:

>> imformats
EXT ISA INFO READ WRITE ALPHA DESCRIPTION
bmp isbmp imbmpinfo readbmp writebmp 0 Windows Bitmap (BMP)
gif isgif imgifinfo readgif writegif 0 Graphics Interchange Format (GIF)
jpg jpeg isjpg imjpginfo readjpg writejpg 0 Joint Photographic Experts Group (JPEG)
pbm ispbm impnminfo readpnm writepnm 0 Portable Bitmap (PBM)
pcx ispcx impcxinfo readpcx writepcx 0 Windows Paintbrush (PCX)
pgm ispgm impnminfo readpnm writepnm 0 Portable Graymap (PGM)
png ispng impnginfo readpng writepng 1 Portable Network Graphics (PNG)
pnm ispnm impnminfo readpnm writepnm 0 Portable Any Map (PNM)
ppm isppm impnminfo readpnm writepnm 0 Portable Pixmap (PPM)

هنگام خواندن تصاویر، یک مشکل قابل حل این است که دستور imread داده تصویر را با دیتاتایپ uint8 بر می گرداند که باید قبل از استفاده به دابل تبدیل شود. بنابراین به جای استفاده مستقیم از imread از تابع M-file زیر برای خواندن و تبدیل تصاویر استفاده می کنیم:
function Img = getimage(Filename)
%GETIMAGE Read an image given a filename
% V = GETIMAGE(FILENAME) where FILENAME is an image file. The image is
% returned either as an MxN double matrix for a grayscale image or as an
% MxNx3 double matrix for a color image, with elements in [0,1].

% Read the file
[Img,Map,Alpha] = imread(Filename);
Img = double(Img);

if ~isempty(Map) % Convert indexed image to RGB
Img = Img + 1;
Img = reshape(cat(3,Map(Img,1),Map(Img,2),Map(Img,3)),size(Img,1),size(Img,2),3);
else
Img = Img/255; % Rescale to [0,1] end
getimage.m را برای استفاده از این M-function ذخیره کنید. اگر تصویر baboon.png در دایرکتوری فعلی (یا جایی در مسیر جستجوی متلب) قرار دارد، می توانید آن را با دستور MyImage = getimage(‹baboon.png›) بخوانید. همچنین می توانید از مسیر های جزئی استفاده کنید، برای مثال اگر تصویر در /images/ قرار دارد با getimage(‹images/baboon.png›).
برای نوشتن یک تصویر grayscale یا RGB از imwrite(MyImage,›myimage.png›);
استفاده کنید. دقت کنید که MyImage یک ماتریس دابل با عناصر [۰,۱] است، اگر بطور نامناسبی اسکیل شود فایل ذخیره شده احتمالا خالی خواهد بود.
هنگام نوشتن فایل های تصویر به شدت توصیه می شود که از فرمت PNG استفاده کنید. از آنجایی که این فرمت از فشرده سازی lossless بهره می برد انتخاب قابل اعتمادی است، از truecolor RGB پشتیبانی کرده و به خوبی فشرده سازی می کند. از فرمت های دیگر با احتیاط استفاده کنید.

عملیات اصلی
در زیر برخی عملیات اصلی بر روی یک تصویر grayscale نمایش داده شده است. دستوراتی که نیازمند جعبه ابزار پردازش تصویر هستند با [Image Toolbox] مشخص شده اند.

Statistics
uMax = max(u(:)); محاسبه مقدار ماکزیمم
uMin = min(u(:)); مینیمم
uPower = sum(u(:).^2); توان
uAvg = mean(u(:)); میانگین
uVar = var(u(:)); واریانس
uMed = median(u(:)); مدین
hist(u(:),linspace(0,1,256)); پلات هیستوگرام

% Basic manipulations
uClip = min(max(u,0),1); % Clip elements to [0,1] uPad = u([1,1:end,end],[1,1:end,end]); % Pad image with one-pixel margin
uPad = padarray(u,[k,k],›replicate›); % Pad image with k-pixel margin [Image Toolbox] uCrop = u(RowStart:RowEnd,ColStart:ColEnd); % Crop image
uFlip = flipud(u); % Flip in the up/down direction
uFlip = fliplr(u); % Flip left/right
uResize = imresize(u,ScaleFactor); % Interpolate image [Image Toolbox] uRot = rot90(u,k); % Rotate by k*90 degrees with integer k
uRot = imrotate(u,Angle); % Rotate by Angle degrees [Image Toolbox] uc = (u – min(u(:))/(max(u(:)) – min(u(:))); % Stretch contrast to [0,1] uq = round(u*(K-1))/(K-1); % Quantize to K graylevels {0,1/K,2/K,…,1}

Simulating noise
uNoisy = u + randn(size(u))*sigma;
اضافه نمودن نویز گوسی سفید با انحراف معیار استاندارد سیگما
uNoisy = u; uNoisy(rand(size(u)) < p) = round(rand(size(u))); نویز نمک و فلفل Debugging any(~isfinite(u(:))) % Check if any elements are infinite or NaN nnz(u > 0.5)
شمردن اینکه چند المان در یک سری شروط برقرارند
(توجه: در هر آرایه سینتکس u(:) به معنی باز کردن u به یک بردار ستونی است. برای مثال اگر u = [1,5;0,2] آنگاه u(:) برابر است با [۱;۰;۵;۲].)
به عنوان مثال، قدرت سیگنال تصویر در محاسبه نسبت سیگنال به نویز (SNR) و پیک نسبت سیگنال به نویز (PSNR) استفاده می شود. اگر تصویر بدون نویز uclean و تصویر آلوده به نویز u باشد:
% Compute SNR
snr = -10*log10( sum((uclean(:) – u(:)).^2) / sum(uclean(:).^2) );

محاسبه PSNR: (که ماکزیمم مقدار ممکن برای uclean 1 است)
psnr = -10*log10( mean((uclean(:) – u(:)).^2) );
به norm دقت کنید: norm(v) در بردار v ، sqrt(sum(v.^2)) را محاسبه می کند ولی norm (A) در ماتریس A ، نرم L2 را محاسبه می کند،
norm(A) = sqrt(max(eig(A›*A)))
بنابراین norm (A) قطعا sqrt(sum(A(:).^2)) نیست. اشتباه متداول این است که در جایی که باید از norm(A(:)) استفاده کنید، norm (A) را به کار برید.

فیلترهای خطی
فیلتر خطی تکنیک اساسی پردازش سیگنال است. برای معرفی مختصر باید گفت فیلتر خطی عملیاتی است که بر روی تک تک پیکسل های xm,n یک تصویر انجام می شود یک تابع خطی بر روی پیکسل و همسایه های آن برای محاسبه یک مقدار جدید برای پیکسل (ym,n) اعمال می شود.
فیلتر خطی دو بعدی به فرم کلی زیر است:
ym,n = ∑j∑k hj,k xm−j,n−k

که درآن X ورودی، y خروجی و h پاسخ ضربه ای فیلتر است. انتخاب های متفاوت h منجر به فیلتر هایی می شود که مثلا لبه یابی کرده، تصویر را smooth یا sharp می کنند. سمت راست معادله بالا به طور مختصر با h x نمایش داده شده و به آن کانولوشن h و x می گویند.
فیلتر حوزه مکان
فیلتر خطی دو بعدی در MATLAB با conv2 پیاده سازی می شود. متاسفانه، conv2 تنها می تواند فیلتر را در نزدیکی مرزهای تصویر با zero-padding انجام دهد که به این معنی است که نتایج فیلتر معمولا برای پیکسل های نزدیک به مرز نامناسب است. برای اصلاح آن، می توانیم تصویر ورودی را لایه گذاری کرده و هنگام فراخوانی conv2 از گزینه ‹valid› استفاده کنیم. M-function زیر این کار را انجام می دهد:
function x = conv2padded(varargin)
%CONV2PADDED Two-dimensional convolution with padding.
% Y = CONV2PADDED(X,H) applies 2D filter H to X with constant
% extension padding.
%
% Y = CONV2PADDED(H1,H2,X) first applies 1D filter H1 along the rows
% and then applies 1D filter H2 along the columns.
%
% If X is a 3D array, filtering is done separately on each channel.
if nargin == 2 % Function was called as «conv2padded(x,h)»
x = varargin{1};
h = varargin{2};
top = ceil(size(h,1)/2)-1;
bottom = floor(size(h,1)/2);
left = ceil(size(h,2)/2)-1;
right = floor(size(h,2)/2);
elseif nargin == 3 % Function was called as «conv2padded(h1,h2,x)»
h1 = varargin{1};
h2 = varargin{2};
x = varargin{3};
top = ceil(length(h1)/2)-1;
bottom = floor(length(h1)/2);
left = ceil(length(h2)/2)-1;
right = floor(length(h2)/2);
else
error(‹Wrong number of arguments.›);
end
% Pad the input image
xPadded = x([ones(1,top),1:size(x,1),size(x,1)+zeros(1,bottom)],…
[ones(1,left),1:size(x,2),size(x,2)+zeros(1,right)],:);

% Since conv2 cannot handle 3D inputs, do filtering channel by channel
for p = 1:size(x,3)
if nargin == 2
x(:,:,p) = conv2(xPadded(:,:,p),h,›valid›); % Call conv2
else
x(:,:,p) = conv2(h1,h2,xPadded(:,:,p),›valid›); % Call conv2
end
end
conv2padded.m را ذخیره کنید تا بتوانید از این M-function استفاده نمایید. در ادامه چند مثال آورده شده است:
یک فیلتر smoothing سبک
h = [0,1,0;
۱,۴,۱;
۰,۱,۰];
h = h/sum(h(:)); % Normalize the filter
uSmooth = conv2padded(u,h);

% A sharpening filter
h = [0,-1,0;
-۱,۸,-۱;
۰,-۱,۰];
h = h/sum(h(:)); % Normalize the filter
uSharp = conv2padded(u,h);

لبه یاب سوبل
hx = [1,0,-1;
۲,۰,-۲;
۱,۰,-۱];
hy = rot90(hx,-1);
u_x = conv2padded(u,hx);
u_y = conv2padded(u,hy);
EdgeStrength = sqrt(u_x.^2 + u_y.^2);
% Moving average
WindowSize = 5;
h1 = ones(WindowSize,1)/WindowSize;
uSmooth = conv2padded(h1,h1,u);

% Gaussian filtering
sigma = 3.5;
FilterRadius = ceil(4*sigma);
گوسی با سیگمای ۴
h1 = exp(-(-FilterRadius:FilterRadius).^2/(2*sigma^2));
h1 = h1/sum(h1);
نرمال سازی فیلتر
uSmooth = conv2padded(h1,h1,u);
می گویند فیلتر دو بعدی h جدایی پذیر است اگر بتوان آن را بصورت ترکیب خارجی دو فیلتر یک بعدی h1 و h2 نوشت یعنی : h = h1(:)*h2(:).
محاسبه h1 و h2 سریع تر از h است، همان طور که در بالا برای فیلتر گوسی و فیلتر average انجام شده است. در حقیقت فیلتر های سوبل hx و hy نیز جدایی پذیر هستند.

فیلتر حوزه فوریه
فیلتر حوزه مکان با conv2 عملیاتی است که از نظر محاسباتی بسیار هزینه بر است. conv2 برای یک فیلتر K×K در یک تصویر M×N ، O(MNK2) عملیات جمع و ضرب انجام می دهد یا با فرض M-N-K، به اندازه O(N4) ، هزینه عملیاتی در بر دارد.
برای فیلتر های بزرگ، فیلتر در حوزه فوریه سریع تر است زیرا هزینه عملیاتی تا O(N2 log N) کاهش می یابد. با استفاده از خاصیت کانولوشن-ضرب تبدیل فوریه، کانولوشن هم ارز است با محاسبه:
% Compute y = h*x with periodic boundary extension
[k1,k2] = size(h);
hpad = zeros(size(x));
hpad([end+1-floor(k1/2):end,1:ceil(k1/2)], …
[end+1-floor(k2/2):end,1:ceil(k2/2)]) = h;
y = real(ifft2(fft2(hpad).*fft2(x)));
نتیجه معادل است با conv2padded(x,h) به جز در نزدیکی مرزها که محاسبه فوق از بسط دوره ای مرزها استفاده می کند.
فیلترینگ مبتنی بر فوریه را می توان با بسط متقارن مرزها انعکاس ورودی در هر جهت نیز انجام داد:
% Compute y = h*x with symmetric boundary extension
xSym = [x,fliplr(x)]; % Symmetrize horizontally
xSym = [xSym;flipud(xSym)]; % Symmetrize vertically
[k1,k2] = size(h);
hpad = zeros(size(xSym));
hpad([end+1-floor(k1/2):end,1:ceil(k1/2)], …
[end+1-floor(k2/2):end,1:ceil(k2/2)]) = h;
y = real(ifft2(fft2(hpad).*fft2(xSym)));
y = y(1:size(y,1)/2,1:size(y,2)/2);
(توجه: روشی که کارایی بهتری دارد فیلترینگ افزودن هم پوشانی FFT است. جعبه ابزار پردازش سیگنال افزودن هم پوشانی FFT را در یک بعد با دستور fftfilt انجام می دهد.

فیلترهای غیر خطی
فیلتر غیر خطی عملیاتی است که در آن هر پیکسل فیلتر شده ym,n یک تابع غیر خطی از xm,n و همسایگان آن است. در اینجا به طور خلاصه در مورد چند نوع از فیلترهای غیر خطی صحبت می شود.

Order statistic filters
اگر Image Toolbox را داشته باشید، فیلتر های order statistic با دستورات ordfilt2 یا medfilt2 انجام می شوند. این نوع فیلتر ها ارزشهای پیکسل را در یک همسایگی مرتب می کنند و kامین ارزش بزرگ را انتخاب می کنند. فیلتر های min، max و مدین برخی از انواع این نوع فیلترینگ هستند.

Morphological filters
اگر Image Toolbox را داشته باشید دستور bwmorph عملیات متعدد مورفولیژیکی بر روی تصاویر باینری مانند باز کردن و بستن انجام می دهند. به علاوه دستوراتی برای موفولوژی در تصاویر سیاه و سفید نیز وجود دارند: imerode، imdilate و imtophat.

فیلتر خود را بسازید
اغلب مواقع نیاز به استفاده از فیلتر جدیدی داریم که متلب آن را آماده ندارد. کد زیر نمونه از پیاده سازی فیلترها می باشد.
[M,N] = size(x);
y = zeros(size(x));
r = 1; % Adjust for desired window size
for n = 1+r:N-r
for m = 1+r:M-r
% Extract a window of size (2r+1)x(2r+1) around (m,n)
w = x(m+(-r:r),n+(-r:r));
% … write the filter here …
y(m,n) = result;
end
end
(توجه: حلقه ها در متلب در نسخه های قدیمی مانند نسخه ۵ به قبل بسیار کند بودند در نسخه های جدید سرعت اجرای حلقه ها افزایش پیدا کرده با این وجود این طور ترجیح داده می شود که تا حد ممکن کمتر از آنها استفاده گردد.)
به عنوان مثال alpha-trimmed mean filter کمترین و بیشترین مقدار d/2 در پنجره را نادیده می گیرد و از مقادیر باقیمانده یعنی (۲r+1)2−d میانگین می گیرد. این فیلتر موازنه ای میان فیلتر مدین و فیلتر متوسط گیر است. فیلتر alpha-trimmed mean در نمونه می تواند بصورت زیر پیاده سازی شود:
% The alpha-trimmed mean filter
w = sort(w(:));
y(m,n) = mean(w(1+d/2:end-d/2)); % Compute the result y(m,n)
به عنوان یک مثال دیگر فیلتر بایلترال بصورت زیر است:

ym,n = ∑j,k hj,k,m,n xm−j,n−k
——————————————–
∑j,k hj,k,m,n
که در آن
hj,k,m,n = e−(j2+k2)/(2σs2) e−(xm−j,n−k − xm,n)2/(2σd2).
این فیلتر را می توان بصورت زیر پیاده سازی کرد:
% The bilateral filter
[k,j] = meshgrid(-r:r,-r:r);
h = exp( -(j.^2 + k.^2)/(2*sigma_s^2) ) .* …
exp( -(w – w(r+1,r+1)).^2/(2*sigma_d^2) );
y(m,n) = h(:)›*w(:) / sum(h(:));
اگر Image Toolbox را ندارید از تمپلیت می توان به عنوان جایگزین برای فیلتر هایی که وجود ندارند استفاده کرد، با این وجود به سرعت پیاده سازی های Image Toolbox نمی رسند.
Filter Code
medfilt2 y(m,n) = median(w(:));
ordfilt2 w = sort(w(:));
y(m,n) = w(k); % Select the kth largest element
imdilate % Define a structure element as a (2r+1)x(2r+1) array
SE = [0,1,0;1,1,1;0,1,0];
y(m,n) = max(w(SE));
تا به اینجا عملیاتی را که معمولا در پردازش تصویر مفید هستند پوشش دادیم. در مقالات بعدی آموزش های بیشتری در زمینه پردازش تصویر در متلب ارائه خواهد شد.

2 دیدگاه

  1. با سلام
    از اطلاعات کاربردی و مفید شما متشکرم …

  2. سلام
    مطالبتون عالی وکاربردی بود .بی نهایت سپاسگزارم
    لطفابرنامه ای که نوشتید( فیلترحوزه مکانی)روش اجرا گرفتن ان را درمتلب۲۰۱۵و همچنین عکس که خوانده نمی شود راخواهشا توضیح دهید

دیدگاهتان را ثبت کنید

آدرس ایمیل شما منتشر نخواهد شدعلامتدارها لازمند *

*

x

شاید بپسندید

بهترین نرم افزارهای رایگان

TeamViewer 10 ابزار کنترل رایانه از راه دور www.snipca.com/14665 آن چه که ...