Matlab Image Processing

https://initpara.wordpress.com/2010/02/05/matlab-image-processing-232y3pcqwxodx-11/

 

 

 

New Updated News:

– Scroll below: Learn Matlab function for active contours.

– Updates at Below: Learn Matlab Image Processing – Edge detection Function with example

– Visit two more links for  Matlab Image Processing with examples:

http://knol.google.com/k/mayur-vora/learn-matlab-image-processing-with/232y3pcqwxodx/12#

http://knol.google.com/k/learn-matlab-image-processing-with-example-part-2#

———————————————————————————————————-

Active Contours in Matlab (Image Processing):

* Inbuilt Matlab functions have been made use of in implementing the
  below code.

close all
clear all
clc

% read the input image
inpImage =imread('coin_1.jpg');
% size of image
[rows cols dims] = size(inpImage);

if dims==3
    inpImage=double(rgb2gray(inpImage));
else
    inpImage=double(inpImage);
end

% Gaussian filter parameter
sigma=1.2;
% Gaussian filter
G=fspecial('gaussian',15,sigma);
% Gaussian smoothed image
inpImage=conv2(inpImage,G,'same');


% gradient of image
[gradIX,gradIY]=gradient(inpImage);
absGradI=sqrt(gradIX.^2+gradIY.^2);


% higher dimensional embedding function phi whose zero level set is our
% contour
% radius of circle - initial embedding function
% radius=min(floor(0.45*rows),floor(0.45*cols));
[u,v] = meshgrid(1:cols, 1:rows);
phi = ((u-cols/2)/(floor(0.45*cols))).^2+((v-rows/2)/(floor(0.45*rows))).^2-1;

% edge-stopping function
g = 1./(1+absGradI.^2);
% gradient of edge-stopping function
[gx,gy]=gradient(g);

% gradient descent step size
dt=.4;

% number of iterations after which we reinitialize the surface
num_reinit=10;

phiOld=zeros(rows,cols);

% number of iterations
iter=0;

while(sum(sum(abs(phi-phiOld)))~=0)
    % gradient of phi
    [gradPhiX gradPhiY]=gradient(phi);
    % magnitude of gradient of phi
    absGradPhi=sqrt(gradPhiX.^2+gradPhiY.^2);
    % normalized gradient of phi - eliminating singularities
    normGradPhiX=gradPhiX./(absGradPhi+(absGradPhi==0));
    normGradPhiY=gradPhiY./(absGradPhi+(absGradPhi==0));

    [divXnormGradPhiX divYnormGradPhiX]=gradient(normGradPhiX);
    [divXnormGradPhiY divYnormGradPhiY]=gradient(normGradPhiY);
    % curvature is the divergence of normalized gradient of phi
    K = divXnormGradPhiX + divYnormGradPhiY;
    % dPhiBydT
    dPhiBydT =( g.*K.*absGradPhi + g.*absGradPhi + (gx.*gradPhiX+gy.*gradPhiY) );
    phiOld=phi;
    % level set evolution equation
    phi = phi + ( dt * dPhiBydT );
    iter=iter+1;
    if mod(iter,num_reinit)==0
        % reinitialize the embedding function after num_reinit iterations
        phi=sign(phi);
        phi = double((phi > 0).*(bwdist(phi < 0)) - (phi < 0).*(bwdist(phi > 0)));
    end
    if mod(iter,10)==0
        pause(0.05)
        iter
        imagesc(inpImage)
        colormap(gray)
        hold on
        contour(phi,[0 0],'r')
        %         close all
        %         surf(phi)
        %         pause
    end
end

Image Thresholding

This demo shows how a image looks like after thresholding. The percentage of the thresholding means the threshold level between the maximum and minimum indesity of the initial image. Thresholding is a way to get rid of the effect of noise and to improve the signal-noise ratio. That is, it is a way to keep the significant imformation of the image while get rid of the unimportant part (under the condition that you choose a plausible thresholding level). In the Canny edge detector part, you will see that, before thinning, we first do some thresholding. You are encouraged to do thinning without thresholding and to see what is the advantage of thresholding.

lena.gif	10% threshold

20% threshold	30% threshold

Image Thresholding Matlab Codes

This program show the effect of thresholding. The output are four subfigures shown in the same figure:

• Subfigure 1: The initial “lena”
• Subfigure 2: Threshold level is one alfa
• Subfigure 3: Threshold level is two alfa
• Subfigure 4: Threshold level is three alfa

The MATLAB codes:
%%%%%%%%%%%%% The main.m file %%%%%%%%%%%%%%
clear;
% Threshold level parameter alfa:
alfa=0.1;% less than 1/3

[x,map]=gifread('lena.gif');
ix=ind2gray(x,map);
I_max=max(max(ix));
I_min=min(min(ix));
level1=alfa*(I_max-I_min)+I_min;
level2=2*level1;
level3=3*level1;
thix1=max(ix,level1.*ones(size(ix)));
thix2=max(ix,level2.*ones(size(ix)));
thix3=max(ix,level3.*ones(size(ix)));
figure(1);colormap(gray);
subplot(2,2,1);imagesc(ix);title('lena');
subplot(2,2,2);imagesc(thix1);title('threshold one alfa');
subplot(2,2,3);imagesc(thix2);title('threshold two alfa');
subplot(2,2,4);imagesc(thix3);title('threshold three alfa');
%%%%%%%%%%%%% End of the main.m file %%%%%%%%%%%%%%

Gaussian function 

These demos show the basic effects of the (2D) Gaussian filter: smoothing the image and wiping off the noise. Generally speaking, for a noise-affected image, smoothing it by Gaussian function is the first thing to do before any other further processing, such as edge detection. The effectiveness of the gaussian function is different for different choices of the standard deviation sigma of the Gaussian filter. You can see this from the following demos.

Smoothing nonnoisy image

lena.gif	filtered with sigma = 3	filtered with sigma = 1

Noise cancelling

noisy lena	filtered with sigma = 3	filtered with sigma =1

 

(Noise is generated by matlab function 0.3*randn(512))

Gaussian filter study matlab codes

This program show the effect of Gaussian filter. The output are four subfigures shown in the same figure:

• Subfigure 1: The initial noise free “lena”
• Subfigure 2: The noisy “lena”
• Subfigure 3: Filtered the initial “lena”
• Subfigure 4: Filtered the noisy “lena”

  The matlab codes:
  %%%%%%%%%%%%% The main.m file  %%%%%%%%%%%%%%%
  clear;
  % Parameters of the Gaussian filter:
  n1=10;sigma1=3;n2=10;sigma2=3;theta1=0;
  % The amplitude of the noise:
  noise=0.1;
  
  [w,map]=gifread('lena.gif');               
  x=ind2gray(w,map);
  filter1=d2gauss(n1,sigma1,n2,sigma2,theta);
  x_rand=noise*randn(size(x));
  y=x+x_rand;
  f1=conv2(x,filter1,'same');
  rf1=conv2(y,filter1,'same');
  figure(1);
  subplot(2,2,1);imagesc(x);
  subplot(2,2,2);imagesc(y);
  subplot(2,2,3);imagesc(f1);
  subplot(2,2,4);imagesc(rf1);
  colormap(gray);
  %%%%%%%%%%%%%% End of the main.m file %%%%%%%%%%%%%%%
  
  
  %%%%%%% The functions used in the main.m file %%%%%%%
  % Function "d2gauss.m":
  % This function returns a 2D Gaussian filter with size n1*n2; theta is 
  % the angle that the filter rotated counter clockwise; and sigma1 and sigma2
  % are the standard deviation of the gaussian functions.
  function h = d2gauss(n1,std1,n2,std2,theta)
  r=[cos(theta) -sin(theta);
     sin(theta)  cos(theta)];
  for i = 1 : n2 
      for j = 1 : n1
          u = r * [j-(n1+1)/2 i-(n2+1)/2]';
          h(i,j) = gauss(u(1),std1)*gauss(u(2),std2);
      end
  end
  h = h / sqrt(sum(sum(h.*h)));
  
  % Function "gauss.m":
  function y = gauss(x,std)
  y = exp(-x^2/(2*std^2)) / (std*sqrt(2*pi));
  %%%%%%%%%%%%%% end of the functions %%%%%%%%%%%%%%%%

Canny edge detector 

There are some results of applying Canny edge detector to real image (The black and white image “lena.gif” we used here was obtained by translating from a color lena.tiff using matlab. So it might not be the standard BW “lena”.) The thresholding parameter alfa is fix as 0.1. The size of the filters is also fixed as 10*10.

These images are all gray images though they might seem a little strange in your browser. To see them more clearly, just click these images and you will find the difference especially from the “result images”, that is, the titled “after thinning” ones. The safe way to see the correct display of these images is to grab these images and show them by “xv” or “matlab”. While, you are encouraged to use the given matlab codes and get these images in matlab by yourself. Try to change the parameters to get more sense about how these parameters affect the edge detection.

The results of choosing the standard deviation sigma of the edge detectors as 3.

lena.gif	vertical edges	horizontal edges

norm of the gradient	after thresholding	after thinning

The results of choosing the standard deviation sigma of the edge detectors as 1.

lena.gif	vertical edges	horizontal edges

norm of the gradient	after thresholding	after thinning

Canny edge detector algorithm matlab codes

This part gives the algorithm of Canny edge detector. The outputs are six subfigures shown in the same figure:

• Subfigure 1: The initial “lena”
• Subfigure 2: Edge detection along X-axis direction
• Subfigure 3: Edge detection along Y-axis direction
• Subfigure 4: The Norm of the image gradient
• Subfigure 5: The Norm of the gradient after thresholding
• Subfigure 6: The edges detected by thinning

The matlab codes:
%%%%%%%%%%%%% The main.m file  %%%%%%%%%%%%%%%
clear;
% The algorithm parameters:
% 1. Parameters of edge detecting filters:
%    X-axis direction filter:
     Nx1=10;Sigmax1=1;Nx2=10;Sigmax2=1;Theta1=pi/2;
%    Y-axis direction filter:
     Ny1=10;Sigmay1=1;Ny2=10;Sigmay2=1;Theta2=0;
% 2. The thresholding parameter alfa:
     alfa=0.1;
     
% Get the initial image lena.gif
[x,map]=gifread('lena.gif');               
w=ind2gray(x,map);
figure(1);colormap(gray);
subplot(3,2,1);
imagesc(w,200);
title('Image: lena.gif');

% X-axis direction edge detection
subplot(3,2,2);
filterx=d2dgauss(Nx1,Sigmax1,Nx2,Sigmax2,Theta1);
Ix= conv2(w,filterx,'same');
imagesc(Ix);
title('Ix');

% Y-axis direction edge detection
subplot(3,2,3)
filtery=d2dgauss(Ny1,Sigmay1,Ny2,Sigmay2,Theta2);
Iy=conv2(w,filtery,'same'); 
imagesc(Iy);
title('Iy');

% Norm of the gradient (Combining the X and Y directional derivatives)
subplot(3,2,4);
NVI=sqrt(Ix.*Ix+Iy.*Iy);
imagesc(NVI);
title('Norm of Gradient');

% Thresholding
I_max=max(max(NVI));
I_min=min(min(NVI));
level=alfa*(I_max-I_min)+I_min;
subplot(3,2,5);
Ibw=max(NVI,level.*ones(size(NVI)));
imagesc(Ibw);
title('After Thresholding');

% Thinning (Using interpolation to find the pixels where the norms of 
% gradient are local maximum.)
subplot(3,2,6);
[n,m]=size(Ibw);
for i=2:n-1,
for j=2:m-1,
 if Ibw(i,j) > level,
 X=[-1,0,+1;-1,0,+1;-1,0,+1];
 Y=[-1,-1,-1;0,0,0;+1,+1,+1];
 Z=[Ibw(i-1,j-1),Ibw(i-1,j),Ibw(i-1,j+1);
    Ibw(i,j-1),Ibw(i,j),Ibw(i,j+1);
    Ibw(i+1,j-1),Ibw(i+1,j),Ibw(i+1,j+1)];
 XI=[Ix(i,j)/NVI(i,j), -Ix(i,j)/NVI(i,j)];
 YI=[Iy(i,j)/NVI(i,j), -Iy(i,j)/NVI(i,j)];
 ZI=interp2(X,Y,Z,XI,YI);
  if Ibw(i,j) >= ZI(1) & Ibw(i,j) >= ZI(2)
  I_temp(i,j)=I_max;
  else
  I_temp(i,j)=I_min;
  end
 else
 I_temp(i,j)=I_min;
 end
end
end
imagesc(I_temp);
title('After Thinning');
colormap(gray);
%%%%%%%%%%%%%% End of the main.m file %%%%%%%%%%%%%%%


%%%%%%% The functions used in the main.m file %%%%%%%
% Function "d2dgauss.m":
% This function returns a 2D edge detector (first order derivative
% of 2D Gaussian function) with size n1*n2; theta is the angle that
% the detector rotated counter clockwise; and sigma1 and sigma2 are the
% standard deviation of the gaussian functions.
function h = d2dgauss(n1,sigma1,n2,sigma2,theta)
r=[cos(theta) -sin(theta);
   sin(theta)  cos(theta)];
for i = 1 : n2 
    for j = 1 : n1
        u = r * [j-(n1+1)/2 i-(n2+1)/2]';
        h(i,j) = gauss(u(1),sigma1)*dgauss(u(2),sigma2);
    end
end
h = h / sqrt(sum(sum(abs(h).*abs(h))));

% Function "gauss.m":
function y = gauss(x,std)
y = exp(-x^2/(2*std^2)) / (std*sqrt(2*pi));

% Function "dgauss.m"(first order derivative of gauss function):
function y = dgauss(x,std)
y = -x * gauss(x,std) / std^2;
%%%%%%%%%%%%%% end of the functions %%%%%%%%%%%%%

——————————————————————————————————————————————-

MATLAN IMAGE PROCESSING : EDGE DETECTION FUNCTION 

Marr/Hildreth and Canny Edge Detection with m-code

Recall the wagon wheel test image used in the resampling example:

Marr/Hildreth edge detection is based on the zero-crossings of the Laplacian of the Gaussian operator applied to the image for various values of sigma, the standard deviation of the Gaussian. What follows is a mosaic of zero-crossings for four choices of sigma computed using the Matlab image processing toolbox. The top left is sigma=1, the top right is sigma=2, the bottom left is sigma=3 and the bottom right is sigma=4. (Matlab Laplacian of Gaussian edge detection normally selects a threshold so that only zero-crossings of sufficient strength are shown. Here, the threshold is forced to be zero so that all zero-crossings are reported, as is required by the Marr/Hildreth theory of edge detection.)

As most commonly implemented, Canny edge detection is based on extrema of the first derivative of the Gaussian operator applied to the image for various values of sigma, the standard deviation of the Gaussian. What follows is a mosaic of edge points for four choices of sigma computed using the Matlab image processing toolbox. The top left is sigma=1, the top right is sigma=2, the bottom left is sigma=3 and the bottom right is sigma=4. The Canny method uses two thresholds to link edge points. The Matlab implementation can estimate both thesholds automatically. For this example, I found the Matlab estimated thresholds to be somewhat conservative. Therefore, for the mosaic, I reduced the thresholds to 75 percent of their automatically estimated values.

You should observe that zero-crossings in Marr/Hildreth edge detection always form connected, closed contours (or leave the edge of the image). This comes, however, at the expense of localization, especially for larger values of sigma. Arguably, Canny edge detection does a better job of localization. Alas, with Canny edge detection, the edge segments can become disconnected

For above example, m-code is below:

%    Image Processing Toolbox Version 5.4(R2007a)

A = imread('/ai/woodham/public_html/cpsc505/images/wheel.tif');

% Marr/Hildreth edge detection
% with threshold forced to zero
MH1 = edge(A,'log',0,1.0);
MH2 = edge(A,'log',0,2.0);
MH3 = edge(A,'log',0,3.0);
MH4 = edge(A,'log',0,4.0);

% form mosaic
EFGH = [ MH1 MH2; MH3 MH4];

%% show mosaic in Matlab Figure window
%log = figure('Name','Marr/Hildreth: UL: s=1  UR: s=2  BL: s=3 BR: s=4');
%iptsetpref('ImshowBorder','tight');
%imshow(EFGH,'InitialMagnification',100);

% Canny edge detection
[C1, Ct1] = edge(A,'canny',[],1.0);
[C2, Ct2] = edge(A,'canny',[],2.0);
[C3, Ct3] = edge(A,'canny',[],3.0);
[C4, Ct4] = edge(A,'canny',[],4.0);

% Recompute lowering both automatically computed
% thresholds by fraction k
k = 0.75
C1 = edge(A,'canny',k*Ct1,1.0);
C2 = edge(A,'canny',k*Ct2,2.0);
C3 = edge(A,'canny',k*Ct3,3.0);
C4 = edge(A,'canny',k*Ct4,4.0);

% form mosaic
ABCD= [ C1 C2; C3 C4 ];

% show mosaic in Matlab Figure window
%canny = figure('Name','Canny: UL: s=1  UR: s=2  BL: s=3 BR: s=4');
%iptsetpref('ImshowBorder','tight');
%imshow(ABCD,'InitialMagnification',100);

% write results to file
% Note: Matlab no longer reads/writes GIF files, owing to licence
%       restrictions.  Translation from TIF to GIF was done
%       manually (with xv) for inclusion in the web example page
imwrite(ABCD,'/ai/woodham/World/cpsc505/images/canny.tif','tif
imwrite(EFGH,'/ai/woodham/World/cpsc505/images/log.tif','tif