Basic 2D Signals

Contents

• Start
• Define Important Variables
• Rectangular Function
• Take Fourier Transform
• Pyramid Function
• Fourier Transform of pyramid
• Cylinder
• Fourier Transform of Cylinder
• Cone
• Fourier Transform of Cone
• Airy PSF
• Fourier Transform of Airy PSF
• Gaussian
• Fourier Transform of Gaussian
• Peak
• Fourier Transform of Peak
• Exponential Decay
• Fourier Transform of Exponential Dacay

Start

close all;
clear all;
clc;

Define Important Variables

N = 64;         % Size of a sunction
D = N/2;        % to indicate origin at the center of the function
a = 8;          % radius for cylindrical function
w = 0.4;        % decay rate for exponential function
y = repmat(1:N,N,1);
x = y';
r = sqrt((D-x).^2+(D-y).^2);    % definition of radius
sig = 5;                        % Variance for gaussian function

Rectangular Function

re = zeros(N);
re(D-a:D+a-1,D-a:D+a-1) = w*ones(2*a);
mesh(re); title('2D Rectangular Function');

Take Fourier Transform

R = fft2(re);
R = fftshift(R);
figure;mesh(abs(R)); title('Fourier Transform of Rectangular function');

Pyramid Function

p = conv2(re,re,'same');
figure; mesh(p); title('2D Pyramidal Function');

Fourier Transform of pyramid

P = fft2(p);
P = fftshift(P);
figure; mesh(abs(P));title('Fourier Transform of pyramidal function');

Cylinder

c = w*double((r.^2)<a^2);
figure; mesh(c); title('2D Cylindrical Function');

Fourier Transform of Cylinder

C = fft2(c);
C = fftshift(C);
figure; mesh(abs(C));title('Fourier Transform of Cylindrical function');

Cone

co = conv2(c,c,'same');
figure; mesh(co); title('2D Conical Function');

Fourier Transform of Cone

CO = fft2(co);
CO = fftshift(CO);
figure; mesh(abs(CO)); title('Fourier Transform of Conical function');

Airy PSF

psf = (airy(w*r/2).^2)/pi;
figure; mesh(psf); title('2D Point Spread Function');

Fourier Transform of Airy PSF

PSF = fft2(psf);
PSF = fftshift(PSF);
figure; mesh(abs(PSF)); title('Fourier Transform of airy PSF function');

Gaussian

g = (2*pi*sig^2)*exp(-((r.^2))./(2*sig^2));
figure;mesh(g); title('2D Gaussian Function');

Fourier Transform of Gaussian

G = fft2(g);
G = fftshift(G);
figure; mesh(abs(G)); title('Fourier Transform of Gaussian function');

Peak

pk = 1./r;
pk(isinf(pk)) = 1;
figure; mesh(pk); title('2D Peak Function');

Fourier Transform of Peak

Pk = fft2(pk);
Pk = fftshift(Pk);
figure; mesh(abs(Pk)); title('Fourier Transform of peak function');

Exponential Decay

exd = exp(-w*r);
figure; mesh(exd); title('2D Exponential Decay Function');

Fourier Transform of Exponential Dacay

EXD = fft2(exd);
EXD = fftshift(EXD);
figure; mesh(abs(EXD));
title('Fourier Transform of Exponential Dacay function');

Published with MATLAB® 7.8

Morphological Operations

This program will demonstrate the various morphological operations used in image processing. For further information please read: Gonzalez, Rafael C, Digital image processing - 3rd ed. - New Delhi Pearson Education 2009

Contents

• Start
• Read the image
• Threshold
• Defining Structuring Element
• Image Dilation
• Image Erosion
• Closing Operation
• Opening Operation
• Region Filling
• Boundary Detection

Start

clear all;
close all;
clc;

Read the image

im = imread('coins.png');
figure; imshow(im);title('original Image');

Thresholding

this step will convert the given gray cycle image into binary image

x = im>128;
figure;imshow(x); title('Binary Image');

Defining Structuring Element

In mathematical morphology, a structuring element (s.e.) is a shape, used to probe or interact with a given image, with the purpose of drawing conclusions on how this shape fits or misses the shapes in the image. It is typically used in morphological operations, such as dilation, erosion, opening, and closing, as well as the hit-or-miss transform. (Source: Wikipedia)

Smat1 = [0 1 0; 1 1 1; 0 1 0];
Smat2 = [1 1 1; 1 0 1; 1 1 1];
SE = strel(Smat1);
figure;
subplot 121;  imshow(Smat1); title('Element1');
subplot 122;  imshow(Smat2); title('Element2');

Image Dilation

The dilation operation usually uses a structuring element for probing and expanding the shapes contained in the input image.

xdilate = imdilate(x,SE);
figure;imshow(xdilate); title('Dilated Image');

Image Erosion

Visit here: http://homepages.inf.ed.ac.uk/rbf/HIPR2/erode.htm

xerode = imerode(x,SE);
figure;imshow(xerode); title('Eroded Image');

Closing Operation

Opening operation is same as close(A,B) = erode(dilate(A,B),B) It removes noise and irregularities inside the object

xclose = imerode(imdilate(x,SE),SE);
figure;imshow(xclose); title('Closed Image');

Opening Operation

Opening operation is same as open(A,B) = dilate(erode(A,B),B) It removes background noise as well as small objects from the scene

xopen = imdilate(imerode(x,SE),SE);
figure;imshow(xopen); title('Opened Image');

Region Filling

Visit below link for more information http://www.cis.rit.edu/class/simg782/lectures/lecture_04/lec782_05_04.pdf

B=[0,1,0;1,1,1;0,1,0];
maxit = 1000;
x = xclose;
sa=size(x);
X=zeros(sa(1),sa(2));
[m start] = min(min(x==1));
X(start)=1;
Y=zeros(sa(1),sa(2));
xc= ~x;
count = 1;
while ((min(X ~= Y) == 0) & (count < maxit))
    count=count+1;
    Y=X;
    X=imdilate(Y,B) & xc;
end
figure;imshow(X); title('Region Filling');

Boundary Detection

we can detect boundaries of the objects using this technique boundary = A - (A eroded by SE);

xhat = imerode(X,SE);
boundary = ~(X-xhat);
figure; imshow(boundary); title('Boundary Detection');

(Part II: Coming soon...)

Template matching using Normalized Cross Correlation

This program demonstrate the implementation of conventional cross correlation and normalized cross correlation metric to find the similarity score between template and the image portion.This program can be used for image registration to align the given images according to correlated pixels. Also this image demonstrate the use of correlation techniques to match the object. This method is very crude that it may fail most of the time if the object selected to find is a bright image.

Contents

• Start
• Read the image
• Read the Template
• Finding correlation map
• Finding Normalized Cross correlation map
• Display the object in image

Start

clear all;
close all;
clc;

Read the image

image = imread('dollar.tif');
x = double(image)/255;
imshow(x); pause(1);

Read the Template

template = imread('temp.tif');
temp = double(rgb2gray(template))/255;
imshow(temp); pause(1);

Finding correlation map

Cxt = conv2(x,rot90(conj(temp)),'same'); % Cross correlation
imshow(Cxt/max(max(Cxt))); pause(1)

Finding Normalised Cross correlation map

h = ones(size(temp));
Cgg = conv2(x,h,'same'); % finding the energy of signal
NCC = Cxt./Cgg;
imshow(NCC/max(max(NCC))); pause(0.5);

Display the object in image

Locate the maxima [m,n] from the figure and then display the box around the object

[a b] = max(max(NCC));
[c d] = max(NCC);
m = d(b);
n = b;
clear a b c d;
[a b] = size(x);
[c d] = size(temp);
c = round(c/2);
for i=1:a
    for j=1:b
        if (i==m-c && j>n-c && j<n+c)||(i==m+c && j>n-c && j<n+c)...
                || (j==n-c && i>m-c && i<m+c)||(j==n+c && i>m-c && i<m+c)
            x(i,j) = 0;
            x(i+1,j+1) = 0;
        end
    end
end
imshow(x);

Connected component Labelling

This program demonstrate the segmentation process by classifying the connected pixels in the image. This program can be used for automatic coins counting machine, if images of coins taken with Black background and the coins are places at a safe distance. to try this, take an image of coins placing them on black surface and use this image as input for this program

Contents

• Start
• Read the Image
• Convert it to binary image
• Connected Component Labelling Algorithm
• Label Conflicts Compensation
• Renaming the labels
• Display of the results

Start

clear all;
close all;
clc;

Read the Image

I = rgb2gray(imread('coin.jpg')); % read graycycle image
title('Original Image');
imshow(I);

Convert it to binary image

ima = I>128; %Thresholding
ima = bwareaopen(ima,30); %Neglect the connected region <30 pixels
figure; imshow(ima);

Connected Component Labelling Algorithm

[a b] = size(ima);
label = 0;
D = [];
for i=1:a
    for j = 1:b
        if ima(i,j) == 0
            Im(i,j) = 0;
        else
            switch ima(i-1,j)+10*ima(i, j-1)
                case 0
                    label = label+1;
                    Im(i,j) = label;
                case 1
                    Im(i,j) = Im(i-1, j);
                case 10
                    Im(i,j) = Im(i ,j-1);
                case 11
                    if Im(i-1,j) == Im(i, j-1)
                        Im(i,j) = Im(i, j-1);
                    else
                        Im(i,j) = Im(i-1, j);
                        Im(Im == Im(i-1,j)) = Im(i, j-1);
                    end
            end
        end
    end
end

Label Conflicts Compensation

for i = 1:2:length(D)
    m = minmax(D(i:i+1));
    Im(Im==m(2)) = m(1);
    d((i+1)/2) = D(i);
end

Renaming the labels

label = 1;
for i = 1:max(max(Im))
    if sum(sum(Im==i))>10
        Im(Im==i) = label;
        label = label+1;
    end
end

Display of the results

figure;imshow(10*uint8(Im));
t = ['Number of connected regions = ' num2str(max(max(Im)))];
title(t);
disp(['Number of coins in the image are: ', num2str(max(max(Im)))]);

Number of coins in the image are: 9

Image Processing by point processing

This program will demonstrate various image enhancement techniques using gray level manipulation

Contents

• Start
• Read the Image
• Digital Negative
• Log Transformation
• Power Law Transformation
• Contrast Stretching
• Gray Level Slicing

Start

clear all;
close all;
clc;

Read the Image

x1 = imread('coins.png');
x = double(x1)/255;
imshow(x);

Digital Negative

In this method: f(x,y) = L-Im(x,y), where L is the largest gray level in the image

x_neg = 1-x;
figure; imshow(x_neg)

Log Transformation

In this method: f(x,y) = c*log(Im(x,y)) where c is any constant and Im > 0
for c = 0.4

c = 0.4;
x_log = c*log(1+x);
figure; imshow(abs(x_log));

For c = 3

c = 3;
x_log = c*log(1+x);
figure; imshow(abs(x_log));

Power Law Transformation

In this method: f(x,y) = c*Im(x,y)^g where c and g are positive constants
for c = 0.5, g = 2

c = 0.5;
g = 2;
x_pow = c*x.^g;
figure; imshow(abs(x_pow));

for 1.5 = 2, g = 0.5

c = 1.5;
g = 0.5;
x_pow = c*x.^g;
figure; imshow(abs(x_pow));

Contrast Stretching

This program demonstrate piecewise continuous gray level transformation here we user three regions as: 0 to a1 transformed with the slope of m1 and a2 to 255 with slope of m1

m1 = 2.5;
m3 = 0.8;
a1 = 50;
a2 = 120;
c1 = 255 - 255*m3;
m2 = (m3*a2+c1-m1*a1)/(a2-a1);
c2 = m1*a1-m2*a1;
[row col] = size(x1);
x_contrast = ones(row,col);
for i = 1:row
    for j = 1:col
        if x1(i,j) < a1
            x_contrast(i,j) = x1(i,j)*m1;
        elseif x1(i,j) > a2
            x_contrast(i,j) = x1(i,j)*m3+c1;
        else
            x_contrast(i,j) = x1(i,j)*m2;
        end
    end
end
figure;imshow(uint8(x_contrast));
y = [(0:a1-1)*m1 ((a1:a2-1))*m2+c2 (a2:255)*m3+c1];
figure; plot(y);

Gray Level Slicing

This program demonstrate the gray level slicing operation on image
Without Background:

x_GLS = ones(row,col);
a = 48;
b = 80;
c = 190;
for i = 1:row
    for j = 1:col
        if x1(i,j) < a || x1(i,j) > b
            x_GLS(i,j) = c;
        else
            x_GLS(i,j) = 0;
        end
    end
end
figure;imshow(uint8(x_GLS));
y = [(0:a1-1)*0 (ones(1,a2-a1+1))*c (a2:255)*0];
figure; plot(y);

With Background:

c = 190;
for i = 1:row
    for j = 1:col
        if x1(i,j) < a || x1(i,j) >b
            x_GLS(i,j) = c;
        else
            x_GLS(i,j) = x1(i,j);
        end
    end
end
figure;imshow(uint8(x_GLS));
y = [(0:a1-1) (ones(1,a2-a1+1))*c (a2:255)];
figure; plot(y);

Spatial Domain Filtering

This program demonstrate the image filtering in spatial domain using various types of box filters

Contents

• Start
• Choose your mask
• Read your image
• Filter the image
• Display the result

Start

clear all;
close all;
clc;

Choose your mask

hLP = ones(3,3)/9;
hLap = [0 -1 0; -1 4 -1; 0 -1 0]; %Laplacian
hPreH = [-1 -1 -1; 0 0 0; 1 1 1]; %Prewitt
hPreV = [1 0 -1; 1 0 -1; 1 0 -1]; %Prewitt
hSobV = [1 0 -1; 2 0 -2; 1 0 -1]; %Sobel
hSobH = [-1 -2 -1; 0 0 0; 1 2 1]; %Sobel
hHP = [-1 -1 -1; -1 8 -1; 1 -1 -1]; %Sharpening
hHB = [-1 -1 -1; -1 8.2 -1; 1 -1 -1]; % High boost
hUnsharp = [-1 -1 -1; -1 1 -1; -1 -1 -1]/9; %Unsharp Masking

Read your image

x = imread('coins.png');
x = double(x)/255;

Filter the image

imLP = conv2(x,hLP,'same');
imLap = conv2(x,hLap,'same');
imPreH = conv2(x,hPreH,'same');
imPreV = conv2(x,hPreV,'same');
imPre = sqrt(imPreH.^2+imPreV.^2);
imSobH = conv2(x,hSobH,'same');
imSobV = conv2(x,hSobV,'same');
imSob = sqrt(imSobH.^2+imSobV.^2);
imHP = conv2(x,hHP,'same');
imHB = conv2(x,hHB,'same');
imUnsharp = conv2(x,hUnsharp,'same');

Display the result

Original Image

imshow(x);

Low Pass Image

figure; imshow(abs(imLP));

Laplacian Image

figure; imshow(abs(imLap));

Prewitt Operated Image

figure; imshow(abs(imPre));

Sobel Operated Image

figure; imshow(abs(imSob));

High Pass Image

figure; imshow(abs(imHP));

High Boost Filtering

figure; imshow(abs(imHB));

Unsharp Masking

figure; imshow(abs(imUnsharp));