FortranからOpenCVを叩く 基本データ編

Fortranは2003になってISO_C_BINDINGモジュールが追加されました。これにより、C言語ライブラリとの連携がよりやりやすくなりました。
他にも関数ポインタとかもありますが、これは別の機会に。あるいは他の人が書いているのでは。

CのヘッダをFortranに翻訳する流れにさらっと触れて、OpenCVの基本的な関数を試し、最終的にはHaar-Like特徴量によるCascade分類器を試したいと思います。

結構長くなりましたので、記事を２つに分けました。
今日の内容は

• CのヘッダをFortranに翻訳
• OpenCVの画像データ型の取り扱い

です。
明日はその続きの

• OpenCVの基本的な画像処理機能
• 物体検出

としたいと思います。

検証環境

• Ubuntu 16.04
• gfortran 5.4.0 20160609
• OpenCV 3.4.3 (ソースからビルドで/usr/localにインストール)

参考資料

OpenCVドキュメント

ここのページでMain Modules以下にある適当なモジュールを選んでリンクをたどり、その先でC structures and operations とか C API といったところをたどると、C言語での実装の参考資料が出てきます。
また、インストールされたライブラリにおいて、たとえばinclude/opencv2/core/core_c.hやtypes_c.hといったように、_c.hで終わっているヘッダはC言語用ですので、これを見れば結構参考になります。

正直これまでOpenCVを使うときはPythonからのnumpy.ndarrayとか、C++でのcv::Matとかなので、低級言語でのAPIはあまり馴染みがないです。

実装例の前例

手前味噌ですが、昔Githubにモジュールを公開しました。

なお、CのヘッダファイルをFortranのモジュールに自動変換する何かがあれば、このプロジェクトは意味がなくなります。せいぜいFortranからみたときの機能の取捨選択・整理をすれば自動変換以上の価値が生まれるかな？（明日の記事でふれますが自動変換とは別の理由で意味がなくなりそう）
あと、このforcvリポジトリではIplImageベースの管理です。今回はCvMatで行いたいと思います。

CのヘッダをFortranのインターフェースへ

大雑把にいうとCとFortranは同じライブラリを使えます。あとはそのCによる表現であるところのヘッダをFortranのインターフェースに書き直してやるだけです。
ライブラリから出てきたデータをFortranでどう取り扱うかは後ほど。

構造体

たとえば画像のサイズを保持する構造体struct CvSizeはtypes_c.hの中で以下のように書かれています。

typedef struct CvSize
{
    int width;
    int height;

#ifdef __cplusplus
    CvSize(int w = 0, int h = 0): width(w), height(h) {}
    template<typename _Tp>
    CvSize(const cv::Size_<_Tp>& sz): width(cv::saturate_cast<int>(sz.width)), height(
cv::saturate_cast<int>(sz.height)) {}
    template<typename _Tp>
    operator cv::Size_<_Tp>() const { return cv::Size_<_Tp>(cv::saturate_cast<_Tp>(wid
th), cv::saturate_cast<_Tp>(height)); }
#endif
}
CvSize;

すごいややこしく見えますが、OpenCVのCヘッダは場合によってはC++などからも使われますし、MicrosoftのVCコンパイラにおいては特殊な修飾語が付いたりします。なので、OpenCV（というか様々な環境で使われること前提のライブラリ）というのはヘッダが冗長になりがちですが、今回のFortranではGCCで見られるであろう部分しか使いません。なので、上記の構造体は

    type,bind(C) :: CvSize
      integer(c_int) :: width
      integer(c_int) :: height
    end type

で終わりです。

関数

まずはcvCreateMatを見てみます。これはメモリ上に領域を確保し、ポインタを返す関数です。

include/opencv2/core/core_c.h

CVAPI(CvMat*)  cvCreateMat( int rows, int cols, int type );

このintは値渡しなので、Fortranでinterfaceを書くときはinteger(c_int), valueとvalueを指定します。
一方、戻り値のCVAPI(CvMat*)ですが、これはgccだと結局CvMat*になります。CVAPI()はC++やVisualStudioでビルドした時に必要に応じて修飾語をつけるためにあります。
これはFortranから見るとただのポインタなので、type(c_ptr)で受け取ります。

    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
    end interface

引数に整数を3つ渡していますが、値渡しなのでFortranではvalueをつけてやります。

一方でcvSetZeroを見てみます。これは行列をすべて０で埋めるサブルーチンです。

CVAPI(void)  cvSetZero( CvArr* arr );

void戻りなので、Fortranではsubroutineで定義します。
一方引数のCvArr*ですが、これはC言語なりに実現したポリモーフィズムです。というのは、　OpenCVにおいて、Cで画像データを保持する構造体として、CvMatとIplImageがあります。一方で画像処理の関数の多くは別にどっちの構造体でもよい話です。構造体が複数種類あるからと言って関数も分けていたら、ただでさえ多いOpenCVの関数群が2倍3倍と膨れ上がってしまうので、ポインタだけ受け取って構造体のフラグから内部で細かな処理を分けるということをしています。C言語なりに実装したオブジェクト指向ですね。
　ただ、この時に使っているポインタはvoid*で、このままだとヘッダを見た人もこれに何を渡しているか分からなくなります。そこでOpenCVではCvArrをvoidで定義して、CvArr*と書くことで「これはOpenCVでの配列データのポインタを渡せばいいんだな」と察しがつくようにしています。

include/opencv2/core/types_c.h

typedef void CvArr;

さて、Fortranでのインターフェースですが、CvMatのポインタ渡しとして

type(c_ptr), value :: arr

type(CvMat) :: arr

の2通りが考えられます。それぞれ以下のような欠点・利点があります。

• type(c_ptr), value :: arr : ポインタの値渡し。

  • FortranでCvMat構造体をローカル変数として管理していた場合、この関数とのやり取りにおいてはc_f_pointerによるポインタ割当てや、c_locによるCポインタの取得が必要。
  • 逆にFortranはポインタの受け渡しだけに徹していた場合、CポインタとFortran構造体の間で変換する必要がない。
  • type(c_ptr)には構造体の種類の情報はないので、void*によるポリモーフィズムが再現できる。

• type(CvMat) :: arr : 値の参照渡し。

  • FortranでCポインタの取得が必要ない。
  • 逆に、Cポインタで得られたデータを渡すときは、c_f_pointerでFortran内の変数に割り当てる必要がある。
  • CではIplImageとCvMatの両方を受け取れていたものが、Fortranからの呼び出しではtype(CvMat)しか受け取れない。

以上を考えるとポインタの値渡しを採用して、cvSetZeroは

    interface
        subroutine cvSetZero(arr) bind(C, name="cvSetZero")
            import :: c_ptr
            type(c_ptr),value :: arr
        end subroutine cvSetZero
    end interface

というように書くことができます。

type(c_ptr),valueはやっぱりFortranの参照渡しが使いこなせていないように見えますが、void型ポインタをやり取りするライブラリだとこうするしか無いのかなと。

実例

データ置き場の生成

まず画像に触れる前に、コンパイルやメモリ割り当てがちゃんとできることを確かめます。

基本的な構造体：coreモジュールのtypes_c.h

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

これらはinclude/opencv2/core/types_c.hに定義されている、OpenCV中で使われている基本的なデータ型とか定数です。
CV_MAKETYPEはsubroutineにせず、結果をCV_MAKETYPE(CV_8U, 3)の結果をCV_8UC3として定数に定義してもいいのですが、今回はサブルーチンで都度計算させました。どうせ最初に1回だけしか呼び出さないし。

行列生成サブルーチン：coreモジュールのcore_c.h

core_c.hにはメモリの割当等、画像を取扱う上で必須になる関数群があります。

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
    end interface
end module cv_core_c

これらはOpenCVを使うならまず間違いなく使われている関数群です。

main関数

test.f90

program test1
    use cv_types_c
    use cv_core_c
    implicit none 
    type(c_ptr) :: p_img
    type(CvMat),pointer :: img
    integer(c_int) :: W, H
    W = 640
    H = 480
    p_img = cvCreateMat(H, W, CV_MAKETYPE(CV_64F, 1))
    call c_f_pointer(p_img, img)
    write(*, *) img%cols, img%rows
    write(*, *) img%data
    write(*, *) img%type, img%step
    call cvReleaseMat(p_img)
end program test1

コンパイルと実行

$ gfortran types_c.f90 core_c.f90 test1.f90 -L/usr/local/lib -lopencv_core

OpenCVのリンクオプションは長くなりがちなので、CMakeなどに任せるのが吉です。まあ今回は手で書いてしまいます。

$ ./a.out                                                                    
         640         480
      140497251598464
  1111638022        5120

ちゃんと構造体の適切なところに適切なデータが入ったようですね。stepの5120は８バイトのデータが１行640個あるものだと思われます。

画素値の設定と表示

わざわざ画像をこういうプログラムで操作する目的は画素値の操作です。あと編集した後すぐに確認したいところ。

必要な関数を定義します。

画素値の設定：coreモジュール

とりあえず、画像全体に同じ画素値を設定する関数です。

core_c.f90のinterfaceに追加

        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

ここで１ピクセルあたりの画素値はスカラー値とは限らず、RGBの３次元だったりRGBAの４次元だったりします。これらを統一的に扱う構造体がCvScalarです。

types_c.f90に追加

    type,bind(C) :: CvScalar
        real(c_double) :: val(4)
    end type

この構造体１つで、１つの画素のRGBAを保持できます。なかのメンバ変数はdouble型ですが、例えばucharの画像に使った時は適当にキャストされるようです。使わないvalの要素（たとえばRGB画像に対してAチャンネルなど）は無視されるようです。

画像描写の処理：highguiモジュール

QtもしくはGTKのGUIを使い、画像を表示したりキー入力を受け付けたりするモジュールです。

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

cvShowImageで画像を表示するウィンドウが生成されます。このままだとプログラムが先に進んでしまいますが、cvWaitKeyでプログラムを指定ミリ秒数経過するか、何かキーが押されるまで停止します。0を渡せば無制限に待ちます¹。
cvWaitKeyの戻り値で押したキーのコードが得られますので、簡単なアプリもできます。
cvDestroyAllWindows()で生成したWindowをすべて閉じます。

Main関数

画像配列メモリ割り当て、データ代入、表示、データ開放。

test2.f90

program test1
    use cv_types_c
    use cv_highgui_c
    use cv_core_c
    implicit none 
    type(c_ptr) :: p_img
    type(CvMat),pointer :: img
    type(CvScalar) :: v
    integer(c_int) :: W, H, ret
    W = 320
    H = 240
    p_img = cvCreateMat(H, W, CV_MAKETYPE(CV_64F, 3))
    call c_f_pointer(p_img, img)
    v%val = 0.d0
    v%val(1) = 1.d0
    call cvSet(p_img, v, c_null_ptr)
    call cvShowImage("hoge", p_img)
    ret = cvWaitKey(0)
    call cvDestroyAllWindows()
    call cvReleaseMat(p_img)
end program test1

コンパイルと実行

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

すると下図のように青１色の画像ができるはずです。

画像データといえば1byte整数というイメージですが、OpenCVでは画素値として浮動小数点も扱えます。このとき、0.d0が0に、1.d0が255に対応し、この範囲外の値はこの範囲にクリッピングされます。

データ直接編集　浮動小数点

画素値を直接いじってみます。画像は結局配列なので、配列が得意なFortranはさぞかし相性が良いことでしょう。

今回はこれまでに紹介した関数だけで大丈夫です。まずはFortranが得意な倍精度浮動小数点の行列を扱います。

main関数

画像配列用メモリ割り当て、円の範囲に代入、縦長の範囲に代入、表示、削除。

test4.f90

program test
    use cv_types_c
    use cv_highgui_c
    use cv_core_c
    implicit none 
    type(c_ptr) :: p_img
    type(CvMat),pointer :: img
    integer(c_int) :: W, H, ret
    double precision, pointer :: imgdata(:, :, :)
    integer :: i, j
    integer,allocatable :: mask(:, :)
    
    W = 80
    H = 120
    p_img = cvCreateMat(H, W, CV_MAKETYPE(CV_64F, 3))
    call c_f_pointer(p_img, img)
    call c_f_pointer(img%data, imgdata, (/3, W, H/))
    imgdata = 0.d0
    allocate(mask(W, H))
    forall(i = 1:W, j = 1:H)mask(i, j) = (i - 40)**2 + (j - 40)**2
    where(mask < 900) imgdata(3, :, :) = 1.d0
    imgdata(1, 20:40, :) = 1.d0
    call cvShowImage("hoge", p_img)
    ret = cvWaitKey(0)
    call cvDestroyAllWindows()
    call cvReleaseMat(p_img)

end program test

Cポインタのp_imgからFortranポインタへのimg、そしてその中にある配列データへのCポインタdataからFortran配列ポインタへのimgdataへの２回のポインタ変換が必要になります。
また、配列の添字はC言語やnumpyだとX[H][W][ch]という順番ですが、Fortranでは逆にX(ch, W, H)という順序になります。
maskは、添字の数字ベースで中心座標(40, 40)からの距離を入れて、whereで画素を代入するピクセルをマスクするために導入しました。もっといいやり方はあるんじゃないかな、と思います。
とりあえず、 where(mask < 900) imgdata(3, :, :) = 1.d0で半径30の円を赤色で描写し、その上から imgdata(1, 20:40, :) = 1.d0で縦長の帯を青色で描写します。

コンパイルと実行

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

狙い通りです。重なったところは紫色（ピンク色？）になりました。

データ直接編集　符号なし８ビット整数

今度は画素が８ビット整数の画像を取り扱ってみます。
まず、C言語で符号なし８ビット整数はucharで扱うのがセオリーです。同様のことをやってみます。

!前略
    character(len=1), pointer :: imgdata(:, :, :)
!中略
    imgdata(1, 20:40, :) = char(255)
!後略

これは代入はできます。ですが、大小比較や算術演算ができません。

では、１バイト整数を使ってみます。

  integer(1), pointer :: imgdata(:, :, :)
  imgdata(2, :, :) = 255

これはコンパイル時に次のようなエラーになります。

!         imgdata(2, :, :) = 255
!                       1
!Error: Arithmetic overflow converting INTEGER(4) to INTEGER(1) at (1). This check can be disabled with the option ‘-fno-range-check’

なので、次のようにする必要があります。

  integer(1), pointer :: imgdata(:, :, :)
  imgdata(2, :, :) = -127

これでとりあえず255の画素で設定できました。ですが、この方法では問題があります。

次のような画像を用意します。

forall(i = 1:W) imgdata(1, i, :) = i - 1

　この画像は左端が(0, 0, 0)で右端が(255, 0, 0)の横方向グラデーションとなっています。この画像で、青の画素値が１００より大きいピクセルについて緑色を重ねたいと思います。その範囲はx軸で100 ~ 255なので、真ん中より若干左側から右端にかけて、半分以上の部分で色が変わるはずです。

この操作を直感的に書くと次のようになります。

    where(imgdata(1, :, :) > 100) imgdata(2, :, :) = -127

ですが、この結果は次の図のようになります。

ちょうど真ん中で画素値が127から128, 129・・・となるのですが、これは符号あり整数だと
127, -127, -126・・・となってしまっています。このため、大小関係が正しく評価できず、画素の操作が思惑通りにいきませんでした。
この問題の回避のため、Fortranが0から255を正の整数として扱える4バイト整数に変換します。やり方は２通り考えられます。
ちなみにucharのデータをそのまま入れたinteger*1を、オーバーフローせずにintegerなどに変換する方法として、

（２バイト以上の整数） = iand(255, uchar)

というのがあります。基本的に式の中でも利用できるこれを使っていきます。

• 1. 演算のたびに変換する

    where(iand(255, imgdata(1, :, :)) > 100) imgdata(2, :, :) = -127

• 1. 普通の４バイト整数に一旦変換して、そっちを操作する

    integer, allocatable :: imgdata_f(:, :, :)
!中略
    allocate(imgdata_f(3, W, H))
    imgdata_f = iand(255, imgdata)
    where(imgdata_f(1, :, :) > 100) imgdata_f(2, :, :) = 255
    imgdata = iand(255, imgdata_f)

これでやると、狙い通りの画像編集ができました。

個人的に直感的な表記に近い２）の方法をとりたいと思います。

結論のmain関数

test7.f90

program test
    use cv_types_c
    use cv_highgui_c
    use cv_core_c
    implicit none 
    type(c_ptr) :: p_img
    type(CvMat),pointer :: img
    type(CvScalar) :: v
    integer(c_int) :: W, H, ret
    integer(1), pointer :: imgdata(:, :, :)
    integer :: i, j
    integer, allocatable :: imgdata_f(:, :, :)
    W = 256
    H = 120
    p_img = cvCreateMat(H, W, CV_MAKETYPE(CV_8U, 3))
    call c_f_pointer(p_img, img)
    call cvSetZero(p_img)
    call c_f_pointer(img%data, imgdata, (/3, W, H/))
    forall(i = 1:W) imgdata(1, i, :) = i - 1
    allocate(imgdata_f(3, W, H))
    imgdata_f = iand(255, imgdata)
    where(imgdata_f(1, :, :) > 100) imgdata_f(2, :, :) = 255
    imgdata = iand(255, imgdata_f)

    call cvShowImage("hoge", p_img)
    ret = cvWaitKey(0)
    call cvDestroyAllWindows()
    call cvReleaseMat(p_img)
end program test

コンパイルと実行

$ gfortran highgui_c.f90 types_c.f90 core_c.f90 test7.f90 -L/usr/local/lib -lopencv_core -lopencv_highgui
$ ./a.out

なぜ今日の記事が長くなったのか

普通、画像は8bitの符号なし整数で1つの画素・チャンネルの強度を表します。C言語などではuchar、numpyではuint8が該当します。ところが、Fotranのintegerは符号ありです。そもそもFortranの仕様で符号なし整数はありません。

配列のポインタを画像処理関数に渡す関係上、変数サイズは揃えたいので、integer(2)などの配列を割り当てて代用するのは避けたいです。
またこういうプログラム上で画像を取り扱うときは、チャンネル強度を直接使った操作もあります。（たとえば2で割って暗くする、など）。そういうわけで、数字としての操作ができないcharacter(len=1)は使えないですし、符号付き整数のinteger*1は先頭ビットが1のとき（負のとき）の変換がiandを使うぶんだけややこしくなりますので、使いたくないです。

コンパイラ拡張だったり内部表現だったりで等価なデータ表現を実現できるかもしれませんが、環境依存が強そうなので避けた方が無難かと思います。コーディングも面倒な事になりそうだし。

幸いOpenCVの関数のほとんどは画素を浮動小数点で表したfloat型のCvMatでも動作しますので、Fortranでの処理は基本的に浮動小数点で読み込んでdouble precisionで、もしくはucharで読み込んでinteger*4に変換して行い、ここぞというところでucharに変換するべきでしょう。

感想

とはいえ、画像配列への直接の代入でループ構文を書かなくてよいので、データ操作の部分はすっきりしているかなと思います。そこはまあ狙い通りでした。

続き

「FortranからOpenCVを叩く 画像処理関数編」に続きます。

1. なお、cvWaitKey(0)で待っている時にwindowをマウスで[x]を押して閉じてしまうと、Ctrl+Cするしかなくなります。 ↩