Математическая морфология.

Электронный математический и медико-биологический журнал. - Т. 14. -

Вып. 2. - 2015. - URL:

http://www.smolensk.ru/user/sgma/MMORPH/TITL.HTM

http://www.smolensk.ru/user/sgma/MMORPH/N-46-html/TITL-46.htm

http://www.smolensk.ru/user/sgma/MMORPH/N-46-html/cont.htm

 

 

МИНИСТЕРСТВО  ОБРАЗОВАНИЯ и НАУКИ  РОССИЙСКОЙ  ФЕДЕРАЦИИ

__________________________________________________________

 

Филиал федерального государственного бюджетного образовательного учреждения высшего профессионального образования 

«Национальный исследовательский университет «МЭИ» в г. Смоленске

 

 

УДК: 56.07

 

 

Кафедра ЭЛЕКТРОНИКИ И МИКРОПРОЦЕССОРНОЙ ТЕХНИКИ

 

Направление подготовки: 11.04.04 - Электроника и наноэлектроника

 

 

 

 

МАГИСТЕРСКАЯ   ДИССЕРТАЦИЯ

 

(shalatkina.docx)

 

Магистерская программа:

«Промышленная электроника и микропроцессорная техника»

 

Тема:

РАЗРАБОТКА АВТОМАТИЗОРОАННОЙ СИСТЕМЫ

 

ОСТЕОМЕТРИЧЕСКИХ ИЗМЕРЕНИЙ

 

Студент

                 ПЭ-13М                                         Шалаткина А.Н.

 

                              группа                                                       подпись                                      фамилия  И.О.             

 

 

Научный руководитель

к.т.н., доцент                                                         Троицкий Ю.В.

 

уч. степень, должность                             подпись                                      фамилия И.О.           

 

 

 

Консультант

 

 

уч. степень, должность                             подпись                                      фамилия И.О.           

 

 

 

 

 

Зав. кафедрой

 

 

д.т.н., доцент                                                        Якименко И.В.

уч. степень, должность                            подпись                                      фамилия И.О.

 

 

Место выполнения магистерской диссертации

     филиал «МЭИ» в гмоленске, кафедра ЭиМТ

 

 

Смоленск 2015 г.

 


АННОТАЦИЯ

Автор работы – Шалаткина Анастасия Николаевна.

Тема — разработка автоматизороанной системы остеометрических измерений.

В соответствии с техническим заданием был произведен обзор литературы по исследуемой теме. После анализа существующих методов была разработана модель автоматизированной системы для остеометрических измерений. На основе этой модели в среде Matlab создан алгоритм для определения длины, ширины и угла скрученности. Проведены измерения реального остеометрического объекта с последующей оценкой точности автоматизированных измерений.  Расчетно-пояснительная записка оформлена с помощью текстового редактора Microsoft Word 2010 и графического редактора Microsoft Visio 2010. Моделирование произведено в системе компьютерной математики Matlab.

Расчетно-пояснительная записка содержит: 80 листов, 60 рисунков, 5 таблиц, 1 приложение, 17 источников литературы.

SUMMARY

The report was written by Shalatkina Anastasiya Nikolaevna.

The theme is the research of method automation osteometric dimensions.

According to the theme and the task the review of the literature was written. After analyzing of existing methods model of automated system for osteometric dimensions was developed. The algorithm for determination the length, width and angle of twist  was created based on Matlab’s model. The real osteometric object was measured with evaluation of the accuracy of automated measurements. The report was written using text editor Microsoft Word 2010 and graphic editor Microsoft Visio 2010. Modeling was carried out using computer-aided design system MATLAB.

The report contains: 80 pages, 60 figures, 5 tables, 1 appendices, 17 items of the literature.

ОГЛАВЛЕНИЕ

1 АНАЛИЗ СУЩЕСТВУЮЩИХ СПОСОБОВ ОБРАБОТКИ ВИДЕОИЗОБРАЖЕНИЙ ПРИ АВТОМАТИЗАЦИИ ОСТЕОМЕТРИЧЕСКИХ ИЗМЕРЕНИЙ.. 7

1.1 Существующие способы осуществления остеометрических измерений  7

1.2 Существующие способы обработки видеоизображений. 15

2 МЕТОД И АЛГОРИТМ ИЗМЕРЕНИЯ ОСТЕОМЕТРИЧЕСКИХ ОБЪЕКТОВ   28

2.1 Модель оптико-физической измерительной системы.. 28

3 МЕТОД И АЛГОРИТМ ОБРАБОТКИ ИЗОБРАЖЕНИЙ СОЗДАННЫЙ НА ОСНОВЕ БИНАРНОГО СПОСОБА.. 35

3.1 Пороговая обработка изображения на основе метода Отсу. 35

3.2 Пороговая обработка цветных изображений. 39

4 АЛГОРИТМ ДЛЯ ОПРЕДЕЛЕНИЯ ПАРАМЕТРОВ В СРЕДЕ MATLAB.. 57

5 ДОКАЗАТЕЛЬСТВО РАБОТОСПОСОБНОСТИ АЛГОРИТМА ДЛЯ ОПРЕДЕЛЕНИЯ ПАРАМЕТРОВ ОБЪЕКТА.. 65

5.1 Результаты измерения параметров остеометрического объекта с использованием ручных средств. 65

5.2 Результаты измерения параметров остеометрического объекта с использованием среды Matlab. 67

5.3 Оценка точности автоматизированных измерений и сравнение с ручными измерениями. 74

ПРИЛОЖЕНИЕ. 84

ПРИЛОЖЕНИЕ А.. 85

ПРИЛОЖЕНИЕ B.. 87

ВВЕДЕНИЕ

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

·        Необходимость сократить время, затрачиваемое человеком на выполнение своих обязанностей и, вследствие этого, увеличить объем выпущенной продукции, измеренных объектов, протестированного оборудования и т.д.;

·        Невысокая точность, увеличить которую можно путем автоматизации процесса;

·        Работа, связанная с хрупкими объектами, способными получить повреждения и потерять свою ценность;

·        Работа, которую человек не способен выполнять на должном уровне по причине несовершенства данных от природы анализаторов (слух, зрение, осязание, мелкая моторика и т.д.);

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

Одной из областей, нуждающихся в автоматизации, является остеометрия. Остеометрия  это совокупность методов измерения костей.

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

Целью является сокращение времени остеометрических измерений и повышение точности путем создания автоматизированной системы остеометрических измерений с использованием оптико-физической измерительной системы, основанной на техническом зрении.

Для достижения этой цели были поставлены следующие основные задачи:

·        провести анализ существующих способов остеометрических измерений и способов обработки изображения;

·        разработать модель автоматизированной системы для остеометрических измерений;

·        на основе разработанной модели создать в среде Matlab алгоритм для определения длины, ширины и угла скрученности кости;

·        провести измерения реального остеометрического объекта и на их основе доказать адекватность работы разработанной системы;

·        оценить точность проведенных автоматизированных измерений и сравнить с ранее существовавшими.

Научная новизна полученных результатов состоит в предложенной  модели установки, с помощью которой определяются остеометрические параметры.

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

По теме диссертационной работы опубликовано 7 (статей и тезисов докладов) на различные научно-технических конференциях.

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

Также рассмотрены способы обработки изображений при автоматизации остеометрических измерений. Широкий спектр алгоритмов цифровой обработки изображений обработки (фильтрации) получены для простейших стационарных моделей объектов и способны улучшать качество, выделять границы, обозначать крупноразмерные объекты и получать бинарные изображения, необходимые для осуществления остеометрических измерений.

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

В третьей главе описаны существующие методы и алгоритмы обработки изображений, которые основаны на бинарном способе. Рассмотрены методы пороговой обработки изображений методом Отсу, Саувола, Кристиана, также рассмотрен метод обработки с использованием адаптивного порога, представлен способ обработки цветных изображений. Предложен оптимальный алгоритм для определения параметров объекта в среде Matlab, приведена блок-схема для определения параметров.

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

В заключении сделаны основные выводы по результатам проведенной работы.

АНАЛИЗ СУЩЕСТВУЮЩИХ СПОСОБОВ ОБРАБОТКИ ВИДЕОИЗОБРАЖЕНИЙ ПРИ АВТОМАТИЗАЦИИ ОСТЕОМЕТРИЧЕСКИХ ИЗМЕРЕНИЙ

1.1 Существующие способы осуществления остеометрических измерений

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

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

Измерительный штатив (Рисунок 1.1). Измерительный штатив состоит из горизонтальной доски, двух вертикальных стенок и подвижной доски. Вертикальные стенки укрепляются под прямым углом друг к другу и к горизонтальной поверхности. Так же как и горизонтальная доска, они должны иметь идеально ровную поверхность. Размеры их могут быть разные. Рекомендуется для длины горизонтальной доски 65 см, для ширины — 22 см. Высота их — 8 см. Но при некоторых измерениях иногда может потребоваться инструмент несколько большего размера.

 

Рисунок 1.1 – Измерительный штатив. Вид сбоку

 

 

Горизонтальная доска прибора расчерчивается таким образом, чтобы и в продольном, и в поперечном направлении можно было определять размеры с точностью до 1 мм, т. е. так же, как расчерчена миллиметровая бумага. Наиболее простой способ достичь этого и заключается в аккуратном наклеивании миллиметровой бумаги на горизонтальную поверхность или в укреплении на ней миллиметровки любым другим способом (нужно следить только, чтобы она плотно прилегала к горизонтальной доске). Нулевая отметка должна соответствовать углу между вертикальными стенками прибора. От этой точки и в продольном и в поперечном направлении миллиметровка должна быть пронумерована по сантиметрам с целью облегчить чтение результатов измерения. Для продольного направления лучше сделать это на той стороне горизонтальной доски, которая свободна от вертикальной стенки, для поперечного — рекомендуется нанести разбивку на сантиметры два или три раза (Рисунок 1.2).

Рисунок 1.2 – Измерительный штатив. Вид сверху

 

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

Измерительный штатив служит для измерения всех длинных костей конечностей — их длины и в ряде случаев некоторых поперечных диаметров. Кость кладется на горизонтальную доску, прижимается либо к продольной, либо к поперечной вертикальной доске и фиксируется в требуемом положении. Затем с помощью подвижной доски фиксируется искомый размер, и его величина прочитывается на миллиметровой измерительной шкале.

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

Измерительный штатив Рида (Рисунок 1.3). Этот штатив в принципе мало отличается от предыдущего и создан на его основе, но он снабжен рядом дополнительных приспособлений, позволяющих фиксировать взаимно пересекающиеся прямые и плоскости, и поэтому пригоден для измерения угловых соотношений. Рекомендуется конструировать его шире, чем предыдущий штатив, так, чтобы его ширина была не менее 25 см. Одна из вертикальных стенок прибора — продольная — не должна доводиться до конца по длине прибора и обычно делается короче на 5 см с тем, чтобы обеспечить подвижность специального приспособления, укрепляемого на этом конце.

 

 

Рисунок 1.3 – Измерительный штатив Рида и периграф. Измеряется изгиб диафиза бедренной кости

 

 

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

Рисунок 1.4 – Измерительный штатив Рида. Вид сверху. Измеряется угол шейки бедренной кости

 

Поперечная стенка также несет на себе аналогичную перекладину, но уже неподвижно укрепленную на стойках. Обе перекладины находятся на одинаковой высоте над горизонтальной доской. На неподвижной перекладине, находящейся над поперечной стенкой прибора, и на эллипсоидном кольце, укреплены подвижные клеммы, между которыми натягиваются нитки. Одна нитка натягивается между двумя клеммами на кольце — другая между клеммами на неподвижной перекладине и определенной точкой на подвижной перекладине, на которой сама нить закрепляется неподвижно. Перекладина при этом поднята над горизонтальной доской. Нитки должны быть туго натянуты, для чего используются маленькие гири. Если измеряется, скажем, угол шейки и головки бедренной кости с телом бедренной кости, то одна длинная нитка, идущая от подвижной перекладины к неподвижной, натягивается вдоль тела кости по его продольной оси, вторая короткая, соединяющая клеммы на кольце,— по оси шейки. Это достигается при помощи подвижных клемм. В последнем случае клеммы на кольце позволяют менять положение нитки на 360°. Угол между ними, соответствующий искомому углу, определяется транспортиром.

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

Прибор состоит из основания, на котором укреплены две стойки (при известной технической смекалке они могут быть заменены одной). Одна из этих стоек градуирована и представляет собой миллиметровую шкалу. Между стоек укреплена подвижная муфта, на которой в свою очередь укреплена изогнутая игла. Игла может быть изогнута в виде полукруга или дважды под прямым углом. Неподвижное положение муфты внизу соответствует нулевому делению шкалы. Подъем муфты означает подъем кончика иглы ровно на столько миллиметров, на сколько поднята муфта. Таким образом, с помощью кончика иглы можно фиксировать положение определенных точек или плоскостей по вертикали, а с помощью измерительной шкалы — определить расстояние между ними.

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

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

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

Рисунок 1.5

Рисунок 1.5 – Параллелограф и держатель. Измеряется угол скрученности большой берцовой кости

 

На эпифизах кости укрепляются две стальные иглы. После этого кость закрепляется в держателе в вертикальном положении. Вертикальное положение кости, что означает вертикальное положение продольной ее оси, может быть определено как визуально, так и при помощи отвеса. Под держатель закладывается на ровную поверхность, на которой он стоит, чистый лист бумаги. При помощи параллелографа на этом листе фиксируется положение концов обеих игл, укрепленных на эпифизах кости. Затем они соединяются прямыми, и угол, образованный этими прямыми, измеряется транспортиром [1; 2].

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

1.2 Существующие способы обработки видеоизображений

При цифровой обработке изображения используется его представление в памяти в виде двумерного цифрового массива (матрицы) , где  - число строк, а  – число столбцов массива. Обработка изображения заключается в выполнении преобразования указанной матрицы , каждый элемент которой может принимать  значений (уровней) квантования. Каждое из  значений будет называться яркостью данного, исходного изображения в заданной точке  (пикселе). Обычно 1<L<24. 24-битовый размер изображения способен обеспечить разрешение слабых изменений в интенсивности излучения. Дальнейшее увеличение количества уровней квантования L не приведет к заметному улучшению качества, а уменьшение битового размера приведет к потере информации и возникновению ложных контуров в изображении. В результате преобразования формируется набор числовых характеристик матрицы или новое, обработанное изображение. Преобразование может касаться значений элементов  или их координат (i,j) (индексов), выполняться над матрицей в целом, группой элементов или над каждым элементом в отдельности [3].

В настоящее время существует группа алгоритмов обработки  изображений: поэлементные преобразования, выделение контуров, обработка скользящим окном (Рисунок 1.6). Названные алгоритмы приводят к улучшению качества при восприятии изображения человеком наблюдателем.

В процессе развития теории и практики работы с изображениями под обработкой изображений принимаются не только преобразования, приводящие к улучшению качества видения изображения, но и те, которые приводят исходное изображение к виду, удобному для дальнейшей их обработки с помощью машин (обеспечение машинного (технического) зрения). В этих случаях интерес вызывают изображения не сами по себе, а та информация, которая в них содержится. Извлечение этой информации и составляет цель обработки видеоинформации (изображений), которое производится автоматически, т. е. без участия оператора. Такая обработка изображений содержит сегментацию изображений, пороговую обработку, проводимую в интересах обнаружения объектов [4].

Рисунок 1.6 – Алгоритмы обработки изображений

 

Вид цифровой обработки изображений, заключающийся в выполнении одного и того же выбранного функционального преобразования для каждого элемента матрицы вне зависимости от его положения и значений других (соседних) элементов, получил название поэлементного преобразования изображений, т. е. перевод значений каждого элемента  f  в новое значение  g  в соответствии с выбранной функциональной зависимостью [5]:

1.1

Размеры входного и выходного изображения здесь совпадают. При практической реализации алгоритмов поэлементного преобразования необходимо вычислять каждое значение преобразованного элемента в соответствии с конкретным видом функции. Алгоритмы поэлементной обработки, как правило, используются для улучшения качества восприятия изображений оператором-наблюдателем (Рисунок 1.7) [6].

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

 

а)

 

б)

Рисунок 1.7– Результаты использования поэлементной обработки изображения: а) изображение до обработки; б) изображение после обработки

 

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

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

 

Обработка изображения выделением контуров

Исследованиями психологов установлено, что в плане зрения распознавания объектов на изображении наиболее информативными являются не значения яркостей объектов, а характеристики их границ – контуров. Другими словами, основная информация заключена не в яркости отдельных областей, а в их очертаниях. Задача выделения контуров состоит в построении изображения именно границ объектов и очертаний однородных областей. Контуром изображения называется совокупность его пикселей, в окрестности которых наблюдается скачкообразное изменение функции яркости.

В соответствии с алгоритмом обработки изображения выделения контуров вычисление производится в два этапа. Сначала изображение обрабатывается  (проводится двумерная свертка) двумя двумерными фильтрами. Импульсные характеристики этих фильтров соответствуют «маскам» («окнам»):

и

1.2

На втором шаге проводится квадратурная обработка:

1.3

Итогом, таковой является изображение с обозначенными контурами (Рисунок 1.8).

Таким образом, операция выделения контуров заключается в выполнении нелинейной свертки изображения и «окна» в двух противоположных направлениях с последующей квадратурной обработкой для получения контурного изображения.

 

 

а)

 

б)

Рисунок 1.8 – Результаты использования обработки изображения выделением контуров: а) изображение до обработки; б) изображение после обработки

Сегментация изображения

Метод сегментации (препарирования) изображений заключается в  приведении их к виду, далекому от естественного, но удобному для визуальной интерпретации или дальнейшего машинного анализа. Сегментация состоит в принятии решения о наличии объекта для каждой точки изображения. Смысл этого уровня состоит в автоматическом разбиении изображения на содержательно интерпретируемые области. Важным в прикладном отношении значением этого уровня обработки информации является подготовка к процессу обнаружения (выделения) объектов (целей), различных как по своим яркостным и геометрическим свойствам, так и по физической интерпретации [6-9].

Математическая постановка задачи сегментации. Пусть f(x,y) – функция яркости анализируемого изображения; X - конечное подмножество точек плоскости, на котором определена функция f(x,y); S = {S1,S2,…,Sk} – разбиение X на К непустых связных подмножеств Si, i = 1,2,...,К; LP – предикат, определенный на множестве S и принимающий истинные значения тогда и только тогда, когда любая пара точек из каждого подмножества Si удовлетворяет некоторому критерию однородности.

Сегментацией   изображения   f(x,y)   по   предикату   LP   называется   разбиение S = {S1,S2,…,Sk}, удовлетворяющее условиям:

а)

б)  

в)

г)

Условие (а) означает, что каждая точка изображения должна быть отнесена к некоторой области, (б) означает, что области  Si должны быть связными, (в) определяет вид (тип) однородности получаемых областей, (г) выражает свойство «максимальности» областей разбиения. Подразумевается, что разбиение S существует и оно единственно. Предикат LP в (б) называется предикатом однородности и его истинное или ложное значение зависит от свойств функции f(x,y).

Предикат может быть определен как

1.4

где  m=1, 2,…, M,   M – число точек в области Si, либо как

1.5

где     – произвольные точки из области Si;

Т – величина некоторого наперед заданного порога.

Таким образом, сегментацию можно рассматривать как оператор вида

 при  

1.6

где  – функции, определяющие исходное и сегментированное изображения соответственно; - метка (имя) i-й области.

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

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

1.7

Таким образом, пороговая обработка является одним из основных подходов к решению задачи сегментации изображений. Это объясняется ее «относительной простотой в вычислительном аспекте, а также тем, что наиболее характерным признаком объектов в реальных сценах является перепад яркости.

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

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

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

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

Выбор глобального порога непосредственно по гистограмме яркостей. Введем следующие обозначения. Сегментации подвергается изображение f(x,у); x=1,…M; y=1,…N с динамическим диапазоном L уровней яркости. Каждому уровню i соответствует ni точек изображения. Общее количество точек равно n=M×N=n1+n2+...+nL. Гистограмма яркостей рассматривается как оценка распределения вероятностей:

, i=1, L; рi≥0,

1.8

Осуществляется сегментация на два класса точек по пороговому уровню t p0={(х,у): f(x,у) t*} — класс точек фона; p1={(x,у) : f(x,y)>t*} —класс точек объекта:

1.9

Самым распространенным методом использования гистограмм яркостей для пороговой обработки изображений является так называемый метод мод [21, 39-41]. Метод вытекает из предположения, что изображение содержит известное число однородных по яркости классов точек, т. е. таких классов, распределения вероятностей яркости которых унимодальны. Кроме того, считается, что граничные участки между замкнутыми областями занимают сравнительно небольшую площадь изображения. Поэтому на гистограмме им должны соответствовать межмодовые впадины, в пределах которых устанавливаются пороги сегментации (Рисунок 1.9).

Исходя из этого изображение f(x,у) преобразуется в бинарное S(x,у). Такое преобразование осуществляется в первую очередь для того, чтобы сократить информационную избыточность изображения, оставить в нем только ту информацию, которая нужна для решения конкретной задачи. В бинарном изображении должны быть сохранены необходимые  детали (например, очертания изображенных объектов) и исключены второстепенные:

1.10

Рисунок 1.9 – Одномерный срез перепада яркости и гистограмма яркостей

 

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

Для автоматического определения порога сегментации может быть использован итерационный алгоритм, в соответствии с которым вначале выбирается начальная оценка порога Т, затем производится пороговая обработка изображения с порогом Т, в результате чего образуется две группы пикселей: G1 и G2. После этого вычисляются средние значения яркостей μ0 и μ1 по областям G1 и G2. На основе μ0 и μ1 определяется новое значение порога  (рисунок 1.10).

Рисунок 1.10 – Полимодальная  гистограмма яркостей изображения объекта и фона

 

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

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

Метод обработки изображений скользящим окном

Практически любое изображение зашумлено, а значит, возникает проблема борьбы с шумами. Одним из наиболее распространенных способов борьбы с шумами является фильтрация, которая с самого начала была применена к обработке одномерных сигналов. Не составляет исключения и борьба с шумами на зашумленных изображениях, но с особенностями, связанными с двумерностью видеоинформации. Указанная особенность приводит к тому, что вместо одномерного фильтра, предназначенного для фильтрации зашумленных сигналов одной переменной (обычно этой переменной является время), для фильтрации изображений чаще всего применяются двумерные фильтры, представляющие собой апертуры (маски) различной конфигурации на плоскости. Каждому элементу апертуры соответствует (присваивается исследователем) число, в дальнейшем называемое весовым множителем [6-9].

Фильтрация осуществляется перемещением слева – направо (или сверху – вниз) апертуры фильтра (маски) на один пиксель. При каждом положении апертуры производятся однотипные операции, а именно перемножение весовых множителей с соответствующими значениями яркостей исходного изображения и суммированием полученных результатов. Полученное значение делится на заранее заданное число (нормирующий множитель), которое и присваивается центральному элементу апертуры (это и есть «выход» фильтра). Размеры апертуры берутся таким образом, чтобы центральный элемент апертуры определялся однозначно. Наиболее распространенные размеры апертуры 3×3 и 5×5.

Указанная процедура имеет вид:

1.11

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

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

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

 

а)

 

б)

 Рисунок 1.11– Результаты использования обработки изображений скользящим окном: а) изображение до обработки; б) изображение после обработки

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

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

Именно поэтому необходима автоматизированная система измерений, способная сократить время работы, повысить точность и уменьшить погрешность полученных результатов. Для реализации необходимой системы был проведен анализ алгоритмов цифровой обработки.  На сегодняшний день широкий спектр алгоритмов цифровой обработки изображений обработки (фильтрации) получены для простейших стационарных моделей объектов и способны улучшать качество, выделять границы, обозначать крупноразмерные объекты и получать бинарные изображения, необходимые для осуществления остеометрических измерений.

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

2 МЕТОД И АЛГОРИТМ ИЗМЕРЕНИЯ ОСТЕОМЕТРИЧЕСКИХ ОБЪЕКТОВ

2.1 Модель оптико-физической измерительной системы

Выше был описан остеометрический инструментарий, однако все измерения данными приборами не точны, достаточно продолжительны и обладают большой стоимостью. Один из вариантов оптимизировать данный процесс – проведение остеометрических измерений с использованием оптико-физической измерительной системы, основанной на обработке фотоизображений. Модель этой системы представлена на Рисунок 2.1.

Рисунок 2.1 – Модель оптико-физической измерительной системы

 

1 – камера для съёмки горизонтальной плоскости

2 – штатив

3 – лапка с подвижной иголкой

4 – подвижный шарнир

5 – иголка

6 – лапка, фиксирующая кость

7 – фронтальный экран

8 – горизонтальный экран

9 – модель кости

 

Примечание: в установке предусмотрена камера, фотографирующая фронтальную плоскость (на рисунке не представлена). Камера 1 свободно перемещается вдоль горизонтальной оси штатива 2.

Принцип работы установки.        

Исследуемая кость помещается в фиксирующую лапку 6. Затем, путём перемещения лапок 3 и вращения подвижных шарниров 4, достигается такое положение иголок, при котором они будут совмещены с эпифизами кости. Проводится съёмка кости в двух плоскостях, двумя камерами. Полученное изображение передаётся на одноплатный компьютер, с целью дальнейшей обработки изображения в программе Matlab, сохранения и систематизации полученных данных [9; 10].

Таким образом, функциональная схема разработанной модели (Рисунок 2.2) выглядит следующим образом:

 

 

Рисунок 2.2 – Функциональная схема оптико-физической измерительной системы

 

Особое влияние на точность полученных измерений оказывает правильно выбранная цифровая камера.

Принцип работы цифрового фотоаппарата (Рисунок 2.3). Изображение преломляется через систему оптики, преобразуется в цифровую информацию на матрице, от разрешающей способности которой и будет зависеть качество снимка. Затем перекодированное изображение в цифровом виде сохраняется на сменном носителе информации. Информацию в виде изображения можно редактировать, перезаписывать и отправлять на другие носители данных [11].

 

 

Рисунок 2.3  – Устройство цифрового фотоаппарата

 

Корпус современного цифрового фотоаппарата значительно тоньше обычного пленочного и имеет место для ЖК экрана, встроенного в корпус, либо выдвижного, и слоты для карт памяти.

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

В профессиональных цифровых фотоаппаратах объектив практически ничем не отличается от аналоговых фотокамер. Он также состоит из линз и набора зеркал и имеет те же механические функции.

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

Микропроцессор (Рисунок 2.4) отвечает за все функции работы цифровой камеры. Все рычаги управления камеры ведут к процессору в котором зашита программная оболочка (прошивка), которая отвечает за действия фотокамеры: работа видоискателя, автофокус, программные сцены съемки, настройки и функции, электрический привод выдвижного объектива, работа фотовспышки [12].

Рисунок 2.4 – Структурная схема цифрового фотоаппарата

 

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

Полученное изображение сохраняется в памяти фотоаппарата в виде информации на внутренней, либо внешней памяти. Фотоаппараты имеют разъемы для карт памяти SD, MMC, CF, XD-Picture и др., а также разъемы для подключения к другим источникам хранения информации компьютеру, HDD сменным носителям и т.п [13].

Применительно к предлагаемой модели в качестве цифровой камеры можно использовать ультразум фотоаппарат Sony Cyber-shot DSC-HX300, основные параметры которой приведены в Таблица 2.1. Использование камеры с подобными параметрами – оптимальный вариант для установки.

Таблица 2.1 – Основные параметры фотокамеры Sony Cyber-shot DSC-HX300

Матрица

Тип матрицы

CMOS

Разрешение матрицы (Мп) 

20.4

Разрешение матрицы (точек)

5184 x 3888

Физический размер матрицы

1/2.3 дюйма

Скорость съемки (кадр/сек) 

10

Максимальная чувствительность ISO

80 - 3200 ISO, Auto ISO, ISO6400, ISO12800

Диапазон выдержки 

30 - 1/4000 с

Объектив

Фокусное расстояние объектива (мм) 

24 - 1200

Диафрагма

F2.8 - F6.3

Оптический зум 

50

Цифровой зум 

2

Стабилизация изображения 

да

Память

Тип карты памяти

SD, SDHC, SDXC, MS Duo, MS Pro Duo

Интерфейсы

USB-порт

да

Видеовыход

micro-HDMI

Корпус

Высота (мм)

93

Ширина (мм)

130

Толщина (мм)

103

Вес (г)

650

 

В качестве одноплатного компьютера выбран Cubietruck Cubieboard3 (Рисунок 2.5). Cubietruck (так же известный под именем Cubieboard3) представляет собой новый миникомпьютер компании Cubiesteam, имеющий процессор с двумя ядрами ARM Cortex-A7, 2 ГБ ОЗУ, модули Gigabit Ethernet, WiFi и Bluetooth [14].

 

 

Рисунок 2.5 – Cubietruck Cubieboard3

 

При размерах платы 110 × 80 компьютер имеет множество портов ввода/вывода, а также место для установки кожуха, где может разместиться жесткий диск 2.5" или твердотельный накопитель.

Система оснащена слотом для карт microSD, двумя портами USB, портами HDMI и VGA, разъемом для наушников, интерфейсом SATA 2.0 и дополнительными выводами, которые можно использовать для программирования, отладки или соединения с периферийными устройствами.

Плата требовательна к мощности источника питания. Рекомендуемый блок питания 5В 2А, работать от USB порт компьютера будет всего несколько секнуд, потом выключится. Если не использовать жесткий диск, то устойчиво работает от блока питания 1А.

Характеристики:
- Процессор Allwinner A20 с двумя ядрами ARM Cortex-A7, графический ускоритель ARM Mali400 MP2 с поддержкой OpenGL ES 2.0/1.1

- ОЗУ 2 ГБ DDR3 480 МГц

- Встроенный разъем для подключения к мониторам HDMI и VGA 1080p

- 10 Мб / 100 Мб / 1 Гб Ethernet

- Модули WiFi и Bluetooth со встроенной на плате антенной

- Интерфейс SATA 2.0 с поддержкой жестких дисков 2.5" (для винчестеров 3.5" требуется еще одна линия питания 12 В)

- Память хранения данных: NAND + MicroSD или TSD + MicroSD или 2 × MicroSD

- 2 × USB HOST, OTG, Toslink (SPDIF Optical), ИК-порт, 4 светодиода, 1 разъем для наушников

- Питание: 5 В 2.5 А с подключенным жестким диском, поддержка литиевого аккумулятора и часов реального времени

- 54 дополнительных вывода, включая I2S, I2C, SPI, CVBS, два АЦП низкого разрешения, UART, PS2, два канала ШИМ, TS/CSI, IRDA, линейный вход, вход микрофона, 4 ТВ-входа

- Размеры платы: 110 × 80 × 1.4 мм.

Одним из самых главных достоинств является поддержка двух систем Google Android и Ubuntu Linux. Логика Matlab стабильно работает на основе системы Linux [1].

Посредством встроенной на компьютере антенны осуществляется передача данных на компьютер, планшет с соответствующим программным обеспечением.

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

3 МЕТОД И АЛГОРИТМ ОБРАБОТКИ ИЗОБРАЖЕНИЙ СОЗДАННЫЙ НА ОСНОВЕ БИНАРНОГО СПОСОБА

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

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

Синтаксис:

level=graythresh(I)

Описание:

Функция level=graythresh(I) вычисляет значение глобального порога, который используется для преобразования яркостей (интенсивностей) изображения в бинарное изображение в функции im2bw.

Значения нормализованных интенсивностей находятся в диапазоне [0, 1].

Функция graythresh использует метод Отсу, который выбирает порог путем минимизации различных вариантов черных и белых пикселов.

Многомерные массив автоматически преобразуется в двумерный с использованием функции reshape. Функция graythresh игнорирует все ненулевые мнимые части массива I.

Требования к исходным данным:

Исходное изображение I может быть представлено неразреженным массивом с форматом представления данных uint8, uint16 или double. Возвращаемое значение имеет формат представления данных double

Для того, чтобы решить эту проблему, существует несколько не очень хитрых способов. Что бы решить проблему бинаризации используют адаптивный порог: т.е. такой порог, который зависит от локальных (в окне, например, размером 19x19 пикселей) особенностей изображения. Для примера возьмём три изображения (Рисунок 3.1):

Рисунок 3.1 – Исходные изображения

 

Одни из таких методов был предложен Ниблэком (англ. Niblack), основная идея которого заключается в создании пороговой поверхности, базирующейся на математическом ожидании и среднеквадратическом отклонении в некотором окне:

T=μ+kσ

3.1

где k - константа (обычно -0.2).

Математическое ожидание в окне можно найти следующим образом: 

μ=(1/N)*∑xi,

3.2

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

σ=((1/N)*∑⎷(Xi-M)2)1/2

3.3

 Саувола (англ. Sauvola) улучшил метод Ниблэка, введя гипотезу о значениях цвета текстовых пикселей (т.е. на изображении имеется текст): текст близок к 0, а все нетекстовые пиксели близки к 255.

T=μ(1−k(1−σR)),

3.4

где k берётся равным 0.5, а R описывает стандартное отклонение динамического диапазона серого изображения и берётся равным 128. Данный метод будет хорошо справляться, где более или менее соблюдаются условия гипотезы. При «тёмном фоне», он будет выдавать плохие результаты [15].

 

Рисунок 3.2 ‑ Метод Саувола

 

 Решая эту проблема Кристиан (англ. Christian) модифицировал формулу до вида: 

T=(1−k)μ+kM+kσR(μ−M)

3.5

где M — минимальное значение цвета в сером изображении, а R — устанавливается как максимальное значение среднеквадратического отклонения, посчитанные для всего изображения (Рисунок 3.3).

Рисунок 3.3– Метод Кристиана

 

 В предыдущем способе есть недостаток: параметры R и M считаются для всего изображения, из-за чего появились артефакты. Связанно это с тем, что используются глобальные характеристики изображения. Дабы улучшить бинаризацию можно использовать два окна: одно больше другого (Рисунок 3.4).

 

 

 

Рисунок 3.4 ‑ Двухоконное изображение

 

 В малом окне считается математическое ожидание и отклонение, а в большом параметры Rσ и M. Такой подход позволяет улучшить качество бинаризации:

 

 

Рисунок 2.12 - Пример бинаризации с использованием метода двухоконного изображения

 

Поиск локального порога можно произвести по следующей формуле: 

T=(1−α1)μ+α2(σRσ(μ−M)+αM)

3.6

где α1=k1(σRσ)γ, α2=k2(σRσ)γ, а α1, γ, k1 и k2 — положительные константы.

При γ = 2 параметры α1, k1 и k2рекомендуется брать в диапазоне 0.1-0.2, 0.15-0.25 и 0.01-0.05.

Мы уже рассматривали один из способов бинаризации, в случае, когда имеется неоднородность освещения. Суть метода достаточно проста:

1.     Используя Гауссов фильтр с большим размером окна сильно размываем изображение.

2.     Вычитаем из изображения фон (абсолютное значение).

3.     Инвертируем результат вычитания.

4.     Бинаризуем с использованием метода Отсу ().

а) б)

 

 

Рисунок 3.5 ‑ Пример бинаризации: а) исходное изображение; б) бинаризированное изображение

3.2 Пороговая обработка цветных изображений

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

Для этого сначала считаем и визуализируем некоторое цветное изображение (Рисунок 3.6).

L=imread('fabric.png');

L=double(L)./255;

figure, imshow(L);

title('Исходное изображение');

 

Рисунок 3.6 ‑ Исходное изображение

 

Теперь сформируем и визуализируем значения интенсивности исходного изображения и построим их гистограмму (Рисунок 3.7).

I=L(:,:,1)+L(:,:,2)+L(:,:,3);

figure, subplot(211);imshow(I);title('Интенсивность');

subplot(212); imhist(I);

imwrite(I,'I.bmp');

 

Рисунок 3.7 ‑ Интенсивность и гистограмма интенсивности

 

Сформируем нормированное значение r составляющей изображения и визуализируем ее. Также построим гистограмму полученных данных (Рисунок 3.8).

r=L(:,:,1)./(I+eps);

figure, subplot(211);imshow(r);title('Составляющая r');

subplot(212); imhist(r);

imwrite(r,'r.bmp');

Рисунок 3.8 ‑ r-составляющая и ее гистограмма

 

Аналогично сформируем и визуализируем g цветовую составляющую исходного изображения и построим ее гистограмму (Рисунок 3.9).

g=L(:,:,2)./(I+eps);figure,

subplot(211);

imshow(g);

title('Составляющая g');

subplot(212);

imhist(g);

imwrite(g,'g.bmp');

 

Рисунок 3.9 – g-составляющая и ее гистограмма

 

Далее сформируем и визуализируем b цветовую составляющую исходного изображения и построим ее гистограмму (Рисунок 3.10).

b=L(:,:,3)./(I+eps);

figure, subplot(211);imshow(b);

title('Составляющая b');

subplot(212); imhist(b);

imwrite(b,'g.bmp');

 

 

Рисунок 3.10 ‑ Составляющая b и ее гистограмма

 

Также необходимо сформировать значения насыщенности S изображения.

[N M]=size(r);

for i=1:N;

    disp(i);

    for j=1:M;

        MINIMUM=min([r(i,j) g(i,j) b(i,j)]);

S(i,j)=1-3*MINIMUM;

    end;

end;

S(S<0)=0;

Визуализируем эти данные и построим их гистограмму (Рисунок 1.1Рисунок 3.11).

figure, subplot(211);imshow(S);title('Насыщенность S');

subplot(212); imhist(S);

imwrite(S,'S.bmp');

 

Рисунок 3.11 ‑ Гистограмма насыщенности

 

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

F=acos((2*r-g-b)./(sqrt(6).*sqrt( (r-1/3).^2+(g-1/3).^2+(b-1/3).^2)));

figure, subplot(211);imshow(F);title('Цветовой фон F');

subplot(212); imhist(F);

imwrite(F,'F.bmp');

 

Рисунок 3.12 ‑Гистограмма цветового фона F

 

Таким образом, нами сформированы наборы данных (r-, g-, b-цветовые составляющие, насыщенность S, интенсивность I и фон F), которые будут использованы при проведении пороговых операций над исходным изображением.

Далее на основе анализа гистограммы полученных данных выберем пороговое значение и проведем бинаризацию изображения:

1)    Для изображения значений интенсивности исходного изображения (Рисунок 3.13). Пороговое значение равно среднеарифметическому значению уровней интенсивности.

а)

б)

в)

г)

Рисунок 3.13– Поэтапная обработка: а) исходное изображение; б) изображение интенсивности; в) гистограмма изображения интенсивности; г) результирующее изображение

 

2)    Для r цветовой составляющей изображения (Рисунок 3.14). Гистограмма r цветовой составляющей изображения имеет несколько ярко выраженных мод, которые и будут определять пороговые значения для обработки этого изображения. Это позволит выделить несколько типов объектов.

а) Пороговое значение определяется диапазоном [0, 0.31).

а)

б)

в)

г)

Рисунок 3.14 ‑ Поэтапная обработка: а) исходное изображение; б) изображение цветной составляющей r; в) гистограмма изображения цветной составляющей r; г) результирующее изображение

б) Пороговое значение определяется диапазоном [0.31, 0.39) (Рисунок 3.15, Рисунок 3.16).

а)

б)

Рисунок 3.15 – Поэтапная обработка: а) исходное изображение; б) изображение цветной составляющей r

в)

г)

Рисунок 3.16 ‑ Поэтапная обработка: в) гистограмма изображения цветной составляющей r; г) результирующее изображение

 

3)    Для g цветовой составляющей изображения. Также как и в предыдущем случае, гистограмма распределения этой цветовой составляющей имеет несколько ярко выраженных мод, которые будут основой для выбора порогового значения.

а) Пороговое значение определяется диапазоном [0, 0.2). (Рисунок 3.17, )

а)

б)

Рисунок 3.17 – Поэтапная обработка: а) исходное изображение; б) изображение цветной составляющей g;

 

в)

г)

Рисунок 3.18 – Поэтапная обработка: в) гистограмма изображения цветной составляющей g; г) результирующее изображение

б) Пороговое значение определяется диапазоном [0.2, 0.25). (Рисунок 3.19)

а)

б)

в)

г)

Рисунок 3.19 – Поэтапная обработка: а) исходное изображение; б) изображение цветной составляющей g; в) гистограмма изображения цветной составляющей g; г) результирующее изображение

 

4)    Для b цветовой составляющей изображения. Гистограмма этой цветовой составляющей также имеет несколько мод, анализ которых будет служить основой для выбора пороговых значений.

а) Пороговое значение определяется диапазоном [0, 0.35).(Рисунок 3.20; Рисунок 3.21).

а) б)

Рисунок 3.20 – Поэтапная обработка: а) исходное изображение; б)изображение цветной составляющей b

в) г)

Рисунок 3.21 ‑Поэтапная обработка: в) гистограмма изображения цветной составляющей b; г) результирующее изображение

б) Пороговое значение определяется диапазоном [0.35, 0,4). (Рисунок 3.22; Рисунок 3.23)

а)

б)

Рисунок 3.22 ‑ Поэтапная обработка: а) исходное изображение; б) изображение цветной составляющей b

в)

г)

Рисунок 3.23 ‑ Поэтапная обработка: в) гистограмма изображения цветной составляющей b; г) результирующее изображение

 

в) Пороговое значение определяется диапазоном [0.4, 0,5).(Рисунок 3.24, Рисунок 3.24)

а)

б)

Рисунок 3.25– Поэтапная обработка: а) исходное изображение; б) изображение цветной составляющей b

 

в)

г)

Рисунок 3.26– Поэтапная обработка в) гистограмма изображения цветной составляющей b; г) результирующее изображение

 

5) Для значений насыщенности S изображения. В этом случае пороговое значение можно выбирать как по гистограмме, так и экспериментальным путем. Пусть пороговое значение будет определяться диапазоном [0, 0.25), (Рисунок 3.27; Рисунок 3.28).

а)

б)

Рисунок 3.27 – Поэтапная обработка: а) исходное изображение; б) изображение насыщенности s

в)

г)

Рисунок 3.28 ‑ Поэтапная обработка: в) гистограмма изображения насыщенности s; г) результирующее изображение

 

6) Для значений цветового фона F изображения. В этом случае гистограмма имеет две моды, на основе которых будут определены пороговые значения.

а) Пороговое значение определяется диапазоном [0, 0.25), (Рисунок 3.29).

а)

б)

в)

г)

Рисунок 3.29 ‑ Поэтапная обработка: а) исходное изображение; б) изображение цветового фона F; в) гистограмма изображения цветового фона F; г) результирующее изображение

 

б) Пороговое значение определяется диапазоном [0.6, 1), (Рисунок 3.30).

а)

б)

в)

г)

Рисунок 3.30– Поэтапная обработка: а) исходное изображение; б) изображение цветового фона F; в) гистограмма изображения цветового фона F; г) результирующее изображение

 

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

clear;

L=imread('fabric.png');

L=double(L)./255;

figure,imshow(L);title('Исходное изображение');

R=imread('r.bmp');

R=R(:,:,1);

R=double(R)./255;

G=imread('g.bmp');

G=G(:,:,1);

G=double(G)./255;

B=imread('b.bmp');

B=B(:,:,1);

B=double(B)./255;

F=imread('F.bmp');

F=F(:,:,1);

F=double(F)./255;

I=imread('I.bmp');

I=I(:,:,1);

I=double(I)./255;

S=imread('S.bmp');

S=S(:,:,1);

S=double(S)./255;

[N M]=size(R);

%Lser=mean(mean(R));

for i=1:N;

    disp(i);

    for j=1:M;

        if  ((R(i,j)-0)<.4)&((G(i,j)-0)>.35)&((B(i,j)-0)<.4)&

        ((F(i,j)-0)>.1)&((S(i,j)-0)>.1)&((I(i,j)-0)>.1);

            L(i,j,1)=0;L(i,j,2)=1;L(i,j,3)=0;           

        end;

    end;

end;

figure,imshow(L);title('Результирующее изображение');

Результирующее изображение представлено на Рисунок 3.31.

Рисунок 3.31 ‑ Результирующее изображение

 

Тексты программ, которые использовались при подготовке материала.

1) Программа, которая реализует формирование исходных данных для проведения пороговых операций.

clear;

%L=imread('peppers.png');

%L=imread('tissue.png');

%L=imread('greens.jpg');

L=imread('fabric.png');

L=double(L)./255;

figure, imshow(L);title('Исходное изображение');

 

I=L(:,:,1)+L(:,:,2)+L(:,:,3);

figure, subplot(211);imshow(I);title('Интенсивность');

subplot(212); imhist(I);

imwrite(I,'I.bmp');

 

r=L(:,:,1)./(I+eps);

figure, subplot(211);imshow(r);title('Составляющая r');

subplot(212); imhist(r);

imwrite(r,'r.bmp');

 

g=L(:,:,2)./(I+eps);

figure, subplot(211);imshow(g);title('Составляющая g');

subplot(212); imhist(g);

imwrite(g,'g.bmp');

 

b=L(:,:,3)./(I+eps);

figure, subplot(211);imshow(b);title('Составляющая b');

subplot(212); imhist(b);

imwrite(b,'b.bmp');

 

[N M]=size(r);

for i=1:N;

    disp(i);

    for j=1:M;

        MINIMUM=min([r(i,j) g(i,j) b(i,j)]);

S(i,j)=1-3*MINIMUM;

    end;

end;

S(S<0)=0;

figure, subplot(211);imshow(S);title('Насыщенность S');

subplot(212); imhist(S);

imwrite(S,'S.bmp');

 

F=acos((2*r-g-b)./(sqrt(6).*sqrt( (r-1/3).^2+(g-1/3).^2+(b-1/3).^2)));

figure, subplot(211);imshow(F);title('Цветовой фон F');

subplot(212); imhist(F);

imwrite(F,'F.bmp');

2) Программа, которая реализует пороговую обработку на основе значений цветовых составляющих, интенсивности, насыщенности и значений фона.

clear;

L=imread('fabric.png');

figure, subplot(221);imshow(L);title('Исходное изображение');

L=imread('F.bmp');

L=L(:,:,1);

L=double(L)./255;

subplot(222);imshow(L);title('Изображение цветового фона F');

subplot(223);imhist(L);grid;title('Гистограмма изобр-ия цветового фона F');

 [N M]=size(L);

Lser=mean(mean(L));

for i=1:N;

    disp(i);

    for j=1:M;

        if  L(i,j)>0.6;

            L(i,j)=1;

        else L(i,j)=0;

        end;

    end;

end;

subplot(224);imshow(L);title('Результирующее изображение');

Выводы

1. Широкий спектр алгоритмов цифровой обработки изображений обработки (фильтрации) получены для простейших стационарных моделей объектов и способны улучшать качество, выделять границы, обозначать крупноразмерные объекты и получать бинарные изображения, необходимые для обеспечения реализации остеометрических измерений.

2. Ни один из методов не даст нам полностью автоматизированный процесс бинаризации, но при совместном использовании обоих методов задача будет решена.

4 АЛГОРИТМ ДЛЯ ОПРЕДЕЛЕНИЯ ПАРАМЕТРОВ В СРЕДЕ MATLAB

Каждый объект изображения можно охарактеризовать набором признаков, которые могут служить основой для их анализа и распознавания. В Matlab для вычисления признаков объектов используется функция regionprops [9].

Все типы признаков объектов, которые вычисляются с помощью функции regionprops, условно можно разбить на пять групп:

1.     признаки типов изображений объекта;

2.     признаки размеров объекта;

3.     признаки определения площади объекта;

4.     признаки формы объекта;

5.     другие признаки.

Остановимся подробнее на первых двух признаках.

Сначала считаем некоторое исходное изображение (Рисунок 4.1), которое будем использовать и визуализируем его.

L=imread('test_image.bmp');

figure, imshow(L);

Рисунок 4.1 – Исходное изображение

 

Далее проиллюстрируем разницу между различными типами изображений. Чтобы приступить к выполнению этого задания сначала необходимо выполнить поиск объектов на изображении и вычисление их признаков:

[L num]=bwlabel(L,8);

Сформируем изображение самого объекта (Рисунок 4.2).

STATS=regionprops (L,'Image');

Image=STATS(1).Image;     

figure, imshow(Image);

 

Рисунок 4.2 ‑ Изображение самого объекта

 

Сформируем изображение с “залитыми” дырами (Рисунок 4.3)

STATS=regionprops(L,'FilledImage');

FilledImage=STATS(1).FilledImage;

figure, imshow(FilledImage);

 

Рисунок 4.3– Изображение с “залитыми” дырами

 

Сформируем изображение “залитого” пикселями объекта (Рисунок 4.4) выпуклого многоугольника, в который вписан объект:

STATS=regionprops(L,'ConvexImage');

ConvexImage=STATS(1).ConvexImage;

figure, imshow(ConvexImage);

Рисунок 4.4 – Изображение “залитого” пикселями объекта выпуклого многоугольника, в который вписан объект

 

Довольно часто возникает необходимость описания габаритных размеров объекта различными способами. Для этого используются следующие опции – “BoundingBox”, ”ConvexHull”, ”Extrema”,”PixelList”. Рассмотрим использование этих опций на примере исходного изображения[4,5].

1.     Получение координат прямоугольника, ограничивающего объект.

STATS=regionprops(L,'BoundingBox');

S

S =

                    BoundingBox: [31.5000 37.5000 205 181]

Результат представлен в виде массива [x y width height] , где (x,y) – координаты левого верхнего угла прямоугольника, width – ширина, height – высота прямоугольника.

Используя код

line([BoundingBox (1) BoundingBox (1)+ BoundingBox (3) BoundingBox (1)+ BoundingBox (3)

     BoundingBox (1) BoundingBox (1)],[ BoundingBox (2)  BoundingBox (2)

     BoundingBox (2)+ BoundingBox (4)  BoundingBox (2)+ BoundingBox (4) BoundingBox (2)]);

Прямоугольник, который ограничивает объект, для наглядности можно визуализировать (Рисунок 4.5).

Рисунок 4.5 – Прямоугольник, ограничивающий объект

 

2.     Получение координат выпуклого многоугольника, в который вписан объект. Для этого используется опция “ConvexHull”.

STATS=regionprops(L,'ConvexHull',8);

ConvexHull=STATS.ConvexHull;

Результат представляется в виде матрицы, каждая строка которой содержит (x,y) координаты вершин описанного выше многоугольника.

Для наглядности визуализируем этот многоугольник (Рисунок 4.6) на основании полученных координат[7].

line(ConvexHull(:,1),ConvexHull(:,2))

Рисунок 4.6 – Выпуклый многоугольник, в который вписан объект

3.     Опция “Extrema” позволяет получить экстремальные координаты объекта. Результат представлен матрицей, которая содержит все экстремальные координаты объекта (Рисунок 4.7).

STATS=regionprops(L,'Extrema');

Для наглядности отметим эти точки на изображении.

Рисунок 4.7 – Экстремальные точки изображения

 

Будем пользоваться первым способом измерения размеров объекта. Суть измерения длины кости сводится к нахождению размеров (в пикселях) прямоугольника, ограничивающего исследуемый объект, и умножению этого числа на заранее известный размер пикселя, который составляет 0,02636 см.

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

Рисунок 4.8 – Изображение кости со спицами.

Чтобы вычислить угол необходимо знать координаты экстремумов, далее при помощи тригонометрии вычисляем угол и выводим его значение. Используя код:

level=graythresh(L);

 bw=im2bw(L,level);

 STATS=regionprops(bw,'FilledImage');

Image=STATS(1).FilledImage;

 figure, imshow(Image);

 imwrite(Image,'Image.bmp');

STATS=regionprops(Image,'Extrema');

Extrems=STATS(1).Extrema;

 

l1=sqrt((Extrems(3,1)-Extrems(1,1))^2+(Extrems(3,2)-Extrems(1,2))^2);

l2=sqrt((Extrems(5,1)-Extrems(3,1))^2+(Extrems(5,2)-Extrems(3,2))^2);

l3=sqrt((Extrems(7,1)-Extrems(5,1))^2+(Extrems(7,2)-Extrems(5,2))^2);

l4=sqrt((Extrems(5,1)-Extrems(1,1))^2+(Extrems(5,2)-Extrems(1,2))^2);

l5=sqrt((Extrems(7,1)-Extrems(3,1))^2+(Extrems(7,2)-Extrems(3,2))^2);

alfa=acos((l2^2+l4^2-l1^2)/(2*l2*l4))*180/pi;

betta=acos((l2^2+l5^2-l3^2)/(2*l2*l5))*180/pi;

gamma=180-(alfa+betta);

Представленный выше алгоритм позволяет вычислить угол скрученности кости.

Исходя из этого, можно составить блок-схему алгоритма для определения длины, ширины и угла скрученности остеометрического объекта (рисунок 2.13).

Рисунок 4.9 ‑  Блок-схема алгоритма для определения параметров объекта

 

Выводы: предложенная модель оптико-физической измерительной системы поможет оптимизировать процесс остеометрических измерений, ускорив его за счет автоматизации процессов. В частности, зрение человека заменено оптическим зрением фотокамеры Sony Cyber-shot DSC-HX300, которая передает качественное, с высоким разрешением, изображение измеряемого остеометрического объекта на одноплатный компьютер Cubietruck Cubieboard3. Одноплатный компьютер, в свою очередь, одновременно заменяет собой ручные средства измерений и проводит процесс исследования объекта, определяет его параметры, позволяет передать полученные данные для дальнейшей систематизации и хранения. Сам процесс измерений реализуется в качестве алгоритма в среде Matlab и сводится к бинаризации фотоизображения, определению по нему размеров (в пикселях) прямоугольника, ограничивающего измеряемый объект, и умножению числа пикселей на известный коэффицинт k = 0,02636 см.

5 ДОКАЗАТЕЛЬСТВО РАБОТОСПОСОБНОСТИ АЛГОРИТМА ДЛЯ ОПРЕДЕЛЕНИЯ ПАРАМЕТРОВ ОБЪЕКТА

5.1 Результаты измерения параметров остеометрического объекта с использованием ручных средств

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

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

Согласно Рекомендации по метрологии Государственной системы обеспечения единства измерений «Измерения прямые однократные. Оценивание погрешностей и неопределенности результата измерений» Р 50.2.038-2004 прямые однократные измерения возможны лишь при определенных условиях:

1.     Объем априорной информации об объекте такой, что определение измеряемой величины не вызывает сомнений;

2.     Изучен метод измерения, его погрешности либо заранее устранены, либо оценены;

3.     Средства измерений исправны, а их метрологические характеристики соответствуют установленным нормам.

Все требуемые условия соблюдены. Ниже представлены результаты измерений длины, ширины и угла скрученности конкретных остеометрических объектов с помощью ручного средства измерений - измерительного штатива. Так как погрешности не избежать, необходимо провести измерения трижды, чтобы определить ошибку [16].

Пример 1. Левая бедренная кость (Рисунок 5.1). Результаты измерений представлены в Таблица 5.1.

 

Рисунок 5.1– Левая бедренная кость

 

Таблица 5.1 – Результаты измерений

1

2

3

L,мм

136

140

138

W,мм

95

98

96

γ,

22

21

23

 

Пример 2. Левая М. берцовая кость (Рисунок 5.2). Результаты измерений представлены в Таблица 5.2.

 

Рисунок 5.2– Левая М. берцовая кость

Таблица 5.2– Результаты измерений

1

2

3

L,мм

117

121

122

W,мм

14

17

16

 

5.2  Результаты измерения параметров остеометрического объекта с использованием среды Matlab

Опираясь на блок-схему алгоритма для определения параметров остеометрического объекта (рисунок 2.12) и способ нахождения параметров, основанный на  определении размеров ограничивающего прямоугольника, составим текст алгоритма в программной среде Matlab и проверим его работоспособность на реальных остеометрических объектах.

Пример 1. Левая бедренная кость.

Бедренная кость имеет три параметра измерений, это длина, ширина и угл скрученности.

Сначала обращаемся к фотоизображению измеряемого объекта, загружаем его для дальнейшей работы. Затем визуализируем его (рисунок 3.4 а). С помощью команды im2bw проводим бинаризацию изображения с определенным пороговым коэффициентом [4,5]. В текущем примере его точным значением не задавались, так как это не являлось основной задачей. Визуализируем бинаризированное изображение (рисунок 3.4 б). Если иметь в виду, что в дальнейшем этот алгоритм будет использован для работы с множеством разнообразных костей, в которых могут быть сколы, дырки и другие разрушения, необходимо предусмотреть команду для получения изображения с «залитыми» дырами. К тому же не исключена возможность появления выпавших пикселей за предел порогового значения, которые также необходимо ликвидировать. После формирования изображения с «залитыми дырами» визуализируем его (рисунок 3.4 в). Следующим этапом является определение ограничивающего прямоугольника, который строится по экстремальным точкам объекта, используя команду regionprops. Определив прямоугольник, программа выдает его координаты, высоту, ширину. Заметим, что длина и ширина  высчитывается в пикселях. Для того, чтобы выразить ее в сантиметрах, нужно умножить полученное число пикселей на размер одного пикселя, который равен k=0.02636см.

а)

б)

в)

Рисунок 5.3 – Левая бедренная кость: а) исходное изображение; б) бинаризированное изображение; в) изображение с «залитыми» дырами

 

Таким образом длина кости составляет 135,49 мм, ширина 95,65 мм.

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

а)

б)

в)

Рисунок 5.4 ‑ Левая бедренная кость: а) исходное изображение; б) бинаризиванное изображение; в) изображение с «залитыми» дырами

 

Таким образом, угол скрученности кости равен 21,746.

 

Ниже приведен текст алгоритма в программной среде Matlab:

clc

clear all;

 

L=imread('itog.bmp');

L=double(L)./255;

figure(1), imshow(L);

title('Исходное изображение');

I=L(:,:,1)+L(:,:,2)+L(:,:,3);

figure(10), subplot(211);imshow(I);title('Интенсивность');

subplot(212); imhist(I);

imwrite(I,'I.bmp');

%выделение красной составляющей

r=L(:,:,1)./(I+eps);

figure(20), subplot(211);imshow(r);title('Составляющая r');

subplot(212); imhist(r);

imwrite(r,'r.bmp');

%выделение зелёной составляющей

g=L(:,:,2)./(I+eps);

figure(30), subplot(211);imshow(g);title('Составляющая g');

subplot(212); imhist(g);

imwrite(g,'g.bmp');

%выделение синей составляющей

b=L(:,:,3)./(I+eps);

figure(40), subplot(211);imshow(b);title('Составляющая b');

subplot(212); imhist(b);

imwrite(b,'b.bmp');

 

%N - количество пикселей по высоте

%M - количество пикселей по ширине

[N M]=size(r);

 

%Выделение составляющей насыщенность

for i=1:N;

    for j=1:M;

        MINIMUM=min([r(i,j) g(i,j) b(i,j)]);

        S(i,j)=1-3*MINIMUM;

    end;

end;

S(S<0)=0;

figure(50), subplot(211);imshow(S);title('Насыщенность S');

subplot(212); imhist(S);

imwrite(S,'S.bmp');

 

%Выделение составляющей фона

F=real(acos((2*r-g-b)./(sqrt(6).*sqrt( (r-1/3).^2+(g-1/3).^2+(b-1/3).^2))));

figure(60), subplot(211);imshow(F);title('Цветовой фон F');

subplot(212); imhist(F);

imwrite(F,'F.bmp');

%clear;

% L=imread('F.bmp');

% L=double(L)./255;

%figure,imshow(L);title('Исходное изображение');

% Запись красной, зелёной и синей составляющей в матрицы R, G и B

% И последующая нормировка значений в них

R=imread('r.bmp');

R=R(:,:,1);

R=double(R)./255;

G=imread('g.bmp');

G=G(:,:,1);

G=double(G)./255;

B=imread('b.bmp');

B=B(:,:,1);

B=double(B)./255;

% Запись составляющих фона, интенсивности и насыщенности в матрицы F, I и S

% И последующая нормировка

F=imread('F.bmp');

F=F(:,:,1);

F=double(F)./255;

I=imread('I.bmp');

I=I(:,:,1);

I=double(I)./255;

S=imread('S.bmp');

S=S(:,:,1);

S=double(S)./255;

 

[N M]=size(R);

%Пороговая обработка цветного изображения

for i=1:N;

    for j=1:M;

         if  (I(i,j)<0.2)

            L(i,j,1)=1;L(i,j,2)=1;L(i,j,3)=1;

        else L(i,j,1)=0;L(i,j,2)=0;L(i,j,3)=0;

        end;

    end;

end;

figure(2),imshow(L);title('Результирующее изображение');

% Бинаризация изображения

level=graythresh(L);

 bw=im2bw(L,level);

 figure(3), imshow(bw);

% Заполнение пробелов

 STATS=regionprops(bw,'FilledImage');

Image=STATS(1).FilledImage;

 figure(4), imshow(Image);

 imwrite(Image,'hg.bmp');

 

%Запись ширины и высоты изображения в переменную Boundings

% STATS=regionprops(Image,'BoundingBox');

% Boundings=STATS(1).BoundingBox;

% Length=Boundings(1,4); %длина кости

% Width=Boundings(1,3); %ширина кости

%Вычисление угла скрученности по экстремумам изображения (по спицам)

% Выделение экстремумов на изображении

STATS=regionprops(Image,'Extrema');

Extrems=STATS(1).Extrema;

l1=sqrt((Extrems(3,1)-Extrems(1,1))^2+(Extrems(3,2)-Extrems(1,2))^2);

l2=sqrt((Extrems(5,1)-Extrems(3,1))^2+(Extrems(5,2)-Extrems(3,2))^2);

l3=sqrt((Extrems(7,1)-Extrems(5,1))^2+(Extrems(7,2)-Extrems(5,2))^2);

l4=sqrt((Extrems(5,1)-Extrems(1,1))^2+(Extrems(5,2)-Extrems(1,2))^2);

l5=sqrt((Extrems(7,1)-Extrems(3,1))^2+(Extrems(7,2)-Extrems(3,2))^2);

alfa=acos((l2^2+l4^2-l1^2)/(2*l2*l4))*180/pi;

betta=acos((l2^2+l5^2-l3^2)/(2*l2*l5))*180/pi;

gamma=180-(alfa+betta);

 

Пример 2. Рассмотрим левую берцовую кость. Ниже представлены результаты обработки изображения.

На Рисунок 5.5 представлено исходное изображение, по которому был произведен расчет длины и ширины кости.

а)

б)

в)

Рисунок 5.5 – Левая М. берцовая кость: а) исходное изображение; б) бинаризированное изображение; в) изображение с «залитыми» дырами

 

Таким образом, длина кости 116,775 мм, ширина составила 14,893мм.

5.3 Оценка точности автоматизированных измерений и сравнение с ручными измерениями

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

Методы оценивания характеристик погрешности результатов измерений регламентируются документом ГОСТ 8.207-76 «Государственная система обеспечения единства измерений» [6].

Все погрешности подразделяются на три вида: систематические, случайные и приборные (Рисунок 5.6).

Рисунок 5.6 – Виды погрешностей

 

Систематические погрешности обусловлены факторами, действующими одинаковым образом при повторении одноплановых измерений. Например:

– неисправностью приборов или их неправильным использованием;

– несовершенством методики измерения или неучётом постоянно действующих факторов, влияющих на исследуемое явление;

– применением приближённых («упрощённых») формул в случае косвенных измерений.

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

– несовершенство органов чувств экспериментатора;

– несовершенство измерительного прибора;

– случайные неконтролируемые изменения внешних воздействий – температуры, влажности, давления; – которые практически невозможно учесть.

Приборными называются погрешности, обусловленные принципиальным несовершенством и ограниченной точностью изготовления и градуировки технических устройств, используемых для измерений. Приборная погрешность определяется для всей совокупности конкретного типа приборов путём сравнения показаний приборов исследуемой партии с показаниями эталонного прибора. За значение приборной погрешности принимается предельно допустимая для данного класса приборов погрешность, которая указывается в паспортных данных прибора.

Задача обработки результатов измерений заключается в том, чтобы определить границы, в которых заключено истинное значение измеряемой величины – доверительный интервал [17].

Принята следующая форма записи результата измерений какой-либо величины:

5.1

Где [x] – наиболее вероятное значение измеряемой величины x, а Δx – рассчитанная погрешность эксперимента [8].

Сравним погрешности измерений для определения длины кости.

Пример 1. В ходе проведенных измерений остеометрического объекта (Рисунок 5.1) с помощью измерительного штатива были получен набор значений (Таблица 5.1):

136 мм, 141 мм, 140 мм.

За наиболее вероятное значение искомой величины [x] принимается среднее арифметическое значение результатов измерений:

5.2

Погрешность измерений Δx оценивается следующим образом. Вычисляются частные отклонения отдельных измерений Δxi:

,

.

5.3

Оценивается абсолютная погрешность измерений Δx:

.

5.4

Полезно также определить относительную погрешность измерений ε:

.

5.5

Относительная погрешность наглядно характеризует точность проведенных измерений.

Что касается приборной погрешности, то она определяется классом точности для стрелочных приборов. А при определении приборной погрешности некоторых «простейших» приборов, не имеющих паспорта, можно ориентироваться на цену наименьшего деления шкалы: она обычно согласована с классом точности прибора. Так, при измерениях длины линейкой с «миллиметровой шкалой» приборную погрешность принимают равной 1 мм (Таблица 5.3).

Таблица 5.3– Приборная погрешность

Вид прибора

Значение
погрешности, мм

1

Миллиметровые линейки

1

2

Штангенциркули (с числом делений нониуса – 10)

0,1

3

Штангенциркули (с числом делений нониуса – 20)

0,05

4

Микрометры

0,01

 

Таким образом:

.

5.6

Для определения доверительного интервала необходимо учесть как случайную погрешность измерения, так и погрешность, вносимую приборами.

Результирующую погрешность будем называть полной погрешностью эксперимента

Для оценки полной погрешности эксперимента при прямых измерениях складывают погрешности измерений и приборные погрешности:

5.7

Полная относительная погрешность равна:

.

5.8

Таким образом, результат измерений в виде доверительного интервала:

5.9

Аналогично рассчитаем величины погрешностей для двух других примеров. Получаем:

Пример 2. Исходя из данных Таблица 5.2:

117 мм, 121 мм, 122 мм.

.

.

5.10

Теперь оценим точность проведенных автоматизированных измерений с помощью алгоритма математической среды Matlab.

Для определения погрешностей, необходимо предусмотреть допущение. Так как при реализации алгоритма в математической среде Matlab не изменяются условия проведения измерений, ракурс фотографирования, а также текст программы и условия ее работы, то значения, получаемые в ходе измерения одного объекта так же не изменятся, независимо от числа повторений опыта. Поэтому, чтобы получить набор значений для оценки точности измерений, необходимо найти ряд из пиксельных значений длины кости. Добиться этого можно путем изменения значения порогового коэффициента бинаризации. Его изменение на 0,1 приводит к изменению длины приблизительно на 5-6 пикселей. А это, в свою очередь, влечет за собой изменение конечного результата. Так как размер пикселя очень мал, можно пользоваться этим допущением и проводить дальнейшую обработку результатов.

Пример 1. Возьмем ряд пиксельных значений длины бедренной кости: 514, 517 и 521. Этому ряду соответствует набор величин:

135,49 мм, 136,28 мм, 137,33 мм.

.

.

5.11

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

5.12

Пример 2. Возьмем ряд пиксельных значений длины для малой берцовой кости: 443, 439 и 446. Этому ряду соответствует набор величин:

116,77 мм, 115,72 мм, 117,57 мм.

.

.

.

5.13

Вычислив погрешности ручных измерений  и измерений с помощью автоматизированной системы, проведем сравнение результатов, представленных в Таблица 5.4

Таблица 5.4– Сводные данные погрешностей измерений

         Способ измерения

Параметр

1

2

ручной

автоматизированный

ручной

автоматизированный

Длина
L, мм

136

135,49

117

116,77

141

136,28

121

115,72

140

137,33

122

117,57

Среднее значение
[
x], мм

139

136,37

120

116,69

Погрешность измерений Dxизм, мм

2

0,64

2

0,64

Относительная погрешность измерений
eизм, %

1,439

0,47

1,667

0,53

Приборная погрешность
Dxпр, мм

1

-

1

-

Относительная приборная погрешность
eпр, %

0,719

-

0,833

-

Полная погрешность
Dx, мм

3

0,64

3

0,64

Полная относительная погрешность
e, %

2,16

0,47

2,5

0,53

 

Выводы: результаты реализации алгоритма для определения параметров остеометрического объекта в программной среде Matlab на конкретных примерах подтвердили его работоспособность. Оценка точности и погрешностей измерений показала, что автоматизированный способ измерения дает более точный результат, до сотых долей миллиметра, и приблизительно в семь раз меньшие погрешности, что доказывает предпочтительность использования алгоритма Matlab.

ЗАКЛЮЧЕНИЕ

В работе рассмотрено достижение цели исследования метода автоматизации остеометрических измерений.

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

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

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

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


 

СПИСОК ИСПОЛЬЗОВАННЫХ ИСТОЧНИКОВ

1.     Алексеев В.П. Остеометрия. Методика антропологических исследований// М.: Наука,1966.

2.     Петри Э. Ю. Антропология. Соматическая антропология // СПб.,1895-1897г.

3.     Шандров Б.В., Чудаков А.Д. Технические средства автоматизации // М.:Академия, 2010, 368с.

4.     Алпатов Б.А., Бабаян П.В., Балашов О.Е., Степашкин А.И. Методы автоматического обнаружения и сопровождения объектов. Обработка изображений и управление. М., Радиотехника, 2008. 176 с

5. Гансалес Р., Вудс Р. Цифровая обработка изображений // М.: Техносфера, 2005.1012 с.

6.     Яне Б. Цифровая обработка изображений. М.: Техносфера, 2007. 583 с.

7. Бакут П. А., Колмогоров Г.С., Ворновицкий И. Э. Сегментация изображений: методы пороговой обработки. // Зарубежная радиоэлектроника. № 10, 1987. С. 6–24

8.         Журавель И. М. Краткий курс теории обработки изображений. Анализ признаков объектов. http://matlab.exponenta.ru/imageprocess/book2/60.php

9. Сойфер В. А. Методы компьютерной обработки изображений. М., ФИЗМАТЛИТ, 2003. 784 с.

10.       Гансалес Р., Вудс Р., Эддинс С. Цифровая обработка изображений в среде MATLAB // М.: Техносфера, 2006. 615 с.

11.  Трубникова Т.А., Гудинов К.К., Двуреченский С.А., Гусев В.П. Цифровая Фотоаппаратура: учебное пособие. Санкт-Петербург: СПБГУКиТ, 2010, 72 с.

12.       Тымкул В.М, Тымкул Л.В. Оптико-электронные приборы и системы. Теория и методы энергетического расчета: учебное пособие. Новосибирск: СГГА., 2005. 215 с.

13.       Якушенков Ю.Г. Теория и расчёт оптико-электронных приборов: учебник. М.: Логос, 2011. 567 с.

14.  Sun-store.ru. Мини-компьютер Cubietruck 8gb. http://sun-store.ru/2128--mini-kompyuter-cubietruck-8gb

15.       Дьяконов В.П. MATLAB R2007/2008/2009 для радиоинженеров // М.: МДК-Пресс, 2010, 976с.

16.       ГОСТ 8.207-76 «Государственная система обеспечения единства измерений»

17.       Ефимова А.И., Зотеев А.В., Склянкин А.А. Общий физический практикум физического факультета МГУ. Погрешности эксперимента: Учебно-
ме­тоди­­­ческое пособие // М.: МГУ, Физический факультет, 2012. – 39 с.

18.       Шалаткина А. Н., Борисова О. И., Корниенко К. Ю. Систематизация и автоматизация остеометрических измерений.  - Математическая морфология. Электронный математический и медико-биологический журнал. - Т. 12. - Вып. 4. - 2013. - URL:

http://www.smolensk.ru/user/sgma/MMORPH/N-40-html/shalatkina/shalatkina.htm

http://www.smolensk.ru/user/sgma/MMORPH/TITL.HTM

http://www.smolensk.ru/user/sgma/MMORPH/N-40-html/TITL-40.htm

http://www.smolensk.ru/user/sgma/MMORPH/N-40-html/cont.htm

 

 

Поступила в редакцию 25.06.2015.