УДК 517.98+621.39
Вестник СПбГУ. Сер. 10, 2007, вып. 1
М. К. Чобану, С. А. Караказъяп
ТРЕХКАНАЛЬНЫЕ МНОГОСКОРОСТНЫЕ СИСТЕМЫ
1. Введение. В работе рассмотрены вопросы реализации методов сжатия изображений с помощью двумерного неразделимого (несепарабелыюго) вейвлет-преобразования. Кодирующая часть основана на переработанном и оптимизированном варианте метода пространственно ориентированных иерархических деревьев (SPIHT).
Базовый кодер изображений можно представить в виде трех основных этапов -декоррелирующее преобразование, процедура квантования и энтропийное кодирование.
До появления кодеров на базе вейвлет-преобразования в качестве декоррелятора использовалось дискретное косинусное преобразование (DCT) [1]. Стандарт JPEG работает с DCT, разбивая при этом изображение па блоки размером 8x8 и выполняя DCT раздельно для каждого блока.
В последнее время самые эффективные из кодеров применяют вейвлет-де-композицию сигнала. Стандартизован кодер JPEG2000, работающий с вейвлет-пре-образованием.
Существует также большое количество нестаидартизованных, но зачастую не менее эффективных алгоритмов. Одним из них является SPIHT. Реализуемый в нем подход уже можно считать классическим. На его базе развиваются новые алгоритмы, и в оценках эффективности любого нового вейвлет-кодера обязательно можно найти сравнительные тесты со SPIHTom.
В настоящее время принято делить кодеры, работающие с вейвлет-декомпозицией, на два класса: межполосные (inter-band) и внутриполосные (intra-band). Алгоритм SPIHT относится к классу межполосных кодеров, т. е. в своей работе он использует избыточность, связанную с корреляцией между уровнями декомпозиции.
2. Сжатие с потерями и без потерь. При сжатии изображений с потерями эти потери должны быть минимально заметны для человека. Сжатие с потерями позволяет достичь сжатия в десятки и даже сотни раз при удовлетворительном качестве. Для большинства изображений, сжатых алгоритмами, рассмотренными в настоящей работе, ошибки при сжатии в 8 раз практически невозможно заметить невооруженным глазом, тогда как классические методы сжатия без потерь, такие как алгоритмы LZ*, позволяют достичь в лучшем случае сжатия в 2-3 раза.
Целью сжатия изображений во многих случаях является компактное хранение данных и/или их последующая передача по каналам связи с максимальной скоростью. При этом в системах, где пропускная способность канала ограничена, а скорость передачи крайне важна, полезным оказывается свойство шкалируемости качества для рассматриваемого семейства алгоритмов. Сначала надо передавать копию изображения с меньшим качеством, но с сохранением самой необходимой информации об изображении, а затем постепенно «уточнять» менее важные детали - это так называемая возможность постепенной передачи (progressive transmission).
Вместе с тем во многих медицинских приложениях сжатие должно обеспечивать возможность анализа больших массивов изображений (наборов трехмерных данных,
Работа выполнена при финансовой поддержке Российского фонда фундаментальных исследований (грант № 06-01-00457) и Японского общества JSPS (грант № 06-07-91751-ЯФ-а).
© М. К. Чобану, С. А. Караказьян, 2007
последовательностей изображений, каталогов) в интерактивном режиме, для исследования зависимых от контекста деталей изображения, а также количественного анализа результатов измерений. Поэтому в данной области существуют жесткие требования к качеству обрабатываемых изображений: недопустимо отбрасывать любые детали при обработке медицинских данных. Неточная их передача может повлиять на диагноз, что нередко влечет за собой тяжелые последствия. Например, детали изображений, полученные с помощью проективной радиографии, могут отображать повреждения внутренних органов, причем различить их можно только по еле заметным изменениям яркости изображения.
Концепция постепенной передачи изображений здесь также крайне важна, так как позволяет быстро определить по полученным ранее фрагментам с малым разрешением важные для постановки диагноза области изображения.
3. Метод пространственно-упорядоченных иерархических деревьев.
3.1. Основы метода. Одним из важных свойств вейвлет-преобразования является его способность перераспределять энергию сигнала, концентрируя ее в малом числе каналов. Исходный сигнал разбивается на субполосы, и энергия концентрируется в одной субполосе, а именно - в низкочастотной (LL - Low-Low, НЧ). Эта компактность энергии ведет к эффективному применению скалярных квантователей. Для квантования областей с низкой энергией отводится меньшее количество бит, что обеспечивает сжатие. Однако такой подход не учитывает остаточную структуру, сохраняющуюся в вейвлет-коэффициентах, в особенности высокочастотных (ВЧ) субполос. В ВЧ субполосах имеются обычно большие области с нулевой или малой энергией. Области с высокой энергией повторяют от субполосы к субполосе свои очертания и местоположение, так как они появляются вокруг контуров в исходном изображении там, где вейвлет-преобразование на высоких НЧ уровнях вейвлет-разложения не может адекватно представить сигнал, что приводит к «утечке» части энергии в ВЧ субполосы. Медленно изменяющиеся гладкие области исходного изображения хорошо описывают НЧ вейвлет-базисы, что вызывает «упаковку» энергии в малом числе коэффициентов НЧ области. Этот процесс примерно повторяется на всех уровнях декомпозиции, что и дает визуальную «похожесть» различных субполос. Априорное знание о том, что изображение состоит из гладких областей, текстур и контуров, помогает учитывать эту межполосную структуру.
Впервые идею использования такого свойства вейвлет-разложения предложили А. Лыоис и Г. Ноулес [2]. Их алгоритм базировался на построении пространственных деревьев над полем вейвлет-коэффициентов.
Дерево выстраивается из выбранной точки, или корня, который может быть расположен в любой субполосе. Корень дает четыре потомка в том же субдиапазоне на предыдущем уровне декомпозиции. Идеи, лежащие в основе алгоритма Лыоиса и Но-улеса, легли в основу многих современных кодеров изображения.
Г. М. Шапиро в 1993 г. [3] разработал очень эффективный и с малой вычислительной сложностью алгоритм для сжатия изображения EZW, который основывался на следующих принципах:
1) частичное упорядочивание компонент изображения по модулю с передачей порядка при помощи алгоритма формирования подмножеств, при этом декодер релизует аналогичный алгоритм;
2) передача битовых плоскостей в соответствии с сортировкой;
3) использование подобия компонент изображения после вейвлет-преобразования.
Рис. 1. Пример вейвлет-декомпозиции (Lena256Gray.bmp, фильтр - 2Ыог2.2, 2 уровня).
В 1996 г. А. Сайд и В. А. Перлманн [4, 5] разработали улучшенный метод кодирования изображений при помощи пространственно упорядоченных иерархических деревьев. Данные алгоритмы являются кодерами битовых плоскостей.
Алгоритмы, работающие с битовыми плоскостями, обладают свойством саморегулируемости. Это означает, что битрейт может быть изменен (уменьшен) после фактического кодирования простым укорачиванием битового потока. Существует несколько типов саморегулируемости битовых потоков в зависимости от того, что происходит с восстанавливаемым изображением при приеме очередной части битового потока.
Если получение большего количества бит приводит к тому, что восстановленное изображение имеет меньшее искажение, то данный битовый поток является саморегулирующимся с последовательным улучшением качества.
Если значительное количество бит приводит изображение большего размера за счет дополнительных высокочастотных субдиапазонов, то такой битовый поток является саморегулирующимся с последовательным увеличением разрешения.
После вейвлет-преобразования (рис. 1) и квантования задача энтропийного кодирования - конвертировать поле кодовых значений в битовый поток и обратно. Одним из таких алгоритмов является ЯРШТ.
3.2. Кодирование битовых плоскостей. Алгоритм за один проход кодирует одну битовую плоскость. Каждое кодовое значение у(пг,пс) хранится в знакопеременном бинарном формате, знак кодового значения - в единственном бите, а модуль - как число с фиксированным количеством бит. Пусть 6 - индекс битовой позиции в модуле кодового значения. Тогда бит с индексом 6 = 0- наименьший значащий.
Пусть и3{пг,пс) - знаковый бит у(пг,пс) и ут(пг,пс, 6) - 6-тый бит модуля ь(пг,пс), тогда 6-тые биты всех кодовых значений ут(пг,пс, 6), где пг = 0,1,..., ЛГг, пс = 0,1,..., Агг, формируют битовую плоскость кодовых значений.
Работа алгоритма начинается с поиска индекса 67П5р, наиболее значащего ненулевой битовой плоскости. Данное число передается декодеру. Битовые плоскости над Ьт$р известны декодеру, так как они равны 0. Алгоритм начинает кодирование наиболее значащего бита с наиболее значащей ненулевой битовой плоскости, т. е. 6 = Ьтзр, и продолжает работу до наименее значащей битовой плоскости, т. е. 6 = 0. Для каждой битовой плоскости кодер помещает кодируемые биты в битовый поток. Важно передать
меньше бит, чем их было в битовой плоскости изначально, иначе ire произойдет компрессии данных. Данный подход обеспечивает саморегулируемость алгоритма за счет последовательного уточнения кодовых значений по мере передачи битовых плоскостей.
3.3. Значимость. Пусть bs(nr, пс) является индексом наиболее значимого ненулевого бита в кодовом значении v(nr,nc) или bs(nr,nc) = —1, если v(nr,nc) = 0. Такая последовательность называется полем значимости. После вейвлет-преобразований существует статистическая зависимость между значениями поля значимости. SPIHT использует такую зависимость для уменьшения объема передаваемых данных. Для этого в алгоритме обработка каждой битовой плоскости разбита на две стадии: определения значимости (significance pass) и уточняющая (refinement pass).
При кодировании в битовой плоскости b кодовое значение значащее, если bs(nr,nc) >= b. В стадии определения значимости кодер передает декодеру, какие кодовые значения впервые стали значимыми. В битовой плоскости Ь существуют кодовые значения.такие, что bs(nr,nc) = b. Они имеют модуль в промежутке [s, 2s), где s - уровень значимости, s = 2Ь.
На данный момент декодер уже знает, какие кодовые слова имеют bs(nr,nc) > b, потому что они были найдены в предыдущих битовых плоскостях. На уточняющей стадии кодер передает 6-тый бит каждого кодового значения, которое было значимым на предыдущих битовых плоскостях. Таким образом, после обеих стадий для битовой плоскости b декодер может определить значение всех битов для этой плоскости и для более значащих битовых плоскостей.
3-4- Построение пространственно цпорядоченных иерархических деревьев. SPIHT-алгоритм группирует кодовые значения в множества, образованные на пространственно упорядоченных иерархических деревьях в вей влет-деком позиции. Пространственные деревья формируются при помощи соотношения «отец/сын», распространяющееся сквозь уровни декомпозиции. Каждое кодовое значение, за исключением трех субдиапазонов наивысшей частоты и тех, которые находятся в низкочастотном субдиапазоне, имеют четырех «детей».
«Дитя» - это 2 х 2-блок кодовых значений в следующем субдиапазоне с более высокой частотой с теми же ориентацией и пространственным расположением. Предполагается, что число строк и столбцов изображения должно быть кратно 2i+1, где L - число уровней декомпозиции. Следовательно, кодовое значение с координатами (пг,пс) имеет четырех «детей» с координатами (2пг,2пс), (2пг,2пс + 1), (2пг + 1,2тгс), (2п,- + 1, 2пс + 1). Чтобы понять, почему используется такая структура, необходимо обратить внимание на визуальную корреляцию амплитуды между поддиапазонами различных уровней вейвлет-декомпозиции на рис. 1.
В используемом алгоритме кодовые значения в верхнем левом квадранте самого последнего уровня декомпозиции не являются корнями деревьев. В других квадрантах данного уровня декомпозиции они представляют собой корни иерархических деревьев в соответствии с определениями, приведенными выше. Это изменение немного упрощает алгоритм, потому что выражения, описывающие отношения «отец/сын», теперь одинаковы для всех поддиапазонов.
Пространственные деревья формируются путем рекурсивного применения отношений «отец/сын» и получения «деревьев-потомков».
Алгоритм и в кодере, и в декодере группирует незначащие кодовые значения в три типа множеств: А - множество «внуков», В - множество «правнуков» и индивидуальное множество (незначащие пиксели).
Рис. 2. Два уровня вейвлет-декомпозиции с применением трехка-иального неразделимого фильтра.
В случае обнаружения алгоритмом значимости множества типа А с корнем в II данное множество разбивается на пять новых множеств: множество типа В с корнем в Я и четыре индивидуальных множества, содержащих «детей» по одному в каждом. Множество типа В с корнем в 7? разбивается на четыре множества типа А с корнями
в «детях» В..
Алгоритм начинает работу с начальной установки набора множеств. Во время инициализации кодер и декодер каждое кодовое значение в низкочастотном субдиапазоне помещают в отдельное индивидуальное множество. Остальные кодовые значения являются членами множеств «потомков». Для каждого значения в низкочастотном суб-диапазоне за исключением верхнего левого квадранта формируется множество типа Л с корнем в этом кодовом значении. Таким образом, после инициализации каждое кодовое значение находится только в одном множестве. Во время работы алгоритма кодовое значение может находиться только в одном множестве. Отметим, что декодер, так же как и кодер, организует кодовые значения в виде множеств. Он уточняет изначально неизвестные кодовые значения.
4. Применение иерархического алгоритма для неразделимых решеток и банков фильтров. Ранее алгоритм БРШТ описывался применительно к вейвлет-разложениям, полученным с применением разделимых фильтров [6, 7]. Свертка с импульсной характеристикой двумерного фильтра этого класса аналогична последовательной свертке с одномерными фильтрами по строкам и столбцам. Данное свойство позволяет использовать разделимую решетку децимации, описываемую диагональной
ф 2 ■ Современные методы построения неразделимых
банков фильтров [8, 9] дают возможность синтезировать системы с некоторыми дополнительными свойствами, например такими как абсолютная симметрия коэффициентов импульсной характеристики и др. Двухканальные лифтинговые схемы существенно повышают скорость выполнения вейвлет-декомпозиции, снижая количество операций.
На рис. 2 представлен пример использования неразделимого банка фильтров. Нечетные уровни декомпозиции имеют ромбическую форму. Это свойство затрудняет практическое применение данного класса фильтров. Такие структуры сложнее хранить в памяти, визуализировать, существенно затруднен их анализ.
Работа двухканальных систем рассматривалась в [4]. Ниже впервые изучается трех-канальпая неразделимая многоскоростная система.
4-1■ Сдвиг. На нечетных уровнях декомпозиции единственно возможный способ группировки отсчетов после децимации - это ромбическая форма. Она не искажает сиг-
матрицей децимации М =
Рис. 3. Пример вейвлет-декомпозиции с использованием неразделимых фильтров, выполненный по предложенной методике.
нал и делает возможным дальнейший анализ. Если, например, сдвинуть («схлопнуть») отсчеты продецимированного сигнала, то он будет искаженным, а его дальнейший анализ неправомерным, так как нарушены геометрические соотношения и относительное расположение коэффициентов. Однако такая структура обладает компактностью.
В настоящей работе предложим довольно простое и эффективное решение этой проблемы. Перед вейвлет-декомпозицией из исходных фильтров дополнительно формируем два фильтра для нечетных и четных строк. При такой организации сигнал визуально искажается только на нечетных уровнях декомпозиции, там, где необходимо ромбическое представление. На четных же уровнях сигнал сохраняет форму и обрабатывается стандартными фильтрами, если выполняется равенство м2 = Уд2 = 21, приводящее к разделимой децимации на четных уровнях.
4-2. Трехканалъные неразделимые системы. Для случая трехканальных
неразделимых систем выбрана матрица децимации М
2 1
-1 -2
. Модуль опре-
делителя т = |скЧМ| = 3 равен числу каналов и отношению частот следования отсчетов исходного и продецимированного сигналов. При этом выполняется необходимое свойство М2 = ЗТ2, что необходимо для построения эффективных систем кодирования многомерных сигналов [10, 11].
Для практического применения данного класса банков фильтров нужно преобразовать внешний вид субполос к виду, приемлемому для обработки алгоритмом ЭРШТ. Для возможности применения алгоритма ЯРШТ субполосы, получаемые в результате трехканальной фильтрации, преобразуются в прямоугольники, а также располагаются в виде древовидной мозаики для последующего обхода этого разложения алгоритмом БРШТ. После подобной компоновки субполос необходимо предисказить импульсную характеристику фильтра.
Следует заметить, что в результате подобного преобразования в итоге получается три предискаженных фильтра из одного базового - один для первых столбцов, другой -для вторых, третий - для третьих. После проведения всех преобразований происходит компоновка всех полученных субполос разложения в виде, приемлемом для алгоритма БРШТ.
На полученном разложении алгоритм 8РШТ строит свою работу аналогично разделимому варианту, используя подобие ВЧ субполос. Такой подход позволил не только применять БРПГТ, но и дал еще целый ряд преимуществ: сделал более доступным ненулевое расширение сигнала при фильтрации (например, периодическое, реализация которого на ромбических структурах технически сложная задача); существенно снизил объем требуемой памяти, упростил визуализацию и обработку вейвлет-разложения (рис. 3).
5. Синтез трехканальных неразделимых фильтров. Хорошо известно, что для обеспечения качества изображения фильтры должны удовлетворять так называемому правилу сумм (условию обнуления моментов), которое обеспечивает соответствующий порядок аппроксимации (см. [12]) и является необходимым условием гладкости масштабирующих функций (см., например, [13]). Метод построения таких фильтров был разработан М. А. Скопиной [9]. Ею был найден полифазный критерий обнуления моментов, который эквивалентен правилу сумм, но, в отличие от него и других ранее известных критериев, позволил выписать явно все фильтры, обладающие этим свойством, для произвольной матрицы децимации. Более того, в [9] для произвольного интерполяционного фильтра даны явные формулы для двойственного фильтра, который удовлетворяет правилу сумм того же порядка, что и исходный. На базе полученных в [9] результатов нами разработан метод симметричных/антисимметричных двойственных фильтров, удовлетворяющих правилу сумм порядка п, для данной матрицы децимации с нечетным определителем. В настоящей работе реализуем этот метод " 2 1
для матрицы М
-1-2
и п — 2. Из [14] следует, что для нахождения требуемых
фильтров можно взять следующие порождающие матричные полифазные функции:
а(:г, у) = —2 вт 2тт(2х — у) + эт 2ж(х — 2г/),
V =
1~Шсг4 1 - 1/18сг2 + 1/Зг ст 1 — 1/18 сг2 — 1/3 г,
-1 + 1/18 а2 1/3 г а 1
л/3
1/2 1/2
1/2 -1/2
1 - 1/18ст2 + 1/3 г
-2 + 1/9ст2 1 + 2/9сг2 — у^ст4 — 2/3г а + 1/27г ст3 2/3 г а 3 - 2/9 сг2 + 2/Зг" <т -1/27I а3
1
1 + 2/9 а2
1/18сг2 - 1/3?: СГ -¡^СГ4 +2/3 ¿сг - 1 /27 г <т3
162
-3 + 2/9 сг2 + 2/3 г ст - 1/27 г ст3
По данным полифазным функциям можно найти импульсные характеристики (маски) фильтров банка анализа (БА):
С 0 =
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 1 81 1/18
0 0 0 0 0 0 0 0 0 0 7 162
0 0 0 0 0 0 0 0 0 -1/3 1/18
0 0 0 0 0 0 1 108 1/18 0 0 1/6
0 0 0 0 0 0 0 5 108 31 36 1/3 0
0 0 0 0 0 0 -1/6 1/18 133 144 1/18 -1/6
0 0 . 0 1 324 1 72 0 0 1/3 31 36 5 108 0
0 0 0 0 13 648 — 1/18 1/6 0 0 1/18 1 108
0 0 0 0 1 72 5 108 1/18 -1/3 0 0 0
1 2592 0 0 0 0 -1/18 7 162 0 0 0 0
0 1 324 0 0 0 0 1/18 1 81 0 0 0
0 0 1 108 0 0 0 0 0 0 0 0
0 0 0 1 81 0 0 0 0 0 0 0
0 0 0 0 1 162 0 0 0 0 0 0
о о
о о
0 о
1 /18 О
1
102 О
108 1/18 0 0 О О О О О 0 О О
X
72
13 648
О 0 О
о о
__I
324
о о о о о о о о о
ООО ООО О О
" 108 о о о о о о о о о о о о о о
324
о о о о о о о о о о о о о
2592
о о о о о о о о о о о о
С1 =
0 0 0 0 0 0 -1/18 0 0
0 0 0 0 0 0 0 1/18 0
0 0 0 0 0 0 0 0 1 72
0 0 0 -1/18 1/2 0 0 0 0
0 0 0 0 31 36 0 0 0 0
0 0 0 0 1/2 -1/18 0 0 0
1 72 0 0 0 0 0 0 0 0
0 1/18 0 0 0 0 0 0 0
0 0 1/18 0 0 0 0 0 0
0 0 0 1/3 0
0 0 -1/2 0 - -1/6
С2 = 0 0 0 0 0
1/6 0 1/2 0 0
0 -1/3 0 0 0
Аналогично определяются маски банка синтеза (БС). Частотные характеристики фильтров БА и БС показаны на рис. 4.
в / Рис. Фильтры банков анализа и синтеза.
а, с, е - фильтры БА СО, БА С1, БА С2 соответственно; 6, й, / - фильтры БС ОО, БС 01, БС Б2 соответственно.
6. Заключение. В ходе работы описанный принцип фильтрации был реализован в среде Ма^аЬ. Фильтр преобразуется в специальной процедуре, т. е. разработчик фильтров избавлен от перегруппировки коэффициентов.
В настоящее время такой подход используется для анализа синтезируемых филь-
тров, проведения ряда тестов и делает доступным использование неразделимых фильтров в системах сжатия изображений.
Summary
Tchobanou М. К., Karakazjan S. A. Three channel multirate systems.
The method of synthesis and implementation of image compression techniques based on two-dimensional nonseparable wavelet transform is considered. Synthesis of the filters is made by application of the method symmetrical/antisymmetrical dual filter for a given sum rule of order n for given decimation matrix with an odd determinant. Coding part is based on an improved and optimized version of set partitioning into hierarchical trees.
Литература
1. Зубарев Ю. Б., Дворкович В. П. Цифровая обработка телевизионных и компьютерных изображений. М.:' Изд. Междунар. центра науч. и технич. информации, 1997. 212 с.
2. Lewis A., Knowles G. Image compression using the 2-d wavelet transform // IEEE Trans. Image Proc. 1992. Vol. 2. P. 244-250.
3. Shapiro J. M. Embedded image coding using zerotrees of wavelets coefficients // IEEE Trans. Signal Proc. 1993. Vol. 41. P. 3445-3462.
4. Said A., Pearlman W. A. A new fast and efficient image codec based on set partitioning in hierarchical trees // IEEE Trans. Circ.., Syst. for Video Technol. 1996. Vol. 6. P. 243-250.
5. Wheeler F. W., Pearlman W. A. Spiht image compression without lists // IEEE Intern. Conf. on Acoustics, Speech and Signal Processing (ICASSP). 2000. Vol. 6. P. 5-9.
6. Черников А. В., Чобану M. К. Современный метод сжатия изображений на базе вейвлет-преобразования и иерархического алгоритма кодирования // Цифровая обработка сигналов. 2005. № 3. С. 40-59.
7. Чобану М. К. Многомерные многоскоростные системы и многомерные вейвлет-функции. Ч. II. Синтез // Вестн. Моск. энергетич. ин-та. 2003. Т. 3. С. 69-78.
8. Батлук А. В., Чобану М. К. Исследование применения банков фильтров для сжатия изображений // Цифровая обработка сигналов. 2005. № 4. С. 29-40.
9. Skopina М. On construction of multivariate wavelets with vanishing moments // Appl. and Сотр. Harmonic Analysis. 2006. Vol. 20, N 3. P. 375-390.
10. Чобану M. К. Многомерные неразделимые матрицы децимации // Труды 13-го Междунар. конгресса Академии информатизации ITS-2005. М.: СТАНКИН, 2005. Т. 1. С. 67-70.
11. Tchobanou М. Parameterization of multidimensional decimation matrices // Proc. the 2005 Intern, workshop on spectral methods and multirate signal processing, SMMSP-2005. Riga, Latvia, 2005. P. 7-10.
12. Jia R. Q. Approximation properties of multivariate wavelets//Math. Сотр. 1998. Vol.67. P. 647-655.
13. Han B. Analysis and construction of optimal multivariate biorthogonal wavelets with compact support // SIAM J. Math. Anal. 1999. Vol. 31, N 2. P. 274-304.
14. Каракаэьян С. А. О построении двумерных симметричных/антисимметричных всплеск-функций // Вестн. С.-Петерб. ун-та. Сер. 10: Прикладная математика, информатика, процессы управления. 2006. Вып. 4. С. 61-69.
Статья представлена к публикации членом редколлегии С. В. Чистяковым.
Статья принята к печати 18 сентября 2006 г.