グラフ読み取り

PyQtでグラフを読み取れるアプリを作った。
グラフの画像をCtrl+Vでペーストし、その座標をクリックすると画素値からグラフの値を教えてくれる。

最初に右クリックで左下と右上の座標を教える必要がる。
なんでこんなんがいるかって、世の中紙でデータを保存する人がまだいるということである。

(OpenCVは使っていない)

from PyQt4 import QtGui, QtCore, QtCore
import sys
import cv2
import numpy as np

def convert_array_to_qimg(cv_img):
    """
    this is based on https://stackoverflow.com/questions/18676888/how-to-configure-color-when-convert-numpy-array-to-qimage
    """
    height, width, bytesPerComponent = cv_img.shape
    bytesPerLine = bytesPerComponent * width;
    #cv2.cvtColor(cv_img, cv2.CV_BGR2RGB)
    return QtGui.QImage(cv_img.data, width, height, bytesPerLine, QtGui.QImage.Format_RGB888)

class GraphCoordinate:
    def __init__(self):
        pass

    def setBottomPixelCoordinate(self, x1, y1):
        """
        x1 and y1 are pixel
        """
        self.x1 = x1
        self.y1 = y1

    def setBottomGraphCoordinate(self, X1, Y1):
        self.X1 = X1
        self.Y1 = Y1
        print "(X1, Y1) = (%.2f, %.2f)" % (self.X1, self.Y1)


    def setTopPixelCoordinate(self, x2, y2):
        self.x2 = x2
        self.y2 = y2

    def setTopGraphCoordinate(self, X2, Y2):
        self.X2 = X2
        self.Y2 = Y2
        print "(X2, Y2) = (%.2f, %.2f)" % (self.X2, self.Y2)

    def getGraphCoordinate(self, x, y):
        w = x - self.x1 
        h = self.y1 - y
        W = self.x2 - self.x1 
        H = self.y1 - self.y2 
        X = self.X1 + w * (self.X2 - self.X1) / float(W)
        Y = self.Y1 + h * (self.Y2 - self.Y1) / float(H)
        return X, Y 



class MainWindow(QtGui.QMainWindow):
    def __init__(self):
        super(MainWindow, self).__init__()
        self.image_frame = QtGui.QLabel()
        fname = r"sample08.png"
        img = cv2.imread(fname)
        #img = QtGui.QImage(fname)
        qimg = convert_array_to_qimg(img)
        self.image_frame.setPixmap(QtGui.QPixmap.fromImage(qimg))
        self.clip = QtGui.QApplication.clipboard()
        
        self.setCentralWidget(self.image_frame)
        self.move(30,30)

        self.counter = 0
        self.coord = GraphCoordinate()
 
    def keyPressEvent(self, e):

        if (e.modifiers() & QtCore.Qt.ControlModifier):
            #selected = self.table.selectedRanges()
            if e.key() == QtCore.Qt.Key_V:#past
                qimg = self.clip.image()
                self.displayImage(qimg)

            elif e.key() == QtCore.Qt.Key_C: #copy
                pass

    def mousePressEvent(self, e):
        #print e.x(), e.y()
        if e.button() == QtCore.Qt.RightButton:
            if self.counter % 2 == 0:
                text, ok = QtGui.QInputDialog.getText(self, 'Input Dialog', 
                                                        'Enter X1, Y1 (comma separated):')
                X1, Y1 = [float(v) for v in text.split(',')]
                if ok:
                    self.counter += 1
                    self.coord.setBottomPixelCoordinate(e.x(), e.y())
                    self.coord.setBottomGraphCoordinate(X1, Y1)


            else:
                text, ok = QtGui.QInputDialog.getText(self, 'Input Dialog', 
                                                        'Enter X2, Y2 (comma separated):')
                X2, Y2 = [float(v) for v in text.split(',')]
                if ok:
                    self.coord.setTopPixelCoordinate(e.x(), e.y())
                    self.coord.setTopGraphCoordinate(X2, Y2)
                    self.counter = 0

        elif e.button() == QtCore.Qt.LeftButton:
            try:
                x, y = self.coord.getGraphCoordinate(e.x(), e.y())
                print "%.5f %.5f" % (x, y)

            except AttributeError as e:
                print "please set the graph coordinate by right clicks."


    def displayImage(self, qimg):
        self.image_frame.setPixmap(QtGui.QPixmap.fromImage(qimg))

 
def main(args):
    app = QtGui.QApplication(sys.argv)
    form = MainWindow()
    form.show()
     
    sys.exit(app.exec_())
 
if __name__ == "__main__":
    main(sys.argv)

今さら聞けないRichardson-Lucy入門

デコンボリューションがとても魅力的だ。これを使えばレンズの収差やボケ、手ブレによる画像劣化が除去できるという夢の技術だ。

いろんな論文を読むとかならずといっていいほど以下の2つの手法が挙げられている。

1.Wiener Filter

Wikipediaに導出方法が載っている。

2.Richardson-Lucy Deconvolution

こちらもWikipediaに載っている。Wikipediaを読んでもわからなかったので、勉強した内容を少し解説してみる。

わかりやすいように1次元で話を進める。

をピクセルiの真の値(全くボケていないときの値)

をピクセルiで観測された値

とする。全くボケていないときはこの2つが一致する。でも今回はPSF(Point Spread Function)によってボケる。ボケがあるとピクセルjへ行く光は広がりをもってしまうので、ピクセルiへ少しおすそ分けしてしまう。
本来jへいくはずの光のボケでのih量を

とする。Richardson-Lucyではこのp(i,j)は既知だとするイメージは下の図。

ここでZ(i,j)なるものを導入する。Z(i,j)は期待値がP(i,j)λjのポアッソン分布に従いう確率変数。

これはピクセルiに届く光が確率的に変化し、期待値はp(i, j)λjということを意味する。

Wikipekiaにこの量は出てこないけど、元論文などを読むと出てくる。
もともとデコンボリューションはCTとかで研究されていて、CTは原子の崩壊(ガンマ線?)を捉えるので、観測される値が確率的になる、という経緯がある。
(一様な確率で起こる事象がある期間で観測される回数はポアッソン分布になる。詳しくはWikipediaをみてください。)

さて、点が1点だけなら復元も簡単かもしれんけど、周りの点の広がりが重なって観測されるので問題が難しくなる。例えば下の図のように重なるのでピクセルiで値は, 周辺のピクセルからの寄与を足しあわせたのものになる。

ここでzは確率的な値なので、観測結果と期待値から予測することで話を前に進めよう。
ピクセルiで観測された値をつかうと、その内ピクセルjからきた値は、割合で考えると下のような式になる。右辺の分母がiへくる光の総和で、分子がその内jからきたものを表す。

zはもともとjにあったものが周辺のピクセルiに散らばったものとかんがえることができるので、集めてやると真の値λjが分かる。真のzはわからないので推定値を使うとλを推定できてる。

この式にさっきのzの推定値を代入すると

このλjは真の値なので、知り得ない。そこで前回推定した値を使って逐次的に求めることにすれば

というふうになる。

実際に実装するときはWikipediaのようにコンボリューションの関数を使って書いたほうが楽だった。


(2014/11/15　式間違ってましたので修正しました。前に書いてたやつではPSFが対称じゃないとうまくいかなかったorz)

参考
https://en.wikipedia.org/wiki/Richardson%E2%80%93Lucy_deconvolution
http://people.csail.mit.edu/dgreensp/may/deconvlucy.html

Juliaで自作VectorとArray{Float64,1}の速度比較

Juliaの勉強を始めた。
Juliaは数値計算の関数を大量に有しているけど、可変長のArrayをもとにしていて、3x3の行列や3x1ベクトルの単純なものだと自分で定義したほうが高速なのだろうか?という疑問が湧いたので実験してみた。
もしかしたらJuliaは裏でものすごいことをやっていてるかもしれんし。

3x1のベクトルVecを以下のように定義した。

immutable Vec
    x::Float64
    y::Float64
    z::Float64
end
dot(a::Vec, b::Vec) = (a.x*b.x + a.y*b.y + a.z*b.z)
cross(a::Vec, b::Vec) = Vec(a.y*b.z - a.z*b.y, a.z*b.x-a.x*b.z, a.x*b.y-a.y*b.x)

ただ外積やらなんやらを計算する関数を作った。

function cross_loop(N::Int, a::Vec, b::Vec)
    tmp = 0.0
    for i in 1:N
        v = cross(a,b)
        tmp += dot(v,v)
    end 
    tmp/N
end

function cross_loop(N::Int, a::Array{Float64,1}, b::Array{Float64,1})
    tmp = 0.
    for i in 1:N
        v = Base.cross(a,b)
        tmp += Base.dot(v,v) 
    end
    tmp / N
end

Arrayの方はデフォルトで入っているBase.dotやcrossをつかった。

N=10^7くらいにして計算
速度は@timeで計測した。1回目はコンパイルに時間を取られるので2回目の時間を使用。

N = 10^7
a = Vec(1,2,3)
b = Vec(4,5,6)
v = [1.,2.,3.]
u = [4.,5.,6.]

@time cross_loop(N, a, b)
@time cross_loop(N, u, v)

結果はダントツで自作のほうが速い。

elapsed time: 0.013432361 seconds (96 bytes allocated)
elapsed time: 2.079223219 seconds (1280000096 bytes allocated, 43.49% gc time)

ちなみにVecをimmutableでなくtypeで定義すると、それだけでかなり遅くなる。GCが働いているようだ。

elapsed time: 0.419364716 seconds (320000096 bytes allocated, 53.07% gc time)

結論:3次元に特化したベクトル計算では自分で定義してやったほうが速い。