Коррекция дисторсии без «шахматной доски» / Хабр

Классическая калибровка камеры в компьютерном зрении
Альтернатива шахматной доске
Оптический метод коррекции дисторсии
Коррекция дисторсии по таблице из Zemax
— Пример: расчёт по таблице
— Дополнение: смещение сенсора
Коррекция по графику дисторсии
— Дополнение: неопределённость фокусного расстояния
Расчёт угловых полей
— Пример: расчёт угловых полей для масштабированного изображения
— Пример: расчёт предельных угловых полей
Программная реализация оптического метода
— Пример кода для исправления дисторсии
— Примеры исправленных изображений
— Режимы масштабирования
— Подготовка данных
— Расчёт коэффициента масштаба
— Построение карты ремаппинга

Привет, Хабр!

Меня зовут Андрей, я – специалист по оптическим системам, оптик и инженер-конструктор в одном лице.

Эта статья написана в соавторстве с коллегами из отдела нейросетей и компьютерного зрения:

Павлом @WIIC
Андреем @locoderr

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

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

Классическая калибровка камеры в компьютерном зрении

Автор раздела: Андрей @locoderr

При необходимости исправить дисторсию любой программист компьютерного зрения, вероятно, воспользуется стандартным решением – проведёт калибровку с помощью шахматной доски. Этот подход давно изучен и реализован в различных библиотеках (например, OpenCV) и его можно использовать прямо из коробки.

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

Если калибровка камеры для вас новая тема, рекомендую прочитать статьи от @Ilya12c:

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

Интересный факт

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

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

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

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

Альтернатива шахматной доске

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

Зачем нужен альтернативный метод, если есть упомянутое выше готовое решение из OpenCV? Всё просто – оптический метод точнее корректирует дисторсию, поскольку опирается напрямую на параметры объектива, а не на исследование изображения, создаваемое объективом.

Оптический метод коррекции дисторсии имеет ряд преимуществ перед калибровкой:

Более точное исправление дисторсии (в случае использования таблицы коэффициентов).
Возможность исправить дисторсию там, где калибровка не вытягивает (например, на краях fisheye-объектива).
Возможность «заочно», «удалённо» или заранее исправить дисторсию без работы с самим объективом. Например, исправить её в тот момент, когда объектив ещё только едет от продавца.
Легко играться с кадрированием и масштабированием (и даже со сменой сенсоров!) без потери исправления дисторсии.

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

Почему же при таком сильном наборе преимуществ этот способ практически не используется?

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

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

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

«Перспективная» дисторсия. Картинка из статьи про калибровку камеры

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

Оптический метод коррекции дисторсии

Информация о дисторсии объектива приводится в виде графика зависимости дисторсии от угла поля зрения. Такой график присутствует в спецификации любого серьёзного объектива.

График дисторсии в Zemax строится так:

Через объектив трассируется множество главных лучей с постепенным увеличением угла поля зрения ω: 0°, 0.50°, 1°, 1.50° и так далее.
Для каждого главного луча определяется высота его падения в плоскости изображения относительно оптической оси (Real Height).
Также для каждого главного луча рассчитывается идеальная высота падения, как если бы в объективе полностью отсутствовала дисторсия (Reference Height): Reference Height = tg ω⋅f′
Для каждого угла поля зрения из первого пункта вычисляется значение дисторсии по формуле:

$Distortion=\frac{RealHeight-Ref. Height}{Ref.Height}\cdot\text{100\%}$

График дисторсии получается интерполяцией этих точек.

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

Коррекция дисторсии по таблице из Zemax

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

Такую таблицу программист может получить двумя способами:

Если объектив рассчитывали в той же организации, то можно запросить таблицу у оптика-расчётчика. Такой вариант часто встречается на крупных оптических предприятиях, но для IT такой сценарий – редкость.
Если объектив покупной – можно попробовать запросить таблицу у продавца или изготовителя объектива.

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

видимый диапазон (Visible): 0,55 – 0,587 мкм
ближний инфракрасный диапазон (SWIR, Short-Wave InfraRed): 1,2 мкм
средний инфракрасный диапазон (MWIR, Middle-Wave InfraRed): 4 мкм
дальний инфракрасный диапазон (LWIR, Long-Wave InfraRed) 10 мкм

Для примера исправим дисторсию для объектива, таблица для которого известна:

Широкоугольный объектив F3 с фокусным расстоянием 35 мм и его график дисторсии

Перед началом работы необходимо проверить таблицу:

Тип дисторсии должен быть F-Tan(Theta). Это означает, что Ref. Height вычислены в классическом виде, через параксиальное фокусное расстояние объектива и тангенс полевого угла.
Убедиться, что используется таблица для нужной длины волны (середина оптического диапазона). Использование другой длины волны не внесёт заметной ошибки, но если есть возможность, лучше использовать центральную, особенно если в объективе заметен хроматизм увеличения.

Для исправления дисторсии из таблицы берутся столбцы Angle, Real Height и Ref. Height. Парадоксально, но столбец Distortion не используется – на него стоит смотреть только для справок. Второй и третий столбцы игнорируются – это данные для астигматизма (сдвиги астигматических фокусов от параксиальной плоскости изображения), и они нужны для другой задачи.

Если визуализировать столбцы Real Height и Ref. Height на модели идеальной линзы, то легко заметить, насколько сильна бочкообразная дисторсия в объективе:

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

Для этого на сенсор накладывается система координат с началом в точке пересечения сенсора с оптической осью (см. далее примечание про смещение сенсора):

Коррекция дисторсии через сдвиг координат. Real Height – красные, Ref. Height – серые — Коррекция дисторсии через сдвиг координат. *Real* *Height* – красные, *Ref. Height* – серые

Сдвиг координат производится по такому алгоритму:

В спецификации находится размер пикселя p – либо берётся напрямую, либо вычисляется делением длины и ширины оптической области сенсора на соответствующее разрешение в пикселях. Второй способ более точен, поскольку учитывает промежутки между микролинзами, которые некоторые производители не включают в размер пикселя.
Для каждого пикселя определяются его координаты X_Real и Y_Real(в пикселях).
X_Real и Y_Realпересчитываются в миллиметры или микрометры умножением координаты в пикселях на размер пикселя p.
По теореме Пифагора определяется Real Height этого пикселя.
По таблице из Zemax находится ближайшее значение Real Height.
Интерполяцией рассчитывается Ref. Height для Real Height пикселя.
Находится коэффициент дисторсии для этого пикселя: ∆d = Ref. Height / Real Height.
Находятся X_Ref и Y_Ref в миллиметрах по формулам: X_Ref = ∆d⋅ X_Real, Y_Ref = ∆d⋅ Y_Real
X_Ref и Y_Ref пересчитываются обратно в количество пикселей, делением координаты в миллиметрах на размер пикселя, с округлением до целого числа, поскольку координата в пикселях не может быть дробной.
Изображение с каждого пикселя сенсора (исходные координаты X_Real и Y_Real) переносится на пиксель, который ближе всего к координатам X_Ref и Y_Ref.

Интересный факт

Координаты в пикселях – всегда целые числа. Но размеры объектов в пикселях могут быть дробными. Получается это просто:

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

Примечание. Помните, что в алгоритме центр системы координат – в точке пересечения оптической оси и сенсора (чаще всего – в середине сенсора), а в самом сенсоре внутренняя кодировка пикселей начинается с одного из углов сенсора.

Пример расчёта представлен ниже, а пример реализации алгоритма в коде – в последней главе.

Пример: расчёт по таблице

Задача:

Рассчитать коррекцию дисторсии для таблицы выше, если с объективом используется Full HD сенсор с разрешением 1920×1080 и размером пикселя 10 мкм.

Решение:

Алгоритм идентичен для каждого пикселя. Для примера рассчитаем коррекцию для пикселя с координатами: –705, +395.

При отсутствии смещения сенсора в камере (см. ниже) центр системы координат по умолчанию расположен в центре кадра:

Иллюстрация к примеру. Для наглядности размеры не в масштабе

Размер пикселя: 10 мкм
Координаты пикселя: X_Real = –705 px; Y_Real = +395 px
Координаты пикселя в мм: X_Real = –7,05 мм ; Y_Real = +3,95 мм
Вычисляем Real Height: √[(–7,05)²+ 3,95²] = 8,08115 мм
Находим ближайшее табличное значение Real Height: 7,931 мм и видим, что Ref. Height для этого значения – 8,08037 мм
Интерполяцией находим Ref. Height для 8,08115 мм: Ref. Height = 8,23985 мм
Находим коэффициент дисторсии для этого Real Height: ∆d = 8,23985/8,08115 = 1,01964
Находим X_Ref и Y_Ref в миллиметрах: X_Ref = –7,05⋅1,01964 = –7,1885 мм, Y_Ref = 3,95⋅1,01964 = 4,0276 мм
Переводим найденные значения в пиксели: X_Ref = –7,1885/0,01 = –718,85 = –719 px, Y_Ref = 4,027 мм/0,01 = 402,76 = 403 px

Ближайшие к 8,08115 значения Real Height для интерполяции — Ближайшие к 8,08115 значения *Real Height* для интерполяции

Ответ:

Для исправления дисторсии изображение с пикселя –705, +395 необходимо перенести на пиксель –719, +403.

Дополнение: смещение сенсора

Обычно считается, что оптическая ось пересекает сенсор в его центре. Но бывают камеры, в которых сенсор конструктивно смещён относительно крепёжной резьбы и, соответственно, оптической оси объектива. Например, сенсоры семейства Python имеют смещение центра светочувствительной области относительно габаритов самого сенсора, и это же смещение без исправления присутствует в сконструированных на его базе камерах Mars Vision Datum Technology:

При столкновении с такой камерой необходимо разместить начало системы координат не в центре сенсора, а в центре крепёжной резьбы (на рисунке выше – center of package), потому что именно туда приходит оптическая ось объектива.

Коррекция по графику дисторсии

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

График дисторсии объектива из примера выше в высоком разрешении

Алгоритм:

Оцифровываем график любым онлайн-сервисом, либо вручную.
Дискретизируем график по полю и получаем столбцы Y Angle и Distortion. Дисторсия зависит от третьей степени угла поля зрения, потому шаг дискретизации для края поля необходимо брать максимально мелким, в то время как для центра поля, где дисторсия мала, шаг можно укрупнить.
Через тангенс для каждого Y Angle вычисляем соответствующие ему Ref. Height по формуле: Ref. Height = f′⋅tg(Y Angle), где f′ – фокусное расстояние объектива.
Вычисляем Real Height, выразив его из формулы в начале статьи. Distortion необходимо подставлять с учётом знака, поскольку знак дисторсии определяет её форму на изображении (подушкообразная или бочкообразная):
$RealHeight=\frac{Distortion\cdot Ref.Height}{{100\%}}+Ref.Height$
Далее по значениям Real Height и Ref. Height проводим коррекцию дисторсии, описанную выше.

Дополнение: неопределённость фокусного расстояния

Посмотрите на типичный объектив для технического зрения (даже не fisheye):

Типичный широкоугольный тепловизионный объектив

Внимание на приписку. Из-за сильной дисторсии в спецификации указаны два фокусных расстояния: реальное (real 11,5 mm) и вспомогательное условно-эквивалентное (10 mm), которое было бы, если бы указанное в спецификации поле зрения формировал объектив без дисторсии. Такая подмена фокусного расстояния – нестандартное для оптики решение, которое легко запутает неспециалиста при вычислении реальных угловых полей.

Настоящее фокусное расстояние объектива – 11,5 мм. Но даже на самом объективе может быть выгравировано более «удобное» 10 мм. Будьте бдительны – всегда используйте в расчётах настоящее фокусное расстояние!

Примечание. Полям зрения (FOV – field of view) в спецификации доверять нельзя – иногда их указывают неверно. Всегда стоит перепроверять угловые поля зрения по формулам арктангенсов. Проблема в том, что как раз при дисторсии эти формулы перестают работать.

Расчёт угловых полей

В ряде задач технического зрения (чаще всего – для измерения изображений) критически необходимо знать точные угловые поля объектива. И вот тут неопытный разработчик ступает на очень тонкий лёд.

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

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

сенсор камеры и область просмотра (например, окно отображения в GUI) имеют одинаковое разрешение (например 1280×1024 пикселей)
изображение с сенсора переносится в область просмотра «пиксель в пиксель»
полудлина и полуширина сенсора равны двум строкам из таблицы:

Величины Real Height, выбранные для габаритов сенсора — Величины *Real Height*, выбранные для габаритов сенсора

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

Коррекция дисторсии только сдвигом пикселей. Серая область выходит за границы области просмотра

Расчёт дисторсии для края кадра по таблице

$L=\sqrt{17,59^{2}+24,09^{2}}=29,83$ $L'=L\cdot\frac{41,7113}{29,803}=41,75$

Спойлер. Строки для длины и ширины сенсора специально выбраны так, чтобы диагональ изображения без дисторсии (41,75) была наиболее близка к максимальному Ref. Height (41.7113) из таблицы. Это позволило лучше всего визуализировать дисторсию выбранного объектива на иллюстрации.

Примечание. Из-за нелинейного характера дисторсии форма дуг будет иной. На рисунке формирование граничных линий по трём точкам сделано для наглядности. Реальные формы границ откорректированного изображения можно посмотреть в статье от @Dmitry_Po – они чем-то похожи на гиперболы:

Разумеется, такие потери краёв нерациональны. Чтобы свести теряемую площадь изображения к минимуму, используют два приёма:

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

Оба приёма могут применяться как независимо, так и совместно.

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

Разные способы совмещения ортоскопического изображения с сенсором

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

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

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

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

Сформулируем разницу:

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

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

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

Пример: расчёт угловых полей для масштабированного изображения

Задача:

Рассчитать горизонтальное, вертикальное и диагональное поля зрения для широкоугольного объектива из примера выше (f′ = 35 мм) в случае масштабирования изображения по высоте.

Решение:

Рассмотрим, как меняется изображение при масштабировании по высоте:

Прямоугольная рамка – светочувствительная область сенсора

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

вертикальное угловое поле =
$2⋅arctg(\frac{AO}{f'})=2⋅arctg(\frac{17.59}{35})=53,3656°≈53,37°$
горизонтальное угловое поле =
$2⋅arctg(\frac{BO}{f'})=2⋅arctg(\frac{24.09}{35})=69,0782°≈69,08°$
диагональное угловое поле =
$2⋅arctg(\frac{CO}{f'})=2⋅arctg(\frac{29,83}{35})=80,88096°≈80,88°$

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

Коэффициент масштабирования по высоте = 17,59/19,4 = 0,9067

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

коэффициент масштабирования по высоте = 17,59/19,4 = 0,9067
коэффициент масштабирования по ширине = 24,09/29,37 = 0,8478
коэффициент масштабирования по диагонали = 29,83/41,75 = 0,7145

Когда мы вписывали изображение в область просмотра, мы умножали отрезки AO, BO и CO на коэффициент масштаба. Чтобы вернуться к физическим размерам изображения (которое объектив строит в плоскости сенсора), необходимо поделить эти размеры на коэффициент масштаба. Итого:

вертикальное угловое поле с учётом масштабирования:
$2\cdot arctg\left(\frac{17,59}{35}\cdot\frac{1}{0,9067}\right)=59,998°\approx60°$
горизонтальное угловое поле с учётом масштабирования:
$2\cdot arctg\left(\frac{24,09}{35}\cdot\frac{1}{0,9067}\right)=74,405°\approx 74,4°$
диагональное угловое поле с учётом масштабирования:
$2\cdot arctg\left(\frac{29,83}{35}\cdot\frac{1}{0,9067}\right)=86,4562°\approx 86,5°$

Ответ:

Угловые поля зрения с учётом масштабирования изображения равны:

вертикальное угловое поле = 60°
горизонтальное угловое поле = 74,4°
диагональное угловое поле = 86,5°

Обратите внимание – для расчёта полей из таблицы Zemax берётся только одна пара цифр: 17,59 и 19,4, при помощи которых определяется коэффициент масштаба. Далее таблица дисторсии не используется, все остальные углы высчитываются только по коэффициенту масштаба.

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

Пример: расчёт предельных угловых полей

Если необходимо рассчитать предельные физические угловые поля, которые ещё можно увидеть сквозь пару «объектив + сенсор», то они рассчитываются только по таблице / размеру изображения. Коэффициент масштабирования и размеры области просмотра не учитываются, поскольку они никак не связаны с физическими параметрами оптической системы – фокусным расстоянием объектива и размерами светочувствительной области сенсора.

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

Предельные угловые поля для этого объектива и этого сенсора:

вертикальное угловое поле = 2⋅arctg(AO/f′) = 2⋅arctg(19,4/35) = 57,998° ≈ 58°
горизонтальное угловое поле = 2⋅arctg(BO/f′) = 2⋅arctg(29,37/35) = 80,0029° ≈ 80°
диагональное угловое поле = 2⋅arctg(CO/f′) = 2⋅arctg(41,75/35) = 100,0522° ≈ 100°

Те же самые значения мы видим в таблице:

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

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

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

Программная реализация оптического метода

Автор раздела: Павел@WIIC

Реализация алгоритма выполнена на Python с использованием библиотек NumPy, SciPy и OpenCV. Функция принимает изображение, таблицу дисторсии из Zemax и размер пикселя сенсора, а также режим масштабирования – способ вписывания подушкообразного ортоскопического изображения в область просмотра.

Пример кода для исправления дисторсии

Местоположение: github

Код

import numpy as np
import cv2
from scipy.interpolate import interp1d
import os
from enum import Enum


class ScalingMode(Enum):
    NONE = "none"      # Без масштабирования
    HEIGHT = "height"  # Масштабирование по высоте
    WIDTH = "width"    # Масштабирование по ширине
    FULL = "full"      # Полное масштабирование


def correct_distortion(image_path: str, 
                       table_data: list, 
                       sensor_pixel_size: float, 
                       scaling_mode = ScalingMode.NONE, 
                       output_path: str = None
                       ) -> np.ndarray[np.uint8]:
    
    img_array = cv2.imread(image_path, cv2.IMREAD_GRAYSCALE)
    height, width = img_array.shape

    cx = width / 2.0
    cy = height / 2.0

    # мкм -> мм
    pixel_size_mm = sensor_pixel_size / 1000.0

    real_heights_mm = np.array([row[1] for row in table_data])
    ref_heights_mm = np.array([row[2] for row in table_data])

    # мм -> пиксели
    real_px = real_heights_mm / pixel_size_mm
    ref_px = ref_heights_mm / pixel_size_mm

    # ref -> real: для построения карты ремаппинга
    interp_ref_to_real = interp1d(ref_px, real_px, kind="linear", fill_value="extrapolate")

    # real -> ref: для расчёта коэффициента масштаба
    interp_real_to_ref = interp1d(real_px, ref_px, kind="linear", fill_value="extrapolate")

    # Коэффициент масштаба в зависимости от режима
    if scaling_mode == ScalingMode.NONE:
        scale = 1.0
    elif scaling_mode == ScalingMode.HEIGHT:
        scale = cy / float(interp_real_to_ref(cy))
    elif scaling_mode == ScalingMode.WIDTH:
        scale = cx / float(interp_real_to_ref(cx))
    elif scaling_mode == ScalingMode.FULL:
        r_corner = np.sqrt(cx ** 2 + cy ** 2)
        scale = r_corner / float(interp_real_to_ref(r_corner))

    # Сетка координат выходного изображения
    y_out, x_out = np.mgrid[0:height, 0:width]
    x_rel = x_out - cx
    y_rel = y_out - cy

    r_out = np.sqrt(x_rel ** 2 + y_rel ** 2)
    theta = np.arctan2(y_rel, x_rel)

    # Радиус в системе ортоскопического изображения
    r_ref = r_out / scale

    # Радиус в исходном изображении
    r_real = interp_ref_to_real(r_ref).astype(np.float64)

    # Обратно в декартовы координаты
    x_src = cx + r_real * np.cos(theta)
    y_src = cy + r_real * np.sin(theta)

    map_x = x_src.astype(np.float32)
    map_y = y_src.astype(np.float32)

    # Пиксели вне границ сенсора определяем нулём
    corrected = cv2.remap(img_array, map_x, map_y, cv2.INTER_LINEAR, 
                          borderMode=cv2.BORDER_CONSTANT, borderValue=0)

    if output_path is None:
        base, ext = os.path.splitext(image_path)
        output_path = f"{base}_corrected_{scaling_mode.value}{ext}"

    cv2.imwrite(output_path, corrected)
    print(f"Сохранено: {output_path}")
    return corrected


def visualize_correction(original_path: str, results: list[tuple[ScalingMode, np.ndarray]]):
    import matplotlib.pyplot as plt
    
    mode_titles = {
        ScalingMode.NONE: "Без масштабирования",
        ScalingMode.HEIGHT: "По высоте",
        ScalingMode.WIDTH: "По ширине",
        ScalingMode.FULL: "Полное\nмасштабирование",
    }

    n = 1 + len(results)
    fig, axes = plt.subplots(1, n, figsize=(5 * n, 5))
    axes: list[plt.Axes] # для удобства, на самом деле это numpy.ndarray

    original = cv2.imread(original_path, cv2.IMREAD_GRAYSCALE)
    axes[0].imshow(original, cmap="gray")
    axes[0].set_title("Оригинал")
    axes[0].axis("off")

    for i, (mode, corrected) in enumerate(results):
        axes[i + 1].imshow(corrected, cmap="gray")
        axes[i + 1].set_title(mode_titles.get(mode, mode.value))
        axes[i + 1].axis("off")

    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    import data as dt

    SENSOR_PIXEL_SIZE = 12  # мкм
    IMAGE_PATH = "42.jpg"

    ext = os.path.splitext(IMAGE_PATH)[1].lower()

    if not os.path.exists(IMAGE_PATH):
        print("Изображение не найдено")
    else:
        results: list[tuple[ScalingMode, np.ndarray]] = []

        for mode in ScalingMode:
            corrected = correct_distortion(
                image_path=IMAGE_PATH,
                table_data=dt.table,
                sensor_pixel_size=SENSOR_PIXEL_SIZE,
                scaling_mode=mode)
            
            results.append((mode, corrected))

        visualize_correction(IMAGE_PATH, results)

Примеры исправленных изображений

Режимы масштабирования

Для задания режима введём перечисление ScalingMode:

class ScalingMode(Enum):
    NONE = "none"      # Без масштабирования
    HEIGHT = "height"  # Масштабирование по высоте
    WIDTH = "width"    # Масштабирование по ширине
    FULL = "full"      # Полное масштабирование

Подготовка данных

Первым шагом таблица переводится из миллиметров в пиксели делением на размер пикселя:

real_px = real_heights_mm / pixel_size_mm

ref_px = ref_heights_mm / pixel_size_mm

Затем строятся два интерполятора по значениям таблицы:

ref → real: для каждого радиуса в выходном (ортоскопическом) изображении интерполятор находит соответствующий радиус во входном (искажённом) изображении. Используется при построении карты ремаппинга.
real → ref: для заданного радиуса во входном изображении интерполятор находит соответствующий радиус Ref. Height в ортоскопическом изображении. Используется только для расчёта коэффициента масштаба.

Расчёт коэффициента масштаба

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

После вычисления scale алгоритм единообразен для всех режимов: радиус r_out каждого пикселя выходного изображения пересчитывается в радиус ортоскопического изображения r_ref = r_out / scale, а затем в радиус источника r_real в искажённом входном изображении.

Для режима без масштабирования (NONE) scale = 1: каждый пиксель выходного изображения напрямую соответствует ортоскопической позиции.

Для масштабирования по высоте (HEIGHT):

Берётся радиус до середины верхнего края сенсора: r_real = cy (половина высоты в пикселях)
По интерполятору real → ref определяется соответствующий Ref. Height: r_ref_top = interp_real_to_ref(cy)
Коэффициент масштаба: scale = cy / r_ref_top

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

Для масштабирования по ширине (WIDTH) расчёт аналогичен, но используется половина ширины: scale = cx / interp_real_to_ref(cx). Ортоскопическое изображение вписывается по горизонтали; верх и низ могут быть обрезаны.

Для полного масштабирования (FULL):

Вычисляется половина диагонали сенсора: r_corner = sqrt(cx² + cy²)
По интерполятору real → ref определяется Ref. Height для углового пикселя: r_ref_corner = interp_real_to_ref(r_corner)
Коэффициент масштаба: scale = r_corner / r_ref_corner

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

Примечание. Для объектива с бочкообразной дисторсией всегда ref > real, поэтому interp_real_to_ref(cy) > cy и scale < 1 – ортоскопическое изображение сжимается, чтобы уместиться в область просмотра.

Построение карты ремаппинга

После определения scale для каждого пикселя выходного изображения с координатами (x, y) выполняется следующее:

Вычисляются координаты относительно центра оптической оси: x_rel = x − cx, y_rel = y − cy
Определяется радиус пикселя в выходном изображении: r_out = √(x_rel² + y_rel²)
Вычисляется угол в полярной системе: θ = arctan2(y_rel, x_rel)
Радиус в ортоскопическом изображении: r_ref = r_out / scale
По интерполятору ref → real находится радиус источника: r_real = interp_ref_to_real(r_ref)
Источник пересчитывается обратно в декартовы координаты: x_src = cx + r_real×cos θ, y_src = cy + r_real×sin θ

Координаты x_src и y_src образуют карту ремаппинга, которая передаётся в cv2.remap. Функция для каждого пикселя выходного изображения берёт значение из рассчитанной позиции источника с билинейной интерполяцией. Пиксели, чьи координаты источника выходят за границы входного изображения, заполняются константой, в моём случае – нулём.

Спасибо за внимание.

Цикл статей по оптике для программистов технического зрения:

Коррекция дисторсии без «шахматной доски»