Scale-Invariant Feature Detection and Matching Overview The goal of this assignm
ID: 3599813 • Letter: S
Question
Scale-Invariant Feature Detection and Matching
Overview
The goal of this assignment is to design a scale-invariant feature detection and matching algorithm using techniques described in lecture slides and Szeliski Chapter 4.1. For this purpose, you will implement a version of the well-known SIFT pipeline. Your implementation will first find keypoints (or corners) on a pair of images and then match them. You can find test images and the starter code here.
There is also an additional written part of the assignment on PCA.
Details
You will implement the three major steps of a feature matching algorithm:
Keypoint detection in detect_features.m. This function is a state of the art scale-invariant keypoint detector that uses differences of Gaussians (DoG) and already done for you. You will use this function but also implement a detect_corners.m function as a baseline. The detect_corners.mfunction will implement the uniscale Harris corner detector (see lecture slides and Szeliski 4.1.1).
Local feature description in SIFTDescriptor.m (see see lecture slides and Szeliski 4.1.2)
Feature matching in SIFTSimpleMatcher.m (see see lecture slides and Szeliski 4.1.3)
1. Keypoint detection (detect_features.m and detect_corners.m)
The detect_features.m function (already implemented for you in the starter code) is a scale-invariant keypoint detector that makes use of differences of Gaussians (DoG). It outputs a set of keypoints each associated with a "best" scale (in the Gaussian pyramid which is also an output of the function). The output of this function is then given in the next step as input to the SIFT descriptor. Note that a Gaussian pyramid is a set of Gaussian filtered and downsampled versions of a given image at increasing levels of scale.
Although you don't need to do anything about the detect_features.m function, except to understand how to use it, you will implement aside the Harris corner detector (detect_corners.m) as a baseline in the way described in the lecture slides and Szeliski 4.1.1. Note that detect_features.m function is a scale-invariant keypoint detector whereas your Harris corner function detect_corners.m is not expected to be scale invariant, hence simpler to implement. A few recommendations for the implementation: Use color and/or grayscale gradient. To compute gradient, you can use the Prewitt or Sobel masks available in Matlab. While computing the correlation matrix, use Gaussian weighting over a neighborhood (try 3x3, 5x5, 7x7). You do not need to worry about scale invariance for your baseline Harris corner detector. So you can process all pixels at the given scale of the images. To avoid spurious keypoint detections near the boundaries, simply suppress the gradients / corners near the edges of the image. Apply a non-maximum suppression procedure once you've thresholded the cornerness score. Here are some helpful functions: BWLABEL and the newer BWCONNCOMP will find connected components in thresholded binary image. You could, for instance, take the maximum value within each component. COLFILT can be used to run a max() operator on each sliding window. You could use this to ensure that every interest point is at a local maximum of cornerness.You may try different options as corner response function (given in the lecture slides).
You will implement detect_corners.m function from scratch. We strongly recommend you to make the two functions detect_features.m and detect_corners.m compatible in terms of output and input so that the matching pipeline would work seamlessly for both cases. This is important since you can then directly use the SIFTDescriptor.m function that you'll implement in the next step in order to find SIFT descriptors on the Harris corners. Hence your Harris corner detector will generate keypoints just like detect_features.m but at only one fixed scale, that is, the scale of the given images. This means that although your detect_corners.m function won't deal with scale, it will still produce a scale value (fixed as the given scale) and a Gaussian image pyramid (though only one level, i.e., the given image, will be relevant and the others will be redundant).
2. Local feature description (SIFTDescriptor.m)
You will implement the SIFT features as described in the lecture slides and Szeliski 4.1.2. Edit the provided SIFTDescriptor.m to produce a SIFT keypoint descriptor for each keypoint. See also Lowe's related paper provided. Note that the keypoint location and scale are provided for you, as is the Gaussian pyramid by the already implemented detect_features.m function. Run the provided EvaluateSIFTDescriptor.m to check and verify your implementation.
If you want to get your matching pipeline working quickly (and maybe to help debug the other algorithm stages), you might want to start with normalized patches as your local feature, that is, your feature vector for each interest point will simply be composed of image intensities in the neighborhood (in this case intensity normalization can help by subtracting from each pixel's intensity the mean intensity of the neighborhood and then dividing the result by the variance of the intensities in the patch).
3. Feature matching (SIFTSimpleMatcher.m)
At this step, you will match the keypoints (or corners) detected on a given pair of images, using detect_features.m (or detect_corners.m) and SIFTDescriptor.m methods that you have already implemented. Edit SIFTSimpleMatcher.m to calculate the Euclidean distance (or simply SSD, sum of squared distances) between a given SIFT descriptor from image 1 and all SIFT descriptors from image 2. Then use this to determine if there is a "good" match via ratio test: if the distance to the closest vector is significantly smaller (by a factor which is given) than the distance to the second-closest, we call it a match. You're free to try and use other thresholding and test schemes for matching if you get better results (you have to explain it clearly in your report).
The output of the function is an array where each row holds the indices of one pair of matching descriptors.
Run the provided EvaluateSIFTMatcher.m to check your implementation.
Using the Starter Code:
The starter code includes detect_features.m (already complete) along with SIFTDescriptor.m and SIFTSimpleMatcher.m (to be completed by you) as well as some visualization and evaluation functions to test your implementations. There exists also a TestMatcher.m function that you can use to run the whole matching pipeline. You will test your implementation on a number of image pairs provided, where each image pair poses a different type of challenge for the matching problem. Note that EvaluateSIFTMatcher.m and EvaluateSIFTDescriptor.m verify your SIFT description and matching implementations based on some reference ground-truth data. Make sure that your implementations pass the evaluation step.
Potentially useful MATLAB functions: imfilter(), fspecial(), bwconncomp(), colfilt(), bwlabel(), sort().
Explanation / Answer
%% Important Variables
% a : Input image
% kpmag : keypoints magnitude
% kpori : keypoints orientation
% kpd : key point descriptors
% kp : keypoints
% kpl : keypoint locations
% Extension of 2 for above variable indicates,it is used for second image for matching
% mp : matching percentage
clc;
close all;
clear all;
fprintf('Some of the images available for checking are: ');
fprintf('1.lena.jpg(512x512) 2.lena1.jpg(256x256) 3.peppers256.png(256x256) 4.testfi.png(256x256) 5.baboon.bmp(256x256) 6.testfile.jpg(color)(1920x1200) ');
image=input('Choose and enter the image name from above : ','s');
a=imread(image);
imshow(a);
title('Selected image');
[row,col,plane]=size(a);
if plane==3
a=rgb2gray(a);
end
a=im2double(a);
original=a;
store1=[];
store2=[];
store3=[];
tic
%% 1st octave generation
k2=0;
[m,n]=size(a);
a(m:m+6,n:n+6)=0;
clear c;
for k1=0:3
k=sqrt(2);
sigma=(k^(k1+(2*k2)))*1.6;
for x=-3:3
for y=-3:3
h(x+4,y+4)=(1/((2*pi)*((k*sigma)*(k*sigma))))*exp(-((x*x)+(y*y))/(2*(k*k)*(sigma*sigma)));
end
end
for i=1:m
for j=1:n
t=a(i:i+6,j:j+6)'.*h;
c(i,j)=sum(sum(t));
end
end
store1=[store1 c];
end
clear a;
a=imresize(original,1/2);
%% 2nd level pyramid generation
k2=1;
[m,n]=size(a);
a(m:m+6,n:n+6)=0;
clear c;
for k1=0:3
k=sqrt(2);
sigma=(k^(k1+(2*k2)))*1.6;
for x=-3:3
for y=-3:3
h(x+4,y+4)=(1/((2*pi)*((k*sigma)*(k*sigma))))*exp(-((x*x)+(y*y))/(2*(k*k)*(sigma*sigma)));
end
end
for i=1:m
for j=1:n
t=a(i:i+6,j:j+6)'.*h;
c(i,j)=sum(sum(t));
end
end
store2=[store2 c];
end
clear a;
a=imresize(original,1/4);
%% 3rd level pyramid generation
k2=2;
[m,n]=size(a);
a(m:m+6,n:n+6)=0;
clear c;
for k1=0:3
k=sqrt(2);
sigma=(k^(k1+(2*k2)))*1.6;
for x=-3:3
for y=-3:3
h(x+4,y+4)=(1/((2*pi)*((k*sigma)*(k*sigma))))*exp(-((x*x)+(y*y))/(2*(k*k)*(sigma*sigma)));
end
end
for i=1:m
for j=1:n
t=a(i:i+6,j:j+6)'.*h;
c(i,j)=sum(sum(t));
end
end
store3=[store3 c];
end
%% Obtaining key point from the image
i1=store1(1:row,1:col)-store1(1:row,col+1:2*col);
i2=store1(1:row,col+1:2*col)-store1(1:row,2*col+1:3*col);
i3=store1(1:row,2*col+1:3*col)-store1(1:row,3*col+1:4*col);
[m,n]=size(i2);
kp=[];
kpl=[];
for i=2:m-1
for j=2:n-1
x=i1(i-1:i+1,j-1:j+1);
y=i2(i-1:i+1,j-1:j+1);
z=i3(i-1:i+1,j-1:j+1);
y(1:4)=y(1:4);
y(5:8)=y(6:9);
mx=max(max(x));
mz=max(max(z));
mix=min(min(x));
miz=min(min(z));
my=max(max(y));
miy=min(min(y));
if (i2(i,j)>mz && i2(i,j)>my) || (i2(i,j)<miz && i2(i,j)<miy)
kp=[kp i2(i,j)];
kpl=[kpl i j];
end
end
end
%% Key points plotting on to the image
for i=1:2:length(kpl);
k1=kpl(i);
j1=kpl(i+1);
i2(k1,j1)=1;
end
figure, imshow(i2);
title('Image with key points mapped onto it');
%% Magnitude and orientation assignment to the key points
for i=1:m-1
for j=1:n-1
mag(i,j)=sqrt(((i2(i+1,j)-i2(i,j))^2)+((i2(i,j+1)-i2(i,j))^2));
oric(i,j)=atan2(((i2(i+1,j)-i2(i,j))),(i2(i,j+1)-i2(i,j)))*(180/pi);
end
end
%% Forming key point neighbourhooods
kpmag=[];
kpori=[];
for x1=1:2:length(kpl)
k1=kpl(x1);
j1=kpl(x1+1);
if k1 > 2 && j1 > 2 && k1 < m-2 && j1 < n-2
p1=mag(k1-2:k1+2,j1-2:j1+2);
q1=oric(k1-2:k1+2,j1-2:j1+2);
else
continue;
end
%% Finding orientation and magnitude for the key point
[m1,n1]=size(p1);
magcounts=[];
for x=0:10:359
magcount=0;
for i=1:m1
for j=1:n1
ch1=-180+x;
ch2=-171+x;
if ch1<0 || ch2<0
if abs(q1(i,j))<abs(ch1) && abs(q1(i,j))>=abs(ch2)
ori(i,j)=(ch1+ch2+1)/2;
magcount=magcount+p1(i,j);
end
else
if abs(q1(i,j))>abs(ch1) && abs(q1(i,j))<=abs(ch2)
ori(i,j)=(ch1+ch2+1)/2;
magcount=magcount+p1(i,j);
end
end
end
end
magcounts=[magcounts magcount];
end
[maxvm maxvp]=max(magcounts);
kmag=maxvm;
kori=(((maxvp*10)+((maxvp-1)*10))/2)-180;
kpmag=[kpmag kmag];
kpori=[kpori kori];
% maxstore=[];
% for i=1:length(magcounts)
% if magcounts(i)>=0.8*maxvm
% maxstore=[maxstore magcounts(i) i];
% end
% end
%
% if maxstore > 2
% kmag=maxstore(1:2:length(maxstore));
% maxvp1=maxstore(2:2:length(maxstore));
% temp=(countl((2*maxvp1)-1)+countl(2*maxvp1)+1)/2;
% kori=temp;
% end
end
%% Forming key point Descriptors
kpd=[];
%% Forming key point neighbourhooods
for x1=1:2:length(kpl)
k1=kpl(x1);
j1=kpl(x1+1);
if k1 > 7 && j1 > 7 && k1 < m-8 && j1 < n-8
p2=mag(k1-7:k1+8,j1-7:j1+8);
q2=oric(k1-7:k1+8,j1-7:j1+8);
else
continue;
end
kpmagd=[];
kporid=[];
%% Dividing into 4x4 blocks
for k1=1:4
for j1=1:4
p1=p2(1+(k1-1)*4:k1*4,1+(j1-1)*4:j1*4);
q1=q2(1+(k1-1)*4:k1*4,1+(j1-1)*4:j1*4);
%% Finding orientation and magnitude for the key point
[m1,n1]=size(p1);
magcounts=[];
for x=0:45:359
magcount=0;
for i=1:m1
for j=1:n1
ch1=-180+x;
ch2=-180+45+x;
if ch1<0 || ch2<0
if abs(q1(i,j))<abs(ch1) && abs(q1(i,j))>=abs(ch2)
ori(i,j)=(ch1+ch2+1)/2;
magcount=magcount+p1(i,j);
end
else
if abs(q1(i,j))>abs(ch1) && abs(q1(i,j))<=abs(ch2)
ori(i,j)=(ch1+ch2+1)/2;
magcount=magcount+p1(i,j);
end
end
end
end
magcounts=[magcounts magcount];
end
kpmagd=[kpmagd magcounts];
end
end
kpd=[kpd kpmagd];
end
fprintf(' Time taken for calculating the SIFT keys and their desccriptors is :%f ',toc);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 2ng image for comparision
a=original;
fprintf('Calculation of key points has been finished for the given image ');
fprintf('If you want to check the robustness of the program by changing some of the properties.. you can check them here ');
fprintf('Note: Run only one condition at a time Select ');
fprintf('1.Change intensity 2.Change scale i.e., sigma 3.Rotate image(Select only for 2,3,4,5 images mentioned above) ');
ev=input(' Enter the choice ');
if ev==1
disp('Select');
disp('1.Increase intensity 2.Decrease Intensity ');
eis=input('Enter choice ');
if eis==1
p1=input('Enter the level of increase in intensity ');
a=a*p1;
figure, imshow(a);
title('Image after changing the intensity');
sigmae=1.6;
elseif eis==2
p1=input('Enter the level of decrease in intensity ');
a=a/p1;
figure, imshow(a);
title('Image after changing the intensity');
sigmae=1.6;
else
disp('Wrong choice');
end
elseif ev==2
sigmae=input('Enter sigma value :');
elseif ev==3
thetai=input('Enter the angle at which image is to be rotated ');
theta=thetai*(pi/180);
b=imrotate(a,thetai);
fprintf('Rotated image is ');
figure, imshow(b);
title('Rotate image')
[m,n]=size(original);
old=[];
map=[];
for i=1:m
for j=1:n
i1=(cos(theta)*(i-m/2))-(sin(theta)*(j-n/2))+m/2;
j1=(sin(theta)*(i-m/2))+(cos(theta)*(j-n/2))+n/2;
old((i-1)*m+j,1:2)=[i j];
map((i-1)*m+j,1:2)=[round(i1) round(j1)];
end
end
[m1,n1]=size(b);
marginr=floor((m1-m)/2);
marginc=floor((n1-n)/2);
if thetai==90
for i=1:m
for j=1:n
new((i-1)*m+j,1:2)=[map((i-1)*m+j,1)+marginr+1 map((i-1)*m+j,2)+marginc];
end
end
elseif thetai==180
for i=1:m
for j=1:n
new((i-1)*m+j,1:2)=[map((i-1)*m+j,1)+marginr+1 map((i-1)*m+j,2)+marginc+1];
end
end
else
for i=1:m
for j=1:n
new((i-1)*m+j,1:2)=[map((i-1)*m+j,1)+marginr+1 map((i-1)*m+j,2)+marginc+1];
end
end
end
clear a;
for i=1:m
for j=1:n
a(i,j)=b(new((i-1)*m+j,1),new((i-1)*m+j,2));
end
end
figure, imshow(a);
title('straitened version of rotated image using coordinate mapping');
sigmae=1.6;
else
disp('entered wrong choice');
disp('No changes in parameters is done');
sigmae=1.6;
end
[row,col]=size(a);
store1=[];
store2=[];
store3=[];
tic
%% 1st octave generation
k2=0;
[m,n]=size(a);
a(m:m+6,n:n+6)=0;
clear c;
for k1=0:3
k=sqrt(2);
sigma=(k^(k1+(2*k2)))*sigmae;
for x=-3:3
for y=-3:3
h(x+4,y+4)=(1/((2*pi)*((k*sigma)*(k*sigma))))*exp(-((x*x)+(y*y))/(2*(k*k)*(sigma*sigma)));
end
end
for i=1:m
for j=1:n
t=a(i:i+6,j:j+6)'.*h;
c(i,j)=sum(sum(t));
end
end
store1=[store1 c];
end
clear a;
a=imresize(original,1/2);
%% 2nd level pyramid generation
k2=1;
[m,n]=size(a);
a(m:m+6,n:n+6)=0;
a=im2double(a);
clear c;
for k1=0:3
k=sqrt(2);
sigma=(k^(k1+(2*k2)))*sigmae;
for x=-3:3
for y=-3:3
h(x+4,y+4)=(1/((2*pi)*((k*sigma)*(k*sigma))))*exp(-((x*x)+(y*y))/(2*(k*k)*(sigma*sigma)));
end
end
for i=1:m
for j=1:n
t=a(i:i+6,j:j+6)'.*h;
c(i,j)=sum(sum(t));
end
end
store2=[store2 c];
end
clear a;
a=imresize(original,1/4);
%% 3rd level pyramid generation
k2=2;
[m,n]=size(a);
a(m:m+6,n:n+6)=0;
a=im2double(a);
clear c;
for k1=0:3
k=sqrt(2);
sigma=(k^(k1+(2*k2)))*sigmae;
for x=-3:3
for y=-3:3
h(x+4,y+4)=(1/((2*pi)*((k*sigma)*(k*sigma))))*exp(-((x*x)+(y*y))/(2*(k*k)*(sigma*sigma)));
end
end
for i=1:m
for j=1:n
t=a(i:i+6,j:j+6)'.*h;
c(i,j)=sum(sum(t));
end
end
store3=[store3 c];
end
%% Obtaining key point from the image
i1=store1(1:row,1:col)-store1(1:row,col+1:2*col);
i2=store1(1:row,col+1:2*col)-store1(1:row,2*col+1:3*col);
i3=store1(1:row,2*col+1:3*col)-store1(1:row,3*col+1:4*col);
if ev==2
figure, imshow(i2);
title('Image after change in the sigma value')
end
[m,n]=size(i2);
kp2=[];
kpl2=[];
for i=2:m-1
for j=2:n-1
y=[];
x=i1(i-1:i+1,j-1:j+1);
y=i2(i-1:i+1,j-1:j+1);
z=i3(i-1:i+1,j-1:j+1);
y(1:4)=y(1:4);
y(5:8)=y(6:9);
mx=max(max(x));
mz=max(max(z));
mix=min(min(x));
miz=min(min(z));
my=max(max(y));
miy=min(min(y));
if (i2(i,j)>mz && i2(i,j)>my) || (i2(i,j)<miz && i2(i,j)<miy)
kp2=[kp2 i2(i,j)];
kpl2=[kpl2 i j];
end
end
end
%% Key points plotting on to the image
for i=1:2:length(kpl2);
i1=kpl2(i);
j1=kpl2(i+1);
i2(i1,j1)=1;
end
figure, imshow(i2);
title('2nd image with key points mapped onto it');
%% Magnitude and orientation assignment to the key points
for i=1:m-1
for j=1:n-1
mag(i,j)=sqrt(((i2(i+1,j)-i2(i,j))^2)+((i2(i,j+1)-i2(i,j))^2));
oric(i,j)=atan2(((i2(i+1,j)-i2(i,j))),(i2(i,j+1)-i2(i,j)))*(180/pi);
end
end
%% Forming key point neighbourhooods
kpmag2=[];
kpori2=[];
for x1=1:2:length(kpl2)
k1=kpl2(x1);
j1=kpl2(x1+1);
if k1 > 2 && j1 > 2 && k1 < m-2 && j1 < n-2
p1=mag(k1-2:k1+2,j1-2:j1+2);
q1=oric(k1-2:k1+2,j1-2:j1+2);
else
continue;
end
%% Finding orientation and magnitude for the key point
[m1,n1]=size(p1);
magcounts=[];
for x=0:10:359
magcount=0;
for i=1:m1
for j=1:n1
ch1=-180+x;
ch2=-171+x;
if ch1<0 || ch2<0
if abs(q1(i,j))<abs(ch1) && abs(q1(i,j))>=abs(ch2)
ori(i,j)=(ch1+ch2+1)/2;
magcount=magcount+p1(i,j);
end
else
if abs(q1(i,j))>abs(ch1) && abs(q1(i,j))<=abs(ch2)
ori(i,j)=(ch1+ch2+1)/2;
magcount=magcount+p1(i,j);
end
end
end
end
magcounts=[magcounts magcount];
end
[maxvm maxvp]=max(magcounts);
kmag=maxvm;
kori=(((maxvp*10)+((maxvp-1)*10))/2)-180;
kpmag2=[kpmag2 kmag];
kpori2=[kpori2 kori];
% maxstore=[];
% for i=1:length(magcounts)
% if magcounts(i)>=0.8*maxvm
% maxstore=[maxstore magcounts(i) i];
% end
% end
%
% if maxstore > 2
% kmag=maxstore(1:2:length(maxstore));
% maxvp1=maxstore(2:2:length(maxstore));
% temp=(countl((2*maxvp1)-1)+countl(2*maxvp1)+1)/2;
% kori=temp;
% end
end
%% Forming key point Descriptors
kpd2=[];
for x1=1:2:length(kpl2)
k1=kpl2(x1);
j1=kpl2(x1+1);
if k1 > 7 && j1 > 7 && k1 < m-8 && j1 < n-8
p2=mag(k1-7:k1+8,j1-7:j1+8);
q2=oric(k1-7:k1+8,j1-7:j1+8);
else
continue;
end
kpmagd=[];
%% Dividing into 4x4 blocks
for k1=1:4
for j1=1:4
p1=p2(1+(k1-1)*4:k1*4,1+(j1-1)*4:j1*4);
q1=q2(1+(k1-1)*4:k1*4,1+(j1-1)*4:j1*4);
%% Finding orientation and magnitude for the key point
[m1,n1]=size(p1);
magcounts=[];
for x=0:45:359
magcount=0;
for i=1:m1
for j=1:n1
ch1=-180+x;
ch2=-180+45+x;
if ch1<0 || ch2<0
if abs(q1(i,j))<abs(ch1) && abs(q1(i,j))>=abs(ch2)
ori(i,j)=(ch1+ch2+1)/2;
magcount=magcount+p1(i,j);
end
else
if abs(q1(i,j))>abs(ch1) && abs(q1(i,j))<=abs(ch2)
ori(i,j)=(ch1+ch2+1)/2;
magcount=magcount+p1(i,j);
end
end
end
end
magcounts=[magcounts magcount];
end
kpmagd=[kpmagd magcounts];
end
end
kpd2=[kpd2 kpmagd];
end
fprintf(' Time taken for calculating the SIFT keys and their desccriptors for 2nd image is :%f ',toc);
%% Two images key point comparision
tic
count=0;
for i=1:2:length(kpl)
for j=1:2:length(kpl2)
if (kpl(i)==kpl2(j)) && (kpl(i+1)==kpl2(j+1))
count=count+1;
break;
end
end
end
mp=(count/length(kp))*100;
fprintf('Time taken for calculating the matching percentage is :%f ',toc);
fprintf('Matching percentage between 2 images by key point location is :%f ',mp);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% END OF THE PROGRAM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%