鱼眼图像提取有效区域（Python实现） 改进的线扫描法

以下节选自论文。

线扫描法

原理很直观，它的思路是将原始图像灰度化，然后在灰度图上，分别使用水平和垂直两对直线，从上下左右四个方向夹逼圆形有效区域，获得四个切点的坐标，之后利用这两对坐标，计算出圆心坐标和半径大小。在扫描线移动的过程中，只要同一扫描线上的像素最大亮度值与最小亮度值之差大于设定的阈值，便认为该扫描线已经与圆形有效区域相切，最大亮度值的像素点即为切点（这是论文原话，但我个人觉得表述不太完整，应该是在直线逼近圆的过程中，第一个满足“同一扫描线上的像素最大亮度值与最小亮度值之差大于设定的阈值”这一条件的点）。   

假设两对切点的图像坐标值分别为：

则圆心和半径的图像坐标值可以这样得到：

改进的线扫描法

传统的线扫描法仅仅从上下左右四个方向逼近圆形区域，如果圆形区域刚好在这四个顶点附近存在较多噪声点，必然会导致该方法产生误差，针对这一情况，本文使用成对的不同倾斜角的直线，按照与传统线扫描法相同的判断方法：在扫描线移动的过程中，只要同一扫描线上的像素最大亮度值与最小亮度值之差大于设定的阈值，便认为该扫描线已经与圆形有效区域相切，最大亮度值的像素点即为切点。

从各个方向对圆形区域夹逼，求取切点坐标，并对这些切点进行筛选，剔除一些无效切点，利用剩余的切点坐标进行圆拟合，进而计算出圆心坐标和半径。

首先，令直线对的倾斜角从0开始，以Π/2n递增（n为设定的常数，越大越精确），分别从圆形区域两侧逼近，从而找到各对切点的坐标值。同一条直线所获得的两个切点称为一对切点。

连接每对切点，其距离的一半为这对切点对整个圆的半径的估测值。求所有估测值的平均值。由于所求切点接近圆的内部（论文中是外部，这一点跟我所观察的结果不一样），故将估测半径小于平均值的点对剔除。

之后，按照下面的方法（Kasa圆拟合法）进行圆拟合，得到圆心和半径。

代码

config.py

就是一些全局变量。

class Config(object):
    pi=3.1415926535
    n=10
    #图片大小
    l=0
    w=0
    #灰度阈值，用于切点判断
    thre=50
    #路径
    path='D:/TheMoth/fisheye.jpg'

main.py

from PIL import Image
from torchvision.transforms import ToTensor,ToPILImage
import matplotlib.pyplot as plt
import numpy as np
import math
from config import Config
import cv2

#设置
opt = Config()
#下面四个数组临时存储切点对，其中(sx[i],sy[i])与(tx[i],ty[i])表示同一直线所切的一对切点
sx=[]
sy=[]
tx=[]
ty=[]
#切点对之间的距离
dis=[]
#存储最终的切点集
xx=[]
yy=[]

#绘制图像，同时获取图像长宽
def open_img():
    '''打开图像'''
    img=Image.open(opt.path)
    #绘制图像
    plt.imshow(img)
    #将图片转化为numpy数组
    data=np.asarray(img)
    #维度转置
    data=np.transpose(data,[2,0,1])
    #获取图像的长宽
    opt.l=len(data[0])
    opt.w=len(data[0][0])
    return data

#获取灰度图像
def open_cvmg():    
    src=cv2.imread(opt.path, cv2.IMREAD_GRAYSCALE)
    data=np.asarray(src)
    return data
pic_one_channel=open_cvmg()

直线扫描

在判断切点的时候，需要扫描整条直线，为保证扫描的精确性，采用计算机图形学中的DDA算法扫描直线。

#DDA算法扫描,(x0,y0)(x1,y1)为起始点坐标
def scanline(x0,y0,x1,y1):
    if(x0==x1 and y0==y1):
        return 0,0,0
    dx=x1-x0
    dy=y1-y0
    x=float(x0)
    y=float(y0)    
    if(abs(dx)>abs(dy)):
        eps=abs(dx)
    else:
        eps=abs(dy)    
    xIncre=float(dx)/float(eps)
    yIncre=float(dy)/float(eps)
    
    maxc=0       #最大的灰度值
    minc=500     #最小的灰度值
    is_cut=0     #是否是切点
    #切点坐标
    r=0
    c=0
    
    for k in range(0,eps+1):
    
        #是否出界
        if(int(x)>=opt.l or int(y)>=opt.w):
            continue
        tem=pic_one_channel[int(x)][int(y)]
        if maxc<tem:
            maxc=tem
            r=int(x)
            c=int(y)
        minc=min(minc,tem)
        
        x+=xIncre
        y+=yIncre
    if(maxc-minc>=opt.thre):
        is_cut=1
    return is_cut,r,c

垂直，水平线用DDA不太方便，这里单独扫描。

#水平&垂直线扫描
def prescan():
    l=opt.l
    w=opt.w
    #从上到下
    for i in range(l):
        is_cut=0
        minc=500
        maxc=0
        r=0
        c=0
        for j in range(w):
            tem=pic_one_channel[int(i)][int(j)]
            if maxc<tem:
                maxc=tem
                r=int(i)
                c=int(j)
            minc=min(minc,tem)
        if(maxc-minc>=opt.thre):
            is_cut=1
        if(is_cut):
            sx.append(r)
            sy.append(c)
            break
        
    #从下到上
    for i in range(l-1,0,-1):
        is_cut=0
        minc=255
        maxc=0
        r=0
        c=0
        for j in range(w):
            tem=pic_one_channel[int(i)][int(j)]
            if maxc<tem:
                maxc=tem
                r=int(i)
                c=int(j)
            minc=min(minc,tem)
        if(maxc-minc>=opt.thre):
            is_cut=1
        if(is_cut):
            tx.append(r)
            ty.append(c)
            break
    
    #从左到右
    for j in range(w):
        is_cut=0
        minc=255
        maxc=0
        r=0
        c=0
        for i in range(l):
            tem=pic_one_channel[int(i)][int(j)]
            if maxc<tem:
                maxc=tem
                r=int(i)
                c=int(j)
            minc=min(minc,tem)
        if(maxc-minc>=opt.thre):
            is_cut=1
        if(is_cut):
            sx.append(r)
            sy.append(c)
            break    
    
    #从右到左
    for j in range(w-1,0,-1):
        is_cut=0
        minc=255
        maxc=0
        r=0
        c=0
        for i in range(l):
            tem=pic_one_channel[int(i)][int(j)]
            if maxc<tem:
                maxc=tem
                r=int(i)
                c=int(j)
            minc=min(minc,tem)
        if(maxc-minc>=opt.thre):
            is_cut=1
        if(is_cut):
            tx.append(r)
            ty.append(c)
            break  

绘制圆形

#参数方程法绘制圆形，r为半径，(u,v)为圆心
def draw_circle(r,u,v):
    theta = np.arange(0, 2*np.pi, 0.01)
    x = u + r * np.cos(theta)
    y = v + r * np.sin(theta)
    #这里圆是红色的
    plt.plot(x,y,'r')

主函数

if __name__ == '__main__':
    
    #读取彩色图像
    open_img()
    #先扫描水平&垂直的直线
    prescan()
    
    l=opt.l
    w=opt.w
    th=0
    
    #直线每次变化的角度
    step=opt.pi/(opt.n*2)
    
    for i in range(opt.n*2):
        
        #斜率
        k=math.tan(th)
        #(fr,fc)为扫描过程中，第一个满足切点条件的点
        #(lr,lc)为扫描过程中，最后一个满足切点条件的点
        #这里斜率相同的两条直线并不是分别从两端向圆逼近，而是穿过整个圆。第一个与最后一个满足切点条件的点为切点
        fr=-1
        fc=-1
        lr=-1
        lc=-1
        if(k<10 and k != 0):
            if(k>0):
                #平移直线
                for x0 in range(l,0,-1):
                    if(int(k*(l-x0))<=w):
                        is_cut,r,c=scanline(x0,0,l,int(k*(l-x0)))
                    elif(int(w/k+x0)<=l):
                        is_cut,r,c=scanline(x0,0,int(w/k+x0),w)
                    if(is_cut and fr==-1):
                        fr=r
                        fc=c
                    if(is_cut):
                        lr=r
                        lc=c
                for y0 in range(w):
                    if(int(k*l+y0)<=w):
                        is_cut,r,c=scanline(0,y0,l,int(k*l+y0))
                    elif(int((w-y0)/k)<=l):
                        is_cut,r,c=scanline(0,y0,int((w-y0)/k),w)
                    if(is_cut and fr==-1):
                        fr=r
                        fc=c
                    if(is_cut):
                        lr=r
                        lc=c
            if(k<0):
                for x0 in range(l):
                    if(int(-k*x0)<=w):
                        is_cut,r,c=scanline(x0,0,0,int(-k*x0))
                    elif(int(w/k+x0)<=l):
                        is_cut,r,c=scanline(x0,0,int(w/k+x0),w)
                    if(is_cut and fr==-1):
                        fr=r
                        fc=c
                    if(is_cut):
                        lr=r
                        lc=c
                for y0 in range(w):
                    if(int(y0-k*l<=w)):
                        is_cut,r,c=scanline(0,int(y0-k*l),l,y0)
                    elif(int((w-y0)/k+l)<=l):
                        is_cut,r,c=scanline(int((w-y0)/k+l),w,l,y0)
                    if(is_cut and fr==-1):
                        fr=r
                        fc=c
                    if(is_cut):
                        lr=r
                        lc=c
        if lc!=-1:
            sx.append(fr)
            sy.append(fc)
            tx.append(lr)
            ty.append(lc)
        th+=step
    
    mid_dis=0
    n=len(sx)
    
    #求距离与平均距离
    for i in range(n):
        d=math.sqrt((sx[i]-tx[i])*(sx[i]-tx[i])+(sy[i]-ty[i])*(sy[i]-ty[i]))
        dis.append(d)
    for i in range(n):
        mid_dis+=dis[i]
    mid_dis/=n
    
    #保留距离大于平均距离的点对
    for i in range(n):
        if(dis[i]>=mid_dis):
            xx.append(sx[i])
            xx.append(tx[i])
            yy.append(sy[i])
            yy.append(ty[i])
    
    #圆拟合法
    m=len(xx)
    a=np.ones((m,3))
    b=np.ones((m,1))
    for i in range(m):
        a[i][0]=xx[i]
        a[i][1]=yy[i]
        a[i][2]=1
    for i in range(m):
        b[i][0]=xx[i]*xx[i]+yy[i]*yy[i]
    p=np.matmul(np.linalg.pinv(a),b)
    p1=p[0][0]
    p2=p[1][0]
    p3=p[2][0]
    u0=p1/2
    v0=p2/2
    R=math.sqrt((p1*p1+p2*p2)/4+p3)
    
    #将圆绘制在原图上
    draw_circle(R, v0, u0)

提取效果

红色区域为提取区域，可见效果还不错。