HSPで画像処理する際の覚書

# 概要
　HSP(Hot Soup Processor)は手軽に使えるスクリプト言語ですが、「処理速度が遅い」などの[致命的な欠陥](http://qiita.com/YSRKEN/items/0c059abfb8685051cb16)があります。
　とはいえ、ちょっとしたGUIツールを作る際には手軽で便利なせいか、この2016年でも未だに生き残っています。
　一方、画像処理には通常、C言語などの高速な言語が使われることが多いです。一見相矛盾する存在ですが、なまじ**[HSPではOpenCVが使える](http://www.onionsoft.net/hsp/ref/docs/hspcv.txt)**ため、HSPだけでも案外回せてしまうのです。
　そこで今回の記事では、HSPで画像処理をしたい人向けのtips集を書いてみました。

※当記事は[Hot Soup Processor Advent Calendar 2016](http://qiita.com/advent-calendar/2016/hsp)の10日目です
※文中、時間計測には自作PCのWindows 10(64bit)＋Core i7-4790K(定格)を使用しました

# 目次
- 画像処理命令を概観する
- VRAMを直接操作してみる
- 画像間のPSNRを取得してみる
- 画像のハッシュを取得してみる
- 画像拡大/縮小アルゴリズムの質について
- 画像をフーリエ変換してみる

# 画像処理命令を概観する
　そもそもHSPでは、「ウィンドウID」によって番号付けられた「画面」が存在します。
　各画面には、ウィンドウとして**表示されるもの**(screen・bgscr命令で作成)と**表示されないもの**(buffer命令で作成)が存在します。以下の説明では便宜上、前者を「実画面」、後者を「仮想画面」、両者を合わせて「標準画面」と称します。
　ヘルプファイルには「screenで初期化されたウィンドウIDをbuffer命令で再初期化することはできません」とありますが、試したところ**buffer命令で初期化したウィンドウIDもscreen命令で使いまわせません**のでご注意ください。つまり、ウィンドウIDは、実画面同士と、仮想画面同士との間でのみ使いませます。
　そして初期化した画面に対しては、lineやcircleやboxfなどの描画命令で絵を描くことができます。絵を描く対象となるウィンドウIDは**描画先**、描く際の起点となる座標を**カレントポジション**と呼びます。この2つはどちらも(暗黙の)グローバル変数ですので、意識してないと忘れがちですから注意が必要です。
　実画面に描画する際、redraw命令を使うと「**一通り描き終わってから表示を更新する**」ことができます。これを行わないと、1操作毎に画面に結果が反映されるので処理が遅くなります。ちなみに、「gsel命令などで非表示にした実画面」と「仮想画面」に対してはredraw命令を使う必要はありません。
　画像処理で最もよく使うのはgcopy命令でしょう。要するに画像を切り貼りする命令ですが、gmode命令によって「透明色付きコピー」や「ピクセルアルファブレンドコピー」などを行わせることができます。切り貼りと同時に画像を伸縮させたい場合はgzoom命令を使いますが、こちらはgmode命令による細かな制御ができないといった問題があります。

　上記(基本命令)だけでも色々できますが、HSPには標準で強力な画像処理ライブラリが複数付属しています。
　まず[Artlet2D](http://sprocket.babyblue.jp/html/hsp_artlet2d.htm)。GDI+の機能を利用していますので、**アルファ値付きの画像を正常に扱うことができる**(標準画面は24bitRGB)上、「ブラシ」を別途定義することで、幅を持った線やグラデーション入りの線などを簡単に表現することができます。また、Artlet2D用に作った仮想イメージは標準画面に書き出したり、逆に読み込んだりすることが可能です。
　次に[HSPCV](http://www.onionsoft.net/hsp/ref/docs/hspcv.txt)。HSPでOpenCVの機能の一部を使うことができるプラグインですが、一部とはいえその威力はなかなかのもの。「様々なアルゴリズムで拡大縮小・回転・ぼかし・二値化・**パターンマッチング**」「**顔認識**」「**Webカメラなどの画像を取得・保存**」などの機能を使うことができます。こちらにも独自のバッファ構造(CVバッファ)が存在しますが、標準画面と読み書き可能です。
　以下、HSPで特殊な画像処理をする際のtipsを書いてみますが、速度面の問題から、**先にこれらのライブラリで実行できないか**をよく勘案しましょう。場合によっては、HSPでこねくり回さず、**画像全体のデータをC言語などで作成したDLLに投げた方が早い**こともありますので……。

# VRAMを直接操作してみる
　標準画面の色を取得/色を書き込みするには、基本的にはpget/pset命令を使います。
　ただ、大量に取得する場合など、それだけだと処理が遅く感じる場合はよくあります。
　そこでよく使われるのが、標準画面の「中身」とも言えるVRAMデータを直接叩く手法です。
　VRAMデータはタダのバイナリ列ですので、そこをpeek/poke命令で読み書きできます。
　ただし、VRAMデータは直感に反し、「BGRBGRBGR...」の順で書き込まれます。
　また、書き込み方向も「左下→右下→左上→右上」といった風になります。
(※WindowsのBMPデータも左下スタートなので、24bitBMPならそのままmemcpy命令等でコピペできる)
　なお、VRAMを書き換えた結果を実画面に反映(表示)する際は、redraw命令を挟む必要があります。
　VRAM操作でpget/psetを代替するとこんな感じ↓。

```vram.hsp
#module vram
	#deffunc pgetv int vx_, int vy_
		// 1行毎に横幅が4の倍数バイトになるよう調整されているので注意
		// 例えばginfo_winx=98なら、vwidth=296バイトが横幅のバイト数である
		vwidth = ((ginfo_winx * 3) + 3) & 0xfffffffc
		// 上から下ではなく、下から上にデータが入っているので注意
		vpos = ((ginfo_winy - 1 - vy_) * vwidth) + (vx_ * 3)
		// データはBGRの順で入っている。66はVRAMデータを示す定数だが、
		// なぜ定数定義が標準で存在しないのか理解に苦しむ
		mref vram, 66
		vr = peek(vram, vpos + 2)
		vg = peek(vram, vpos + 1)
		vb = peek(vram, vpos)
		color vr, vg, vb
	return
	#deffunc psetv int vx_, int vy_
		vwidth = ((ginfo_winx * 3) + 3) & 0xfffffffc
		vpos = ((ginfo_winy - 1 - vy_) * vwidth) + (vx_ * 3)
		mref vram, 66
		poke vram, vpos + 2, ginfo_r
		poke vram, vpos + 1, ginfo_g
		poke vram, vpos, ginfo_b
	return
#global
```

　ちなみにベンチすると、次のような結果になりました。
(表中の※は、画面サイズが変わらないと変化しない変数vwidthを事前計算した場合)

|手法|時間[ms]|速度比(標準命令を1とする)|
|----|--------|-------------------------|
|標準命令-pget|4443|1|
|VRAM-pget|3756|1.18|
|VRAM-pget(※)|3302|1.35|
|標準命令-pset|1219|1|
|VRAM-pset|2321|0.52|
|VRAM-pset(※)|1941|0.63|

```vram2.hsp
#uselib "winmm.dll"
#cfunc global timeGetTime "timeGetTime"
#func global timeBeginPeriod "timeBeginPeriod" int
#func global timeEndPeriod "timeEndPeriod" int

#module vram
// 中略
#global

loop_max = 100:x_size = 200:y_size = 200:buffer 1, 640, 480

timeBeginPeriod 1

// 空ループの速度を計測
start_time = timeGetTime()
for k, 0, loop_max:for y, 0, x_size:for x, 0, y_size
next:next:next
diff_time0 = timeGetTime() - start_time

// 標準命令-pget
start_time = timeGetTime()
for k, 0, loop_max:for y, 0, x_size:for x, 0, y_size
	pget x, y
next:next:next
diff_time1 = timeGetTime() - start_time - diff_time0

// VRAM-pget
start_time = timeGetTime()
for k, 0, loop_max:for y, 0, x_size:for x, 0, y_size
	pgetv x, y
next:next:next
diff_time2 = timeGetTime() - start_time - diff_time0

// 標準命令-pset
start_time = timeGetTime()
for k, 0, loop_max:for y, 0, x_size:for x, 0, y_size
	pset x, y
next:next:next
diff_time3 = timeGetTime() - start_time - diff_time0

// VRAM-pset
start_time = timeGetTime()
for k, 0, loop_max:for y, 0, x_size:for x, 0, y_size
	psetv x, y
next:next:next
diff_time4 = timeGetTime() - start_time - diff_time0

timeEndPeriod 1

gsel 0
mes diff_time1:mes diff_time2:mes diff_time3:mes diff_time4
```

# 画像間のPSNRを取得してみる
　普通に考えれば、VRAMを直接叩いて数値を読み出し、ゴリゴリ演算するのが最速なはずです。
　腹立たしいことに、HSPでは**for文よりrepeat文の方が1.88倍ほど速い**のでそちらを使って測ってみると、次のような結果になりました(単位は[ms])。

|画像サイズ|peekで読む部分|差の二乗を計算して足す部分|合計|
|----------|--------------|--------------------------|----|
|640x480|160|60|220|
|1920x1080|400|140|540|

```psnr1.hsp
#include "hspmath.as"
#uselib "winmm.dll"
#cfunc global timeGetTime "timeGetTime"
#func global timeBeginPeriod "timeBeginPeriod" int
#func global timeEndPeriod "timeEndPeriod" int

// 通常
buffer 1 :picload "a.png"
buffer 2 :picload "b.png"
timeBeginPeriod 1 :start_time = timeGetTime()
	px = ginfo_winx :py = ginfo_winy
	vwidth = ((px * 3) + 3) & 0xfffffffc
	gsel 1:mref vram1, 66
	gsel 2:mref vram2, 66
	sum_r = 0.0 :sum_g = 0.0 :sum_b = 0.0
	repeat py
		// 比較する座標の順番はPSNRに関係しない
		vpos_b = cnt * vwidth
		vpos_g = vpos_b + 1
		vpos_r = vpos_b + 2
		repeat px
			diff_b = peek(vram1, vpos_b) - peek(vram2, vpos_b)
			diff_g = peek(vram1, vpos_g) - peek(vram2, vpos_g)
			diff_r = peek(vram1, vpos_r) - peek(vram2, vpos_r)
			sum_b += diff_b * diff_b
			sum_g += diff_g * diff_g
			sum_r += diff_r * diff_r
			vpos_b += 3
			vpos_g += 3
			vpos_r += 3
		loop
	loop
	if(sum_r != 0.0) :psnr_r = str(-logf(255.0 / sqrt(sum_r)) * 20 / M_LN10) :else :psnr_r = "∞"
	if(sum_g != 0.0) :psnr_g = str(-logf(255.0 / sqrt(sum_g)) * 20 / M_LN10) :else :psnr_g = "∞"
	if(sum_b != 0.0) :psnr_b = str(-logf(255.0 / sqrt(sum_b)) * 20 / M_LN10) :else :psnr_b = "∞"
diff_time = timeGetTime() - start_time :timeEndPeriod 1
gsel 0
mes diff_time
mes psnr_r
mes psnr_g
mes psnr_b
```

　ただ実は、**RGB値に対するPSNR**に限った場合、HSPCVを使えば**高速化することが可能**です。
　ポイントとしては、HSPCVにはCVバッファ間でコピーするcvcopy命令があり、演算オプションに**差の絶対値を返すCVCOPY_DIFがある**ということ。つまり、一旦2枚の画像をCVバッファにコピーし、差の絶対値を取った後に標準画面に書き出し、その後VRAMを叩けばいいのです。

```psnr2.hsp
#include "hspcv.as"
#include "hspmath.as"
#uselib "winmm.dll"
#cfunc global timeGetTime "timeGetTime"
#func global timeBeginPeriod "timeBeginPeriod" int
#func global timeEndPeriod "timeEndPeriod" int

// OpenCV利用
buffer 1 :picload "a.png"
buffer 2 :picload "b.png"
timeBeginPeriod 1 :start_time = timeGetTime()
	gsel 1
	px = ginfo_winx :py = ginfo_winy
	cvbuffer 1, px, py :gsel 1 :cvputimg 1
	cvbuffer 2, px, py :gsel 2 :cvputimg 2
	cvcopy 1,0,0,2,CVCOPY_DIF
	buffer 3, px, py :cvgetimg 2, 0 :mref vram, 66
	sum_r = 0.0 :sum_g = 0.0 :sum_b = 0.0
	vwidth = ((px * 3) + 3) & 0xfffffffc
	repeat py
		// 比較する順番はPSNRに関係しない
		vpos_b = cnt * vwidth
		vpos_g = vpos_b + 1
		vpos_r = vpos_b + 2
		repeat px
			diff_b = peek(vram, vpos_b)
			diff_g = peek(vram, vpos_g)
			diff_r = peek(vram, vpos_r)
			sum_b += diff_b * diff_b
			sum_g += diff_g * diff_g
			sum_r += diff_r * diff_r
			vpos_b += 3
			vpos_g += 3
			vpos_r += 3
		loop
	loop
	if(sum_r != 0.0) :psnr_r = str(logf(255.0 / sqrt(sum_r/px/py)) * 20 / M_LN10) :else :psnr_r = "∞"
	if(sum_g != 0.0) :psnr_g = str(logf(255.0 / sqrt(sum_g/px/py)) * 20 / M_LN10) :else :psnr_g = "∞"
	if(sum_b != 0.0) :psnr_b = str(logf(255.0 / sqrt(sum_b/px/py)) * 20 / M_LN10) :else :psnr_b = "∞"
diff_time = timeGetTime() - start_time :timeEndPeriod 1
gsel 0
mes diff_time
mes psnr_r
mes psnr_g
mes psnr_b
```

　とはいえ、VRAMで読む作業が残ってしまうのもまた事実。時間計測したところ、次のような結果になりました(単位は[ms])。

|画像サイズ|cvcopyする部分|peekで読む部分|差の二乗を計算して足す部分|合計|
|----------|--------------|--------------|--------------------------|----|
|640x480|2|100|58|160|
|1920x1080|10|270|150|430|

# 画像のハッシュを取得してみる
　『画像にハッシュが存在するのか』と思われるかもしれませんが、実は存在します。
　SHA-256など、通常の一方向性ハッシュ関数だと「僅かにデータが異なるとハッシュ値が大きく異なる」といった特性があります。これは利点でもありますが同時に欠点でもあり、画像のように「似たデータなら似た結果を返して欲しい」場合には向きません。
　このことから、画像の特徴量取得のために様々な手法が開発されていますが、その中でも簡単かつ有効な手段の一つに「Average Hash」があります。手順としては次の通り。

- 画像を8x8など手頃なサイズに縮小(画素数がnなら、n[bit]のハッシュ値になる)
- 画像をグレースケール化
- 画素値の平均を計算
- 各画素を見ていって、画素の明るさ＜画像値の平均なら0、そうでないなら1としてビットを割り当てる

　これをHSPのみで実装するとこんな感じ。1920x1080の画像相手でも0.15[ms/回]ほどですので、その高速性が窺えるでしょう。
　ちなみに参考資料にあるように16x16で計算すると1.97[ms/回]ほどですが、lena.jpg相手に計算すると、参考資料のと比較して14bitしか違いませんでした(細かな部分の差が結構あるのに)。

```ahash.hsp
#module
	// ハッシュに使用する縮小画像のサイズを指定
	#const ahash_x 8
	#const ahash_y 8
	#const vwidth (((ahash_x * 3) + 3) & 0xfffffffc)
	// 64bitなので、ハッシュ値は32bitづつ格納する
	#deffunc ahash array hash_
		// 元画像の情報を取得する
		sel_win_id = ginfo_sel
		wx = ginfo_winx :wy = ginfo_winy
		// 縮小を掛ける
		buffer 99, ahash_x, ahash_y
		gzoom ahash_x, ahash_y, sel_win_id, 0, 0, wx, wy, 1
		// 画像を読み込み、同時にグレースケール化・平均値計算を行う
		mref vram, 66
		ddim clr, ahash_y, ahash_x
		average = 0.0
		repeat ahash_y
			y = cnt
			vpos_b = (ahash_y - cnt - 1) * vwidth
			vpos_g = vpos_b + 1
			vpos_r = vpos_b + 2
			repeat ahash_x
				temp = 0.299 * peek(vram, vpos_r) + 0.587 * peek(vram, vpos_g) + 0.114 * peek(vram, vpos_b)
				clr(y, cnt) = temp
				average += temp
				vpos_b += 3
				vpos_g += 3
				vpos_r += 3
			loop
		loop
		average = average / ahash_x / ahash_y
		// 32bitづつ格納する
		p = 0 :hash = 0
		repeat 4
			y = cnt
			repeat ahash_x
				if(p != 0) :hash = hash << 1
				if(average <= clr(y, cnt)) :hash = hash | 1
				p++
			loop
		loop
		hash_(0) = hash
		p = 0 :hash = 0
		repeat 4, 4
			y = cnt
			repeat ahash_x
				if(p != 0) :hash = hash << 1
				if(average <= clr(y, cnt)) :hash = hash | 1
				p++
			loop
		loop
		hash_(1) = hash
		gsel sel_win_id
	return
#global
```

# 画像拡大/縮小アルゴリズムの質について
　HSPでは、標準命令による拡大縮小・Artlet2Dによる拡大縮小・HSPCVによる拡大縮小などのルーチンが使用できます。
　しかし、その「質」について考えてみたことはあるでしょうか？
　例えば、画像の横幅が1/2になったとします。元画像が[CZP](http://wazalabo.com/scaling4.html)のようにみっちり高周波が詰まったものだとすると、**縮小後は外側1/2の部分は高周波を表しきれずにグレー化する**はずです。逆に、何らかの模様が外側1/2の部分に存在した場合、一般的な画像では**潰しきれない高周波成分がモアレとして現れる**ことになります。
　そこで、HSPにおける各種画像リサイズ命令を比較してみることにしました。[CZPの元画像はここにあります](http://image-d.isp.jp/commentary/czp/czp512.bmp)。
![CZP-all.png](https://qiita-image-store.s3.amazonaws.com/0/29802/1afc963e-7d55-ace7-e643-439ad9c0c73a.png)
　……こうして見ると、OpenCVバイキュービックの粗さが目立ちます。適当なリサイズソフトで「バイキュービック」指定で縮小しても、そこまで縞が目立つことはありませんでした。後で調べてみたところ、[単なるバイキュービックだとシャープ化(高周波強調)の処理が入る](http://www7a.biglobe.ne.jp/~fairytale/article/program/graphics.html#bi_cubic)ので、[事前にローパスフィルタでぼかすことでエイリアシングを軽減する](http://nodamushi.hatenablog.com/entry/20110901/1314928825)作戦があるようです。
　ちなみに画像を比較してみたところ、HSPCVにおいてCZPを50％縮小した場合、**『最近傍法＝バイキュービック、バイリニア≒モアレ軽減』**の関係が成り立ちました(40％縮小や60％縮小などではこの関係は成り立ちません)。
　また、本来OpenCVでは[Lanczosアルゴリズムによる縮小が可能](http://opencv.jp/opencv-2svn/cpp/geometric_image_transformations.html)ですが、2007年に追加された都合からか、**HSPCVに採用されているOpenCVはなんと[バージョン1.0](http://www.onionsoft.net/hsp/ref/docs/hspcv.txt)！**　道理で[Lanczosアルゴリズムが使えない](http://opencv.jp/opencv-1.0.0/document/opencvref_cv_sampling.html)わけだ。

# 画像をフーリエ変換してみる
　そもそもフーリエ変換やDCTなどは、計算量を食う処理ですのでHSPには向きません。
　ただ、[なまじDLL(APIなど)を叩けてしまう](http://www.geocities.co.jp/Playtown-Yoyo/6130/hsp/hsp-win32api.html)ことから、HSPでフーリエ変換を扱うことはできます。
　例えば、FFTライブラリとして著名なFFTWには2次元FFTを行うルーチンがありますので、それを流用すると簡単に画像をフーリエ変換することができます。

```fftw.hsp
;FFTWのAPIを呼び出すための宣言
#uselib "libfftw3-3.dll"
#cfunc fftw_malloc "fftw_malloc" int
#cfunc fftw_plan_dft_2d "fftw_plan_dft_2d" int, int, sptr, sptr, int, int
#func fftw_execute "fftw_execute" int
#func fftw_free "fftw_free" sptr
;FFTWのAPIを呼び出すための定数定義
#const fftw_complex_size 16
#const FFTW_FORWARD -1
#const FFTW_ESTIMATE 64

;対象となる画像を読み込む
screen 0, 640, 480 :title "変換前" :picload "lena.jpg"
wx = ginfo_winx :wy = ginfo_winy
;入力用バッファを確保する
in_ptr = fftw_malloc(fftw_complex_size * wx * wy)
dupptr in_data, in_ptr, fftw_complex_size  * wx * wy, vartype("double")
;画像のY成分を取り出して代入する
repeat wy
	y = cnt
	repeat wx
		x = cnt
		p = y * wx + x
		pget x, y
		in_data(p * 2) = 0.299 * ginfo_r + 0.587 * ginfo_g + 0.114 * ginfo_b
		in_data(p*2+1) = 0.0
	loop
loop
;出力用バッファを確保する
out_ptr = fftw_malloc(fftw_complex_size * wx * wy)
dupptr out_data, out_ptr, fftw_complex_size * wx * wy, vartype("double")
;FFTを行う
p = fftw_plan_dft_2d(wy, wx, in_ptr, out_ptr, FFTW_FORWARD, FFTW_ESTIMATE)
fftw_execute p
;結果を別の配列に取り出しつつ、スケーリングのために最大値・最小値を取得する
ddim power, wx, wy
power_min = 9999.0 :power_max = -9999.0
repeat wy
	y = cnt
	repeat wx
		x = cnt
		;転地処理を行うので、座標をズラしている
		p = ((y + wy / 2) \ wy) * wx + ((x + wx / 2) \ wx)
		;実部と虚部を取得し、絶対値を取りつつ対数処理している
		re = out_data(p * 2)
		im = out_data(p*2+1)
		norm = re * re + im * im
		if(norm != 0.0) :norm = logf(norm) / 2
		;結果(パワースペクトルの対数値)を配列に代入しつつ、最大値・最小値を取得
		power(x, y) = norm
		if(power_min > norm) :power_min = norm
		if(power_max < norm) :power_max = norm
	loop
loop
power_diff = power_max - power_min
;スケーリングしながら、画像として書き出す
screen 1, wx, wy
title "変換後"
redraw 0
repeat wy
	y = cnt
	repeat wx
		x = cnt
		clr = int((power(x, y) - power_min) / power_diff * 255)
		color clr, clr, clr
		pset x, y
	loop
loop
redraw 1
```
以下に実行例を示します(流石にlena画像ぐらいはググってダウンロードしてください)。
![fftw_.png](https://qiita-image-store.s3.amazonaws.com/0/29802/b4823ca2-d7d6-4cb4-5742-8c1c7035b48a.png)


# 参考資料
- [小ワザ/VRAMを直接操作してみる](http://wiki.hsp.moe/%E5%B0%8F%E3%83%AF%E3%82%B6%EF%BC%8FVRAM%E3%82%92%E7%9B%B4%E6%8E%A5%E6%93%8D%E4%BD%9C%E3%81%97%E3%81%A6%E3%81%BF%E3%82%8B.html)
- [簡単な画像の類似度計算手法「Average Hash」](http://www.toyamaguchi.com/blog/archives/20120731_average_hash/)
- [縮小・拡大処理について　その４～CZPを用いたリサイズ品質評価手法の提案～ | 技ラボ](http://wazalabo.com/scaling4.html)
- [FFTW@wiki](https://www32.atwiki.jp/amaeda/)