HW6, computer problem
EECS 203A, UCI, Fall 2004

by Nasser Abbasi

FIrst start by reading my Mathematica pakcage. This contains Mathematica functions I have written for basic image processing and filtering

Clear["Global`*"];
<< nma.m
<< ImageProcessing`

\(\(\(
\)\(g2D[x_, 
        y_, σ_] := \(1\/\(\(\@\(2\ π\)\) σ\)\) E\^\(-\(\((x\^2 + y\^2)\)\/\(2\ σ\^2\)\)\)

    Block[{$DisplayFunction = 
          Identity}, 

p1 = 
        Plot3D[g2D[x, y, 40], {x, \(-240\), 240}, {y, \(-320\), 320}, 
          PlotPoints -> 40, PlotRange -> All, 
          PlotLabel -> "standard deviation=40"]; 

p2 = 
        Plot3D[g2D[x, y, 5], {x, \(-240\), 240}, {y, \(-320\), 320}, 
          PlotPoints -> 40, PlotRange -> All, 
          PlotLabel -> "standard deviation=5"];
]

    Show[GraphicsArray[{p1, 
          p2}]]



    \)\)\)

Now, sample a 2D continuos gaussian to obtain a discrte version of Gaussian 2D, to use as a filtering window to convolve the image with

nma`cd
Directory[]

nRow = 480;  nCol = 640;
f = BinaryReadList["triangle.raw"];
f =  Partition[f, nCol];
Print["Image has dimensions ",  
        Dimensions[f]];
fsaved = Chop[N[f]];
Print["Max[f]=", Max[f]];
Print["Min[f]=", Min[f]];

Now create the 2D gaussian filter, have it centers at the center of the image

σ = 5; 
 nRowG = 480;  nColG = 640;
h = 
        Table[g2D[x, y, σ], {x, - (nRowG/2), 
            (nRowG/2) - 1}, {y, - (nColG/2), 
            (nColG/2) - 1}];
hsaved = Chop[N[h]];
Print["Min value of gaussian, std=5, is ", 
        Min[h]];
Print["Max value of gaussian, std=5, is ", 
        Max[h]];

Now generate 'g', the degraded image

h = hsaved;

    f = fsaved;


    G = DFT[h] DFT[f];

    G = nma`centerImage[G];

    g = 
        InverseDiscreteFourierTransform[G];
  

  

    Max[Re[g]]

    Min[Re[g]]

    t = Re[g];

    t = (255/Max[t])*t;


    Export["triangleDegraded.raw", 
        Round[N[t]], "Byte"];


    ArrayPlot[Round[Re[g]]]

\(\(Export["triangleDegraded.raw", 
        t, "Byte"];\)

 (*file = 
        OpenWrite["triangleDegraded.raw", 
          BinaryFormat -> True]; 
BinaryWrite[file, 
        Flatten[t], "Character8"]; 
Close["triangleDegraded.raw"];*) 
\)
g = nma`centerImage[Re[g]];
G = DFT[g];

Print["degraded image dimensions=",  
        Dimensions[g]];

Now center the images, and generate their DFT

f = nma`centerImage[f];
h = 
        nma`centerImage[h];

F = DFT[f];
H = DFT[h];

\(Block[{$DisplayFunction = Identity}, \ 
p0 = 
        ArrayPlot[f, 
          PlotLabel -> "f, input image 480x640"]; 
p1 = 
        ArrayPlot[h, 
          PlotLabel -> "h, 32x32 gaussian"]; 
p2 = 
        ArrayPlot[g, 
          PlotLabel -> "g, 480x640 degraded image"]; 

p01 = 
        ArrayPlot[2\ Log[Abs[F]], 
          PlotLabel -> "2 Log[Abs[F]"]; 
p11 = 
        ArrayPlot[2\ Log[Abs[H]], 
          PlotLabel -> "2 Log[Abs[H]"]; 
p21 = 
        ArrayPlot[2\ Log[Abs[G]], 
          PlotLabel -> "2 Log[Abs[G]"]; 

p02 = 
        ArrayPlot[Arg[F], 
          PlotLabel -> "Phase F"]; 
p12 = 
        ArrayPlot[Arg[H], 
          PlotLabel -> "Phase H"]; 
p22 = 
        ArrayPlot[Arg[G], 
          PlotLabel -> "Phase G"];

]\)
Show[GraphicsArray[{p0, p1, p2}]];
Show[GraphicsArray[{p01, p11, p21}]];
Show[
        GraphicsArray[{p02, p12, 
            p22}]];

\(\(\(
\)\( (*Show[
        GraphicsArray[\(ListPlot3D[#1, PlotRange -> All, 
                DisplayFunction -> Identity] &\) /@ {Abs[F1], 
              Arg[F1]}]]; 
Show[
        GraphicsArray[\(ListPlot3D[#1, PlotRange -> All, 
                DisplayFunction -> Identity] &\) /@ {Abs[H1], 
              Arg[H1]}]]; 
Show[
        GraphicsArray[\(ListPlot3D[#1, PlotRange -> All, 
                DisplayFunction -> Identity] &\) /@ {Abs[G1], 
              Arg[G1]}]];*) \)\(
\)\(
\)\(\(Print["Size of F=", Dimensions[F]];\)

    \(Print["Max element in Re[F]=", 
        Max[Re[F]], \ "  Min Element in Re[F]=", 
        Min[Re[F]]];\)

    \(Print["Max element in Im[F]=", 
        Max[Im[F]], \ "  Min Element in Im[F]=", 
        Min[Im[F]]];\)

    \(Print["Max element in Abs[F]=", 
        Max[Abs[F]], \ "  Min Element in Abs[F]=", 
        Min[Abs[F]]];\)


    \(Print["Size of H=", Dimensions[H]];\)

    \(Print["Max element in Re[H]=", 
        Max[Re[H]], \ "  Min Element in Re[H]=", 
        Min[Re[H]]];\)

    \(Print["Max element in Im[H]=", 
        Max[Im[H]], \ "  Min Element in Im[H]=", 
        Min[Im[H]]];\)

    \(Print["Max element in Abs[H]=", 
        Max[Abs[H]], \ "  Min Element in Abs[H]=", 
        Min[Abs[H]]];\)


    \(Print["Size of G=", Dimensions[G]];\)

    \(Print["Max element in Re[G]=", 
        Max[Re[G]], \ "  Min Element in Re[G]=", 
        Min[Re[G]]];\)

    \(Print["Max element in Im[G]=", 
        Max[Im[G]], \ "  Min Element in Im[G]=", 
        Min[Im[G]]];\)

    \(Print["Max element in Abs[G]=", 
        Max[Abs[G]], \ "  Min Element in Abs[G]=", 
        Min[Abs[G]]];\)

    \)\)\)

Now generate restored image

{nRow, nCol} = Dimensions[G];
{hRow, hCol} = 
        Dimensions[
          H];

  
GG =  G/ H;
Print["Result of G/H has dimensions = ", 
      Dimensions[GG]]

ArrayPlot[10 Log[Abs[GG]], 
        PlotLabel -> "Log[Abs[RestoredImage]]"];
ArrayPlot[Arg[GG], 
        PlotLabel -> "Phase RestoredImage"];
m = nma`centerImage[GG];

ff = 
        InverseDiscreteFourierTransform[m];

    
ff = 
        nma`centerImage[Re[ff]];

    
Close["triangleRestore.raw"];
file = 
        OpenWrite["triangleRestore.raw", 
          BinaryFormat -> True];
BinaryWrite[file,Floor[ff],"Character8"];
Close["triangleRestore.raw"];

ArrayPlot[ff, PlotLabel -> "restored Image"];

\(Block[{$DisplayFunction = Identity}, \ 
p0 = 
        ArrayPlot[f, 
          PlotLabel -> "f, input image 480x640"]; 
p1 = 
        ArrayPlot[g, 
          PlotLabel -> "g, degraded image"]; 
p2 = 
        ArrayPlot[ff, 
          PlotLabel -> "ff, restored  image"];

]\)
Show[
        GraphicsArray[{p0, p1, p2}]];