CC BY
41
ВОССТАНОВЛЕНИЕ ИНФОРМАЦИИ ОБ ОБЪЕКТЕ В ЦИФРОВОЙ ГОЛОГРАФИИ НА ОСНОВЕ УРАВНЕНИЯ ПЕРЕНОСА
С.А. Балтийский, И.П. Гуров
Записанное голографическим способом изображение объекта содержит информацию об амплитуде и фазе объектной волны. Методы цифровой голографии позволяют выделить и представить эту информацию в количественной форме. Реальные голографические системы имеют ограниченное разрешение, регистрируемые голограммы искажаются влиянием шума, что вызывает неоднозначность восстановления полной фазы объектной волны. Рассматриваемые в статье методы цифровой голографии основаны на использовании уравнения переноса интенсивности и представляют возможные подходы к решению проблемы развертывания фазы объектной волны.
Введение
Методы голографии, как известно, основываются на записи и последующем восстановлении полной информации об оптическом поле [1]. При записи голограммы регистрируется картина интерференции объектной и опорной волн в некоторой плоскости. Для восстановления исходного оптического поля голограмма освещается волной, аналогичной опорной.
Цифровая голография обеспечивает восстановление информации об объектной волне при цифровой обработке зарегистрированного распределения интенсивности в плоскости голограммы, представленного в дискретной форме [2]. При этом оптическое поле характеризуется матрицей комплексных чисел, каждое из которых содержит информацию об амплитуде и фазе поля в соответствующей точке. Голограмма представляет собой зарегистрированное распределение интенсивности
(1)
где 0(х,у) и R(x,y) являются комплексным представлением исследуемой объектной волны и опорной волны, соответственно. При пространственной дискретизации значения /(.г, у) определяются для дискретных координат, и матрица отсчетов содержит вещественные значения, точность которых ограничивается разрядностью цифрового представления в регистрирующем устройстве (фотокамере или сканере). Матрица отсчетов может быть представлена в виде матрицы целых чисел без потери информации. Это позволяет хранить и передавать голографические изображения в виде растровых файлов-изображений стандартных форматов.
Восстановление изображения в цифровой голографии осуществляется в форме количественного расчета произведения значений интенсивности и комплексной амплитуды опорной волны, а именно
(2)
После подстановки правой части (1) в (2) результат восстановления поля в плоскости голограммы можно представить следующим образом:
R(x,y)l(xiy)=\R(x,yfo{x,y)+R2(x,yp(x,y)+
+ R(x, уЩх, у)|2 + R{x, ур(х, yf, (3)
где первое слагаемое в правой части уравнения является точной копией исходного поля 0\х, v) = \0{х, v)j ехр(/ф0 (л\ у)), остальные слагаемые являются неинформативными и
устраняются специальными методами цифровой обработки (см., например, [2] - [4]). В результате оказывается возможным выделить информацию об амплитуде и фазе объектной волны в количественной форме.
Восстановление информации о фазе оптического поля
Полученная при восстановлении оптического поля в плоскости голограммы матрица комплексных чисел 0{х,у) содержит информацию об амплитуде и фазе исходного поля в каждой точке. Вещественная матрица значений интенсивности, полученная возведением в квадрат амплитуды в каждой точке, образует «фотографическое» изображение объекта
А(х,у) = \0(х,у}2. (4)
Матрица значений фазы вычисляется как
ц>{х,у)^аг^(0(х,у)). (5)
При определенных условиях вычисленные значения фазы (5) позволяют восстановить форму поверхности объекта, поскольку значения полной фазы объектной волны однозначно связаны с длиной оптическог о пути для каждой точки поверхности.
Построенное в (4) амплитудное изображение объекта является четким только в том случае, если записываемое оптическое поле было предварительно сфокусировано в плоскости записи голограммы. Значения фазы (5) являются пригодными для реконструкции формы поверхности, если они получены из данных при точной фокусировке.
В отличие от обычного фотографического изображения, отсутствие фокусировки которого, как правило, не может быть скомпенсировано после записи математическими методами, содержащаяся в голограмме информация является достаточной для виртуальной фокусировки цифровыми методами на этапе восстановления изображения из голограммы. Для этого, в соответствии со скалярной теорией дифракции, требуется рассчитать распространение оптического поля из плоскости голограммы до параллельной ей плоскости изображения, находящейся на фокусном расстоянии от плоскости голограммы вблизи «поверхности» мнимого изображения объекта (см. рис. 1),
Рис. 1. Относительное расположение плоскости фокусировки и восстановленного мнимого изображения объекта (1), плоскости голограммы (2), исходного объекта во время записи(3)
Фокусное расстояние может быть вычислено исходя из параметров установки записи и расстояния до записываемого объекта, либо подобрано при наблюдении изображения объекта, полученного из (4).
Расчёт распространения оптического поля обычно основывается на применении приближения Френеля для интеграла дифракции Релея-Зоммерфельда. Восстановлен-
мое дифрагированное поле О в плоскости изображения на расстоянии гт% от плоскости голограммы при учете параксиального приближения может быть записано в следующем виде [5], [6]:
О7(х>у)- —— expí/—г |х
^ У) 1кг Ч Ь )
х £ £ *(и>у)ехр(^((* - и)2 + (у - V)2 . (6)
Интеграл свертки в (6) обладает слишком большой вычислительной сложностью для непосредственного применения в расчётах распространения оптического поля. В целях оптимизации вычислений интеграл свертки может быть приведен к виду преобразования Френеля от интенсивности голограммы, что позволяет численно рассчитать его с помощью двумерного преобразования Фурье (см., например, [7]), в том числе с использованием алгоритмов быстрого преобразования Фурье (БПФ).
Фокусировка достигается при достаточно малом расстоянии от всех точек поверхности объекта до плоскости фокусировки. В случае если этот критерий не может быть удовлетворён для плоскости фокусировки из-за большой разницы высот в различных точках объекта, требуется использовать набор независимых плоскостей фокусировки для разных участков объекта. При этом матрица значений фазы может быть преобразована в матрицу относительных разностей хода с1{х,у) исходя из очевидного соотношения разности фаз и разности хода волн
Дф = Ш . (7)
При условии, что плоскость восстановления находится достаточно близко к поверхности объекта, можно считать, что расстояние отсчитывается перпендикулярно этой плоскости. Значение оптической разности хода й при такой интерпретации соответствует высоте точки объекта, находящейся напротив рассматриваемой точки плоскости, относительно некоторой общей начальной точки объекта, высота в которой принимается равной нулю.
Для голо графических установок, осуществляющих запись в отражённых лучах, все расстояния должны быть поделены на два, так как оптическая разность хода составлена из длины пути к поверхности объекта и обратного пути от поверхности к плоскости восстановления. Таким образом, полученная матрица ¿(х,у) является по сути картой высот объекта, которая определяет трёхмерную форму его поверхности.
Рис. 2. Соответствие значений фазы и разности хода в пикселах плоскости восстановления. 1 — полутоновая картина развёрнутой фазы в пикселах плоскости восстановления, 2 — плоскость восстановления оптического поля, 3 — оптическая разность хода, соответствующая значению фазы, 4 — поверхность объекта, аппроксимированная площадками с заданной высотой (пунктирная линия), 5 — истинная поверхность объекта (сплошная линия)
На рис, 2 иллюстрируется соответствие оптической разности хода высоте соответствующих точек на поверхности объекта при выборе самой верхней точки рельефа в качестве начальной точки.
Поскольку обрабатываемое поле представлено в дискретной форме, каждой точке матрицы d{x,y) соответствует прямоугольный «пиксел» размерами Дц х Дт]. На поверхности объекта ему соответствует некоторая область, более сложная из-за отклонений формы объекта, получаемая проекцией прямоугольника на поверхность объекта.
Форма и положение поверхности в пространстве могут быть аппроксимированы прямоугольной площадкой с теми же размерами Д£ х Дг|, параллельной плоскости изображения и отстоящей от неё на расстояние, вычисленное из значения фазы. Набор таких площадок, построенных для всех точек плоскости, аппроксимирует форму объекта в пределах этой плоскости (см. пример на рис. 3).
Рис. 3. Трёхмерное изображение формы поверхности объекта по карте высот: карта высот, представленная полутоновыми уровнями яркости (а), аппроксимация трёхмерной формы поверхности объекта набором прямоугольных площадок (б), сглаженная поверхность(в)
Основными препятствиями для получения точных данных о форме поверхности объекта являются аберрации объектива и необходимость развёртывания фазы.
Влияние аберраций проявляется в форме зависимости оптической разности хода, привносимой линзами объектива записывающей установки в фазу для различных точек, от расположения этих точек, так что при некоторой расфокусировке восстановленная поверхность оказывается вогнутой. Наиболее простым, но при этом надёжным способом компенсации аберраций является восстановление формы поверхности из голограммы плоского образца, которая после внесения аберраций непосредственно описывает величину искажения для каждой точки матрицы высот [8, 9].
Необходимость развёртывания фазы следует из того, что поточечный способ построения фазы как аргумента комплексного значения поля приводит к свёртыванию значений фазы в диапазон (- я, я], как это иллюстрируется на рис. 4.
Восстановление развёрнутой фазы требует использования взаимосвязи между различными точками оцифрованного комплексного поля. Данная задача не имеет точного аналитического решения или хотя бы точных и устойчивых критериев для проверки правильности результата. Приближённое численное решение может быть найдено с помощью одного из известных эвристических алгоритмов, реализующих разные методики и основанные на различных критериях. При этом невозможно выделить один универсальный алгоритм, который подходил бы для обработки любого набора входных данных.
Рис, 4. Одномерное представление изменяющихся значений фазы; линейные изменения фазы (например, для падающей под фиксированным углом плоской волны) (а) и квадратично изменяющиеся значения фазы (б); сплошной линией показана свёрнутая фаза, пунктирной — значения исходной (или развёрнутой) фазы; на оси ординат отложены значения фазы, нормированные на величину п
Проблема развёртывания фазы обусловлена ограниченным разрешением регистрирующей системы и, следовательно, оцифрованных значений свёрнутой фазы, влиянием шумов, а также возможным наличием сингулярных точек. Интенсивность оптического поля в точках сингулярности равна нулю, и значение фазы не определено, являясь значением аргумента комплексного числа с нулевым модулем. Сложность представляют не сами такие точки, а их окрестности, так как в них неизбежно нарушается непрерывность восстановленной фазы при движении вокруг точки, поскольку процедура развёртывания фазы заключается в интегрировании приращения фазы вдоль координат. Значение фазы присваивается новым точкам на основании уже известных значений в соседних точках с использованием некоторых дополнительных критериев.
Применение уравнения переноса интенсивности к задаче развёртывания фазы
Алгоритмы развертывания фазы можно подразделить на две основные группы. В алгоритмах первой группы используются значения комплексного поля только в той же плоскости, в которой происходит восстановление; в частном случае принимается в расчет только фаза комплексного числа, т. е. входными данными для алгоритма является вещественная матрица свёрнутой фазы. Такие алгоритмы являются наиболее распространёнными, поскольку пригодны для решения задач интерферометрии, где доступная информация ограничивается интерференционной картиной, по смыслу аналогичной картине свёрнутой фазы. Алгоритмы, применяемые для обработки распределений фаз интерференционных картин достаточно высокого качества, позволяют использовать матрицу вещественных или комплексных чисел в готовом виде и абстрагироваться от способа получения этой информации.
Алгоритмы данной группы могут использовать глобальные либо локальные критерии оптимальности, что означает наличие, либо отсутствие, критериев качества результата развёртывания в целом; использовать или не использовать информацию об интенсивности в точке в качестве локального критерия или для выявления точек сингулярности (см., например, [4]) и т. д.
Алгоритмы второй группы используют информацию за пределами основной плоскости восстановления для уточнения значений развёрнутой фазы. Такие алгоритмы представляются предпочтительными для цифровой голографии, так как позволяют более полно использовать имеющуюся в распоряжении информацию о восстановленном
полном оптическом поле. Именно к этой категории относятся алгоритмы развёртывания фазы на основании уравнения переноса интенсивности.
Исходными данными для подобных алгоритмов является информация об интенсивности электромагнитного поля /(*,>>, 2), заданная вместе со своей первой производной 1(х, у, г) в плоскости, в которой осуществляется развёртывание фазы. В данном случае такой плоскостью является плоскость изображения, г = . Первая производная интенсивности может быть получена с помощью разностной аппроксимации [11]:
^ 1
—1\х,у,2тъ)*—[1 (х, у, + Дг)- /{*» у, 2;П1е)], (8)
где Д2 мало, Для этого требуются значения интенсивности оптического поля уже в двух близкорасположенных плоскостях, одна из которых совпадает с плоскостью изображения, а другая отстоит на малое расстояние Дг в какую-либо сторону.
Может быть использована также аппроксимация первой производной более высокого порядка на основе вычисления центральной разности в форме
£ 4х' У• 2<те)« ~ И*' У> г1тВ + Лг/2)-/(х, у, 2,т0 - Дг / 2)]. (9)
Применение формулы (9) предполагает работу с тремя плоскостями, одна из которых является плоскостью развёртывания фазы, а две другие, расположенные на расстоянии Дг/2 по разные стороны от неё, предоставляют значения интенсивности оптического поля. При дальнейшем рассмотрении будем опускать индекс 2, считая, что
В идеальном случае исходных данных в виде значений интенсивности в двух плоскостях достаточно для восстановления фазы (с точностью до граничных значений), и данные могут быть получены без использования интерферометрических установок прямой записью оптического поля для двух близких фокусных расстояний. Наличие в системе фокусирующих линз (или их аналога) при этом является обязательным, в отличие от голо графической записи, при которой они нужны лишь для изменения углового размера объекта для обеспечения соответствия условиям приближений, принимаемых при обработке голограмм. Поэтому в данном случае влияние аберраций оптической системы на значения восстановленной фазы неизбежно. Требование когерентности излучения источника остаётся в силе, однако интерференция с опорной волной при записи не требуется. Это свойство автоматически устраняет неинформативные составляющие поля в уравнении (3), которые возникают при традиционном восстановлении изображения из голограммы.
В случае обработки данных методами цифровой голографии значения интенсивности во второй плоскости могут быть получены численным расчётом распространения оптического поля на малое расстояние Д2 в соответствии со скалярной теорией дифракции по аналогии с процедурой получения сфокусированного изображения.
К трудностям применения данного подхода относится высокий уровень шума при формировании картины интенсивности и цифровой фокусировке. Отчасти появление шума связано с влиянием неинформативных слагаемых поля, например, сильно расфокусированного действительного изображения, наложенного поверх мнимого изображения. Запись нескольких голограмм с шаговым изменением разности хода между опорной и объектной волной позволяет подавить эти слагаемые, но значительно усложняет процесс записи и ужесточает требования к стабильности установки [4].
Значения интенсивности, используемые в методе уравнения переноса, могут быть выведены из параксиального дифференциального уравнения (см. [12] - [14])
г
. д V2 , г — + — + к дг 2 к
\
0{х,у) = О,
(10)
которое разделяет производные по пространственным осям, и оператор в левой части
2 г)2 д2
подразумевается двумерным, V = + .
Уравнение переноса интенсивности может быть записано в виде
VШ—к^
дг
(И) (12)
а уравнение переноса фазы - в форме
¿-/V2/ ~-|-(У/)2 -/2(Уср)2 + 1к212 =к21г~. I ^ 02
Уравнение (11) содержит достаточную информацию для восстановления фазы при наличии всех необходимых краевых условий.
Рассмотрим вначале одномерный случай. При этом уравнение переноса интенсивности (11) принимает вид
А
дх
скр
дх
д!
/гг
дг
(13)
Аналитическое решение уравнения (13) получается двойным интегрированием по .г в виде
дх
д!{м>)
■«о дг
<Ы.
(14)
Из этого соотношения видно, что для численного восстановления развёрнутой фазы достаточно задать значения фазы в двух первых точках линии и далее рекуррентно находить значения в каждой следующей точке.
Другим возможным подходом является задание в некоторой точке Л'0 значения
"дф(х)~
фазы ф(х0) и её производной
дх
. В принципе, требуется задать только произ-
водную, так как значение развёрнутой фазы в любом случае находится лишь с точностью до константы. Для правильного задания производной начальная точка выбирается гак, чтобы первое слагаемое под внешним интегралом уравнения (14), содержащее производную в качестве множителя, было мало по сравнению со вторым слагаемым.
Для получения численного аналога уравнения (14) раскроем скобки в уравнении (13), при этом получим
(15)
дх дх дх2 дг
Далее подставим разностное выражение для каждой из частных производных первого или второго порядка. Выражение для выбирается из (8) или (9) в зависимости от
используемой конфигурации плоскостей, а значения для первых производных и ^
строятся по аналогии с (8). Вторая производная использует разностную аппроксимацию
ф(х,у) к [ф(.х + Дг, у) + ф(х - Дх, у)- 2ф(х,у)].
Дх
(16)
Результат подстановки разностных аппроксимаций в уравнение (15) является довольно громоздким в аналитическом описании- По существу, реализуется разностная схема в виде четырёхточечного шаблона, показанного на рис. 5.
Значение развёрнутой фазы в очередной точке на оси х строится на основании уже вычисленного значения развёрнутой фазы в двух соседних точках и значений интенсивности в трёх точках для различных значений координат х и г. Начальными условиями являются значения в первых двух точках, либо в первой и в последней точке.
Рис, 6. 7-точечный шаблон разностной схемы, соответствующей уравнению (18)
Рис. 5. Четырёхточечный шаблон разностной схемы, соответствующей уравнению (15)
Практический интерес представляет двумерный случай, при котором возможно построить трёхмерную поверхность объекта, используя матрицу значений развёрнутой фазы.
Раскроем оператор в левой части уравнения переноса интенсивности (11) в форме
(я А ( Я N
у(Мр)= — Уф + / —.Уф
I дх
дх
<н_
дх
дх
ду.
д ц> а2ф
дх2 дхду
'дГ
ду;
61 +■-
ау
/ (еУ
+ ду2
дхду
Л ду]
=(гх + /; & +(/; + 4 ++ 1<р%, + 2/Ф;. (17)
Приравняв результат правой части уравнения (11), получим следующее дифференциальное уравнение в частных производных второго порядка относительно ф :
й+к + (/;++(18) Разностные аппроксимации для производных Гх, 1'у,1'х, ф'т и ф^ основы-
ваются на формулах (8), (9) и (16). Для смешанной производной ф" запишем:
ф(*, у) - Ь(х + Ах, у) + ф(дг, у + Ау)- ф(*, у) - ф(* + ДХ, у + Ду)]. (19) охоу ДеДу
Аналитический результат подстановки разностных аппроксимаций в уравнение (18) описывает трёхмерную разностную схему с 7-точечным шаблоном, изображённым на рис. 6, которая позволяет получить значение развёрнутой фазы следующей точки ф(* + Ах,у + Ду) с использованием значений в пяти точках с меньшими координатами в той же плоскости, значение развёрнутой фазы в которых уже известно. При этом возможно получить значения ф(х, у) во всей плоскости изображения при условии, что начальные (граничные) условия заданы в достаточном объёме.
Самым удобным для вычисления является вариант граничных условий, при котором значения развёрнутой фазы заданы в два ряда вдоль края плоскости с минимальным значением координаты г ив два столбца вдоль смежной стороны с минимальными значениями у (см. рис. 7, а). При этом значения всех точек могут быть вычислены последовательно и в явном виде, начиная от угла с минимальными координатами. Аналогичный результат получится при задании значений одного ряда и одного столбца вместе с первыми производными в тех же точках. При отсутствии точек сингулярности внутри области и достаточном разрешении полученное решение является однозначным
[15].
Рис. 7, Начальные условия для решения уравнения (18) в виде двойного ряда точек вдоль двух смежных сторон (а) и в виде одного ряда точек вдоль границы области (б)
Более точный для расчёта вариант задания начальных условий заключается в определении значений разЕёрнутой фазы на всём периметре рассматриваемого участка плоскости (рис. 7, б). Значения на внутреннем участке плоскости могут быть рассчитаны, например, с помощью релаксационного алгоритма для уравнений второго порядка. Такой способ предпочтительнее в случае, если внутри области присутствуют большие скачки фазы или сильный шум, который может проявляться на особенно высоких точках рельефа объекта, «выпадающих» из фокусировки.
Точность разностной аппроксимации может оказаться недостаточной, тогда наблюдается «срыв» значений на всём ещё не обработанном прямоугольнике плоскости. В любом случае, если области «плохих» значений фазы образуют замкнутый контур, фаза внутри такого контура может быть восстановлена корректно в пределах этого контура, но независимо от фазы вне его, и «состыковать» такие участки поверхности возможно в интерактивном режиме при участии оператора, использующего априорные данные об объекте.
Оба варианта соответствуют задаче, при которой развёрнутая фаза уже известна на прямоугольном контуре, ограничивающем сложный объект, например, если объёмная фигура расположена на относительно плоской подложке, для которой восстановление практически неизменной фазы не представляет труда. Развёртывание фазы для произвольного объекта является значительно более сложной задачей.
Из-за требований к начальным условиям данный алгоритм не может применяться самостоятельно для восстановления развёрнутой фазы на некоторой матрице и поэтому не подходит для обработки пары записанных неголографическим способом изображений [12].
В случае записи и восстановления данных голографическим методом может использоваться дополнительная информация, в том числе на участках непрерывности неразвернутой фазы, а также результаты применения «обычных» алгоритмов развёртывания фазы в плоскости. Обеспечение совместной работы алгоритмов первого вида (на плоскости) и второго вида (на основании уравнения переноса) с целью более надёжного восстановления трёхмерной формы поверхности объекта из развёрнутой фазы представляет предмет дополнительных исследований.
Помимо рассмотренного выше алгоритма существуют другие методы восстановления развёрнутой фазы на основе уравнения переноса интенсивности. Следует отметить метод [14], в котором осуществляется независимое вычисление величины и направления градиента фазы в каждой точке рассматриваемой области, и затем выполняется восстановление развёрнутой фазы во всей плоскости на основе вычисленного градиента и произвольно выбранного начального значения. Метод избавляет от необходимости вычислять большой объём начальных условий, но в то же время накладывает существенные ограничения на характеристики обрабатываемого оптического поля.
Заключение
Восстановление значений развернутой фазы объектной волны, количественно характеризующих исследуемый объект с интерферометрической точностью, является сложной проблемой, прежде всего, ввиду того, что многие объекты, исследуемые методами голографии, имеют негладкие поверхности с локальными отклонениями рельефа, превышающими значение длины волны оптического излучения. При ограниченной разрешающей способности регистрирующей камеры, когда разность фаз в соседних пикселах превышает 2тг, невозможно обеспечить однозначность развертывания фазы при использовании излучения с одной длиной волны. Вместе с тем, значительный объем информации, содержащейся в оптических полях, позволяет найти подходы к решению рассматриваемой проблемы. Один из возможных подходов состоит в использовании метода уравнений переноса, рассмотренного в настоящей статье.
Литература
1. Оптическая голография /Под ред. Г. Колфилда. В 2-х т. М.: Мир, 1982.
2. Балтийский С.А., Гуров И.П., Де Никола С., Коппола Д., Ферраро П. Современные методы цифровой голографии. В кн.: Проблемы когерентной и нелинейной оптики /Под ред. И.П. Гурова и С.А. Козлова. СПб: СПбГУ ИТМО, 2004. С. 91-117.
3. Kreis Th., Jiiptner W, Suppression of the dc term in digital holography //Opt. Eng. 1997 V 36. P. 2357-2360.
4. Yamaguchi I., Inomoto O., Kato J. Surface shape measurement by phase-shifting digital holography // Proc. Fringe'2001, 4th Int. Workshop on Automatic Processing of Fringe Patterns. Bremen, 2001. P. 365-372.
5. Schnars U., Juptner W. Direct recording of holograms by a CCD target and numerical reconstruction //Appl. Opt. 1994. V. 33. P. 179-181.
6. Goodman J.W. Introduction to Fourier Optics. McGraw-Hill: New York, 1996.
7. Kreis Th., Juptner W. Principles of digital holography //Proc. Fringe* 1997, Int. Workshop on Automatic Processing of Fringe Patterns. Akademie-Verlag. 1997. P. 353-363.
8. De Nicola S., Ferraro P., FinizioA,, Pierattini G. Wave front reconstruction of Fresnel off-ax is holograms with compensation of aberrations by means of phase-shifting digital holography //Opt. Las. Eng. 2002. V.37. P. 331-340.
9. Baitiysky S., Gurov L, De Nicola S,, Ferraro P., Finizio A., Coppola G. Characterization of microelectromechanical systems by digital holography method //Imag. Sei. Journ. 2006. V. 54. N. 2. P.103-110.
10. Jüngling R. Phase unwrapping speckle wave fronts /German Aerospace Center DLR-IB. Report 554-00/14. 2000.
11. Самарский A.A. Введение в теорию разностных схем. M.: Наука, 1971.
12. Teague M. Deterministic phase retrieval: a Green's function solution // J. Opt. Soc. Am. 1983. V. 73. N. 11. P. 1434-1441.
13. Teague M. Image formation in terms of the transport equation // J. Opt. Soc. Am. A. 1985. V. 2. N. 11. P. 2019-2026.
14. Kolenovic E. Correlation between intensity and phase in monochromatic light // J. Opt. Soc. Am. A. 2005. V. 22. N. 5. P. 899-906.
15. Gureyev T., Roberts A., Nugent K. Partially coherent fields, the transport-of-intensity equation, and phase uniqueness//J. Opt. Soc. Am. A. 1995, V. 12. N. 9. P. 1942-1946,
