汚いコードさらし。

一応コードも載せておく。

# -*- encoding : utf-8 -*-

from PIL import Image, ImageDraw
from random import randint

img = Image.open('test3.bmp') # 画像ファイルの読み込み
bmp = list(img.getdata()) # 画像の情報を取得

x, y = img.size  # 画像サイズを取得

# imgは2次元配列。座標を表す。画面左上が原点。
img = []
for i in range(x):
    for j in range(y):
        if bmp[(i*x)+j] == 255:
            img.append((j,i))

# ノイズ。ロバスト推定の有効性を確認するため。
'''for i in range(50):
    img.append((randint(0,100), randint(0,100)))
    img.append((randint(400,500), randint(400, 500)))
'''


xSum = 0.0
xSum2 = 0.0
ySum = 0.0
xySum = 0.0
for i in img:
    xSum += i[0]
    xSum2 += i[0]**2
    ySum += i[1]
    xySum += (i[0] * i[1])

# 最小二乗法で直線のパラメータ推定を行う
a = (1.0/(xSum2*(len(img)+1) - xSum*xSum)) * ((len(img)+1) * xySum - xSum * ySum)
b = (1.0/(xSum2*(len(img)+1) - xSum*xSum)) * (-xSum * xySum + xSum2 * ySum)
print a, b

zahyo = []
for i in img:
    zahyo.append((i[0], a*i[0] + b))

# ロバスト推定で直線のパラメータ推定を行う
# TukeyのBiweight推定法
W = 1000.0  # 許容誤差
wList = []
for i in img:
    d = i[1] - (a*i[0] + b)
    if -W < d and d < W:
        wList.append((1-((d/W)**2))**2)
    else:
        wList.append(0)
wxSum = 0.0
wxSum2 = 0.0
wySum = 0.0
wxySum = 0.0
wSum = 0.0
count = 0
while count < len(img):
    wxSum += wList[count]*img[count][0]
    wxSum2 += wList[count]*(img[count][0]**2)
    wySum += wList[count]*img[count][1]
    wxySum += wList[count]*img[count][0]*img[count][1]
    wSum += wList[count]
    count += 1

aa = (1.0/(wxSum2*wSum - wxSum*wxSum)) * (wSum * wxySum - wxSum * wySum)
bb = (1.0/(wxSum2*wSum - wxSum*wxSum)) * (-wxSum * wxySum + wxSum2 * wySum)
print aa, bb

zahyo2 = []
for i in img:
    zahyo2.append((i[0], aa*i[0] + bb))

f = open("original", "w")
for i in img:
    f.write(str(i[0]) + " " + str(-i[1]) + "\n")
f.close()

f = open("last-square-method", "w")
for i in zahyo:
    f.write(str(i[0]) + " " + str(-i[1]) + "\n")
f.close()

f = open("robust", "w")
for i in zahyo2:
    f.write(str(i[0]) + " " + str(-i[1]) + "\n")
f.close()

#座標の比較(同一x座標から得られるy座標の違い)
#for i in range(len(zahyo)):
#    print img[i], zahyo[i], zahyo2[i]