如何消除数独方格中的凸缺陷？

Question

我正在做一个有趣的项目：使用OpenCV（如Google goggles等）从输入图像解决数独问题。我已经完成了任务，但最后发现了一个小问题，所以来这里求助。

我使用OpenCV 2.3.1的Python API进行编程。

以下是我的做法:

1. 读取图像
2. 查找轮廓
3. 选择最大面积的轮廓，（也要略微等同于正方形）。

4. 查找角点。

   例如下面给出的:

   

   (请注意，绿线正确地与数独的真实边界重合，因此可以正确地扭曲数独。查看下一张图片)

5. 将图像扭曲成完美的正方形

   例如图像:

   

6. 执行OCR（我使用了Simple Digit Recognition OCR in OpenCV-Python中提供的方法）

这种方法效果很好。

问题：

检查this image.

对这张图片执行第4步骤会得到以下结果：

红色线条是数独边界的真实轮廓，绿色线条是近似轮廓，将成为变形图像的轮廓。当然，在数独的顶部边缘，绿线和红线之间存在差异。因此，在变形时，我无法获得数独的原始边界。

我的问题是：如何在数独的正确边界上进行变形，即红线，或者如何消除红线和绿线之间的差异？OpenCV中是否有任何方法可以实现这一点？

- Abid Rahman K

---

根据你所说的，你是基于角点进行检测的，而红线和绿线在这一点上达成了一致。我不太了解OpenCV，但我可以推测你可能想在这些角点之间检测出直线，并根据它们进行扭曲处理。 - Danica

也许可以强制连接角点的线与图像中的粗黑色像素重合。也就是说，不要让绿色线条只是在角点之间找到一条直线，而是强制它们穿过粗黑色像素。我认为这会使你的问题更加困难，而且我不知道有哪些OpenCV内置函数能立即对你有用。 - ely

@ Dougal：我认为绘制的绿线是红线的近似直线，因此它是连接这些角点之间的线。当我根据绿线进行变形时，在变形图像的顶部会出现弯曲的红线。（希望你能理解，我的解释可能有点不好） - Abid Rahman K

@ EMS：我认为所画的红线恰好在数独边界上。但问题是，如何将图像精确地扭曲到数独边界上（我的意思是，问题出在扭曲上，即将那些弯曲的边界转换为精确的正方形，就像我在第二张图片中展示的那样）。 - Abid Rahman K

Answer 1

我有一个可行的解决方案，但你需要自己将其翻译为OpenCV。它是用Mathematica编写的。

第一步是通过将每个像素除以闭运算的结果来调整图像的亮度：

src = ColorConvert[Import["http://davemark.com/images/sudoku.jpg"], "Grayscale"];
white = Closing[src, DiskMatrix[5]];
srcAdjusted = Image[ImageData[src]/ImageData[white]]

下一步是找到数独区域，这样我就可以忽略（掩盖）背景。为此，我使用连通组件分析，并选择具有最大凸面积的组件：

components = 
  ComponentMeasurements[
    ColorNegate@Binarize[srcAdjusted], {"ConvexArea", "Mask"}][[All, 
    2]];
largestComponent = Image[SortBy[components, First][[-1, 2]]]

通过对该图像进行填充，我可以得到数独网格的掩模：

mask = FillingTransform[largestComponent]

现在，我可以使用二阶导数滤波器在两个单独的图像中找到垂直和水平线：

lY = ImageMultiply[MorphologicalBinarize[GaussianFilter[srcAdjusted, 3, {2, 0}], {0.02, 0.05}], mask];
lX = ImageMultiply[MorphologicalBinarize[GaussianFilter[srcAdjusted, 3, {0, 2}], {0.02, 0.05}], mask];

我再次使用连通组件分析从这些图像中提取网格线。网格线比数字要长得多，因此我可以使用卡尺长度仅选择网格线连接的组件。通过按位置对它们进行排序，我可以获得每个垂直/水平网格线的2x10掩码图像：

verticalGridLineMasks = 
  SortBy[ComponentMeasurements[
      lX, {"CaliperLength", "Centroid", "Mask"}, # > 100 &][[All, 
      2]], #[[2, 1]] &][[All, 3]];
horizontalGridLineMasks = 
  SortBy[ComponentMeasurements[
      lY, {"CaliperLength", "Centroid", "Mask"}, # > 100 &][[All, 
      2]], #[[2, 2]] &][[All, 3]];

接下来，我会取每一对垂直/水平网格线，进行膨胀操作，计算像素级别的交点，并计算结果的中心点。这些点就是网格线的交点：

centerOfGravity[l_] := 
 ComponentMeasurements[Image[l], "Centroid"][[1, 2]]
gridCenters = 
  Table[centerOfGravity[
    ImageData[Dilation[Image[h], DiskMatrix[2]]]*
     ImageData[Dilation[Image[v], DiskMatrix[2]]]], {h, 
    horizontalGridLineMasks}, {v, verticalGridLineMasks}];

最后一步是定义两个插值函数来通过这些点进行X/Y映射，并使用这些函数转换图像：

fnX = ListInterpolation[gridCenters[[All, All, 1]]];
fnY = ListInterpolation[gridCenters[[All, All, 2]]];
transformed = 
 ImageTransformation[
  srcAdjusted, {fnX @@ Reverse[#], fnY @@ Reverse[#]} &, {9*50, 9*50},
   PlotRange -> {{1, 10}, {1, 10}}, DataRange -> Full]

所有操作都是基本图像处理函数，因此在OpenCV中也应该可以实现。基于样条的图像转换可能会更难一些，但我认为你真的不需要它。可能在每个单独的单元格上使用您现在使用的透视变换将给出足够好的结果。

- Niki

---

哦，我的天啊！！！！！那太棒了。这真的真的太好了。我会尝试在OpenCV中实现它。希望你能帮我解释一些函数和专业术语的细节……谢谢。 - Abid Rahman K

@arkiaz：我不是OpenCV专家，但如果可以的话，我会尽力帮忙。 - Niki

请问“closing”函数是用来做什么的？我的意思是背后发生了什么？在文档中，它说closing可以去除椒盐噪声？那么closing是低通滤波器吗？ - Abid Rahman K

太棒了！您是从哪里得到将图像亮度归一化的想法，通过除以关闭来实现？我正在尝试提高此方法的速度，因为移动电话上的浮点除法非常缓慢。您有什么建议吗？@AbidRahmanK - 1''

@1*: 我认为它被称为“白色图像调整”。不要问我在哪里读到的，这是一个标准的图像处理工具。这个想法背后的模型很简单：从（Lambertian）表面反射的光量只是表面亮度乘以同一位置上白色物体反射的光量。估计同一位置上白色物体的表面亮度，将实际亮度除以它，就可以得到表面的亮度。 - Niki

Answer 2

Nikie的回答解决了我的问题，但他的答案是用Mathematica写的。因此我想我应该在这里给出它的OpenCV版本。但在实现后，我发现OpenCV代码比nikie的Mathematica代码要长得多。而且，我在OpenCV中找不到nikie使用的插值方法（虽然可以使用scipy来完成，我会在合适的时候提及）。

1.图像预处理（闭运算）

import cv2
import numpy as np

img = cv2.imread('dave.jpg')
img = cv2.GaussianBlur(img,(5,5),0)
gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
mask = np.zeros((gray.shape),np.uint8)
kernel1 = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(11,11))

close = cv2.morphologyEx(gray,cv2.MORPH_CLOSE,kernel1)
div = np.float32(gray)/(close)
res = np.uint8(cv2.normalize(div,div,0,255,cv2.NORM_MINMAX))
res2 = cv2.cvtColor(res,cv2.COLOR_GRAY2BGR)

结果：

2. 查找数独方格并创建遮罩图像

thresh = cv2.adaptiveThreshold(res,255,0,1,19,2)
contour,hier = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)

max_area = 0
best_cnt = None
for cnt in contour:
    area = cv2.contourArea(cnt)
    if area > 1000:
        if area > max_area:
            max_area = area
            best_cnt = cnt

cv2.drawContours(mask,[best_cnt],0,255,-1)
cv2.drawContours(mask,[best_cnt],0,0,2)

res = cv2.bitwise_and(res,mask)

结果：

3. 寻找垂直线

kernelx = cv2.getStructuringElement(cv2.MORPH_RECT,(2,10))

dx = cv2.Sobel(res,cv2.CV_16S,1,0)
dx = cv2.convertScaleAbs(dx)
cv2.normalize(dx,dx,0,255,cv2.NORM_MINMAX)
ret,close = cv2.threshold(dx,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
close = cv2.morphologyEx(close,cv2.MORPH_DILATE,kernelx,iterations = 1)

contour, hier = cv2.findContours(close,cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)
for cnt in contour:
    x,y,w,h = cv2.boundingRect(cnt)
    if h/w > 5:
        cv2.drawContours(close,[cnt],0,255,-1)
    else:
        cv2.drawContours(close,[cnt],0,0,-1)
close = cv2.morphologyEx(close,cv2.MORPH_CLOSE,None,iterations = 2)
closex = close.copy()

结果：

4. 寻找水平线

kernely = cv2.getStructuringElement(cv2.MORPH_RECT,(10,2))
dy = cv2.Sobel(res,cv2.CV_16S,0,2)
dy = cv2.convertScaleAbs(dy)
cv2.normalize(dy,dy,0,255,cv2.NORM_MINMAX)
ret,close = cv2.threshold(dy,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
close = cv2.morphologyEx(close,cv2.MORPH_DILATE,kernely)

contour, hier = cv2.findContours(close,cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)
for cnt in contour:
    x,y,w,h = cv2.boundingRect(cnt)
    if w/h > 5:
        cv2.drawContours(close,[cnt],0,255,-1)
    else:
        cv2.drawContours(close,[cnt],0,0,-1)

close = cv2.morphologyEx(close,cv2.MORPH_DILATE,None,iterations = 2)
closey = close.copy()

结果：

当然，这个不是很好。

5. 寻找网格点

res = cv2.bitwise_and(closex,closey)

结果：

6. 纠正缺陷

在这里，nikie做了某种内插，关于这个我不太了解。我也找不到OpenCV中对应的函数（可能有，但我不知道）。

查看这篇SOF，它解释了如何使用SciPy进行操作，但我不想使用它：OpenCV中的图像转换

因此，在这里，我取每个子正方形的四个角，并对每个角应用透视变换。

为此，首先我们找到质心。

contour, hier = cv2.findContours(res,cv2.RETR_LIST,cv2.CHAIN_APPROX_SIMPLE)
centroids = []
for cnt in contour:
    mom = cv2.moments(cnt)
    (x,y) = int(mom['m10']/mom['m00']), int(mom['m01']/mom['m00'])
    cv2.circle(img,(x,y),4,(0,255,0),-1)
    centroids.append((x,y))

但是得到的质心不会被排序。查看下面的图像以查看它们的顺序： 因此，我们将它们从左到右，从上到下进行排序。

centroids = np.array(centroids,dtype = np.float32)
c = centroids.reshape((100,2))
c2 = c[np.argsort(c[:,1])]

b = np.vstack([c2[i*10:(i+1)*10][np.argsort(c2[i*10:(i+1)*10,0])] for i in xrange(10)])
bm = b.reshape((10,10,2))

现在看一下它们的顺序：

最后，我们应用变换并创建一个大小为450x450的新图像。

output = np.zeros((450,450,3),np.uint8)
for i,j in enumerate(b):
    ri = i/10
    ci = i%10
    if ci != 9 and ri!=9:
        src = bm[ri:ri+2, ci:ci+2 , :].reshape((4,2))
        dst = np.array( [ [ci*50,ri*50],[(ci+1)*50-1,ri*50],[ci*50,(ri+1)*50-1],[(ci+1)*50-1,(ri+1)*50-1] ], np.float32)
        retval = cv2.getPerspectiveTransform(src,dst)
        warp = cv2.warpPerspective(res2,retval,(450,450))
        output[ri*50:(ri+1)*50-1 , ci*50:(ci+1)*50-1] = warp[ri*50:(ri+1)*50-1 , ci*50:(ci+1)*50-1].copy()

结果：

结果与nikie的几乎相同，但是代码长度较大。或许还有更好的方法，但在此之前，这个方法可以正常工作。

敬礼 ARK.

- Abid Rahman K

---

“我宁愿我的应用程序崩溃，也不愿得到错误的答案。”<- 我也完全同意这个观点。 - Viktor Sehr

谢谢，Nikie已经给出了真正的答案。但那是在Mathematica中，所以我只是将其转换为OpenCV。因此，我认为真正的答案已经得到了足够的赞同。 - Abid Rahman K

啊，我没看到你也发了这个问题 :) - Viktor Sehr

我得到了错误： output[ri*50:(ri+1)50-1 , ci50:(ci+1)50-1] = warp[ri50:(ri+1)50-1 , ci50:(ci+1)*50-1].copy TypeError: long() argument must be a string or a number, not 'builtin_function_or_method' - user898678

好的。我发现我的错误：在那一行的末尾，我只写了“.copy”，实际上应该是“.copy()”;-) - user898678

Answer 3

你可以尝试使用基于网格的建模方法来处理任意扭曲。由于数独本身已经是一个网格，所以这应该不难实现。

因此，你可以尝试检测每个3x3子区域的边界，然后单独对每个区域进行扭曲。如果检测成功，它会给出更好的近似结果。

- sietschie

Answer 4

我认为这是一篇很棒的文章，ARK提供了一个很好的解决方案，并且解释得非常清楚易懂。

我曾经也遇到过类似的问题，而且已经把整个东西构建出来了。虽然有些变化（例如xrange改成了range，在cv2.findContours中的参数），但这应该可以直接使用（Python 3.5，Anaconda）。

这是以上元素的综合，加入了一些缺失的代码（例如点的标记）。

'''

https://dev59.com/b2kw5IYBdhLWcg3wBV7j

'''

import cv2
import numpy as np

img = cv2.imread('test.png')

winname="raw image"
cv2.namedWindow(winname)
cv2.imshow(winname, img)
cv2.moveWindow(winname, 100,100)


img = cv2.GaussianBlur(img,(5,5),0)

winname="blurred"
cv2.namedWindow(winname)
cv2.imshow(winname, img)
cv2.moveWindow(winname, 100,150)

gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
mask = np.zeros((gray.shape),np.uint8)
kernel1 = cv2.getStructuringElement(cv2.MORPH_ELLIPSE,(11,11))

winname="gray"
cv2.namedWindow(winname)
cv2.imshow(winname, gray)
cv2.moveWindow(winname, 100,200)

close = cv2.morphologyEx(gray,cv2.MORPH_CLOSE,kernel1)
div = np.float32(gray)/(close)
res = np.uint8(cv2.normalize(div,div,0,255,cv2.NORM_MINMAX))
res2 = cv2.cvtColor(res,cv2.COLOR_GRAY2BGR)

winname="res2"
cv2.namedWindow(winname)
cv2.imshow(winname, res2)
cv2.moveWindow(winname, 100,250)

 #find elements
thresh = cv2.adaptiveThreshold(res,255,0,1,19,2)
img_c, contour,hier = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)

max_area = 0
best_cnt = None
for cnt in contour:
    area = cv2.contourArea(cnt)
    if area > 1000:
        if area > max_area:
            max_area = area
            best_cnt = cnt

cv2.drawContours(mask,[best_cnt],0,255,-1)
cv2.drawContours(mask,[best_cnt],0,0,2)

res = cv2.bitwise_and(res,mask)

winname="puzzle only"
cv2.namedWindow(winname)
cv2.imshow(winname, res)
cv2.moveWindow(winname, 100,300)

# vertical lines
kernelx = cv2.getStructuringElement(cv2.MORPH_RECT,(2,10))

dx = cv2.Sobel(res,cv2.CV_16S,1,0)
dx = cv2.convertScaleAbs(dx)
cv2.normalize(dx,dx,0,255,cv2.NORM_MINMAX)
ret,close = cv2.threshold(dx,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
close = cv2.morphologyEx(close,cv2.MORPH_DILATE,kernelx,iterations = 1)

img_d, contour, hier = cv2.findContours(close,cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)
for cnt in contour:
    x,y,w,h = cv2.boundingRect(cnt)
    if h/w > 5:
        cv2.drawContours(close,[cnt],0,255,-1)
    else:
        cv2.drawContours(close,[cnt],0,0,-1)
close = cv2.morphologyEx(close,cv2.MORPH_CLOSE,None,iterations = 2)
closex = close.copy()

winname="vertical lines"
cv2.namedWindow(winname)
cv2.imshow(winname, img_d)
cv2.moveWindow(winname, 100,350)

# find horizontal lines
kernely = cv2.getStructuringElement(cv2.MORPH_RECT,(10,2))
dy = cv2.Sobel(res,cv2.CV_16S,0,2)
dy = cv2.convertScaleAbs(dy)
cv2.normalize(dy,dy,0,255,cv2.NORM_MINMAX)
ret,close = cv2.threshold(dy,0,255,cv2.THRESH_BINARY+cv2.THRESH_OTSU)
close = cv2.morphologyEx(close,cv2.MORPH_DILATE,kernely)

img_e, contour, hier = cv2.findContours(close,cv2.RETR_EXTERNAL,cv2.CHAIN_APPROX_SIMPLE)

for cnt in contour:
    x,y,w,h = cv2.boundingRect(cnt)
    if w/h > 5:
        cv2.drawContours(close,[cnt],0,255,-1)
    else:
        cv2.drawContours(close,[cnt],0,0,-1)

close = cv2.morphologyEx(close,cv2.MORPH_DILATE,None,iterations = 2)
closey = close.copy()

winname="horizontal lines"
cv2.namedWindow(winname)
cv2.imshow(winname, img_e)
cv2.moveWindow(winname, 100,400)


# intersection of these two gives dots
res = cv2.bitwise_and(closex,closey)

winname="intersections"
cv2.namedWindow(winname)
cv2.imshow(winname, res)
cv2.moveWindow(winname, 100,450)

# text blue
textcolor=(0,255,0)
# points green
pointcolor=(255,0,0)

# find centroids and sort
img_f, contour, hier = cv2.findContours(res,cv2.RETR_LIST,cv2.CHAIN_APPROX_SIMPLE)
centroids = []
for cnt in contour:
    mom = cv2.moments(cnt)
    (x,y) = int(mom['m10']/mom['m00']), int(mom['m01']/mom['m00'])
    cv2.circle(img,(x,y),4,(0,255,0),-1)
    centroids.append((x,y))

# sorting
centroids = np.array(centroids,dtype = np.float32)
c = centroids.reshape((100,2))
c2 = c[np.argsort(c[:,1])]

b = np.vstack([c2[i*10:(i+1)*10][np.argsort(c2[i*10:(i+1)*10,0])] for i in range(10)])
bm = b.reshape((10,10,2))

# make copy
labeled_in_order=res2.copy()

for index, pt in enumerate(b):
    cv2.putText(labeled_in_order,str(index),tuple(pt),cv2.FONT_HERSHEY_DUPLEX, 0.75, textcolor)
    cv2.circle(labeled_in_order, tuple(pt), 5, pointcolor)

winname="labeled in order"
cv2.namedWindow(winname)
cv2.imshow(winname, labeled_in_order)
cv2.moveWindow(winname, 100,500)

# create final

output = np.zeros((450,450,3),np.uint8)
for i,j in enumerate(b):
    ri = int(i/10) # row index
    ci = i%10 # column index
    if ci != 9 and ri!=9:
        src = bm[ri:ri+2, ci:ci+2 , :].reshape((4,2))
        dst = np.array( [ [ci*50,ri*50],[(ci+1)*50-1,ri*50],[ci*50,(ri+1)*50-1],[(ci+1)*50-1,(ri+1)*50-1] ], np.float32)
        retval = cv2.getPerspectiveTransform(src,dst)
        warp = cv2.warpPerspective(res2,retval,(450,450))
        output[ri*50:(ri+1)*50-1 , ci*50:(ci+1)*50-1] = warp[ri*50:(ri+1)*50-1 , ci*50:(ci+1)*50-1].copy()

winname="final"
cv2.namedWindow(winname)
cv2.imshow(winname, output)
cv2.moveWindow(winname, 600,100)

cv2.waitKey(0)
cv2.destroyAllWindows()

- asylumax

Answer 5

我想补充说明，以上方法只适用于数独板直立时，否则高度/宽度（或反之）比率测试很可能会失败，您将无法检测到数独的边缘。（另外，如果图像边框不垂直于图像边界的线，则Sobel操作（dx和dy）仍将按照两个轴的边缘处理线条。）

为了能够检测出直线，您应该在轮廓或基于像素的分析上进行处理，例如轮廓面积/边界矩形面积、左上角和右下角点……

编辑：我成功地通过应用线性回归并检查误差来检查一组轮廓是否构成一条直线。然而，当线的斜率太大（即>1000）或非常接近于0时，线性回归表现不佳。因此，在进行线性回归之前应用上面的比率测试（在最受欢迎的答案中）是合理的，并且对我有用。

- Ali Eren Çelik

Answer 6

为了去除未检测到的角落，我使用了伽马校正，伽马值为0.8。

红色圆圈是为了显示缺失的角落而绘制的。

代码如下：

gamma = 0.8
invGamma = 1/gamma
table = np.array([((i / 255.0) ** invGamma) * 255
                  for i in np.arange(0, 256)]).astype("uint8")
cv2.LUT(img, table, img)

这是针对阿比德·拉赫曼的答案所遗漏的一些角点进行补充。

- Vardan Agarwal