Real time AR tag detection.

algorithms/AR tag detection/scripts/ARTagDetection.m

@@ -0,0 +1,104 @@
+clear;
+clc;
+
+%% INPUT IMAGE AND EDGES 
+
+N = 7; % no of cubes in a row  in artag 
+% cam = webcam;
+% Ir = snapshot(cam);
+% preview(cam);
+% figure, imshow(Ir), title('Input image');
+% clear('cam');
+
+Ir = imread('sample1.jpg');
+I1 = rgb2gray(Ir);
+figure, imshow(Ir);
+
+T_Gaus= imgaussfilt(I1);
+I = imbinarize(T_Gaus,0.5);
+%figure,imshow(I), title('Binarized image after applying Gaussian filter');
+
+BW = edge(I,'Sobel');    % edge detection 
+%figure, imshow(BW), title('Edge detection using sobel method');
+se1 = strel('disk', 5);
+se2 = strel('disk', 4);
+BW = imdilate(BW,se1);
+%figure, imshow(BW) , title('After applying dilation filter');
+BW = imerode(BW,se2);
+%figure, imshow(BW), title('After applying erosion filter');
+F = imfill(BW,'holes');   % fill the image 
+%figure, imshow(F), title('Filling the gaps in the edges');
+
+%Fill a gap 
+se = strel('disk',2);
+F = imclose(F,se);
+F = bwareaopen(F, 5000);
+figure, imshow(F), title('Removing small boundaries');
+
+%% boundries 
+B1 = bwboundaries(F, 8, 'noholes'); % Boundary detection
+B_size = size(B1);
+
+
+% figure, imshow(Ir);
+% hold on
+% for k = 1:B_size(1,1)
+% plot(B1{k}(:,2),B1{k}(:,1),'b*')
+% end
+% hold off
+
+%%
+
+D = zeros(N,N,B_size(1,1));   % ar tag data
+ID  = zeros(B_size(1,1));         % tag id 
+
+for k = 1:B_size(1,1)
+    BB = B1{k};
+ ps = dpsimplify(BB,10); %Douglas-Peucker Algorithm
+
+        ps_size = size(ps);
+        ps_size
+        final_corners = zeros(4,2);
+
+             if(  ps_size(1) == 5)                                              %If it is a polygon with 4 corners, then detect as quad
+            for k1=1:ps_size(1)-1
+                for kk=1:2 
+                    final_corners(k1,kk) = ps(k1,kk); %Only four corners
+                end
+            end
+
+            fc = final_corners; 
+            Ir = insertShape(Ir,'FilledPolygon',[fc(1,2) fc(1,1) fc(2,2) fc(2,1) fc(3,2) fc(3,1) fc(4,2) fc(4,1)],'Color','green','Opacity',0.6);
+
+%% HOMOGRAPHY
+             %Reference marker points
+
+             quad_pts(1,:) = [1, 1];
+             quad_pts(2,:) = [600, 1];
+             quad_pts(3,:) = [600, 600];
+             quad_pts(4,:) = [1, 600];
+
+             %Corner points
+             final_pts = [final_corners(:,2), final_corners(:,1)];
+
+
+             %Estimate homography with the 2 sets of four points
+             H = fitgeotrans( final_pts,quad_pts,'projective');
+
+             %Warp the marker to image plane
+             RA = imref2d([quad_pts(3,1) quad_pts(3,2)], [1 quad_pts(3,1)-1], [1 quad_pts(3,1)-1]);
+             [warp,r] = imwarp(I1, H, 'OutputView', RA); 
+
+
+              th = graythresh(warp);
+              markBin = imbinarize(warp, th);
+              se3 = strel('square', 1);
+              markBin = imerode(markBin,se3);
+
+              % decoding the ar tag 
+              [D(:,:,k),ID(k),img] = DECODE_AR_TAG(markBin,N);
+
+              figure,imshow(img), title('AR tag after decoding');              
+             end
+end
+figure,imshow(Ir);

algorithms/AR tag detection/scripts/DECODE_AR_TAG.m

@@ -0,0 +1,39 @@
+function [D,ID,img] = DECODE_AR_TAG(img,N)
+
+D = zeros(N,N);
+ID = 0;
+
+[av,ah] = size(img);
+x = 1;
+
+for i = 1:N
+
+   y = 1;
+    for j = 1:N 
+        D(j,i) = 1;
+        p1 = int16(ah/(2*N));
+        p2 = int16(av/(2*N));
+        D(j,i) = img(y+p2,x+p1 );
+
+        ID = ID + i*j*D(j,i);
+
+        if D(j,i)
+
+         img(y+int16(av/(2*N))-10:y+int16(av/(2*N))+10 , x+int16(ah/(2*N))-10:x+int16(ah/(2*N))+10) = zeros(21,21);
+
+        else 
+
+        img(y+int16(av/(2*N))-10:y+int16(av/(2*N))+10 , x+int16(ah/(2*N))-10:x+int16(ah/(2*N))+10) = ones(21,21);
+
+        end
+
+         y = y + int16(av/N);
+    end
+
+    x = x + int16(ah/N);
+
+end
+
+
+end
+

algorithms/AR tag detection/scripts/dpsimplify.m

@@ -0,0 +1,243 @@
+function [ps,ix] = dpsimplify(p,tol)
+
+% Recursive Douglas-Peucker Polyline Simplification, Simplify
+%
+% [ps,ix] = dpsimplify(p,tol)
+%
+% dpsimplify uses the recursive Douglas-Peucker line simplification 
+% algorithm to reduce the number of vertices in a piecewise linear curve 
+% according to a specified tolerance. The algorithm is also know as
+% Iterative Endpoint Fit. It works also for polylines and polygons
+% in higher dimensions.
+%
+% In case of nans (missing vertex coordinates) dpsimplify assumes that 
+% nans separate polylines. As such, dpsimplify treats each line
+% separately.
+%
+% For additional information on the algorithm follow this link
+% http://en.wikipedia.org/wiki/Ramer-Douglas-Peucker_algorithm
+%
+% Input arguments
+%
+%     p     polyline n*d matrix with n vertices in d 
+%           dimensions.
+%     tol   tolerance (maximal euclidean distance allowed 
+%           between the new line and a vertex)
+%
+% Output arguments
+%
+%     ps    simplified line
+%     ix    linear index of the vertices retained in p (ps = p(ix))
+%
+% Examples
+%
+% 1. Simplify line 
+%
+%     tol    = 1;
+%     x      = 1:0.1:8*pi;
+%     y      = sin(x) + randn(size(x))*0.1;
+%     p      = [x' y'];
+%     ps     = dpsimplify(p,tol);
+%
+%     plot(p(:,1),p(:,2),'k')
+%     hold on
+%     plot(ps(:,1),ps(:,2),'r','LineWidth',2);
+%     legend('original polyline','simplified')
+%
+% 2. Reduce polyline so that only knickpoints remain by 
+%    choosing a very low tolerance
+%
+%     p = [(1:10)' [1 2 3 2 4 6 7 8 5 2]'];
+%     p2 = dpsimplify(p,eps);
+%     plot(p(:,1),p(:,2),'k+--')
+%     hold on
+%     plot(p2(:,1),p2(:,2),'ro','MarkerSize',10);
+%     legend('original line','knickpoints')
+%
+% 3. Simplify a 3d-curve
+% 
+%     x = sin(1:0.01:20)'; 
+%     y = cos(1:0.01:20)'; 
+%     z = x.*y.*(1:0.01:20)';
+%     ps = dpsimplify([x y z],0.1);
+%     plot3(x,y,z);
+%     hold on
+%     plot3(ps(:,1),ps(:,2),ps(:,3),'k*-');
+%
+%
+%
+% Author: Wolfgang Schwanghart, 13. July, 2010.
+% w.schwanghart[at]unibas.ch
+
+
+if nargin == 0
+    help dpsimplify
+    return
+end
+
+error(nargoutchk(2, 2, nargin))
+
+% error checking
+if ~isscalar(tol) || tol<0;
+    error('tol must be a positive scalar')
+end
+
+
+% nr of dimensions
+nrvertices    = size(p,1); 
+dims    = size(p,2);
+
+% anonymous function for starting point and end point comparision
+% using a relative tolerance test
+compare = @(a,b) abs(a-b)/max(abs(a),abs(b)) <= eps;
+
+% what happens, when there are NaNs?
+% NaNs divide polylines.
+Inan      = any(isnan(p),2);
+% any NaN at all?
+Inanp     = any(Inan);
+
+% if there is only one vertex
+if nrvertices == 1 || isempty(p);
+    ps = p;
+    ix = 1;
+
+% if there are two 
+elseif nrvertices == 2 && ~Inanp;
+    % when the line has no vertices (except end and start point of the
+    % line) check if the distance between both is less than the tolerance.
+    % If so, return the center.
+    if dims == 2;
+        d    = hypot(p(1,1)-p(2,1),p(1,2)-p(2,2));
+    else
+        d    = sqrt(sum((p(1,:)-p(2,:)).^2));
+    end
+
+    if d <= tol;
+        ps = sum(p,1)/2;
+        ix = 1;
+    else
+        ps = p;
+        ix = [1;2];
+    end
+
+elseif Inanp;
+
+    % case: there are nans in the p array
+    % --> find start and end indices of contiguous non-nan data
+    Inan = ~Inan;
+    sIX = strfind(Inan',[0 1])' + 1; 
+    eIX = strfind(Inan',[1 0])'; 
+
+    if Inan(end)==true;
+        eIX = [eIX;nrvertices];
+    end
+
+    if Inan(1);
+        sIX = [1;sIX];
+    end
+
+    % calculate length of non-nan components
+    lIX = eIX-sIX+1;   
+    % put each component into a single cell
+    c   = mat2cell(p(Inan,:),lIX,dims);
+
+    % now call dpsimplify again inside cellfun. 
+    if nargout == 2;
+        [ps,ix]   = cellfun(@(x) dpsimplify(x,tol),c,'uniformoutput',false);
+        ix        = cellfun(@(x,six) x+six-1,ix,num2cell(sIX),'uniformoutput',false);
+    else
+        ps   = cellfun(@(x) dpsimplify(x,tol),c,'uniformoutput',false);
+    end
+
+    % write the data from a cell array back to a matrix
+    ps = cellfun(@(x) [x;nan(1,dims)],ps,'uniformoutput',false);    
+    ps = cell2mat(ps);
+    ps(end,:) = [];
+
+    % ix wanted? write ix to a matrix, too.
+    if nargout == 2;
+        ix = cell2mat(ix);
+    end
+
+
+else
+
+
+% if there are no nans than start the recursive algorithm
+ixe     = size(p,1);
+ixs     = 1;
+
+% logical vector for the vertices to be retained
+I   = true(ixe,1);
+
+% call recursive function
+p   = simplifyrec(p,tol,ixs,ixe);
+ps  = p(I,:);
+
+% if desired return the index of retained vertices
+if nargout == 2;
+    ix  = find(I);
+end
+
+end
+
+% _________________________________________________________
+function p  = simplifyrec(p,tol,ixs,ixe)
+
+    % check if startpoint and endpoint are the same 
+    % better comparison needed which included a tolerance eps
+
+    c1 = num2cell(p(ixs,:));
+    c2 = num2cell(p(ixe,:));   
+
+    % same start and endpoint with tolerance
+    sameSE = all(cell2mat(cellfun(compare,c1(:),c2(:),'UniformOutput',false)));
+
+
+    if sameSE; 
+        % calculate the shortest distance of all vertices between ixs and
+        % ixe to ixs only
+        if dims == 2;
+            d    = hypot(p(ixs,1)-p(ixs+1:ixe-1,1),p(ixs,2)-p(ixs+1:ixe-1,2));
+        else
+            d    = sqrt(sum(bsxfun(@minus,p(ixs,:),p(ixs+1:ixe-1,:)).^2,2));
+        end
+    else    
+        % calculate shortest distance of all points to the line from ixs to ixe
+        % subtract starting point from other locations
+        pt = bsxfun(@minus,p(ixs+1:ixe,:),p(ixs,:));
+
+        % end point
+        a = pt(end,:)';
+
+        beta = (a' * pt')./(a'*a);
+        b    = pt-bsxfun(@times,beta,a)';
+        if dims == 2;
+            % if line in 2D use the numerical more robust hypot function
+            d    = hypot(b(:,1),b(:,2));
+        else
+            d    = sqrt(sum(b.^2,2));
+        end
+    end
+
+    % identify maximum distance and get the linear index of its location
+    [dmax,ixc] = max(d);
+    ixc  = ixs + ixc; 
+
+    % if the maximum distance is smaller than the tolerance remove vertices
+    % between ixs and ixe
+    if dmax <= tol;
+        if ixs ~= ixe-1;
+            I(ixs+1:ixe-1) = false;
+        end
+    % if not, call simplifyrec for the segments between ixs and ixc (ixc
+    % and ixe)
+    else   
+        p   = simplifyrec(p,tol,ixs,ixc);
+        p   = simplifyrec(p,tol,ixc,ixe);
+
+    end
+
+end
+end