Алгоритм перепроецирования растровых изображений средствами программного комплекса «Топографический калькулятор» Текст научной статьи по специальности «Электротехника, электронная техника, информационные технологии»

УДК 911.2:551.48 Е. А. Паниди

АЛГОРИТМ ПЕРЕПРОЕЦИРОВАНИя РАСТРОВЬІХ ИЗОБРАЖЕНИЙ СРЕДСТВАМИ ПРОГРАММНОГО КОМПЛЕКСА «ТОПОГРАФИЧЕСКИЙ КАЛЬКУЛЯТОР»

В данной статье рассмотрена методика перепроецирования растровых изображений реализованная в программном комплексе «Топографический калькулятор» («Topographic Calculator») [1, 2].

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

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

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

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

Формат привязки World File соответствует второму типу. Это открытый формат и он поддерживается большинством современных ГИС. По этому в «Топографическом калькуляторе» для задания привязки растровых изображений был выбран именно этот формат.

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

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

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

© Е. А. Паниди, 2008

А расширение должно состоять из трех букв, первая из которых соответствует первой букве расширения файла изображения, вторая — последней букве расширения файла изображения, а третья — буква «w» (World). [3, 4]

Пример содержания файла привязки:

8,000000 разрешение растра по оси X

0,000000 параметры поворота матрицы, которые обычно равны нулю

0,000000

-8,000000 разрешение растра по оси У

6329621,756784 координата X центра левого верхнего пикселя

6660578,042448 координата У центра левого верхнего пикселя

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

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

В ГИС, как правило, системы координат разделяются на два типа «географические» и «проекционные» (в данном случае применена терминология ArcGIS). По сути своей первые — это обычные геодезические системы координат на поверхностях различных эллипсоидов и сфер, координаты в которых измеряются в градусах долготы и широты, вторые — плоские прямоугольные системы координат различных картографических проекций, координаты в которых обычно измеряются в метрах. Однако, отображать данные приходится на плоскости монитора компьютера, не зависимо от того, в какой системе координат они хранятся и в каких единицах измеряются их координаты. В связи с этим, геодезические системы координат в ГИС, вероятно, можно понимать как особый класс проекций, в которых координаты измеряются в градусах. Причем при отображении на мониторе, такая система координат выглядит как плоская прямоугольная: по оси абсцисс измеряются долготы, по оси ординат — широты. Земная поверхность в такой проекции выглядит в целом аналогично нормальным цилиндрическим проекциям.

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

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

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

Процесс перепроецирования растрового изображения, в соответствии с предлагаемым алгоритмом, можно разделить на три основных этапа:

• первый — вычисление размеров матрицы выходного изображения,

• второй — определение разрешения результирующего изображения,

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

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

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

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

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

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

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

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

Длина дуги параллели может быть определена по следующей формуле [5]:

S = r • l,

где l — абсолютное значение разности долгот начальной и конечной точек дуги,

а

r = N • cos ф — радиус параллели, ф — широта параллели, N =--------------------------------радиус

,--------- 1 - e sinfy "

„ [ TbY „

кривизны первого вертикала в данной точке, e = AJ1 - I — I — первый эксцентриситет

Земного эллипсоида, а и b — соответственно большая и малая полуоси Земного эллипсоида.

Длина дуги меридиана определяется следующим образом [6]:

S = a • (1 - e)2.

ч Л 3 2 45 4 А 1 , . „ . ^ (3 2 15 4

(ф2 - ф,) • I 1 + — e2 + — е4 + ... I-• (sin 2ф2 - sin 2ф.) • I — e2 + — е4 + ...

I 4 64 | 2 14 16

+ — • (sin 4ф2 - sin 4ф.) • I-15- e4 + ...

4 2 1 I 64

где a — большая полуось Земного эллипсоида, e — первый эксцентриситет Земного эллипсоида, ф1 и ф2 — широты начальной и конечной точек дуги

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

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

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

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

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

• сохранить значение цвета данного пикселя,

• от пиксельных координат центра данного пикселя перейти к проекционным координатам этой точки в системе координат исходной проекции,

• от проекционных координат точки в системе координат исходной проекции перейти к проекционным координатам в системе координат требуемой проекции,

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

• присвоить сохраненное значение цвета пикселю, на территорию которого попадает данная точка.

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

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

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

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

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

• считать значение цвета пикселя, на территорию которого попадает данная точка,

• присвоить считанное значение цвета текущему пикселю результирующего изображения.

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

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

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

Summary

Panidi E. A. Algorithm of the raster imagery Reprojecting by the instruments of program complex «Topographical calculator»

There is technique of the raster imagery reprojecting realized by author in program complex «Topographical calculator» is considered in the article. The reference format of images in any coordinate systems is disassembled. The general technique of transformation of the raster image from one projection to another is covered and necessary attending calculations are quoted. Features of the imagery reprojecting into geographical coordinate systems are considered.

Литература

1. Программный комплекс «Топографический калькулятор» («Topographic Calculator»). Автор — Паниди Е. А. 2. Афанасьев В. А., Павлова О. А., Паниди Е. А, Щербаков В. М. Автоматизированная система генерации номенклатуры российских топографических карт и пересчёта координат между различными системами // Сб. ст. «Теория и практика эколого-географических исследований». СПб. 2005. 3. Материалы сайта компании DATA + http://www.dataplus.ru. 4. Материалы сайта http://www.gis-lab.info. 5. Бугаевский Л.М. Математическая картография. М. 1998. 6. Закатов П. С. Курс высшей геодезии. М. 1976.