Вернуться в оглавление        Предыдущий параграф       Следующий параграф

 

3.2.2. Анализ изображения при помощи вейвлет-базиса

 

            Следующий пример демонстрирует использование вейвлет-преобразования, для сжатия двумерного изображения. Файл swphoto.eps содержит сканированное изображение.

 

Загружаем пакет.

In[1]:= Needs["Wavelets`Wavelets`"]

 

Слежующая команда считывает информацию в файле и преобразует её в 8-битовую числовую матрицу.

In[2]:= data = Reverse[Partition[ToExpression /@
(StringJoin[Join[{"16^^"}, #]]& /@
Partition[Characters[StringJoin @@ Drop[Drop[
ReadList[ToFileName[{"Wavelets", "Data"}, "swphoto.eps"], String],
22], -3]], 2]), 250]];

 

            Проверяем размер данных.

In[3]:= Dimensions[data]

Out[3]=

 

            Даем дополнительные строки наверху и дополняем тремя  дополнительными столбцами, два слева и один справа так, чтобы данные имели бы размерность, вида 2N  .

In[4]:= (data = Transpose[Join[Table[#[[1]], {3}], #,
Table[#[[-1]], {3}]]&[Transpose[Drop[data, 21]]]];

Dimensions[data])

Out[4]=

 

            Мы построим график из файла, содержащий портрет ребенка.

In[5]:= ListDensityPlot[data, Mesh -> False]

            Мы используем наименее асимметричный фильтр порядка 4. Здесь мы изменили начальный индекс от 0 до  так, чтобы продолжение от границы было симметричным.

In[6]:= (s4 = LeastAsymmetricFilter[4];
s4 = MapAt[-3&, s4, {2}])

Out[6]=

 

            Данные преобразованы, используя WaveletTransform.

In[7]:= wtdata = WaveletTransform[data, s4, 3,
BoundaryCondition -> Fixed];

 

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

In[8]:= Dimensions /@ wtdata

Out[8]=

 

Мы можем построить график коэффициентов одномерного вейвлет-преобразования в фазовом пространстве используя PhaseSpacePlot или графики коэффициентов  по  для каждого уровня  и расположить вместе графики для различных уровней используя PlotCoefficients. Для двумерных вейвлет-коэффициентов, мы располагаем матрицы различных уровней в большую матрицу размера первоначальных данных и строим их графики,  используя PlotCoefficients2D. Отметим, что мы не может использовать PhaseSpacePlot для двумерных коэффициентов. Двумерный эквивалент функции ShowBasisPosition есть  ShowBasisPosition2D, которая показывает как пространство приближения   расщеплено в пространство деталей и пространство приближения более низкого уровня разрешения.

 

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

In[9]:= ShowBasisPosition2D[wtdata, AspectRatio -> Automatic]

 

            Используя функцию ListDensityPlot  строим графики вейвлет-коэффициентов. Отметим, что остаточный тренд нормализован отлично от деталей, чтобы сделать график более заметным.

In[10]:= PlotCoefficients2D[wtdata]

 

            Мы можем также построить график каждой матрицы отдельно. Это - остаточная или "усредненная" часть изображения.

In[11]:= ListDensityPlot[wtdata[[1]], Mesh -> False]

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

 

In[12]:= Show[GraphicsArray[ListDensityPlot[
#, Mesh -> False, DisplayFunction -> Identity]& /@
wtdata[[2]]]]

            Далее детали на следующем уровне.

In[13]:= Show[GraphicsArray[ListDensityPlot[
#, Mesh -> False, DisplayFunction -> Identity]& /@
wtdata[[3]]]]

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

 

            Если мы возьмём самый тонкий уровень , тогда самый грубый уровень . Это дает тренд в диадических точках.

In[14]:= temp1 = InverseWaveletTransform[
wtdata-MapAt[#-#&, wtdata, {1}], s4,
BoundaryCondition -> Fixed];

 

            Это - график плотности .

In[15]:= ListDensityPlot[temp1, Mesh -> False]

 

 

            Аналогично детали можно получить, оставляя соответствующие коэффициенты деталей и делая обратное преобразование. Здесь мы полагаем . Мы можем также получить, скажем   , оставляя только коэффициенты в положении  из wtdata и остальные - обнулить.

In[16]:= InverseWaveletTransform[
wtdata-MapAt[#-#&, wtdata, {2}], s4,
BoundaryCondition -> Fixed];

 

            Это - график плотности деталей на самом грубом уровне.

In[17]:= ListDensityPlot[%, Mesh -> False]

 

            Сумма тренда   изображена ранее и детали на том же самом уровне изображенные выше дают тренд на следующем более тонком уровне показанный здесь. Несмотря на то, что детали на уровнях  6 и 7 неудовлетворительны, воспроизведение разумно. Мы можем также получить прямо используя обратное преобразование на wtdata со всеми нулевыми коэффициентами деталей выше уровня 5.

In[18]:= ListDensityPlot[temp1+%%, Mesh -> False]

 

            Здесь мы делаем пакетное вейвлет-преобразование изображения.

In[19]:= wpdata = WaveletPacketTransform[data, s4, 3,
BoundaryCondition -> Fixed];

 

            По сравнению с вейвлет-преобразованием, пакетное вейвлет-преобразование раскладывает некоторые подпространства далее.

In[20]:= ShowBasisPosition2D[wpdata, AspectRatio -> 1]

 

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

In[21]:= PlotCoefficients2D[wpdata]

 

            Мы вычисляем кумулятивную энергию пакетных вейвлет-коэффициентов и просто вейвлет- коэффициентов.

In[22]:= (wtE=CumulativeEnergy[wtdata, 3500]/
(Flatten[wtdata].Flatten[wtdata]);
wpE=CumulativeEnergy[wpdata, 3500]/
(Flatten[wpdata].Flatten[wpdata]));

 

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

In[23]:= Show[{ListPlot[Transpose[{Table[256 256/i, {i, 1000, 3500}],
Sqrt[1-Take[wtE,{1000, 3500}]]}], PlotJoined -> True,
DisplayFunction -> Identity, PlotStyle -> Dashing[{0.02}],
PlotLabel -> "relative error vs. compression ratio"],
ListPlot[Transpose[{Table[256 256/i, {i, 1000, 3500}],
Sqrt[1-Take[wpE, {1000, 3500}]]}], PlotJoined -> True,
DisplayFunction -> Identity]},
DisplayFunction -> $DisplayFunction]

 

            Из вышеупомянутой кривой мы видим, что, если мы хотим сохранить погрешность ниже 5 %, степень сжатия должна быть ниже приблизительно 35. Мы хотим сохранять 1985 наибольших коэффициентов. Это соответствует 3 % полных коэффициентов и степени сжатия примерно 33.

In[24]:= InverseWaveletPacketTransform[
Compress[wpdata, 1985], s4,
BoundaryCondition -> Fixed];

 

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

In[25]:= ListDensityPlot[%, Mesh -> False]

            Это дает погрешности сжатия.

In[26]:= Sqrt[#.#&[Flatten[data-%%]]/
(Flatten[data].Flatten[data])]

Out[26]=

 

            Это - вычисленная погрешность, используя кумулятивную энергию. Несоответствие с вышеупомянутым результатом связано с граничными эффектами.

In[27]:= Sqrt[1-wpE[[1985]]]

Out[27]=

 

Вернуться в оглавление        Предыдущий параграф       Следующий параграф