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

昨日の[FortranからOpenCVを叩く 基本データ編](https://qiita.com/sage-git/items/4d9fc49e416105097ca0)　からの続きです。

# 現在のヘッダファイル

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

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

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

```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インターフェース

```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インターフェース

```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`という構造体があります。

```fortran:type_c.f90に追加
    type,bind(C) :: CvPoint
        integer(c_int) :: x
        integer(c_int) :: y
    end type
```


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

```fortran: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関数

```fortran: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
```

![Screenshot 2018-11-12 21:45:27.png](https://qiita-image-store.s3.amazonaws.com/0/258019/ff7de843-ed8c-fa86-4ed7-a06b91e02132.png)

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

## 画像処理関数

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

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

```fortran
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モジュール：色変更とイコライザ

```fortran: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関数

```fortran: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
```

### コンパイルと実行

```bash
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つの画像が順番に見れるはずです。

* 元画像

<img src="https://qiita-image-store.s3.amazonaws.com/0/258019/33b0836e-0240-3ef8-2572-6e7be9ea677f.jpeg" width="200px">

* グレースケール

<img src="https://qiita-image-store.s3.amazonaws.com/0/258019/0048af40-412b-73f1-9f65-b92fdf245985.png" width="200px">

* イコライズ

<img src="https://qiita-image-store.s3.amazonaws.com/0/258019/689cba78-b107-3cc9-b7ae-961f9ae4b282.png" width="200px">

`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: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`のためのインターフェース。

```fortran: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`型に保存していますので、その構造体も用意します。

```fortran: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`のインターフェースを用意します。

```fortran: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関数

```fortran: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
```

### コンパイルと実行

```bash
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

### 結果

<img src="https://qiita-image-store.s3.amazonaws.com/0/258019/07ca162b-6758-af96-7cc0-c36a7ace1acc.jpeg" width="200px">

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

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

# 今後のOpenCV

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

さらに、[4.0.0-alphaの公式ページ](https://www.opencv.org/opencv-4-0-0-alpha.html)には、

> 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++がいいです。