Гранулометрия

Обработка изображений применяется во многих сферах деятельности человека. Развитие информационных технологий способствует повышению качественного уровня анализа данных. Также прогресс компьютерной техники влияет на быстродействие и достоверность такой обработки. Приведем несколько наглядных примеров из тех областей, где применяется обработка изображений и, в частности, гранулометрия.

В первую очередь, это медицина. Гранулометрия в медицине применяется для подсчета количества определенных объектов на изображениях срезов биологических тканей. Например, позволяет автоматически выделять различные объекты изображения и оценивать их процентный состав в зависимости от площади; также применяется при анализе крови (построении и анализе эритроцитарной гистограммы, подсчете тромбоцитов и т.д.). Это может быть также анализ снимков микроскопа (подсчет бактерий, анализ хромосом и т.д.). Похожая область - это анализ фармацевтических препаратов при их изготовлении.

Вторая задача, где применяется гранулометрия - это анализ состава материалов в промышленности. Это может быть количественный анализ зерен в зернистых структурах (сталях, сплавах и т.п.) и получение их распределения по размерам и форме; оценка параметров включений и дефектов (размеры, форма, распределение). Причем все эти измерения, в большинстве случаев производятся автоматически, т.е. частицы сначала выделяются, измеряются и классифицируются. Ниже приведены примеры двух изображений, которые могут использоваться при гранулометрических исследованиях.

Под термином гранулометрия подразумевается измерение и классификация объектов изображения по одному из имеющихся параметров, например, размеру. Рассмотрим более детально решение этой задачи на конкретном примере. Результатом обработки является гистограмма, показывающая процентное содержание объектов каждого класса. Содержание

• Считывание изображения.

• Улучшение контрастности изображения.

• Определение яркости поверхности на улучшенном изображении.

• Вычисление первой производной распределения.

• Формирование объектов (снежинок) с учетом вычисленного радиуса.

Считывание изображения. Считываем файл 'snowflakes.png', который представляет собой изображение хлопьев снега. I=imread('snowflakes.png'); figure, imshow(I)

Улучшение контрастности изображения. Эффективность решения любой задачи по обработке изображений, в большой мере, зависит от качества исходных данных. В данном случае, когда необходимо идентифицировать объекты и анализировать их размер, то высокое качество исходного изображения является одним из необходимых условий. Для решения задачи повышения контрастности изображения в системе Matlab можно использовать функцию ADAPTHISTEQ. Эта функция реализует контрастно-ограниченную адаптивную эквализацию (выравнивание) гистограммы. Для коррекции яркостей изображения используется функция IMADJUST, которая позволяет управлять диапазоном яркостей изображения. claheI=adapthisteq(I, 'Divisions', [10 10]); claheI=imadjust(claheI); imshow(claheI);

Определение яркости поверхности на улучшенном изображении. При решении задач гранулометрии результат измерения интенсивностей элементов исследуемого изображения представляется в виде функции. Эта функция привязана к размерам исходного изображения. Каждому объекту изображения присваиваются определенные параметры, основной из которых - размер. На основании этих параметров производится построение гистограммы распределения объектов по определенных признаках. Вторым параметром, на основании которого осуществляется распределение, является интенсивность (яркость) объектов изображения. for counter=0:22 remain=imopen(claheI, strel('disk', counter)); intensity_area(counter+1)=sum(remain(:)); end figure,plot(intensity_area, 'm - *'), grid on; title('Функция зависимости суммы пикселей раскрытого изображения от их радиуса'); xlabel('радиус раскрытия (в пикселях)'); ylabel('значение суммы пикселей раскрытых объектов (интенсивность)');

Вычисление первой производной распределения. Распределение объектов изображения, полученное в результате гранулометрической обработки, характеризуется также рядом параметров. Одним из них является первая производная, которая вычисляется при помощи функции DIFF. intensity_area_prime=diff(intensity_area); plot(intensity_area_prime, 'm-*'), grid on; title('Granulometry (Size Distribution) of Snowflakes'); set(gca, 'xtick', [0 2 4 6 8 10 12 14 16 18 20 22]); xlabel('радиус снежинок (в пикселях)'); ylabel('зависимость суммы пикселей снежинок от радиуса');

Формирование объектов (снежинок) с учетом вычисленного радиуса Отметим минимум, который встречается на полученном графе. Этот минимум характеризует наиболее мелкие объекты (снежинки) на изображении. Более отрицательная минимальная точка свидетельствует о повышении кумулятивной интенсивности снежинок в зависимости от их радиуса. open5=imopen(claheI, strel('disk', 5)); open6=imopen(claheI, strel('disk', 6)); rad5=imsubtract(open5, open6); imshow(rad5, []);

Идентификация округлых предметов

Рассмотрим задачу, основной целью которой будет идентификация округлых объектов с помощью функции bwboundaries.

Содержание

• Шаг 1: Считывание изображения.

• Шаг 2: Пороговая обработка изображения.

• Шаг 3: Устранение шума.

• Шаг 4: Поиск границ объектов изображения.

• Шаг 5: Определение округлости объектов.

Шаг 1: Считывание изображения.

Считаем данные из pills_etc.png.

RGB=imread('pillsetc.png');

imshow(RGB);

Шаг 2: Пороговая обработка изображения.

Преобразуем исследуемое изображение в бинарное и таким образом подготовим его для применения

функции bwboundaries, которая реализует выделение границ объектов.

I=rgb2gray(RGB);

threshold=graythresh(I);

bw=im2bw(I, threshold);

imshow(bw)

Шаг 3: Устранение шума.

С использованием морфологических функций устраним объекты, которые не являются объектами нашего интереса, т.е. шум.

% удаление всех объектов, содержащих меньше чем 30 пикселей

bw=bwareaopen(bw,30);

% заполнение пустот

se=strel('disk', 2);

bw=imclose(bw, se);

bw=imfill(bw,'holes');

imshow(bw)

Шаг 4: Поиск границ объектов изображения.

При решении этого вопроса будем рассматривать только внешние границы.

Опция 'noholes' приводит к ускорению обработки с помощью функции bwboundaries

и повышает достоверность выделения границ объектов изображения.

[B, L]=bwboundaries(bw, 'noholes');

% отображение матрицы меток и извлечение границ

imshow(label2rgb(L, @jet, [.5 .5 .5]))

hold on

for k=1:length(B)

boundary=B{k};

plot(boundary(:, 2), boundary(:, 1), 'w', 'LineWidth', 2)

end

Шаг 5: Определение округлости объектов. Измерим площадь каждого объекта и его периметр. Эти данные будут использованы при вычислении характеристики округлостей объектов: metric=4*pi*area/perimeter^2. Параметр metric является достоверным только в случае идеальной окружности и является меньшим при других формах исследуемого объекта. Кроме того, при решении задач подобного рода, можно контролировать установку порога бинаризации, что приведет к улучшению качества обработки. При обработке исследуемых изображений использовался порог, равный 0.94. Используем функцию regionprops для получения измерений площади для всех объектов. Отметим, что матрица меток, возвращаемая функцией bwboundaries, может повторно использоваться в функции regionprops.

stats=regionprops(L, 'Area', 'Centroid');

threshold=0.94;

% окружность границ

for k=1:length(B)

% получение координат границ (X, Y), соответствующих метке 'k'

boundary=B{k};

% вычисление измерений на основе периметра объектов

delta_sq=diff(boundary).^2;

perimeter=sum(sqrt(sum(delta_sq, 2)));

% получение вычисленной площади, соответствующей метке 'k'

area=stats(k).Area;

% вычисление характеристики округлости metric

metric=4*pi*area/perimeter^2;

% отображение результатов

metric_string=sprintf('%2.2f', metric);

% маркирование объектов

if metric>threshold

centroid=stats(k).Centroid;

plot(centroid(1), centroid(2), 'ko');

end

text(boundary(1, 2)-35, boundary(1, 1)+13, metric_string,'Color', 'y',...

'FontSize', 14, 'FontWeight', 'bold');

end

title(['Metrics closer to 1 indicate that ',...

'the object is approximately round']);

Измерение углов пересечения

Часто в машинном зрении применяются автоматические измерения различных геометрических параметров с применением технологии захвата и обработки изображений. Измерение угла и определение точки пересечения двух элементов некоторой конструкции выполняется с использованием функции bwtraceboundary, которая реализует также автоматическое отслеживание границ.

Содержание

• Шаг 1: Считывание изображения.

• Шаг 2: Определение области интереса.

• Шаг 3: Пороговая обработка изображения.

• Шаг 4: Поиск начальных точек каждой границы.

• Шаг 5: Отслеживание границ.

• Шаг 6: Подгонка линий границ.

• Шаг 7: Определение угла пересечения.

• Шаг 8: Поиск точки пересечения.

• Шаг 9: Отображение результатов.

Шаг 1: Считывание изображения.

Считаем изображение из графического файла gantrycrane.jpg и определим область интереса

при решении этой задачи. Это изображение портала подъемного крана.

RGB=imread('gantrycrane.png');

imshow(RGB);

text(size(RGB, 2), size(RGB, 1)+15, 'Image courtesy of Jeff Mather', ...

'FontSize', 7, 'HorizontalAlignment', 'right');

line([300 328], [85 103], 'color', [1 1 0]);

line([268 255], [85 140], 'color', [1 1 0]);

text(150, 72, 'Измерение углов', 'Color', 'y', ...

'FontWeight', 'bold');

Шаг 2: Определение области интереса. Вырежем изображение, содержащее только те элементы конструкции крана, угол между которыми мы хотим определить. Далее необходимо определить контур двух металлических конструкций.

% Координаты области интереса можно задать вручную,

% далее отображаем нужную область с помощью imview

start_row=34;

start_col=208;

cropRGB=RGB(start_row:163, start_col:400, :);

imshow(cropRGB)

offsetX=start_col-1;

offsetY=start_row-1;

Шаг 3: Пороговая обработка изображения. Преобразуем исследуемое изображение в бинарное. Далее эти бинарные данные будут использованы для получения координат границ с помощью функции bwtraceboundary.

I=rgb2gray(cropRGB);

threshold=graythresh(I);

BW=im2bw(I, threshold);

BW=~BW; % дополнение изображения (объект интереса должен быть белым)

imshow(BW)

Шаг 4: Поиск начальных точек каждой границы.

Функция bwtraceboundary требует описания определенных точек контура.

Эти точки являются стартовыми при отслеживании границ объектов.

Функция bwtraceboundary также требует описание исследуемого изображения,

представленное в специальной форме.

dim=size(BW);

% горизонтальное направление

col1=4;

row1=min(find(BW(:, col1)));

% угловое направление

row2=12;

col2=min(find(BW(row2, :)));

Шаг 5: Отслеживание границ.

Функция bwtraceboundary позволяет реализовать также процедуру определения координат граничных точек. Для увеличения достоверности вычисленных углов и точек при расчетах нужно использовать максимальное число исследуемых граничных точек. Это число определяется экспериментально. При проведении процедуры определения координат границ также необходимо указывать направление отслеживания.

boundary1=bwtraceboundary(BW, [row1, col1], 'N', 8, 70);

% установим направление отслеживания - против часовой стрелки, сверху вниз

boundary2=bwtraceboundary(BW, [row2, col2], 'E', 8, 90,'counter');

imshow(RGB); hold on;

% отображение результатов на исходном изображении

plot(offsetX+boundary1(:, 2), offsetY+boundary1(:, 1), 'g', 'LineWidth', 2);

plot(offsetX+boundary2(:, 2), offsetY+boundary2(:, 1), 'g', 'LineWidth', 2);

Шаг 6: Подгонка линий границ. Из набора пар координат (X,Y), полученных на прежнем шаге, нужно отбросить заведомо ложные, которые не принадлежат искомой линии. Для проведения этой процедуры каждой точке контура необходимо на основании некоторого критерия присвоить соответствующий весовой коэффициент. Уравнение линии y=[x 1]*[a; b]. Параметры 'a' и 'b' определяются на основании критерия минимальной квадратической ошибки с использованием функции polyfit.

ab1=polyfit(boundary1(:, 2), boundary1(:, 1), 1);

ab2=polyfit(boundary2(:, 2), boundary2(:, 1), 1);

Шаг 7: Определение угла пересечения.

Используем функцию dot для поиска угла.

vect1=[1 ab1(1)]; % создание вектора на основании линейных выражений

vect2=[1 ab2(1)];

dp=dot(vect1, vect2);

% вычисление длины вектора

length1=sqrt(sum(vect1.^2));

length2=sqrt(sum(vect2.^2));

% вычисление угла

angle=180-acos(dp/(length1*length2))*180/pi

angle=

129.5860

Шаг 8: Поиск точки пересечения.

Решение системы двух выражений в порядке получения координат (X,Y) точек пересечения.

intersection=[1 ,-ab1(1); 1, -ab2(1)] \ [ab1(2); ab2(2)];

% получение результата на исходном,

% т.е. не вырезанном изображении.

intersection=intersection+[offsetY; offsetX]

intersection=

143.0453

295.7910

Шаг 9: Отображение результатов.

inter_x=intersection(2);

inter_y=intersection(1);

% отображение точек пересечения

plot(inter_x, inter_y, 'yx', 'LineWidth', 2);

text(inter_x-60, inter_y-30, [sprintf('%1.3f', angle), '{\circ}'],...

'Color', 'y', 'FontSize', 14, 'FontWeight', 'bold');

interString=sprintf('(%2.1f, %2.1f)', inter_x, inter_y);

text(inter_x-10, inter_y+20, interString,...

'Color', 'y', 'FontSize', 14, 'FontWeight', 'bold');

Измерение радиуса части мотка ленты

Рассмотрим задачу определения радиуса круглого объекта. В качестве примера возьмем катушку ленты. Решение задачи будет усложняться тем, что на изображении представлено только часть исследуемого объекта. При решении этой задачи будем использовать функцию bwtraceboundary.

Содержание

• Шаг 1: Считывание изображения.

• Шаг 2: Пороговая обработка изображения.

• Шаг 3: Определение положения начальной точки границы объекта.

• Шаг 4: Отслеживание границ.

• Шаг 5: Подгонка окружности на выделенные границы.

Шаг 1: Считывание изображения.

Считаем изображение из файла tape.png.

RGB=imread('tape.png');

imshow(RGB);

text(15,15,'Измерение радиуса',...

'FontWeight','bold','Color','y');

Шаг 2: Пороговая обработка изображения.

Преобразуем изображение в бинарное для последующего получения координат границ с помощью функции bwtraceboundary. I=rgb2gray(RGB); threshold=graythresh(I); BW=im2bw(I, threshold); imshow(BW)

Шаг 3: Определение положения начальной точки границы объекта. Для реализации функции bwtraceboundary нужно указать начальную точку определения границы. Эта точка используется для начала отслеживания границы. dim=size(BW); col=round(dim(2)/2)-90; row=min(find(BW(:, col))); Шаг 4: Отслеживание границ. Результатом работы функции bwtraceboundary является поиск (X, Y) - координат точек границ исследуемого объекта. Для более точного определения радиуса исследуемого объекта нужно использовать максимальное число точек, которые принадлежат границе. connectivity=8; num_points=180; contour=bwtraceboundary(BW, [row, col], 'N', connectivity, num_points); imshow(RGB); hold on; plot(contour(:, 2), contour(:, 1), 'g', 'LineWidth', 2);

Шаг 5: Подгонка окружности на выделенные границы.

Отобразим основное уравнение окружности:

(x-xc)^2+(y-yc)^2=radius^2, где (xc, yc) - центр окружности.

Это же уравнение можно представить в другом виде через параметры a, b, c

x^2+y^2+a*x+b*y+c=0, где a=-2*xc, b=-2*yc, c=xc^2+yc^2-radius^2

Параметры a, b, c нужны для вычисления радиуса.

x=contour(:, 2);

y=contour(:, 1);

% при вычислении радиуса будет

% использоваться метод наименьших квадратов

abc=[x y ones(length(x), 1)]\[-(x.^2+y.^2)];

a=abc(1); b=abc(2); c=abc(3);

% вычислим положение центра и радиус

xc=-a/2;

yc=-b/2;

radius=sqrt((xc^2+yc^2)-c)

% отобразим вычисленный центр

plot(xc, yc, 'yx', 'LineWidth', 2);

% отобразим полную окружность

theta=0:0.01:2*pi;

% используем параметрическое представление

% окружности для получения координат точек

Xfit=radius*cos(theta)+xc;

Yfit=radius*sin(theta)+yc;

plot(Xfit, Yfit);

message=sprintf('Вычисленный радиус равен %2.3f пикселей', radius);

text(15, 15, message, 'Color', 'y', 'FontWeight', 'bold');

radius=

80.7567

В работе приведен пример решения задачи определения радиуса окружности мотка ленты. Следует отметить, что представленная задача находит свое применение в реальных системах управления различными автоматизированными технологическими процессами. Такое внедрение компьютерных технологий приводит к повышению эффективности производства. Приведенные выше подходы также нашли свое применение в различных системах дистанционного контроля. Однако точность измерений во многом зависит от качества используемой аппаратуры.

Обнаружение объектов с помощью сегментации изображений

Если контраст объект-фон является достаточным, тогда обнаружить такой объект не представляет труда. Применим средства выделения границ и морфологические функции для обнаружения некоторого медицинского объекта.

Ключевые слова:	Обнаружение границ, структурный элемент, эрозия, морфологические операторы, сегментация
Ключевые функции:	edge, strel, imdilate, imerode, imfill, imclearborder

Демонстрационный пример Демонстрационный пример включает следующие шаги:

• Шаг 1: Считывание изображения.

• Шаг 2: Детектирование всего объекта.

• Шаг 3: Заполнение промежутков.

• Шаг 4: Морфологическое расширение изображений.

• Шаг 5: Заполнение внутренних промежутков.

• Шаг 6: Операции со связными объектами на границе.

• Шаг 7: Сглаживание объекта.

Шаг 1: Считывание изображения. Считаем некоторое медицинское изображение из файла 'cell.tif'. I=imread('cell.tif'); figure, imshow(I), title('original image');

Шаг 2: Детектирование всего объекта. На изображении представлено два объекта и только один из них видно полностью. Наша задача заключается в выделении этого объекта. Другими словами этот объект необходимо сегментировать. Качество сегментации в большой мере зависит от того насколько сильно объект отличается от фона. Поэтому усиление контрастности изображения приводит к увеличению качества детектирования. Одним их путей детектирования является применение оператора Собеля с использованием некоторого порогового значения. В данном случае пороговое значение определим с помощью функции graythresh. Создадим некоторую градиентную маску и применим ее в функции выделения края. BWs=edge(I, 'sobel', (graythresh(I)*.1)); figure, imshow(BWs), title('binary gradient mask');

Шаг 3: Заполнение промежутков. Бинарная градиентная маска отображает линии с высоким контрастом на изображении. Сравнивая исходное и обработанное изображения, можно сделать вывод, что эти линии вполне точно описывают контуры объекта интереса. se90=strel('line', 3, 90); se0=strel('line', 3, 0); Шаг 4: Морфологическое расширение изображений. Бинарная градиентная маска при выполнении морфологической операции расширения использует последовательно вертикальные и горизонтальные структурные элементы. Для выполнения операции морфологического расширения используется функция imdilate. BWsdil=imdilate(BWs, [se90 se0]); figure, imshow(BWsdil), title('dilated gradient mask');

Шаг 5: Заполнение внутренних промежутков. Применение морфологической операции расширения для обработки контуров с помощью градиентной маски дает вполне приемлемый результат, который заключается в заполнении отверстий в середине исследуемого объекта. Для реализации данного шага используется функция imfill. BWdfill=imfill(BWsdil, 'holes'); figure, imshow(BWdfill); title('binary image with filled holes');

Шаг 6: Операции со связными объектами на границе. На предыдущем этапе была проведена операция сегментации. Однако эта операция затронула не только исследуемый объект. Каждый объект, который имеет связные границы на изображении, может быть удален с помощью функции imclearborder. Связность в функции imclearborder устанавливается на уровне 4 по диагонали. BWnobord=imclearborder(BWdfill, 4); figure, imshow(BWnobord), title('cleared border image');

Шаг 7: Сглаживание объекта. И наконец, чтобы проведенная обработка выглядела более натурально, необходимо провести сглаживание с помощью морфологической операции эрозии с использованием специального структурного элемента. Создадим этот элемент с использованием функции strel. seD=strel('diamond', 1); BWfinal=imerode(BWnobord, seD); BWfinal=imerode(BWfinal, seD); figure, imshow(BWfinal), title('segmented image');

Существуют также альтернативные методы отображения сегментированных объектов. Приведем пример с использованием функции bwperim. BWoutline=bwperim(BWfinal); Segout=I; Segout(BWoutline)=255; figure, imshow(Segout), title('outlined original image');

Реконструкция изображений по их проекционным данным

В этом примере рассмотрим применение функций radon, iradon, fanbeam и ifanbeam для получения проекций изображений и обратную реконструкцию данных на основе их проекций. В функциях radon и iradon для получения проекций используется параллельно-лучевая геометрия, а в функциях fanbeam и ifanbeam - веерно-лучевая геометрия. Для того, чтобы сравнить параллельно-лучевую и веерно-лучевую геометрии, создадим некоторые синтезированные проекции для каждой геометрии и затем применим их для восстановления исходного изображения.

Практической задачей, где применяется реконструкция изображений, является рентгеновская томография, где проекции формируются на основании измерения ослабления радиации при прохождении через разные физические среды. В результате на изображении отображаются данные о структуре среза исследуемого объекта. Проекционные данные формируются и хранятся в специальной медицинской аппаратуре. На основании этих данных существует возможность реконструкции изображения исследуемых органов. В среде Matlab моделирование восстановления вида объекта по его томографическим проекциям можно реализовать с помощью функций iradon и ifanbeam.

Функция iradon применяется для реконструкции изображений из параллельно-лучевых проекций. В параллельно-лучевой геометрии каждая проекция формируется в результате линейной комбинации изображений, представляющих некоторый исследуемый орган под разными углами. Функция ifanbeam применяется для реконструкции данных на основании их веерно-лучевых проекций, которые получены из одного источника излучения, но зафиксированы несколькими датчиками.

Рассмотрим применение пакета Image Processing Toolbox для моделирования описанных выше процедур при использовании обеих геометрий.

Содержание

• Создание фантома головы.

• Параллельные лучи - вычисление синтезированных проекций.

• Параллельные лучи - реконструкция фантома головы на основании проекционных данных.

• Веерные лучи - вычисление синтезированных проекций.

• Веерные лучи - реконструкция фантома головы на основании проекционных данных.

Создание фантома головы.

Тестовое изображение фантома головы Шеппа-Логана может быть сгенерировано с помощью функции phantom. Изображение фантома иллюстрирует много особенностей, которые в действительности свойственны реальным томографическим изображениям головы человека. Размещенные в центре фантома объекты в виде эллипса моделируют реальные внутренние особенности строения головы человека.

P=phantom(256);

imshow(P)

Параллельные лучи - вычисление синтезированных проекций.

Рассмотрим вычисление синтезированных проекций на основе параллельно-лучевой геометрии с разным числом углов. Для обработки этих данных используется функция radon, а результат представляется в виде матрицы, в которой каждый столбец представляет собой преобразование Радона для одного значения угла.

theta1=0:10:170;

[R1, xp]=radon(P, theta1);

num_angles_R1=size(R1, 2)

num_angles_R1=

18

theta2=0:5:175;

[R2, xp]=radon(P, theta2);

num_angles_R2=size(R2, 2)

num_angles_R2=

36

theta3=0:2:178;

[R3, xp]=radon(P, theta3);

num_angles_R3=size(R3, 2)

num_angles_R3=

90

Отметим, что для каждого угла проекция вычисляется на основании N точек вдоль оси xp, где N является константой,

которая зависит от расстояния до всех точек изображения, спроектированных под разными углами.

N_R1=size(R1, 1)

N_R2=size(R2, 1)

N_R3=size(R3, 1)

N_R1=

367

N_R2=

367

N_R3=

367

Таким образом, когда мы используем меньший фантом головы, при вычислении проекций требуется меньшее число точек вдоль оси xp.

P_128=phantom(128);

[R_128, xp_128]=radon(P_128, theta1);

N_128=size(R_128, 1)

N_128=

185

Отобразим проекционные данные R3. Некоторые свойства изображения исходного фантома видны на изображении R3.

Первый столбец R3 соответствует проекциям нулевого угла, которые интегрированы в вертикальном направлении.

Последний столбец соответствует проекциям под 90 градусов, которые просуммированы в горизонтальном направлении.

figure, imagesc(theta3, xp, R3)

colormap(hot)

colorbar

xlabel('Параллельный поворот угла - \theta');

ylabel(' Параллельное расположение датчиков x\ ');

На изображении представлена зависимость местонахождения сенсоров от поворота угла.