这次还是一个卫星项目，问题是这样的：

在一个栅格图像上，有一个点的位置上有一个向量，需要计算的是这个向量所在的直线所经过的像元有哪些。

问题看上去简单，但似乎实现起来还有点麻烦。

还好我们有 PIL 可以用，问题就简化成只要找到这条直线的两个端点就行了，思路是这样的：

1. 这两个端点必然在栅格图像最外圈上，只要找到这上面的两个点就行了；
2. 首先提取最外圈所有像素的坐标，然后计算每个点与给定点之间的向量，计算该向量与所给向量之前相互平行的一个度量；
3. 找到最为相互平行的两个度量所对应的点，那就基本大功告成了。

最后利用 PIL，把找到的两个点在栅格图像上画出直线，找到直线经过的像素点就完成了。

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

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


def vector_across_image(p=(10, 30), u=-1, v=2, size=100, isPlot=False):
    """vector line across a square image
    
    A point on the image with a direction like u,v components, make a vector,
    which across the image, result as a line of pixels
    
    Keyword Arguments:
        p {tuple} -- position of the point on the image (default: {(10, 30)})
        u {number} -- u component (default: {-1})
        v {number} -- v component (default: {2})
        size {number} -- image size (default: {100})
        isPlot {bool} -- view the result or not (default: {False})
    
    Returns:
        points -- pixels on the line
    """
    # make a square
    x, y = np.mgrid[:size, :size]
    # take the outline of the square
    outline = zip(x[0, :], y[0, :]) + zip(x[-1, :], y[-1, :]) \
            + zip(x[:, 0], y[:, 0]) + zip(x[:, -1], y[:, -1])
    outline = np.array(outline)

    # find the two interception points on the outline, i.e. p0 and p1
    # delta = [(_[0]-p[0])*v - (_[1]-p[1])*u for _ in outline]
    delta = (outline[:, 0] - p[0]) * v - (outline[:, 1] -p[1]) * u
    delta = np.abs(delta)
    ascding = np.argsort(delta)
    p0 = outline[ascding[0]]
    p1 = outline[ascding[1]]

    # get the line using PIL
    image = Image.new('1', (size, size))
    draw = ImageDraw.Draw(image)
    draw.line((p0[0], p0[1], p1[0], p1[1]), fill='white')
    color = list(image.getdata())
    color = np.reshape(color, (size, size))
    px = y[color==255].ravel()  # !!! x-->y
    py = x[color==255].ravel()  # !!! y-->x
    points = zip(px, py)

    if isPlot:
        plt.subplot(211)
        plt.plot(px, py, 'rs')
        plt.plot(p[0], p[1], 'ko', ms=20)
        plt.plot([p0[0], p1[0]], [p0[1], p1[1]], 'b-')
        plt.axis('scaled')
        plt.axis([0, size, 0, size])
        plt.subplot(212)
        plt.imshow(color)
        plt.show()

    return points, color


def test():
    points, color = vector_across_image(p=(10, 30), u=-1, v=2, size=100, isPlot=True)
    # print color


if __name__ == '__main__':
    test()
    # from timeit import timeit
    # t = timeit('vector_across_image()', 'from __main__ import vector_across_image', number=140*140)
    # print(t)