Министерство образования и науки Российской Федерации

УДК

ГРНТИ

Инв. №

УТВЕРЖДЕНО:

Исполнитель:

Государственное учебно-научное учреждение

Факультет вычислительной математики и

кибернетики Московского государственного

университета имени М.В.Ломоносова

От имени Руководителя организации

/_/

М.П.

НАУЧНО-ТЕХНИЧЕСКИЙ

ОТЧЕТ о выполнении 1 этапа Государственного контракта № 14.740.11.0996 от 23 мая 2011 г.

Исполнитель: Государственное учебно-научное учреждение Факультет вычислительной математики и кибернетики Московского государственного университета имени М.В.Ломоносова Программа (мероприятие): Федеральная целевая программа «Научные и научнопедагогические кадры инновационной России» на 2009-2013 гг., в рамках реализации мероприятия № 1.2.2 Проведение научных исследований научными группами под руководством кандидатов наук.

Проект: Математические методы анализа и обработки стохастических изображений.

Руководитель проекта:

/Шестаков Олег Владимирович (подпись) Москва 2011 г.

СПИСОК ОСНОВНЫХ ИСПОЛНИТЕЛЕЙ

по Государственному контракту 14.740.11.0996 от 23 мая 2011 на выполнение поисковых научно-исследовательских работ для государственных нужд Организация-Исполнитель: Государственное учебно-научное учреждение Факультет вычислительной математики и кибернетики Московского государственного университета имени М.В.Ломоносова Руководитель темы:

Шестаков О. В.

кандидат физикоматематических наук, подпись, дата доцент Исполнители темы:

Мизотин М. М.

кандидат физикоматематических наук, без подпись, дата ученого звания Насонов А. В.

без ученой степени, без ученого звания подпись, дата Сорокин Д. В.

без ученой степени, без ученого звания подпись, дата Черноморец А. А.

без ученой степени, без ученого звания подпись, дата Семашко А. С.

без ученой степени, без ученого звания подпись, дата Сторожилова М. В.

без ученой степени, без ученого звания подпись, дата Сергеев В. В.

без ученой степени, без кандидат физикоматематических наук, без подпись, дата ученого звания без ученой степени, без без ученой степени, без Отчет 76 с., 1 ч., 1 рис., 0 табл., 55 источн., 0 прил.

стохастические изображения, линейный однородный фильтр, подавление шума, регуляризирующие методы В отчете представлены результаты исследований, выполненных по 1 этапу Государственного контракта № 14.740.11.0996 "Математические методы анализа и обработки стохастических изображений." (шифр "2011-1.2.2-111от 23 мая 2011 по направлению "Проведение научных исследований научными группами под руководством кандидатов наук в следующих областях:- математика;- механика" в рамках мероприятия 1.2.2 "Проведение научных исследований научными группами под руководством кандидатов наук.", мероприятия 1.2 "Проведение научных исследований научными группами под руководством докторов наук и кандидатов наук", направления 1 "Стимулирование закрепления молодежи в сфере науки, образования и высоких технологий." федеральной целевой программы "Научные и научнопедагогические кадры инновационной России" на 2009-2013 годы.

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

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

Для выполнения работ первого этапа были использованы материалы электронных библиотек издательств, подписка на которые осуществляется через РФФИ. Для проведения вычислительных экспериментов задействована вычислительная и программная база кафедры математической статистики и лаборатории математических методов обработки изображений факультета ВМК МГУ имени М.В. Ломоносова.

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

-аналитический обзор;

-новые постановки задач анализа и обработки стохастических сигналов и изображений;

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

СОДЕРЖАНИЕ

Введение

Проведение I этапа исследований по проблеме «Математические методы анализа и обработки стохастических изображений».................. оптимального направления исследований при разработке методов анализа и обработки стохастических сигналов и изображений, рассматриваемых в проекте

Разработка плана проведения теоретических исследований и вычислительных экспериментов на ЭВМ

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

Заключение

Список использованных источников

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Ожидаемые результаты На первом этапе будут получены следующие результаты:

Аналитический обзор.

Новые постановки задач анализа и обработки стохастических сигналов и изображений.

пропущенных через линейный однородный фильтр.

На втором этапе будут получены следующие результаты:

сферических преобразований Радона.

Метод регуляризации для реконструкции и оценки погрешностей Устойчивый метод определения поверхностного натяжения жидкости методом лежащей капли.

На третьем этапе будут получены следующие результаты:

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

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

использующие системы ортогональных многочленов.

На четвертом этапе будут получены следующие результаты:

Методы построения оценок погрешностей при пороговой обработке вейвлет-коэффициентов.

Метод повышения качества изображений на основе метода суперразрешения изображений.

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

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

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

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

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

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

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

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

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

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

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

приравниваются к нулю (режектируются) при реконструкции сигнала.

Критерии выбора порогов могут быть различны. В общем случае пороги могут выбираться адаптивно на основе статистического анализа сигнала.

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

Пусть {g m }m — ортонормированный базис гильбертова пространства H. Рассмотрим класс сигналов, которые разложимы по первым N векторам базиса.

Пусть f — сигнал длины N. Это значит, что данный сигнал f определяется набором чисел f [0], f [1],, f [ N 1]. Обозначим через W возникающие помехи или шум, которые являются реализацией некоторого случайного процесса, тогда искаженный сигнал будет определяться как = {g m }0m< N соответственно Скалярное произведение левой и правой частей равенства X = f W с g m дает Предположим, что W — гауссовский белый шум с нулевым средним значением и дисперсией 2.

диагонального оператора D в базисе. Диагональный оператор оценивает каждую компоненту сигнала f [m] по соответствующей зашумленной компоненте X [m] с помощью функций d m (x) :

Запишем функции d m в виде где a[m] зависит от X [m].

Обозначим через r ( D, f ) риск оценки F, вычисляемый по формуле Для линейного оператора D (т.е. когда a[m] не зависит от X [m] ) это выражение принимает наименьшее значение при и минимальный риск равен В случае, когда D нелинейный, а коэффициенты a[m] принимают значения из множества {0,1}, минимум выражения rinf ( f ) достигается при При этом минимальный риск равен В обоих случаях оператор D не может использоваться практически, так как a[m] зависит от | f [ m] |, значение которого не известно.

Функции d m будем называть пороговыми, если обработка сигнала этими функциями зависит от сравнения амплитуды сигнала с некоторой константой (порогом).

Жесткой пороговой обработкой (с жестким порогом T) называется процедура оценивания сигнала, при которой пороговые функции d m определяются как Риск жесткой пороговой обработки равен Следовательно, Мягкой пороговой обработкой (с мягким порогом T) называется процедура оценивания сигнала, при которой пороговые функции d m определяются как Введем функцию где X - нормально распределенная случайная величина со средним значением и дисперсией 1, а 1A - индикатор события A.

Так как f [m] - константа, то зашумленные компоненты сигнала X [m] = f [m] W [m] являются гауссовскими случайными величинами со средним значением f [m] и дисперсией 2. Следовательно, риск мягкой пороговой оценки F с порогом T определяется формулой Отметим, что при мягкой пороговой обработке пороговые функции обладают свойством непрерывности, а при жесткой — пороговые функции разрывны. Отсюда вытекают и различия в полученных восстановленных сигналах.

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

Метод пороговой обработки, как показано в [6] и [10], является близким к оптимальному (в смысле минимизации риска) для широкого класса сигналов, зашумленных аддитивным гауссовским белым шумом.

Одним из важнейших показателей качества аудиоаппаратуры, в том числе и звуковых карт, является отношение сигнал/шум (ОСШ, англ. SNR, Signal-to-Noise Ratio). Эта величина равна отношению мощности полезного сигнала к мощности шума. Чем больше это отношение, тем менее заметен шум.

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

1.2 Методы выбора порога 1.2.1 Универсальный порог В 1993 году Донохо (D. Donoho) и Джонстон (I. Johnstone) предложили использовать порог, основанный на следующих результатах.

Пусть случайная величина имеет нормальное распределение N (a, 2 ), и m — медиана абсолютного отклонения от ее математического ожидания, тогда где C3/4 — 3/4-квантиль стандартного нормального распределения.

Приближенное значение квантили C3/4 приведено в стандартных таблицах и обычно берется равным 0,6745. Квантиль C3/4 является медианой случайной величины | z |, где z N (0,1).

При наличии выборки устойчивой оценкой стандартного отклонения в случае нормального распределения является (с некоторой погрешностью округления квантили C3/4 ) величина где MAD (Median Absolute Deviation) — выборочное медианное абсолютное отклонение (см. [7]).

стандартное нормальное распределение и Более того, если Для нашей модели zi N (0, 2 ), поэтому события An и Bn (t ) будут иметь вид:

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

Порог T = 2 ln n во многом решает эту проблему. Согласно лемме 3.2, с большой вероятностью происходит удаление основного шума, так как вероятность An стремится к нулю, а оставшийся шум в обработанном экспоненциально убывает. Неизвестный параметр можно заменить его оценкой.

Эти рассуждения объясняют выбор порога названного универсальным.

При выборе универсального порога T U риск и мягкой, и жесткой Риск rt ( f ) мягкой или жесткой пороговой оценки удовлетворяет при всех N 4 неравенству 1.2.2 SURE-порог Риск обработки сигнала можно уменьшить, выбирая изменяющийся адаптивный порог, т.е. зависящий от зашумленных данных X. Несмещенные оценки риска принято называть SURE-оценки, по аббревиатуре слов Stein Unbiased Risk Estimator. Порог, минимизирующий такую оценку при заданных зашумленных данных, также называют SURE-порогом.

Исследуем влияние порога на риск. Обозначим через rt ( f, T ) риск мягкой пороговой обработки с порогом T, который зависит от неизвестных значений f [m], и через rt ( f, T ) — оценку этого риска, зависящую от наблюдаемых зашумленных данных X. Напомним, что разложения сигнала X по базису. Заметим, что если | X [m] | T, то dT ( X [m]) = 0, следовательно вклад в риск m -го слагаемого равен | f [m] |2.

Это число можно оценить разностью | X [m] |2 2, так как равен Следовательно, риск мягкой пороговой обработки можно оценить следующим образом:

Определенная выше оценка риска rt ( f, T ) является SURE-оценкой: при мягкой пороговой обработке оценка риска rt ( f, T ) не смещена, т.е.

пересортированного (ранжированного) вектора. Так, что | X [k ] || X [k 1] | Пусть l — такой индекс, что принадлежит.

4. Из всех точек QRS-комплекса, определенных в пункте 3, выбрать такие, которые образуют один и тот же QRS-комплекс.

С этой целью с параметрами 2 = 0,1 и t =0,0005 сек., определенными эмпирическим путем, применить следующее правило: выбрать ближайшие по значению вейвлет-коэффициенты d j,n и d j,m, соответствующие отсчетам n и m, и если то n -ый и m -ый отсчеты принадлежат одному и тому же QRSкомплексу, иначе — разным.

5. Вычислить расстояние между соседними QRS-комплексами и определить число обнаруженных QRS-комплексов.

Анализ ритмограммы Алгоритм обнаружения QRS-комплексов использутеся для построении так называемой ритмограммы. Ритмограмма RR -интервалов как мера вариабельности сердечного ритма (ВСР) получила в последние годы достаточно широкое распространение в диагностике сердечно-сосудистых заболеваний (ССЗ) (см. [26] и [27]). (Под RR интервалами подразумеваются временные промежутки между пиками соседних R -зубцов.) Благодаря развитию компьютерной техники стал возможным анализ десятков тысяч записей и выявление общих закономерностей. В частности, оказалось, что существует зависимость между ССЗ и спектральными характеристиками ритмограммы.

Существует ряд рекомендаций относительно построения ритмограммы.

Они устанавливают минимальную продолжительность измерений, при которой, во-первых, возможно вычислить спектр для определенного диапазона частот, и, во-вторых, спектральные характеристики будут постоянны [26]. Однако в ритмограмме обычно присутствуют так называемые эктопические импульсы - очень короткие или очень длинные интервалы между биениями. Обычно сначала идет короткий RR интервал, а сразу после него - длинный. Эти импульсы имеют иную природу, в отличие от нормальных (или, как их еще называют, синусовых) интервалов [28].

Наличие всего нескольких эктопических импульсов оказывает существенное влияние на спектр ритмограммы, поэтому обычно их стараются удалить. Стоит отметить, что ритмограмма относится к сигналам с неравномерным квантованием времени, так как длина отсчета времени равна длине очередного интервала. Это необходимо учитывать при вычислении спектральных характеристик.

Весьма популярный метод удаления эктопических импульсов состоит в следующем:

1. Удалить все интервалы, короче 300мс и длиннее 2000мс.

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

3. Удалить интервалы, отличающиеся по длине более чем на 20% от среднего значения, где среднее вычисляется по пяти предыдущим синусовым интервалам.

Критические значения вычислены по базе ЭКГ Массачусетского технологического института (MIT-BIH Normal Sinus Rhythm Database), насчитывающей десятки тысяч записей. Эти значения являются некоторым усреднением, и поэтому вполне может случиться так, что какие-то эктопические значения не будут отсеяны и останутся в сигнале. При таком подходе фактически не учитывается форма (структура) ритмограммы и индивидуальные особенности (амплитуда колебаний). В предлагаемом методе присутствует адаптация под конкретный сигнал.

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

На практике для спектрального анализа обычно используются ритмограммы двух типов: короткие (5 минут) и длинные (24 часа). Частотные характеристики пятиминутных ритмограмм можно считать постоянными при выполнении ряда условий, тогда как измерения, полученные за 24 часа, таким свойством не обладают. При этом длинные записи несут более полную спектральную информацию. Далее будут рассматриваться только короткие ритмограммы.

При работе с пятиминутными ритмограммами выделяют 2 диапазона частот: диапазон высоких частот ВЧ ( 0,15 0.4 Гц) и диапазон низких частот НЧ ( 0.04 0.15 Гц). Анализ данных показал, что в этих диапазонах часто можно четко выделить два биоритма: волну Майера (Mayer wave) в НЧ области и дыхательную волну (Respiratory Sinus Arrhythmia wave) в ВЧ области. Волна Майера соответствует колебаниям кровяного давления, а дыхательная волна отражает влияние дыхания на сердечный ритм.

Таким образом, в спектре ритмограммы имеются два подъема с максимумами, по одному в НЧ и ВЧ диапазонах. Это наталкивает на мысль о представлении ритмограммы в виде суммы двух синусоид с круговыми частотами, соответствующими пикам спектра в НЧ и ВЧ областях. Таким образом, получаем представление:

RR (t ) - ритмограмма;

1 = 2f1 и 2 = 2f 2 — круговые частоты, соответствующие частотам A1 и A2 — амплитуды, а 1 и 2 — фазы синусоид;

(t ) — шум, нормально распределенные случайные величины с нулевым средним и дисперсией 2 ;

C — некоторая константа, выполняющая роль сдвига графика ритмограммы относительно нуля.

Время меняется дискретно: t {t1, t2,, t N |0 = t1 < t2 < < t N }. Следует отметить, что данная модель является весьма грубой, так как форма реальной ритмограммы несколько сложнее, чем сумма двух синусоид. Однако она вполне достаточна, ведь стоит задача оценки амплитуды синусовой части ритмограммы, а синусовая часть является довольно «гладкой» в смысле малых изменений в соседних точках, в то время, как эктопические импульсы соответствуют скачкам на графике ритмограммы.

Идея предлагаемого метода состоит в следующем:

1. Оценить параметры 1, 2, 1, 2, A1 и A2 по спектру с помощью робастного варианта регрессии.

4. Отсеять импульсы, не попадающие в доверительные интервалы для RRh (t ) A1 sin (1t 1 ) A2 sin (2t 2 ) будет выполнено:

При построении доверительных интервалов предполагается, что дисперсия шума 2 известна. Тогда нормальная случайная величина с нулевым средним и дисперсией 2 с вероятностью 1 попадает в интервал, 1 1, где (x) суть функция распределения стандартной Доверительные интервалы для разностей RR (t ) будут выглядеть так:

Однако необходимо оценить. Для этого можно использовать разности RRi = RR (ti ) RR (ti 1 ) функции RR (t ), исключив некоторым образом из этих разностей влияние полезного сигнала. Предлагается делать это в соответствии с моделью: в качестве шума рассматривать значения Коэффициенты RR h оцениваются методом робастной регрессии.

В качестве оценки для медианное отклонение от медианы * :

Эта оценка является наиболее робастной, ее пороговая точка (раздел 5) равна 0.5. Кроме того, * будет состоятельной оценкой стандартного отклонения нормальной случайной величины. Таким образом, осталось вычислить 1, 2, 1, 2, A1 и A2.

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

наименьших квадратов строится следующим образом. Необходимо оценить y = ( y1, y2,, yN )T, причем имеет место связь:

где xij — некоторые известные постоянные коэффициенты, ui — независимые одинаково распределенные случайные величины, i = 1,, N, j = 1,, p.

Эта задача сводится к решению задачи:

Если матрица X имеет ранг p, то решение системы будет выглядеть так:

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

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

Ее производная равна Вблизи нуля функция имеет квадратичный порядок роста, а начиная с некоторого значения - нулевой. В качестве c берется значение c = 0.85, которое было подобрано опытным путем. Полезным свойством данной функции является то, что очень большие по модулю значения аргумента не оказывают влияния на сумму, так как значение функции в этих точках равно нулю. Это полностью согласуется с практикой, которая говорит, что очень длинные ( 2000 мс) и очень короткие интервалы ( 300 мс) однозначно являются эктопическими. Именно такие интервалы дают большие абсолютные значения разностей.

В робастном варианте регрессии возникает задача оценки масштаба, т.к. функции и практически всегда имеют разный порядок роста для разных значений аргумента. Соответственно, параметр масштаба s позволяет выполнить нормировку данных для эффективного деления на большие и маленькие значения. Нужно решать систему Значения s и находятся итеративным путем, но, так как функции x являются линейными относительно j, то решение находится за один шаг. В качестве начального значения параметра масштаба ~ можно взять абсолютное медианное отклонение от медианы (AMO) которое обладает хорошими теоретическими свойствами, в частности, AMO имеет пороговую точку * = 0.5. На интуитивном уровне понятие пороговой точки можно объяснить как максимальную долю паразитных импульсов, с которыми может справиться оценка.

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

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

Такой вид линейной регрессии называется компенсированным дискретным преобразованием Фурье (КДПФ) (date-compensated discrete Fourier transform) [32]. Дело в том, что неравномерное квантование времени, а именно такое ортогональность и нормировку тригонометрических векторов X 2, X 3 и единичного вектора X1 :

неправильно, причем в отдельных случаях отличия от реальных значений амплитуд весьма существенные [33].

КДПФ решает эту проблему. К трем векторам применяется процедура ортогонализации Грама-Шмидта:

Коэффициенты a1, a 2 и a3 находятся из условия ортонормированности векторов x1, x 2 и x 3. Здесь частота фиксирована, и p = 3. Поэтому коэффициенты регрессии выглядят так:

винзоризованных остатков r, получим, соотвественно, векторы коэффициентов 0 и. Далее находим значение амплитуды для частоты :

максимумы в НЧ и ВЧ:

[0,2 ). (Здесь supp(G) - носитель G ( x, y ), а supp( g ) - носитель g (s ) ).

{( x, y ) : x 2 y 2 a}, используя R f (s ) только при s [a, a]. Предполагая, что и положив G ( x, y ) = 2G( x/, y/ ), мы могли бы аппроксимировать f ( x, y ) при {( x, y ) : x 2 y 2 a} сколь угодно точно, так как свертку выбирая соответствующее > 0.

Однако g (s ) и G ( x, y ) не могут одновременно иметь компактный носитель по следующей причине. Предположим, например, что G ( x, y ) имеет компактный носитель, тогда хорошо известно, что G(1, 2 ) можно аналитически продолжить до целой функции в C 2, поэтому для любого [0,2 ) функция G( cos, sin ) представляет сужение на R целой функции в C. Следовательно, для того чтобы функция g (s ) имела компактный носитель, необходимо, чтобы функция G ( cos, sin ) представляла собой сужение целой функции на R. Но это невозможно, так как предполагается, что Следовательно, функция g (s ) не может иметь компактный носитель, и локальная реконструкция невозможна.

Хотя компактность носителя функции g (s ) невозможно обеспечить, если компактный носитель имеет G ( x, y ), можно обеспечить быстрое убывание g (s ) на бесконечности. Этого можно достичь требованием, чтобы функция G ( cos, sin ) была как можно более гладкой. Поскольку G ( cos, sin ) бесконечно дифференцируема везде за исключением точки = 0, для достижения желаемой гладкости необходимо, чтобы G(1, 2 ) имела ноль достаточно высокого порядка в точке (0,0). Это означает наличие у G ( x, y ) достаточно большого числа нулевых моментов.

Приведенные выше рассуждения показывают, что самое большее, на что можно надеяться, это "почти" локально восстановить f * G( x, y), где G ( x, y ) имеет достаточно большое количество нулевых моментов. Выбирая G ( x, y ) соответствующим образом, мы приходим к вейвлет-преобразованию функции f ( x, y ).

условию Определим четную функцию ( s) = (, s) с помощью следующего соотношения:

(в силу круговой симметрии эта функция не зависит от ), тогда функция (s) представляет собой допустимую вейвлет-функцию, и для вычислить, применив оператор обратного проецирования к вейвлетпреобразованиям проекций R f (s ), используя допустимую вейвлет-функцию для каждого [0,2 ).

Справедливо также следующее утверждение.

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

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

Если функция h(s ) удовлетворяет следующим условиям:

1. h(s ) имеет компактный носитель, производная h( ), и Первая формула показывает скорость убывания функции I 1[h]( s) для больших значений s, а вторая формула дает усредненную оценку убывания и показывает, что "энергия" I 1[h]( s) мала вдали от начала координат.

кратномасштабном анализе томографического изображения. Аппроксимацию функции изображения можно представить в виде а (x) и (x) — масштабирующая функция и вейвлет-функция, ассоциированные с одномерным кратномасштабным анализом (см [2]). При этом Пусть ( x, y ) и [l ] ( x, y), l = 1,2,3, - соответственно масштабирующая анализом. Определим для каждого [0,2 ) функции o (s ) и l ] ( s ), l = 1,2,3, с помощью следующих соотношений:

Важным обстоятельством является тот факт, что если вейвлет-функция (x) имеет компактный носитель и достаточно большое число нулевых моментов, то функции l ] ( s ), l = 1,2,3, [0,2 ), существенно отличны от нуля только в области примерно такого же радиуса, что и радиус носителя фиксированных m и n коэффициенты b[jl,m,n можно вычислить достаточно точно, используя значения R f (s ) при [0,2 ) и s из интервала радиуса 2 j T с центром в точке 2 j (m cos n sin ), где T - радиус области, в которой функции [l ] ( x, y) существенно отличны от нуля. Эти значения R f (s ) проходящим в круге радиуса 2 j T c центром в точке (2 j m,2 j n). Радиус круга становится меньше с ростом j, т.е. для восстановления более мелких деталей изображения требуется более локальная информация, чем для обстоятельство позволяет значительно снизить дозу облучения в процессе получения проекционных данных.

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

Классический подход к преодолению этой трудности заключается в применении линейной регуляризации по Тихонову (см., например, [24]).

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

Для того чтобы решить эту проблему, был разработан метод нелинейной регуляризации, основанный на пороговой обработке вейвлет-коэффициентов (см. [8]).

следующих соотношений:

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

убедиться в том, что Таким образом, функции o представляют из себя сдвиги и сжатия "базовых" функций o0,0,0 (, s ) и 0,0,0 (, s ), как и в случае вейвлет-базисов. Также можно убедиться в том, что (см. [4] и [8]) [jl,]m,n (, s ), растут со скоростью 2 j/2. Это и характеризует некорректность задачи реконструкции.

В практической ситуации проекционные данные регистрируются для дискретных отсчетов { r }r =1 и {st }tN 1 по переменным и s. Не ограничивая томографическое изображение, является круг единичного радиуса с центром в начале координат. Таким образом, [0, ) и s [1,1]. Рассмотрим проекционные данные в виде где прямая L(, s), как и раньше, задается уравнением а nr,t - независимые случайные величины, имеющие нормальное распределение с нулевым средним и дисперсией 2. Для регуляризации метода реконструкции воспользуемся методом пороговой обработки с универсальным порогом.

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

Из такого определения следует, что коэффициенты распределены нормально со средними, приблизительно равными и дисперсиями, приблизительно равными Последовательность функций образует в пространстве L2 (T ), где T = [0, ) R, устойчивый базис (базис Рисса) (см., например, [8]). Как показано в [4], это означает, что эмпирические коэффициенты можно считать практически независимыми.

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

Увеличение уровня шума указывает на то, что использование является адекватным. Вместо этого следует использовать свой порог на каждом масштабе. В частности, если для описания изображения на каждом j для каждого вида вейвлет-функций используется 2 j 2 j масштабе коэффициентов (как, например, в [8]), то, пользуясь методикой получения универсального порога, мы приходим к порогам регуляризованному варианту формулы для представления функции f :

обработки с порогами T [0] и T j[l ].

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

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

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

Таким образом, рассматривается следующая модель:

где индекс i обозначает номер отсчета измеряемого сигнала, X i наблюдаемые данные, K - некоторое линейное преобразование, f истинная (незашумленная) функция сигнала, а i - случайные погрешности измерения. Будем предполагать, что все i независимы и имеют одинаковое гауссово распределение с нулевым средним и дисперсией 2.

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

для любого x0 и любого a > 0, и рассматриваем метод представления сигнала, предложенный в работе [34] и получивший название вейглетвейвлет разложение (Vaguelette-Wavelet Примерами однородных линейных преобразований служат оператор интегрирования, преобразование Абеля и некоторые виды операторов свертки.

В работе [4] предложен метод вейвлет-вейглет разложения сигнала, идея которого заключается в представлении функции сигнала f L2 ( R) в виде ряда из сдвигов и растяжений некоторой вейвлет-функции :

ортонормированный базис в L2 (I R) ). Индекс j называется масштабом, а индекс k - сдвигом. Функция должна удовлетворять определенным требованиям (см. [11]), однако ее можно выбрать таким образом, чтобы она обладала некоторыми полезными свойствами, например имела компактный носитель, была дифференцируемой нужное число раз и имела заданное число M нулевых моментов, т.е.

Мы предполагаем выполнение некоторых условий гладкости. Будем считать, что функция Kf L2 ( R) задана на конечном отрезке [a, b] и равномерно регулярна по Липшицу с некоторым параметром > 0, т.е. (см.

[11]) существует константа L > 0 и полином Py степени n = такой, что для любого y [a, b] и любого x I R непрерывно дифференцируема ( M ), имеет M нулевых моментов и быстро убывает на бесконечности вместе со своими производными, т.е. для всех 0 k M и любого m N найдется константа C m, что при всех то найдется такая константа C > 0, что преобразование K обладает тем свойством, что функция Kf оказывается более гладкой, чем функция f. Таким образом, ограничения на гладкость f могут быть менее жесткими.

Поскольку наблюдается не сигнал f, а его линейное преобразование Kf, коэффициенты разложения вычислить напрямую нельзя. В работе [4] предлагается использовать последовательность функций j,k (получивших название «вейглеты») таких, что представляют собой сдвиги и растяжения некоторой функции (см. [4]).

При этом семейство { j,k } уже не обладает свойством ортонормированности, однако образует устойчивый базис, т.е. существуют такие константы для всех квадратично суммируемых последовательностей {c j,k }.

Иногда это свойство называют «почти ортогональностью».

В работе [34] предложен альтернативный метод, получивший название вейглет-вейвлет разложения сигнала. В этом методе в ряд по функциям j,k раскладывается не функция сигнала f, а ее линейное преобразование Kf :

Здесь для удобства мы не изменили обозначения для функций j,k, но в общем случае они отличаются от функций j,k в вейвлет-вейглет разложении.

Заметим, что если преобразование K однородно. Действительно, пусть K однородно с показателем. Положим g = Kf, тогда Функция f представляется в виде ряда где u j,k = K 1 j,k / j,k (по аналогии с предыдущим методом функции u j,k также называются «вейглетами»). Последовательность {u j,k } не образует ортонормированную систему, однако если выполнены некоторые условия гладкости на K * и K 1, а именно: существуют такие константы Ai > 0, то последовательность {u j,k } образует устойчивый базис.

При практической реализации метода вместо ряда функция Kf аппроксимируется конечной суммой следующего вида где 0,0 - так называемая масштабирующая функция, которая фактически описывает среднее значение измеряемых данных (см. [4]).

Соответственно, функция сигнала задается формулой Поскольку в измерениях присутствует шум, необходимо использовать некоторые процедуры для его удаления. Мы рассмотрим процедуру мягкой пороговой обработки с порогом T j = 2ln 2 j. В результате функция сигнала f оценивается следующим образом:

где X 0,0 представляет собой зашумленный коэффициент при K 10,0, а X W,k - зашумленные вейвлет-коэффициенты. Коэффициент X 0,0 обычно не подвергается пороговой обработке, так как входит в ту часть оценки, которая описывает среднее значение измеряемых данных, и дисперсия шума в этом коэффициенте много меньше его абсолютного значения, поэтому мы не рассматриваем его влияние, полагая, что X 0,0 = Kf, 0,0, т.е, что X 0,0 не содержит шума.

Заключение На первом этапе были получены следующие результаты:

Аналитический обзор.

Новые постановки задач анализа и обработки стохастических пропущенных через линейный однородный фильтр.

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

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

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

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

Список использованных источников 1. Berenstein C.A., Walnut D. Wavelets and local tomography//Wavelets in Medicine and Biology, A. Aldroubi and M. Unser (eds.). CRC Press. 1996.

P. 231–261.

2. Boggess A., Narkowich F. A First Course in Wavelets with Fourier Analysis. Prentice Hall, 2001.

3. Chang S.G., Yu B., Vetterli M. Adaptive Wavelet Thresholding for image Denoising and Compression//IEEE Trans. Image Processing. 2000. 9. N 9.

P. 1532–1546.

4. Donoho D. Nonlinear solution of linear inverse problems by waveletvaguelette decomposition//Appl. Comput. Harm. Anal. 1995. 2. P. 101–126.

5. Donoho D., Johnstone I. Adapting to unknown smoothness via wavelet shrinkage//J. Am. Statist. Assoc. 1995. 90. P. 1200–1224.

shrinkage//Biometrika. 1994. 81. P. 425–455.

7. Hoaglin D., Tukey J., Mosteller F. Understanding Robust and Exploratory Data Analysis. John Wiley and Sons, 1983.

8. Kolaczyk E.D. A wavelet shrinkage approach to tomographic image reconstruction//J. Am. Statist. Assoc. 1996. 91. P. 1079–1090.

9. Kozakevicius A.J., Rodrigues C.R., Nunes R.C., Guerra Gilho R. Adaptive BioMedical Engineering. 2005. P. 109–114.

10.Mallat S. A wavelet tour of signal processing, Second Edition. Academic Press, 1999.

11.Mallat S. Zero-crossings of a wavelet transform//IEEE Trans. Inform.

Theory. 1991. 37. N 4. P. 1019–1033.

12.Mallat S., Zhong S. Wavelet transform maxima and multiscale edges//Wavelets and their Applications, M.B. Ruskai et al. (eds.). Jones and Bartlett Publishers. 1992. P. 67–104.

13.Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical Recipes in C. Cambridge University Press, 1988.

14.Powell K.J., Sapatinas T., Bailey T.C., Krzanowski W.J. Application of wavelets to the Pre-processing of underwater sounds//Statistics and Computing. 1995. 5. P. 265–273.

15.Rashid-Farokhi F., Lin K.J.R., Berenstein C.A., Walnut D. Wavelet-based Multiresolution Local Tomography//IEEE Ttrans. Image Processing. 1997.

6. N 10. P. 1412–1430.

16.Stein C. Estimation of the mean of a multivariate normal distribution//Ann.

Statist. 1981. 9. N 6. P. 1135–1151.

17.Vidacovic B., Ruggeri F. BAMS method: theory and simulations//Sankhya.

Ser. B. 2001. 63. N 2. P. 234–249.

18.Vidacovic B. Statistical Modeling by Wavelets. Wiley and Sons, 1999.

19.Астафьева Н.М. Вейвлет-анализ: основы теории и примеры применения//УФН. 1996. 166. C. 11–45.

20.Баскаков С.И. Радиотехнические цепи и сигналы. М.: Высшая школа, 2000.

21.Добеши И. Десять лекций по вейвлетам. Ижевск: НИЦ "Регулярная и хаотическая динамика"$,$ 2001.

22.Лидбеттер М., Линдгрен Г., Ротсен Х. Экстремумы случайных последовательностей и процессов. М.: Мир, 1989.

23.Свешников А.Г., Тихонов А.Н. Теория функций комплексной переменной. М.: Наука, 1974.

24.Тихонов А.Н., Арсенин В.Я. Методы решения некорректных задач. М.:

Наука, 1974.

25.Троицкий И.Н. Статистическая теория томографии. М.: Радио и связь, 1989.

26.Malik M. Heart rate variability // European Heart Journal, volume 17, number 3, 1996, pp. 354–381.

27.Azuaje F., Clifford G., McSharry P. Advanced methods and tools for ECG data analysis. Artech House Publishers, 2005.

28.Clifford G., McSharry P., Tarassenko L. Characterizing artefact in the normal human 24-hour RR time series to aid identification and artificial replication of circadian variations in human beat to beat heart rate using a simple threshold // Computers in Cardiology, 2002, pp. 129–132.

29.Ивченко Г.И., Медведев Ю.И. Математическая статистика. Высш. шк., 1984.

30.Хьюбер П. Робастность в статистике. Мир, 1984.

31.Andrews D. A robust method for multiple linear regression // Technometrics, volume 16, number 4, 1974, pp. 523–531.

32.Ferraz-Mello S. Estimation of Periods from Unequally Spaced Observations // The As\-tro\-no\-mi\-cal Jour\-nal, volume 86, number 4, 1981, pp. 619– 33.Foster G. The cleanest Fourier spectrum // The Astronomical Journal, volume 109, number 4, 1995, pp. 1889–1902.

34.Abramovich F., Silverman B. W. Wavelet Decomposition Approaches to Statistical Inverse Problems // Biometrika, 1998. Vol. 85. No. 1. P. 115–129.

35.Захарова Т.В., Шестаков О.В. Вейвлет-анализ и его приложения. М.:

МАКС Пресс. 2009. 152 С.

36.Маркин А.В., Шестаков О.В. Отсев эктопических импульсов из ритмограммы с использованием робастных оценок//Информатика и ее применения. 2008. Т 2. № 2. С. 47-54.

37.Кудрявцев А.А., Шестаков О.В. Асимптотика оценки риска при вейглет-вейвлет разложении наблюдаемого сигнала//T-Comm — Телекоммуникации и Транспорт. 2011. №2 С. 54-57.

38.Bashforth F, Adams JC. An attempt to test the theory of capillary action.

Cambridge, 1892.

39.Blaisdell BE. The Physical Properties of Fluid Interfaces of Large Radius of Curvature, J. Math Phys 1940; 19:186.

40.Tawde NR, Parvatikar KG. Unstable pendant drops in relation to dropweight method of surface tension. Indian J Phys 1958;32:174.

41.Fordham S. On the Calculation of Surface Tension from Measurements of.

Pendant Drops. Proc R Soc Lond 1948;194A:1.

42.Mills OS. Tables for use in the measurement of interfacial tensions between liquids with small density differences. Br J Appl Phys 1953;4:24.

43. Paddy JF. In: Matijevic E, editor. Surface Tension. Part II. The Measurement of Surface Tension, vol.1. New York: Wiley; 1968.

44.Hartland S, Hartley RW. Axisymmetric fluid–liquid interfaces. Amsterdam:

Elsevier; 1976.

45.Maze C., Burnet G. A Non-linear Regression Method for Calculating the Surface Tension and Contact Angle from the Shape of a Sessile Drop // Surf.

Sci. 1969. V. 13. P. 451.

46.Huh C, Reed RL. J Colloid Interface Sci vol. 91, 1983, p. 472.

47.Boyce JF, Schrch S, Rotenberg Y, Neumann AW. The measurement of surface and interfacial tension by the axisymmetric drop technique. Colloids Surf. Vol. 9,, 1984, p. 307.

48.The determination of interfacial tension by video image processing of pendant fluid drops / Anastasiadis SH, Chen JK, Koberstein JT, et al. J Colloid Interface Sci, vol. 119, p. 55, 1986.

49.Hoggs RV. An introduction to robust estimation. NewYork: Academic Press;

1979.

50.Rotenberg Y, Boruvka L, Neumann AW. Determination of surface tension and contact angle from the shapes of axisymmetric fluid interfaces. J Colloid Interface Sci, vol. 93, 1983, p. 169.

51.Automation of axisymmetric drop shape analysis for measurements of interfacial tensions and contact angles / Cheng P, Li D, Boruvka L, Rotenberg Y, Neumann AW. Colloids Surf. vol. 43, 1990, p. 151.

52.Cheng P, Neumann AW. Computational evaluation of axisymmetric drop shape analysis-profile (ADSA-P). Colloids Surf vol. 62, 1992, p. 297.

53.del Ro OI, Neumann AW. Axisymmetric drop shape analysis:

Computational methods for the measurement of interfacial properties from the shape and dimensions of pendant and sessile drops. J Colloid Interface Sci. vol. 196, 1997, p. 136.

54.Dormand, J.R., Prince, P.J., A family of embedded Runge-Kutta formulae.

Journal of Computational and Applied Mathematics 6 (1): 19–26 (1980).

55.Canny, J., A Computational Approach To Edge Detection // IEEE Trans.

Pattern Analysis and Machine Intelligence, 8(6):679–698, 1986.