Download Solutions for Homework 4 - Digital Image Processing | ECE 533 and more Assignments Digital Signal Processing in PDF only on Docsity!
Digital Image Processing ECE 533
Solutions Assignment 4
Department of Electrical and Computing Engineering, University of New Mexico.
Professor Majeed Hayat, [email protected]
March 16, 2008
1 Image enhancement using intensity transformations
- By inspecting the given figures we can see that the transformation is simply a polynomial transformation of degree γ for both the positive and negative intensity transformations. Let us call s( kp )the positive transformation and s( kn )for the neg- ative transformation. Therefore, the class of functions, parameterized by γ, γ > 0, that maps the normalized intensity levels rk ∈ { 0 , (^) L^1 − 1 ,... , 1 } into the normalized intensity levels sk = rk are:
s( kp ) = rγk s( kn ) = (1 − rk)γ^.
- Below we present a MATLAB function that implements both the intensity trans- formation as well as the contrast stretching for positive and negative images. In Figs.1 and 2 we show some results obtained after applying the transformations to an image. Note that: i) the gamma correction enhances the white intensities when γ < 1 and it enhances the black intensities when γ > 1 despite of considering the
(^00) 0.2 0.4 0.6 0.8 1
rk
n/n^0
a) γ=0.5 pos.
(^00) 0.2 0.4 0.6 0.8 1 0.005^ 0.
0.015^ 0.
0.025^ 0.
0.035^ 0.
sk
n/n^0
γ=0.5 neg.
(^00) 0.2 0.4 0.6 0.8 1 0.005^ 0.
0.015^ 0.
0.025^ 0.
0.035^ 0.
sk
n/n^0
b) c) γ=2 pos.
(^00) 0.2 0.4 0.6 0.8 1 0.005^ 0.
0.015^ 0.
0.025^ 0.
0.035^ 0.
sk
n/n^0
γ=2 neg.
(^00) 0.2 0.4 0.6 0.8 1 0.002^ 0.
0.006^ 0.
0.012 0. 0.014^ 0.
0.018^ 0.
sk
n/n^0
d) e)
Figure 1: Intensity transformations: positive and negative gamma corrections: a) Original image and its histogram. Positive, b), and negative, c), images corrected with γ = 0.5. Positive, d), and negative, e), images corrected with γ = 2.
positive or the negative image; ii) the histogram stretching produces a binary image when the slope α is infinite.
function ImgT=IntensityTransformation(Img,Type,Gamma,Inv) % % IntensityTransformation: function to modify the intensity of the image % Img according to the polynomial function: ImgT(x,y)=Img(x,y)^Gamma % % Inputs:
% Gamma: For type the ’Gammma’ is a scalar with the transformation % factor. For type ’CStretch’, Gamma=[m Alpha] with the threshold % value m and the slope Alpha of the correction. Alpha can be equal % to inf. For the type ’BinTrans’ is a 2x2 matrix with [rmin rmax; sL % sH], i.e., the range and the binary values. For the type ’LevBoost’ % is a vector of 1x3: [rmin rmax sH], i.e., the range and the % value at the boosted range % Inv: create the negative of the image (optional). 0: positive % image. 1: negative image. % Output: % ImgT: tranformed image
switch Type case ’Gamma ’ % If invertion parameter is not used if (nargin==3) ImgT=Img.^Gamma; end % If invertion parameter is used if (nargin==4) ImgT=(1Inv+(-1)^InvImg).^Gamma; end case ’CStretch’ m=Gamma(1); Alpha=Gamma(2); % If invertion parameter is not used if (nargin==3) if isfinite(Alpha) ImgT=1./(1+(m./(Img+eps)).^Alpha); else % infinity slope ImgT=double(Img>=m); end end % If invertion parameter is used if (nargin==4) if isfinite(Alpha) ImgT=1./(1+(m./(1Inv+(-1)^InvImg+eps)).^Alpha); else % infinity slope if Inv==
ImgT=double(Img>=m); else ImgT=double(Img<=m); end end end case ’BinTrans’ % Flag with ones in the desired range Gamma=Gamma/255; Flags=((Gamma(1,1)<=Img) & (Img<=Gamma(1,2))); % Create the new image negating the flag and adding ImgT=Gamma(2,1)double(~Flags)+Gamma(2,2)double(Flags); case ’LevBoost’ Gamma=Gamma(:)’/255; % Flag with ones in the desired range Flags=((Gamma(1,1)<=Img) & (Img<=Gamma(1,2))); % Create the new image negating the flag and adding ImgT=Img.double(not(Flags))+Gamma(1,3)double(Flags); end
2 Gray-level slicing
The function IntensityTransformation presented above implements also the gray-level slicing transformations. In Fig. 3 we show the results obtained after applying the required transformations to a sample image.
3 Bit-level slicing
Bit-level slicing is a technique used to partition the intensities of an image into n levels that are not gray-levels. The bit-level slicing technique highlights the contribution made by specific bits to the total image intensity. It can be said that only the five most significant bits (MSB) contain visually significant data while the three less significant bit (LSB) planes contribute the more subtle details. Recall from binary arithmetics that shifting i bits to the left (right) and padding
(^00 50 100 150 200 )
rk
n/n^0
(^00 50 100 150 200 )
rk
n/n^0
(^00 50 100 150 200 )
rk
n/n^0
(^00 50 100 150 200 )
rk
n/n^0
a) b)
(^00 64 128 )
rk
n/n^0
(^00 64 128 )
rk
n/n^0
c)
Figure 4: Bit-level slicing transformations. a) Effect of setting the i = 3 LSB of the intensity values of the image to zero. b) Effect of setting the i = 3 MSB of the intensity values of the image to zero. c) Effect of inverting the i = 3 MSB of the intensity values of the image.
same. Finally, if we invert the higher order bit planes then the result is that we are interchang- ing all the pixels values between one bit plane and its complement. If we invert the first i MSB, then all the intensity values that belong to the interval [n∗ 28 −i, (n+1)∗ 28 −i−1], n = 0 , 1 ,... , 2 i^ − 1 will be mapped to the interval [k ∗ 28 −i, (k + 1) ∗ 28 −i^ − 1], k = 2i^ − 1 − n and viceversa. Note that when we mantain the 2^8 −i^ LSB unchanged then all the inten- sity values keep the same positions within the mapped interval. Hence, the number of
intensities the number will not change, however the histogram will look reflected if we look at the level of the intervals determined by high-order bit planes, but it will remain unchanged within the intervals. The following code implement the three bit-plane operations commented above. In Fig. 4 the results of the processing are shown for the case of modifying i = 3 bits. It is easy to see that the mask uint8(bin2dec(’00001000’)*ones(M,N)) in conjunction with the and operation will activate the third order bit-plane of an image. clear all; clc; Img=imread(’IRImage2.gif’); N=size(Img,2); M=size(Img,1);
% Set the 3 LSB to zero Mask=uint8(bin2dec(’11111000’)*ones(M,N)); ImgT=bitand(Img,Mask);
% Set the 3 MSB to zero Mask=uint8(bin2dec(’00011111’)*ones(M,N)); ImgT=bitand(Img,Mask);
% Invert the 3 MSB ImgT=Img; for i=1:M for j=1:N for k=6: % Get the 3 MSB and negate them ImgT(i,j)=bitset(ImgT(i,j),k,not(bitget(ImgT(i,j),k))); end end end
4 Histogram equalization and histogram specification
- Histogram equalization is a mapping or redistribution of the histogram components on a new intensity scale. To obtain an uniform histogram it is required that pixel intensities mapped form, say, L groups of N/L pixels per group, where L is the number of allowed discrete gray-levels and N is the total number of pixels of the
- The code given in the assignment implements equation (1) in the following steps: (i) it creates the normalized histogram of the image (hist1); (ii) it produces the cumulated sum (CDF1) to generate the summation in (1), note that in this step if some of the intensities is missing in the image, then they are not mapped into the range; and (iii) it assigns the new intensities.
- The following code is an example on how to perform the histogram specification.
clear all; close all; img1 = imread(’phobos.jpg’);
%Transform to Uniform Distribution [hist1, bins1] = hist(double(img1(:)),0:255); hist1 = hist1./length(img1(:)); T = cumsum(hist1); img1eq = zeros(size(img1)); for i=0: img1eq(find(img1==i)) = T(i+1); end [hist1eq, bins1eq] = hist(double(255*img1eq(:)),0:255); hist1eq = hist1eq./length(img1(:)); S = cumsum(hist1eq);
%Specify New Histogram load HistogramV7.mat
%Compute New CDF’s from Specified Histogram (Iterative) G = (cumsum(z)/length(img1(:))); Ginv = zeros(size(G)); for k=1: dff = -1; m = 0; while(dff < 0) m = m+1; dff = G(m) - S(k); end Ginv(k) = m-1; end
img1mt = zeros(size(img1eq)); ieq = floor(255*img1eq); for i=0:
img1mt(find(ieq==i)) = Ginv(i+1); end
[hist1mt, bins1mt] = hist(double(img1mt(:)),0:255); img1mt = img1mt/255;
The main difference between histogram specification and histogram equalization is that in histogram specification one can freely specify the shape of the histogram for the processed image, while the histogram equalization seeks to produce a processed image with an uniform histogram.
5 Spatial filtering
Table 1: Summary of the performance metrics achieved by the LF and the MF for different mask sizes. The original image in Fig. 5a) was corrupted with Gaussian noise (GN) of zero mean and variance 0.01. In addition the image was corrupted with Salt & Pepper noise (SPN) affecting the 5% of the pixels. Metric 3 × 3 GN 9 × 9 GN 15 × 15 GN 3 × 3 SPN 9 × 9 SPN 15 × 15 SPN CP U LFtime 0.0700 0.0200 0.0200 0.0100 0.0100 0. CP U M Ftime 0.0800 0.2500 0.6400 0.0300 0.2200 0. M SELF^ 0.0019 0.0033 0.0056 0.0025 0.0035 0. M SEM F^ 0.0019 0.0014 0.0028 0.0002 0.0008 0. ρLF^ 0.1020 0.0399 0.0314 0.0909 0.0403 0. ρM F^ 0.1137 0.0407 0.0311 0.0391 0.0300 0. QLF^ 0.3931 0.5527 0.4902 0.3629 0.5256 0. QM F^ 0.3413 0.5276 0.4994 0.8997 0.7484 0.
To prove that a median filter (MF) is nonlinear consider to images I 1 (x, y) and I 2 (x, y). If we apply a median filter to each image we obtain: O 1 (x, y){I 1 } and O 2 (x, y){I 2 }. Now, consider the combined image I(x, y) = I 1 (x, y) + I 2 (x, y) and apply a median filter to it: O(x, y){I} = O(x, y){I 1 + I 2 } 6 = O 1 (x, y){I 1 } + O 2 (x, y){I 2 } because given an N × N mask we have to sort (rank) the intensities of the pixels within the mask and pick as the output of the (x, y) pixels the intensity value located in the middle of the rank. So, if
we add the images first then the ranking operation will not produce necessarily the same intensity value in the middle of the rank as the sum of the intensity values obtained after ranking individually the masks. So, superposition does not hold and the MF is nonlinear. We tested and compared the performance of the MF and an average spatial linear filter (LF) under two scenarios. In Fig. 5 we show the filtered images obtained by each filter, and in Table 1 we list the summary of the performance metrics achieved by them. Clearly, the MF takes a larger time than the LF to produce the result. Such situation is expected due to the complexity of the ranking operation. In addition, the results obtained show that for the evaluated image the MF outperforms the LF in almost every performance metric but the computing time. Finally, a visual evaluation also coincides with the objective results. The following MATLAB code is an example of the required implementation.
clear all; clc; kf=1; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % IMAGE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Read image and calculate its FFT-2D % The images available are: IRImage.gif, IRImage2.gif, and IRImage3.gif Img2Use=2; % Select the image to use switch Img2Use case 1 Img=double(imread(’IRImage1.gif’))/255; case 2 Img=double(imread(’IRImage2.gif’))/255; case 3 Img=double(imread(’IRImage3.gif’))/255; end N=size(Img,2); M=size(Img,1); for i=1: switch i case 1 ImgN=imnoise(Img,’Gaussian’,0,0.01); case 2 ImgN=imnoise(Img,’salt & pepper’); end
Nv=[3 9 15]; for k=1:length(Nv) % Compute the cputime t=cputime; ImgLF=filter2(fspecial(’average’,Nv(k)),ImgN); DtLF(i,k)=cputime-t; % Compute performance metrics MSELF(i,k)=mean2((Img-ImgLF).^2); RhoLF(i,k)=Roughness(ImgLF); QLF(i,k)=Qindex(Img,ImgLF); % Compute the cputime t=cputime; ImgMF=medfilt2(ImgN,[Nv(k) Nv(k)]); DtMF(i,k)=cputime-t; % Compute performance metrics MSEMF(i,k)=mean2((Img-ImgMF).^2); RhoMF(i,k)=Roughness(ImgMF); QMF(i,k)=Qindex(Img,ImgMF); end end [DtLF(1,:) DtLF(2,:); DtMF(1,:) DtMF(2,:); ... MSELF(1,:) MSELF(2,:); MSEMF(1,:) MSEMF(2,:); ... RhoLF(1,:) RhoLF(2,:); RhoMF(1,:) RhoMF(2,:); ... QLF(1,:) QLF(2,:); QMF(1,:) QMF(2,:)]
6 Laplacian spatial filtering
The following MATLAB code is an example of the required implementation.
clear all; close all; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Laplacian Filtering img1 = im2double(imread(’moon.jpg’)); lap = [1 1 1; 1 -8 1; 1 1 1;]; img2 = conv2(img1, lap, ’same’); img3 = [img2 - min(img2(:))] ./ max(img2(:) - min(img2(:))); img4 = img1 - img2; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
img2 = [img2 - min(img2(:))] ./ max(img2(:) - min(img2(:)));
7 Change detection using image subtraction
Scene 1 Scene 2
Difference Mask Detected Motion
Figure 6: Detection of movement by subtracting images. The following MATLAB code is an example of the required implementation. In Fig. 6 we show the resulting images.
clear all; close all; imga = im2double(imread(’scene1.jpg’)); imgb = im2double(imread(’scene2.jpg’)); img1 = imga(10:size(imga,1)-9, 10:size(imgb,2)-9); img2 = imgb(10:size(imga,1)-9, 10:size(imgb,2)-9);
idff = abs(img1 - img2); idx = find(idff>.2); mask = zeros(size(idff)); mask(idx) = 1; red = img1; green = img1; blue = img1; red(idx) = 1; green(idx) = 0; blue(idx) = 0; img3 = cat(3, cat(3,red,green), blue);
figure;
subplot(2,2,1); imshow(img1); title(’Scene 1’); subplot(2,2,2); imshow(img2); title(’Scene 2’); subplot(2,2,3); imshow(mask); title(’Difference Mask’); subplot(2,2,4); imshow(img3); title(’Detected Motion’);
