GSL Shell:2次元補間(interp2d)

2016-06-04 :  PCクリニック
本文の前に、
-・・・ -・-
現時点での blogramのランクインカテゴリは、
7、2、0、0、1、 0、0、0、0、0(40)で、換算ポイント 88pt 。
「Perl」「化学業界」「硝子業界」「FM COCOLO」「e-radio」「グルコサミン」、
「Python」「Firefox」bg値変動。「C言語」「FM青森」変化無し。
・-・ - -・

さて、本文。

4ヶ月前(16-02-22)の記事「GSL Shell 学習:General Data Tables ( GDT )
で、
1次元データの補間(Interpolation:インターポレーション)について書いた。

今回は、2次元データの補間(画像データに限らず)について。


「GSL - GNU Scientific Library」には、関数があるようダ。
でも、Ver.2.0 からでしょうか。

Gsl for Windows”は Vre.1.8 なので、含まれていない。

“GSL Shell”にある“libgsl-0.dll”にも無い?

本家“GNU Scientific Library”には(当然)あるが、
これは、ソース・コードのみ。

つまり、“Lua”から、簡単に利用できるライブラリは見つからない?


画像データに限定すると、“OpenCV”ではどうか?
cvResize」でしょうか?

でも、いろいろ試してみたが結果がオカシイ???


自作するか。

あちこち探したところ、

Cubic interpolation - Paulinternet.nl
を見つけた。

比較的易しいコード。“Java”と“C++”の例が載っている。

このアルゴリズムは?
次のページに載っていた。
第3章 第5節 Catmull-Romスプライン曲線
   Catmull-Romスプライン曲線は通過点だけから
  滑らかな曲線を定義する方法の一つです。
   Edwin CatmullとRaphael Romによって開発されました。
  Ferguson / Coons 曲線のように速度を与える必要がなく、
  制御点だけで曲線を決めることができるので、
  こちらの方がよく使われます。
  ・・・・・
  ・・・・・
これは、good


そう云うことで、
頑張って、自作した:序でに処理時間も測ってみた。
local time = require'time'

DblA = |x| ffi.cast( 'double*', x )
Int = |x| ffi.cast( 'int32_t', x )

My, Mx = 6, 10 -- 入力区間数
Myy, Mxx = My+1, Mx+1 -- 入力行列数
Ry, Rx = 5, 20 -- 拡大倍率(整数)
Ny, Nx = Ry*My, Rx*Mx -- 出力区間数
Nyy, Nxx = Ny+1, Nx+1 -- 出力行列数

-- 入力データ ----------------------------------
M = matrix.new( Myy, Mxx )
for y=1,Myy do
for x=1,Mxx do
M[y][x] = 10*(x-1) - (x-1)^2
end
end

-- 出力領域:横 Rx==10 倍拡大、縦 Ry==3 倍拡大
N = matrix.new( Nyy, Nxx )

---------------------------------------------------------------------------------
-- a x^3 + b x^2 + c x + d <--- 4点:p0, p1, p2, p3 :[ p1 ~ p2 ] 間補間
function C_abcd( p0, p1, p2, p3 )
local a = -1/2*p0+3/2*p1-3/2*p2+1/2*p3
local b = p0-5/2*p1+2*p2-1/2*p3
local c = -1/2*p0+1/2*p2
local d = p1
return a,b,c,d
end
function C_val( a, b, c, d, x ) -- x : 0.0 ~ 1.0
return a*x^3+b*x^2+c*x+d
end

function Row( y ) -- 1行分処理:入力 M → 出力 N y:0 ~ My-1
local Mptr = DblA( Int(M.data)+8*Mxx*y )
local Nptr = DblA( Int(N.data)+8*Nxx*Ry*y )
local a, b, c, d = C_abcd( 2*Mptr[0]-Mptr[1], Mptr[0], Mptr[1], Mptr[2] )
for k=0,Rx-1 do Nptr[k] = C_val( a, b, c, d, k/Rx ) end
for x=1,Mx-2 do
a, b, c, d = C_abcd( Mptr[x-1], Mptr[x], Mptr[x+1], Mptr[x+2] )
for k=0,Rx-1 do Nptr[x*Rx+k] = C_val( a, b, c, d, k/Rx ) end
end
a, b, c, d = C_abcd( Mptr[Mx-2], Mptr[Mx-1], Mptr[Mx], 2*Mptr[Mx]-Mptr[Mx-1] )
for k=0,Rx do Nptr[(Mx-1)*Rx+k] = C_val( a, b, c, d, k/Rx ) end
end

function Col( x ) -- 1列分処理:入力 N → 出力 N x:0 ~ Nx-1
local Nptr = DblA( Int(N.data)+8*x )
local a, b, c, d = C_abcd( 2*Nptr[0]-Nptr[Nxx*Ry], Nptr[0], Nptr[Nxx*Ry], Nptr[Nxx*2*Ry] )
for k=1,Ry-1 do Nptr[Nxx*k] = C_val( a, b, c, d, k/Ry ) end
for y=1,My-2 do
a, b, c, d = C_abcd( Nptr[Nxx*Ry*(y-1)], Nptr[Nxx*Ry*y], Nptr[Nxx*Ry*(y+1)], Nptr[Nxx*Ry*(y+2)] )
for k=1,Ry-1 do Nptr[Nxx*(Ry*y+k)] = C_val( a, b, c, d, k/Ry ) end
end
a, b, c, d = C_abcd( Nptr[Nxx*Ry*(My-2)], Nptr[Nxx*Ry*(My-1)], Nptr[Nxx*Ry*My], 2*Nptr[Nxx*Ry*My]-Nptr[Nxx*Ry*(My-1)] )
for k=1,Ry-1 do Nptr[Nxx*(My-1+k)] = C_val( a, b, c, d, k/Ry ) end
end

st=time.time() ---------
for y=0,My-1 do Row(y) end ---------
for x=0,Nx-1 do Col(x) end ---------
et=time.time() ---------
print( et-st )

------------------------------------------------ 1行分について、グラフ化 -------------
x = matrix.new( Mxx, 1 ); for k=1,Mxx do x[k]=Rx*(k-1) end -- サンプルデータのX座標値
X = matrix.new( Nxx, 1 ); for k=1,Nxx do X[k]=k-1 end -- 結果データのX座標値

p = graph.plot( 'Graph' ) -- グラフタイトル
g = graph.xyline( x, M[1] )
p:addline( g, 0x0000FFFF, {{'marker', size=8}} )
gg= graph.xyline( X, N[1] )
p:addline( gg )
p:show()
これで、行けそう。
処理時間も、大丈夫。


但し、(画像の様に)等間隔メッシュの場合のみ。
上記の以前のブログの1次元の場合の、
 { 0, 2, 3, 4, 8 }
みたいに不等間隔なデータは扱え無い。


それについては、別途考えよう。


本日はここまで。


Lua ( GSL Shell ) 学習は続く。


見ていただいた序でとは厚かましい限りですが、
お帰りに投票して頂けるとなお嬉しいです。 ⇒ blogram投票ボタン


160502
関連記事
スポンサーサイト

コメントの投稿

管理者にだけ表示を許可する

おきてがみ/blogram
blogram投票ボタン



おきてがみ

最新記事
カレンダー
09 | 2017/10 | 11
1 2 3 4 5 6 7
8 9 10 11 12 13 14
15 16 17 18 19 20 21
22 23 24 25 26 27 28
29 30 31 - - - -
月別アーカイブ
カテゴリ
最新コメント
検索フォーム
リンク
プロフィール

<紙>

Author:<紙>
ようこそ。
「パソコンヲタクの雑記帳」
もろもろなことを綴っています。
パソコン ヲタクってねくら?
画像は kami でなく kani です。

カウンター(fc2、i2i) /Google Analytics


i2i(from 2010-08-24)
Total =
Today  =  
Yesterday=
アンチエイジング

Google Analytics
ブックマーク