日曜技術者のメモ

趣味でやった事のメモ書きです。

Vivado IP Integratorでap_fifoとfifoを簡単に接続する

ap_fifofifoを接続するときの小ネタ。

ap_fifofifoと接続されない

Vivado HLSでap_fifoを多用する私ですが、Vivado IP Integrator上でFIFO Generatorと接続しようとしてもそのままでは繋がりません。
f:id:ginnyu-tei:20170717120816j:plain:w400
↑の様に接続候補に現れない。

ap_fifoは何故かfullやempty信号が負論理になっているので反転する必要があります。

f:id:ginnyu-tei:20170717121804j:plain:w400
↑この様な形で配線が必要。

この配線を毎回するのは面倒な上にせっかくIP IntegratorでFIFO I/Fがまとめてられているのに毎度ばらして配線するのは見た目が悪いので変換するIPを作成した。
これを使うと↓の様に配線がすっきりする。
f:id:ginnyu-tei:20170717123204j:plain:w400

ap_fifo2fifo

まずIPにするRTLだが↓の様にFULLを反転するコードを用意する。
これでIPを作成すれば良い。

module ap_fifo2fifo(
    WR_DATA_o,
    FULL_N_o,
    WR_EN_o,
    WR_DATA_i,
    FULL_i,
    WR_EN_i
);

parameter WIDTH = 32;

//ap_fifo
input [WIDTH-1:0] WR_DATA_o;
output FULL_N_o;
input WR_EN_o;

//fifo
output [WIDTH-1:0] WR_DATA_i;
input FULL_i;
output WR_EN_i;


assign WR_DATA_i = WR_DATA_o;
assign FULL_N_o = ~FULL_i;
assign WR_EN_i = WR_EN_o;

endmodule

このままIPを生成すると各ポートがそのままになるのでI/F定義が必要になる。
f:id:ginnyu-tei:20170717124414j:plain:w400

din側はInterface Definitionを"acc_fifo_write_rtl"のslaveにして各ポートを割り当てる。
f:id:ginnyu-tei:20170717124525j:plain:w400
f:id:ginnyu-tei:20170717124718j:plain:w400
dout側はInterface Definitionを"fifo_write_rtl"のmasterにして各ポートを割り当てる。
f:id:ginnyu-tei:20170717124838j:plain:w400 f:id:ginnyu-tei:20170717124851j:plain:w400

fifo2ap_fifo

ap_fifo2fifoもemptyを反転するRTLを用意してIPを作成する。

module fifo2ap_fifo(
    RD_DATA_i,
    EMPTY_i,
    RD_EN_i,
    RD_DATA_o,
    EMPTY_N_o,
    RD_EN_o
);

parameter WIDTH = 32;

//fifo
input [WIDTH-1:0] RD_DATA_i;
input EMPTY_i;
output RD_EN_i;

//ap_fifo
output [WIDTH-1:0] RD_DATA_o;
output EMPTY_N_o;
input RD_EN_o;

assign RD_DATA_o = RD_DATA_i;
assign EMPTY_N_o = ~EMPTY_i;
assign RD_EN_i = RD_EN_o;

endmodule

こちらもI/F定義をやっていく。
f:id:ginnyu-tei:20170717130240j:plain:w400

din側はInterface Definitionを"fifo_read_rtl"のmasterにして各ポートを割り当てる。 f:id:ginnyu-tei:20170717130343j:plain:w400
f:id:ginnyu-tei:20170717130354j:plain:w400
dout側はInterface Definitionを"acc_fifo_read_rtl"のslaveにして各ポートを割り当てる。
f:id:ginnyu-tei:20170717130425j:plain:w400
f:id:ginnyu-tei:20170717130439j:plain:w400

後はVivadoのプロジェクト上から各IPをインスタンスすれば簡単に配線できる。

VivadoHLSで射影変換を実装してみた-3-

↓の続きで今回が多分射影変換の最終回

se.hatenablog.jp

原因追求

前回射影変換をお手軽に実装したがかなり遅かったのでまずは原因を確かめる。
と、言っても原因は↓のAXI HP3の波形を見れば明らか。
f:id:ginnyu-tei:20170716121548j:plain:w400
Write/Read共にシングルアクセスなので1pix読み書きに時間がかかっている。なのでピクセルクロックの3倍で動作させても遅いのである。

しかし、射影変換は座標を変換するのでアドレスが飛び飛びになる。これ自体はどうしようもないので他の方法で速度アップを図る。

クロックアップ(150MHz→250MHz)

まず、手抜き手軽な対応としてバス周りのクロックアップをやってみる。homography_dmaとHP3バス周りを250MHzまで上げてみた。
ちなみにhomography_dmaは300MHzまで合成できたが、PSからPLへ供給できるクロックが250MHzまでだったので250MHzにした。
その結果が約33fpsまで上がった。
f:id:ginnyu-tei:20170716211315j:plain:w400

ただし、AXI HP3の波形はあまり変わっていない。
f:id:ginnyu-tei:20170716215304j:plain:w400

ライトアクセス改善

次にライトアクセスを改善してみる。
上でも書いたとおり射影変換のアドレスが飛び飛びになるが今回の実装では、画像データライトは連続している。
そこで、homography_dmaを改造してライトはバースト転送にする。

f:id:ginnyu-tei:20170716214224j:plain
↑の様にリードしたデータを一旦FIFOにプッシュして別モジュールからバーストライトを実行する。
これでDRAMへの負荷軽減を期待している。

では、コードを記載していく。
まずはhomography_dma_write_fifoここはo_dataのm_axiをap_fifo変えるだけ。この様にI/Fを手軽に変更できるのも高位合成の良い所。

void homography_dma_write_fifo(
        unsigned int *i_data,
        unsigned int *o_data,
        unsigned int *addr
        ){
#pragma HLS INTERFACE ap_fifo depth=786432 port=addr
#pragma HLS INTERFACE m_axi depth=786432 port=i_data
#pragma HLS INTERFACE ap_fifo depth=786432 port=o_data

    unsigned int tmp;
    for (int i = 0 ; i < 1024*768; i++){
#pragma HLS PIPELINE
        tmp = *(addr + i);
        if (tmp > 786432){
            o_data[i] = 0;
        } else {
            o_data[i] = i_data[tmp];
        }
    }
}

次にfifo2axi_xgaのコード

void fifo2axi_xga(
        unsigned int *i_data,
        unsigned int *o_data
        ){
#pragma HLS INTERFACE ap_fifo depth=786432 port=i_data
#pragma HLS INTERFACE m_axi depth=786432 port=o_data

    for (int i = 0 ; i < 1024*768; i++){
        o_data[i] = i_data[i];
    }
}

ただfifoのデータをaxiに書き込むだけである。ここはバーストlengthを指定できないか色々やってみたが、この書き方が16バースト固定で一番速度が早かった。
Vivado HLSでCosimした結果が↓
f:id:ginnyu-tei:20170716223136j:plain:w400

実際にはFIFOにライトする速度の方が早いのでfifoにある程度データが溜まってからライトを開始するようにすれば効率よくAXIライト出来ると思う。
(これ書いてて思いついたが、ap_fifoのempty信号にAlmost Fullを使ってバーストlength分溜まったら信号を出す様にすれば必ず1バースト分は連続で出力できる。でも、1フーレムのデータ数と1バーストの転送量の関係によってはFIFOに半端なデータが残りそう)

これで必要なモジュールができたのでIP Integratorで配線した結果が↓
f:id:ginnyu-tei:20170717001902j:plain:w400

ブロックがかなり増えたのでCreate Hierarchyで階層化した。
この状態でZedboardでhomography_dma_write_fifoのap_doneを観察した結果が↓
f:id:ginnyu-tei:20170717004256j:plain:w400
約47FPSとかなり改善した。
実際はap_doneの後fifo2axi_xgaがまだ動いているのがパイプライン的な動作になっているのでFPSには影響しない(はず)。
(例えパイプライン動作でなくとも残っているは1バースト+α位で多少増えるだけ)
AXI HP3の波形を確認するとリード側はあまり変化していないが、ライト側はバースト転送になっているのがわかる。
f:id:ginnyu-tei:20170717010613j:plain:w400

さらなる高速化

以上で簡単にできそうな改良をして2倍位早くなった。
これ以上早くするアイディアはいくつかあるが実装に時間がかかりそうなのでメモだけしておく。

  1. ブロックRAMを使う
    フレームデータをDRAMでなくブロックRAMに格納する。これならRead/Writeが1サイクルで可能なので高速化できる。でも、ブロックRAMが多いFPGAじゃないと足りなくなる。
  2. 先にアドレス計算をしておく
    一度パラメータを決めたらしばらく変えなくて良いアプリケーション向けだが、先に1フレーム分座標を計算してバースト転送出来るように並び替えをする。
  3. 変換する範囲をなるだけ狭くする
    今回の様にフレーム分の一部のみ変換する場合はフレーム全て変換せず、変換前/変換後の各4座標から変換に必要な範囲だけを転送するようにする。
  4. ブロックサイズでリードする。
    変換する次のアドレスは今のアドレスの近くにある可能性が高いので8×8や16×16のサイズでリードしてデータを貯めておく。

とりあえず、この位思いついた。実装は3以外大変そうなので誰かこれを試してみてほしい。
また思いついたら書き足す。

Zynqのベアメタルアプリケーションで画像をメモリに書き込む

Zynqでベアメタルアプリケーション(OSなし)を走らせる際
テスト入力として画像データをメモリ上に展開したくなる時がある。
その際GIMPを使えば簡単に画像データをソースコードに変換できるので
blogに書いておく。

GIMP

まず、GIMPをインストール&起動する。
そしてメモリに展開したい画像を開く

開いた後、ファイル→名前をつけてエクスポートを選ぶ、
そして、右下のエクスポートされた全ての画像から
「Cソースコードヘッダ」を選んで適当に名前をつけてエクスポートする。
f:id:ginnyu-tei:20170709005400j:plain

(かなり大きなヘッダファイルが出力されるので開く時は注意!)

出力されたファイルを確認すると以下の様なコードが出力される。

static unsigned int width = 1024;
static unsigned int height = 768;

/*  Call this macro repeatedly.  After each use, the pixel data can be extracted  */

#define HEADER_PIXEL(data,pixel) {\
pixel[0] = (((data[0] - 33) << 2) | ((data[1] - 33) >> 4)); \
pixel[1] = ((((data[1] - 33) & 0xF) << 4) | ((data[2] - 33) >> 2)); \
pixel[2] = ((((data[2] - 33) & 0x3) << 6) | ((data[3] - 33))); \
data += 4; \
}
static char *header_data =
    "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
    "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
    "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
    "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!"
~以下画像データが続く~

SDK

ヘッダファイルができたらXilinx SDKでヘッダファイルをincludeする。
HEADER_PIXELを呼び出す毎にラスタースキャンでRGB色が取得できるので
それを直接メモリに書き込む。

#include <stdio.h>
#include "platform.h"
#include "xil_printf.h"

#include "test.h" //GIMPでエクスポートしたヘッダファイル


int main()
{
    int i,j;
    char pixel[3]; //RGB格納配列
    unsigned int R,G,B;

    volatile unsigned int *Addr;
    volatile unsigned int *reg_base;

    Addr = (unsigned int*)0x10000000;

    init_platform();

    Xil_DCacheDisable();


    print("Hello World\n\r");

    for (j = 0 ; j < 768;j++){
        for (i = 0 ; i < 1024;i++){
            HEADER_PIXEL(header_data,pixel); //1pix読み込む
            R = pixel[0] & 0xFF;
            G = pixel[1] & 0xFF;
            B = pixel[2] & 0xFF;
            Addr[j*1024+i]= (R << 16) | (G << 8) | (B); //メモリへ書き込み
            //Addr[j*1024+i]=j*1024+i;
        }
    }
~以下省略~

この手法はjtagで画像をFPGAに転送するので大きな画像の場合
少し時間がかかるのが難点だが手軽に使えるのが利点。

VivadoHLSで射影変換を実装してみた-2-

前回の続き

se.hatenablog.jp

Vivado HLS

射影変換は以下図のように2モジュールに分けて実装をした。
f:id:ginnyu-tei:20170708185251j:plain:w200

fifo使ってモジュールを小分けにすると
実装が楽になったり速度が出やすくなるので良く使います。
ここら辺の話は需要があればblogかどっかで纏めるかも。

homography

まずはhomographyのC++コード

void homography(
        int a,
        int b,
        int c,
        int d,
        int e,
        int f,
        int g,
        int h,
        unsigned int *addr)
{
#pragma HLS INTERFACE s_axilite port=a bundle=reg
#pragma HLS INTERFACE s_axilite port=b bundle=reg
#pragma HLS INTERFACE s_axilite port=c bundle=reg
#pragma HLS INTERFACE s_axilite port=d bundle=reg
#pragma HLS INTERFACE s_axilite port=e bundle=reg
#pragma HLS INTERFACE s_axilite port=f bundle=reg
#pragma HLS INTERFACE s_axilite port=g bundle=reg
#pragma HLS INTERFACE s_axilite port=h bundle=reg

   #pragma HLS INTERFACE ap_fifo port=addr

    unsigned int u,v;

    for (int y = 0 ; y < 768 ; y++){
        for (int x = 0;x < 1024; x++){
#pragma HLS PIPELINE
            u = (x*a + y*b + c) / (x*g + y*h + 65536);
            v = (x*d + y*e + f) / (x*g + y*h + 65536);
            int tmp = v*1024+u;
            if ((v >= 0) & (v < 768) & (u >= 0) & (u < 1024)){
                *addr = tmp;
            } else {
                *addr = 786433; //786432 + 1
            }
        }
    }
}

↑の図では書いていなかったがlパラメータa~hはAXI LiteのレジスタにしてCPUから設定。
射影変換の式はC++コードそのままで、変換後座標を画像左上から順に変化させ
変換前座標を計算し、そこからメモリアドレスを計算してFIFOに入れる。
変換後座標値の値によっては範囲外の変換前座標になる事があるので
その時は範囲外という意味で1024*768+1の値をFIFOに入れる。
ここではパイプラインオプションを入れているので毎サイクル計算結果が出力される。

modelsimでシミュレーションしてみると毎サイクル出てくるのが確認できる。
f:id:ginnyu-tei:20170708190530j:plain:w200

そして回路規模
f:id:ginnyu-tei:20170708203559j:plain
積和除算があるのでDSPはそこそこ使ってる。

homography_dma

次はhomography_dmaのコード

void homography_dma(
        unsigned int *i_data,
        unsigned int *o_data,
        unsigned int *addr
        ){
#pragma HLS INTERFACE ap_fifo port=addr
#pragma HLS INTERFACE m_axi depth=786432 port=i_data
#pragma HLS INTERFACE m_axi depth=786432 port=o_data

    unsigned int tmp;
    for (int i = 0 ; i < 1024*768; i++){
#pragma HLS PIPELINE
        tmp = *addr;
        if (tmp > 786432){
            o_data[i] = 0;
        } else {
            o_data[i] = i_data[tmp];
        }
    }
}

homographyで計算したアドレスを使ってo_dataから出力するだけ。
元画像の範囲外座標は黒色を出力する。
なお、今回は実装していないがi_data[tmp]だけリードせず、
例えば2x2pixを読み込みしてバイリニアで補完したり、
4x4pixでバイキュービックで補完すると画質が良くなるはずである。
(この際はhomographyから1次元アドレスではなく座標値を送る用変更が必要)

↓がSim結果
AXIリード・ライト共にシングルアクセスである。
f:id:ginnyu-tei:20170708214451j:plain:w200

↓回路規模
f:id:ginnyu-tei:20170708214833j:plain

Vivado

これでモジュールが完成したので実機で試す。
TOPの配線は全てIP Integratorで済ます。
↓が全て配線したTOP
f:id:ginnyu-tei:20170708215611j:plain:w200

↓概略図
f:id:ginnyu-tei:20170708221823j:plain:w200

まず、0x10000000にCPUから画像データを書き込む。
その後、0x10000000からhomography_dmaモジュールが画像を読みこみ0x11000000に書き出す。
アナログRGBモジュールは0x10000000にある画像データを出力する。
なお、上図には書いていないが書くモジュールのap_startはVIOで制御している。
テスト稼働時はVIOを使うと非常に楽である。

本当は射影変換とアナログRGBモジュールのFPSが異なるので
フレームの同期が取れていないが今回はそのままにしておく。
静止画であれば気にならないと思う。

SDK

以下がPSで動かすプログラムである。
テンプレートにあるHello Worldを改変して使っている。

#include <stdio.h>
#include "platform.h"
#include "xil_printf.h"

#include "test.h"

int main()
{
    int i,j;
    char pixel[3];
    unsigned int R,G,B;

    volatile unsigned int *Addr;
    volatile unsigned int *reg_base;

    Addr = (unsigned int*)0x10000000;
    reg_base   =  (unsigned int*)0x43C00000;


    init_platform();

    Xil_DCacheDisable();

    print("Hello World\n\r");

    for (j = 0 ; j < 768;j++){
        for (i = 0 ; i < 1024;i++){
            HEADER_PIXEL(header_data,pixel);
            R = pixel[0] & 0xFF;
            G = pixel[1] & 0xFF;
            B = pixel[2] & 0xFF;
            Addr[j*1024+i]= (R << 16) | (G << 8) | (B);
            //Addr[j*1024+i]=j*1024+i;
        }
    }

    /射影変換パラメータ設定
    reg_base[4]  = 0x2C6C6;
    reg_base[6]  = 0xF6EE;
    reg_base[8]  = 0xFCD32CF7;
    reg_base[10] = 0x248A;
    reg_base[12] = 0x2C492;
    reg_base[14] = 0xFE1BDF57;
    reg_base[16] = 0x49;
    reg_base[18] = 0x2A;

    cleanup_platform();
    return 0;
}

ここで画像データをCPUに書き込むが、
その画像データはGIMPで作成している。 これについては別途blogに書く予定。

Zedboardでテスト

Zedboardで動かした結果が↓
f:id:ginnyu-tei:20170705210623j:plain:w200
前回のソフトで動作させた結果とほぼ同じ形に変形している。

この状態でフレームレートを測定した。
フレームレートの測定はhomography_dmaのap_doneを
Pmodに出力しオシロスコープで間隔を測定した。
↓がオシロスコープの画面
f:id:ginnyu-tei:20170708225509j:plain:w200
ap_doneは1パルスしか出なくオシロで取り逃す可能性がある
そこでap_done毎に信号を反転するモジュールを入れている。
その結果、約24fps出ている事が分かった。
Twitterでは約11fpsとツイートしたが(↓)あれは測定ミスでした。)

以上で射影変換の実装&実機テストは完了である。
今回は射影変換をお手軽に実装するという脳内テーマで
スタートしたので速度や規模はあまり気にしていなかった。
(無事動いたのでOKとする)
とは言え、homography_dmaは150MHz(ピクセルクロックの3倍)で動作させているのに
24フーレムしか出ていないので、次回は原因についてと高速化について言及する。

VivadoHLSで射影変換を実装してみた-1-

台形補正のアルゴリズムを調べてたら射影変換でやるみたいなのでFPGAに実装してみた。

OpenCV

まずはOpenCVでさくっとテスト

#include <opencv2/opencv.hpp>
#include <opencv2/highgui/highgui.hpp>
int main(void)
{
    //cv::Mat src_img;
    //src_img = cv::imread("test.jpg", 1);
    IplImage *pImg;
    pImg = cvLoadImage("test.jpg");
    //IplImage -> Mat
    cv::Mat mat(pImg);
    cv::Mat src_img;
    src_img = pImg;

    cv::Mat dst_img = src_img.clone();

    if (src_img.empty()){
        return -1;
    }

    const cv::Point2f src_pt[]={
        cv::Point2f(256, 128),
        cv::Point2f(768, 128),
        cv::Point2f(768, 640),
        cv::Point2f(256, 640) };

    const cv::Point2f dst_pt[]={
        cv::Point2f(356, 228),
        cv::Point2f(768, 228),
        cv::Point2f(668, 640),
        cv::Point2f(256, 540)};

        //パラメータ計算
    cv::Mat homography_matrix = cv::getPerspectiveTransform(src_pt,dst_pt);
        //射影変換
    cv::warpPerspective( src_img, dst_img, homography_matrix,src_img.size());

    cv::namedWindow("Image", CV_WINDOW_AUTOSIZE | CV_WINDOW_FREERATIO);
    cv::imshow("Image", dst_img);
    cv::waitKey(0);
}

f:id:ginnyu-tei:20170705220240j:plain f:id:ginnyu-tei:20170705220131j:plain こんな風に変換前の4点座標と変換後の4点座標を指定することで
座標に合うように画像を変形してくれます。

C++で実装

C++での実装は↓のページを丸パクリ参考にしました。
mf-atelier.sakura.ne.jp

C++で実装する際、ついでに16ビットシフトして整数演算化しました。
多分小数点そこまでなくても問題ないと思う

void homography(
    int x,
    int y,
    int a,
    int b,
    int c,
    int d,
    int e,
    int f,
    int g,
    int h,
    int *u,
    int *v)
{
    *u = (x*a + y*b + c) / (x*g + y*h + 65536);
    *v = (x*d + y*e + f) / (x*g + y*h + 65536);
}


int main(void)
{

    double a, b, c, d, e, f, g, h;
    int i_a, i_b, i_c, i_d, i_e, i_f, i_g, i_h;
    int u, v;

    int    naOrig[4][2] = { { 356, 228 }, { 768, 228 }, { 668, 640 }, { 256, 540 } };
    int    naTran[4][2] = { { 256, 128 }, { 768, 128 }, { 768, 640 }, { 256, 640 } };

    //cv::Mat src_img;
    //src_img = cv::imread("test.jpg", 1);

    IplImage *pImg;
    pImg = cvLoadImage("test.jpg");
    //IplImage -> Mat
    cv::Mat mat(pImg);
    cv::Mat src_img;
    src_img = pImg;

    cv::Mat dst_img = src_img.clone();

    if (src_img.empty()){
        return -1;
    }

    homography_param
        (
        naOrig,
        naTran,
        &a,
        &b,
        &c,
        &d,
        &e,
        &f,
        &g,
        &h
        );

    std::cout <<
        "a=" << a << std::endl <<
        "b=" << b << std::endl <<
        "c=" << c << std::endl <<
        "d=" << d << std::endl <<
        "e=" << e << std::endl <<
        "f=" << f << std::endl <<
        "g=" << g << std::endl <<
        "h=" << h <<
        std::endl;

    i_a = (int)(a * 65536);
    i_b = (int)(b * 65536);
    i_c = (int)(c * 65536);
    i_d = (int)(d * 65536);
    i_e = (int)(e * 65536);
    i_f = (int)(f * 65536);
    i_g = (int)(g * 65536);
    i_h = (int)(h * 65536);

    std::cout <<
        "i_a=" << i_a << std::endl <<
        "i_b=" << i_b << std::endl <<
        "i_c=" << i_c << std::endl <<
        "i_d=" << i_d << std::endl <<
        "i_e=" << i_e << std::endl <<
        "i_f=" << i_f << std::endl <<
        "i_g=" << i_g << std::endl <<
        "i_h=" << i_h <<
        std::endl;

    cv::Vec3b *src = src_img.ptr<cv::Vec3b>(0);
    cv::Vec3b *dst = dst_img.ptr<cv::Vec3b>(0);

    for (int y = 0; y < 768; y++){
        for (int x = 0; x < 1024; x++){
            homography(x, y, i_a, i_b, i_c, i_d, i_e, i_f, i_g, i_h, &u, &v);
            if ((u >= 0) & (u < 1024) &
                (v >= 0) & (v < 768)){
                dst[y * 1024 + x] = src[v * 1024 + u];
            }
            else {
                dst[y * 1024 + x] = cv::Vec3b(0, 0, 0);
            }
        }
    }
    cv::namedWindow("Image", CV_WINDOW_AUTOSIZE | CV_WINDOW_FREERATIO);
    cv::imshow("Image", dst_img);
    cv::waitKey(0);
}

パラメータ変換式は参考元のコードをほぼそのまま使用した。 (関数名をhomography_paramに変更)
逆行列計算はOpenCVを使って↓の用に実装

 // 逆行列
    //matinv(8, dATA, dATA_I);
    cv::Mat matA = (cv::Mat_<double>(8, 8));

    for (int j = 0; j < 8; j++){
        for (int i = 0; i < 8; i++){
            matA.at<double>(i, j) = dATA[i][j];
        }
    }

    matA = matA.inv();

    for (int j = 0; j < 8; j++){
        for (int i = 0; i < 8; i++){
            dATA_I[i][j] = matA.at<double>(i, j);
        }
    }

パラメータをどうやって計算するか知りたかったので数式を見たら
8次元連立一次方程式で自分で解くのは諦めました・・・

それと大事な事なんですが、ここでは変換後座標から変換前座標を計算したいので
変換前座標と変換後座標を入れ替えて計算しています。

こうする事で変換後画像内に変換前画像が当てはまらず
ピクセルが抜けてしまう事が防げます。
その為、ピクセルを補完をしなくても良くなります。
(正確には一番近い画素を変換後座標にコピーしているので
ニアレストネイバーを使ってる言える。)

余談 射影変換の逆計算方法が初め分からず逆行列の計算を
するのかとげんなりしていましたが以下ページを見て解決しました。

detail.chiebukuro.yahoo.co.jp

長くなったのでVivado HLSのコードは次回

多層パーセプトロンをPythonで実装してみた-手書き数字認識-

前回(↓)の続き。

se.hatenablog.jp

簡単な分類で動作を確認できたので次は数字認識をやってみる。

多層パーセプトロンをPythonで実装してみた-1- - 日曜技術者のメモ
で書いたクラスは出力層が二個以上に対応していなかったので修正して使った。
ブログは修正済み。

数字の分類では出力層にソフトマックス関数を使って出力を
[-1,1]ではなく確立にするのが普通の様だが今回は出力層をそのままにして
学習データは正しい数字は0.5、その他の数字は0として学習した。

sklearn.datasets

sklearn.datasetsに手書き文字があったのでまずはそれでやってみる。

The Digit Dataset — scikit-learn 0.18 documentation

8x8の数字なのでかなり荒い。

入力層は8pix×8pix=64個
出力層は0から9の10個
学習率0.03
600文字学習を繰り返す。

実装コード

def digit_sepalate():
  #テストデータ作成
  iteration_num = 100000
  sample_num = 600

  digits = load_digits()

  x = digits.data
  y = digits.target
  x_max = x.max()

  x /= x.max()
  x *= 0.5

  iLayer = nn.inputLayer(64,10,10)
  oLayer = nn.outputLayer(64,10,1,0.05)

  inputLayerY = np.zeros(64)
  dut_output = np.zeros(10)

  x1_array = np.arange(0,8,1)
  x2_array = np.arange(0,8,1)
  pos_x1,pos_x2 = np.meshgrid(x1_array,x2_array)
  val_y = np.identity(40)


  plt.ion()

  match_list = []
  cnt_list = []

  shuffle_list = list(range(sample_num))

  plt.pause(0.05)

  for i in range(iteration_num):
    match = 0
    for j in range(sample_num):
      inputLayerY = x[j]

      iLayer.foward(inputLayerY)
      oLayer.foward(iLayer.getY())

      exp_y = np.zeros(10) - 1


      dut_output = oLayer.getY()

      exp_y[y[j]] = 1.0

      oLayer.backprop(exp_y)

      oLayer.updateWeight()

      if(dut_output.argmax() == digits.target[j]):
         match += 1
      
    #次の学習時に順番をシャッフルする
    random.shuffle(shuffle_list)
    
    if (1):
      plt.clf()
      match_list.append(match/sample_num)
      cnt_list.append(i)

      for i in range(10):
        plt.subplot(3,5,6+i)
        plt.xlim(0,7)
        plt.ylim(0,7)
        plt.pcolor(pos_x1,pos_x2,oLayer._weight[i].reshape(8,8))

      plt.subplot(3,1,1)
      plt.plot(cnt_list,match_list)
      plt.pause(0.05)


  while True:
    plt.pause(0.05)

↓は実行結果
f:id:ginnyu-tei:20161101215515j:plain:w400
10個のヒートマップは各出力層の重み値を出力している。
うまくいけば数字に見えるが・・・微妙。
数回の学習でかなりの精度になった。

MNIST

次は有名なMNISTでやってみる。

データは以下にあるがsklearnで関数が容易されているのでそっちを使った。
MNIST handwritten digit database, Yann LeCun, Corinna Cortes and Chris Burges

入力層は28pix×28pix=784個
出力層は0から9の10個
学習率0.03
1000文字学習を繰り返す。

↓は実行結果
f:id:ginnyu-tei:20161101220631j:plain:w400

10回の学習で正解率91%!!
幾ら何でも学習速すぎないかと思っていたが
これは学習した数字で学習率を見てるので速いのは当たり前だった

次は学習する数字と正解率を評価する数字を別々にしてみる。
正解率は100文字中の正解数

1000文字学習
f:id:ginnyu-tei:20161101224349j:plain:w400

5000文字学習
f:id:ginnyu-tei:20161101224406j:plain:w400

10000文字学習
f:id:ginnyu-tei:20161101224421j:plain:w400

50000文字学習
f:id:ginnyu-tei:20161101225229j:plain:w400

さっきよりは学習に時間がかかっているが問題なく分類できてる
重みのヒートマップも0,3,8はそれっぽい。

多層パーセプトロンをPythonで実装してみた-2-

前回(↓)の続き。

se.hatenablog.jp

多層パーセプトロン各層の実装ができたので実際に動かしてみる。
パーセプトロンの数や層の構成は↓を参考にした。

googlecloudplatform-japan.blogspot.jp

簡単な分類

直線分類

まず↓をやってみる。

f:id:ginnyu-tei:20161019194257g:plain:w400

実装は↓

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import random as random
from sklearn.datasets import make_blobs
from sklearn.datasets import make_circles

import NNLayer as nn

def linearlySeparable():
  
  #学習数、ループ回数
  iteration_num = 10
  sample_num = 100
  
  #学習結果表示グラフ設定
  x1_array = np.arange(-8,8,0.4)
  x2_array = np.arange(8,-8,-0.4)
  pos_x1,pos_x2 = np.meshgrid(x1_array,x2_array)
  val_y = np.identity(40)

  shuffle_list = list(range(sample_num))
  match_list = []
  cnt_list = []
  inputLayerY = np.zeros(2)

  #NN作成
  iLayer = nn.inputLayer(2,2,1)
  oLayer = nn.outputLayer(2,1,1,0.03)

  #テストデータ作成
  centers = [(3,3),(-3,-3)]
  #centers = [(-3,-3),(3,3)]
  #centers = [(0,-3),(0,3)]
  #centers = [(-3,0),(3,0)]
  x_exp,y_exp = make_blobs(n_samples= sample_num,centers=centers,random_state=42,cluster_std=1.5)

  #yの値を[0,1]から[-1,1]にする
  y_exp = [x if x != 0 else -1 for x in y_exp]
  
  #初期値で一度グラフ描画
  plt.ion()

  for index_x1 in range(len(x1_array)):
    for index_x2 in range(len(x2_array)):

      #X1
      inputLayerY[0] = x1_array[index_x1]
      #X2
      inputLayerY[1] = x2_array[index_x2]

      iLayer.foward(inputLayerY)
      oLayer.foward(iLayer.getY())

      val_y[index_x1][index_x2] = oLayer.getY()
  #グラフ表示
  plt.clf()
  plt.subplot(2,1,1)
  plt.xlim(-8,8)
  plt.ylim(-8,8)
  plt.pcolor(pos_x1,pos_x2,val_y)
  plt.colorbar()
  plt.scatter(x_exp[:,0],x_exp[:,1],c=y_exp,s=80)
  plt.subplot(2,1,2)
  plt.plot(cnt_list,match_list)
  plt.pause(0.05)

  
  for i in range(iteration_num):
    match = 0
    #学習フェーズ
    for j in shuffle_list:
      
      #X1
      inputLayerY[0] = x_exp[j][0]
      #X2
      inputLayerY[1] = x_exp[j][1]

      iLayer.foward(inputLayerY)
      oLayer.foward(iLayer.getY())

      if (((y_exp[j] < 0) & (oLayer.getY() < 0)) |
          ((y_exp[j] > 0) & (oLayer.getY() > 0))):
        match += 1

      oLayer.backprop(y_exp[j])
      oLayer.updateWeight()

    #次の学習時に順番をシャッフルする
    random.shuffle(shuffle_list)
    
    match_list.append(match/sample_num)
    cnt_list.append(i)
    #学習結果(ヒートマップ)
    for index_x1 in range(len(x1_array)):
      for index_x2 in range(len(x2_array)):

        #X1
        inputLayerY[0] = x1_array[index_x1]
        #X2
        inputLayerY[1] = x2_array[index_x2]

        iLayer.foward(inputLayerY)
        oLayer.foward(iLayer.getY())

        val_y[index_x2][index_x1] = oLayer.getY()

    #結果グラフ表示
    plt.clf()
    plt.subplot(2,1,1)
    plt.xlim(-8,8)
    plt.ylim(-8,8)
    plt.pcolor(pos_x1,pos_x2,val_y)
    plt.colorbar()
    plt.scatter(x_exp[:,0],x_exp[:,1],c=y_exp,s=80)
    plt.subplot(2,1,2)
    plt.plot(cnt_list,match_list)
    plt.pause(0.05)

  while True:
    plt.pause(0.05)

if __name__ == '__main__': 
  linearlySeparable()

オンライン学習、学習率は0.03固定で実行
実行結果は↓
f:id:ginnyu-tei:20161019194635j:plain:w400

悪いときでも100データ×10回学習すればほぼ100パーセントになった。

円形分類

次は隠れ層を追加する。
f:id:ginnyu-tei:20161019195229g:plain:w400

実装は直線分類とほぼ同じ。
異なるのはパーセプトロンの構成とループ数とテストデータ作成

テストデータ

  #テストデータ作成
  x_exp,y_exp = make_circles(n_samples = sample_num,factor=0.3,noise=0.05)

  x_exp *= 3

  #yの値を[0,1]から[-1,1]にする
  y_exp = [x if x != 0 else -1 for x in y_exp]

3倍しているのはそのままだと円が小さかったので適当に大きくしている。
パーセプトロン構成

  #NN作成
  iLayer  = nn.inputLayer(2,4,4)
  h0Layer = nn.hiddenLayer(2,4,1,0.03,True)
  oLayer  = nn.outputLayer(4,1,1,0.03)

結果は↓
f:id:ginnyu-tei:20161019195458j:plain:w400

100データ×50回学習すればほぼ100パーセントになった。
ただ、初期値がランダムの為か、たまに↓みたいに収束しない時がある。
f:id:ginnyu-tei:20161019195709j:plain:w400

渦巻き分類

もっと隠れ層を追加する。
f:id:ginnyu-tei:20161019200151g:plain:w400

sklearnに渦巻きのデータセットがなさそうだったので数式を実装した。

テストデータ

  x_exp = np.zeros((sample_num,2))
  y_exp = np.zeros(sample_num)
  for i in range (sample_num):
    theta = random.uniform(-3*math.pi,3*math.pi)
    if (theta > 0):
      x_exp[i][0] = 0.5*math.cos(theta)*theta
      x_exp[i][1] = 0.5*math.sin(theta)*theta
      y_exp[i] = 0.5
    else:
      x_exp[i][0] = 0.5*math.cos(theta)*theta
      x_exp[i][1] = 0.5*math.sin(theta)*-theta
      y_exp[i] = -0.5

テストデータは[-1,1]ではなく[-0.5,0.5]にしている。
これは[-1,1]でなかなか収束しない時に↓を見たら、

バックプロパゲーション - Wikipedia

tanhの時は出力を-0.65848~0.6584にすると良いと
書いてあったのでそうすると収束率がかなり改善した。

パーセプトロン構成

  #NN作成
  iLayer  = nn.inputLayer(2,8,8)
  h0Layer = nn.hiddenLayer(2,8,8,0.03,True)
  h1Layer = nn.hiddenLayer(8,8,5,0.03,True)
  h2Layer = nn.hiddenLayer(8,5,1,0.03,True)
  oLayer  = nn.outputLayer(5,1,1,0.03)
  

結果は↓
f:id:ginnyu-tei:20161019201319j:plain:w400

データ数が100だと渦巻きがかなりまばらになったので
これはデータ数を200にしている。
200データ×300回くらいでほぼ100パーセントになる。

長くなったので今回はここまで
次回こそはMINSTやるよ!