Устройство предварительной обработки изображений на ПЛИС

д., можно и с помощью свободно распространяемых инструментов - таких, как C-to-Verilog генератор. Кроме того, в OpenCV есть готовая реализация функции, вычисляющей контуры по алгоритму Канни.
Так же, есть возможность использовать интрументарий пакета MATLAB. В него входит HDL Coder - генератор некоторого произвольного описания аппаратуры для некоей гипотетической релаизации, или для конкретных ПЛИС от Xilinx, Altera и т.д. Так же, Image Processing Toolbox содержит удобный и обширный инструментраий и набор функций для цифровой обработки изображений, в том числе - включающий функцию нахождения границ, в том числе и по методу Канни.
Из этих двух вариантов мы выбираем Matlab, который в сочетании с расширениями Image Processing Toolbox и HDL Coder позволяет производить проектирование полностью интегрировано, в одной среде, и достаточно эффективно. Далее, предполагается что работа с полученным описанием аппаратуры будет производиться уже в некотором специализированном САПР. Одновременно, до этого будет выбрана ПЛИС и архитекутра устройства.
Функциональное описание алгоритма
В качестве исходных данных для проектирования у нас выступает выбранный алгоритм выделения контуров. Описание алгоритма производится в MATLAB на встроенном языке в виде так называемого m-файла. Используется библиотечная функция edge, реализующая оптимизированное вычисление границ, и тестовое окружение. Полные тексты программ приведены в приложении А.
Алгоритм состоит из следующих шагов:
применение ко всем элементам фильтра Гаусса, с окном 5х5 (входные данные - исходное изображение, выходные данные - отфильтрованное изображение в таком же формате);
вычисления для всех элементов производных посредством модифицированного оператора Собеля (входные данные - результаты предыдущего этапа, выходные данные - производные для каждого элемента);
вычисление градиента, модуля и угла (входные данные - производные для каждого элемента изображения, выходные данные - значения модуля и угла градиента для каждого элемента);
подавление не-максимумов градиента (входные данные - значения модуля и угла градиента для каждого элемента, выходные данные - значения максимумов градиентов);
двойная пороговая фильтрация (входные данные - значения максимумов градиентов, выходные данные - результирующий рисунок контуров).
Укрупнённое структурное описание алгоритма приведено на Рис. 1.

Рис. 1
После отладки программы и настройки фильтра Канни получаем следующий результат (Рис. 2). Исходное и отфильтрованное изображения приведены в приложении Б.
Рис. 2
Первые два этапа алгоритма, и последний, выполняются практически полностью независимо, на основе исходных данных предыдущих этапов, и могут реализовываться параллельно. Только вычисление градиента требует информации о состоянии соседних ячеек памяти изображения.
Разработка архитектуры устройства
Внешне устройство представляет собой вычислительное ядро с некоторыми высокопроизводительными интерфейсами ввода-вывода, позволяющими загружать и выгружать большие объёмы информации (Рис. 3).
Рис. 3
Пропускная способность вычислений (независимо от задержки) и скорость обработки информации должны совпадать. Допустим, нам нужно обрабатывать изображения размером 484х272 (просто пример, размер может быть выбран другой). Тогда одно изображение будет весить 1028Кбит, примерно один килобайт. Желательно стараться оргнизовывать вычисления и коммуникации так, чтобы они были одинаковыми по производительности. Однако, эти оценки лучше делать на этапе симуляции.
Так же, вычисления должны производиться по возможности в параллели, для оптимизации скорости выполнения. Поэтому, после генерации кода необходимо оценить, насколько это условие выполнятся, и если нет - произвести перенастройку генератора кода или же исправлять вручную.
Вычисления организуются в соответствии с алгоритмом, указанным на Рис. 1, отдельно для разных элементов. При этом, вычисления скользящей маской реализуются одновременными вычислениями для всех положений маски самостоятельными блоками.
Таким образом, на первом этапе для каждого элемента из 1053184 необходимо осуществить 24 операции умножения и операции сложения. Исходя из этого реализация фильтрации полностью параллельно является слишком трудоёмкой, и её желательно конвейеризовать.
Остальные операции менее трудоёмкие, и их можно выполнять в большей степени в параллели.
Создание HDL-кода
HDL-описание реализации алгоритма генерируем автоматически при помощи пакета HDL coder (Рис. 4).
Рис. 4
Генерация HDL кода для ПЛИС и ASIC происходит в несколько шагов:
создание проекта, добавление кода MATLAB;
оптимизация модели для достижения требуемых показателей занимаемой на кристалле площади и скорости;
генерация HDL кода;
проверка сгенерированного кода при помощи HDL Verifier.
В процессе происходит автоматическая конвертация кода MATLAB из арифметики с плавающей точкой в арифметику с фиксированной точкой и генерация синтезируемого VHDL или Verilog кода.
HDL Coder предоставляет библиотеки готовых к использованию логических элементов, таких как счетчики и таймеры, которые написаны на MATLAB.
Компромисс между занимаемой на кристалле площадью и скоростью работы достигается оптимизацией HDL код за счет распределенной конвееризации, потоковой обработки и совместного использования ресурсов.
В MATLAB имеется имеется возможность улучшенной оптимизацией циклов, потоковая обработка циклов и развертывание циклов, для проектов MATLAB содержащих циклы for или матричные операции. Постоянные массивы и матричные переменные в коде MATLAB можно расположить в блоки RAM.
После ряда итераций и улучшений получаем приемлемые результаты по занимаемым на кристалле ресурсам.
Выбор ПЛИС
После генерации алгоритма, мы получаем отчёт о необходимых для его реализации ресурсов (Рис. 4). В соответствии с этими данными, и полученным ранее обзором существующих ПЛИС (Таблица 1), и будет происходить выбор конкретного устройства для реализации алгоритма.
Рис. 4
Исходя из имеющихся данных, наиболее дешёвым, но всё же удовлетворяющим требованиям вариантом, видится ПЛИС производства Xilinx Spartan6 (Рис. 5).

Рис. 5

Выводы по разделу
В разделе, на основе полученных ранее данных, было произведено системное проектирование устройства на базе ПЛИС, реализующее алгоритм выделения границ по методу Канни.
Был выбран и обоснован выбор маршрута проектирования, заключающийся в следующих шагах: создание функционального последовательного описания алгоритма функционирования устройства, создание архитектурного описания устройства, генерация из функционального последовательного описания кода описания аппратуры, оценка реализуемости и выбор материальной базы для реализации полученного аппаратного описания.
Алгоритм был реализован и отлажен средствами Matlab с использованием библиотечных функций. Была произведена настройка и оценка различных параметров фильтра.
Была произведен генерация из последовательного описания на языке Matlab в hdl-описаие для реализации в ПЛИС. На основе данных об использовании ресурсов была выбрана оптимально подходящая по своим параметрам и наименее дорогая ПЛИС.

ЗАКЛЮЧЕНИЕ
В работе был проведён подробный анализ методов и средств реализации алгоритмов выделение контуров изображений. Было установлено, что имеется некоторый набор вариантов технической базы (ПК общего назначения с использованием CPU, использование GPU-сопроцессоров, использование ПЛИС, DSP, интегральных схем). Было установлено, что ПЛИС среди этих вариантов реализации позволяет осуществлять потоковую обработку информации, имеет высокую производительность, при этом может быть перепрограммирована.
Был проведён анализ методов проектирования для ПЛИС и осуществлён обзор рынка производителей. Были установлены основные маршруты проектирования, инструментарии и возможные устройства.
В ходе обзора методов и средств реализации было установлено, что часто используется нисходящий маршрут проектирования от функционального описания. Был произведён поиск имеющихся средств, позволяющих описыать алгоритмы и использовать их на ПЛИС. Для цифровой обработки это средства MATLAB, библиотеки OpenCV и OpenCL.
Был произведён патентный поиск и исследование предметной области, с целью определения перспективных и наиболее эффективных методов и средств обработки изображений.
Был подробно разобран алгоритм Канни, и произведено системное проектирование устройства, включающее в себя выбор архитектуры, описание алгоритма высокоуровневыми средствами, и расссмотрение их комбинирования в процессе проектирования. Было произведено моделирование в MATLAB, генерация кода в Verilog.

СПИСОК ЛИТЕРАТУРЫ1. Fykse E. Performance Comparison of GPU , DSP and FPGA implementations of image processing and computer vision algorithms in embedded systems Egil Fykse. 2013. №11 June 25.
2. Hussain W. Basics of field programmable Lecture, 2011. P. 1–11.
3. Brown S., Rose J. Architecture of FPGAs and CPLDs: A tutorial // IEEE Des. Test Comput. 1996. Vol. 13, № 2. P. 42–57.
4. Официальный сайт Xilinx. Ссылка: http://www.xilinx.com/.
5. Официальный сайт Altera. Ссылка: https://www.altera.com/.
6. Oruklu E. System-on-Chip Design Using High-Level Synthesis Tools // Circuits Syst. 2012. Vol. 03, № 01. P. 1–9.
7. Минь Н.К. Методы и алгоритмы обработки изображений в системах технического зрения промышленных робототехнических комплексов. Московский Государственный Технический Университет имени Н.Э. Баумана, 2015. 136 p.
8. Aresh K.N., Ahender M.M. FPGA Based Implementation of Image Edge Detection using Canny Edge Detection Algorithm // Int. J. Sci. Engeneering Techology Res. 2014. Vol. 03, № 29. P. 5840–5844.
9. Engineering T. Hardware Implementation of Edge Detection Algorithm-A Review 1,2 // Int. J. Eng. Dev. Res. 2014. Vol. 3, № 2. P. 55–56.
10. Intel Corporation B., OpCV. The OpenCV Tutorials 2.3. 2011. P. 351.
11. Munshi A. OpenCL 1.2 Specification // Version 1.2. 2012. P. 380.
12. Toolbox I.P. Image Processing Toolbox // Image Process. 2004. P. 1–4.
13. Евгений С. Способ помехоустойчивого градиентного выделения контуров объектов на цифровых изображениях: pat. RU 2 459 251 C2 USA. Россия: Александрович, Самойлин Евгений, 2009. P. 3.
14. Быков В., Новак В., Александр Р.А. Способ выделения локальных объектов на цифровых изображениях поверхности: pat. RU 2 459 251 C2 USA. Россия: Закрытое акционерное общество “Нанотехнология МДТ,” 2010. P. 5.
15. Гданский Николай, Юлия М. Градиентный способ выделения контуров объектов на матрице полутонового растрового изображения: pat. RU 2325044 C1 USA. Россия: Московский государственный университет инженерной экологии, 2007. P. 5.
16. Томакова Римма , Филист Сергей К.Н., Ольга Ш., Александр К. Способ автоматической сегментации полутоновых сложноструктурированных растровых изображений: pat. RU 2580074 C1 USA. Россия, 2016.
17. Мирошниченко, Сергей Труфанов М., Анциферов А. Способ выделения контуров объектов изображения и устройство для его реализации: pat. RU 2383925 C2 USA. Россия: Государственное образовательное учреждение высшего профессионального образования Курский государственный технический университет, 2007. P. 4.
18. Белим С.В., Кутлунин П.Е. Выделение контуров на изображениях с помощью алгоритма кластеризации. 2015. P. 119–124.
19. Жарких А.А. Двухэтапный алгоритм выделения контуров на изображении. 2007. P. 202–205.
20. Станислав Ч. « Синтез нейроподобных помехоустойчивых детекторов границ ». Санкт-Петербургский Государственный Университет Аэрокосмического Приборостроения, 2012. 121 p.
21. А.Власов И.Ц. Модификация алгоритма канни применительно к обработке рентгенографических изображений // Вестник науки Сибири. 201 3 . № 4. 2013. Vol. 4, № 10. P. 120–127.
22. Fei-Fie. Canny Edge Detection. 2015. P. 1–25.
Приложение А.
Тексты программ.
Файл edge.m - реализация функции выделения границ.
function [eout,thresh,gv_45,gh_135] = edge(varargin)

[a,method,thresh,sigma,thinning,H,kx,ky] = parse_inputs(varargin{:});

% Check that the user specified a valid number of output arguments
if ~any(strcmp(method,{'sobel','roberts','prewitt'})) && (nargout>2)
error(message('images:edge:tooManyOutputs'))
end

% Transform to a double precision intensity image if necessary
if ~isa(a,'double') && ~isa(a,'single')
a = im2single(a);
end

[m,n] = size(a);

if strcmp(method,'canny')
% Magic numbers
PercentOfPixelsNotEdges = .7; % Used for selecting thresholds
ThresholdRatio = .4; % Low thresh is this fraction of the high.

% Calculate gradients using a derivative of Gaussian filter
[dx, dy] = smoothGradient(a, sigma);

% Calculate Magnitude of Gradient
magGrad = hypot(dx, dy);

% Normalize for threshold selection
magmax = max(magGrad(:));
if magmax > 0
magGrad = magGrad / magmax;
end

% Determine Hysteresis Thresholds
[lowThresh, highThresh] = selectThresholds(thresh, magGrad, PercentOfPixelsNotEdges, ThresholdRatio, mfilename);

% Perform Non-Maximum Suppression Thining and Hysteresis Thresholding of Edge
% Strength
e = thinAndThreshold(dx, dy, magGrad, lowThresh, highThresh);
thresh = [lowThresh highThresh];

elseif strcmp(method,'canny_old')
% Magic numbers
GaussianDieOff = .0001;
PercentOfPixelsNotEdges = .7; % Used for selecting thresholds
ThresholdRatio = .4; % Low thresh is this fraction of the high.

% Design the filters - a gaussian and its derivative

pw = 1:30; % possible widths
ssq = sigma^2;
width = find(exp(-(pw.*pw)/(2*ssq))>GaussianDieOff,1,'last');
if isempty(width)
width = 1; % the user entered a really small sigma
end

t = (-width:width);
gau = exp(-(t.*t)/(2*ssq))/(2*pi*ssq); % the gaussian 1D filter

% Find the directional derivative of 2D Gaussian (along X-axis)
% Since the result is symmetric along X, we can get the derivative along
% Y-axis simply by transposing the result for X direction.
[x,y]=meshgrid(-width:width,-width:width);
dgau2D=-x.*exp(-(x.*x+y.*y)/(2*ssq))/(pi*ssq);

% Convolve the filters with the image in each direction
% The canny edge detector first requires convolution with
% 2D Gaussian, and then with the derivative of a Gaussian.
% Since Gaussian filter is separable, for smoothing, we can use
% two 1D convolutions in order to achieve the effect of convolving
% with 2D Gaussian. We convolve along rows and then columns.

%smooth the image out
aSmooth=imfilter(a,gau,'conv','replicate'); % run the filter across rows
aSmooth=imfilter(aSmooth,gau','conv','replicate'); % and then across columns

%apply directional derivatives
ax = imfilter(aSmooth, dgau2D, 'conv','replicate');
ay = imfilter(aSmooth, dgau2D', 'conv','replicate');

mag = sqrt((ax.*ax) + (ay.*ay));
magmax = max(mag(:));
if magmax>0
mag = mag / magmax; % normalize
end

% Select the thresholds
if isempty(thresh)
counts=imhist(mag, 64);
highThresh = find(cumsum(counts) > PercentOfPixelsNotEdges*m*n,...
1,'first') / 64;
lowThresh = ThresholdRatio*highThresh;
thresh = [lowThresh highThresh];
elseif length(thresh)==1
highThresh = thresh;
if thresh>=1
error(message('images:edge:thresholdMustBeLessThanOne'))
end
lowThresh = ThresholdRatio*thresh;
thresh = [lowThresh highThresh];
elseif length(thresh)==2
lowThresh = thresh(1);
highThresh = thresh(2);
if (lowThresh >= highThresh) || (highThresh >= 1)
error(message('images:edge:thresholdOutOfRange'))
end
end

% The next step is to do the non-maximum suppression. We will accrue
% indices which specify ON pixels in strong edgemap The array e will become
% the weak edge map.
e = cannyFindLocalMaxima(ax,ay,mag,lowThresh);

if ~isempty(e)
[rstrong,cstrong] = find(mag>highThresh & e);

if ~isempty(rstrong) % result is all zeros if idxStrong is empty
e = bwselect(e, cstrong, rstrong, 8);
e = bwmorph(e, 'thin', 1); % Thin double (or triple) pixel wide contours
end
end

elseif any(strcmp(method, {'log','zerocross'}))
% The output edge map:
e = false(m,n);

rr = 2:m-1; cc=2:n-1;

% We don't use image blocks here
if isempty(H),
fsize = ceil(sigma*3) * 2 + 1; % choose an odd fsize > 6*sigma;
op = fspecial('log',fsize,sigma);
else
op = H;
end

op = op - sum(op(:))/numel(op); % make the op to sum to zero
b = imfilter(a,op,'replicate');

if isempty(thresh)
thresh = .75*mean2(abs(b));
end

% Look for the zero crossings: +-, -+ and their transposes
% We arbitrarily choose the edge to be the negative point
[rx,cx] = find( b(rr,cc) < 0 & b(rr,cc+1) > 0 ...
& abs( b(rr,cc)-b(rr,cc+1) ) > thresh ); % [- +]
e((rx+1) + cx*m) = 1;
[rx,cx] = find( b(rr,cc-1) > 0 & b(rr,cc) < 0 ...
& abs( b(rr,cc-1)-b(rr,cc) ) > thresh ); % [+ -]
e((rx+1) + cx*m) = 1;
[rx,cx] = find( b(rr,cc) < 0 & b(rr+1,cc) > 0 ...
& abs( b(rr,cc)-b(rr+1,cc) ) > thresh); % [- +]'
e((rx+1) + cx*m) = 1;
[rx,cx] = find( b(rr-1,cc) > 0 & b(rr,cc) < 0 ...
& abs( b(rr-1,cc)-b(rr,cc) ) > thresh); % [+ -]'
e((rx+1) + cx*m) = 1;

% Most likely this covers all of the cases. Just check to see if there
% are any points where the LoG was precisely zero:
[rz,cz] = find( b(rr,cc)==0 );
if ~isempty(rz)
% Look for the zero crossings: +0-, -0+ and their transposes
% The edge lies on the Zero point
zero = (rz+1) + cz*m; % Linear index for zero points
zz = (b(zero-1) < 0 & b(zero+1) > 0 ...
& abs( b(zero-1)-b(zero+1) ) > 2*thresh); % [- 0 +]'
e(zero(zz)) = 1;
zz = (b(zero-1) > 0 & b(zero+1) < 0 ...
& abs( b(zero-1)-b(zero+1) ) > 2*thresh); % [+ 0 -]'
e(zero(zz)) = 1;
zz = (b(zero-m) < 0 & b(zero+m) > 0 ...
& abs( b(zero-m)-b(zero+m) ) > 2*thresh); % [- 0 +]
e(zero(zz)) = 1;
zz = (b(zero-m) > 0 & b(zero+m) < 0 ...
& abs( b(zero-m)-b(zero+m) ) > 2*thresh); % [+ 0 -]
e(zero(zz)) = 1;
end

else % one of the easy methods (roberts,sobel,prewitt)

if strcmp(method,'sobel')
op = fspecial('sobel')/8; % Sobel approximation to derivative
x_mask = op'; % gradient in the X direction
y_mask = op;

scale = 4; % for calculating the automatic threshold
offset = [0 0 0 0]; % offsets used in the computation of the threshold

elseif strcmp(method,'prewitt')
op = fspecial('prewitt')/6; % Prewitt approximation to derivative
x_mask = op';
y_mask = op;

scale = 4;
offset = [0 0 0 0];

elseif strcmp(method, 'roberts')
x_mask = [1 0; 0 -1]/2; % Roberts approximation to diagonal derivative
y_mask = [0 1;-1 0]/2;

scale = 6;
offset = [-1 1 1 -1];

else
error(message('images:edge:invalidEdgeDetectionMethod', method))
end

% compute the gradient in x and y direction
bx = imfilter(a,x_mask,'replicate');
by = imfilter(a,y_mask,'replicate');

if (nargout > 2) % if gradients are requested
gv_45 = bx;
gh_135 = by;
end

% compute the magnitude
b = kx*bx.*bx + ky*by.*by;

% determine the threshold; see page 514 of "Digital Imaging Processing" by
% William K. Pratt
if isempty(thresh), % Determine cutoff based on RMS estimate of noise
% Mean of the magnitude squared image is a
% value that's roughly proportional to SNR
cutoff = scale*mean2(b);
thresh = sqrt(cutoff);
else % Use relative tolerance specified by the user
cutoff = (thresh).^2;
end

if thinning
e = computeedge(b,bx,by,kx,ky,int8(offset),100*eps,cutoff);
else
e = b > cutoff;
end

end

if nargout==0,
imshow(e);
else
eout = e;
end

if isempty(a)
if nargout==2
if nargin == 2
if strcmp(method,'canny')
thresh = nan(1,2);
else
thresh = nan(1);
end
end
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Local Function : parse_inputs
%
function [I,Method,Thresh,Sigma,Thinning,H,kx,ky] = parse_inputs(varargin)
% OUTPUTS:
% I Image Data
% Method Edge detection method
% Thresh Threshold value
% Sigma standard deviation of Gaussian
% H Filter for Zero-crossing detection
% kx,ky From Directionality vector

narginchk(1,5)

I = varargin{1};

validateattributes(I,{'numeric','logical'},{'real','nonsparse','2d'},mfilename,'I',1);

% Defaults
Method = 'sobel';
Direction = 'both';
Thinning = true;

methods = {'canny','canny_old','prewitt','sobel','marr-hildreth','log','roberts','zerocross'};
directions = {'both','horizontal','vertical'};
options = {'thinning','nothinning'};

% Now parse the nargin-1 remaining input arguments

% First get the strings - we do this because the interpretation of the
% rest of the arguments will depend on the method.
nonstr = []; % ordered indices of non-string arguments
for i = 2:nargin
if ischar(varargin{i})
str = lower(varargin{i});
j = find(strcmp(str,methods));
k = find(strcmp(str,directions));
l = find(strcmp(str,options));
if ~isempty(j)
Method = methods{j(1)};
if strcmp(Method,'marr-hildreth')
error(message('images:removed:syntax','EDGE(I,''marr-hildreth'',...)','EDGE(I,''log'',...)'))
end
elseif ~isempty(k)
Direction = directions{k(1)};
elseif ~isempty(l)
if strcmp(options{l(1)},'thinning')
Thinning = true;
else
Thinning = false;
end
else
error(message('images:edge:invalidInputString', varargin{ i }))
end
else
nonstr = [nonstr i]; %#ok
end
end

% Now get the rest of the arguments
[Thresh,Sigma,H,kx,ky] = images.internal.parseNonStringInputsEdge(varargin,Method,Direction,nonstr);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Local Function : smoothGradient
%
function [GX, GY] = smoothGradient(I, sigma)

% Create an even-length 1-D separable Derivative of Gaussian filter

% Determine filter length
filterExtent = ceil(4*sigma);
x = -filterExtent:filterExtent;

% Create 1-D Gaussian Kernel
c = 1/(sqrt(2*pi)*sigma);
gaussKernel = c * exp(-(x.^2)/(2*sigma^2));

% Normalize to ensure kernel sums to one
gaussKernel = gaussKernel/sum(gaussKernel);

% Create 1-D Derivative of Gaussian Kernel
derivGaussKernel = gradient(gaussKernel);

% Normalize to ensure kernel sums to zero
negVals = derivGaussKernel < 0;
posVals = derivGaussKernel > 0;
derivGaussKernel(posVals) = derivGaussKernel(posVals)/sum(derivGaussKernel(posVals));
derivGaussKernel(negVals) = derivGaussKernel(negVals)/abs(sum(derivGaussKernel(negVals)));

% Compute smoothed numerical gradient of image I along x (horizontal)
% direction. GX corresponds to dG/dx, where G is the Gaussian Smoothed
% version of image I.
GX = imfilter(I, gaussKernel', 'conv', 'replicate');
GX = imfilter(GX, derivGaussKernel, 'conv', 'replicate');

% Compute smoothed numerical gradient of image I along y (vertical)
% direction. GY corresponds to dG/dy, where G is the Gaussian Smoothed
% version of image I.
GY = imfilter(I, gaussKernel, 'conv', 'replicate');
GY = imfilter(GY, derivGaussKernel', 'conv', 'replicate');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Local Function : selectThresholds
%
function [lowThresh, highThresh] = selectThresholds(thresh, magGrad, PercentOfPixelsNotEdges, ThresholdRatio, ~)

[m,n] = size(magGrad);

% Select the thresholds
if isempty(thresh)
counts=imhist(magGrad, 64);
highThresh = find(cumsum(counts) > PercentOfPixelsNotEdges*m*n,...
1,'first') / 64;
lowThresh = ThresholdRatio*highThresh;
elseif length(thresh)==1
highThresh = thresh;
if thresh>=1
error(message('images:edge:thresholdMustBeLessThanOne'))
end
lowThresh = ThresholdRatio*thresh;
elseif length(thresh)==2
lowThresh = thresh(1);
highThresh = thresh(2);
if (lowThresh >= highThresh) || (highThresh >= 1)
error(message('images:edge:thresholdOutOfRange'))
end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Local Function : thinAndThreshold
%
function H = thinAndThreshold(dx, dy, magGrad, lowThresh, highThresh)
% Perform Non-Maximum Suppression Thining and Hysteresis Thresholding of
% Edge Strength

% We will accrue indices which specify ON pixels in strong edgemap
% The array e will become the weak edge map.

E = cannyFindLocalMaxima(dx,dy,magGrad,lowThresh);

if ~isempty(E)
[rstrong,cstrong] = find(magGrad>highThresh & E);

if ~isempty(rstrong) % result is all zeros if idxStrong is empty
H = bwselect(E, cstrong, rstrong, 8);
else
H = false(size(E));
end
else
H = false(size(E));
end


Файл edge_detection.m - высокоуровневое представление функции, в котором задаются коэффициенты и вызывается основнfя функция edge.
function [BW2] = edge_detection(I1)
TR = [0.15 0.2];
BW2 = edge(I1,'Canny',TR);

Файл edge_detection_launch.m - скрипт, осуществляющий запуск функции и демонстрацию её результатов.
I = imread('армата_s.tif');
I1 = I(1:272,1:484);
TR = [0.15 0.2];
BW2 = edge(I1,'Canny',TR);
%BW2 = edge_detection(I1);
imshowpair(I1,BW2,'montage')



Приложение Б.
Исходное и обработанное изображения.












46