Image processing & computer vision definitions

This notebook defines a comprehensive set of definitions useful for image processing and computer vision and will evolve as the course progresses. This is version 1.1. Subsequent updates will be tagged with different version numbers.

Version history

Version 1: Initial version posted to web (2/5/2004).

Version 1.1: Added image scaling functions (2/8/2004).

General definitions

Version number

$DefinitionsVersion  =  1.1;

Set default image file directory

\(\($ImagePath = "/Users/nechyba/public_html/eel6562/course_materials/images/mathematica";\)\)

Turn off annoying spelling warnings

Off[General::spell];
Off[General::spell1];

Load important packages

<< LinearAlgebra`MatrixManipulation`;
<< Statistics`MultinormalDistribution`;
<< Statistics`NormalDistribution`;
<< Statistics`DataManipulation`;
<< Graphics`Graphics`;
<< Graphics`Colors`;
<< Graphics`ImplicitPlot`;
<< Utilities`FilterOptions`;
<< Calculus`FourierTransform`;

Set default options for specific functions

SetOptions[Plot, Frame -> True];
SetOptions[ListPlot, Frame -> True,  PlotJoined -> True,  
        Axes -> {True, False}];
SetOptions[Fourier, 
        FourierParameters -> {1, -1}];
SetOptions[InverseFourier, 
        FourierParameters -> {1, -1}];
SetOptions[FourierTransform,  
        FourierParameters -> {0, -2 π}];
SetOptions[FourierTransform,  
        FourierParameters -> {0, -2 π}];
SetOptions[InverseFourierTransform,  
        FourierParameters -> {0, -2 π}];

Basic image processing

Basic image display and conversion

Matrix to Mathematica Image object

ImageObject[m_, size_:  1, max_:  255,  
          min_:  0] := 
Module[{sz, dimx, dimy, 
            cf}, 
If[size > 20, sz = 20, 
            sz = size]; 
{dimy, dimx} = 
            Take[Dimensions[m], {1, 2}]; 
cf = 
            If[MatrixQ[m],  GrayLevel, 
              RGBColor]; 
Graphics[
            Raster[m, {{0, 0}, {dimx, dimy}}, {min, max}, 
              ColorFunction -> cf], 
            ImageSize -> {Round[sz*dimx], Round[sz*dimy]}, 
            PlotRange -> {{0, dimx}, {0, dimy}},  
            AspectRatio -> Automatic]
];

Display Mathematica image object

ShowImage[m_,  size_:  1, max_:  255,  min_:  0,  
          opts___] := 
Show[
          ImageObject[m, size, max,  min],  opts];

Read image file to matrix

ReadImage[file_] := Import[file][[1, 1]];

Save matrix to image file

SaveImage[file_,  m_,  
          format_: "PPM"] := 
Export[file,  
          ImageObject[m],  format];

Color to grayscale

ColorToGrayScale[m_] :=  
        Map[((#[[1]] + #[[2]] + #[[3]])/3. ) &,  
          m,  {2}];

Grayscale to color (data conversion)

GrayScaleToColor[m_] := Map[{#, #, #} &,  m, {2}];

Image scaling

Scale one image channel to new dimensions ((Subscript[x, n]), (Subscript[y, n]))

ScaleOneChannel[
          mc_, {xn_,  yn_}] := 
Module[{xo, yo, f1, 
            f2}, 
{yo, xo} = 
            mc  //  Dimensions; 
f1 = 
            Partition[
              Flatten[Table[{y, x,  mc[[y, x]]}, {y, 1, yo}, {x, 1, 
                    xo}]], 3]; 
f2 = 
            Interpolation[f1, InterpolationOrder -> 1]; 
          Round[Table[
              f2[y, x], {y, 1, yo, (yo - 1)/(yn - 1)}, {x, 1, 
                xo, (xo - 1)/(xn - 
                      1)}]]
];

Scale image to new dimensions ((Subscript[x, n]), (Subscript[y, n]))

ScaleImage[
          m_, {xn_, yn_}] := 
Module[{r, g, 
            b}, 
If[
            MatrixQ[m], 
ScaleOneChannel[
              m, {xn, yn}], 
r = 
              ScaleOneChannel[
                Map[#[[1]] &, m, {2}], {xn, yn}]; 
g = 
              ScaleOneChannel[
                Map[#[[2]] &, m, {2}], {xn, yn}]; 
b = 
              ScaleOneChannel[
                Map[#[[3]] &, m, {2}], {xn, 
                  yn}]; 
Transpose  /@  ({r, g, 
                    b}  //  
                  Transpose)
]
];

Scale image by a factor α such that ((Subscript[x, n]), 
          (Subscript[y, n])) = α ((Subscript[x, old]), (Subscript[y, old]))

ScaleImage[m_, 
          fac_] := 
Module[{x, y}, 
If[
            MatrixQ[m], 
{y, x} = 
              Dimensions[m], 
{y, x} = 
                Drop[Dimensions[
                    m], -1];
]; 
ScaleImage[m, 
            Round[fac*{x, y}]]
];

Simple image creation

Create single color image

SingleColorImage[rgb_,  {dimx_, dimy_}] := 
Table[
          rgb, {dimy}, {dimx}];

Clip single pixel/channel image values

Clip[v_, max_:  255,  min_:  0] :=  
        If[Length[v] == 0, 
          Min[Max[v, min], max], Min[Max[#, min], max] &  /@  
            v];

Create Gaussian noise image

GaussianNoiseImage[{μ_,  σ_},  {dimx_, dimy_}, 
          max_:  255,  
          min_:  0] := 
Module[{m, 
            offset}, 
offset = 
            max/2  //  N; 
m  =  
            Table[Random[
                  NormalDistribution[μ, σ]], {dimx}, {dimy}] + 
              offset; 
Map[Clip[#, max] &,  
            m,  {2}]
];

Histogram modifications

Support function for histogram plotting

HistogramListPlot[data_, color_:  Black, 
          opts___] := 
Module[{g1}, 
g1 = Line[{{#[[1]], 0}, {#[[1]], #[[2]]}}] &  /@  data; 
          Show[Graphics[{color, g1}], opts]
];

Plot image histogram

ShowHistogram[m_,  max_:  255,  
          opts___] := 
Module[{mi, vals, 
            histogram}, 
mi = 
            Round[m]; 
vals = 
            Table[i, {i, 0, max}]; 
histogram = 
            CategoryCounts[Flatten[mi], 
              vals]; 
HistogramListPlot[{vals, 
                histogram}  //  Transpose,  Blue,  Frame -> True, 
            PlotRange -> All, opts]
];

Histogram equalization

HistogramEqualization[m_, 
          max_:  255] := 
Module[{mi, histogram, step1, 
            map}, 
mi = 
            Round[m]; 
histogram = 
            CategoryCounts[Flatten[mi], 
              Table[i, {i, 0, max}]]; 
step1 = 
            CumulativeSums[histogram]; 
map = 
            Round[step1/Last[step1]*max]; 
Map[
            map[[# + 1]] &, mi, {2}]
];

Histogram stretch

HistogramStretch[m_,  min_:  0,  
          max_:  255] := 
Module[{imin, 
            imax}, 
imin  =  Min[m];  
          imax  =  
            Max[m]; 
Map[((# - imin)/(imax - 
                          imin)*(max - min) + min) &, 
            m, {2}]
];

Histogram slide

HistogramSlide[m_,  offset_,  max_:  255] := 
Map[
          Clip[#, max] &, m + offset, {2}];

Pixel modifications (apply function to individual pixels)

Gray scale: \(TraditionalForm\`f : ℛ\^1 -> ℛ\^1\)such that \(TraditionalForm\`\(i\^*\) = f(i)\)

Color: \(TraditionalForm\`f : ℛ\^3 -> ℛ\^3\)such that \(TraditionalForm\`\((\(r\^*\), \(g\^*\), \(b\^*\))\) = f(r, g, b)\)

PixelModification[m_,  f_,  max_:  255] := 
Round[
          Map[Clip[#, max] &, Map[f,  m, {2}], {2}]];

2D convolution with masks

Mask operators

Smoothing

MeanFilter[n_] := Table[1, {n}, {n}]/(n*n)  //  N;

GaussianFilter[
          n_, σ_:  1] := 
Module[{nc, 
            gf1}, 
nc = (n - 1)/
              2; 
gf1 =  
            Table[Exp[-((x - nc)^2  +  (y - 
                                nc)^2)/(2  σ^2)], {x, 0, 
                  n - 1}, {y, 0, n - 1}]  //  
              N; 
gf1/(Plus  @@  
                Flatten[gf1]) 
];

Edges

PrewittRow[] := {{-1, -1, -1}, {0, 0, 0}, {1, 1, 
            1}};
PrewittColumn[] := PrewittRow[] // Transpose;

SobelRow[] := {{-1, -2, -1}, {0, 0, 0}, {1, 2, 1}}  //  
          N;
SobelColumn[] := SobelRow[]  //  Transpose;

Laplace1[] := {{0, -1, 0}, {-1, 4, -1}, {0, -1, 
            0}};
Laplace2[] := {{1, -2, 1}, {-2, 4, -2}, {1, -2, 
            1}};
Laplace3[] := {{-1, -1, -1}, {-1, 
            8, -1}, {-1, -1, -1}};

Texture example

(L = Level, E = Edge, S = Spot, R = Ripple) TexureMasks: (EL, RL, ES, RR, SL, EE, ER, SR)

L5T[] := {1, 4, 6, 4, 1};
E5T[] := {-1, -2, 0, 2, 1};
S5T[] := {-1, 0, 2, 0, -1};
R5T[] := {1, -4, 6, -4, 1};

    
L5E5[] := Outer[Times, L5T[], E5T[]];
E5L5[] := Outer[Times, E5T[], L5T[]];
EL[] := (L5E5[]  +  E5L5[])/2;

    
L5R5[] := Outer[Times, L5T[], R5T[]];
R5L5[] := Outer[Times, R5T[], L5T[]];
RL[] := (L5R5[]  +  R5L5[])/2;

    
E5S5[] := Outer[Times, E5T[], S5T[]];
S5E5[] := Outer[Times, S5T[], E5T[]];
ES[] := (E5S5[]  +  S5E5[])/2;

    
SS[] := Outer[Times, S5T[], S5T[]];
RR[] := Outer[Times, R5T[], R5T[]];

    
L5S5[] := Outer[Times, L5T[], S5T[]];
S5L5[] := Outer[Times, S5T[], L5T[]];
SL[] := (L5S5[]  +  S5L5[])/2;

    
EE[] := Outer[Times, E5T[], E5T[]];

    
E5R5[] := Outer[Times, E5T[], R5T[]];
R5E5[] := Outer[Times, R5T[], E5T[]];
ER[] := (E5R5[]  +  R5E5[])/2;

    
S5R5[] := Outer[Times, S5T[], R5T[]];
R5S5[] := Outer[Times, R5T[], S5T[]];
SR[] := (S5R5[]  +  R5S5[])/2;

2D convolution

Pad image with border rows so that 2D convolution does not shrink image

PadImage[
          m_] := 
Module[{mout}, 
mout = 
            Join[{First[m]}, m]; 
mout = 
            Join[mout, {Last[mout]}]; 
mout  =  
            Transpose[mout]; 
mout = 
            Join[{First[mout]}, mout]; 
mout = 
            Join[mout, {Last[mout]}]; 
mout  //  
            Transpose
];

2D convolution of a single channel

MaskOneChannel[mc_,  
          kernel_] := 
Module[{m2}, 
m2 = 
            Nest[PadImage, 
              mc, (First[Dimensions[kernel]] - 1)/
                2]; 
Round[
            Abs[ListConvolve[kernel, m2]]]
];

2D convolution of an image w/kernel (image dimensions unaffected for square masks with odd-numbered dimensions)

Mask[m_,  
          kernel_] := 
Module[{r, g, 
            b}, 
If[
            MatrixQ[m], 
MaskOneChannel[m, 
              kernel], 
r = 
              MaskOneChannel[Map[#[[1]] &, m, {2}], 
                kernel]; 
g = 
              MaskOneChannel[Map[#[[2]] &, m, {2}], 
                kernel]; 
b = 
              MaskOneChannel[Map[#[[3]] &, m, {2}], 
                kernel]; 
Transpose  /@  ({r, g, 
                    b}  //  
                  Transpose)
]
];

Median filtering

Filter one channel with (2n+1) x (2n+1) median filter

MedianFilterOneChannel[mc_,  
          n_] := 
Module[{ker, m1, 
            m2}, 
ker = 
            Table[0, {2*n + 1}, {2*n + 1}]; 
m1 = 
            Nest[PadImage, mc, n]; 
m2 = 
            ListConvolve[ker, m1, {-1, 1},  0,  Plus, 
              List]; 
Map[Median[Flatten[#]] &,  
            m2, {2}]
];

Filter image with (2n+1) x (2n+1) median filter

MedianFilter[m_, 
          n_] := 
Module[{r, g, b}, 
If[
            MatrixQ[m], 
MedianFilterOneChannel[m, 
              n], 
r = 
              MedianFilterOneChannel[Map[#[[1]] &, m, {2}], 
                n]; 
g = 
              MedianFilterOneChannel[Map[#[[2]] &, m, {2}], 
                n]; 
b = 
              MedianFilterOneChannel[Map[#[[3]] &, m, {2}], 
                n]; 
Transpose  /@  ({r, g, b}  //  
                  Transpose)
]
];

Image thresholding

Simple thresholding

ThresholdImage[m_, t_,  
          max_:  255] := 
Module[{m1}, 
m1 = 
            If[MatrixQ[m], m, ColorToGrayScale[m]]; 
Map[
            If[#  <  t,  0,  max] &,  
            m1,  {2}]
];

Gaussian thresholding (defined only for color images)

Mahalanobis distance from point to Gaussian distribution

Mahalanobis[
          x_, μ_, Σi_] := 
(x - μ) . Σi . (x - μ);

Gaussian threshold an image

GaussianThresholdImage[m_, t_, {μ_, Σ_}, 
          max_:  255] := 
Module[{Σi}, 
Σi = 
            Inverse[Σ]; 
Map[
            If[Mahalanobis[#, μ, Σi]  >  t,  0,  
                max] &,  m,  {2}]
];

Create Gaussian distance image

GaussianDistanceImage[
          m_, μ_, Σ_] := 
Module[{Σi}, 
Σi = 
            Inverse[Σ]; 
Map[
            Mahalanobis[#, μ, Σi] &, 
            m, {2}]
];

Hough transform (straight lines)

Hough transform helper function

ConvertEdgePixel[{xc_, yc_}, θt_,  {ρmax_, 
            dρ_}] := 
Module[{t1, 
            t2}, 
t1 = (xc Cos[#]  +  
                    yc Sin[#]) &  /@  θt; 
t2 = Ceiling[(# + ρmax)/dρ] &  /@  
              t1; 
Table[{t2[[i]], i}, {i, 1, 
              Length[t2]}]
];

From binary image, generate a sorted list of scaled line parameters

HoughEngine[mb_, 
          n_: {100, 100}] := 
Module[{ep,  nθ,  
            nρ,  θt, ρmax, dρ,  ydim, xdim, 
            ha1}, 
{nθ,  nρ} = 
            n; 
{ydim, xdim} = 
            Dimensions[mb]; 
ρmax = 
            Sqrt[xdim^2  +  ydim^2]  //  
              N; 
dρ  =  
            2. *ρmax/nρ; 
θt  =  
            Table[θ, {θ, 
                  0, π - π/nθ, π/nθ}]  //  
              N; 
ep = 
            Position[mb, 
              255]; 
ha1 = ConvertEdgePixel[#, θt, {ρmax, dρ}] &  /@  ep; 
Reverse[
            Sort[Frequencies[
                Join  @@  ha1]]]
];

From binary image, generate Hough accumulator matrix

HoughAccumulator[mb_, 
          n_: {100, 100}] := 
Module[{tmp, 
            ha}, 
tmp  =  
            HoughEngine[mb, n]; 
ha = 
            Table[0, {n[[2]]}, {n[[1]]}]; 
Fold[
            ReplacePart[#1, #2[[1]], #2[[2]]] &,  ha,  
            tmp]
];

From binary image, returns the (ρ,  θ) parameters for the top (Subscript[n, lines]), where ρ  =  
        x cos(θ)  +  y sin(θ)

StraightLines[mb_,  nlines_:  1,  
          n_: {100, 100}] := 
Module[{nθ,  
            nρ, ρmax, dρ,  ydim, 
            xdim}, 
{nθ,  nρ} = 
            n; 
{ydim, xdim} = 
            Dimensions[mb]; 
ρmax = 
            Sqrt[xdim^2  +  ydim^2]  //  
              N; 
dρ  =  
            2. *ρmax/
                nρ; 
{#[[1]] dρ  - ρmax - dρ/2,  (#[[2]] - 1) π/nθ} &  /@  
              Take[Last  /@  HoughEngine[mb, n], {1, nlines}]  // 
            N
];

Visualize detected straight lines over some image

ShowStraightLines[m_, sl_, 
        opts___] := 
Module[{xdim, ydim, 
          g1}, 
xdim = 
          Dimensions[m] //  Last; 
ydim = 
          Dimensions[m] // First; 
g1 = 
          ImplicitPlot[
            y Cos[sl[[2]]]  + x Sin[sl[[2]]]  ==  
              sl[[1]], {x, 0, xdim},  
            PlotStyle -> {Red, Thickness[ .01]},  
            PlotRange -> {{0, xdim}, {0, ydim}}, 
            DisplayFunction -> Identity]; 
g2 = 
          ShowImage[m, 1, 255, 0, 
            DisplayFunction -> Identity]; 
Show[g2, 
          g1, DisplayFunction -> $DisplayFunction, 
          opts]
]

Visualize Hough accumulator

ShowHoughAccumulator[m_,  ha_, 
          opts___] := 
Module[{ρmax}, 
ρmax  =  
            Sqrt[Dimensions[m][[1]]^2 + Dimensions[
                        m][[1]]^2]  //  
              N; 
ListDensityPlot[ha,  Mesh -> False, 
            ColorFunction -> Hue, 
            MeshRange -> {{0, π}, {-ρmax, ρmax}},  
            PlotRange -> All, opts]
];

Blob finding ("Connected Component Labeling")

Identify foreground neighbors of pixel (x, y) in image

NonZeroNeighbors[
        m_, {x_, 
          y_}] := 
DeleteCases[{m[[y - 1]][[x - 
                1]], m[[y - 1]][[x]], m[[y]][[x - 
                1]]}, 0]

First pass of blob finding algorithm

BlobPass1[
          mc_] := 
Module[{d1, m2, m3, label, nzn, 
            new}, 
d1 = 
            Dimensions[mc]; 
m2 = 
            PadRight[PadLeft[mc, d1 + {1, 1}], 
              d1 + {2, 2}]; 
label = 
            1; 
Do[
new = 
              If[m2[[y]][[x]] > 0, 
nzn = 
                  NonZeroNeighbors[m2, {x, y}]; 
If[
                  Length[nzn] > 0, 
Min[
                    nzn], 
label = label + 1], 
                0]; 
m2 = ReplacePart[m2, new, {y, x}], {y, 
              2, d1[[1]] + 1}, {x, 2, 
              d1[[2]] + 1}]; 
m3 = 
            TakeMatrix[
              m2, {2, 2}, {d1[[1]] + 1, 
                d1[[2]] + 1}]; 
Map[
            If[# == 0, #, # - 1] &, 
            m3, {2}]
];

Unify adjacencies of labels after first pass of blob finding algorithm: support functions

PartialAdjacency[
          l_] := (Union  /@  (Flatten  /@  
                Table[l[[First  /@  Position[l, i]]], {i, 1,  
                    Max[l]}]));

AdjacencyRules[l_] := 
        Table[l[[i]] -> Last[l], {i, 1, 
            Length[l] - 1}];

Second pass of blob finding algorithm

BlobPass2[
          mc_] := 
Module[{tmp1, tmp2, tmp3, 
            tmp4}, 
tmp1 = 
            ListConvolve[{{0, 0}, {0, 0}}, mc, {-1, 1},  0,  Plus, 
              List]; 
tmp2 = 
            Map[Sort[Union[DeleteCases[Flatten[#], 0]]] &, 
              tmp1, {2}]; 
tmp3 = 
            Flatten[Map[If[Length[#] < 2, Null, #] &, tmp2, {2}], 
              1]; 
tmp4 =  
            AdjacencyRules  /@  
              FixedPoint[PartialAdjacency, 
                Union[DeleteCases[tmp3, Null]]]; 
mc  /.  
            Flatten[Union[tmp4]]
];

Label all blobs with unique label in a binary image

BlobFind[mc_] := BlobPass2[BlobPass1[mc]];

Count the number of blobs in a labeled blob image (after BlobFind[ ])

BlobCount[mb_] :=  Length[Union[Flatten[mb]]] - 1;

Eliminate all but the largest blob from a labeled blob image (after BlobFind[ ])

KeepLargestBlob[
          mb_] := 
Module[{tmp1}, 
tmp1 = 
            Last  /@  Drop[
                Reverse[
                  Sort[Frequencies[DeleteCases[mb  //  Flatten, 0]]]], 
                1]; 
mb  /. ((# -> 
                      0) &  /@  
                tmp1)
];

Display a blob image with different random color assigned to each blob

BlobImage[
          mb_] := 
Module[{levels, 
            colors}, 
levels = 
            Drop[Sort[Union[Flatten[mb]]], 1]; 
colors =  
            Join[{0 -> {0, 0, 0}}, 
              Table[levels[[i]] -> {Random[Integer, {50, 255}], 
                    Random[Integer, {50, 255}], 
                    Random[Integer, {50, 255}]}, {i, 1, 
                  Length[levels]}]]; 
mb /. 
            colors];

Morphological operators for binary images

Dilation

Dilation[mc_, kernel_, 
          max_:  255] := 
Module[{m2, d1, 
            d2}, 
d1 = 
            Dimensions[mc]; 
d2 = 
            Dimensions[kernel]/2  //  Floor; 
m2 = 
            PadRight[PadLeft[mc, d1 + d2], 
              d1 + 2*d2]; 
ListConvolve[kernel, 
              m2, {-1, 1},  0,  BitAnd,  Max]  /.  
            1 -> max
];

Erosion

Erosion[mc_, kernel_,  
          max_:  255] := 
Module[{m2, d1, 
            d2}, 
d1 = 
            Dimensions[mc]; 
d2 = 
            Dimensions[kernel]/2  //  Floor; 
m2 = 
            PadRight[PadLeft[mc, d1 + d2], 
              d1 + 2*d2]; 
ListCorrelate[kernel, 
              m2, {1, -1},  0, BitOr[BitAnd[#1, #2], 1 - #1] &,  
              Min]  /.  1 -> max
];

Closing

Closing[mc_, kernel_, max_:  255] := 
Erosion[
          Dilation[mc, kernel, max], kernel, max];

Opening

Opening[mc_, kernel_, max_:  255] := 
Dilation[
          Erosion[mc, kernel, max], kernel, max];

Discrete Fourier Transform (DFT)

Visualization functions

Fix time variable for sampled signals

FixTime[signal_, 
          sampleRate_] := 
Module[{time, 
            data}, 
time = 
            Table[i/sampleRate, {i, 1, 
                Length[signal]}]; 
data  =  {time, 
                signal} // Transpose
];

Plot discrete-time signal as sequence of amplitude lines (large end point)

DiscreteListPlot[signal_, color_:  Black, 
          opts___] := 
Module[{data, g1, 
            g2}, 
If[
            Length[signal[[1]]]  ==  2, 
data = 
              signal, 
data = 
              FixTime[signal, 
                1]
]; 
g1 = Line[{{#[[1]], 0}, {#[[1]], #[[2]]}}] &  /@  
              data; 
g2 = {PointSize[ .015], 
                  Point[#]} &  /@  data; 
Show[
            Graphics[{color, g1, g2}], 
            opts]
];

Plot discrete-time signal as sequence of amplitude lines (small end point)

DiscreteListPlot2[signal_, color_:  Black, 
          opts___] := 
Module[{data, g1, 
            g2}, 
If[
            Length[signal[[1]]]  ==  2, 
data = 
              signal, 
data = 
              FixTime[signal, 
                1]
]; 
g1 = Line[{{#[[1]], 0}, {#[[1]], #[[2]]}}] &  /@  
              data; 
g2 = {PointSize[ .005], 
                  Point[#]} &  /@  data; 
Show[
            Graphics[{color, g1, g2}], 
            opts]
];

Fourier rearrange

FourierRearrange[fft] rearranges the result of an FFT (Mathematica function Fourier[ ]), to a easier to visualize the spectrum represented by 'fft'. The length of the list 'fft' is assumed to be even; if it is not, the last term in 'fft' is dropped.

FourierRearrange[fft_,  
          samplingRate_] := 
Module[{fft2, n, 
            n2}, 
If[EvenQ[Length[fft]],  
            fft2  =  fft,  
            fft2  =  Drop[fft, -1]];  
n = 
            Length[fft]; 
          n2 = n/2; 
{Table[ 
                samplingRate i/n, {i, -n2 + 1, n2}], 
              Take[fft2, {n2 + 2, n2*2}]~Join~Take[fft2, {1, n2 + 1}]} // 
            Transpose
];

Filtering definitions

Filter[n, sampleRate, {lowFreq, highFreq}] generates an ideal discrete filter for an FFT list of length 'n' (sampled at 'sampleRate'), with a bandpass from 'lowFreq' to 'highFreq' (in Hz).

Filter[n_, 
          samplingRate_, {lowF_, 
            highF_}] := 
Module[{coeff, l, 
            h}, 
coeff = 
            Table[0, {n}]; 
l = Floor[lowF*n/samplingRate]; 
          h = Ceiling[highF*n/samplingRate]; 
If[
            l  ==  0,  l  =  1;  
            coeff = ReplacePart[coeff, 1, 1]]; 
coeff = 
            ReplacePart[coeff, 1, Table[{i}, {i, l + 1, h + 1}]]; 
          coeff = ReplacePart[coeff, 1, 
              Table[{i}, {i, n + 1 - h, 
                  n + 1 - l}]]
];

Sampler (CT to DT converter)

Sampler[xc_, fs_, 
          range_: {0, 5}] := 
Table[{n, xc[n/fs]}, {n, 
            range[[1]]*fs, range[[2]]*fs}];

Discrete Fourier Transform

DFT[
          x_] := 
Module[{x2,  k}, 
If[
            Length[First[x]]  >  1,  x2  =  Last  /@  x,  
            x2  =  x]; 
{Table[
                k, {k, 0, Length[x2] - 1}], Fourier[x2]}  //  
            Transpose];

Discrete Wavelet Transform (DWT)

Generate filter coefficients

Generate Debauchies filter coefficients of order n (n = even)

PDB[n_, y_] := 
      Sum[Binomial[n - 1 + k, k] y^k, {k, 0, n - 1}]

DaubechiesFilter[n_, 
        prec_:  20] := 
Module[{poly, r, small, z, 
            l}, 
poly = 
            PDB[n/2, 1/2 - 1/(4  z) - z/4]; 
r = 
            z /. N[Solve[poly == 0, z], 
                prec]; 
small = 
            Select[r, 
              Abs[#] < 
                  1 &]; 
poly = (z + 1)^(n/
                    2) Times @@ (z - small); 
l = 
            CoefficientList[poly, z] // Reverse; 
l/
              Sqrt[l . l]  //  Chop
] /; 
        EvenQ[n]  &&  n ≥ 2

Generate Haar filter coefficients

HaarFilter[] := {1, 1}/Sqrt[2]  //  N;

Generate high-pass filter coefficients from low-pass filter coefficients

HighPassFilter[filter_] := 
        Table[(-1)^i, {i, 0, Length[filter] - 1}]*
          Reverse[filter];

Generate inverse filter coefficients from low-pass filter coefficients

InverseFilter[
          filter_] := 
Module[{n, lpf, 
              hpf}, 
n = 
              Length[filter]; 
lpf = 
              PadRight[filter, 2  n - 2]; 
hpf = 
              PadRight[HighPassFilter[filter], 
                2  n - 2]; 
(Partition[
                    Flatten[
                      Table[{RotateRight[lpf, i], RotateRight[hpf, i]}, {i, 
                          0, n - 2, 2}]], 2  n - 2]  //  
                  Transpose)[[n - 1]]
]  /;  
          EvenQ[Length[filter]]  &&  Length[filter] > 0;

Decimation functions

Decimation for DWT (1D)

Decimate1D[data_] := First  /@  Partition[data, 2];

Decimation for inverse DWT (1D)

InverseDecimate1D[data_] := 
        Last  /@  Partition[data, 2];

1D Discrete Wavelet Transform (DWT)

1D DWT at specified level

Wavelet1DLevel[data_, filter_, 
          level_] := 
Module[{lpf, hpf, nd, d1, 
              d2}, 
lpf = filter; 
hpf = 
              HighPassFilter[lpf]; 
nd = 
              Length[data]; 
If[
              level == 
                1, 
{ListCorrelate[lpf, 
                      data, {1, 1}]  //  Decimate1D, 
                  ListCorrelate[hpf, data, {1, 1}]  //  Decimate1D}  //  
                Flatten, 
d2 = 
                Drop[data, 
                  nd/(2^(level - 1))]; 
d1 = 
                Drop[data, -Length[d2]]; 
If[
                Length[d1] < 
                  Length[lpf], 
data, 
w1 = {ListCorrelate[lpf, d1, {1, 1}]  //  Decimate1D, 
                    ListCorrelate[hpf, d1, {1, 1}]  //  
                      Decimate1D}; 
Join[Flatten[w1], 
                  d2]
]
]
]  /;  EvenQ[Length[filter]];

1D DWT

Wavelet1D[data_, filter_, 
          levels_:  
            Complete] := 
Module[{l2}, 
If[levels  ==  Complete,  
            l2 = Floor[Log[2, Length[data]/Length[filter]]] + 1, , 
            l2 = levels]; 
Fold[
            Wavelet1DLevel[#1, filter, #2] &, data, 
            Table[i, {i, 1, l2}]]
];

Inverse 1D DWT at single level

InverseWavelet1DLevel[data_, filter_, 
          level_] := 
Module[{lpf, hpf, nd, d1, 
              d2}, 
lpf = 
              InverseFilter[filter]; 
hpf = 
              HighPassFilter[lpf]; 
nd = 
              Length[data]; 
d2 = 
              Drop[data, nd/(2^(level - 1))]; 
d1 = 
              Drop[data, -Length[d2]]; 
d1  =  
              Transpose[Partition[d1, Length[d1]/2]]  //  
                Flatten; 
If[
              Length[d1] < 
                Length[lpf], 
data, 
w1 = {ListCorrelate[lpf, d1, {-1, -1}]  //  InverseDecimate1D, 
                  ListCorrelate[hpf, d1, {-1, -1}]  //  
                    InverseDecimate1D}; 
Join[
                Flatten[Transpose[w1]], 
                d2]
]
]  /;  
          EvenQ[Length[filter]];

Inverse 1D DWT

InverseWavelet1D[data_, filter_, 
          levels_:  
            Complete] := 
Module[{l2}, 
If[levels  ==  Complete,  
            l2 = Floor[Log[2, Length[data]/Length[filter]]] + 1, , 
            l2 = levels]; 
Fold[
            InverseWavelet1DLevel[#1, filter, #2] &, data, 
            Table[i, {i, l2, 1, -1}]]
];

2D Discrete Wavelet Transform (DWT)

2D DWT for one level (operates on whole input matrix)

WaveletOneLevel2D[m_, 
          filter_] := 
Module[{w1, 
              w2}, 
w1 = Wavelet1DLevel[#, filter, 
                    1] &  /@  
                m; 
w2 = Wavelet1DLevel[#, filter, 
                    1] &  /@  Transpose[w1]; 
w2  //  
              Transpose
] /;  
          EvenQ[Length[filter]];

2D DWT at specified level

This function assumes a square matrix with dimensions n ×n, where n =  (2^i),  i = {2,  3,  4,   ... }

Wavelet2DLevel[m_, filter_,  
          level_] := 
Module[{dim, n, ll, hl, lh, hh, 
              llw}, 
If[level == 1, 
              WaveletOneLevel2D[m, filter], 
dim = 
                Length[m]; 
n = 
                dim/(2^(level - 1)); 
ll = 
                TakeMatrix[m, {1, 1}, {n, n}]; 
lh = 
                TakeMatrix[m, {1, n + 1}, {n, dim}]; 
hl = 
                TakeMatrix[m, {n + 1, 1}, {dim, n}]; 
hh = 
                TakeMatrix[
                  m, {n + 1, n + 1}, {dim, dim}]; 
If[
                Length[ll] < 
                  Length[filter], 
m, 
llw = 
                  WaveletOneLevel2D[ll, 
                    filter]; 
BlockMatrix[{{llw, lh}, {hl, 
                      hh}}]
]
]
]  /;  EvenQ[Length[filter]];

2D DWT (uncompiled, slow)

Wavelet2D[data_, filter_, levels_] := 
Fold[
          Wavelet2DLevel[#1, filter, #2] &, data, 
          Table[i, {i, 1, levels}]];

2D DWT (compiled, fast)

WaveletTransform2D  =  
        Compile[{{m, _Real, 2}, {filter, _Real, 1}, {level, _Integer}}, 
          Wavelet2D[m, filter, level]];

Inverse 2D DWT for one level (operates on whole input matrix)

InverseWaveletOneLevel2D[m_, 
          filter_] := 
Module[{w1, 
              w2}, 
w1 = InverseWavelet1DLevel[#, filter, 
                    1] &  /@  
                Transpose[
                  m]; 
w2 = InverseWavelet1DLevel[#, 
                    filter, 1] &  /@  
                Transpose[w1]
] /;  
          EvenQ[Length[filter]];

Inverse 2D DWT at specified level

This function assumes a square matrix with dimensions n ×n, where n =  (2^i),  i = {2,  3,  4,   ... }

InverseWavelet2DLevel[m_, filter_,  
          level_] := 
Module[{dim, n, ll, hl, lh, hh, 
              llw}, 
If[level == 1, 
              InverseWaveletOneLevel2D[m, filter], 
dim = 
                Length[m]; 
n = 
                dim/(2^(level - 1)); 
ll = 
                TakeMatrix[m, {1, 1}, {n, n}]; 
lh = 
                TakeMatrix[m, {1, n + 1}, {n, dim}]; 
hl = 
                TakeMatrix[m, {n + 1, 1}, {dim, n}]; 
hh = 
                TakeMatrix[
                  m, {n + 1, n + 1}, {dim, dim}]; 
If[
                Length[ll] < 
                  Length[filter], 
m, 
llw = 
                  InverseWaveletOneLevel2D[ll, 
                    filter]; 
BlockMatrix[{{llw, lh}, {hl, 
                      hh}}]
]
]
]  /;  EvenQ[Length[filter]];

Inverse 2D DWT (uncompiled, slow)

InverseWavelet2D[data_, filter_, levels_] := 
Fold[
          InverseWavelet2DLevel[#1, filter, #2] &, data, 
          Table[i, {i, levels, 1, -1}]];

Inverse 2D DWT (compiled, fast)

InverseWaveletTransform2D  =  
        Compile[{{m, _Real, 2}, {filter, _Real, 1}, {level, _Integer}}, 
          InverseWavelet2D[m, filter, level]];

Visualization functions

Extract bands at specified level from 2D DWT

This function assumes a square matrix with dimensions n ×n, where n =  (2^i),  i = {2,  3,  4,   ... }

ExtractBands[mw_, 
          level_] := 
Module[{n1, n2, ll, lh, hl, 
            hh}, 
n2 = 
            Length[mw]/(2^(level - 1)); 
n1 = 
            n2/2; 
ll = 
            TakeMatrix[mw, {1, 1}, {n1, n1}]; 
lh = 
            TakeMatrix[mw, {1, n1 + 1}, {n1, n2}]; 
hl = 
            TakeMatrix[mw, {n1 + 1, 1}, {n2, n1}]; 
hh = 
            TakeMatrix[
              mw, {n1 + 1, n1 + 1}, {n2, n2}]; 
{{ll, 
              lh}, {hl, hh}}
];

Helper functions (visualize signed wavelet coefficients, or magnitudes only)

Magnitude[a_] := Abs[a];
Signed[a_] := a;

Construct visualization for one level of 2D DWT

VisualizeOneLevel[mw_, level_, g_:  Null,  f_:  Signed,  
          opts___] := 
Module[{bands, g1, g2, 
            opts2}, 
bands = 
            ExtractBands[mw, level]; 
opts2 = 
            FilterOptions[{PlotRange}, opts]; 
g1 = 
            Map[ListDensityPlot[f[#], Mesh -> False, 
                  DisplayFunction -> Identity, Frame -> False, 
                  opts2] &, bands, {2}]; 
g2 = 
            If[g == Null, g1, , 
              ReplacePart[g1, g, {1, 1}]]; 
opts2 = 
            FilterOptions[{GraphicsSpacing}, opts]; 
Show[
            GraphicsArray[g2], opts2, GraphicsSpacing -> -0.02,  
            DisplayFunction -> 
              Identity]
];

Display 2D DWT (w/specified number of levels)

The output of this function is controlled by three main things: (1) levels = number of levels in 2D DWT; (2) f = {Signed, Magnitude} - display signed wavelet coefficients or magnitudes of coefficients only; (3) the PlotRange -> Automatic option (default) gives more detail, while the PlotRange -> All option adjusts output range to not clip coefficients (less detail).

ShowWaveletTransform2D[mw_, levels_, f_:  Signed,  
          opts___] := 
Module[{gl, dsp, 
            opts2}, 
gl = 
            Drop[FoldList[VisualizeOneLevel[mw, #2, #1, f, opts] &, Null, 
                Table[i, {i, levels, 1, -1}]], 
              1]; 
dsp = 
            Table[DisplayFunction -> Identity, {Length[gl] - 1}]~
              Join~{DisplayFunction -> $DisplayFunction}; 
opts2 = 
            FilterOptions[{ImageSize}, 
              opts]; 
Show[#[[1]], opts2, #[[2]], 
                ImageSize -> 500] &  /@  
            Transpose[{gl, dsp}]
];

Visualize 2D DWT levels

VisualizePyramidOneLevel[ndim_, level_, g_:  Null, 
          opts___] := 
Module[{m, v, bands, g1, g2, 
            opts2}, 
v = 
            level/(Log[2, ndim]  - 1)  +   .6; 
If[
            v > 1, v = v - 1]; 
m = 
            Table[v, {ndim}, {ndim}]; 
bands = 
            ExtractBands[m, level]; 
opts2 = 
            FilterOptions[{PlotRange}, opts]; 
g1 = 
            Map[ListDensityPlot[#, Mesh -> False, 
                  DisplayFunction -> Identity, Frame -> False, 
                  PlotRange -> {0, 1}, ColorFunction -> Hue] &, 
              bands, {2}]; 
g2 = 
            If[g == Null, g1, , 
              ReplacePart[g1, g, {1, 1}]]; 
opts2 = 
            FilterOptions[{GraphicsSpacing}, opts]; 
Show[
            GraphicsArray[g2], opts2, GraphicsSpacing -> 0.05,  
            DisplayFunction -> 
              Identity]
];

ShowPyramid2D[levels_,  
          opts___] := 
Module[{gl, dsp, 
            opts2}, 
gl = 
            Drop[FoldList[
                VisualizePyramidOneLevel[2^(levels + 1), #2, #1, opts] &, 
                Null, Table[i, {i, levels, 1, -1}]], 
              1]; 
dsp = 
            Table[DisplayFunction -> Identity, {Length[gl] - 1}]~
              Join~{DisplayFunction -> $DisplayFunction}; 
opts2 = 
            FilterOptions[{ImageSize}, 
              opts]; 
Show[#[[1]], opts2, #[[2]], 
                ImageSize -> 500] &  /@  
            Transpose[{gl, dsp}]
];

Visualize 1D DWT bands

ConstructOneBand1D[{{r1_, r2_}, h_}, 
          colors_] := 
{colors[[h]],  
          Rectangle[{r1, 0}, {r2, 1}]};

VisualizeBands1D[level_,  
          opts___] := 
Module[{ht, rt, m, 
            cols}, 
cols = {Green, Red, Blue, Magenta, 
              Cyan, Orange, Yellow, Pink, Purple, 
              Black}; 
rt  =  
            Partition[{0}~Join~Table[2^(5 - i), {i,  level, 0, -1}], 
              2, 1]; 
ht  = {2 level}~Join~
              Reverse[Table[
                  2  i + 1,  {i, 0,  
                    level - 1}]]; 
m = {rt, ht}  //  
              Transpose; 
Show[
            Graphics[ConstructOneBand1D[#, cols] &  /@  m],  
            AspectRatio -> Automatic, 
            opts]
];