%Image Perception

%Chapter 3

%Image Perception

% Line spread function

clear

% generating line spread function

B = [zeros(1,2) [0:-0.02:-0.12] [-0.14:0.18:1.0]] ;

C = [B(1:length(B)-1) fliplr(B) ]

C = C/sum(C) ;

len1 = length(C)-1 ; % L is odd

plot(C)

y = [1:200] ;

% generating step function

L = [zeros(1,100) 100*ones(1,115)] ;

x = conv(L,C) ;

x = x(len1/2+1:len1/2+length(L)-15) ;

length(x)

plot(y,L(1:200),y,x) ;

print -dtiff plot.tiff

% Mach Band

clear

% generating line spread function

B = [zeros(1,2) [0:-0.03:-0.18] [-0.20:0.20:1.0]] ;

C = [B(1:length(B)-1) fliplr(B) ]

C = C/sum(C) ;

len1 = length(C)-1 ; % L is odd

plot(C)

y = [1:600] ;

% generating step function

p = ones(1,100) ;

L = [zeros(1,100) 10*p 20*p 30*p 40*p 50*p 50*p] ;

x = conv(L,C) ;

x = x(len1/2+1:len1/2+length(L)-100) ;

length(x)

plot(y,L(1:600),y,x) ;

%Plotting Spatial Frequency Response of HVS

% MATLAB script to generate frequency response of the HVS

r =[0:0.1:30];

r1 = 0.2 + 0.45*r ;

r2 = exp(-0.18*r) ;

H = r1.*r2 ;

plot(r,H)

xlabel(‘Spatial frequency in cycles/degree’)

ylabel(‘Normalized spatial frequency response’),grid

print -dtiff plot.tiff

r =[0:1:50];

r1 = 0.2 + 0.45*r ;

r2 = exp(-0.18*r) ;

H1D = r1.*r2 ;

H2D = zeros(61,61);

for p= -30:30

for q = -30:30

r= round(sqrt(p*p+q*q)) ;

H2D(p+31,q+31) = H1D(r+1) ;

end

mesh([-30:30], [-30:30],H2D) ;

maxf =30 ;

r =[0:1:50];

r1 = 0.2 + 0.45*r ;

r2 = exp(-0.18*r) ;

H1D = r1.*r2 ;

H2D = zeros(2*maxf+1,2*maxf+1);

for p= -maxf:maxf

for q = -maxf:maxf

r= round(sqrt(p*p+q*q)) ;

H2D(p+maxf+1,q+maxf+1) = H1D(r+1) ;

end

mesh([-maxf:maxf], [-maxf:maxf],H2D) ;

xlabel(‘Horizontal frequency in cycles/deg’)

ylabel(‘Vertical frequency in cycles/deg’)

zlabel(‘Spatial frequency response’), grid

print -dtiff plot.tiff

% Additive Color Mixing

% MATLAB script to demonstrate effect of color addition

clear

r = zeros(256,256) ; % Size of the image

g = zeros(256,256) ;

b = zeros(256,256) ;

x1=92; y1=80; % centre of the red circle

x2=92; y2=175;

x3=175; y3=128;

radius = 70 ;

for x= 0:255

for y = 0:255

p = x+1 ; q = y+1 ;

d1 = sqrt((x1-x)^2 + (y1-y)^2) ; % distance from the center of the red circle

d2 = sqrt((x2-x)^2 + (y2-y)^2) ; % distance from the center of the green circle

d3 = sqrt((x3-x)^2 + (y3-y)^2) ; % distance from the center of the blue circle

r(p,q) = d1 < radius ; % red circle

g(p,q) = d2 < radius ; % green circle

b(p,q) = d3 < radius; % green circle

if ((d1>radius) & (d2>radius) & (d3>radius))

% if ((d1>=radius) & (d2>=radius) & (d3>=radius))

r(p,q) = 1 ; g(p,q)=1 ; b(p,q)=1 ; % white circle

end

im = zeros(256,256,3) ;

im(:,:,1) = double(r) ;

im(:,:,2) = double(g) ;

im(:,:,3) = double(b) ;

imshow(im)

imwrite(grating1,’F:\grat1.tif’,’tif’) ;

print -dtiff plot.tiff

% Subtractive Color Mixing

% MATLAB script to demonstrate effect of color subtraction

clear

m = zeros(256,256) ; % Size of the image

y = zeros(256,256) ;

c = zeros(256,256) ;

r = zeros(256,256) ;

g = zeros(256,256) ;

b = zeros(256,256) ;

x1=92; y1=80; % centre of the magenta circle

x2=92; y2=175; % centre of the yellow circle

x3=175; y3=128; % centre of the cyan circle

radius = 70 ;

for x= 0:255

for y = 0:255

p = x+1 ; q = y+1 ;

d1 = sqrt((x1-x)^2 + (y1-y)^2) ; % distance from the center of the magenta circle

d2 = sqrt((x2-x)^2 + (y2-y)^2) ; % distance from the center of the yellow circle

d3 = sqrt((x3-x)^2 + (y3-y)^2) ; % distance from the center of the cyan circle

m(p,q) = d1 < radius ; % magenta circle

y(p,q) = d2 < radius ; % yellow circle

c(p,q) = d3 < radius; % cyan circle

% The color assignment is not correct

if ((m(p,q)==0) & (y(p,q)==0) & (c(p,q)==0))

r(p,q)=1 ; g(p,q)=1 ; b(p,q)= 1 ; % white pixel

elseif ((m(p,q)==0) & (y(p,q)==0) & (c(p,q)==1))

r(p,q)=0 ; g(p,q)=1 ; b(p,q)= 1 ; % cyan pixel

elseif ((m(p,q)==0) & (y(p,q)==1) & (c(p,q)==0))

r(p,q)=1 ; g(p,q)=1 ; b(p,q)= 0 ; % yellow pixel

elseif ((m(p,q)==0) & (y(p,q)==1) & (c(p,q)==1))

r(p,q)=0 ; g(p,q)=1 ; b(p,q)= 0 ; % green pixel

elseif ((m(p,q)==1) & (y(p,q)==0) & (c(p,q)==0))

r(p,q)=1 ; g(p,q)=0 ; b(p,q)= 1 ; % magenta pixel

elseif ((m(p,q)==1) & (y(p,q)==0) & (c(p,q)==1))

r(p,q)=0 ; g(p,q)=0 ; b(p,q)= 1 ; % blue pixel

elseif ((m(p,q)==1) & (y(p,q)==1) & (c(p,q)==0))

r(p,q)=1 ; g(p,q)=0 ; b(p,q)= 0 ; % red pixel

elseif ((m(p,q)==1) & (y(p,q)==1) & (c(p,q)==1))

r(p,q)=0 ; g(p,q)=0 ; b(p,q)= 0 ; % black pixel

end

im = zeros(256,256,3) ;

im(:,:,1) = double(r) ;

im(:,:,2) = double(g) ;

im(:,:,3) = double(b) ;

imshow(im)

imwrite(im,’F:\im.tif’,’tif’) ;

print -dtiff plot.tiff

% Subtractive Color Mixing

% MATLAB script to demonstrate effect of color subtraction

clear

r = zeros(200,256) ; % Size of the image

g = zeros(200,256) ;

b = zeros(200,256) ;

x1=100; y1=100; % centre of the red circle

x2=100; y2=155;

radius = 55 ;

for x= 0:199

for y = 0:255

p = x+1 ; q = y+1 ;

d1 = (x1-x)^2 + (y1-y)^2 ;

outside = 1 ; % the point is outside all circles by default

if (sqrt(d1) < radius)

r(p,q) = 1 ; g(p,q)=1 ; % yellow circle

outside = 0 ; % The point is outside red circle

end

d2 = (x2-x)^2 + (y2-y)^2 ;

if (sqrt(d2) < radius)

b(p,q) = 1 ; % blue circle

outside = 0 ; % The point is outside red circle

if (sqrt(d1) < radius)

r(p,q) = 1 ; g(p,q)=1 ; b(p,q)=1 ; % green circle

end

if outside == 1

r(p,q) = 1 ; g(p,q)=1 ; b(p,q)=1 ; % white circle

end

im = zeros(200,256,3) ;

im(:,:,1) = double(r) ;

im(:,:,2) = double(g) ;

im(:,:,3) = double(b) ;

imshow(im)

imwrite(grating1,’F:\grat1.tif’,’tif’) ;

print -dtiff plot.tiff

% Nonuniform Color Space Experiments

% MATLAB script to demonstrate that RGB color space is non-uniform

% First we generate two colors with a distance of 10 in the blue region

clear

r = zeros(256,256) ; g = zeros(256,256) ; b = zeros(256,256) ; % Size of the image

x1=128; y1=128; % center of the circles

r1 = 35 ; r2 = 70 ;

for x= 0:255

for y = 0:255

p = x+1 ; q = y+1 ;

d = (x1-x)^2 + (y1-y)^2 ;

% Initializing different pixels

if (sqrt(d) < r1)

% g(p,q) = 1 ; % red circle

r(p,q) = 0.20 ; g(p,q)=0.6 ; b(p,q)=0.2 ;

elseif (sqrt(d) < r2)

% g(p,q) = 0.98 ; % The point is inside outer circle

r(p,q) = 0.20 ; g(p,q)=0.62 ; b(p,q)=0.2 ;

else

r(p,q)=1 ; g(p,q)=1; b(p,q)=1 ; % white background

end

im = zeros(256,256,3) ;

im(:,:,1) = double(r) ; im(:,:,2) = double(g) ; im(:,:,3) = double(b) ;

imshow(im)

imwrite(im,’F:\image1.tif’,’tif’) ;

% We now generate two colors with a distance of 10 in the blue region

clear

r = zeros(256,256) ; g = zeros(256,256) ; b = zeros(256,256) ; % Size of the image

x1=128; y1=128; % center of the circles

r1 = 35 ; r2 = 70 ;

for x= 0:255

for y = 0:255

p = x+1 ; q = y+1 ;

d = (x1-x)^2 + (y1-y)^2 ;

% Initializing different pixels

if (sqrt(d) < r1)

r(p,q) = 0.20 ; g(p,q)=0.2 ; b(p,q)=0.6 ;

elseif (sqrt(d) < r2)

r(p,q) = 0.20 ; g(p,q)=0.2 ; b(p,q)=0.62 ; % The point is inside outer circle

else

r(p,q) = 1 ; g(p,q)=1; b(p,q)=1 ; % white background

end

im = zeros(256,256,3) ;

im(:,:,1) = double(r) ; im(:,:,2) = double(g) ; im(:,:,3) = double(b) ;

imshow(im)

imwrite(im,’F:\image2.tif’,’tif’) ;

r = 0.2, g=0.6, 0.62, b=0.2 r = 0.2, g=0.2, b=0.6, 0.62

r=51, g=153,158, b=51 r=51, g=51, b=153,158

R=0, G=1,0.98, B=0 R=0, G=0, B=1,0.98 R=1,0.98, G=B=0

Chapter 4

Multimedia Data Acquisition

% Audio signal and its spectrum

clear;

%fname = 'F:\data\audio\Julien.raw' ; % Name of the Input file

fname = 'F:\imageresearch\Nelly3.raw' ; % Name of the Input file

fid=fopen(fname,'r');

no_samples=inf;

SigOrig=fread(fid,no_samples,'int');

SigOrig=SigOrig';

fclose(fid);

N=length(SigOrig);

plot(SigOrig) ;

Sig = SigOrig/max(abs(SigOrig)) ;

wavwrite(sig,fs,nb,'f:\test1_enhance.wav'); % Write the output signal as a wav file

% Audio signal and its spectrum

infile = 'F:\data\audio\test_audio1.wav' ; % Name of the Input file

% Name of the Input file

[x, Fs, bits]=wavread(infile);

%Plotting the input signal

plot([1:length(x)]/Fs,x) ;

title ('original speech signal')

xlabel('Time (in seconds)'); ylabel('Normalized Amplitude');

print -dtiff plot.tiff

hann_wind = round(2048*0.8) ;

[Pd,freq] = psd(x,2048,Fs,[],hann_wind); % power spectral density of the input signal

plot(freq/1000,10*log10(Pd))

xlabel('Frequency (in KHz)'); ylabel('Power (in dB)');

title('Power Spectrum of the Audio signal');

print -dtiff plot.tiff

% Audio signal and its spectrum

infile = 'F:\data\audio\GSS22KM8B.wav' ; % Name of the Input file

% Name of the Input file

[x, Fs, bits]=wavread(infile);

%Plotting the input signal

plot([1:length(x)]/Fs,x) ;

title ('Speech signal')

xlabel('Time (in seconds)'); ylabel('Normalized Amplitude');

print -dtiff plot.tiff

blocklen = 512 ;

hann_wind = round(blocklen*0.8) ;

[Pd,freq] = psd(x,blocklen,Fs,[],hann_wind); % power spectral density of the input signal

plot(freq/1000,10*log10(Pd))

xlabel('Frequency (in KHz)'); ylabel('Power (in dB)');

title('Power Spectrum of the Audio signal');

print -dtiff plot.tiff

% Audio signal and its spectrum

%Lowpass filtering

%MATLAB Program for Lowpass filtering

infile = 'F:\data\audio\bell1.wav' ;

% Name of the Input file

[x, Fs, bits]=wavread(infile);

%Plotting the input signal

plot([1:length(x)]/1000,128*x) ;

title ('original speech signal')

xlabel('Samples (x1000)'); ylabel('Sample Values');

print -dtiff plot.tiff

%Calculating Power spectral density

hann_wind = round(2048*0.8) ;

[Px,F] = psd(x,2048,Fs,[],hann_wind); % power spectral density of the input signal

plot(F/1000,10*log10(Px))

xlabel('Frequency (in KHz)'); ylabel('Power Spectral Density (in dB)');

print -dtiff plot.tiff

% Reconstructing sampled audio signal

% MATLAB script to demonstrate effect of aliasing

interval = 1/8000 ; %sampling interval in sec

t =[0:interval/30:0.001];

s = size(t) ;

samp = [0:s(2)-1]/30 ;

f = 5*cos(2*pi*4500*t) ;

f1 = 5*cos(2*pi*3500*t) ;

%plot(t,f,t,f1)

plot(samp,f,samp,f1)

xlabel(‘Sampling instances’)

ylabel(‘Signal values’),grid

print -dtiff plot.tiff

% Reconstructing sampled image

% MATLAB script to demonstrate effect of aliasing

x =[0:0.1:5];

y =[0:0.1:5];

grating1 = zeros(51,51);

grating2 = zeros(51,51);

grating3 = zeros(51,51);

for p= 1:51

for q = 1:51

grating1(p,q) = 127*cos(2*pi*(4*x(p)+6*y(q))) ;

grating2(p,q) = 127*cos(2*pi*(4*x(p)-4*y(q))) ; % The correct image

grating3(p,q) = 127*cos(2*pi*(4*x(p)+4*y(q))) ;

end

imshow(grating1)

imwrite(grating1,’F:\grat1.tif’,’tif’) ;

imwrite(grating2,’F:\grat2.tif’,’tif’) ;

imwrite(grating3,’F:\grat3.tif’,’tif’) ;

diff = grating1-grating2 ;

sdiff1 = sum(sum(abs(diff))) ;

diff = grating1-grating3 ;

sdiff2 = sum(sum(abs(diff))) ;

4x+6y 4x-4y 4x+4y

x =[0:0.05:5];

y =[0:0.05:5];

grating1 = zeros(101,101);

grating2 = zeros(101,101);

grating3 = zeros(101,101);

for p= 1:101

for q = 1:101

grating1(p,q) = 127*cos(2*pi*(4*x(p)+6*y(q))) ;

grating2(p,q) = 127*cos(2*pi*(4*x(p)-4*y(q))) ;

grating3(p,q) = 127*cos(2*pi*(4*x(p)+4*y(q))) ;

end

imshow(grating)

imwrite(grating1,’F:\grat1.tif’,’tif’) ;

imwrite(grating2,’F:\grat2.tif’,’tif’) ;

imwrite(grating3,’F:\grat3.tif’,’tif’) ;

print -dtiff plot.tiff

4x+6y 4x-4y 4x+4y

% Anti-aliasing Filter

Wp=3200; Ws=4000; Gp=-2; Gs=-20 ;

%Ideal Filter

mag0 = [ones(1,4001) zeros(1,4000)] ;

%Butterworth Filter

[n, Wc] = buttord(Wp,Ws,-Gp,-Gs,’s’) ;

[num,den] = butter(n,Wc,’s’) ;

w=[0:1:8000]; w = w’;

[mag1,phase1,w] = bode(num,den,w);

plot(w,mag1)

xlabel(‘Frequency (in Hz)’)

ylabel(‘Filter gain’),grid

print -dtiff plot.tiff

% Anti-aliasing Filter

Wp=3200; Ws=4000; Gp=-2; Gs=-20 ;

%Ideal Filter

mag0 = [ones(1,4001) zeros(1,4000)] ;

phase0 = [ones(1,4001) zeros(1,4000)] ;

%Butterworth Filter

[n, Wc] = buttord(Wp,Ws,-Gp,-Gs,’s’) ;

[num,den] = butter(n,Wc,’s’) ;

w=[0:1:8000]; w = w’;

[mag1,phase1,w] = bode(num,den,w);

%Chebyshev-1

[n, Wp] = cheb1ord(Wp,Ws,-Gp,-Gs,’s’) ;

[num,den] = cheby1(n,-Gp,Wp,’s’) ;

w=[0:1:8000]; w = w’;

[mag2,phase2,w] = bode(num,den,w);

%Chebyshev-2

[n, Ws] = cheb2ord(Wp,Ws,-Gp,-Gs,’s’) ;

[num,den] = cheby2(n,-Gs,Ws,’s’) ;

w=[0:1:8000]; w = w’;

[mag3,phase3,w] = bode(num,den,w);

plot(w,mag0,w,mag1,w,mag2,w,mag3)

xlabel(‘Frequency in Hz’)

ylabel(‘Amplitude Response’),grid

print -dtiff plot.tiff

plot(w,phase0,w,phase1,w,phase2,w,phase3)

xlabel(‘Frequency in Hz’)

ylabel(‘Phase Response’),grid

print -dtiff plot.tiff

% Verification

clear;

w = [0:10:8000];

num = 4.909e045 ;

den = [1, 2.714e004, 3.682e008, 3.298e012, 2.171e016, 1.107e020, 4.493e023, 1.47e027, 3.874e030, 8.129e033, 1.322e037, 1.579e040, 1.245e043, 4.909e045] ;

[mag1,phase1] = bode(num,den,w) ;

num = 2.742e016 ;

den = [1, 2261, 1.536e007, 2.272e010, 4.817e013, 2.742e016 ] ;

[mag2,phase2] = bode(num,den,w) ;

num = [1962, - 1.36e-008, 1.196e011, - 0.2839, 1.457e018] ;

den = [1, 1.267e004, 7.84e007, 3.215e011, 7.672e014, 1.457e018] ;

[mag3,phase3] = bode(num,den,w) ;

plot(w,mag1,w,mag2,w,mag3);

grid;

xlabel('w (frequency)')

ylabel('Filter Gain')

%title('Amplitude Spectrum')

plot(w,phase1,w,phase2,w,phase3);

grid;

xlabel('w (frequency)')

ylabel('Phase of H(w)')

title('Phase Spectrum')

%print -dtiff plot.tiff

2-D Filtering

% Digitization of a sinusoidal signal

clear

% generating a sinusoidal function

t = 0:0.008:0.8 ;

g = sin(2*pi*t) ;

plot(10*t,g) ;

xlabel(‘Sampling Time’)

ylabel(‘Signal Value’)

print -dtiff plot.tiff

clf

t1 = 0:0.1:0.8 ;

g1 = sin(2*pi*t1) ;

g1 = [0 0.5878 0.9511 0.9511 0.5878 0.0000 -0.5878 -0.9511 -0.9511] ;

qg1 = [0.125 0.625 0.875 0.875 0.625 0.125 -0.625 -0.875 -0.875] ;

e = g1 – qg1 ;

stem([0:8],e,'filled'), grid

xlabel(‘Sampling Time’)

ylabel(‘Quantization Error’)

print -dtiff plot.tiff

% Quality versus Bits/sample

infile = 'F:\imageresearch\lena.tiff';

frame = imread(infile) ; % frame1 is a 2-D array

imshow(frame)

framed = double(frame) ;

%imshow(frame6)

imwrite(frame6,’F:\research\lena6.tif’,’tif’) ;

frame6 = zeros(512,512) ;

frame5 = zeros(512,512) ;

frame4 = zeros(512,512) ;

frame3 = zeros(512,512) ;

frame2 = zeros(512,512) ;

frame1 = zeros(512,512) ;

frame6 = uint8(4*floor(framed/4)+2) ;

frame5 = uint8(8*floor(framed/8)+4) ;

frame4 = uint8(16*floor(framed/16)+8) ;

frame3 = uint8(32*floor(framed/32)+16) ;

frame2 = uint8(64*floor(framed/64)+32) ;

frame1 = uint8(128*floor(framed/128)+64) ;

imwrite(frame6,’F:\imageresearch\lena6.tif’,’tif’) ;

imwrite(frame5,’F:\imageresearch\lena5.tif’,’tif’) ;

imwrite(frame4,’F:\imageresearch\lena4.tif’,’tif’) ;

imwrite(frame3,’F:\imageresearch\lena3.tif’,’tif’) ;

imwrite(frame2,’F:\imageresearch\lena2.tif’,’tif’) ;

imwrite(frame1,’F:\imageresearch\lena1.tif’,’tif’) ;

(a)

(b) (c) (d)

(e) (f) (g)

(a) 8 bit, 256 levels, (b) 6 bit, 64 levels, (c) 5 bit, 32 levels, (d) 4 bit, 16 levels,

(e) 3 bit, 8 levels, (f) 2 bit, 4 levels, (g) 1 bit, 2 levels.

Chapter 5

Image Transforms

% Basis Functions of DFT

% ylabel(['k= ' int2str(m)])

clear

x = 0:7 ;

wr = zeros(8,8) ;

wi = zeros(8,8) ;

for m=0:7

v = zeros(1,8) ;

v(m+1)= 1 ;

w = ifft(v) ;

wr(m+1,:) = real(w)*sqrt(8) ;

wi(m+1,:) = imag(w)*sqrt(8) ;

end

% 8-point real basis functions

subplot(8,1,1),stem(x,wr(1,:),'filled'), grid

ylabel(‘k=0’)

subplot(8,1,2),stem(x,wr(2,:),'filled'), grid

ylabel(‘k=1’)

subplot(8,1,3),stem(x,wr(3,:),'filled'), grid

ylabel(‘k=2’)

subplot(8,1,4),stem(x,wr(4,:),'filled'), grid

ylabel(‘k=3’)

subplot(8,1,5),stem(x,wr(5,:),'filled'), grid

ylabel(‘k=4’)

subplot(8,1,6),stem(x,wr(6,:),'filled'), grid

ylabel(‘k=5’)

subplot(8,1,7),stem(x,wr(7,:),'filled'), grid