Scilab SIP - eiichiromomma/CVMLAB GitHub Wiki
(Scilab) SIP
Scilabで画像処理。Scilab Image Processing toolboxを使う。
※SIVPの方がオススメ
準備
インストール
Windowsの場合インストーラを使って入れるだけ。
使い方
README.txtには
exec loader.sce
しろと書いてあるが、ディレクトリを移動しないとファイルが見付からないと怒られる。 [toolboxes]メニューで[siptoolbox]を選ぶのが手っ取り早い。
ロードに失敗
なぜかCORE_RL_magic_.dllのロードだけ失敗する。
link(ImageMagickPath + 'CORE_RL_magick_.dll');
link failed for dll C:\PROGRA~1\SCILAB~1.1\contrib\siptoolbox\lib\CORE_RL_magick_.dll
link(ImageMagickPath + 'CORE_RL_magick_.dll');
!--error 236
link: the shared archive was not loaded
at line 21 of exec file called by :
exec("C:\PROGRA~1\SCILAB~1.1\contrib\siptoolbox/loader.sce")
in execstr instruction called by :
execstr(toolboxes(1))
while executing a callback
のようなエラーが出る。 原因は不明だが、dllのある場所まで移動してdll指定するとちゃんとロードされる。
loader.sceの変更
C:\Program Files\scilab-4.1.1\contrib\siptoolboxのloader.sceの
link(ImageMagickPath + 'CORE_RL_magick_.dll');
の行を以下のように変更する。
//link(ImageMagickPath + 'CORE_RL_magick_.dll');
chdir(ImageMagickPath);
link('CORE_RL_magick_.dll');
で、再びメニューからsiptoolboxを呼ぶ
デモ
--------------------------------------------
SIP - Scilab Image Processing Toolbox loaded
Enter "exec(SIPDEMO);" for demonstrations.
と表示されれば成功
exec(SIPDEMO)
とするとデモを見られる。
##SIPで用意されている関数について
メモ
- angle:横、縦の微分画像から画素毎の勾配を画素値とする画像を返す(筈)。なぜか未定義と言われる。
- bwborder:2値画像について4あるいは8近傍の境界を抽出。
- bwdist:2値画像についての境界からの距離画像の作成。形状の概要の作成に使える。
- bwlabel:2値画像についてのラベリング処理。ラベル画像とラベルの最大値を返す。
- curvature:輪郭の各点における曲率
- drawline:画像に線を引く
- edilate:膨張処理
- erosion:収縮処理。edilateと組合せてノイズ除去。
ユーザガイド
- SIP Documentation and Contribs のpre-doc.pdfを追ってみる。
インストール
(省略)
とりあえず使えるようになる
スタックの拡張
SIPはメモリを大量に使うので予めスタックを拡張しておく。
stacksize(3e7)
それか.scilabファイルに
stacksize(3e7);
を記述しておく。
画像のロード
何でも良いのだがSIPの画像を使ってみる。
chdir(SIPDIR+'images')
mat=imread('ararauna.jpg)');
chdirはディレクトリの移動。imreadは画像の読み込み。
[mat,map]=imread('ararauna.png)');
とすればカラーマップも取得できる。
画像の表示
xbasc();imshow(mat);
あるいは
xbasc();imshow(mat,map);
xbascは画像ウィンドのクリア。
画像の保存
imwrite(mat,'/home/druel/test.jpg)')
あるいは
imwrite(mat,map,'/home/druel/test.jpg)')
/home/druelは各自の環境の置き換えて、保存可能な場所を指定する。
正規化
dmat=mat*0.6;
xset('window',2);xbasc();imshow(dmat);
とすると、画像がやや暗くなる。(濃度値の0.6倍)
max(dmat)
とすると、0.6を返す。(SIPにおいても画像は0〜1.0の濃度値になる)
nmat=normal(dmat);
で正規化される。
xset('window',3);xbasc();imshow(nmat);
とすると元の画像と同じ画像を得る。(元の画像のmaxが1.0の場合)
画像上の軸について
左上が原点になる。
image=imread(SIPDIR+'images/ararauna.jpg)');
[r c]=size(image);
画像を読んでr,cに行、列数を取得。
image(10,:,1)=1;
image(:,10,2)=1;
10行目に赤のライン、10列目に緑のラインを引く。
基本的な処理
※SIPの画像ではなく、一般的な画像の話になる。
2画像の合成(平均画像)
mat1=imread('image1.jpg)');mat2=imread('image2.jpg');
mat_mean=(mat1+mat2)/2;
xbasc(); imshow(mat_mean);
画像image1,image2を読んで平均をmat_meanに入れて表示。
画素を二乗した画像
mat=imread('image.jpg)');
pmat=mat.*mat;
MATLABの基本だが、要素同士の乗算は".*"になる事に注意
256を法とする加算、減算
要するに8bit画像として扱う方法(0〜255)
mat1=imread('image1.jpg)');mat2=imread('image2.jpg');
mat_modulo=uint8(mat1*255)-uint8(mat2*255);
でmat_moduloはmat1-mat2の結果の8bit画像になる。
mat_modulo=double(mat_modulo)/255;
SIPで用いるdouble画像にする場合はdoubleにして255で割る。
色の扱い
再びSIPの画像を例に扱う。
グレースケール
gray_imreadを使う。
mat=gray_imread('onca.gif');
xbasc();imshow(mat);
size(mat)
画像は8,16bitのいずれかが利用可能で、0〜1.0に正規化される。
フルカラー
imreadを使う。
mat=imread('ararauna.jpg)');
xbasc();imshow(mat);
size(mat)
mat(:,:,1)が赤、mat(:,:,2)が緑、mat(:,:,3)が青のチャネルを示す。
疑似カラーの概要
再びonca.gifを使う。(gifは256色のカラーテーブルを使っているため)
mat=imread('onca.gif');
として
min(mat)
とすると1を得て、
max(mat)
とすると256を得る。つまり1〜256である点に注意
xbasc();imshow(mat,graycolormap(256));
読み込み時にmapを取得して、imshowに渡すことも出来るが、自由にカラーマップを割り当てられる。 graycolormap(256)は256階調のグレースケールを作成している。
map=hotcolormap(256);
とすると暑苦しいカラーマップになる。これを使うには
xset('window',1);xbasc();imshow(image,map);
とする。
更に、
invmap(256:-1:1,1:3)=map(1:1:256,1:3);
とするとカラーマップが反転する
xbasc();imshow(image,invmap);
で表示。
hotcolormapについて1からの差分をとると、涼しげなカラーマップも作成できる。
xbasc();imshow(image,1-hotcolormap(256));
デフォルトではgraycolormapとhotcolormapしか存在しないが、EnricoというScilab用Toolboxを用意すると他のカラーマップも使えるようになる。
カラーマップの確認法
pseudo=ones(1:256)'*(1:256);
xbasc();imshow(pseudo,hotcolormap(256));
のような感じ。
疑似カラー
画像がカラーテーブルを持っている場合、
image=imread('ararauna.png)');
xbasc();imshow(image,graycolormap(256));
のようにあてずっぽうにカラーマップを割り当てても綺麗なグレー画像にならない。
[image,map]=imread('ararauna.png)');
xbasc();imshow(image,map);
のようにmapを取得する必要がある。
前述の通り、
pseudo=ones(1:256)'*(1:256);
xbasc();imshow(pseudo,map)
とするとカラーマップが見られる。
保存するには
imwrite(image,map,filename);
としてmapも渡す必要がある。
カラーからグレースケールへの変換
mat=imread('ararauna.jpg)');
gmat=im2gray(mat);
xset('window',1);xbasc();imshow(gmat);
im2grayを用いるか、
gmat=gray_imread('ararauna.jpg)');
のように初めからグレーで読む。
しきい値処理
bmat=im2bw(gmat,0.7);
とすると、bmatは0.7をしきい値とした2値画像になる。
フィルタ処理
畳み込み演算
ここでは3x3の演算子をkernelと呼ぶ。
| 1 | 5 | 1 |
|---|---|---|
| 2 | 4 | 3 |
| 8 | 2 | 1 |
のような画像の一部があり、kernelが
| 1 | 1 | 1 |
|---|---|---|
| 1 | 1 | 1 |
| 1 | 1 | 1 |
のように平滑化処理の場合、演算後の画像の中心の画素値は(1x1)+(1x5)+(1x1)+(1x2)+(1x4)+(1x3)+(1x8)+(1x2)+(1x1)=27を画素数9で割った3となる。(この場合割るのは手動)
| 1/9 | 1/9 | 1/9 |
|---|---|---|
| 1/9 | 1/9 | 1/9 |
| 1/9 | 1/9 | 1/9 |
というkernelを作れば良い。
ただし、1チャネル画像のみ対応
一般的な使い方
画像matを取得した場合
mat=gray_imread('ararauna.jpg)');
kern=[1 1 1;1 1 1;1 1 1]
mat_conv=(1/9)*imconv(mat,kern);
xbasc(); imshow(mat_conv);
のようにimconvを使う。
既存のフィルタ
SIPで用意されているフィルタ
mkfilter('名前')
とすればkernelが作られる。
ローパスフィルタ
- mean:平均値フィルタ
- low-pass:meanよりユルい平滑化フィルタ
ラプラシアン(エッジ抽出)
- laplace1:4近傍のラプラシアン
- laplace2:8近傍のラプラシアン
- laplace3:重み付けした8近傍のラプラシアン
非線形フィルタ:メディアン
中間値フィルタ。雑音の除去に適するが処理に時間がかかる。
mat=gray_imread('ararauna.jpg)');
mmat=mogrify(mat,['-median', '1']);
xset('window',1);xbasc();imshow(mat);
xset('window',2);xbasc();imshow(mmat);
として使う。-medianの引数は中心画素からの距離。1なら3x3、2なら5x5になる。 ※ドキュメントのmogrifyの引数は'-median 1'となっているが誤りで、引数が2つ以上の場合は行列として渡す必要がある。
非線形フィルタ:ヒストグラムの平滑化
fmat=mogrify(mat,'-equalize');
xset('window',3);xbasc();imshow(fmat);
でヒストグラムの平滑化になる。
xset('window',4);histplot(64,mat);
xset('window',5);histplot(64,fmat);
とすると、比較が可能。
二次元FFT処理
スペクトルの取得
とりあえず小さいスリットフィルタの空間周波数スペクトルを調べてみる
image=zeros(400,300);
image(180:220,145:155)=1;
として300x400の黒い画像に11x41の窓(値が1)を作る。 表示すると
xset('window',0);xbasc();imshow(image);
となる。
FFTを行ない
IM=fft(image,-1);
共役複素数を要素ごとに乗算する。
spectrum=real((IM).*conj(IM));
xset('window',1);xbasc();imshow(spectrum,[]);
このままだと、四隅が原点(直流)で見辛いので
spectrum2=sip_fftshift(spectrum);
xset('window',2);xbasc();imshow(spectrum2,[]);
とすると、中心が直流成分になる。
また、FFTの結果は数値がデカいので対数にすると概要が掴みやすい。
spectrum3=log(spectrum2+1);
xset('window',3);xbasc();
plot3d1(1:4:400,1:4:300,spectrum3(1:4:400,1:4:300));
とすると3次元表示が出来る。
あるいはしきい値で切って飽和させてしまう手もある。
spectrum3=spectrum2;threshold=1e3;
しきい値thresholdを1000にして
spectrum3(find(spectrum3>threshold))=threshold;
それより大きい値を1000にしてしまう。
xset('window',3);xbasc();
plot3d1(1:4:400,1:4:300,spectrum3(1:4:400,1:4:300));
とすると違いが分かる。
逆変換は
im2=real(fft(IM,1));
xset('window',4);xbasc();imshow(im2)
とすることで元画像に戻せるが、その前にフィルタ処理をすれば、空間周波数レベルでの処理が可能となる。
画像を読んでスペクトルを表示
image=gray_imread(SIPDIR+'images/ararauna.jpg)');
[r,c]=size(image);
xset("window",0);xbasc();imshow(image);
IM=fft(image,-1);
spectrum=real((IM).*conj(IM));
xset("window",1);xbasc();imshow(spectrum,[]);
sp2=sip_fftshift(spectrum);
//to center the spectrum
xset("window",2);xbasc();imshow(sp2,[]);
sp3=log(sp2+1);
xset("window",3);xbasc();imshow(sp3,[]);
円形のフィルタを作成
//design the binary filter
z=zeros(image);
x0=round(r/2);radius=10;y0=round(c/2);
for x=x0-radius:x0+radius
y=round(sqrt(radius^2-(x-x0)^2));
z(x,y0-y:y0+y)=1;
end;
//if you want to inverse the filter
//z=abs(1-z);//complementary filter
xset("window",4);xbasc();imshow(z,[]);
逆FFT
IM2=IM.*sip_fftshift(z);//spectrum modification
//reverse transform
im2=real(fft(IM2,1));
xset("window",5);xbasc();imshow(im2,[]);
のようになる
※sip_fftshiftの対象に注意 IMは中心が高周波なので、zもsip_fftshiftして乗算する必要がある。fft関数も中心を高周波として扱う。
逆にやってしまうと
IM3=sip_fftshift(IM).*z;//spectrum modification
im3=real(fft(IM3,1));
xset("window",6);xbasc();imshow(im3,[]);
のように変な結果になる。
ちなみにフィルタはmkfftfilterなる便利な関数があるので、例にあるようなループでわざわざ作成する必要は無い。 mkfftfilterについては
help mkfftfilter
参照。
Cからの利用
必要が生じたら書く予定。