CC BY
78
УДК 616-073.75:681.32
С. В. Меркурьев, канд. физ.-мат. наук, НИПК «Электрон», Санкт-Петербург
Методы подавления шума цифровых рентгенограмм
Ключевые слова: цифровые рентгенограммы, интенсивность сигнала, визуальная оценка качества, уровень дисперсии шума
В работе рассматриваются различные тенденции в исследованиях методов подавления шума цифровых изображений. Описаны уже получившие практическое применение подходы, а также представляющиеся перспективными способы фильтрации шумов цифровых рентгенограмм. Проведено сравнение некоторых методик подавления шума и показано, что объединение оценки шума, зависящего от интенсивности сигнала, с алгоритмами шумоподавления, в особенности с усредняющими в пространстве исходного изображения фильтрами, обеспечивают хорошие результаты с точки зрения числовых метрик и визуальной оценки качества.
Введение
На данный момент времени цифровая рентгенография является распространенным методом диагностики различных заболеваний. В клинической практике для коррекции и облегчения интерпретации цифровых снимков применяется широкий спектр методов и алгоритмов обработки: подавление шума, повышение резкости, автоматическое выделение зоны интереса и т. д. Важнейшими среди алгоритмов обработки цифровых снимков в рентгенографии являются методы подавления шумов. Подавление шума улучшает визуальное восприятие и увеличивает возможность различения низко-контрастных деталей, что делает анализ снимка рентгенологом более удобным и повышает качество выполнения других алгоритмов обработки изображений.
Настоящая статья является обзорной по методам оценки и подавления шума цифровых снимков. Особый акцент сделан на методах, уже хорошо себя зарекомендовавших на практике либо имеющих большой потенциал применения в обработке цифровых рентгенограмм.
1. Оценка шума цифровых снимков
Многие эффективные методы подавления шума использует информацию об уровне шума на изображении. Поэтому до применения того или иного мето-
да фильтрации шума обычно встает следующая задача: располагая лишь исходным снимком, определить уровень дисперсии шума. Данная задача осложняется тем, что при рассмотрении цифровых снимков и, в частности, снимков проекционной цифровой рентгенографии уровень шума на изображении зависит от интенсивности сигнала.
Под шумом обычно понимают случайное отклонение измеренного значения некоторой физической величины от ее точного значения. В цифровой рентгенографии детектор измеряет интенсивность рентгеновского излучения — каждая ячейка работает как сумматор фотонов, накапливая за время экспонирования в среднем N электронов путем поглощения фотонов. В силу физической природы излучения количество накопленных электронов является величиной случайной, распределенной по закону Пуассона:
- N
р N = п) = ехр(-N) , п > 0.
п!
Случайные флуктуации числа поглощенных фотонов называют квантовым шумом или шумом фотонов. В современных детекторах основным источником шума является шум фотонов. К дополнительным источникам шума относят шумы детекторной системы: шум чтения, тепловой шум, шумы усилителей, шумы квантования и др. Суммарный эффект данных источников шума можно моделировать распределенной по закону Гаусса случайной величиной [1-2]. Общепринятой является модель, согласно которой при линейных электронных схемах дисперсия шума пикселя выходного изображения линейно растет со значением накопленного в данном пикселе сигнала [3]:
ст2!(р) = а!(р) + Ь , (1)
где 1(р) — уровень интенсивности сигнала в пикселе р.
Проблеме оценки и подавления зависящего от сигнала шума цифровых снимков посвящено много публикаций [1, 2, 4]. Так, в работе [4] рассматривается алгоритм оценки шума и исследуются различные схемы подавления шума в рентгенографии. Особый акцент делается на построении алгоритма,
биотехносфера
| № 4(10)/2010
способного работать в реальном времени. Алгоритм оценки шума является непараметрическим, т. е. не учитывает модель (1). В публикации [2] рассматривается двухпараметрический подход к оценке шума цифровых снимков. Шум исходного цифрового изображения (полученного непосредственно с сенсора и не прошедшего через нелинейные преобразования, такие как гамма-коррекция и т. п.) моделируется как аддитивная по отношению к сигналу случайная величина, дисперсия которой зависит от сигнала согласно закону (1). Предлагается способ построения модельных кривых дисперсии шума, учитывающий нелинейности в работе сенсора, приводящие к недоэкспонированию и переэкспонированию, т. е. при нарушении закона (1) на краях динамического диапазона. Исследуются возможные схемы подавления шума, использующие полученные на этапе оценки шума данные.
Суммируя информацию, изложенную по данной теме в различных публикациях, можно заключить, что для построения оценки зависящего от сигнала шума по исходному снимку достаточно выполнить следующую процедуру:
1) получить оценочный снимок с помощью низкочастотной фильтрации (линейной, медианной и др.) исходного снимка и построить снимок шума (разность между исходным и оценочным снимками);
2) отбросить тем или иным способом выделяющиеся значения пикселей снимка шума, соответствующие резким изменениям (границы, одиночные «горячие» пиксели) в исходном снимке;
3) разбить диапазон интенсивности оценочного снимка на интервалы: каждый пиксель оценочного снимка принадлежит определенному интервалу;
4) накопить для каждого интервала пиксели снимка шума, соответствующие пикселям оценочного снимка;
5) вычислить дисперсию шума в каждом интервале интенсивности по накопленным в нем пикселям шумового изображения.
После выполнения данных этапов по полученным интервальным оценкам шума можно построить непараметрическую оценку шума (например, как в работе [4]) — создать интерполяционную табличную функцию, позволяющую каждой яркости исходного снимка поставить в соответствие оценку дисперсии шума. Точками входа в таблицу могут служить, например, интенсивности оценочного снимка. Иначе, в случае параметрической оценки шума, можно построить тем или иным способом (например, методами наименьших квадратов, минимизации функции правдоподобия, направленной оптимизации) оценку параметров а, Ь модели шума (1). Возможен также учет нелинейностей сенсора [2].
Для использования найденных параметров (будь то интерполирующая таблица или параметры а, Ь модели шума) на основе оценочного снимка можно
Рис. 1\ Цифровая рентгенограмма руки
построить карту шума — изображение, каждый пиксель которого оценивает среднеквадратическое отклонение (СКО) шума в соответствующем пикселе исходного снимка. Карта шума дает попиксель-ную оценку шума с достаточной для практического применения точностью.
С другой стороны, при параметрической оценке шума возможно применение преобразования, стабилизирующего дисперсию шума исходного снимка [3], т. е. такого преобразования, которое сводит проблему фильтрации шума, зависящего от сигнала, к задаче подавления аддитивного независимого от сигнала шума постоянной (заданной) дисперсии. В любом случае использование как можно более точной информации об уровне шума на изображении позволяет существенно улучшить качество алгоритмов шумоподавления, что продемонстрировано в ряде работ [1-2].
Пример вычисления параметров модели шума (1) для реального цифрового снимка руки (рис. 1) приведен на рис. 2, демонстрирующем хорошее совпадение теоретической модели (1) с экспериментальными данными.
а(1)
80
70
60
50
40
30
г 'О 'сР
(1
= 7 а1 1 Ь
Ог.фР '■ръ
5000 5500 6000 6500 7000 7500 8000 8500 9000
Рис. 2
Зависимости СКО шума от интенсивности для снимка рис. 1. Параметры модели шума найдены методом наименьших квадратов:
о — интервальные оценки СКО; - — модельная кривая;
а = 1,4в,Ь = -695в,3
№ 4(Ю)/2010~|'
биотехносфера
2. Методы подавления шума цифровых снимков
К настоящему времени разработан и используется весьма разнообразный математический инструментарий, непосредственно относящийся к проблеме подавления шума цифровых изображений. Простое перечисление с указанием характерных свойств тех или иных методов шумоподавления, да еще и с учетом разнообразных улучшений и всевозможных гибридных схем, займет целую книгу и представляется невозможным в рамках статьи. Здесь мы рассмотрим две группы методов шумоподавления:
• методы фильтрации в пространстве преобразования;
• методы фильтрации в пространстве исходного изображения.
Первая группа методов сводится к трансформации исходного изображения в иную форму посредством применения некоторого преобразования. Затем следует обработка коэффициентов преобразования по определенному правилу, например с помощью специальным образом сконструированного нелинейного оператора. Целями такой обработки могут являться уменьшение энергии шума и, возможно, увеличение резкости снимка. Завершает фильтрацию реконструкция изображения с помощью обратного преобразования. К данной группе можно отнести весьма широкий класс методов, основанных на преобразовании Фурье, вейв-лет-преобразовании, комплексном вейвлет-преобразо-вании, курвлет-преобразовании и др. Качество и скорость работы фильтра определяются выбранным преобразованием и используемым правилом обработки.
Ко второй группе методов относится не менее обширный список подходов: анизотропное размытие, вариационные методики, фильтр Ярославского, билатеральный фильтр, регрессионный фильтр, фильтр с оптимальной пространственной адаптацией и др. Как правило, данные методы сводятся к поиску подобия между пикселями или группами пикселей (пат-чами) в самом изображении, возможно с учетом геометрической информации (такой как локальная ориентированность структур изображения), и применению той или иной схемы модификации (усреднению) яркостей этих пикселей, при этом сила модификации зависит от степени подобия пикселей (патчей). Особо отметим изобретенный недавно фильтр нелокальных средних [6], демонстрирующий хорошие результаты в подавлении шума изображений и породивший большой поток публикаций, касающихся различных аспектов его улучшения и практического применения.
3. Методы фильтрации
в пространстве преобразования
Как следует из названия данной группы методов, они основаны на идее анализа снимка после выполнения его преобразования в некоторое иное
представление — переход в пространство преобразования. После выполнения преобразования происходит обработка коэффициентов преобразования (например, обнаружение и подавление шумовых составляющих, подъем резкости посредством усиления полезных коэффициентов и др.). Обработанные коэффициенты подвергаются реконструкции — обратному преобразованию.
В последние три десятилетия глубокое теоретическое обоснование и широкое практическое распространение получили методы анализа сигналов на основе вейвлет-преобразования (ВП). Общим для преобразований Фурье и ВП является разложение сигнала по базису некоторой анализирующей функции (комплексной экспоненты в Фурье-анализе, локальной волне вейвлет-анализа). Отличия существенны: ВП способно осуществить кратно-масштабный локализованный анализ сигнала — не только выделить частоты в сигнале (в вейвлет-анализе чаще говорят не о частотах, а о деталях сигнала определенного размера или масштаба), но и указать, в какие моменты времени (пространства при анализе изображений) эти частоты присутствуют. Поэтому ВП, в отличие от преобразования Фурье, особенно хорошо подходит для анализа нестационарных сигналов. Дискретное ВП относительно легко сконструировать и реализовать через быстрые алгоритмы линейной сложности.
В настоящее время существует большое разнообразие алгоритмов ВП (ортогональные преобразования, преобразования высокой плотности, комплексные преобразования и др.). Выбор преобразования, обладающего специальными свойствами, главным образом и определяет качество и скорость фильтра. В задаче подавления шума изображений особенно важны такие характеристики ВП, как:
• локальность анализа по пространству и масштабу: локальность обеспечивается выбором вейв-лета, энергия которого как можно лучше сконцентрирована по времени и частоте;
• число нулевых моментов вейвлета: при выборе вейвлета с большим числом нулевых моментов энергия кусочно-полиномиального сигнала лучше концентрируется в небольшом числе коэффициентов ВП, отвечающих резким изменениям сигнала;
• инвариантность ВП по отношению к сдвигу сигнала: данное свойство обеспечивает независимость отклика фильтра от положения сингулярности в сигнале, что дает более точную информацию о частотной локализации резких изменений;
• выделение как можно большего числа направлений: чем лучше и больше с помощью преобразования выделяется направлений, тем проще моделировать и обрабатывать различные геометрические структуры в изображении, обычное ВП выделяет четыре направления (вертикальное, горизонтальное и два смешанных в одном субдиапазоне диагональных направлений);
биотехносфера
| № 4(10)/200
• плотность выборки по пространству и масштабу. преобразования более высокой плотности обеспечивают большую избыточность выборки частотно-временного пространства и аппроксимируют свойство инвариантности по отношению к сдвигу.
На практике необходимо соблюдать разумный баланс по соблюдению отмеченных выше свойств ВП, поскольку их совместное выполнение затруднительно (например, чем выше число нулевых моментов, тем больше носитель вейвлета — меньше локальность) и приводит к увеличению времени преобразования (инвариантное по отношению к сдвигу ВП выполняется дольше, чем обычное ортогональное ВП с прореживанием). Комплексные ВП и в особенности комплексные ВП двойной плотности дают хороший баланс по обеспечению перечисленных свойств и обладают некоторыми дополнительными полезными свойствами.
Комбинируя ВП с обработкой коэффициентов посредством специального оператора, можно получить весьма неплохие по качеству и эффективные по скорости выполнения алгоритмы подавления шума и повышения резкости цифровых медицинских снимков. Подходящий для этих целей оператор не должен приводить к существенным потерям важной диагностической информации или трудностям в ее различении на финальном снимке. Для этого оператор конструируют таким образом, чтобы он был адаптивным к уровню шума на изображении и обеспечивал более сильное подчеркивание низкоконтрастных деталей по сравнению с высококонтрастными. Пример подобного оператора, представляющего собой склейку из кусочно-полиномиальных функций (1 и 2 степеней) вейвлет-коэффи-циента ш, показан на рис. 3.
Детализирующие вейвлет-коэффициенты ш подвергаются обработке оператором (2), осуществляющим пороговое подавление шума и повышение резкости деталей изображения. Приведем формулу такого оператора.
О(ш) 120
а = 10, Т = За, Т2 = 4а, Т3 = 8а
100 80 60 40 20 0
Т3
20
40
60
80
100 120 ш
Рис. 3
Пример оператора подавления шума и повышения резкости:
—— — функция 0(и>);---— пороги Тр Т2, Т3
в (ш) =
Щш, |ш| < Т{;
(аш + р)ш, Т < И < Т;
(уш + д)ш, т2* < И < Т; ш, И > Т ■
(2)
Здесь а =
грк грк
Р = Щ - аТ?; У =
1 - Щ ;
А = 1 - уТ£; к = 1, Ь, Ь — число уровней ВП; Т, Т2 , Т —пороговые значения.
Параметр 0 < Их < 1 позволяет регулировать степень подавления коэффициентов ВП (процент подавления шума), а параметр - 1 задает силу повышения резкости деталей изображения. С каждым уровнем вниз по иерархии дерева ВП значения порогов Т масштабируются по следующему закону.
Т = (к)а, I = 1,3, к = 1, Ь,
где — заданные априори множители (например, ах = 3,0, а2 = 4,0, а3 = 8,0); /(А:) = [0, 1] — специфичная для каждого ВП функция, характеризующая уменьшение шума при переходе к следующему, более грубому уровню ВП; а — СКО шума на изображении, которое можно определить, используя методики раздела 2 настоящей статьи.
4. Методы фильтрации в пространстве изображения
К простейшим пространственным методам фильтрации можно отнести линейные фильтры (такие как усреднение по окрестности, гауссово размытие, биномиальные фильтры и т. п.). Такие фильтры просты в реализации и эффективны по скорости выполнения, однако приводят к искажениям полезной информации, в частности к размытию границ объектов и текстур. Более сложную группу образуют нелинейные фильтры (медианный фильтр, фильтр Ярославского, билатеральный фильтр, фильтр нелокальных средних и др.). Откликом таких фильтров в каждом пикселе снимка является значение, нелинейно зависящее от значений пикселей в некоторой окрестности (области усреднения) данного пикселя. Нелинейные фильтры лучше сохраняют полезные детали изображения и могут учитывать контекст каждого пикселя, однако они менее эффективны по скорости и труднее в реализации.
Разработанный в начале 2000-х годов фильтр нелокальных средних (НЛС) [6] в определенном смысле близок билатеральному фильтру. Билатеральный фильтр каждому пикселю из области усреднения текущего пикселя присваивает вес, зависящий от фотометрического и пространственного расстояний между текущим пикселем и пикселем сравнения [5]. Таким образом, билатеральный
№ 4(!0)/2010 |
биотехносфера
фильтр использует принцип избыточности изображений: в любом реальном изображении найдется много пикселей, подобных усредняемому пикселю. Фильтр НЛС использует более тонкий критерий близости пикселей: каждому пикселю из области усреднения текущего пикселя присваивается вес, зависящий от евклидова расстояния между патчем (небольшой окрестностью) данного пикселя и патчем текущего пикселя, т. е. в данном фильтре используется контекстная избыточность изображений. При этом фильтр НЛС не налагает ограничений на размер области усреднения — она может быть и нелокальной, т. е. совпадать со всем снимком (отсюда и название фильтра). Однако на практике столь большая окрестность усреднения приведет к чрезмерным вычислительным затратам и пересглаживанию, поэтому все-таки ограничиваются локальной областью усреднения.
Несмотря на локальность области усреднения, арифметическая сложность фильтра НЛС весьма высока, что поначалу являлось, пожалуй, главным фактором, препятствующим широкому внедрению данного фильтра на практике. Со временем, однако, рост производительности компьютеров и очевидные преимущества в качестве обработки изображений, которые демонстрирует фильтр НЛС по сравнению со многими другими подходами, привели к интенсивным поискам возможных путей ускорения и дальнейшего улучшения качества работы данного фильтра [6]. На сегодня в литературе уже описаны применения фильтра НЛС в обработке 2Б и 3Б МШ-изображений [7], низкотемпературной электронной микроскопии, внедрение в коммерческое ПО для профессиональной цифровой фотографии [8] и др. Весьма перспективной представляется возможность применения фильтра НЛС и при обработке цифровых рентгенограмм, особенно при объединении данного фильтра с оценкой шума по одной из методик раздела 1. Пример такого применения НЛС-подобного фильтра совместно с оценкой зависимого от сигнала шума при обработке рентгеновского снимка приведен в следующем разделе.
5. Сравнение методов
Приведем иллюстративное сравнение следующих подходов к подавлению шума реального цифрового рентгеновского снимка: пороговая обработка в избыточном базисе (фрейме на вейвлетах Хаара) — инвариантном по отношению к сдвигу ВП без прореживания (алгоритм А1); пороговая обработка коэффициентов комплексного ВП двойной плотности [8е1евшк] (алгоритм А2); фильтрация НЛС-по-добным фильтром (алгоритм А3). Для получения числовых метрик качества (РйМК — пиковое отношение сигнал/шум) реальный цифровой снимок (радиологическая тест-таблица) был отфильтрован с помощью алгоритма А3 на 100 % (таким обра-
зом, получен снимок без шума — оценочный снимок). На оценочный снимок программно наложен аддитивный шум с зависящей от яркости дисперсией согласно модели (1) с параметрами а = 1,25, Ь = 1024, таким образом, получен зашумленный снимок. В качестве пороговой функции в методах А1 и А2 использовался оператор вида (2) с параметрами = 0,2 (что соответствует 80%-ному подавлению шума, или уменьшению шума в 5 раз), . = 1,0 (без подъема резкости). При подавлении шума во всех трех методах использовалась приближенная попиксельная оценка шума (карта шума), полученная согласно методике раздела 1.
На рис. 4 представлен уменьшенный фрагмент зашумленного изображения рентгеновской тест-таблицы (рис. 5). Тест-таблица содержит увеличенные фрагменты данного изображения, в которых яркость—контрастность выставлены таким образом, чтобы подчеркнуть особенности работы каждого из алгоритмов А1—А3 при обработке границ объектов. Визуальный анализ показывает, что НЛС-фильтр гораздо лучше справляется с задачей сохранения границ по сравнению с пороговой обработкой А1—А2 в вейвлет-фрейме, которая приводит, особенно на больших уровнях подавления шума (скажем, при подавлении более чем в 2 раза), к весьма заметным артефактам на резких изменениях в яркости.
Тест-таблица на рис. 6 демонстрирует поведение алгоритмов А1, А2 на тонких структурах, характеризующих пространственное разрешение системы. Здесь алгоритм А3 также показывает значительно более высокую способность к передаче мелких структур без привнесения сколько-нибудь заметных артефактов. Для лучшего визуального представления фрагменты были обработаны фильтром повышения резкости.
Рис. 4
Фрагмент изображения рентгеновской тест-таблицы (зашумленный снимок)
биотехносфера
I № 4(10)/200
Рис. 5
Реализация алгоритмов А1, А2 и АЗ на участке тест-таблицы: а — исходный снимок; б—г — обработка
Рис. 6 Реализация алгоритмов А1, А2 и АЗ на штриховой мире: а — исходный сни-| мок; б—г — снимки обработаны по методам А1, А2 и АЗ соответственно
№ 4(10)/20Ю |
биотехносфера
Наконец, ниже приведены числовые метрики качества рассматриваемых здесь алгоритмов, из которой видно преимущество НЛС-фильтра и с количественной точки зрения.
Метрики РЗЫИ для методов А1 —АЗ (при подавлении шума в 5 раз)
AI A2 A3
67,30 69,17 69,64
Заключение
В данной статье рассмотрены некоторые направления в исследованиях по проблематике оценки и подавления шума цифровых изображений на примере проверенных временем, а также изобретенных сравнительно недавно подходов. Приведен метод оценки по исходному изображению шума, зависящего от интенсивности сигнала. Показано, что объединение оценки шума цифровых снимков (нахождение зависимости шума от интенсивности) с нелокальной фильтрацией в пространстве изображения обеспечивает весьма впечатляющие результаты (как с точки зрения числовых метрик, так и при визуальной оценке качества) и имеет большой потенциал применения в обработке цифровых рентгенограмм. Однако можно утверждать, что, несмотря на хорошие результаты, достигнутые в последние десятилетия, проблема улучшения существующих
и разработки новых, универсальных по применимости к широкому классу изображений, а также быстрых алгоритмов подавления шума, особенно в свете намечающегося перехода к технологии больших плоско-панельных детекторов, и по сей день остается актуальной.
| Л и т е р а т у р а |
1. Bosco A., Battiato S., Bruna A. R., Rizzo R. Noise reduction for cfa image sensors exploiting hvs behavior. Sensors 9(3). 2009. P. 1692-1713.
2. Foi A., Trimeche M., Katkovnik V., Egiazarian K. Practical poissonian-gaussian noise modeling and fitting for singleimage raw-data. Image Processing, IEEE Transactions on 17. (October 2008). P. 1737-1754.
3. Яне Б. Цифровая обработка изображений. М.: Техносфера, 2007. 583 с.
4. Hensel M., Lundt B., Pralow T., Grigat R.-R. Robust and Fast Estimation of Signal-Dependent Noise in Medical X-ray Image Sequences. In Handels, H., et al., ed.: Bildverarbeitung fur die Medizin 2006: Algorithmen, Systeme, Anwendun-gen. Springer. 2006. P. 46-50.
5. Tomasi C., Manduchi R. Bilateral filtering for gray and color images. In Proc. of the Sixth Int. Conf. on Computer Vision. 1998. P. 839-846.
6. Buades A., Coll B., Morel J.-M. On image denoising methods. SIAM Multiscale Modeling and Simulation. 2005. 4(2). P. 490-530.
7. Thaipanich T., Kuo C.-C. Jay. An Adaptive Nonlocal Means Scheme for Medical Image Denoising. SPIE Medical Imaging, San Diego, CA, USA, February 2010.
8. DxO. An Introduction To DxO Noise Engine, October 2005. http://www.sironline.com/press/pdf/dolabs/100405/ DxO_Noise_Engine_Introduction.pdf
ИЗДАТЕЛЬСТВО «ПОЛИТЕХНИКА» ПРЕДЛАГАЕТ
Семенова М. В., Ястребов А. С. Нематериальные ресурсы. Оценка и правовое регулирование в процессе коммерциализации. — СПб. : Политехника, 2011. — 444 с. : ил. Цена 451 руб.
ISBN 978-5-7325-0952-6
В книге рассматриваются имущественные (право пользования природными ресурсами, земельными участками, имуществом; авторское право; смежные права), неимущественные (интеллектуальные) и конкурентные (торговая марка, деловая репутация предприятия, фирменные наименование и стиль, корпоративные коммуникации и культура) нематериальные ресурсы хозяйствующих субъектов, их роль и место в структуре капитала. Дана оценка нематериальных ресурсов хозяйствующего субъекта. На основе формирования инвестиционных предпочтений хозяйствующих субъектов в отношении нематериальных ресурсов рассматривается инновационный путь их развития.
Реальную реализацию инвестиционно-инновационной идеи иллюстрирует комплексная программа по разработке проекта и внедрению почтовой пломбы нового поколения с повышенной надежностью защиты почтовых отправлений на сетях почтовой связи ФГУП «Почта России».
Книга рассчитана на научных работников и специалистов, занимающихся нематериальными ресурсами хозяйствующих субъектов.
Книгу можно приобрести в издательстве по адресу: 191023, Санкт-Петербург, Инженерная ул., д. 6, 3-й этаж.
Тел./факс: 312-44-95, 312-53-90, тел.: 571-61-44 E-mail: sales@polytechnics.ru, masha@polytechnics.ru, gfm@polytechnics.spb.ru
биотехносфера
| № 4С10)/2010
