汚いコードさらし。
一応コードも載せておく。
# -*- 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]