суббота, 17 декабря 2011 г.

Открытый софт для научных расчетов

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

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

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

Выбор софта для такой типичной задачки вполне себе наличествует:
На практике же оказалось не так радужно, как хотелось бы. Из того, что есть в дебиане, работает "из коробки" только пакет octave. Зато для него можно установить множество расширений - например, необходимый мне пакет octave-image. Если же расширения нет в дебиан, можно установить из репозитория octave. Например, если мы хотим поставить последнюю версию разширения для работы с изображениями, качаем архив image-1.0.15.tar.gz: и далее в текущей директории выполняем:

aptitude install octave3.2-headers 
sudo octave
pkg install image-1.0.15.tar.gz
exit

Как видим, впечатление octave производит самое приятное. К сожалению, не умеет задействовать несколько ядер процессора, это досадно.

Далее я занялся тем, для чего, собственно, и подбирал софт. А именно, построением 3D модели земных недр по космоснимкам высокого разрешения. Несколько ссылок из предметной области:


На 32-бит хосте обработка космоснимка порядка 60 Мб размером требует около 1 Гб ОЗУ. Если же выполнять передискретизацию (скажем, уменьшить матрицу значений в 3 или 5 раз), то и намного меньше. Для Landsat 7 панхроматический снимок будет имеет вдвое большее разрешение и примерно вчетверо больший размер, так что 60 Мб снимок видимого и ИК диапазонов  соответствует примерно 250 Мб панхроматическому. Значит,  8 Гб ОЗУ на 64 бит хосте должно хватить и на обработку панхромата без передискретизации. Интересно, удастся ли это сделать на 32 бит хосте. Скоро проверю и это, как отлажу расчет на ближнем ИК диапазоне. Update. На 32-бит хосте панхроматический снимок прочитать не удалось - octave сообщает, что невозможно выделить требуемый объем памяти.

Статьи и документация по Octave плюс еще некоторые полезные ссылки:

6 комментариев:

Анонимный комментирует...

А почему бы не использовать специализированный ГИС-софт? Тот же GRASS, к примеру? И преобразования Фурье есть, и 3D модели строит, и скриптуется на ура.

Алексей Печников комментирует...

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

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

На практике же куча скриптов для ArcGIS легко заменяются несколькими строчками в Octave, наверняка и для GRASS будет точно так же.

Анонимный комментирует...

Понял все аргументы, кроме одного -- почему это не известно, как реализован тот или иной алгоритм в GRASS?

Алексей Печников комментирует...

Да вот как пример описание фильтрации в GRASS: r.mfilter - Performs raster map matrix filter.

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

Далее, по примечанием, функция работает лишь для целочисленной исходной матрицы -для вычислений с плавающей точкой нужно использовать функцию r.mfilter.fp., для которой даже и не удалось найти, работает она с float или double числами (что может быть критично, знаю по опыту). Матрица фильтрации, судя по представленным примерам (и наличию делителя), может содержать только нули и единицы, вдобавок можно указать делитель - тогда как куда удобнее оперировать с матрицей дробных чисел (скажем, кольцевой фильтр задать удастся только матрицей не целочисленных коэффициентов). Еще можно нарваться на передискретизацию, когда разрешение региона и слоя не совпадают - для полностью автоматизированной обработки потребуется делать множество проверок этой и других ситуаций просто для того, чтобы GRASS не проявлял самодеятельность. Практически это означает, что все надо вручную проверять - что исходные данные и результаты обработки похожи на то, что мы ожидаем, и обработка выполнена корректно.

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

Анонимный комментирует...

Аха, теперь понял, что вы имели в виду, спасибо за подробные разъяснения.

(ну и просто в догонку: на самом деле модули для БПФ нужно искать среди модулей обработки изображений, а не растров :), например тут: http://grass.osgeo.org/grass64/manuals/html64_user/i.fft.html , кстати для меня это одна из нехороших черт -- не известно, где искать нужный модуль, если не знаешь, как он называется)

Алексей Печников комментирует...

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

В то время как математически кольцевой фильтр задать достаточно просто (это можно записать компактнее, но имхо так нагляднее):
octave:1> r=2;
octave:2> f = fspecial("disk", r) -
[zeros(1,2*r+1); zeros(2*r-1,1), fspecial("disk", r-1), zeros(2*r-1,1); zeros(1,2*r+1)]

0.00000 0.00000 0.07692 0.00000 0.00000
0.00000 0.07692 -0.12308 0.07692 0.00000
0.07692 -0.12308 -0.12308 -0.12308 0.07692
0.00000 0.07692 -0.12308 0.07692 0.00000
0.00000 0.00000 0.07692 0.00000 0.00000

Если фильтр задан в пикселях растра, то отсчет=пиксель и мы просто считаем БПФ растра и фильтра или же сразу выполняем фильтрацию соответствующей функцией (через свертку, по умолчанию, но можно и через БПФ). Отсчетов на пиксель можно взять и более одного, но мы выберем целое число отсчетов (степень двойки, обычно)... тогда как GRASS в географических координатах этого сделать не сможет. Из-за всего этого в GRASS мы получим зашумленную картинку и потеряем мелкие объекты... может, они для нас и не значимы, но каждый раз придется вручную просчитывать реальное пространственное разрешение результирующей картинки, так как мы будем наблюдать мелкие объекты, на самом деле представляющие собой шум. Уверен, глубоко ошибочно делать обработку цифрового изображения иначе, чем в пикселях - одни недостатки и усложнения, преимуществ вообще никаких.


(C) Alexey Pechnikov aka MBG, mobigroup.ru