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にインストール)
参考資料
ここのページで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を見てみます。これはメモリ上に領域を確保し、ポインタを返す関数です。
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での配列データのポインタを渡せばいいんだな」と察しがつくようにしています。
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*によるポリモーフィズムが再現できる。
- Fortranで
-
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
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にはメモリの割当等、画像を取扱う上で必須になる関数群があります。
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関数
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モジュール
とりあえず、画像全体に同じ画素値を設定する関数です。
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です。
type,bind(C) :: CvScalar
real(c_double) :: val(4)
end type
この構造体１つで、１つの画素のRGBAを保持できます。なかのメンバ変数はdouble型ですが、例えばucharの画像に使った時は適当にキャストされるようです。使わないvalの要素（たとえばRGB画像に対してAチャンネルなど）は無視されるようです。
画像描写の処理：highguiモジュール
QtもしくはGTKのGUIを使い、画像を表示したりキー入力を受け付けたりするモジュールです。
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を渡せば無制限に待ちます1。
cvWaitKeyの戻り値で押したキーのコードが得られますので、簡単なアプリもできます。
cvDestroyAllWindows()で生成したWindowをすべて閉じます。
Main関数
画像配列メモリ割り当て、データ代入、表示、データ開放。
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関数
画像配列用メモリ割り当て、円の範囲に代入、縦長の範囲に代入、表示、削除。
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)
というのがあります。基本的に式の中でも利用できるこれを使っていきます。
-
- 演算のたびに変換する
where(iand(255, imgdata(1, :, :)) > 100) imgdata(2, :, :) = -127
-
- 普通の４バイト整数に一旦変換して、そっちを操作する
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関数
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を叩く 画像処理関数編」に続きます。
-
なお、
cvWaitKey(0)で待っている時にwindowをマウスで[x]を押して閉じてしまうと、Ctrl+Cするしかなくなります。 ↩