end0tknr's kipple - web写経開発

太宰府天満宮の狛犬って、妙にカワイイ

opencv.threshold() によるグレースケール画像の二値化

間取り図(画像)にある淡い色を変換(削除)したかったので

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
import cv2
import sys

def main():
    org_file = sys.argv[1]

    # 強制的にグレー画像として読む
    # http://opencv.jp/opencv-2.1/cpp/reading_and_writing_images_and_video.html
    img = cv2.imread(org_file, 0)

    ret, img_new = cv2.threshold(img,
                                 200,               # 閾値
                                 256,               # 画素値の最大値
                                 cv2.THRESH_BINARY) # 2値化type

    cv2.imwrite('org.png', img)
    cv2.imwrite('new.png', img_new)


if __name__ == '__main__':
    main()

↑こう書くと↓こう変換されます

f:id:end0tknr:20171103202830p:plain:w260 f:id:end0tknr:20171103202833p:plain:w260

opencv.medianBlur() で、中央値フィルタ

「間取り図に対して実行することで、ぼかし、その後、減算すれば、細線を楽に消去できるかも」と思いましたが、期待はずれでした。

というより、そりゃそうだ、という感じ

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
import cv2
import sys

def main():
    org_file = sys.argv[1]

    # 強制的にグレー画像として読む
    # http://opencv.jp/opencv-2.1/cpp/reading_and_writing_images_and_video.html
    img = cv2.imread(org_file, 0)

    # 中央値フィルタ
    img_new = cv2.medianBlur(img, ksize=5)
    cv2.imwrite('new.png', img_new)
#    cv2.imwrite('org.png', img)


if __name__ == '__main__':
    main()

↑こう書くと↓こう変換されます

f:id:end0tknr:20171103191712p:plain:w270 f:id:end0tknr:20171103191702p:plain:w270

画像をグレースケールで開き、opencv or numpy or matplotlib.pylab でヒストグラム

画像ファイルを、線画と写真に分類したいと思ったのが、きっかけ

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
from PIL import Image
import cv2
import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb
import numpy as np
import sys

def main():
    org_file = sys.argv[1]

    hist_by_opencv(org_file)     # opencvで算出(高速)
    hist_by_numpy(org_file)      # numpy で算出
    hist_by_matplotlib(org_file) # matplotlib.pylab で算出&描画


def hist_by_opencv(org_file):
    # 強制的にグレー画像として読む
    # http://opencv.jp/opencv-2.1/cpp/reading_and_writing_images_and_video.html
    img = cv2.imread(org_file, 0)

    histr = cv2.calcHist([img],[0],None,[256],[0,256])

    # 線画と写真を判定する為、試しに値のある度数をカウント
#    print bins_size_has_val(histr)

    bins  = range(256)
    
    plb.figure()  # 新規図
    plb.title("by opencv calcHist()")
    # ヒストグラム 描画
    # https://matplotlib.org/examples/pylab_examples/index.html
    plb.step(bins, histr)

    plb.savefig( 'hist_by_opencv.png' )


def hist_by_numpy(org_file):
    # 画像をグレースケールで開き、配列に読み込む
    img = np.array(Image.open(org_file).convert('L'))
    # ヒストグラム算出
    histr, bins = np.histogram(img, bins=256)

    plb.figure()  # 新規図
    plb.title("by numpy histogram()")
    # ヒストグラム 描画
    # https://matplotlib.org/examples/pylab_examples/index.html
    plb.step(bins[:-1], histr)

    plb.savefig( 'hist_by_numpy.png' )

def hist_by_matplotlib(org_file):
    # 画像をグレースケールで開き、配列に読み込む
    img = np.array(Image.open(org_file).convert('L'))

    plb.figure()                    # 新規図
    plb.title("by matplotlib.pylab hist()")
    plb.hist(img.flatten(),256)     # ヒストグラム算出 & 描画
    plb.savefig( 'hist_by_matplotlib.png' )


def bins_size_has_val(histr):
    bins_size = 0
    for val in histr:
        if val[0] > 0:
            bins_size += 1
    return bins_size

if __name__ == '__main__':
    main()

↑こう書くと↓このように3種類、算出&表示されます

f:id:end0tknr:20171103163202j:plain:w300

f:id:end0tknr:20171103163258p:plain:w260 f:id:end0tknr:20171103163308p:plain:w260 f:id:end0tknr:20171103163317p:plain:w260

または↓こう算出&表示されます

f:id:end0tknr:20171103163626p:plain:w300

f:id:end0tknr:20171103163644p:plain:w260 f:id:end0tknr:20171103163700p:plain:w260 f:id:end0tknr:20171103163706p:plain:w260

installing opencv-python and tutorial

opencv-python 3.3.0.10 : Python Package Index

実践コンピュータビジョン サンプルプログラムOpenCVのソースに付属するサンプル( opencv/samples/python )を実行

install opencv-python

# yum install openssl-devel*

$ wget https://www.python.org/ftp/python/2.7.14/Python-2.7.14.tar.xz
$ tar -xvf Python-2.7.14.tar.xz
$ cd Python-2.7.14
$ ./configure
$ make
$ make test
# make install

$ /usr/local/bin/python --version
Python 2.7.14

# curl -kL https://bootstrap.pypa.io/get-pip.py | /usr/local/bin/python
# pip install opencv-python
  :
Installing collected packages: numpy, opencv-python
Successfully installed numpy-1.13.3 opencv-python-3.3.0.10

実践コンピュータビジョン の10章

実践コンピュータビジョン サンプルプログラム

画像の読込みと、形状(size)取得、グレースケール化

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
import cv2

def main():
    # 画像load
    im = cv2.imread('empire.jpg')

    # 画像size
    h,w = im.shape[:2]
    print h,w

    # 別形式で保存
    cv2.imwrite('result.png',im)

    # グレースケール化
    gray = cv2.cvtColor(im,cv2.COLOR_BGR2GRAY)
    cv2.imwrite('result_gray.png',gray)

if __name__ == '__main__':
    main()

↑こう書くと、↓こう変換

f:id:end0tknr:20081110180136j:plain:w250  f:id:end0tknr:20171029081155p:plain:w250

flood fillによる塗りつぶし

  • flood fillは領域判定にも使用されるらしい
  • cv2.floodFill()の仕様は、ちゃんとは理解していません
#!/usr/local/bin/python
# -*- coding: utf-8 -*-
from numpy import *
import cv2

def main():
    # 画像を読み込む
    filename = 'fisherman.jpg'
    im = cv2.imread(filename)
    h,w = im.shape[:2]

    # 塗りつぶしの例
    diff = (6,6,6)
    mask = zeros((h+2,w+2),uint8)
    cv2.floodFill(im,           # 元の画像
                  mask,         # ????
                  (10,10),      # 塗りつぶしの開始座標
                  (0,0,255),    # 塗りつぶしの色(RGBの逆?(B->G->R))
                  diff,         # 塗りつぶしの上限
                  diff)         # 塗りつぶしの下限

    # 結果を保存する
    cv2.imwrite('result_floodFill.jpg',im)

if __name__ == '__main__':
    main()

↑こう書くと、↓こう変換

f:id:end0tknr:20101030192652j:plain:w250 f:id:end0tknr:20171029081900j:plain:w250

OpenCV3では、OpenCV2と異なり、画像サイズの縦横が反転する気がしますが、 気のせいかな?

積分画像作成

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

def main():
    # 画像を読み込む
    im = cv2.imread('empire.jpg')
    gray = cv2.cvtColor(im,cv2.COLOR_BGR2GRAY)

    # 積分画像を計算する
    intim = cv2.integral(gray)

    # 255内に収まるように正規化し、保存
    intim = (255.0*intim) / intim.max()
    cv2.imwrite('result_integral.jpg',intim)

if __name__ == '__main__':
    main()

画像積分の仕様も、利用目的も理解していません...が ↑こう書くと、↓こう変換

f:id:end0tknr:20081110180136j:plain:w250  f:id:end0tknr:20171029082430j:plain:w250

OpenCV3 に付属のサンプルソース

opencv/samples/python 以下に大量にあります

Hough変換による直線検出

#!/usr/local/bin/python

# Python 2/3 compatibility
from __future__ import print_function

import cv2
import numpy as np
import sys
import math


def main():
    fn = sys.argv[1]

    src = cv2.imread(fn)
    dst = cv2.Canny(src, 50, 200)
    cdst = cv2.cvtColor(dst, cv2.COLOR_GRAY2BGR)

    if True: # HoughLinesP
        lines = cv2.HoughLinesP(dst, 1, math.pi/180.0, 40, np.array([]), 50, 10)
        a,b,c = lines.shape
        for i in range(a):
            cv2.line(cdst,
                     (lines[i][0][0], lines[i][0][1]),
                     (lines[i][0][2], lines[i][0][3]),
                     (0, 0, 255), 3, cv2.LINE_AA)

    else:    # HoughLines
        lines = cv2.HoughLines(dst, 1, math.pi/180.0, 50, np.array([]), 0, 0)
        if lines is not None:
            a,b,c = lines.shape
            for i in range(a):
                rho = lines[i][0][0]
                theta = lines[i][0][1]
                a = math.cos(theta)
                b = math.sin(theta)
                x0, y0 = a*rho, b*rho
                pt1 = ( int(x0+1000*(-b)), int(y0+1000*(a)) )
                pt2 = ( int(x0-1000*(-b)), int(y0-1000*(a)) )
                cv2.line(cdst, pt1, pt2, (0, 0, 255), 3, cv2.LINE_AA)

            cv2.imshow("detected lines", cdst)

    cv2.imwrite('houghlines.png',cdst)


if __name__ == '__main__':
    main()

cv2.HoughLinesP()と、cv2.HoughLines()の違いを理解していませんが、 ↑こう書くと、↓こう変換

f:id:end0tknr:20081110180136j:plain:w250  f:id:end0tknr:20171029082831p:plain:w250

install opencv 3.3.1 to centos 7.4 from src

概要

CentOS 7 に OpenCV git 最新版 (3.0.0-rc1) をインストールする - Qiita

上記urlを参考にさせて頂きました。

ただし、

$ make -j
   :
g++: internal compiler error: Killed (program cc1plus)

のようなメモリ不足?によるエラーとなった為、「make」としました。

また、

$ make
  :
In file included from /home/endo/tmp/opencv/modules/videoio/src/cap_ffmpeg.cpp:47:0:
/home/endo/tmp/opencv/modules/videoio/src/cap_ffmpeg_impl.hpp: In function ‘AVStream* icv_add_video_stream_FFMPEG(AVFormatContext*, AVCodecID, int, int, int, double, int)’:
/home/endo/tmp/opencv/modules/videoio/src/cap_ffmpeg_impl.hpp:1573:21: error: ‘CODEC_FLAG_GLOBAL_HEADER’ was not declared in this scope
         c->flags |= CODEC_FLAG_GLOBAL_HEADER;
  :
/home/endo/tmp/opencv/modules/videoio/src/cap_ffmpeg_impl.hpp: In static member function ‘static AVStream* OutputMediaStream_FFMPEG::addVideoStream(AVFormatContext*, AVCodecID, int, int, int, double, AVPixelFormat)’:
/home/endo/tmp/opencv/modules/videoio/src/cap_ffmpeg_impl.hpp:2379:25: error: ‘CODEC_FLAG_GLOBAL_HEADER’ was not declared in this scope
             c->flags |= CODEC_FLAG_GLOBAL_HEADER;
                         ^
make[2]: *** [modules/videoio/CMakeFiles/opencv_videoio.dir/src/cap_ffmpeg.cpp.o] Error 1
make[1]: *** [modules/videoio/CMakeFiles/opencv_videoio.dir/all] Error 2
make: *** [all] Error 2

ようなエラーもありました。ググると、opencv3.3.1 が ffmpeg3.4に対応していないようでしたので ver.3.2.9を使用しています。

詳細

# yum install autoconf automake cmake freetype-devel gcc gcc-c++ \
      git libtool make mercurial nasm pkgconfig zlib-devel
$ git clone git://github.com/yasm/yasm.git
$ cd yasm
$ autoreconf -fiv
$ ./configure
$ make
# make install
$ wget http://www.nasm.us/pub/nasm/releasebuilds/2.13.01/nasm-2.13.01.tar.gz
$ tar -xvf nasm-2.13.01.tar.gz
$ cd nasm-2.13.01
$ ./configure
$ make
# make install
$ git clone git://git.videolan.org/x264
$ cd x264
$ ./configure --enable-static --enable-pic
$ make
# make install
$ git clone git://git.code.sf.net/p/opencore-amr/fdk-aac
$ cd fdk-aac
$ autoreconf -fiv
$ ./configure --disable-shared --with-pic
$ make
# make install
$ wget http://downloads.sourceforge.net/project/lame/lame/3.99/lame-3.99.5.tar.gz
$ tar -xvf lame-3.99.5.tar.gz
$ cd lame-3.99.5
X $ ./configure --disable-shared --enable-nasm --with-pic
$ ./configure --enable-nasm --with-pic
$ make
# make install
$ git clone git://git.opus-codec.org/opus.git
$ cd opus
$ autoreconf -fiv
X $ ./configure --disable-shared --with-pic
$ ./configure --with-pic
$ make
# make install
$ wget http://downloads.xiph.org/releases/ogg/libogg-1.3.2.tar.gz
$ tar xzvf libogg-1.3.2.tar.gz
$ cd libogg-1.3.2
$ ./configure --disable-shared --with-pic
$ make
# make install
$ wget http://downloads.xiph.org/releases/vorbis/libvorbis-1.3.4.tar.gz
$ tar xzvf libvorbis-1.3.4.tar.gz
$ cd libvorbis-1.3.4
X $ ./configure --with-ogg --disable-shared --with-pic
$ ./configure --with-ogg --with-pic
$ make
# make install
$ git clone https://chromium.googlesource.com/webm/libvpx.git
$ cd libvpx
$ ./configure --disable-examples --enable-pic
$ make
# make install
# echo "/usr/local/lib" | sudo tee /etc/ld.so.conf.d/usr-local-lib.conf
# ldconfig
$ wget http://ffmpeg.org/releases/ffmpeg-3.2.9.tar.gz 
$ tar -xvf ffmpeg-3.2.9.tar.gz
$ cd ffmpeg-3.2.9
$ PKG_CONFIG_PATH=/usr/local/lib/pkgconfig ./configure \
  --enable-gpl --enable-nonfree --enable-libfdk_aac \
  --enable-libfreetype --enable-libmp3lame --enable-libopus \
  --enable-libvorbis --enable-libvpx --enable-libx264 \
  --enable-pic --disable-static --enable-shared
$ make
# make install
$ git clone https://github.com/Itseez/opencv.git
$ cd opencv
$ mkdir build
$ cd build
$ PKG_CONFIG_PATH=/usr/local/lib/pkgconfig cmake \
  -DPNG_INCLUDE_DIR=/usr/local/include \
  -DPNG_LIBRARY_RELEASE=/usr/local/lib/libpng16.so \
  ..
$ make
# make install

cmakeのオプションに -DPNG_INCLUDE_DIR と -DPNG_LIBRARY_RELEASE を付けています。 これは、私の環境にlibpng yum install版と、src install版があり、 次のようなエラーとなった為です

$ /usr/local/bin/opencv_createsamples \
  -img window_900.png  -vec window_900.vec \
  -num 1000 -bgcolor 255 -w 50 -h 100
    :
Create training samples from single image applying distortions...
libpng warning: Application built with libpng-1.6.32 but running with 1.5.13

HOG特徴量 + ナイーブベイズ による ジェスチャー(手話)画像認識

先日は、k近傍による 画像認識を行いましたが、今回は、ナイーブベイズを使用。

HOG特徴量 + KNN(k近傍) による ジェスチャー(手話)画像認識 - end0tknr's kipple - 新web写経開発

以下は、更に以前の主成分分析に関するエントリ

固有値、固有ベクトルからの主成分分析 オレオレ入門 (PCA: Principal Component Analysis) - end0tknr's kipple - 新web写経開発

PIL / Pillow , numpy による画像処理入門 その3 - 平板化(1次元化)と主成分分析(principal component analysis; PCA) - end0tknr's kipple - 新web写経開発

ただし、単純にベイズ分類器を使用せず、主成分分析(PCA)により次元数を 50に削減し、実行。

まっ、いつもの

実践コンピュータビジョン サンプルプログラム

の写経です

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

from PIL import Image
import numpy
import os
import bayes
import math

def main():
    # 訓練データとtestデータのload
    features, labels = read_gesture_features_labels('train/')
    test_features, test_labels = read_gesture_features_labels('test/')

    classnames = numpy.unique(labels)
    print "CATEGORIES:",
    print classnames

    # 主成分分析 - 返り値:写像行列(次元の重要度順), 分散, 平均
    V, S, m = my_pca(features)
    # 重要な次元を50コ残す
    V = V[:50]

    # numpy.dot():内積
    features = numpy.array([numpy.dot(V, f - m) for f in features])
    test_features = numpy.array([numpy.dot(V, f - m) for f in test_features])

    # ベイズ分類器
    bc = BayesClassifier()
    blist = [features[numpy.where(labels == c)[0]] for c in classnames]
    bc.train(blist, classnames)
    res = bc.classify(test_features)[0]
    acc = sum(1.0 * (res == test_labels)) / len(test_labels)
    print 'Accuracy:', acc
    print_confusion(res, test_labels, classnames)


def read_gesture_features_labels(path):
    # *.dsift のfile名list
    featlist = [
        os.path.join(path, f) for f in os.listdir(path) if f.endswith('.dsift')
    ]
    # 特徴量をload
    features = []
    for featfile in featlist:
        l, d = read_features_from_file(featfile)
        features.append(d.flatten())
    features = numpy.array(features)
    # ラベルを生成する
    labels = [featfile.split('/')[-1][0] for featfile in featlist]
    return features, numpy.array(labels)


def print_confusion(res, test_labels, classnames):
    n = len(classnames)
    # 混同行列
    class_ind = dict([(classnames[i], i) for i in range(n)])
    confuse = numpy.zeros((n, n))
    for i in range(len(test_labels)):
        confuse[class_ind[res[i]], class_ind[test_labels[i]]] += 1
    print 'Confusion matrix for'
    print classnames
    print confuse


# 特徴量を読込み行列形式で返す
def read_features_from_file(filename):
    f = numpy.loadtxt(filename)
    return f[:, :4], f[:, 4:]  # 特徴点の配置と記述子


# 主成分分析
#   入力:X, 訓練データを平板化した配列を行として格納した行列
#   出力:写像行列(次元の重要度順), 分散, 平均
def my_pca(X):
    # 次元数を取得
    num_data, dim = X.shape
    # データをセンタリング
    mean_X = X.mean(axis=0)
    X = X - mean_X
    if dim > num_data:
        # PCA - 高次元のときはコンパクトな裏技を用いる
        M = numpy.dot(X, X.T)  # 共分散行列
        e, EV = numpy.linalg.eigh(M)  # 固有値と固有ベクトル
        tmp = numpy.dot(X.T, EV).T  # ここがコンパクトな裏技
        V = tmp[::-1]  # 末尾の固有ベクトルほど重要なので、反転する
        S = numpy.sqrt(e)[::-1]  # 固有値の並びも反転する
        for i in range(V.shape[1]):
            V[:, i] /= S
    else:
        # PCA - 低次元なら特異値分解を用いる
        U, S, V = numpy.linalg.svd(X)
        V = V[:num_data]  # 最初のnum_dataの分だけが有用
    # 写像行列と、分散、平均を返す
    return V, S, mean_X


class BayesClassifier(object):
    def __init__(self):
        self.labels = []  # クラスのラベル
        self.mean = []    # クラスの平均
        self.var = []     # クラスの分散
        self.n = 0        # クラスの数

        
    # data (n*dimの配列のリスト)で学習。
    # labelsはオプションで、デフォルトは 0...n-1 
    def train(self, data, labels=None):
        if labels == None:
            labels = range(len(data))
        self.labels = labels
        self.n = len(labels)
        for c in data:
            self.mean.append(numpy.mean(c, axis=0))
            self.var.append(numpy.var(c, axis=0))  # var():分散

    # 各クラスの確率を計算し確率の高いラベルを返すことでpointsを分類
    def classify(self, points):
        # 各クラスの確率を計算する
        est_prob = numpy.array(
            [gauss(m, v, points) for m, v in zip(self.mean, self.var)])
        # 最も確率の高いindex noを求め、クラスのラベルを返す
        ndx = est_prob.argmax(axis=0)
        est_labels = numpy.array([self.labels[n] for n in ndx])
        return est_labels, est_prob


# 独立した平均m分散vの点を、xの行として持つ d次元ガウス分布を評価
def gauss(m, v, x):
    if len(x.shape) == 1:
        n, d = 1, x.shape[0]
    else:
        n, d = x.shape
        
    # 共分散行列を求め、xから平均を引く
    S = numpy.diag(1 / v)  # diag():行列の対角成分
    x = x - m
    # 確率の積,  exp():ネイピア数(e)の?乗, dot():内積
    y = numpy.exp(-0.5 * numpy.diag(numpy.dot(x, numpy.dot(S, x.T))))
    # 正規化して返す, prod():配列の積
    return y * (2 * numpy.pi)**(-d / 2.0) / (numpy.sqrt(numpy.prod(v)) + 1e-6)


if __name__ == '__main__':
    main()

↑こう書くと↓こう表示されます

$ python 0_8.2.1.bayes_g.py 
CATEGORIES: ['A' 'B' 'C' 'F' 'P' 'V']
0_8.2.1.bayes_g.py:87: RuntimeWarning: invalid value encountered in sqrt
  S = numpy.sqrt(e)[::-1]  # 固有値の並びも反転する
0_8.2.1.bayes_g.py:109: FutureWarning: comparison to `None` will result in an elementwise object comparison in the future.
  if labels == None:
Accuracy: 0.744680851064
Confusion matrix for
['A' 'B' 'C' 'F' 'P' 'V']
[[ 23.   1.   3.   2.   0.   0.]
 [  0.  22.   0.   3.   2.   1.]
 [  0.   1.  28.   1.   1.   1.]
 [  1.   0.   0.  22.   0.   0.]
 [  0.   5.   2.   0.  23.   5.]
 [  5.   0.   3.   1.  10.  22.]]

srcではガウス確率分布(正規分布)モデルを使ってベイズ分類器を実装していますが、 未だ私の理解は今ひとつ

画像認識 - ナイーブベイズ(単純ベイズ)分類器 by python

で、ベイズ分類器の実装

参考url

分類器 (ナイーブベイズ)

機械学習ナイーブベイズ分類器のアルゴリズムを理解したのでメモ。そしてPythonで書いてみた。 - Qiita

第3回 ベイジアンフィルタを実装してみよう:機械学習 はじめよう|gihyo.jp … 技術評論社

ナイーブベイズ分類器とは?

ある文書Dがあるとき、その文書がカテゴリCに属する確率は、 ベイズの定理により、次のように表せます。

{ \Large
P(C|D) = \frac{ P(D|C) P(C) }{ P(D) }
}

上記において

P(D) は、カテゴリに依存せず、無視しても P(C|D)の大小は判定できる為、無視。

P(C) は、カテゴリが得られる確率で

{ \Large
P(C) = \frac{ カテゴリCと判定された文書数 }{ 全文書数 }
}

P(D|C) は、カテゴリが与えられた場合、文書Dが得られる確率で

{ \Large
P(D|C) = \frac{ 文書Dの数 }{ カテゴリCに属する文書数 }
}

となるが、実際には以下の近似値を使用する。

{ \Large
P(D|C) ≒ P(W_1|C) \cdot P(W_2|C) \cdots  P(W_n|C)
}
{ \Large
P(W_i|C) = \frac{ 単語 W_i の数}{ カテゴリCにおける語彙数 }
}

例題:ある文書のカテゴリ推定

以下の訓練データがある場合、文書「ボールスポーツ」が属するカテゴリを推定

NO カテゴリ 出現単語
1 サッカー ボール、スポーツ、ワールドカップ、ボール
2 野球 ボール、スポーツ、グローブ、バット
3 テニス ボール、ラケット、コート
4 サッカー ボール、スポーツ

P(C) (各カテゴリが得られる確率)の計算

- サッカー 野球 テニス
P(C) 2/4 1/4 1/4

P(D|C) (あるカテゴリにおいて文書Dが得られる確率)の計算

まず、語彙数をカウント (このとき重複もそのままカウント)

- サッカー 野球 テニス
語彙数 6 4 3

次にカテゴリ内の単語出現率を算出。 今回は文書「ボールスポーツ」の判定ですので「ボール」と「スポーツ」を算出。

P(Wi|C) サッカー 野球 テニス
ボール 3/6 1/4 1/3
スポーツ 2/6 1/4 0/3

よって、

- サッカー 野球 テニス
P(D|C) 6/36 1/8 0/9

P(C|D)≒ P(C) P(D|C) であるから

- サッカー 野球 テニス
P(C|D) (2/4)x(6/36)= 1/12 (1/4)x(1/8)= 1/64 ((1/4)x(0/9)=0/36

以上より、文書「ボールスポーツ」は、サッカーカテゴリに属すると考えられる

ナイーブベイズ分類器 by python による画像分類

実践コンピュータビジョン サンプルプログラム いつもの?↑こちらの写経。今回は8章。

尚、訓練データとテストデータは以前のエントリにある oints_normal.pkl と points_normal_test.pkl を使用しています pythonでのk近傍法( k-nearest neighbor algorithm, k-NN ) - end0tknr's kipple - 新web写経開発

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb
import numpy
import pickle


def main():
    # 訓練データのload
    with open('points_normal.pkl', 'r') as f:
#    with open('points_ring.pkl', 'r') as f:
        class_1 = pickle.load(f)
        class_2 = pickle.load(f)
        labels = pickle.load(f)

    # ベイズ分類器による学習
    bc = BayesClassifier()
    bc.train([class_1, class_2], [1, -1])
    
    # テストデータのload
    with open('points_normal_test.pkl', 'r') as f:
#    with open('points_ring_test.pkl', 'r') as f:
        class_1 = pickle.load(f)
        class_2 = pickle.load(f)
        labels = pickle.load(f)
        

    # いくつかの点についてテスト
    print bc.classify(class_1[:10])[0]
    
    # 点と判別境界を表示
    def classify0(x, y, bc=bc):
        points = numpy.vstack((x, y))
        return bc.classify(points.T)[0]
    plot_2D_boundary([-6, 6, -6, 6], [class_1, class_2], classify0, [1,-1])

    plb.savefig( '8_2_bayes_n.png' )


# 独立した平均m分散vの点を、xの行として持つ d次元ガウス分布を評価
def gauss(m, v, x):
    if len(x.shape) == 1:
        n, d = 1, x.shape[0]
    else:
        n, d = x.shape
        
    # 共分散行列を求め、xから平均を引く
    S = numpy.diag(1 / v)  # diag():行列の対角成分
    x = x - m
    # 確率の積,  exp():ネイピア数(e)の?乗, dot():内積
    y = numpy.exp(-0.5 * numpy.diag(numpy.dot(x, numpy.dot(S, x.T))))
    # 正規化して返す, prod():配列の積
    return y * (2 * numpy.pi)**(-d / 2.0) / (numpy.sqrt(numpy.prod(v)) + 1e-6)

class BayesClassifier(object):
    def __init__(self):
        self.labels = []  # クラスのラベル
        self.mean = []    # クラスの平均
        self.var = []     # クラスの分散
        self.n = 0        # クラスの数

        
    # data (n*dimの配列のリスト)で学習。
    # labelsはオプションで、デフォルトは 0...n-1 
    def train(self, data, labels=None):
        if labels == None:
            labels = range(len(data))
        self.labels = labels
        self.n = len(labels)
        for c in data:
            self.mean.append(numpy.mean(c, axis=0))
            self.var.append(numpy.var(c, axis=0))  # var():分散

    # 各クラスの確率を計算し確率の高いラベルを返すことでpointsを分類
    def classify(self, points):
        # 各クラスの確率を計算する
        est_prob = numpy.array(
            [gauss(m, v, points) for m, v in zip(self.mean, self.var)])
        # 最も確率の高いindex noを求め、クラスのラベルを返す
        ndx = est_prob.argmax(axis=0)
        est_labels = numpy.array([self.labels[n] for n in ndx])
        return est_labels, est_prob


# plot_range:(xmin,xmax,ymin,ymax)、points:クラスの点のリスト、
# decisionfcn:評価関数、
# labels:各クラスについてdecisionfcnが返すラベルのリスト、
# values:表示する判別の輪郭のリスト """
def plot_2D_boundary(plot_range, points, decisionfcn, labels, values=[0]):
    clist = ['b', 'r', 'g', 'k', 'm', 'y']  # クラスの描画色
    # グリッドを評価して判別関数の輪郭を描画する
    x = numpy.arange(plot_range[0], plot_range[1], .1)
    y = numpy.arange(plot_range[2], plot_range[3], .1)
    xx, yy = numpy.meshgrid(x, y)
    xxx, yyy = xx.flatten(), yy.flatten()  # グリッドのx,y座標のリスト
    zz = numpy.array(decisionfcn(xxx, yyy))
    zz = zz.reshape(xx.shape)
    # valuesにそって輪郭を描画する
    plb.contour(xx, yy, zz, values)
    # クラスごとに正しい点には*、間違った点には'o'を描画する
    for i in range(len(points)):
        d = decisionfcn(points[i][:, 0], points[i][:, 1])
        correct_ndx = labels[i] == d
        incorrect_ndx = labels[i] != d
        plb.plot(
            points[i][correct_ndx, 0],
            points[i][correct_ndx, 1],
            '*',
            color=clist[i])
        plb.plot(
            points[i][incorrect_ndx, 0],
            points[i][incorrect_ndx, 1],
            'o',
            color=clist[i])
    plb.axis('equal')

if __name__ == '__main__':
    main()

↑こう書くと↓こう表示されますが、ベイズ分類器の画像への適用はあまり理解できていません。ホント...コードを写経しただけ f:id:end0tknr:20171022145923p:plain

ベイズの定理の証明 - 条件付き確立

単純ベイズ(ナイーブベイズ)分類器やモンティ・ホール問題に入る前の振り返り

ベイズの定理とは?

全ての事象が独立で、互いに相関がない場合、以下が成立

{ \Large
P(B|A) = \frac { P(A∩B) }{P(A)} = \frac { P(A|B) P(B) }{P(A)} … 式1
}

ここで

  • P(B|A) : 事象Aの後に、事象Bが発生する確率
  • P(A∩B) : 事象Aと事象Bが同時に発生する確率

「全ての事象が独立で、互いに相関がない」の仮定が、「単純(ナイーブ)」のゆえん

f:id:end0tknr:20171022074451p:plain

証明

と言っても、

{ \Large {
P(A) =   \frac{A}{U} 、P(B) =  \frac{B}{U} \\
P(A|B) = \frac{A∩B}{B} 、P(B|A) = \frac{A∩B}{A} }
}

上記を式1へ代入するだけ

HOG特徴量 + KNN(k近傍) による ジェスチャー(手話)画像認識

先程のエントリの続きとして、ジェスチャー画像のHOG特徴量を算出し、 その結果を、k近傍することで、手話の画像認識を行います。

画像処理(画像認識) - 「密なSIFT(HOG : Histogram of Oriented Gradients )」算出 - end0tknr's kipple - 新web写経開発

pythonでのk近傍法( k-nearest neighbor algorithm, k-NN ) - end0tknr's kipple - 新web写経開発

前準備 - 訓練データと、確認用データ

実践コンピュータビジョン サンプルプログラム

に記載されている通り、

http://www.idiap.ch/resource/gestures/data/shp_marcel_test.tar.gz

の uniform/ 以下にある {A,B,C,Five,Point,V} の 画像を train と test dirにcopy 。

次にのように一気に行うこともできます

$ cd ~/tmp/VISION/chap8
$ mkdir train
$ mkdir test
$ wget http://www.idiap.ch/resource/gestures/data/shp_marcel_test.tar.gz
$ tar -xvf shp_marcel_test.tar.gz
$ mv Marcel-Test/*/uniform/*uniform[0-2]?.ppm train
$ mv Marcel-Test/*/uniform/*.ppm test

STEP.1 各画像の HOG特徴量を算出

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
from PIL import Image
import os
import dsift
import numpy

## http://www.vlfeat.org/download/vlfeat-0.9.20.tar.gz
SIFT_CMD = '/home/endo/local/vlfeat/sift'

def main():
    d = 'test'  # test用data
    imlist = [os.path.join(d,f) for f in os.listdir(d) if f.endswith('.ppm')]
    d = 'train' # 訓練用data
    imlist += [os.path.join(d,f) for f in os.listdir(d) if f.endswith('.ppm')]

    # 50x50sizeで、「密なSIFT記述子(HOG)」算出
    for filename in imlist:
        featfile = filename[:-3]+'dsift'
        process_image_dsift(filename,featfile,10,5,resize=(50,50))


# 画像から「密なSIFT記述子(HOG)」を抽出し、fileに保存.
#   args size=特徴量size, steps=grid間隔, resize=画像resize用のタプル
#   force_orientation= True →最大の局所勾配の方向で記述子が正規化
#                      False→全て上向き
def process_image_dsift(imagename,
                        resultname,
                        size=20,
                        steps=10,
                        force_orientation=False,
                        resize=None):
    im = Image.open(imagename).convert('L')  # file open後、グレースケール化
    if resize != None:
        im = im.resize(resize)
    m, n = im.size

    if imagename[-3:] != 'pgm':
        im.save('tmp.pgm')  # pgmファイルを作成
        imagename = 'tmp.pgm'

    # frameを作成し一時fileに保存
    scale = size / 3.0
    x, y = numpy.meshgrid(range(steps, m, steps), range(steps, n, steps))
    xx, yy = x.flatten(), y.flatten()
    frame = numpy.array([xx, yy,
                         scale * numpy.ones(xx.shape[0]),
                         numpy.zeros(xx.shape[0])])

    numpy.savetxt('tmp.frame', frame.T, fmt='%03.3f')

    if force_orientation:
        cmmd = str(SIFT_CMD +" "+ imagename + " --output=" + resultname +
                   " --read-frames=tmp.frame --orientations")
    else:
        cmmd = str(SIFT_CMD +" "+ imagename + " --output=" + resultname +
                   " --read-frames=tmp.frame")
    os.system(cmmd)

    print 'processed', imagename, 'to', resultname


if __name__ == '__main__':
    main()

↑こう書くと、test , train dir 以下にある画像のHOGを算出し、 結果を *.dsift ファイルに保存します

STEP.2 k-NNによる画像認識

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
import numpy
import os


def main():
    features, labels = read_gesture_features_labels('train/')
    test_features, test_labels = read_gesture_features_labels('test/')
    classnames = numpy.unique(labels)

    # k近傍法を試す
    k = 1
    knn_classifier = KnnClassifier(labels, features)
    res = numpy.array([
        knn_classifier.classify(test_features[i], k)
        for i in range(len(test_labels))
    ])
    # 精度
    acc = sum(1.0 * (res == test_labels)) / len(test_labels)
    print 'Accuracy:', acc
    print_confusion(res, test_labels, classnames)


# HOGが登録された *.dsift の filename list 作成
def read_gesture_features_labels(path):
    featlist = [
        os.path.join(path, f) for f in os.listdir(path) if f.endswith('.dsift')
    ]

    # 特徴量を読込み
    features = []
    for featfile in featlist:
        l, d = read_features_from_file(featfile)
        features.append(d.flatten())
    features = numpy.array(features)
    
    # file名からラベルを作成
    labels = [featfile.split('/')[-1][0] for featfile in featlist]
#    print labels
    
    return features, numpy.array(labels)


# 特徴量を読込み行列形式で返す
def read_features_from_file(filename):
    f = numpy.loadtxt(filename)
    return f[:, :4], f[:, 4:]  # 特徴点の配置と記述子


def print_confusion(res, test_labels, classnames):
    n = len(classnames)
    # 混同行列
    class_ind = dict([(classnames[i], i) for i in range(n)])
    confuse = numpy.zeros((n, n))
    for i in range(len(test_labels)):
        confuse[class_ind[res[i]], class_ind[test_labels[i]]] += 1
    print 'Confusion matrix for'
    print classnames
    print confuse


class KnnClassifier(object):
  def __init__(self,labels,samples):
    """ 教師データを使って分類器を初期化する """
    self.labels = labels
    self.samples = samples

  def classify(self,point,k=3):
    """ pointを教師データ中のk個の最近傍を使って分類し、
        ラベルを返す """

    # 全教師データへの距離を計算する
    dist = numpy.array([L2dist(point,s) for s in self.samples])
    # ソートする
    ndx = dist.argsort()
    # k個の最近傍を保持するのに辞書を用いる
    votes = {}
    for i in range(k):
      label = self.labels[ndx[i]]
      votes.setdefault(label,0)
      votes[label] += 1

    return max(votes, key=lambda x: votes.get(x))


def L2dist(p1,p2):
    return numpy.sqrt( sum( (p1-p2)**2) )


if __name__ == '__main__':
    main()

↑こう書くと↓こう表示されます

$ python 0_8.1.3_2.knn_g.py 
Accuracy: 0.617021276596
Confusion matrix for
  [   'A'  'B'  'C'  'F'  'P'  'V']
[['A' 23.   2.   2.   1.   0.   2.]
 ['B'  2.  20.   2.   5.   4.   1.]
 ['C'  0.   1.  23.   1.   0.   0.]
 ['F'  0.   1.   0.  20.   0.   0.]
 ['P'  0.   5.   3.   1.  12.   8.]
 ['V'  4.   0.   6.   1.  20.  18.]]

特に右から2列目の結果(12と30)から、P が Vと混同されやすいことが分かります。

画像処理(画像認識) - 「密なSIFT(HOG : Histogram of Oriented Gradients )」算出

oreillyの「実践コンピュータビジョン」写経の続き

実践コンピュータビジョン サンプルプログラム

SIFT特徴量検出の続きでもあります

画像処理 - vlfeatによるSIFT特徴点検出 - end0tknr's kipple - 新web写経開発

SIFT特徴点検出結果の対応点探索 (python) - end0tknr's kipple - 新web写経開発

密なSIFT(HOG : Histogram of Oriented Gradients ) とは

SIFT特徴量は、特徴点の周辺の勾配に対し、勾配ヒストグラムを算出しましたが、 HOGは、一定間隔で勾配ヒストグラムを算出します。

pythonHOGを算出し、その座標を表示

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
from PIL import Image
import numpy
import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb
import os
## http://www.vlfeat.org/download/vlfeat-0.9.20.tar.gz
SIFT_CMD = '/home/endo/local/vlfeat/sift'


def main():
    # 画像から「密なSIFT記述子(HOG)」を抽出し、fileに保存.
    process_image_dsift('empire.jpg', 'empire.sift', 90, 40, True)

    l,d = read_features_from_file('empire.sift')
    im = numpy.array(Image.open('empire.jpg'))
    plot_features(im,l,True)
    plb.savefig( '8_1_2.png' )


# 画像から「密なSIFT記述子(HOG)」を抽出し、fileに保存.
#   args size=特徴量size, steps=grid間隔, resize=画像resize用のタプル
#   force_orientation= True →最大の局所勾配の方向で記述子が正規化
#                      False→全て上向き
def process_image_dsift(imagename,
                        resultname,
                        size=20,
                        steps=10,
                        force_orientation=False,
                        resize=None):
    im = Image.open(imagename).convert('L')  # file open後、グレースケール化
    if resize != None:
        im = im.resize(resize)
    m, n = im.size

    if imagename[-3:] != 'pgm':
        im.save('tmp.pgm')  # pgmファイルを作成
        imagename = 'tmp.pgm'

    # frameを作成し一時fileに保存
    scale = size / 3.0
    x, y = numpy.meshgrid(range(steps, m, steps), range(steps, n, steps))
    xx, yy = x.flatten(), y.flatten()
    frame = numpy.array([xx, yy,
                         scale * numpy.ones(xx.shape[0]),
                         numpy.zeros(xx.shape[0])])

    numpy.savetxt('tmp.frame', frame.T, fmt='%03.3f')

    if force_orientation:
        cmmd = str(SIFT_CMD +" "+ imagename + " --output=" + resultname +
                   " --read-frames=tmp.frame --orientations")
    else:
        cmmd = str(SIFT_CMD +" "+ imagename + " --output=" + resultname +
                   " --read-frames=tmp.frame")
    os.system(cmmd)

    print 'processed', imagename, 'to', resultname


# 特徴量を読込み行列形式で返す
def read_features_from_file(filename):
    f = numpy.loadtxt(filename)
    return f[:, :4], f[:, 4:]  # 特徴点の配置と記述子

# 画像を特徴量とともに描画。
# 入力:im(配列形式の画像),locs(各特徴量の座標とスケール、方向)
def plot_features(im, locs, circle=False):
    def draw_circle(c, r):
        t = numpy.arange(0, 1.01, .01) * 2 * numpy.pi
        x = r * numpy.cos(t) + c[0]
        y = r * numpy.sin(t) + c[1]
        plb.plot(x, y, 'b', linewidth=2)

    plb.imshow(im)

    if circle:
        for p in locs:
            draw_circle(p[:2],p[2])
    else:
        plb.plot(locs[:,0],locs[:,1],'ob')
    plb.axis('off')


if __name__ == '__main__':
    main()

↑こう書くと、↓こう表示されます。(単に座標を表示しているので、全く面白くありませんが) f:id:end0tknr:20171021194701p:plain

emacs 24.5 for win に python 開発環境( コード整形 )を整える

perl tidy と同様のコード整形ツールが必要になった為

python-mode

install時期やきっかけは、すっかり忘れていますが、 .emacsには次のように記載されていました。

(add-to-list 'load-path "c:/emacs-24.5-IME-patched/python-mode.el-6.2.0")
(setq py-install-directory "c:/emacs-24.5-IME-patched/python-mode.el-6.2.0")
(require 'python-mode)

install python 3.6 for win

普段、python 2.7を使用している為、一旦、python2.7をinstallしましたが、 emacspython srcを開いた際、errorが表示されるようになった為、 python3.6をinstall. (原因は調べていません)

autopep8 や yapf 等のコード整形ツールは不使用

pythonの code formatterには、autopep8 や yapf があるようですが 私のemacs環境では、上手く動作しなかった為、使用しませんでした

install elpy (python ide for emacs )

https://elpy.readthedocs.io/en/latest/introduction.html

installは上記urlに記載の通り

(require 'package)
(add-to-list 'package-archives
             '("elpy" . "http://jorgenschaefer.github.io/packages/"))
(elpy-enable)

↑このように .emacs に追記し、「 M-x package-install elpy 」を実行

elpy によるcode format

「 M-x elpy-mode 」で起動後、範囲指定し、 「 M-x elpy-autopep8-fix-code」 or 「 M-x elpy-format-code 」 両者の違いは理解していませんが、上手く動作しているように見えます。

Highlight-Indentation-for-Emacs によるインデントのハイライト

https://github.com/antonj/Highlight-Indentation-for-Emacs

から highlight-indentation.el をsite-listへダウンロードし、 .emacsに「 (require 'highlight-indentation) 」を追記すると、 インデントのグループ?が表示されます。

(その他) 指定した範囲をまとめて字下げ

「 C-c > 」or「 C-c < 」

emacs for win の M-x package install で http://~ Not found エラー

「M-x package install elpy」を実行したら、「http://~ Not found 」というエラー

.emacs には

;; http://emacs-jp.github.io/packages/package-management/package-el.html
(require 'package)
(add-to-list 'package-archives
         '("melpa" . "http://melpa.milkbox.net/packages/"))
(add-to-list 'package-archives
         '("marmalade" . "http://marmalade-repo.org/packages/"))
(add-to-list 'package-archives
             '("elpy" . "http://jorgenschaefer.github.io/packages/"))
(package-initialize)

と設定されているのに、と思ったら、emacsのcacheの問題らしい。

「 M-x package-refresh-contents 」と実行することで解消。

SIFT特徴点検出結果の対応点探索 (python)

正規化相互相関を用いた Harrisコーナー検出結果の対応点探索 (python) - end0tknr's kipple - 新web写経開発

画像処理 - vlfeatによるSIFT特徴点検出 - end0tknr's kipple - 新web写経開発

上記のエントリの続きです。

最も一致する特徴点への距離の比を双方向で比較する がポイントです

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

from PIL import Image
import numpy
import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb
import os

#imname1 = 'crans_1_small.jpg'
#imname2 = 'crans_2_small.jpg'
imname1 = 'sf_view1.jpg'
imname2 = 'sf_view2.jpg'
#imname1 = 'climbing_1_small.jpg'
#imname2 = 'climbing_2_small.jpg'

SIFT_CMD = '/home/endo/local/vlfeat/sift'
PROCESS_IMAGE_PARAMS = "--edge-thresh 10 --peak-thresh 5"

def main():
    im1 = numpy.array(Image.open(imname1).convert('L'))
    im2 = numpy.array(Image.open(imname2).convert('L'))
    # 特徴量検出
    process_image(imname1,imname1 + '.sift')
    l1,d1 = read_features_from_file(imname1 + '.sift')
    process_image(imname2,imname2 + '.sift')
    l2,d2 = read_features_from_file(imname2 + '.sift')

    # 特徴量のマッチング
    matches = match_twosided(d1,d2)

    # マッチング結果の表示
    plb.figure()
    plb.gray()
    plot_matches(im1,im2,l1,l2,matches)
    plb.savefig( '2_3_4.png' )


# 画像の特徴点を検出しfileに保存
def process_image(imagename,resultname,params=PROCESS_IMAGE_PARAMS):
    if imagename[-3:] != 'pgm':
        # pgmファイルを作成する
        im = Image.open(imagename).convert('L')
        imagename = 'tmp.pgm'
        im.save(imagename)

    cmmd = str(SIFT_CMD +" "+imagename+" --output="+resultname+ " "+params)
    os.system(cmmd)
    print 'processed sift to', resultname


# SIFT特徴量を読み込んで行列形式で返す
def read_features_from_file(filename):
    f = numpy.loadtxt(filename)
    return f[:,:4],f[:,4:] # 特徴点の配置と記述子


def match_twosided(desc1,desc2):
    """ 双方向対称バージョンのmatch() """

    matches_12 = match(desc1,desc2)
    matches_21 = match(desc2,desc1)

    ndx_12 = matches_12.nonzero()[0]

    # 対称でないものは除去する
    for n in ndx_12:
        if matches_21[int(matches_12[n])] != n:
            matches_12[n] = 0

    return matches_12

def match(desc1,desc2):
    """ 第1の画像の各記述子について、第2の画像の対応点を求める。
    入力:desc1(第1の画像の記述子)、desc2(第2の画像の記述子)"""

    desc1 = numpy.array([d/numpy.linalg.norm(d) for d in desc1])
    desc2 = numpy.array([d/numpy.linalg.norm(d) for d in desc2])

    dist_ratio = 0.6
    desc1_size = desc1.shape

    matchscores = numpy.zeros(desc1_size[0],'int')
    desc2t = desc2.T # あらかじめ転置行列を計算しておく

    for i in range(desc1_size[0]):
        dotprods = numpy.dot(desc1[i,:],desc2t) # 内積ベクトル
        dotprods = 0.9999*dotprods
        # 第2の画像の特徴点の逆余弦を求め、ソートし、番号を返す
        indx = numpy.argsort(numpy.arccos(dotprods))

        # 最も近い近接点との角度が、2番目に近いもののdist_rasio倍以下か?
        if numpy.arccos(dotprods)[indx[0]]<dist_ratio * numpy.arccos(dotprods)[indx[1]]:
            matchscores[i] = int(indx[0])

    return matchscores

def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):
    """ 対応点を線で結んで画像を表示する
    入力: im1,im2(配列形式の画像)、locs1,locs2(特徴点座標)
    machescores(match()の出力)、
    show_below(対応の下に画像を表示するならTrue)"""

    im3 = appendimages(im1,im2)
    if show_below:
        im3 = numpy.vstack((im3,im3))

    plb.imshow(im3)

    cols1 = im1.shape[1]
    for i,m in enumerate(matchscores):
        if m>0: plb.plot([locs1[i][0],locs2[m][0]+cols1],[locs1[i][1],locs2[m][1]],'c')
    plb.axis('off')


def appendimages(im1,im2):
    """ 2つの画像を左右に並べた画像を返す """
    
    # 行の少ない方を選び、空行を0で埋める
    rows1 = im1.shape[0]
    rows2 = im2.shape[0]

    if rows1 < rows2:
        im1 = numpy.concatenate((im1,numpy.zeros((rows2-rows1,im1.shape[1]))),axis=0)
    elif rows1 > rows2:
        im2 = numpy.concatenate((im2,numpy.zeros((rows1-rows2,im2.shape[1]))),axis=0)
    # 行が同じなら、0で埋める必要はない
    return numpy.concatenate((im1,im2), axis=1)


if __name__ == '__main__':
    main()

↑こう書くと↓こう表示されます f:id:end0tknr:20171021093951p:plain

正規化相互相関を用いた Harrisコーナー検出結果の対応点探索 (python)

前回と以前のエントリの続き。というか、これまで理解不足で、さかのぼっていた感じ

Harrisコーナー検出(再び)と、pythonによる実装 - end0tknr's kipple - 新web写経開発

自己相関と相互相関の違いや、正規化相互相関 - end0tknr's kipple - 新web写経開発

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
from PIL import Image
import numpy
import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb
from scipy.ndimage import filters

im1 = numpy.array(Image.open('sf_view1.jpg').convert('L'))
im2 = numpy.array(Image.open('sf_view2.jpg').convert('L'))
wid = 5

def main():
    # 画像1のHariisコーナー検出
    harrisim = compute_harris_response(im1,5)
    filtered_coords1 = get_harris_points(harrisim,wid+1)
    d1 = get_descriptors(im1,filtered_coords1,wid)
    # 画像2のHariisコーナー検出
    harrisim = compute_harris_response(im2,5)
    filtered_coords2 = get_harris_points(harrisim,wid+1)
    d2 = get_descriptors(im2,filtered_coords2,wid)

    # 画像1,2のコーナーの対応点探索
    matches = match_twosided(d1,d2)

    plb.figure()
    plb.gray()
    plot_matches(im1,im2,filtered_coords1,filtered_coords2,matches[:100])
    plb.savefig( '2_1_1.png' )

# グレースケール画像の画素に対応する Harrisコーナー検出器の応答関数
def compute_harris_response(im,sigma=3):
    # 微分係数
    imx = plb.zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)
    imy = plb.zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)
    # Harris行列の成分
    Wxx = filters.gaussian_filter(imx*imx,sigma)
    Wxy = filters.gaussian_filter(imx*imy,sigma)
    Wyy = filters.gaussian_filter(imy*imy,sigma)

    Wdet = Wxx*Wyy - Wxy**2 # 判別式   det()
    Wtr = Wxx + Wyy         # 対角成分 trace()

    return Wdet / Wtr

# Harris応答画像からコーナーを返す。
# min_distはコーナーや画像境界から分離する最小ピクセル数
def get_harris_points(harrisim,min_dist=10,threshold=0.1):
    # 閾値(threshold)を超えるコーナー候補を探す
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1
    # 候補の座標を取得. 「T」は転置method
    coords = plb.array(harrisim_t.nonzero()).T
    # 候補の値を得る
    candidate_values = [harrisim[c[0],c[1]] for c in coords]
    # 候補をsort
    index = plb.argsort(candidate_values)
    # 許容する点の座標を配列に格納
    allowed_locations = plb.zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
    # 最小距離を考慮しながら、最良の点を得る
    filtered_coords = []
    for i in index:
        if allowed_locations[coords[i,0],coords[i,1]] == 1:
            filtered_coords.append(coords[i])
            allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist),
                              (coords[i,1]-min_dist):(coords[i,1]+min_dist)] = 0

    return filtered_coords

# 各点について、点の周辺で幅 2*wid+1 の近傍ピクセル値を返す。
#(点の最小距離 min_distance > wid を仮定する)
def get_descriptors(image,filtered_coords,wid=5):
    desc = []
    for coords in filtered_coords:
        patch = image[coords[0]-wid:coords[0]+wid+1,
                      coords[1]-wid:coords[1]+wid+1].flatten()
        desc.append(patch)

    return desc


# match()の双方向で一致を調べる
def match_twosided(desc1,desc2,threshold=0.5):

    matches_12 = match(desc1,desc2,threshold)
    matches_21 = match(desc2,desc1,threshold)

    ndx_12 = numpy.where(matches_12 >= 0)[0]

    # 非対称の場合を除去する
    for n in ndx_12:
        if matches_21[matches_12[n]] != n:
            matches_12[n] = -1

    return matches_12

# 【正規化相互相関】を用い、第1の画像の各コーナー点記述子について
#  第2の画像の対応点を選択
def match(desc1,desc2,threshold=0.5):
    n = len(desc1[0])

    # 対応点ごとの距離
    d = -numpy.ones((len(desc1),len(desc2)))
    for i in range(len(desc1)):
        for j in range(len(desc2)):
            d1 = (desc1[i] - numpy.mean(desc1[i])) / numpy.std(desc1[i])
            d2 = (desc2[j] - numpy.mean(desc2[j])) / numpy.std(desc2[j])
            ncc_value = sum(d1 * d2) / (n-1)
            if ncc_value > threshold:
                d[i,j] = ncc_value

    ndx = numpy.argsort(-d)
    matchscores = ndx[:,0]
    return matchscores


# 対応点を線で結んで画像を表示する
# 入力: im1,im2(配列形式の画像)、locs1,locs2(特徴点座標)
# machescores(match()の出力)、
# show_below(対応の下に画像を表示するならTrue)
def plot_matches(im1,im2,locs1,locs2,matchscores,show_below=True):

    im3 = appendimages(im1,im2)
    if show_below:
        im3 = numpy.vstack((im3,im3))

    plb.imshow(im3)

    cols1 = im1.shape[1]
    for i,m in enumerate(matchscores):
        if m>0: plb.plot([locs1[i][1],locs2[m][1]+cols1],
                         [locs1[i][0],locs2[m][0]],'c')
    plb.axis('off')

# 2つの画像を左右に並べる
def appendimages(im1,im2):

    # 行の少ない方を選び、空行を0で埋める
    rows1 = im1.shape[0]
    rows2 = im2.shape[0]

    if rows1 < rows2:
        im1 = numpy.concatenate((im1,numpy.zeros((rows2-rows1,im1.shape[1]))),
                                axis=0)
    elif rows1 > rows2:
        im2 = numpy.concatenate((im2,numpy.zeros((rows1-rows2,im2.shape[1]))),
                                axis=0)
    # 行が同じなら、0で埋める必要なし

    return numpy.concatenate((im1,im2), axis=1)

if __name__ == '__main__':
    main()

↑こう書くと↓こう表示されます。数分待ちますが... f:id:end0tknr:20171003214833p:plain

Harrisコーナー検出(再び)と、pythonによる実装

先日までのエントリで「Harrisコーナー検出は理解できた」と思い 実践コンピュータビジョンのサンプルプログラムを写経しようとしたら、 理解度が不十分だった為、やり直し。

テイラー/マクローリン展開からのHarrisコーナー検出 (その1) - end0tknr's kipple - 新web写経開発

テイラー/マクローリン展開からのHarrisコーナー検出 (その2) - end0tknr's kipple - 新web写経開発

実践コンピュータビジョン サンプルプログラム

Harrisコーナー検出の基礎

まず、ある画像上の点xについて、行列 Miを定義します。

 \Large
M_I = \nabla I \ \nabla I^{T} =
\left[ \begin{array}{c} I_x \\ I_y \end{array} \right]
\left[ \begin{array}{cc} I_x & I_y \end{array} \right] =
\left[ \begin{array}{cc} I_x^2   & I_x I_y \\\
                          I_x I_y & I_y^2  \end{array} \right]

ここで、∇については以前のエントリに記載した通りです。

ベクトル偏微分演算子 : ∇ (ナブラ) - end0tknr's kipple - 新web写経開発

また、∇ I は微分係数Ix, Iyから構成され、次のように求められます。

 \Large
I_x = I \ast G_{σx} \ 、\ I_y = I \ast G_{σy}

( 上記において、Gσx, Gσy, * はそれぞれ、標準偏差σのガウス関数と、畳込み記号 )

また、先程のMiにおいては、階数(rank)=1 で、固有値
\lambda_1 = | \nabla I |^2  と \lambda_2 = 0

次にガウス分布 Gσ 等の重み行列を用いて  \Large \overline{M_{I}} を定義します ( Harris行列 ) 。

 \Large
\overline{M_{I}} = W \ast M_I

 \overline{M_{I}} は、 成分ごとの畳み込みで周辺ピクセルの局所平均となります。

( ガウス関数や畳み込みについては、以前のエントリ参照

画像処理 - エッジ抽出フィルタ(sobel, prewitt)と、平滑化フィルタ(gaussian) - end0tknr's kipple - 新web写経開発 )

 \nabla I の値により  \overline{M_{I}} 固有値は変化し、 次のケースに分けられます。

固有値 判定
λ1, λ2 が大きな正値 コーナーあり
λ1 が大きく, λ2 ≒ 0 エッジあり
λ1 ≒0, λ2 ≒ 0 平面

ただし、コーナー判定において固有値を算出すると、 計算コスト:大 の為、次のような判定式を使用します。

 \Large
  \frac { det( \overline{M_{I}} ) }{ trace( \overline{M_{I}} ) }

( det():行列式、trace():転置 )

pythonによるHarrisコーナー検出

実践コンピュータビジョン サンプルプログラム

相変わらずの上記の写経ですが...

#!/usr/local/bin/python
# -*- coding: utf-8 -*-
from PIL import Image
import numpy
import matplotlib
matplotlib.use('Agg')
import matplotlib.pylab as plb
from scipy.ndimage import filters

def main():
    im = numpy.array(Image.open('empire.jpg').convert('L'))
    # Harrisコーナー検出器の応答関数
    harrisim = compute_harris_response(im)
    #
    filtered_coords = get_harris_points(harrisim,6)
    # 検出されたコーナーの表示
    plot_harris_points(im, filtered_coords)

# グレースケール画像の画素に対応する Harrisコーナー検出器の応答関数
def compute_harris_response(im,sigma=3):
    # 微分係数
    imx = plb.zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (0,1), imx)
    imy = plb.zeros(im.shape)
    filters.gaussian_filter(im, (sigma,sigma), (1,0), imy)
    # Harris行列の成分
    Wxx = filters.gaussian_filter(imx*imx,sigma)
    Wxy = filters.gaussian_filter(imx*imy,sigma)
    Wyy = filters.gaussian_filter(imy*imy,sigma)

    Wdet = Wxx*Wyy - Wxy**2 # 判別式   det()
    Wtr = Wxx + Wyy         # 対角成分 trace()

    return Wdet / Wtr

# Harris応答画像からコーナーを返す。
# min_distはコーナーや画像境界から分離する最小ピクセル数
def get_harris_points(harrisim,min_dist=10,threshold=0.1):
    # 閾値(threshold)を超えるコーナー候補を探す
    corner_threshold = harrisim.max() * threshold
    harrisim_t = (harrisim > corner_threshold) * 1
    # 候補の座標を取得. 「T」は転置method
    coords = plb.array(harrisim_t.nonzero()).T
    # 候補の値を得る
    candidate_values = [harrisim[c[0],c[1]] for c in coords]
    # 候補をsort
    index = plb.argsort(candidate_values)
    # 許容する点の座標を配列に格納
    allowed_locations = plb.zeros(harrisim.shape)
    allowed_locations[min_dist:-min_dist,min_dist:-min_dist] = 1
    # 最小距離を考慮しながら、最良の点を得る
    filtered_coords = []
    for i in index:
        if allowed_locations[coords[i,0],coords[i,1]] == 1:
            filtered_coords.append(coords[i])
            allowed_locations[(coords[i,0]-min_dist):(coords[i,0]+min_dist),
                              (coords[i,1]-min_dist):(coords[i,1]+min_dist)] = 0

    return filtered_coords

def plot_harris_points(image,filtered_coords):
    plb.figure()
    plb.gray()
    plb.imshow(image)
    plb.plot([p[1] for p in filtered_coords],[p[0] for p in filtered_coords],'*')
    plb.axis('off')
    plb.savefig( '2_1.png' )

if __name__ == '__main__':
    main()

↑こう書くと↓こう表示されます f:id:end0tknr:20171003212441p:plain