# SFTW453 Digital Image Processing

AI and Robotics

Nov 6, 2013 (4 years and 6 months ago)

105 views

SFTW453

Digital Image Processing

Assignment 3

Discrete Fourier Transform

&
High
pass
Filtering

D
-
A
7
-
2819
-
0

Cao Han

Abstract

Part I: Discrete Fourier Transform

Implement a two
-
dimensional discrete Fourier Transform with following capabilities:

a
)

Multiply the input image by (
-
1)
x+y

to center the transform for filtering.

b)

Compute the Fourier transforms.

c)

Compute the inverse Fourier transform.

d)

Multiply the result by (
-
1)

x+y

and take the real part.

e)

Compute the Fourier spectrum.

The image

size may be limited to integer powers of 2. You may need to zoom or
shrink an image to the proper size by using the program you developed in Assignment
I.

Part II: High
pass Filtering

a)

Implement the Gaussian high
pass filter in Eq. (4.3
-
7). You must be a
ble to specify
the size,

M x N, of the resulting 2D function. In addition, you must be able to specify
where the 2D

location of the center of the Gaussian function.

b)

Download attached image Fig. 4.11(a) [this image is the same as Fig. 4.18(a)] and

high
pa
ss filter it to obtain Fig. 4.18(c).

Technical discussion

Part A

In this part
, we need

a program to compute the Discrete Fourier Transform (DFT)
,

Inv
erse Discrete Fourier Transform

and Fourier spectrum
.

We need to multiply the
input image by (
-
1)
x+y

to ce
nter the transform for filtering

before transform
.

T
he
method of
Discrete Fourier Transform

T
he
method of Inverse
Discrete Fourier Transform

T
he
method of
Fourier spectrum

Since the method include complex. So, we need to include the complex st
andard
library of C++ (#include <complex>). It can easy to take the real part and
imaginary

part

by calling function real() and I
mag().

I
t need to using the namespace std.

Part B

In the second part,
we need to write a program to impleme
nt the Gaussian
high
pass
filter.

The Gaussian high
pass filter

Basic of filter in the frequency domain

Filtering in the frequency domain is straightforward. It consists of the following steps:

1.

Multiply the input by (
-
1)
x+y

to center the transform.

2.

Compute F(u,v), t
he DFT of the image from (1).

3.

Multiply F(u,v) by a filter function H(u,v).

4.

Compute the inverse DFT of the result in(3).

5.

Obtain the real part of the result in(3).

6.

Multiply the result in (5) by (
-
1)
x+y
.

Discussion

of results:

W
e perform some actions: Fou
rier Transformation, inverse Fourier
Transformation, and therefore we get Fourier spectrum.

I
n part A,
W
e can use the spectrum to
simulate

the image

s change rate of
frequency.
A
nd also, we can use it to perform some special actions to optimize the
imag
es.

In part B, we use Gaussian filter to make the im
age smooth, looking much better
than the source image.

R
esult
:

Part A

Original Image

A
fter
Multiply the input by (
-
1)
x+y

After Fourier
Transformation

Inverse DFT and multiply by (
-
1)
x+y

Fourier spectrum

Part B

Original

Image

15

Program listings

#include "CImg.h"

#include <cmath>

#include <string.h>

#include <complex>

#define PI 3.14159

using namespace std;

using namespace
cimg_library;

complex<double> F[100][100];

complex<double> f[100][100];

complex<double> H[100][100];

// (b) Compute the two
-
dimensional discrete Fourier transform.

void FourierTransform(CImg<double> Image){

printf("Apply Fourier Transform
\
n");

for( int

u = 0; u <= Image.dimx()
-

1; u++){

for ( int v = 0; v <= Image.dimy()
-

1; v++){

complex<double> temp(0.0, 0.0);

for ( int x = 0; x <= Image.dimx()
-

1; x++){

for ( int y = 0; y<= Image.dimy()
-

1; y++){

temp += Image(x,y) * exp(complex
<double>( 0.0, (
-
2 * PI *
((double)u*x/Image.dimx() + (double)v*y/Image.dimy()))));

}

}

F[u][v] = temp / sqrt(Image.dimx() * Image.dimy());

}

}

for ( int i = 0; i < Image.dimx(); i++){

for ( int j = 0; j < Image.dimy(); j++)

Image(i,j)

= real(F[i][j]);

}

CImgDisplay MultiplyDisplay(Image, "Apply Fourier Transform");

system("PAUSE");

}

//(c) Compute the inverse two
-
dimensional discrete Fourier transform.

CImg<double> InverseFourierTransform(CImg<double> Image){

CImg<double> outputI
mage(Image.dimx(),Image.dimy(),1,1,0);

printf("Apply Inverse Fourier Transform
\
n");

int u,v;

for( u = 0; u <= Image.dimx()
-

1; u++){

for ( v = 0; v <= Image.dimy()
-

1; v++){

complex<double> temp(0.0, 0.0);

for ( int x = 0; x <= Image.dimx()

-

1; x++){

for ( int y = 0; y<= Image.dimy()
-

1; y++){

temp += F[x][y] * exp(complex<double>( 0.0, (2 * PI *
((double)u*x/Image.dimx() + (double)v*y/Image.dimy()))));

}

}

f[u][v] = temp / sqrt(Image.dimx() * Image.dimy());

}

}

for ( int x = 0; x <= Image.dimx()
-

1; x++){

for ( int y = 0; y<= Image.dimy()
-

1; y++){

outputImage(x,y)=real(f[x][y]);

}

}

return outputImage;

system("PAUSE");

}

//(a) Multiply the input image by (
-
1)x+y to center the transform
for filtering.

CImg<double> Multiply(CImg<double> Image){

printf("Multiply (
-
1)^(x+y) to center the transform for filtering
\
n");

CImg<double> MultiplyImage(Image.dimx(),Image.dimy(),1,1,0);

for ( int i = 0; i < Image.dimx(); i++){

for ( int j = 0; j
< Image.dimy(); j++)

if ( (i + j)%2 != 0)

MultiplyImage(i,j) =
-
Image(i,j);

else

MultiplyImage(i,j) = Image(i,j);

}

CImgDisplay MultiplyDisplay(MultiplyImage, "After (
-
1)^(x+y)");

system("PAUSE");

return MultiplyImage;

}

//(e) Compute
the Fourier spectrum.

void FourierSpectrum(CImg<double> Image){

printf("The Fourier spectrum
\
n");

CImg<double> FourierSpectrumImage(Image.dimx(),Image.dimy(),1,1,0);

for ( int i = 0; i < Image.dimx() ; i++){

for ( int j = 0; j < Image.dimy() ; j++){

FourierSpectrumImage(i,j) = log(1 + sqrt(pow(real(F[i][j]), 2) +
pow(imag(F[i][j]), 2)));

}

}

CImgDisplay FourierSpectrumDisplay(FourierSpectrumImage,
"FourierSpectrum");

system("PAUSE");

}

//partb:(a)Write a C++ function for the transfer function
of the Gaussian highpass
filter in Eq. (4.4
-
4).

void Gaussian(CImg<double> Image){

printf("Construct Gaussian highpass filter
\
n");

double D0, D;

scanf("%lf", &D0);

for ( int u = 0; u < Image.dimx() ; u++){

for ( int v

= 0; v < Image.dimy() ; v++){

D = sqrt(pow((double)u
-

((double)Image.dimx()/2),2) +
pow((double)v
-

((double)Image.dimy()/2),2));

H[u][v] =(complex<double>)(1,0.0)
-

exp(complex<double>(0.0,
-
(double)((double)pow(D, 2)/ (double)(2* pow(D0,2)))));

}

}

}

void Multiply_HF(CImg<double> Image){

printf("Multiply H(u,v) with F(u,v)
\
n");

for ( int u = 0; u < Image.dimx() ; u++){

for ( int v = 0; v < Image.dimy() ; v++)

F[u][v] *= H[u][v];

}

}

void main(){

char Pic[20];

ge:");

scanf("%s", Pic);

CImg<double> Image(Pic),result;

CImg<double> Result;

complex<double> Transform[50][50];

CImgDisplay main_disp(Image,"Initial Image");

printf("Assignment III highpass Filtering
\
n");

//(a) Multiply the input image by (
-
1)x+y to

center the transform for filtering.

Result = Multiply(Image);

//(b) Compute the two
-
dimensional discrete Fourier transform.

FourierTransform(Result);

//(c) Compute the inverse two
-
dimensional discrete Fourier transform.

result=InverseFourierTrans
form(Image);

//(d) Multiply the result by (
-
1)x+y and take the real part.

Multiply(result);

//(e) Compute the Fourier spectrum.

FourierSpectrum(Image);

//partb:(a) Write a C++ function for the transfer function of the Gaussian highpass
filter in
Eq. (4.4
-
4).

Gaussian(Result);

Multiply_HF(Result);

//

InverseFourierTransform(Result);

for ( int i = 0; i < Image.dimx(); i++){

for ( int j = 0; j < Image.dimy(); j++)

Result(i,j) = real(f[i][j]);

}

Result = Multiply(Result);

FourierSpectrum(R
esult);

}