最近再做一个卫星遥感监测台风的项目，其中有个算法要提取台风中心一定半径的一圈像素然后做一些统计。查了查，原来这个算法还有个名字叫做中圆点算法（Midpoint circle algorithm）。

中圆点算法，说简单点就是设定一个半径画个圆，然后看看哪些像素在这个圆圈上。

其实利用 PIL 可以很简单的实现中圆点算法的功能，找出哪些像素组成了圆圈。

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib
import matplotlib.pyplot as plt
from PIL import Image, ImageDraw


def get_midpoint_circle_pixels(size=100, radius=30, isPlot=False):
    '''get midpoint circle pixels easily with PIL
    
    Keyword Arguments:
        size {number} -- image size (default: {100})
        radius {number} -- circle radius, no more than size/2 (default: {30})
        isPlot {bool} -- show image or not (default: {False})
    
    Returns:
        x -- x coordinate
        y -- y coordinate
        points -- [[x[i], y[i]] for i in range(len(x))]
    '''

    # create new image, size x size pixels, 1 bit per pixel
    image = Image.new('1', (size, size))
    draw = ImageDraw.Draw(image)

    # actually a circle was drawn
    lower_left = size/2 - radius
    upper_right = size/2 + radius
    draw.ellipse((lower_left, lower_left, upper_right, upper_right), outline='white')

    # image.show()

    # get pixel values
    color = list(image.getdata())
    color = np.reshape(color, (size, size))

    # get pixel coordinate with color='white'
    grid_x, grid_y = np.mgrid[:size, :size]
    x = grid_x[color==255].ravel()
    y = grid_y[color==255].ravel()
    points = [[x[_], y[_]] for _ in range(len(x))]

    if isPlot:
        plt.pcolor(color)
        plt.axis('scaled')
        plt.hold(True)
        plt.plot(x, y, 'wo')
        plt.show()

    return x, y, points


if __name__ == '__main__':
    size = 100
    x1, y1, _ = get_midpoint_circle_pixels(size=size, radius=10)
    x2, y2, _ = get_midpoint_circle_pixels(size=size, radius=25)
    x3, y3, _ = get_midpoint_circle_pixels(size=size, radius=40)
    plt.hold(True)
    plt.plot(x1, y1, 'ro')
    plt.plot(x2, y2, 'go')
    plt.plot(x3, y3, 'bo')
    plt.axis('scaled')
    plt.xlim(0, size)
    plt.ylim(0, size)
    plt.show()

调用上面的方法，就得到的一组坐标点，但是这些点并没有按照画圆的方向排列……

x, y, points = get_midpoint_circle_pixels(size=10, radius=3, isPlot=True)
print x
print y
print points

[2 2 2 2 2 3 3 3 3 4 4 5 5 6 6 7 7 7 7 8 8 8 8 8]
[3 4 5 6 7 2 3 7 8 2 8 2 8 2 8 2 3 7 8 3 4 5 6 7]
[[2, 3], [2, 4], [2, 5], [2, 6], [2, 7], [3, 2], [3, 3], [3, 7], [3, 8], [4, 2], [4, 8], [5, 2], [5, 8], [6, 2], [6, 8], [7, 2], [7, 3], [7, 7], [7, 8], [8, 3], [8, 4], [8, 5], [8, 6], [8, 7]]

不过排列起来也并不困难，做法是这样的：

1. 先预定输出排列好的列表为 ring；
2. 从 points 里找到 center top 的那一个点，放入 ring；
3. 然后按照顺时针方向找到下一个点追加到 ring 里面，直到结束；
   • 这个点必然在上一个点周围的 8 个点之中；
   • 按照先上/右/下/左然后右上/右下/左下/左上顺序搜索；
   • 只要从这 8 个点里面找到一个点，当它在 points 里但是尚未在 ring 里，那就找到了；
   • 特别注意，第一次搜索时，只需要找到唯一一个 next 点，这才能确保 ring 的排序方向。

# sort points clockwise
ymax = max(y)
xmid = int(np.median(x))
first_point = [xmid, ymax]
ring = [first_point,]
for _ in range(len(points)):
    x, y = ring[-1]
    # search in the 8 points around: firstly top/right/bottom/left and secondly ur/lr/ll/ul
    possible = [
        [x, y+1], [x+1, y], [x, y-1], [x-1, y],
        [x+1, y+1], [x+1, y-1], [x-1, y-1], [x-1, y+1]]
    for p in possible:
        if (p in points) and (p not in ring):
            ring.append(p)
            # at the first point, find the only one next point: this will insure the direction
            if len(ring) == 2:
                break
x = [ring[_][0] for _ in range(len(ring))]
y = [ring[_][1] for _ in range(len(ring))]
points = [[x[_], y[_]] for _ in range(len(x))]

这样就完成了 x, y, points 的排列。这几句话可以追加到 get_midpoint_circle_pixels() 方法里，直接输出排序结果。最后画个图出来验证一下排序结果是否正确：

size = 50
for r in range(1, 23):
    x, y, points = get_midpoint_circle_pixels(size=size, radius=r)
    plt.plot(x, y, '-s', label=str(r))
    plt.axis('scaled')
    plt.xlim(0, size)
    plt.ylim(0, size)
    # plt.legend(loc='best')
plt.show()

结果是正确的。