<!--親の顔より見た光景-->

日々の発見を残していきます。

【C++】離散フーリエ級数展開

理解するためにすべて手書きで書いた。
乱数で作ったグラフを、
三角関数の足し合わせのみでキレイに再現できている。
f:id:akagi13213:20190203015531p:plain

(あと一息で離散フーリエ変換ができそう)

#include "pch.h"
#include <iostream>
#include <fstream>
#include <sstream>
#include <vector>

#define M_PI    3.14159265358979323846

int main(void)
{
    std::ifstream ifs( "./source.csv" );
    if ( ifs.fail() ) {
        std::cout << "file read error..." << std::endl;
        return -1;
    }

    std::string str;
    std::vector<double> f;
    while ( getline( ifs, str ) ) {
        f.push_back(std::stod( str ));
    }

    double a0 = 0;
    std::vector<double> an;
    std::vector<double> bn;

    double sum = 0;
    int N = f.size();
    for ( int x = 0; x < N; x++ ) {
        sum += f[x];
    }
    a0 = sum / N;

    for ( int n = 1; n < N; n++ ) {
        double a_sum = 0;
        double b_sum = 0;
        
        for ( int x = 0; x < N; x++ ) {
            a_sum += f[x] * cos( 2 * M_PI * n * x  / N );
            b_sum += f[x] * sin( 2 * M_PI * n * x / N );
        }
        
        if ( n == ( N / 2 - 1 ) ) {
            an.push_back( a_sum / N );
            bn.push_back( b_sum / N );
        }
        else {
            an.push_back( 2 * a_sum / N );
            bn.push_back( 2 * b_sum / N );
        }
    }
    
    std::vector<double> fx_fourier;
    for (int x = 0; x < N; x++) {
        double a_sum = 0;
        double b_sum = 0;
        for ( int n = 1; n < N; n++ ) {
            a_sum += an[n - 1] * cos( 2 * M_PI * n * x / N );
            b_sum += bn[n - 1] * sin( 2 * M_PI * n * x / N );
        }
        fx_fourier.push_back( a0 + a_sum + b_sum );
    }

    std::ofstream dist;
    dist.open("dist.csv", std::ios::out);
    for ( int i = 0; i < fx_fourier.size(); i++ ) {
        dist << fx_fourier[i] << std::endl;
    }
    
    return 0;
}

【画像処理】2次微分フィルタ

二次微分フィルタを自作した。

def edge_filter_2(file):
    image = cv2.imread(file, cv2.IMREAD_GRAYSCALE)
    (height, width) = image.shape[:2]

    edge_width = image.copy()
    edge_height = image.copy()
    edge_lap = image.copy()
    for i in range(1, height-1):
        for j in range(1, width-1):
            left = image[i, j-1]
            right = image[i, j+1]
            top = image[i-1, j]
            bottom = image[i+1, j]
            center = -image[i, j] * 2

            edge_width[i, j] = abs(left + right + center)
            edge_height[i, j] = abs(top + bottom + center)
            edge_lap[i, j] = abs(top + bottom + left + right + center*2)
 
    cv2.imwrite('gray.png', image)
    cv2.imwrite('edge2_width.png', edge_width)
    cv2.imwrite('edge2_height.png', edge_height)
    cv2.imwrite('edge2_lap.png', edge_lap)

参考は以下。
画像処理フィルタの1次微分と2次微分の違い | ぱーくん plus idea

適用フィルタは以下。
横方向

0 0 0
1 -2 1
0 0 0

縦方向

0 1 0
0 -2 0
0 1 0

4近傍

0 1 0
1 -4 1
0 1 0

●横方向二次微分
f:id:akagi13213:20190120184119p:plain

●縦方向二次微分
f:id:akagi13213:20190120184130p:plain

●4近傍(ラプラシアンフィルタ)
f:id:akagi13213:20190120184140p:plain

★なぜ1次微分より2次微分のほうがエッジが目立つ?
f:id:akagi13213:20190120184749p:plain
左から、元画像・1次微分(横方向)・2次微分(横方向)

横方向微分の1次微分と2次微分の出力差を見てみる。
仮に、顔の輪郭にフィルタを適用するとして、
フィルタ適用前画素列を、以下とする

32 32 64 64 64 64 255 255

★★1次微分(横方向)

-1/2 0 1/2

このフィルタを適用すると(左端・右端は元画像のままとする)、

32 32 64 64 64 64 255 255

 ↓↓ 1次微分適用

32 16 16 0 0 95 95 255

★★2次微分(横方向)

1 -2 1

このフィルタを適用すると(左端・右端は元画像のままとする)、
※絶対値を画素値として採用

32 32 64 64 64 64 255 255

 ↓↓ 2次微分適用

32 32 -32 0 0 191 -191 255

1次微分より白成分が2倍になっている。

【画像処理】1次微分フィルタ(横、縦)

自作で微分フィルタを作った。

#!/usr/env/bin python
import cv2
import sys

args = sys.argv

def edge_filter(file):
    image = cv2.imread(file, cv2.IMREAD_GRAYSCALE)
    image = cv2.medianBlur(image, 3)
    (height, width) = image.shape[:2]

    edge_width = image.copy()
    edge_height = image.copy()
    for i in range(1, height-1):
        for j in range(1, width-1):
            left = image[i, j-1] / -2
            right = image[i, j+1] / 2
            top = image[i-1, j] / -2
            bottom = image[i+1, j] / 2

            edge_width[i, j] = abs(left + right)
            edge_height[i, j] = abs(top + bottom)

    cv2.imwrite('gray.png', image)
    cv2.imwrite('edge_width.png', edge_width)
    cv2.imwrite('edge_height.png', edge_height)

if __name__ == '__main__':
     edge_filter(args[1])

フィルタは以下を適用した。
ただし、画像の端には適用していない。(コード上の、range(1, height-1) と range(1, width-1))

●横方向微分フィルタ

0 0 0
-1/2 0 1/2
0 0 0

●縦方向微分フィルタ

0 -1/2 0
0 0 0
0 1/2 0

■元画像
f:id:akagi13213:20190120155743p:plain

■横方向微分画像
f:id:akagi13213:20190120155754p:plain

■縦方向微分画像
f:id:akagi13213:20190120155808p:plain

余計なエッジを消すために、メディアンフィルタを適用してから微分フィルタを適用する。

メディアンフィルタ 方向微分
f:id:akagi13213:20190120155836p:plain

メディアンフィルタ 方向微分
f:id:akagi13213:20190120155843p:plain

★★比較
f:id:akagi13213:20190120161812j:plain

縦方向微分は横線を抽出し、
横方向微分は縦線を抽出している。

【plantuml】ピリオドをクラス名として扱わせる

クラス図にヘッダファイルを登場させたくなった。
.hの、.が namespaceとして扱われるようで、思い通り描いてくれない。
f:id:akagi13213:20190113220705p:plain
↑こうなっちゃう。

class goal_functions as "goal_functions.h" {
    +double getGoalPositionDistance(const gemetry_msgs::PoseStamped& global_pose
}

上のソースのように、

class 適当な名前 as "描画させたい名前" {

の形式で記述すれば、ピリオドをクラス名に含ませることが可能。
f:id:akagi13213:20190113220600p:plain

ヘビの安否確認システム

帰省時、蛇は雪国に持っていけないので安否システムを作った。

Python と USBカメラだけでヘビが動いてるか(生きてるか)を、自宅から離れていても確認出来る。

f:id:akagi13213:20181225000054j:image


f:id:akagi13213:20181225000911j:image

少しでも動いたらgmailに画像添付して通知を送ってくれる。監視はとりあえず1分毎。

 

仕組みは簡単で、前回の画像と現在の画像をグレースケール化、走査し、輝度差がしきい値以上なら差分としてカウント。

1分前から変化があった画素は緑にした。上手く蛇のとこだけ検知できてる。うちのヘビは白いので閾値が効きやすい…


最初はAKAZEとかORBとか、特徴点マッチングを試したけど、水槽の定点カメラでは似すぎて厳しかった。

 

https://github.com/akagi132/basilisk