Sudoku Solver or How to solve Sudoku from its pic | Part 3

Hello.

In this post, I will discuss about the first few phases of the work. The system needs lots of improvements, which include technology used and the approach followed. Anyways, it was all part of a learning process.

First the image is converted to greyscale, then thresholding is performed on the image and we obtain a binary image. Thresholding can be used to distinguish background from foreground. Usually in simple images, the threshold value is the mean or median of all the pixel values present in it. Pixels whose value is greater than the threshold is assigned a value 1 or otherwise 0. The situation can also be reversed. However, setting a global thresholding value has a drawback - it cannot handle changing lighting conditions e.g. those happening due to strong illumination gradient or shadows. This problem can be solved using adaptive thresholding (also known as local or dynamic thresholding). It is based on the assumption that the illumination throughout a small region will be uniform. This method selects individual thresholds for each pixel based on its neighborhood.

After we get the binary image, we search for the largest blob (based on the assumption that the grid will be the largest blob). Using hough transform on the grid, we can get an estimate of the amount of rotation. We rotate the image to get an upright puzzle area. We crop out the puzzle area and divide it into 9x9 parts. Using template matching, we identify the numbers. template matching is a very naive method. It got confused between 3 and 8, 1 and 4, 2 and 9, 5 and 6. We solved this problem by finding the number of holes in the figure. For example, 8 and 3 has two and zero holes respectively.

The code is as follows. It is rough at the edges and a lot can be improved.

warning('off') 

%% adaptive threshholding
I=imread('sudoku_1.jpg');
I=imrotate(I,30,'bilinear');
J=mat2gray(I);
K=imfilter(J,fspecial('average',10),'replicate');
C = 0.05;
sIM=K-J-C;
bw=im2bw(sIM,0);
%%imshow(bw)

%% get biggest blob
[label, num] = bwlabel(bw,4);
stats = regionprops (label, 'basic');
mx=max([stats.Area]);
mp = find([stats.Area]==mx);
X=label==mp;
%%imshow(X) 

 

%% crop and get the puzzle area
xy=max(1,floor(stats(mp).BoundingBox));
fI=bw(xy(2):xy(2)+xy(4),xy(1):xy(1)+xy(3));
%%imshow(fI) 

%% Rotation
cleanI=X;
cleanI=edge(cleanI,'canny');
[H,theta,rho] = hough(cleanI);
P = houghpeaks(H,5,'threshold',ceil(0.3*max(H(:))));
lines = houghlines(cleanI,theta,rho,P,'FillGap',5,'MinLength',7);
angl=max([lines.theta]);
fI=imrotate(fI,angl-90,'bilinear');
[label, num] = bwlabel(fI,4);
stats = regionprops (label, 'basic');
mx=max([stats.Area]);
mp = find([stats.Area]==mx);
X=label==mp;
xy=max(1,floor(stats(mp).BoundingBox));
fI=fI(xy(2):xy(2)+xy(4),xy(1):xy(1)+xy(3));
fI=fI-X(xy(2):xy(2)+xy(4),xy(1):xy(1)+xy(3)); 

%% segment and identify
sz=size(fI);
s=round(sz/9);
XY=zeros(9,26,18);
box=zeros(9,9);
perc=round(s(1).*s(2)*0.04);
for i=1:9
    filename=strcat('',num2str(i),'.jpg');
    XY(i,:,:)=im2bw(imread(filename));
end
for i =0:8
    for j=0:8
        %subplot(9,9,(9*i+j+1))
        xm=min(sz(1),(i+1)*s(1));
        ym=min(sz(2),(j+1)*s(2));
        gh=fI(i*s(1)+1:xm,j*s(2)+1:ym);
        gh=bwareaopen(gh,perc);
        [label n]=bwlabel(gh,8);
     
        if n >0
            stats=regionprops(label,'BoundingBox','PixelIdxList');
         
            for region =1:length(stats)
                xy=stats(region).BoundingBox;
                if (xy(3) <10 || xy(4) < 10 )
                    label(stats(region).PixelIdxList) = 0;
                end
            end
            gh = (label~=0);
            [label n]=bwlabel(gh,8);
         
            if n < 1
                %imshow(gh)
                continue
            end
            stats=regionprops(label,'basic');
            mx=max([stats.Area]);
            mp = find([stats.Area]==mx);
            gh=label;
            gh(gh==mp)=1;
            xy=stats(mp).BoundingBox;
            if (xy(3)/xy(4) > 0.75 || xy(3)/xy(4) <0.25)
                %imshow(gh)
                continue
            end
            szx=size(gh);
            gh=gh(xy(2):min(xy(2)+xy(4),szx(1)), xy(1):min(xy(1)+xy(3),szx(2)));
            gh=im2bw(imresize(gh,[26,18]));
            gh=bwareaopen(gh,120);
            mx_corr=0;
            for k=1:9
                temp_(:,:)=XY(k,:,:) ;
                temp=sum(sum(gh.*temp_));
                if mx_corr < temp
                    mx_corr=temp;
                    box(i+1,j+1)=k;
                end
            end
            se1=strel('line',2,0);
            se2=strel('line',2,90);
            gh=imdilate(gh,[se1 se2]);
            gh=imerode(gh,[se1 se2]);
            if (box(i+1,j+1)==8 || box(i+1,j+1)==3)
               nx=regionprops(gh,'EulerNumber');
               if (nx(1).EulerNumber==1)
                   box(i+1,j+1)=3;
               else
                   box(i+1,j+1)=8;
               end
            end
            if (box(i+1,j+1)==4 || box(i+1,j+1)==1)
               nx=regionprops(gh,'EulerNumber');
               if (nx(1).EulerNumber==0)
                   box(i+1,j+1)=4;
               else
                   box(i+1,j+1)=1;
               end
            end
         
            if (box(i+1,j+1)==5 || box(i+1,j+1)==6)
               nx=regionprops(gh,'EulerNumber');
               if (nx(1).EulerNumber<=0)
                   box(i+1,j+1)=6;
               else
                   box(i+1,j+1)=5;
               end
            end
            if (box(i+1,j+1)==9 || box(i+1,j+1)==2)
                nx=regionprops(gh,'EulerNumber');
                if (nx(1).EulerNumber<=0)
                    box(i+1,j+1)=9;
                else
                    box(i+1,j+1)=2;
                end
            end
        end
        %imshow(gh)
    end
end
file_=fopen('out_dip.txt','w');
fprintf(file_,'%d',box');
fclose(file_);
[status result]=system('python do.py');
clc
disp(result)
warning on