111

function sudokuvideo_fn()

%%%%%%%%%% Parameters you need to set: %%%%%%%%%%%%

% The initial threshold. Threshold at "bwthresh" times darker than the
% median in each block. Use the slider bar to adjust this in real time.
blksize = 92;  %blocksize for block processing /

% Specify whether to overlay the solution not. 1 or 0. /
tracking = 1;

% You may need to adjust your video settings here:
% Also, if you use RGB instead of YUY2, you *may* need to adjust
% lines 14-16 of sudokuvideo_fn_trace.m. I'm not sure, my camera doesn't output RGB.
imaqreset
obj = videoinput('winvideo',1,'YUY2_640x480');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
try
    
    %Initialize various parameters, and load in the template data
    load TEMPLATEDATA
    set(obj,'framesperTrigger',10,'TriggerRepeat',Inf);
    start(obj);
    
    A_tmin = 30; % Bounds for the digit pixel area
    A_tmax = 1000;
    digitbox_minarea = 20; % Bounds for the digit bounding box area
    digitbox_maxarea = 25^2;
    
    found = 0; % If a solution is FOUND, then start tracking.
    
    % h = The B/W image. h_g = Green box. h_pts = magenta box
    h = imshow(zeros(480,640));
    hold on;
    h_g = plot(90+[0 460 460 0 0],10+[0 0 460 460 0],'g');
    h_pts  = [];
    
    Plocal_old = []; % Variable to check the positions and values
    
    same = 0;
    samemax = 3; %Consecutive Matches
    
    figure(1);
    hslider = uicontrol('style','slider','units','n','pos',[0 0 1 0.05],'min',1.1,'max',1.4,'val',1.3);
    uicontrol('style','text','units','n','pos',[0 0.05 0.2 0.03],'String','Threshold adjustment:');
    
    while islogging(obj);
        
        bwthresh = get(hslider,'value');
        
        Icam = getdata(obj,1);
        %Icam = imread('sample.bmp'); %<--- For debugging
        
        Icam = Icam(:,:,1);
        
        % I0 represents only the inner green square
        I0 = Icam(10+(1:460),90+(1:460));
        flushdata(obj);
        
        %Block processed threshhold
        makebw2 = @(I) im2bw(I.data,median(double(I.data(:)))/bwthresh/255);
        IBW = ~blockproc(I0,[blksize blksize],makebw2);
        
        % Noise reduction and border elimination
        IBW = imclose(IBW,[1 1; 1 1]);
        I = IBW;
        I = bwareaopen(I,A_tmin);
        I = imclearborder(I);
        
        % Iout is an augmented 640x480 version of I, for display.
        Iout = zeros(480,640);
        Iout(10+(1:460),90+(1:460)) = I;
        
        % Show Iout
        if ~tracking || ~found
            set(h,'Cdata',Iout);
        end
        
        R = regionprops(I,'Area','BoundingBox','PixelList');
        NR = numel(R);
        
        % Find the largest object -> assume it to be the puzzle outline
        maxArea = 0;
        for k = 1:NR
            A(k) = prod(R(k).BoundingBox(3:4));
            if R(k).Area > maxArea
                maxArea = R(k).Area;
                kmax = k;
            end
        end
        figure(1);
        itsok = 0;
        if maxArea > 1000 && A(kmax) > 150^2 % Here we find the corners of the box
            itsok = 1;
            set(h_g,'color','g');
            BBmax = R(kmax).BoundingBox;
            DIAG1 = sum(R(kmax).PixelList,2);
            DIAG2 = diff(R(kmax).PixelList,[],2);
            
            [m,dUL] = min(DIAG1);
            [m,dDR] = max(DIAG1);
            [m,dDL] = min(DIAG2);
            [m,dUR] = max(DIAG2);
            pts = R(kmax).PixelList([dUL dDL dDR dUR dUL],:);
            h_pts = plot(90+pts(:,1),10+pts(:,2),'m');
            
            XYLIMS = [BBmax(1) + [0 BBmax(3)] BBmax(2) + [0 BBmax(4)]];
        end
        
        % If we are in tracking mode, have found a solution, and still have a
        % bounding quadrilateral to work with, then call the overlay function
        if tracking && found && numel(pts) == 10 && polyarea(pts(1:4,1),pts(1:4,2)) > 150^2
            sudokuvideo_fn_trace(h,Icam,Iover,pts);
            try, delete(h_y); end; h_y = [];
            try, delete(h_pts); end; h_pts = [];
            continue
        end
        set(h,'Cdata',Iout);
        drawnow;
        
        % If a decent puzzle outline could not be found, then continue
        if ~itsok
            set(h_g,'color','r');
            continue
        end
        
        
        try, delete(h_y); end; h_y = [];
        try, delete(h_pts); end; h_pts = [];
        if found, continue; end
        
        h_digitcircles = [];
        kgood = zeros(1,NR);
        Pnew = zeros(NR,2);
        % Now we need to identify those blobs which are nearly digit sized,
        % and are inside our puzzle quadrilateral.
        for k = 1:NR
            if R(k).Area < A_tmax && A(k) > digitbox_minarea && A(k) < digitbox_maxarea ...
                    && R(k).BoundingBox(3) < 40 && R(k).BoundingBox(4) < 40 ...
                    && R(k).BoundingBox(3) > 2 && R(k).BoundingBox(4) > 2
                
                Pnew(k,:) = [R(k).BoundingBox(1)+R(k).BoundingBox(3)/2 R(k).BoundingBox(2)+R(k).BoundingBox(4)/2];
                
                
                if inpolygon(Pnew(k,1),Pnew(k,2),pts(:,1),pts(:,2))
                    h_digitcircles(k) = plot(90+Pnew(k,1),10+Pnew(k,2),'ro','markersize',24);
                    kgood(k) = 1;
                else
                    h_digitcircles(k) = plot(90+Pnew(k,1),10+Pnew(k,2),'bo','markersize',24);
                end
                
            end
        end
        [kgoodvals,kgoodlocs] = find(kgood);
        Pnew = Pnew(kgoodlocs,:);
        
        % This is where we send the remaining candidate blobs for
        % identification:
        % Plocal represents the puzzle as an Nx3 matrix: [x, y, digit]
        [Plocal,del] = identifynumbers_fun(pts,Pnew,NT,IBW);
        
        % Remove circles from non-digit locations
        try
            delete(h_digitcircles(kgoodlocs(del)));
            h_digitcircles(kgoodlocs(del)) = 0;
        catch
            keyboard
        end
        drawnow;
        try, delete(nonzeros(h_digitcircles)); end;
        
        % If good numbers were not identified, continue
        if isnan(Plocal)
            continue
        end
        
        % We need to wait until the identified digits are the same for a
        % specified number of times consecutively. If so, then call the
        % solver
        if isequal(Plocal_old, Plocal) && ~isempty(Plocal)
            same = same+1;
            if same == samemax-1
                title([num2str(same+1) ' / ' num2str(samemax)],'FontSize',12);
                M = zeros(9);
                for k = 1:size(Plocal,1)
                    M(Plocal(k,1),Plocal(k,2)) = Plocal(k,3);
                end
                M = M'
                M_sol = drawgraph(M);
                figure;
                imshow(I0);
                % If a solution was found, M_sol is non-empty.
                if ~isempty(M_sol)
                    if ~tracking, return; end
                    found = 1;
                    Iover = ~solution2image(M,M_sol); %Make the solution overlay
                    figure(1);
                else
                    same = 0;
                end
            end
        else
            same = 0;
        end
        
        title([num2str(same+1) ' / ' num2str(samemax)],'FontSize',12);
        Plocal_old = Plocal;
        
    end
catch
    % This attempts to take care of things when the figure is closed
    stop(obj);
    imaqreset
    %keyboard
    figure(1);
    imshow(get(h,'cdata'));
    drawnow;
end



function Msol = drawsudoku(M)

% Basically, show the input M, and solution Msol, in a table

D1 = cell(9);
for m = 1:9
    for n = 1:9
        if M(m,n)
            D1{m,n} = M(m,n);
        end
    end
end

%Call the sudoku solver, return if it fails
Msol = sudoku_solver(M);
if isempty(Msol)
    return
else
    D2 = num2cell(Msol);
end

figure;
set(gcf,'units','pixels','Position',[200 200 800 400]);


h1 = uitable('Data',D1,'FontSize',16,'ColumnWidth',...
    num2cell(repmat(30,1,9)),'columne',false(1,9));
set(h1,'units','norm','position',[.05 .05 .4 .9]);
h2 = uitable('Data',D2,'FontSize',16,'ColumnWidth',num2cell(repmat(30,1,9)));
set(h2,'units','norm','position',[.55 .05 .4 .9]);



function [Plocal,del] = identifynumbers_fun(pts,Pnew,NT,IBW)

% This function uses the vertices of the puzzle quadrilateral "pts"
% to identify which blobs are in digit locations, and then compares those
% blobs to template digits for identification

% No blobs, quick return
if isempty(Pnew)
    Plocal = nan;
    del = [];
    return
end

% Use the vertices to transform the blob coordinates
try
    T = cp2tform(pts(1:4,:),[0.5 0.5; 9.5 0.5; 9.5 9.5; 0.5 9.5],'projective');
catch
    Plocal = nan;
    del = [];
    return
end
Plocal = (tformfwd(T,Pnew));
Plocal = round(2*Plocal)/2;

% "del" represents those blobs which are not in digit-like locations
del = find(sum(Plocal - floor(Plocal) > 0 |  Plocal < 1 | Plocal > 9,2)) ;
Plocal(del,:) = [];
Pnew(del,:) = [];

if any(isnan(Plocal(:))) || isempty(Plocal)
    Plocal = nan;
    del = [];
    return
end


% The actual identification algorithm
try
    Plocal(end,3) = 0;
    
    for k = 1:size(Pnew,1)
        for s = [0 -1 1 -2 2 -3 3 -4 4 -5 5]
            N = bwselect(IBW,Pnew(k,1) + s ,Pnew(k,2));

            if any(N(:))
                break
            end
        end
        if s == 5
            Plocal = nan;
            return
            %continue
        end
        
        [i,j] = find(N);
        N = N(min(i):max(i),min(j):max(j));
        N0 = N;
        
        % Resize to be 20x20
        N = imresize(N,[20 20]);
        
        %for each digit, S(v) represents the degree of matching
        for v = 1:9
            S(v) = sum(sum(N.*NT{v}));
        end
        
        Plocal(k,3) = find(S == max(S),1);
        
        if (Plocal(k,3) == 5 || Plocal(k,3) == 6) && abs(S(5) - S(6)) < 0.1 %If it's a 5 or 6, use the Euler number
            E = regionprops(N,'EulerNumber');
            if ~E(1).EulerNumber
                Plocal(k,3) = 6;
            end
        end
        
    end
catch
    %keyboard
end
Plocal = sortrows(Plocal);



function mysudokusolver(filename)
%% Read in a file
M = dlmread(filename);
  
%% Show the initial state
figure;
set(gcf,'units','pixels','Position',[200 200 800 400]);

D1 = cell(9);
for m = 1:9
    for n = 1:9
        if M(m,n)
            D1{m,n} = M(m,n);
        end
    end
end

%% Call the solver
    Msol = sudoku_solver(M);

%% Show the solution
if isempty(Msol)
    D2 = cell(9);
else
    D2 = num2cell(Msol);
end


h1 = uitable('Data',D1,'FontSize',16,'ColumnWidth',...
    num2cell(repmat(30,1,9)),'columne',false(1,9));
set(h1,'units','norm','position',[.05 .05 .4 .9]);
h2 = uitable('Data',D2,'FontSize',16,'ColumnWidth',num2cell(repmat(30,1,9)));
set(h2,'units','norm','position',[.55 .05 .4 .9]);



%% Read In a File / 
I_cam = imread('sample.bmp');%赋值
imshow(I_cam);%显示sanple
load TEMPLATEDATA

%% Crop the Image (1) / 
hold on
h_g = plot(90+[0 460 460 0 0],10+[0 0 460 460 0],'g');%画绿格子
hold off

%% Crop the Image (2) / 夋憸僼傽僀儖撉傒崬傒

I = I_cam(10+(1:460),90+(1:460));
imshow(I);%覆盖

%% Convert to Black and White / 黑白
makebw = @(I) im2bw(I.data,median(double(I.data(:)))/1.2/255);
I = ~blockproc(I,[92 92],makebw);

imshow(I);

%% Remove Noise / 除燥
I = bwareaopen(I,30);
imshow(I);
%% Clear the border / 去边
I = imclearborder(I);
imshow(I);

%% Find the largest box / 堦斣戝偒偄榞傪扵偡
hold on;
R = regionprops(I,'Area','BoundingBox','PixelList');
NR = numel(R);

maxArea = 0;
for k = 1:NR
    A(k) = prod(R(k).BoundingBox(3:4));
    if R(k).Area > maxArea
        maxArea = R(k).Area;
        kmax = k;
    end
end


BBmax = R(kmax).BoundingBox;
DIAG1 = sum(R(kmax).PixelList,2);
DIAG2 = diff(R(kmax).PixelList,[],2);

[m,dUL] = min(DIAG1);    [m,dDR] = max(DIAG1);
[m,dDL] = min(DIAG2);    [m,dUR] = max(DIAG2);

pts = R(kmax).PixelList([dUL dDL dDR dUR dUL],:);
h_pts = plot(pts(:,1),pts(:,2),'m','linewidth',3);

XYLIMS = [BBmax(1) + [0 BBmax(3)] BBmax(2) + [0 BBmax(4)]];

%% Identify objects inside the box /丂榞撪偺僆僽僕僃僋僩傪扵偡
A_tmin = 30; % Bounds for the digit pixel area
A_tmax = 1000;
digitbox_minarea = 20; % Bounds for the digit bounding box area
digitbox_maxarea = 25^2;

kgood = zeros(1,NR);
Pnew = zeros(NR,2);
        for k = 1:NR
            if R(k).Area < A_tmax && A(k) > digitbox_minarea && A(k) < digitbox_maxarea ...
                    && R(k).BoundingBox(3) < 40 && R(k).BoundingBox(4) < 40 ...
                    && R(k).BoundingBox(3) > 1 && R(k).BoundingBox(4) > 1
                
                Pnew(k,:) = [R(k).BoundingBox(1)+R(k).BoundingBox(3)/2 R(k).BoundingBox(2)+R(k).BoundingBox(4)/2];
                
                
                if inpolygon(Pnew(k,1),Pnew(k,2),pts(:,1),pts(:,2))
                    h_digitcircles(k) = plot(Pnew(k,1),Pnew(k,2),'ro','markersize',24);
                end
                
            end
        end


%% Draw the grid based on the corners / 妏偐傜僌儕僢僪傪嶌惉

T = cp2tform(pts(1:4,:),0.5 + [0 0; 9 0; 9 9; 0 9],'projective');
for n = 0.5 + 0:9, [x,y] = tforminv(T,[n n],[0.5 9.5]); plot(x,y,'g'); end
for n = 0.5 + 0:9, [x,y] = tforminv(T,[0.5 9.5],[n n]); plot(x,y,'g'); end
%% Only keep elements in the boxes / 傑偡偵擖偭偰偄傞梫慺偺傒僉乕僾
T = cp2tform(pts(1:4,:),[0.5 0.5; 9.5 0.5; 9.5 9.5; 0.5 9.5],'projective');
Plocal = (tformfwd(T,Pnew));
Plocal = round(2*Plocal)/2;

del = find(sum(Plocal - floor(Plocal) > 0 |  Plocal < 1 | Plocal > 9,2)) ;
Pnew(del,:) = [];

delete(nonzeros(h_digitcircles(del)));

%% Show the coordinate transforms / 嵗昗曄姺偺寢壥
figure;
T = cp2tform(pts(1:4,:),500*[0 0; 1 0; 1 1; 0 1],'projective');
IT = imtransform(double(I),T);
imshow(IT);
%% Show the template data / 僥儞僾儗乕僩僨乕僞傪昞帵
figure;

for n = 1:9
    subplot(3,3,n),imagesc(NT{n});
end
colormap gray;
%% Calculate the Solution / 僜儕儏乕僔儑儞傪寁嶼偡傞


Plocal = identifynumbers_fun(pts,Pnew,NT,I);
M = zeros(9);
for k = 1:size(Plocal,1)
    M(Plocal(k,2),Plocal(k,1)) = Plocal(k,3);
end
M_sol = drawgraph(M);

%% Generate an image from the solution / 僜儕儏乕僔儑儞偐傜夋憸傪嶌惉
I = solution2image(M,M_sol);
figure; imshow(I);
%% Overlay the solution on the original image / 夋憸傪僆乕僶儗僀

figure(1); 
clf;

T = cp2tform([1 1; 200 1; 200 200; 1 200],[pts(1:4,1),pts(1:4,2)],'projective');
I = imtransform(~I,T,'XData',[1 460], 'YData',[1 460],'XYscale',1);

Imask = zeros(480,640);
Imask(10 + (1:460), 90+(1:460)) = I;
Imask = I_cam .* uint8(~Imask);

h = imshow(cat(3,I_cam,Imask,I_cam));
set(h,'Cdatamapping','direct');

hold on;
plot(90 + pts(:,1), 10 + pts(:,2),'m')



function D = sudoku_solver(D)

onefound = 1;
N = 0;
poss = 1:9;

while onefound
    splits = {};
    onefound = 0;
    for m = 1:9
        rowdata = nonzeros(D(m,:)');
        for n = 1:9
            E = D(m,n);
            if E ~= 0, continue, end
            
            coldata = nonzeros(D(:,n));
            blk = [ceil(m/3) ceil(n/3)]-1;
            blkdata = nonzeros(D(blk(1)*3+[1:3],blk(2)*3+[1:3]));
            
            EE = zeros(1,9);
            RCB = [rowdata; coldata; blkdata(:)];
            EE(RCB) = 1;
            Enew = find(~EE);
            
            if isempty(Enew)
                D = []; return;
            elseif length(Enew) == 1;
                onefound = 1;
                D(m,n) = Enew;
                rowdata = nonzeros(D(m,:)');
            else
                splits{end+1} = [m n Enew];
            end
        end
    end
end

if isempty(splits)
    return
end

splitlength = cellfun(@length,splits);
splits = splits{find(splitlength == min(splitlength),1)};
m = splits(1); n = splits(2);

for test = 3:length(splits)
    D(m,n) = splits(test);
    D0 = sudoku_solver(D);
    if ~isempty(D0)
        D = D0;
        return
    end
end
D = [];



I_cam = imread('sample.bmp');
imshow(I_cam);
load TEMPLATEDATA
hold on
h_g = plot(90+[0 460 460 0 0],10+[0 0 460 460 0],'g');
I = I_cam(10+(1:460),90+(1:460));
imshow(I);
makebw = @(I) im2bw(I.data,median(double(I.data(:)))/1.2/255);
I = ~blockproc(I,[92 92],makebw);

imshow(I);
I = bwareaopen(I,30);
imshow(I);
I = imclearborder(I);
imshow(I);
hold on;
R = regionprops(I,'Area','BoundingBox','PixelList');
NR = numel(R);

maxArea = 0;
for k = 1:NR
    A(k) = prod(R(k).BoundingBox(3:4));
    if R(k).Area > maxArea
        maxArea = R(k).Area;
        kmax = k;
    end
end


BBmax = R(kmax).BoundingBox;
DIAG1 = sum(R(kmax).PixelList,2);
DIAG2 = diff(R(kmax).PixelList,[],2);

[m,dUL] = min(DIAG1);    [m,dDR] = max(DIAG1);
[m,dDL] = min(DIAG2);    [m,dUR] = max(DIAG2);

pts = R(kmax).PixelList([dUL dDL dDR dUR dUL],:);
h_pts = plot(pts(:,1),pts(:,2),'m','linewidth',3);

XYLIMS = [BBmax(1) + [0 BBmax(3)] BBmax(2) + [0 BBmax(4)]];

%% Identify objects inside the box /丂榞撪偺僆僽僕僃僋僩傪扵偡
A_tmin = 30; % Bounds for the digit pixel area
A_tmax = 1000;
digitbox_minarea = 20; % Bounds for the digit bounding box area
digitbox_maxarea = 25^2;

kgood = zeros(1,NR);
Pnew = zeros(NR,2);
        for k = 1:NR
            if R(k).Area < A_tmax && A(k) > digitbox_minarea && A(k) < digitbox_maxarea ...
                    && R(k).BoundingBox(3) < 40 && R(k).BoundingBox(4) < 40 ...
                    && R(k).BoundingBox(3) > 1 && R(k).BoundingBox(4) > 1
                
                Pnew(k,:) = [R(k).BoundingBox(1)+R(k).BoundingBox(3)/2 R(k).BoundingBox(2)+R(k).BoundingBox(4)/2];
                
                
                if inpolygon(Pnew(k,1),Pnew(k,2),pts(:,1),pts(:,2))
                    h_digitcircles(k) = plot(Pnew(k,1),Pnew(k,2),'ro','markersize',24);
                end
                
            end
        end


%% Draw the grid based on the corners / 妏偐傜僌儕僢僪傪嶌惉

T = cp2tform(pts(1:4,:),0.5 + [0 0; 9 0; 9 9; 0 9],'projective');
for n = 0.5 + 0:9, [x,y] = tforminv(T,[n n],[0.5 9.5]); plot(x,y,'g'); end
for n = 0.5 + 0:9, [x,y] = tforminv(T,[0.5 9.5],[n n]); plot(x,y,'g'); end
%% Only keep elements in the boxes / 傑偡偵擖偭偰偄傞梫慺偺傒僉乕僾
T = cp2tform(pts(1:4,:),[0.5 0.5; 9.5 0.5; 9.5 9.5; 0.5 9.5],'projective');
Plocal = (tformfwd(T,Pnew));
Plocal = round(2*Plocal)/2;

del = find(sum(Plocal - floor(Plocal) > 0 |  Plocal < 1 | Plocal > 9,2)) ;
Pnew(del,:) = [];



I_cam = imread('sample.bmp');
imshow(I_cam);
load TEMPLATEDATA
hold on
h_g = plot(90+[0 460 460 0 0],10+[0 0 460 460 0],'g');
I = I_cam(10+(1:460),90+(1:460));
imshow(I);
makebw = @(I) im2bw(I.data,median(double(I.data(:)))/1.2/255);
I = ~blockproc(I,[92 92],makebw);

imshow(I);
I = bwareaopen(I,30);
imshow(I);
I = imclearborder(I);
imshow(I);
hold on;
R = regionprops(I,'Area','BoundingBox','PixelList');
NR = numel(R);

maxArea = 0;
for k = 1:NR
    A(k) = prod(R(k).BoundingBox(3:4));
    if R(k).Area > maxArea
        maxArea = R(k).Area;
        kmax = k;
    end
end


BBmax = R(kmax).BoundingBox;
DIAG1 = sum(R(kmax).PixelList,2);
DIAG2 = diff(R(kmax).PixelList,[],2);

[m,dUL] = min(DIAG1);    [m,dDR] = max(DIAG1);
[m,dDL] = min(DIAG2);    [m,dUR] = max(DIAG2);

pts = R(kmax).PixelList([dUL dDL dDR dUR dUL],:);
h_pts = plot(pts(:,1),pts(:,2),'m','linewidth',3);

XYLIMS = [BBmax(1) + [0 BBmax(3)] BBmax(2) + [0 BBmax(4)]];

%% Identify objects inside the box /丂榞撪偺僆僽僕僃僋僩傪扵偡
A_tmin = 30; % Bounds for the digit pixel area
A_tmax = 1000;
digitbox_minarea = 20; % Bounds for the digit bounding box area
digitbox_maxarea = 25^2;

kgood = zeros(1,NR);
Pnew = zeros(NR,2);
        for k = 1:NR
            if R(k).Area < A_tmax && A(k) > digitbox_minarea && A(k) < digitbox_maxarea ...
                    && R(k).BoundingBox(3) < 40 && R(k).BoundingBox(4) < 40 ...
                    && R(k).BoundingBox(3) > 1 && R(k).BoundingBox(4) > 1
                
                Pnew(k,:) = [R(k).BoundingBox(1)+R(k).BoundingBox(3)/2 R(k).BoundingBox(2)+R(k).BoundingBox(4)/2];
                
                
                if inpolygon(Pnew(k,1),Pnew(k,2),pts(:,1),pts(:,2))
                    h_digitcircles(k) = plot(Pnew(k,1),Pnew(k,2),'ro','markersize',24);
                end
                
            end
        end


%% Draw the grid based on the corners / 妏偐傜僌儕僢僪傪嶌惉

T = cp2tform(pts(1:4,:),0.5 + [0 0; 9 0; 9 9; 0 9],'projective');
for n = 0.5 + 0:9, [x,y] = tforminv(T,[n n],[0.5 9.5]); plot(x,y,'g'); end
for n = 0.5 + 0:9, [x,y] = tforminv(T,[0.5 9.5],[n n]); plot(x,y,'g'); end
%% Only keep elements in the boxes / 傑偡偵擖偭偰偄傞梫慺偺傒僉乕僾
T = cp2tform(pts(1:4,:),[0.5 0.5; 9.5 0.5; 9.5 9.5; 0.5 9.5],'projective');
Plocal = (tformfwd(T,Pnew));
Plocal = round(2*Plocal)/2;

del = find(sum(Plocal - floor(Plocal) > 0 |  Plocal < 1 | Plocal > 9,2)) ;
Pnew(del,:) = [];



clc;
clear all;
close all;
%% 读取图像
root='./data';
img=read_train(root);
%% 提取特征
img_feature=feature_lattice(img);
%% 构造标签
class=10;
numberpclass=500;
ann_label=zeros(class,numberpclass*class);
ann_data=img_feature;
for i=1:class
 for j=numberpclass*(i-1)+1:numberpclass*i
     ann_label(i,j)=1;
 end
end

%% 选定训练集和测试集
k=rand(1,numberpclass*class);  
[m,n]=sort(k);  
ntraindata=4500;
ntestdata=500;
train_data=ann_data(:,n(1:ntraindata));
test_data=ann_data(:,n(ntraindata+1:numberpclass*class));
train_label=ann_label(:,n(1:ntraindata));
test_label=ann_label(:,n(ntraindata+1:numberpclass*class));
%% BP神经网络创建，训练和测试
net=network_train(train_data,train_label);
predict_label=network_test(test_data,net);
%% 正确率计算
[u,v]=find(test_label==1);
label=u';
error=label-predict_label;
accuracy=size(find(error==0),2)/size(label,2)


function feature = feature_lattice(img)
% 输入:黑底白字的二值图像。输出：35维的网格特征
% ======提取特征，转成5*7的特征矢量,把图像中每10*10的点进行划分相加，进行相加成一个点=====%
%======即统计每个小区域中图像象素所占百分比作为特征数据====%
for i=1:length(img);
bw2=im2bw(img{i},graythresh(img{i}));
bw_7050=imresize(bw2,[70,50]);
for cnt=1:7
    for cnt2=1:5
        Atemp=sum(bw_7050(((cnt*10-9):(cnt*10)),((cnt2*10-9):(cnt2*10))));%10*10box
        lett((cnt-1)*5+cnt2)=sum(Atemp);
    end
end
lett=((100-lett)/100);
lett=lett';
feature(:,i)=lett;
end


function out = network_test(test_data,net)
% 输入：测试数据的特征和真值。输出：测试数据的label以及误差图
% BP网络预测

an=sim(net,test_data);
for i=1:length(test_data)
    out(i)=find(an(:,i)==max(an(:,i)));
end

end



function net = network_train(train_data,train_label )
% 输入：训练图像特征和label。输出：训练好的神经网络

% BP网络训练
% 初始化网络结构
layer=25;
net=newff(train_data,train_label,layer);
net.trainParam.epochs=1;
net.trainParam.lr=0.1;
net.trainParam.goal=0.001;
net.trainFcn='trainrp';
% 网络训练
net=train(net,train_data,train_label);
end




function [imglist] = read_train(root)
% ni为读取图片张数，n为文件夹数目
%========读取文件夹========%
out_Files = dir(root);%展开
tempind=0;
imglist=cell(0);
n=length(out_Files);
%========读取文件========%
for i = 1:n;
    if strcmp(out_Files(i).name,'.')|| strcmp(out_Files(i).name,'..')
    else
        rootpath=strcat(root,'/',out_Files(i).name);
        in_filelist=dir(rootpath);
        ni=length(in_filelist);
        for j=1:ni
            if strcmp(in_filelist(j).name,'.')|| strcmp(in_filelist(j).name,'..')|| strcmp(in_filelist(j).name,'Desktop_1.ini')|| strcmp(in_filelist(j).name,'Desktop_2.ini')
            else
                tempind=tempind+1;
                imglist{tempind}=imread(strcat(rootpath,'/',in_filelist(j).name));
            end
        end
    end
end
end