how can i reduce false features points in an image ?
2 views (last 30 days)
Show older comments
Hello , Ineed help please ,Iam working about fatigue detection and i use this SIFT code to detect facial expressions , the problem here is while i run the code i find false features points in the image ,so who can help me to reduce this points , I modified the threshold but it didnt work
thanks in advance
%%%%%%%%%%%%%%%%%%%%%%
% Scale-Space Extrema Detection
tic
% original sigma and the number of actave can be modified. the larger
% sigma0, the more quickly-smooth images
sigma0=sqrt(2);
octave=3;%6*sigma*k^(octave*level)<=min(m,n)/(2^(octave-2))
level=3;
D=cell(1,octave);
for i=1:octave
D(i)=mat2cell(zeros(row*2^(2-i)+2,colum*2^(2-i)+2,level),row*2^(2-i)+2,colum*2^(2-i)+2,level);
end
% first image in first octave is created by interpolating the original one.
temp_img=kron(I1,ones(2));
temp_img=padarray(temp_img,[1,1],'replicate');
figure;
subplot(2,2,1);
imshow(I1);
% %create the DoG pyramid.
for i=1:octave
temp_D=D{i};
for j=1:level
scale=sigma0*sqrt(2)^(1/level)^((i-1)*level+j);
p=(level)*(i-1);
figure(1);
subplot(octave,level,p+j);
f=fspecial('gaussian',[1,floor(6*scale)],scale);
L1=temp_img;
if(i==1&&j==1)
L2=conv2(temp_img,f,'same');
L2=conv2(L2,f','same');
temp_D(:,:,j)=L2-L1;
imshow(uint8(255 * mat2gray(temp_D(:,:,j))));
L1=L2;
else
L2=conv2(temp_img,f,'same');
L2=conv2(L2,f','same');
temp_D(:,:,j)=L2-L1;
L1=L2;
if(j==level)
temp_img=L1(2:end-1,2:end-1);
end
imshow(uint8(255 * mat2gray(temp_D(:,:,j))));
end
end
D{i}=temp_D;
temp_img=temp_img(1:2:end,1:2:end);
temp_img=padarray(temp_img,[1,1],'both','replicate');
end
toc
% %% Keypoint Localistaion
% % search each pixel in the DoG map to find the extreme point
tic
interval=level-1;
number=0;
for i=2:octave+1
number=number+(2^(i-octave)*colum)*(2*row)*interval;
end
extrema=zeros(1,4*number);
flag=1;
for i=1:octave
[m,n,~]=size(D{i});
m=m-2;
n=n-2;
volume=m*n/(4^(i-1));
for k=2:interval
for j=1:volume
% starter=D{i}(x+1,y+1,k);
x=ceil(j/n);
y=mod(j-1,m)+1;
sub=D{i}(x:x+2,y:y+2,k-1:k+1);
large=max(max(max(sub)));
little=min(min(min(sub)));
if(large==D{i}(x+1,y+1,k))
temp=[i,k,j,1];
extrema(flag:(flag+3))=temp;
flag=flag+4;
end
if(little==D{i}(x+1,y+1,k))
temp=[i,k,j,-1];
extrema(flag:(flag+3))=temp;
flag=flag+4;
end
end
end
end
idx= extrema==0;
extrema(idx)=[];
toc
[m,n]=size(I1);
x=floor((extrema(3:4:end)-1)./(n./(2.^(extrema(1:4:end)-2))))+1;
y=mod((extrema(3:4:end)-1),m./(2.^(extrema(1:4:end)-2)))+1;
ry=y./2.^(octave-1-extrema(1:4:end));
rx=x./2.^(octave-1-extrema(1:4:end));
figure(2)
subplot(1,2,2);
imshow(I1)
hold on
plot(ry,rx,'r+');
%% accurate keypoint localization
%eliminate the point with low contrast or poorly localised on an edge
% x:|,y:-- x is for vertial and y is for horizontal
% value comes from the paper.
tic
threshold=0.09;
r=10;
extr_volume=length(extrema)/4;
[m,n]=size(I1);
secondorder_x=conv2([-1,1;-1,1],[-1,1;-1,1]);
secondorder_y=conv2([-1,-1;1,1],[-1,-1;1,1]);
for i=1:octave
for j=1:level
test=D{i}(:,:,j);
temp=-1./conv2(test,secondorder_y,'same').*conv2(test,[-1,-1;1,1],'same');
D{i}(:,:,j)=temp.*conv2(test',[-1,-1;1,1],'same')*0.5+test;
end
end
for i=1:extr_volume
x=floor((extrema(4*(i-1)+3)-1)/(n/(2^(extrema(4*(i-1)+1)-2))))+3;
y=mod((extrema(4*(i-1)+3)-1),m/(2^(extrema(4*(i-1)+1)-2)))+3;
rx=x+1;
ry=y+1;
rz=extrema(4*(i-1)+2);
z=D{extrema(4*(i-1)+1)}(rx,ry,rz);
if(abs(z)<threshold)
extrema(4*(i-1)+4)=0;
end
end
idx=find(extrema==0);
idx=[idx,idx-1,idx-2,idx-3];
extrema(idx)=[];
extr_volume=length(extrema)/4;
x=floor((extrema(3:4:end)-1)./(n./(2.^(extrema(1:4:end)-2))))+1;
y=mod((extrema(3:4:end)-1),m./(2.^(extrema(1:4:end)-2)))+1;
ry=y./2.^(octave-1-extrema(1:4:end));
rx=x./2.^(octave-1-extrema(1:4:end));
figure(2)
subplot(2,2,3);
imshow(I1)
hold on
plot(ry,rx,'g+');
for i=1:extr_volume
x=floor((extrema(4*(i-1)+3)-1)/(n/(2^(extrema(4*(i-1)+1)-2))))+1;
y=mod((extrema(4*(i-1)+3)-1),m/(2^(extrema(4*(i-1)+1)-2)))+1;
rx=x+1;
ry=y+1;
rz=extrema(4*(i-1)+2);
Dxx=D{extrema(4*(i-1)+1)}(rx-1,ry,rz)+D{extrema(4*(i-1)+1)}(rx+1,ry,rz)-2*D{extrema(4*(i-1)+1)}(rx,ry,rz);
Dyy=D{extrema(4*(i-1)+1)}(rx,ry-1,rz)+D{extrema(4*(i-1)+1)}(rx,ry+1,rz)-2*D{extrema(4*(i-1)+1)}(rx,ry,rz);
Dxy=D{extrema(4*(i-1)+1)}(rx-1,ry-1,rz)+D{extrema(4*(i-1)+1)}(rx+1,ry+1,rz)-D{extrema(4*(i-1)+1)}(rx-1,ry+1,rz)-D{extrema(4*(i-1)+1)}(rx+1,ry-1,rz);
deter=Dxx*Dyy-Dxy*Dxy;
R=(Dxx+Dyy)/deter;
R_threshold=(r+1)^2/r;
if(deter<0||R>R_threshold)
extrema(4*(i-1)+4)=0;
end
end
idx=find(extrema==0);
idx=[idx,idx-1,idx-2,idx-3];
extrema(idx)=[];
extr_volume=length(extrema)/4;
x=floor((extrema(3:4:end)-1)./(n./(2.^(extrema(1:4:end)-2))))+1;
y=mod((extrema(3:4:end)-1),m./(2.^(extrema(1:4:end)-2)))+1;
ry=y./2.^(octave-1-extrema(1:4:end));
rx=x./2.^(octave-1-extrema(1:4:end));
figure(2)
subplot(2,2,4);
imshow(I1)
hold on
plot(ry,rx,'b+');
toc
%Orientation Assignment(Multiple orientations assignment)
tic
kpori=zeros(1,36*extr_volume);
minor=zeros(1,36*extr_volume);
f=1;
flag=1;
for i=1:extr_volume
% search in the certain scale
scale=sigma0*sqrt(2)^(1/level)^((extrema(4*(i-1)+1)-1)*level+(extrema(4*(i-1)+2)));
width=2*round(3*1.5*scale);
count=1;
x=floor((extrema(4*(i-1)+3)-1)/(n/(2^(extrema(4*(i-1)+1)-2))))+1;
y=mod((extrema(4*(i-1)+3)-1),m/(2^(extrema(4*(i-1)+1)-2)))+1;
% make sure the point in the searchable area
if(x>(width/2)&&y>(width/2)&&x<(m/2^(extrema(4*(i-1)+1)-2)-width/2-2)&&y<(n/2^(extrema(4*(i-1)+1)-2)-width/2-2))
rx=x+1;
ry=y+1;
rz=extrema(4*(i-1)+2);
reg_volume=width*width;%3? thereom
% make weight matrix
weight=fspecial('gaussian',width,1.5*scale);
% calculate region pixels' magnitude and region orientation
reg_mag=zeros(1,count);
reg_theta=zeros(1,count);
for l=(rx-width/2):(rx+width/2-1)
for k=(ry-width/2):(ry+width/2-1)
reg_mag(count)=sqrt((D{extrema(4*(i-1)+1)}(l+1,k,rz)-D{extrema(4*(i-1)+1)}(l-1,k,rz))^2+(D{extrema(4*(i-1)+1)}(l,k+1,rz)-D{extrema(4*(i-1)+1)}(l,k-1,rz))^2);
reg_theta(count)=atan2((D{extrema(4*(i-1)+1)}(l,k+1,rz)-D{extrema(4*(i-1)+1)}(l,k-1,rz)),(D{extrema(4*(i-1)+1)}(l+1,k,rz)-D{extrema(4*(i-1)+1)}(l-1,k,rz)))*(180/pi);
count=count+1;
end
end
% make histogram
mag_counts=zeros(1,36);
for x=0:10:359
mag_count=0;
for j=1:reg_volume
c1=-180+x;
c2=-171+x;
if(c1<0||c2<0)
if(abs(reg_theta(j))<abs(c1)&&abs(reg_theta(j))>=abs(c2))
mag_count=mag_count+reg_mag(j)*weight(ceil(j/width),mod(j-1,width)+1);
end
else
if(abs(reg_theta(j)>abs(c1)&&abs(reg_theta(j)<=abs(c2))))
mag_count=mag_count+reg_mag(j)*weight(ceil(j/width),mod(j-1,width)+1);
end
end
end
mag_counts(x/10+1)=mag_count;
end
% find the max histogram bar and the ones higher than 80% max
[maxvm,~]=max(mag_counts);
kori=find(mag_counts>=(0.8*maxvm));
kori=(kori*10+(kori-1)*10)./2-180;
kpori(f:(f+length(kori)-1))=kori;
f=f+length(kori);
temp_extrema=[extrema(4*(i-1)+1),extrema(4*(i-1)+2),extrema(4*(i-1)+3),extrema(4*(i-1)+4)];
temp_extrema=padarray(temp_extrema,[0,length(temp_extrema)*(length(kori)-1)],'post','circular');
long=length(temp_extrema);
minor(flag:flag+long-1)=temp_extrema;
flag=flag+long;
end
end
idx= minor==0;
minor(idx)=[];
extrema=minor;
%delete unsearchable points and add minor orientation points
idx= kpori==0;
kpori(idx)=[];
extr_volume=length(extrema)/4;
toc
%keypoint descriptor
tic
d=4;% In David G. Lowe experiment,divide the area into 4*4.
pixel=4;
feature=zeros(d*d*8,extr_volume);
for i=1:extr_volume
descriptor=zeros(1,d*d*8);% feature dimension is 128=4*4*8;
width=d*pixel;
% x,y centeral point and prepare for location rotation
x=floor((extrema(4*(i-1)+3)-1)/(n/(2^(extrema(4*(i-1)+1)-2))))+1;
y=mod((extrema(4*(i-1)+3)-1),m/(2^(extrema(4*(i-1)+1)-2)))+1;
z=extrema(4*(i-1)+2);
if((m/2^(extrema(4*(i-1)+1)-2)-pixel*d*sqrt(2)/2)>x&&x>(pixel*d/2*sqrt(2))&&(n/2^(extrema(4*(i-1)+1)-2)-pixel*d/2*sqrt(2))>y&&y>(pixel*d/2*sqrt(2)))
sub_x=(x-d*pixel/2+1):(x+d*pixel/2);
sub_y=(y-d*pixel/2+1):(y+d*pixel/2);
sub=zeros(2,length(sub_x)*length(sub_y));
j=1;
for p=1:length(sub_x)
for q=1:length(sub_y)
sub(:,j)=[sub_x(p)-x;sub_y(q)-y];
j=j+1;
end
end
distort=[cos(pi*kpori(i)/180),-sin(pi*kpori(i)/180);sin(pi*kpori(i)/180),cos(pi*kpori(i)/180)];
% accordinate after distort
sub_dis=distort*sub;
fix_sub=ceil(sub_dis);
fix_sub=[fix_sub(1,:)+x;fix_sub(2,:)+y];
patch=zeros(1,width*width);
for p=1:length(fix_sub)
patch(p)=D{extrema(4*(i-1)+1)}(fix_sub(1,p),fix_sub(2,p),z);
end
temp_D=(reshape(patch,[width,width]))';
% create weight matrix.
mag_sub=temp_D;
temp_D=padarray(temp_D,[1,1],'replicate','both');
weight=fspecial('gaussian',width,width/1.5);
mag_sub=weight.*mag_sub;
theta_sub=atan((temp_D(2:end-1,3:1:end)-temp_D(2:end-1,1:1:end-2))./(temp_D(3:1:end,2:1:end-1)-temp_D(1:1:end-2,2:1:end-1)))*(180/pi);
% create orientation histogram
for area=1:d*d
cover=pixel*pixel;
ori=zeros(1,cover);
magcounts=zeros(1,8);
for angle=0:45:359
magcount=0;
for p=1:cover;
x=(floor((p-1)/pixel)+1)+pixel*floor((area-1)/d);
y=mod(p-1,pixel)+1+pixel*(mod(area-1,d));
c1=-180+angle;
c2=-180+45+angle;
if(c1<0||c2<0)
if (abs(theta_sub(x,y))<abs(c1)&&abs(theta_sub(x,y))>=abs(c2))
ori(p)=(c1+c2)/2;
magcount=magcount+mag_sub(x,y);
end
else
if(abs(theta_sub(x,y))>abs(c1)&&abs(theta_sub(x,y))<=abs(c2))
ori(p)=(c1+c2)/2;
magcount=magcount+mag_sub(x,y);
end
end
end
magcounts(angle/45+1)=magcount;
end
descriptor((area-1)*8+1:area*8)=magcounts;
end
descriptor=normr(descriptor);
% cap 0.2
for j=1:numel(descriptor)
if(abs(descriptor(j))>0.2)
descriptor(j)=0.2;
end
end
descriptor=normr(descriptor);
else
continue;
end
feature(:,i)=descriptor';
end
index=find(sum(feature));
feature=feature(:,index);
toc
0 Comments
Answers (0)
See Also
Categories
Find more on Feature Detection and Extraction 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!