行列演算ライブラリTiny
ちょっとSPALMで行列演算が使いたくなったので作ってみた。ただしSPALMmoreでないと手直し無しでは動かないと思う。
このライブラリはオブジェクト指向？な感じにできていて、ベクトルやマトリクスの「ハンドル」を取得してから使います。
ハンドルの実体はただのインデックスなんですが、これでSPALM上でも擬似的に複数のインスタンスを作るなんてことが簡単にできます。もちろん取得したハンドルは使用後解放しなくてはなりませんが、（しなくてもいいけど･･･）
下のサンプルは３次曲線を描画するものですが、行列のアフィン変換により曲線を回したり、拡大縮小したり、せん断して表示することができます。
上下・・・曲線のパラメーター変化
左右・・・回転
1、2・・・拡縮
4、5・・・せん断（X方向）
tmp=0
func addTmp(hm1){'行列hm1を一時オブジェクトとして登録'
global tmp[tmp++]=hm1
return hm1
}
func remTmp(){'行列の一時オブジェクト全削除'
while(tmp){removeMat(tmp[--tmp])}
}
func printMat(hm1,x,y){'行列の内容表示'
for(i=0;i<4;++i){
for(j=0;j<4;++j){
text(mAt[hm1][i][j],x+i*60+60,y+j*12,RIGHT)
}}
}
func printVec(hv1,x,y){'ベクトルの内容表示'
for(i=0;i<4;++i){
text(vEc[hv1][i],x+i*60+60,y,RIGHT)
}
}
func vecMat(hv1,hm1){'ベクトルと行列の掛け算hv1に結果は格納される'
tmp=getVec(hv1)
for(i=0;i<4;vEc[tmp][i++]=k){
for(k=j=0;j<4;k+=vEc[hv1][j]*mAt[hm1][i][j++]){
}}
setVec(hv1,tmp)
removeVec(tmp)
return hv1
}
func getMat(k){'k以降で有効な行列のハンドルを返す'
global mAt=0
while(global mAt[++k]){}
mAt[k]=1
@(mAt[k][0],1,0,0,0)
@(mAt[k][1],0,1,0,0)
@(mAt[k][2],0,0,1,0)
@(mAt[k][3],0,0,0,1)
return k
}
func endMat(i){'行列iの使用を終了'
mAt[i]=0
}
func removeMat(hm1){'行列hm1を削除'
for(i=0;i<4;++i){
for(j=0;j<4;++j){
disvar(mAt[hm1][i][j])
}}
disvar(mAt[hm1])
}
func addMat(hm1,hm2,hm3){'行列の足し算 hm1 = hm2 + hm3'
for(i=0;i<4;++i){
for(j=0;j<4;mAt[hm1][i][j]=
mAt[hm2][i][j]+
mAt[hm3][i][j]
,++j){
}}
return hm1
}
func setMat(hm1,hm2){'hm1にhm2の内容をコピー'
for(i=0;i<4;++i){
copy(mAt[hm2][i],0,mAt[hm2][i],0,4)
}
return hm1
}
func mulMat(hm1,hm2,hm3){'行列の掛け算 hm1 = hm2 * hm3'
for(i=0;i<4;++i){
for(j=0;j<4;
mAt[hm1][i][j]=
mAt[hm2][0][j]*mAt[hm3][i][0]+
mAt[hm2][1][j]*mAt[hm3][i][1]+
mAt[hm2][2][j]*mAt[hm3][i][2]+
mAt[hm2][3][j]*mAt[hm3][i][3]
,++j){
}}
return hm1
}
func identMat(hm1){'hm1を単位行列化'
@(mAt[hm1][0],1,0,0,0)
@(mAt[hm1][1],0,1,0,0)
@(mAt[hm1][2],0,0,1,0)
@(mAt[hm1][3],0,0,0,1)
return hm1
}
func rotXMat(hm1,r){'hm1をX軸に関してr度回転させる行列にする'
@(mAt[hm1][0],1,0,0,0)
@(mAt[hm1][1],0,lc=cos(r)/100.0,-ls=sin(r)/100.0,0)
@(mAt[hm1][2],0,ls,lc,0)
@(mAt[hm1][3],0,0,0,1)
return hm1
}
func rotYMat(hm1,r){'hm1をY軸に関してr度回転させる行列にする'
@(mAt[hm1][0],lc=cos(r)/100.0,0,ls=sin(r)/100.0,0)
@(mAt[hm1][1],0,1,0,0)
@(mAt[hm1][2],-ls,0,lc,0)
@(mAt[hm1][3],0,0,0,1)
return hm1
}
func rotZMat(hm1,r){'hm1をZ軸に関してr度回転させる行列にする'
@(mAt[hm1][0],lc=cos(r)/100.0,-ls=sin(r)/100.0,0,0)
@(mAt[hm1][1],ls,lc,0,0)
@(mAt[hm1][2],0,0,1,0)
@(mAt[hm1][3],0,0,0,1)
return hm1
}
func transMat(hm1,x,y,z){'hm1をx,y,z分平行移動'
return
(mAt[hm1][0][3]+=x,
mAt[hm1][1][3]+=y,
mAt[hm1][2][3]+=z,
hm1)
}
func zoomMat(hm1,z){'hm1をXYZ軸に関してz倍'
for(i=0;i<4;mAt[hm1][0][i]*=z,
mAt[hm1][1][i]*=z,
mAt[hm1][2][i]*=z
,++i){
}
return hm1
}
func getVec(i){'i以降で有効なベクトルのハンドルを返す'
glb vEc=0
while(global vEc[++i]){}
vEc[i]=1
@(vEc[i],0,0,0,1)
return i
}
func endVec(i){'ハンドルiのベクトルオブジェクトの使用を終了'
vEc[i]=0
}
func removeVec(i){'ハンドルiのベクトルオブジェクトを削除'
for(j=0;j<4;++j){
disvar(vEc[i][j])
}
disvar(vEc[i])
}
func setVec(hv1,hv2){'hv1にhv2の内容をコピー'
copy(vEc[hv2],0,vEc[hv1],0,4)
return hv1
}
func normalVec(hv1){'hv1を単位化'
return mulVec(hv1,1.0/lenVec(hv1))
}
func lenVec(hv1){'hv1の大きさを返す'
return
sqrt(pow(vEc[hv1][0],2)+pow(vEc[hv1][1],2)+pow(vEc[hv1][2],2)+pow(vEc[hv1][3],2))
}
func mulVec(hv1,m){'hv1の各要素にmをかける'
return
(vEc[hv1][0]*=m,
vEc[hv1][1]*=m,
vEc[hv1][2]*=m,
vEc[hv1][3]*=m,
hv1)
}
func addVec(hv1,hv2){'hv1にhv2を足す'
return
(vEc[hv1][0]+=vEc[hv2][0],
vEc[hv1][1]+=vEc[hv2][1],
vEc[hv1][2]+=vEc[hv2][2],
vEc[hv1][3]+=vEc[hv2][3],
hv1)
}
'ここまでライブラリ'
M=getMat()
Z=0.8
R=0
c=1
X=0
label 1
mAt[identMat(M)][0][1]=X
M=mulMat(getMat(),transMat(addTmp(getMat()),-120,-120,0),addTmp(M))
zoomMat(M,Z)
M=mulMat(getMat(),addTmp(M),rotZMat(addTmp(getMat()),R))
transMat(M,120,120,0)
v0=getVec()
v1=getVec()
vecMat(v0,M)
'setVecElem(v1,0,1)'
vEc[v1][0]=1
vecMat(v1,M)
sx=vEc[v0][0]
sy=vEc[v0][1]
dxx=vEc[v1][0]-sx
dyx=vEc[v1][1]-sy
vEc[v1][0]=0
vEc[v1][1]=1
vecMat(v1,M)
dxy=vEc[v1][0]-sx
dyy=vEc[v1][1]-sy
removeVec(v0)
removeVec(v1)
remTmp()
label 0
lock()
'clear(0,0,240,240)'
clear(0,0,width,height)
a=2*(c-1.0)
b=3*(1.0-c)
x=0
y=0
for(i=0.0;i<1.0;x=nx,y=ny){
nx=239*i=i+0.125
ny=((a*i+b)*i+c)*i*239
line(int(sx+dxx*x+dxy*y),int(sy+dyx*x+dyy*y),int(sx+dxx*nx+dxy*ny),int(sy+dyx*nx+dyy*ny))
}
'text(c_" 倍率:"_setscl(Z,2)_" 角度:"_R_" Xせん断:"_X,0,0,0)'
text(c_" 倍率:"_Z,0,0,0)
text("角度:"_R_" Xせん断:"_X,0,15,0)
unlock(1)
switch(scan|input(0)){
case 0x1000:
c-=0.125 break
case 0x8000:
c+=0.125 break
case 0x2000:
R-=5 goto 1
case 0x4000:
R+=5 goto 1
case 2:
Z-=0.1 goto 1
case 4:
Z+=0.1 goto 1
case 16:
X-=0.1 goto 1
case 32:
X+=0.1 goto 1
case 0:break
default:
removeMat(M)
end
}
goto 0
