FortranからOpenCVを叩く 画像処理関数編

昨日のFortranからOpenCVを叩く 基本データ編　からの続きです。

現在のヘッダファイル

前の記事で紹介したヘッダファイルのFortran版です。以下では、これに必要な関数・変数を追加したり、新しいFortranファイルを追加したりして対応します。

なお、今回は画像データを直接配列に持ってきて弄ることはしないです。

types_c.hのうちの必要な分のFortranインターフェース

types_c.f90

module cv_types_c
    use,intrinsic :: iso_c_binding
    implicit none
    type,bind(C) :: CvMat
        integer(c_int) :: type
        integer(c_int) :: step
        type(c_ptr) :: refcount
        integer(c_int) :: hdr_refcount
        type(c_ptr) :: data
        integer(c_int) :: rows
        integer(c_int) :: cols
    end type
    type,bind(C) :: CvSize
        integer(c_int) :: width
        integer(c_int) :: height
    end type
    type,bind(C) :: CvScalar
        real(c_double) :: val(4)
    end type
    integer(c_int), parameter :: CV_16S = 3
    integer(c_int), parameter :: CV_16U = 2
    integer(c_int), parameter :: CV_32F = 5
    integer(c_int), parameter :: CV_32S = 4
    integer(c_int), parameter :: CV_64F = 6
    integer(c_int), parameter :: CV_8S = 1
    integer(c_int), parameter :: CV_8U = 0
    integer(c_int), parameter :: CV_CN_SHIFT = 3
    integer(c_int), parameter :: CV_DEPTH_MAX = ishft(1, CV_CN_SHIFT)
    integer(c_int), parameter :: CV_DEPTH_MASK = CV_DEPTH_MAX - 1
contains
    function CV_MAKETYPE(depth, cn)
        integer(c_int), intent(in):: depth, cn
        integer(c_int) :: CV_MAKETYPE
        CV_MAKETYPE = iand(depth, CV_DEPTH_MASK) + ishft(cn - 1, CV_CN_SHIFT)
    end function CV_MAKETYPE

end module cv_types_c

core_c.hのうちの必要な分のFortranインターフェース

core_c.f90

module cv_core_c
    use,intrinsic :: iso_c_binding
    use cv_types_c
    implicit none

    interface
        function cvCreateMat(rows, cols, type) bind(C, name="cvCreateMat")
            import :: c_ptr, c_int
            type(c_ptr) :: cvCreateMat
            integer(c_int), value :: rows, cols, type
        end function cvCreateMat
        subroutine cvReleaseMat(mat) bind(C, name="cvReleaseMat")
            import :: c_ptr
            type(c_ptr) :: mat
        end subroutine cvReleaseMat
        subroutine cvSet(arr, val, mask) bind(C, name="cvSet")
            import :: c_ptr, CvScalar
            type(c_ptr),value :: arr, mask
            type(CvScalar),value :: val
        end subroutine cvSet
    end interface
end module cv_core_c

highgui_c.hのうちの必要な分のFortranインターフェース

highgui_c.f90

module cv_highgui_c
    use,intrinsic :: iso_c_binding
    use cv_types_c
    implicit none

    interface
        subroutine cvShowImage(name, image) bind(C, name="cvShowImage")
            import :: c_ptr, CvMat
            character :: name(*)
            type(c_ptr),value :: image
        end subroutine cvShowImage
        function cvWaitKey(delay) bind(C, name="cvWaitKey")
            import :: c_int
            integer(c_int),value :: delay
            integer(c_int) :: cvWaitKey
        end function cvWaitKey
        subroutine cvDestroyAllWindows() bind(C, name="cvDestroyAllWindows")
        end subroutine cvDestroyAllWindows
    end interface

end module cv_highgui_c

描画関数の利用

OpenCVの円を描く関数を使ってみます。

追加の構造体

円の中心など座標のデータのために、OpenCVにはCvPointという構造体があります。

type_c.f90に追加

    type,bind(C) :: CvPoint
        integer(c_int) :: x
        integer(c_int) :: y
    end type

画像処理用モジュール：imgproc

imgproc_c.f90

module cv_imgproc_c
    use,intrinsic :: iso_c_binding
    use cv_types_c
    implicit none
    interface
        subroutine cvCircle(img, center, radius, color, thickness, line_type, shift) bind(C, name="cvCircle")
            import :: c_ptr, c_int, CvPoint, CvScalar
            type(c_ptr), value :: img
            type(CvPoint), value :: center
            type(CvScalar), value :: color
            integer(c_int), value :: radius, thickness, line_type, shift
        end subroutine cvCircle
    end interface
end module cv_imgproc_c

main関数

test3.f90

program test
    use cv_types_c
    use cv_core_c
    use cv_highgui_c
    use cv_imgproc_c
    implicit none 
    type(c_ptr) :: p_img
    type(CvScalar) :: v, color
    type(CvPoint) :: center
    integer(c_int) :: W, H, ret
    W = 320
    H = 240
    p_img = cvCreateMat(H, W, CV_MAKETYPE(CV_64F, 3))
    v%val = 0.d0
    color%val = [0.d0, 1.d0, 0.d0, 0.d0]
    center%x = W/2
    center%y = H/2
    call cvSet(p_img, v, c_null_ptr)
    call cvCircle(p_img, center, 100, color, 10, 1, 0)
    call cvShowImage("hoge", p_img)
    ret = cvWaitKey(0)
    call cvDestroyAllWindows()
    call cvReleaseMat(p_img)
end program test

cvSet(p_img, v, c_null_ptr)はcvSetZero(p_img)でもよいです。

コンパイルと実行

$ gfortran -g highgui_c.f90 types_c.f90 core_c.f90 imgproc_c.f90 test3.f90 -L/usr/local/lib -lopencv_core -lopencv_highgui -lopencv_imgproc
$ ./a.out

狙い通り、黒い背景に緑の丸（幅１０px）を描くことができました。

画像処理関数

モノクロにしてイコライズをかけてみます。画像もせっかくなのでプログラムに作らせたものではなく、lena.jpgを読み込んで見たいと思います。

imgcodecモジュール：画像ファイルIO

module cv_imgcodecs_c
    use,intrinsic :: iso_c_binding
    use cv_types_c
    implicit none
    interface
        function cvLoadImageM( filename, iscolor ) bind(C, name="cvLoadImageM")
            import :: c_ptr, c_int
            character :: filename(*)
            integer(c_int),value :: iscolor
            type(c_ptr) :: cvLoadImageM
        end function cvLoadImageM
    end interface
    integer(c_int), parameter :: CV_LOAD_IMAGE_UNCHANGED = -1
    integer(c_int), parameter :: CV_LOAD_IMAGE_GRAYSCALE = 0
    integer(c_int), parameter :: CV_LOAD_IMAGE_COLOR = 1
    integer(c_int), parameter :: CV_LOAD_IMAGE_ANYDEPTH = 2
    integer(c_int), parameter :: CV_LOAD_IMAGE_ANYCOLOR = 4
    integer(c_int), parameter :: CV_LOAD_IMAGE_IGNORE_ORIENTATION = 128

end module cv_imgcodecs_c

cvLoadImageは画像ファイルを読み込んでIplImageポインタを、後ろにMがついたcvLoadImageMはCvMatポインタを返します。

構造体は初期化しておかないとvoidポインタで渡された先の関数が構造体の種類を識別できなくて困ってしまいますので、nullポインタではなくデータを設定しない構造体を用意します。

imgprocモジュール：色変更とイコライザ

imgproc_c.f90に追加

        subroutine cvCvtColor(src, dst, code) bind(C, name="cvCvtColor")
            import :: c_ptr, c_int
            type(c_ptr), value :: src
            type(c_ptr), value :: dst
            integer(c_int), value :: code
        end subroutine cvCvtColor
        subroutine cvEqualizeHist(src, dst) bind(C, name="cvEqualizeHist")
            import :: c_ptr
            type(c_ptr), value :: src
            type(c_ptr), value :: dst
        end subroutine cvEqualizeHist

これ以外には特に新しいOpenCVの機能は使いません。

main関数

test.f90

program test
    use iso_c_binding
    use cv_types_c
    use cv_highgui_c
    use cv_core_c
    use cv_imgcodecs_c
    use cv_imgproc_c
    implicit none 
    type(c_ptr) :: p_img, p_gimg, p_eqimg
    type(CvMat),pointer :: img
    integer(c_int) :: ret
    character(len=*,kind=c_char), parameter :: filename = "lena.jpg" // char(0)

    p_img = cvLoadImageM(filename, CV_LOAD_IMAGE_COLOR)
    call cvShowImage("hoge", p_img)
    ret = cvWaitKey(0)

    call c_f_pointer(p_img, img)

    p_gimg = cvCreateMat(img%rows, img%cols, CV_MAKETYPE(CV_8U, 1))
    p_eqimg = cvCreateMat(img%rows, img%cols, CV_MAKETYPE(CV_8U, 1))

    call cvCvtColor(p_img, p_gimg, CV_BGR2GRAY)
    call cvShowImage("hoge", p_gimg)
    ret = cvWaitKey(0)
    call cvEqualizeHist(p_gimg, p_eqimg)
    call cvShowImage("hoge", p_eqimg)
    ret = cvWaitKey(0)

    call cvDestroyAllWindows()
    call cvReleaseMat(p_img)
    call cvReleaseMat(p_gimg)
    call cvReleaseMat(p_eqimg)

end program test

コンパイルと実行

gfortran -g highgui_c.f90 types_c.f90 core_c.f90 imgproc_c.f90 imgcodecs_c.f90 test.f90 -L/usr/local/lib -lopencv_core -lopencv_highgui -lopencv_imgcodecs -lopencv_imgproc
./a.out

これを実行すると以下のように3つの画像が順番に見れるはずです。

• 元画像

• グレースケール

• イコライズ

cvLoadImageMのiscolor引数にvalue属性を忘れて悪戦苦闘したことを除いて、とくに難しいところはありません。

物体検出

画像を読み込み、HaarLike特徴量を使ったCascade識別でなにか見つけて、丸で囲って画像を保存したいと思います。

残念なお知らせ

以下の関数はC APIが無い、もしくは使えなかったものです。

• cvLoadHaarClassifierCascade

C++とPythonにおいてはCascade分類器の学習器はXML形式のファイルから読み込むのですが、C APIのcvLoadHaarClassifierCascade(const char * directory,
CvSize orig_window_size)はXMLファイルを読み込むものではありません。ひとつ目の引数の名前がdirectoryであることから分かるように、一定の規則で構成されたフォルダパスを渡す必要がります。
　では、XMLからClassifierに関する構造体を生成するところだけをC++でやるということも考えましたが、どうも前調査の段階ではメモリリークというかいろいろ間違えていたのか、Fortranの方にClassifier構造体を持ってくることができませんでした。
　更に、OpenCVの4.0.0を見てみると、この関数が使えるobjdetectのモジュールのC APIがどっかいきました。なので、Classifier関連の処理（XMLファイルから識別器を読み込み、画像を渡して特徴量の座標を得るまで）はすべてC++にやらせることにします。

• cvSaveImage

これ自体は何の変哲もない、画像をファイルに書き出す関数です。しかし、なぜかC APIから呼び出すとCvMatでもIplImageでも型チェックAssertで失敗します。なので、これもC++でcv::matで読み込んでからCvMatに変換するという関数をC++で書いてしまいます。

C++のラッパー

cpp_wrapper.cpp

#include "opencv2/imgcodecs.hpp"
#include "opencv2/highgui.hpp"
#include "opencv2/imgproc.hpp"
#include "opencv2/objdetect.hpp"
#include <iostream>
#include <vector>

using namespace cv;
using namespace std;

extern "C" {
    CvSeq* cvLoadRunHaarClassifier(const char* filename, CvMat* img);
    int cvImWrite(const char* filename, CvMat* img);
}

CvSeq* cvLoadRunHaarClassifier(const char* filename, CvMat* img){
    CascadeClassifier cascade = CascadeClassifier(String(filename)); 
    vector<Rect> fpoints;
    Mat m = Mat(Size(img->cols, img->rows), img->type, img->data.ptr );
    cascade.detectMultiScale(m, fpoints);

    CvMemStorage* storage = cvCreateMemStorage(0);
    CvSeq* ret = cvCreateSeq(0, sizeof(CvSeq), sizeof(CvRect), storage);
    for(Rect fp: fpoints){
        CvRect cvfp = cvRect(fp.x, fp.y, fp.width, fp.height); 
        int* added = (int*) cvSeqPush(ret, &cvfp);
    }
    cvClearMemStorage(storage);
    cvReleaseMemStorage(&storage);
    return ret;
}

int cvImWrite(const char* filename, CvMat* img){
    Mat m = Mat(Size(img->cols, img->rows), img->type, img->data.ptr );
    imwrite(String(filename), m);

}

FortranとC++で関数名を一致させる必要があるので、FortranのBind(C, name="hogehoge")に相当するextern C{ hogehoge() }を使います。
またC++中ではCvMatもしくはIplImageからcv::Matに変換しますが、このときデータへのポインタのほか、行列の形状の情報（rows, colsなど）が必要になります。このときCvArrを渡していると、voidポインタなので、cv::Matのコンストラクタが対応していない、エディタの補完が効かない、引数に形状の情報が別に必要など不都合がいろいろあります。なので、もういっそCvMatを渡すことにします。

Fortranヘッダー

まずはObjdetectモジュール、というかcvLoadRunHaarClassifierのためのインターフェース。

objdetect_c.f90

module cv_objdetect_c
    use,intrinsic :: iso_c_binding
    use cv_types_c
    implicit none
    interface
        function cvLoadRunHaarClassifier(filename, img) bind(C, name="cvLoadRunHaarClassifier")
            import :: c_ptr, CvMat
            character :: filename(*)
            type(CvMat) :: img
            type(c_ptr) :: cvLoadRunHaarClassifier
        end function
    end interface
end module cv_objdetect_c

付随して、CvSeqに関するデータ型。ついでに、見つけた特徴量をCvRect型に保存していますので、その構造体も用意します。

types_c.f90に追加

    type,bind(C) :: CvRect
        integer(c_int) :: x
        integer(c_int) :: y
        integer(c_int) :: width
        integer(c_int) :: height
    end type

    type,bind(C) :: CvSeq
        integer(c_int) :: flags
        integer(c_int) :: header_size
        type(c_ptr) :: h_prev ! CvSeq*
        type(c_ptr) :: h_next ! CvSeq*
        type(c_ptr) :: v_prev ! CvSeq*
        type(c_ptr) :: v_next ! CvSeq*
        integer(c_int) :: total
        integer(c_int) :: elem_size
        type(c_ptr) :: block_max ! schar*
        type(c_ptr) :: ptr ! schar*
        integer(c_int) :: delta_elems
        type(c_ptr) :: storage ! CvMemStorage*
        type(c_ptr) :: free_blocks ! CvSeqBlock*
        type(c_ptr) :: first ! CvSeqBlock*
    end type

    type,bind(C) :: CvSeqBlock
        type(c_ptr) :: prev ! CvSeqBlock*
        type(c_ptr) :: next ! CvSeqBlock*
        integer(c_int) :: start_index
        integer(c_int) :: count
        type(c_ptr) :: data ! schar*
    end type

cvMemStorageの定義はポインタがメンバーにあるだけなので、Fortranでの定義は要りません。

あとはC++のcvImWriteのインターフェースを用意します。

imgcodecs_c.f90に追加

        function cvImWrite(filename, image) bind(C, name="cvImWrite")
            import :: c_ptr, c_int, CvMat
            character :: filename(*)
            type(CvMat) :: image
            integer(c_int) :: cvImWrite
        end function cvImWrite

CvSeqの使い方

　CvSeqは任意の構造体の配列を管理できます。この構造体には要素数totalがあります。そしてその以下にCvSeqBlockという構造体があります。中身の変数の名前を見ていると、要素1つに1つのCvSeqBlockがあり、それが複数まとまって1つのCvSeqとして管理されているように見えます。
　ということはまずCvSeqのfirstをCvSeqBlockとして取り出し、CvSeqBlock%dataから必要なデータを引っ張ってきて処理、そして次のポインタCvSeqBlock%nextに対して同じ様に行うということをnextがnullになるまでwhileループで繰り返せば良いように思います。

　しかし、やってみるとすべてのCvSeqBlock%dataがCvSeq%first%dataと同じものになってしまいました。
　もしやと思い、CvSeq%first%dataがそのまま配列だと思って、要素数がCvSeq%totalのCvRectの配列として直接Allocateすると、すべてのデータを取り出すことができました。
　なんか使い方を間違っている気がする。

main関数

test.f90

program test
    use cv_types_c
    use cv_highgui_c
    use cv_core_c
    use cv_imgcodecs_c
    use cv_objdetect_c
    use cv_imgproc_c
    implicit none 
    type(c_ptr) :: p_img, storage, p_result
    type(CvSeq),pointer :: results
    type(CvMat),pointer :: img
    integer(c_int) :: ret
    character(len=*), parameter :: filename = "lena.jpg" // c_null_char
    character(len=*), parameter :: cascade_param = "haarcascade_eye.xml" // c_null_char
    type(c_ptr) :: p_seqblock
    type(CvSeqBlock),pointer :: seqblock
    type(CvRect),pointer :: rects(:)
    integer :: i
    type(CvPoint) :: center
    type(CvScalar) :: color
    integer(c_int) :: radius, thickness, line_type, shift
    thickness = 2
    line_type = 1
    shift = 0
    color%val = (/0.0, 255.0, 0.0, 0.0/)

    write(*, *) filename 
    write(*, *) cascade_param
    p_img = cvLoadImageM(filename, CV_LOAD_IMAGE_COLOR)
    call cvShowImage("hoge", p_img)
    ret = cvWaitKey(0)

    call c_f_pointer(p_img, img)
    p_result = cvLoadRunHaarClassifier(cascade_param, img)

    call c_f_pointer(p_result, results)

    write(*, *) "total detected:", results%total
    p_seqblock = results%first
    call c_f_pointer(p_seqblock, seqblock)
    call c_f_pointer(seqblock%data, rects, [results%total])

    do i = 1, results%total
        center%x = rects(i)%x + rects(i)%width/2
        center%y = rects(i)%y + rects(i)%height/2
        radius = (rects(i)%width + rects(i)%height)/4
        call cvCircle(p_img, center, radius, color, thickness, line_type, shift)
    enddo
    call cvShowImage("hoge", p_img)
    ret = cvWaitKey(0)
    ret = cvImWrite("lena_detected.jpg" // C_null_char, img)

    call cvDestroyAllWindows()
    call cvReleaseMat(p_img)

end program test

コンパイルと実行

g++ -std=c++11 -c cpp_wrapper.cpp 
gfortran cpp_wrapper.o highgui_c.f90 types_c.f90 core_c.f90 imgproc_c.f90 imgcodecs_c.f90 objdetect_c.f90 test.f90 -L/usr/local/lib -lopencv_core -lopencv_highgui -lopencv_imgcodecs -lopencv_objdetect -lopencv_imgproc -lstdc++
./a.out

C++で無駄にC++11準拠のfor文を使ったので、-std=c++といったオプションが要ります。また、C++によるプログラムはlibstdc++へのリンクが必要になります。g++でリンクまでやると自動で足してくれるようですが、それ以外のコンパイラだと明示してやらないといけません。

とりあえずlena.jpgから「目」を検出してみます。
目の検出のためのパラメータはhaarcascade_eye.xmlに保存されていますが、このデータはOpenCVのdataディレクトリに入っています。
https://github.com/opencv/opencv/tree/master/data/haarcascades

結果

両目は正しく検出されましたが、鼻の穴が目として誤検出されています。確かにここを拡大すると目のように見えなくもないです。

なんかひどい絵面ですが、これ以上に良いサンプルが見つけられなかった。

今後のOpenCV

OpenCVは2018年末の段階で3.4.3で、4.0.0のpreが出されている状態ですが、C APIについてはメンテナンスを中止するとのことです。
実際、4.0.0のドキュメントを眺めると、3.4.3にはあったobjdetectのC APIが消えています。

さらに、4.0.0-alphaの公式ページには、

The legacy C API from OpenCV 1.x (using CvMat, IplImage, etc.) is partially excluded; the cleanup should mostly be finished by OpenCV 4.0 gold.

とあります。こんだけ苦労して使い方調べたのにもうすぐ使えなくなるというこの感じ。

　今後FortranからOpenCVを活用するには、C++でラッパーを書くか、いっそ別の画像処理用のライブラリを使うのを検討した方がいいかもしれません。たとえば、画像・動画のファイルパスを与えると特徴量の配列を返す関数をC++でOpenCVを使って書き、その結果を使ってFortranで数値計算するという流れが考えられます。

まあ、フロントエンドの言語とライブラリの言語でAPIが用意されてなくて直接呼び出せない時、両方と連携ができる中間層の言語を導入するのはよくあるのかなと。
Python ⇔ C++/CLI ⇔ C# とか。IronPythonが使えなかった。

その他

• ForCVリポジトリ、更新しようかな・・・
• OpenCVを使うならPython、次点でC++がいいです。