Image de-noising using decision tree
Show older comments
%DENOISING OF IMPULSE IMAGES USING DTBTM:
clc;
close all;
x=imread('lena512.bmp');
x =imresize(x,[512 512]);
NOISE_VAR = 0.2;
y=x;
%for z = 1:75
r = randi(255,1,2);
y=pepperOrSalt(y,NOISE_VAR,3,r(1),r(2));
%end
y = double(y);
decision1=0;decision2=0;decision3=0;decision4=0;decision5=0;decision6=0;decision7=0;decision8=0;
f = y;f_m = y;
[R, C] = size(x);
th_ima=20;th_imb=25;th_fma=40;th_fmb=80;th_sma=15;th_smb=60;
for z =1:1
for i = 4:R-1
for j = 4:C-1
tmp =y(i-2:i+2,i-1:i+1,j-2:j+2,j-1:j+1);
w_top_half=[y(i-2,j-2),y(i-2,j-1),y(i-2,j),y(i-2,j+1),y(i-2,j+2),y(i-1,j-2),y(i-1,j-1),y(i-1,j),y(i-1,j+1),y(i-1,j+2),y(i,j-2),y(i,j-1)];
w_bottom_half=[y(i,j+1),y(i,j+2),y(i+1,j-2),y(i+1,j-1),y(i+1,j),y(i+1,j+1),y(i+1,j+2),y(i+2,j-2),y(i+2,j-1),y(i+2,j),y(i+2,j+1),y(i+2,j+2)];
s_top = sort(w_top_half);
s_bottom = sort(w_bottom_half);
s_top_max = s_top(12);
s_top_min = s_top(1);
s_bottom_min = s_bottom(1);
s_bottom_max = s_bottom(12);
top_half_diff = abs(s_top_max-s_top_min);
bottom_half_diff = abs(s_bottom_max-s_bottom_min);
if (top_half_diff>=th_ima || bottom_half_diff>=th_ima)
decision1=1;
else
decision1=0;
end
if(decision1==0)
if(abs(y(i,j)-s_top_max)>=th_imb || abs(y(i,j)-s_top_min)>=th_imb)
im_top_half=1;
else
im_top_half=0;
end
if(abs(y(i,j)-s_bottom_max)>=th_imb || abs(y(i,j)-s_bottom_min)>=th_imb)
im_bottom_half=1;
else
im_bottom_half=0;
end
if(im_top_half==1 || im_bottom_half==1)
decision2 =1;
else
decision2 =0;
end
%else
% y(i,j) = y(i,j);
end
%if(decision2==1)
if(abs(y(i-2,j-2)-y(i,j))>=th_fma || abs(y(i+2,j+2)-y(i,j))>=th_fma || abs(y(i-2,j-2)-y(i+2,j+2))>=th_fma)
fm_e1=0;%no edge
else
fm_e1=1;
end
if(abs(y(i-2,j+2)-y(i,j))>=th_fma || abs(y(i+2,j-2)-y(i,j))>=th_fma || abs(y(i-2,j+2)-y(i+2,j-2))>=th_fma)
fm_e2=0;
else
fm_e2=1;
end
if(abs(y(i-1,j-2)-y(i,j))>=th_fma || abs(y(i-1,j+2)-y(i,j))>=th_fma || abs(y(i-1,j-2)-y(i-1,j+2))>=th_fma)
fm_e3=0;
else
fm_e3=1;
end
if(abs(y(i,j-2)-y(i,j))>=th_fma || abs(y(i,j+2)-y(i,j))>=th_fma || abs(y(i,j-2)-y(i,j+2))>=th_fma)
fm_e4=0;
else
fm_e4=1;
end
if(abs(y(i+1,j-2)-y(i,j))>=th_fma || abs(y(i+1,j+2)-y(i,j))>=th_fma || abs(y(i+1,j-2)-y(i+1,j+2))>=th_fma)
fm_e5=0;
else
fm_e5=1;
end
if(abs(y(i-2,j-1)-y(i,j))>=th_fma || abs(y(i+2,j-1)-y(i,j))>=th_fma || abs(y(i-2,j-1)-y(i+2,j-1))>=th_fma)
fm_e6=0;
else
fm_e6=1;
end
if(abs(y(i-2,j)-y(i,j))>=th_fma || abs(y(i+2,j)-y(i,j))>=th_fma || abs(y(i-2,j)-y(i+2,j))>=th_fma)
fm_e7=0;
else
fm_e7=1;
end
if(abs(y(i-2,j+1)-y(i,j))>=th_fma || abs(y(i+2,j+1)-y(i,j))>=th_fma || abs(y(i-2,j+1)-y(i+2,j+1))>=th_fma)
fm_e8=0;
else
fm_e8=1;
end
if(fm_e1==1 || fm_e2==1 || fm_e3==1 || fm_e4==1 || fm_e5==1 || fm_e6==1 || fm_e7==1 || fm_e8==1)
decision3=0;
else
decision3=1;
end
%else
% y(i,j) = y(i,j);
%end
%if(decision3==1)
s_w=sort(tmp);
w11=s_w(11);
w12=s_w(12);
w_median=s_w(13);
w14=s_w(14);
w15=s_w(15);
w_max=w11+w12+th_sma;
w_min=w14-w15-th_sma;
if(w_max<=w_median+th_smb)
n_max=w_max;
else
n_max=w_median+th_smb;
end
if(w_min>=w_median-th_smb)
n_min=w_min;
else
n_min=w_median-th_smb;
end
if(y(i,j)>=n_max || y(i,j)<=n_min)
decision4=1;
else
decision4=0;
end
%else
% y(i,j) = y(i,j);
%end
if(((decision4==1)||(decision3==1)||((decision2==1)&&(decision1==0))))
d1=abs(y(i-2,j-2)-y(i-1,j+2))+abs(y(i-1,j-2)-y(i,j+2))+abs(y(i,j-2)-y(i+1,j+2))+abs(y(i+1,j-2)-y(i+2,j+2));
d2=abs(y(i-2,j-2)-y(i+2,j-1))+abs(y(i-2,j-1)-y(i+2,j))+abs(y(i-2,j)-y(i+2,j+1))+abs(y(i-2,j+1)-y(i+2,j+2));
d3=(abs(y(i-2,j-1)-y(i+2,j-1)))*2;
d4=(abs(y(i-2,j)-y(i+2,j+1)))*2;
d5=(abs(y(i-2,j+1)-y(i+2,j+1)))*2;
d6=abs(y(i-2,j-1)-y(i+2,j-2))+abs(y(i-2,j)-y(i+2,j-1))+abs(y(i-2,j+1)-y(i+2,j))+abs(y(i-2,j+2)-y(i+2,j+1));
d7=abs(y(i-1,j-2)-y(i-2,j+2))+abs(y(i,j-2)-y(i-1,j+2))+abs(y(i+1,j-2)-y(i,j+2))+abs(y(i+2,j-2)-y(i+1,j+2));
d8=(abs(y(i-1,j-2)-y(i-1,j+2)))*2;
d9=(abs(y(i,j-2)-y(i,j+2)))*2;
d10=(abs(y(i+1,j-2)-y(i+1,j+2)))*2;
d11=(abs(y(i-2,j-2)-y(i+2,j+2)))*2;
d12=(abs(y(i-2,j+2)-y(i+2,j-2)))*2;
if((y(i,j-2)>=w_max || y(i,j-2)<=w_min) && (y(i,j-1)>=w_max ||y(i,j-1)<=w_min) && (y(i,j+1)>=w_max || y(i,j+1)<=w_min )&& (y(i,j+2)>=w_max || y(i,j+2)<=w_min) && (y(i+1,j-2)>=w_max || y(i+1,j-2)<=w_min)&& (y(i+1,j-1)>=w_max || y(i+1,j-1)<=w_min)&& (y(i+1,j)>=w_max || y(i+1,j)<=w_min)&& (y(i+1,j+1)>=w_max || y(i+1,j+1)<=w_min)&& (y(i+1,j+2)>=w_max || y(i+1,j+2)<=w_min)&& (y(i+2,j-2)>=w_max || y(i+2,j-2)<=w_min)&& (y(i+2,j-1)>=w_max || y(i+2,j-1)<=w_min)&& (y(i+2,j)>=w_max || y(i+2,j)<=w_min)&& (y(i+2,j+1)>=w_max || y(i+2,j+1)<=w_min)&& (y(i+2,j+2)>=w_max || y(i+2,j+2)<=w_min))
f(i,j)=((y(i-2,j-2)+y(i-2,j-1)+(y(i-2,j)*2)+y(i-2,j+1)+y(i-2,j+2))+(y(i-1,j-2)+y(i-1,j-1)+(y(i-1,j)*2)+y(i-1,j+1)+y(i-1,j+2)))/8;
else
if(y(i,j-2)>=w_max || y(i,j-2)<=w_min)
d_m=[d2 d3 d4 d5 d6 d8 d10 d11 d12];
elseif(y(i,j+2)>=w_max ||y(i,j+2)<=w_min)
d_m=[d2 d3 d4 d5 d6 d8 d10 d11 d12];
elseif(y(i+1,j-2)>=w_max || y(i+1,j-2)<=w_min)
d_m=[d2 d3 d4 d5 d6 d8 d9 d11 d12];
elseif(y(i+1,j+2)>=w_max || y(i+1,j+2)<=w_min)
d_m=[d2 d3 d4 d5 d6 d8 d9 d11 d12];
elseif(y(i+2,j-2)>=w_max || y(i+2,j-2)<=w_min)
d_m=[d1 d2 d3 d4 d5 d8 d9 d10 d11];
elseif(y(i+2,j-1)>=w_max || y(i+2,j-1)<=w_min)
d_m=[d1 d2 d4 d5 d8 d9 d10 d11 d12];
elseif(y(i+2,j)>=w_max || y(i+2,j)<=w_min)
d_m=[d1 d3 d5 d7 d8 d9 d10 d11];
elseif(y(i+2,j+1)>=w_max || y(i+2,j+1)<=w_min)
d_m=[ d1 d3 d4 d7 d8 d9 d10 d11 d12];
elseif(y(i+2,j+2)>=w_max || y(i+2,j+2)<=w_min)
d_m=[d3 d4 d5 d6 d7 d8 d9 d10 d12];
else
d_m=[d1 d2 d3 d4 d5 d6 d7 d8 d9 d10 d11 d12];
end
s_d=sort(d_m);
d_min=s_d(1);
if(d_min==d1)
y(i,j)=(y(i-2,j-2)+y(i-1,j+2)+y(i-1,j-2)+y(i,j+2)+y(i,j-2)+y(i+1,j+2)+y(i+1,j-2)+y(i+2,j+2))/8;
elseif(d_min==d2)
y(i,j)=(y(i-2,j-2)+y(i+2,j-1)+y(i-2,j-1)+y(i+2,j)+y(i-2,j)+y(i+2,j+1)+y(i-2,j+1)+y(i+2,j+2))/8;
elseif(d_min==d3)
y(i,j)=(y(i-2,j-1)+y(i+2,j-1))/2;
elseif(d_min==d4)
y(i,j)=(y(i-2,j)+y(i+2,j+1))/2;
elseif(d_min==d5)
y(i,j)=(y(i-2,j+1)+y(i+2,j+1))/2;
elseif(d_min==d6)
y(i,j)=(y(i-2,j-1)+y(i+2,j-2)+y(i-2,j)+y(i+2,j-1)+y(i-2,j+1)+y(i+2,j)+y(i-2,j+2)+y(i+2,j+1))/8;
elseif(d_min==d7)
y(i,j)=(y(i-1,j-2)+y(i-2,j+2)+y(i,j-2)+y(i-1,j+2)+y(i+1,j-2)+y(i,j+2)+y(i+2,j-2)+y(i+1,j+2))/8;
elseif(d_min==d8)
y(i,j)=(y(i-1,j-2)+y(i-1,j+2))/2;
elseif(d_min==d9)
y(i,j)=(y(i,j-2)+y(i,j+2))/2;
elseif(d_min==d10)
y(i,j)=(y(i+1,j-2)+y(i+1,j+2))/2;
elseif(d_min==d11)
y(i,j)=(y(i-2,j-2)+y(i+2,j+2))/2;
elseif(d_min==d12)
y(i,j)=(y(i-2,j+2)+y(i+2,j-2))/2;
end
end
m=[y(i,j) y(i-1,j) y(i,j-1) y(i,j+1) y(i+1,j)];
y(i,j) = median(m);
%else
% y(i,j) = y(i,j);
end
%if(f_m(i,j)>=n_max||f_m(i,j)<=n_min)
% if((f_m(i-2,j-2)<n_max)&&(f_m(i-2,j-2)>n_min))
% f_m(i,j)=f_m(i-2,j-2);
% elseif((f_m(i-2,j-1)<n_max)&&(f_m(i-2,j-1)>n_min))
% f_m(i,j)=f_m(i-2,j-1);
% elseif((f_m(i-2,j+1)<n_max)&&(f_m(i-2,j+1)>n_min))
% f_m(i,j)=f_m(i-2,j+1);
% elseif((f_m(i-2,j+2)<n_max)&&(f_m(i-2,j+2)>n_min))
% f_m(i,j)=f_m(i-2,j+2);
% elseif((f_m(i-1,j-2)<n_max)&&(f_m(i-1,j-2)>n_min))
% f_m(i,j)=f_m(i-1,j-2);
%elseif((f_m(i-1,j-1)<n_max)&&(f_m(i-1,j-1)>n_min))
% f_m(i,j)=f_m(i-1,j-1);
%elseif((f_m(i-1,j)<n_max)&&(f_m(i-1,j)>n_min))
% f_m(i,j)=f_m(i-1,j);
% elseif((f_m(i-1,j+1)<n_max)&&(f_m(i-1,j+1)>n_min))
% f_m(i,j)=f_m(i-1,j+1);
%elseif((f_m(i-1,j+2)<n_max)&&(f_m(i-1,j+2)>n_min))
% f_m(i,j)=f_m(i-1,j+2);
%elseif((f_m(i,j-2)<n_max)&&(f_m(i,j-2)>n_min))
% f_m(i,j)=f_m(i,j-2);
% elseif((f_m(i,j-1)<n_max)&&(f_m(i,j-1)>n_min))
% f_m(i,j)=f_m(i,j-1);
% elseif((f_m(i,j+1)<n_max)&&(f_m(i,j+1)>n_min))
% f_m(i,j)=f_m(i,j+1);
% elseif((f_m(i,j+2)<n_max)&&(f_m(i,j+2)>n_min))
% f_m(i,j)=f_m(i,j+2);
%elseif((f_m(i+1,j-2)<n_max)&&(f_m(i+1,j-2)>n_min))
% f_m(i,j)=f_m(i+1,j-2);
%elseif((f_m(i+1,j-2)<n_max)&&(f_m(i+1,j-2)>n_min))
% f_m(i,j)=f_m(i+1,j-2);
%elseif((f_m(i+1,j-1)<n_max)&&(f_m(i+1,j-1)>n_min))
% f_m(i,j)=f_m(i+1,j-1);
%elseif((f_m(i+1,j)<n_max)&&(f_m(i+1,j)>n_min))
% f_m(i,j)=f_m(i+1,j);
%elseif((f_m(i+1,j+1)<n_max)&&(f_m(i+1,j+1)>n_min))
% f_m(i,j)=f_m(i+1,j+1);
%elseif((f_m(i+1,j+2)<n_max)&&(f_m(i+1,j+2)>n_min))
% f_m(i,j)=f_m(i+1,j+2);
% elseif((f_m(i+2,j-2)<n_max)&&(f_m(i+2,j-2)>n_min))
% f_m(i,j)=f_m(i+2,j-2);
% elseif((f_m(i+2,j-1)<n_max)&&(f_m(i+2,j-1)>n_min))
% f_m(i,j)=f_m(i+2,j-1);
% elseif((f_m(i+2,j)<n_max)&&(f_m(i+2,j)>n_min))
% f_m(i,j)=f_m(i+2,j);
% elseif((f_m(i+2,j+1)<n_max)&&(f_m(i+2,j+1)>n_min))
% f_m(i,j)=f_m(i+2,j+1);
% elseif((f_m(i+2,j+2)<n_max)&&(f_m(i+2,j+2)>n_min))
% f_m(i,j)=f_m(i+2,j+2);
%else
% f_m(i,j)=(f_m(i-1,j-1)+(f_m(i-1,j)*2)+f_m(i-1,j+1))/4;
f_m(i,j)=(f_m(i-2,j-2)+f_m(i-2,j-1)+(f_m(i-2,j)*2)+f_m(i-2,j+1)+f_m(i-2,j+2)+f_m(i-1,j-2)+f_m(i-1,j-1)+(f_m(i-1,j)*2)+f_m(i-1,j+1)+f_m(i-1,j+2))/8;
%end
%end
end
end
end
figure(1);imshow(x);figure(2);imshow(f_m,[]);figure(3);imshow(y,[]);
squaredErrorImage = (double(x) - (y)) .^ 2;
mse = sum(sum(squaredErrorImage)) / (R * C);
PSNR = 10 * log10( 255^2 / mse);
r1=rand_index(y);
v=uint8(y);
u=uint8(f_m);
gce1 = global_consistancy_error(y);
jaccardIndex_ac = sum(x(:) & v(:)) / sum(x(:) | v(:));
message = sprintf('The mean square error for denoised image is %.2f.\nThe PSNR = %.2f.\nThe Rand_index = %.2f.\nThe global_consistancy_error = %.2f.\nThe jaccard coefficient = %.2f.\nThe jaccard distance = %.2f.', mse, PSNR, r1,gce1,jaccardIndex_ac,1-jaccardIndex_ac);
msgbox(message);
imwrite(v,'G:\lena_gray.bmp');
squaredErrorImage1 = (double(x) - (f_m)) .^ 2;
mse1 = sum(sum(squaredErrorImage1)) / (R * C);
PSNR1 = 10 * log10( 255^2 / mse1);
r2=rand_index(f_m);
gce2 = global_consistancy_error(f_m);
jaccardIndex_ac2 = sum(x(:) & u(:)) / sum(x(:) | u(:));
message = sprintf('The mean square error for noisy image is %.2f.\nThe PSNR = %.2f.\nThe Rand_index = %.2f.\nThe global_consistancy_error = %.2f.\nThe jaccard coefficient = %.2f.\nThe jaccard distance = %.2f.', mse1, PSNR1, r2, gce2,jaccardIndex_ac2,1-jaccardIndex_ac2);
msgbox(message);
x1=logical(x);
y1=logical(y);
n=jaccard_coefficient(x1,y1);
i have to denoise image but iam getting error like
Index in position 3 exceeds array bounds (must not exceed 1).
Error in mycode (line 22)
tmp =y(i-2:i+2,i-1:i+1,j-2:j+2,j-1:j+1);
>>
here iam using 5 by 5 matrix
i need help
Accepted Answer
Image Analyst
on 3 May 2020
You have the wrong starting and stopping indexes for a 5x5 scanning window:
for i = 4:R-1
for j = 4:C-1
You need to start at the 3rd index and go until the 3rd one to the end:
windowWidth = 5;
halfWindowWidth = ceil(windowWidth / 2); % 3
for i = halfWindowWidth : R - halfWindowWidth + 1
for j = halfWindowWidth : C - halfWindowWidth + 1
4 Comments
krishna kotha
on 4 May 2020
Edited: krishna kotha
on 4 May 2020
Image Analyst
on 4 May 2020
Well, what is j-2 and j+2? When you see a message like that you need to look at your variables and figure out why they're wrong. This is just basic debugging. It must be an integer between 1 and the length of the third dimension, which is 1. So j-2 and j+2 cannot be 1. If j is 3, then j-2 can be 1, but j+2 will be 5 which is larger than 1. I have no idea why your third dimension is only 1 instead of something larger -- you'll have to look into that yourself. I can't delve into hundreds of lines of uncommented code. With no comments, no one knows what you want or what you're doing.
krishna kotha
on 7 May 2020
Image Analyst
on 7 May 2020
What does a 2-D 5x5 mask mean when you're trying to scan a 4-D image? Maybe you should extract the 2-D image out of the 4-D hypervolume.
More Answers (0)
Categories
Find more on Correlation and Convolution in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
