ARToolKit で影をつける。その２

以前の結果が間違っていたので訂正します。
http://d.hatena.ne.jp/fortran66/20090214

glutSolidCube を用いると、座標原点が動かされるのか、影を作る射影行列が正しく働かないようです。
- そこで glVertex3f を用いて描画するようにしました。
射影面の法線ベクトルの第 4 成分には、原点から射影平面への距離を負にしたものが入ります。
光源の位置ベクトルの第 4 成分には、普通の座標点と同じくスケールファクターが入ります。ゆえにそれが 0 の時は、第 1 〜 3 成分 (x､y､z) 方向の無限遠に光源が存在することになり、それが 1 の時は、(x,y,z) 位置に点光源が存在することになるようです。
- OpenGL の 4 元位置ベクトルでは第 4 成分がスケールファクターとなっていて、(x,y,z,w) = (x/w, y/w, z/w, 1.0) というスケール関係にあります。

■影を生成する射影行列の導出

まず射影面の原点を通る場合を考えるとして、以下の三次元空間のベクトルを定義する。
$|l>$ 光源の座標ベクトル。
$|n>$ 原点から射影面までの距離ベクトル。（法線ベクトル）
$|x>$ 影を生じさせる物体の座標ベクトル。
$|s>$ 物体によって生じた影の座標ベクトル。

今求めたいものは、物体の座標を影の座標に変換する射影演算子。これを $M$ とすると、
$|s>=M|x>$ 　　　　　　　　　　　　(1)
と書ける。

影は射影面上にあるので、影のベクトル $|s>$ は射影面の法線ベクトル $|n>$ と直交する。ゆえに、それらの内積は消えて
[tex:=0] 　　　　　　　　　　　　(2)
となる。

影は光源と物体を結ぶ直線状にある。その直線はパラメータ $t$ を用いると、
$|s>=t|x>+(1-t)|l>$ 　　　　　　　　(3)
とあらわせる。

式 (2),(3) より、
[tex:=t+(1-t)=0] 　　　　(4)
となる。

これを満たすパラメータ $t$ は、比例定数 $\alpha$ を導入して、
[tex:t=\alpha ,]　　　　　　　　　　(5)
[tex:1-t=-\alpha ]　　　　　　　　　(6)
とあらわせる。

ここで式 (5),(6) から $\alpha$ を求めると、
[tex:\alpha=(-)^{-1}]　　　　　(7)
となる。

以上の結果から影の座標は、
[tex:|s>=\alpha(|x>-|l>)]　　　(8)
と求まる。

これを形式的に $|x>$ をくくり出した形で書き直すと
[tex:|s>=\alpha(I-|l>]　　　　(9)
となる。ただし、 $I$ は単位行列である。式 (1) と比較すると、射影の演算子 $M$ は
[tex:M=\alpha(I-|l> とあらわせる。

ここで、OpenGL の 4 元ベクトルを用いることを考えると、比例定数[tex:\alpha]を第 4 成分に引き受けさせることができる。[tex:|l>,|x>]の第 4 成分が 1 であることに注意すると、求めるべき第 4 成分も三次元ベクトルの時と同じ形式となることが分かるので、
[tex:M=(I-|l> と書きあらわせる。

また射影面が原点を通らない場合を考えた場合も、法線ベクトル[tex:|n>]の第 4 成分を、原点から平面までの距離に三次元法線ベクトルのノルムを掛けて負にしたものとすることで、式 (11) がそのまま成立することが導ける。

参照:
http://son-son.sakura.ne.jp/programming/opengl_1.html
http://www.opengl.org/resources/code/samples/mjktips/TexShadowReflectLight.html

□追記：式 (11) が原点を通らない平面でもそのまま成立すること

原点から射影平面へ直交する三次元のベクトルは、射影面のベクトルであるが、OpenGL の 4 元ベクトルにするときに、その第 4 成分を三次元ベクトルのノルムの二乗にしておく。これは原点から射影面までの距離の二乗に等しい。

射影面の法線ベクトルを成分にわけて書いたとき
$|n>=|n_x, n_y, n_z, n_w>$ 　　　　　　　　　(A.1)
であるとするとき、あらたに
$|n_0>=|n_x, n_y, n_z, 0>$ ,　　　　　　　　　(A.2)
を定義する。

ここで、定義によって
[tex:n_w=-]　　　　　　　　　(A.3)
の関係がある。

また座標原点を射影面に移したベクトルを $'$ を付けることにして、
$|l'>=|l>-|n_0>$ ,　　　　　　　　　(A.4)
$|x'>=|x>-|n_0>$ ,　　　　　　　　　(A.5)
$|s'>=|s>-|n_0>$ ,　　　　　　　　　(A.6)
を定義する。

さて座標原点を射影面上に持ってきているので、影の生じる点は (8) より
[tex:|s'>=|x'>-|l'>]　　　　　　　　　(A.7)
とあらわせる。

今、
[tex:==-=-n_w-=],　(A.8)　　
[tex:==-=-n_w-=]　 (A.9)
とまとめられる。

これを用いて、
[tex:|s'>=|x'>-|l'>=|x>-|l>-(-)|n_0>] (A.10)
と書き直せる。（ただし、右辺は全体を合わせた時のみ意味がある。）

したがって、
[tex:|s>=|s'>+|n_0>=(|x>-|l>-|n_0>)+|n_0>=|x>-|l>]　　　　　　　　　(A.11)
と求める結果を得る。

解き方が第 4 成分の取り扱いに関して不自然ですっきりしません。OpenGL の第 4 成分を含めた代数がよく分かりません。普通の３次元ベクトルの形から始めて、だんだんと OpenGL の 4 元ベクトルに直して行ったので複雑になったかもしれません。不備なところはチラシの裏の走り書きということで勘弁してｗ

■実行結果

■Fortran のソース

光源（赤玉）が周回運動をするようにしてあります。

MODULE m_ART
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
!---- TYPE definitions ---------------------------------
TYPE, BIND(C) :: t_ARParam
 INTEGER(C_INT) :: ixsize, iysize
 REAL(C_DOUBLE) :: dmat(4, 3)
 REAL(C_DOUBLE) :: dist_factor(4)
END TYPE
!
TYPE, BIND(C) :: t_ARMarkerInfo 
 INTEGER(C_INT) :: iarea
 INTEGER(C_INT) :: id
 INTEGER(C_INT) :: dir
 REAL(C_DOUBLE)  :: cf
 REAL(C_DOUBLE)  :: pos(2)
 REAL(C_DOUBLE)  :: dline(3, 4)
 REAL(C_DOUBLE)  :: vertex(2, 4)
END TYPE
!---- global variables ---------------------------------
INTEGER(C_INT) :: ithresh = 100
INTEGER(C_INT) :: icount = 0
INTEGER(C_INT) :: ipatt_id
REAL(C_DOUBLE) :: patt_width = 80.0_c_double
REAL(C_DOUBLE) :: patt_center(2) = (/0.0_c_double, 0.0_c_double/)
REAL(C_DOUBLE) :: patt_trans(3, 4)
INTEGER(C_INT) :: ixsize, iysize
TYPE(t_ARParam) :: cparam
!---- INTERFACE BLOCK -----------------------------------
INTERFACE
 SUBROUTINE arVideoCapStop() BIND(C, NAME = 'arVideoCapStop')
 END SUBROUTINE arVideoCapStop
! 
 SUBROUTINE arVideoClose() BIND(C, NAME = 'arVideoClose')
 END SUBROUTINE arVideoClose
!
 SUBROUTINE argCleanup() BIND(C, NAME = 'argCleanup')
 END SUBROUTINE argCleanup
! 
 INTEGER(C_INT) FUNCTION arVideoGetImage() BIND(C, NAME = 'arVideoGetImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE arUtilSleep(i) BIND(C, NAME = 'arUtilSleep')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: i 
 END
!
 SUBROUTINE arUtilTimerReset() BIND(C, NAME = 'arUtilTimerReset')
 USE, INTRINSIC :: ISO_C_BINDING 
 END
!
 SUBROUTINE argDrawMode2D() BIND(C, NAME = 'argDrawMode2D')
 END 
!
 SUBROUTINE argDispImage( iptr, i, j ) BIND(C, NAME = 'argDispImage')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER, VALUE :: i, j
 END 
!
 INTEGER(C_INT) FUNCTION arDetectMarker( iptr, i, mark, k ) BIND(C, NAME = 'arDetectMarker')
 USE, INTRINSIC :: ISO_C_BINDING 
 IMPORT :: t_ARMarkerInfo
 INTEGER(C_INTPTR_T), VALUE :: iptr
 INTEGER(C_INT)     , VALUE :: i
 TYPE(t_ARMarkerInfo), POINTER :: mark(:)
 INTEGER(C_INT) :: k
 END 
!
 SUBROUTINE arVideoCapNext() BIND(C, NAME = 'arVideoCapNext')
 END 
!
 SUBROUTINE argSwapBuffers() BIND(C, NAME = 'argSwapBuffers')
 END 
!
 SUBROUTINE arGetTransMat(ptr, d1, d, d2) BIND(C, NAME = 'arGetTransMat')
 USE, INTRINSIC :: ISO_C_BINDING 
 TYPE(C_PTR), VALUE :: ptr
 REAL(C_DOUBLE), VALUE :: d
 REAL(C_DOUBLE):: d1(2), d2(3, 4)
 END 
!
 SUBROUTINE argDrawMode3D() BIND(C, NAME = 'argDrawMode3D')
 END SUBROUTINE argDrawMode3D
 !
 SUBROUTINE argDraw3dCamera(ix, iy) BIND(C, NAME = 'argDraw3dCamera')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER(C_INT), VALUE :: ix, iy 
 END
 !
 SUBROUTINE argConvGlpara(x, y) BIND(C, NAME = 'argConvGlpara')
 USE, INTRINSIC :: ISO_C_BINDING 
 REAL(C_DOUBLE) :: x(4, 3), y(4, 4) 
 END
! 
 SUBROUTINE arVideoCapStart() BIND(C, NAME = 'arVideoCapStart')
 END
!
 SUBROUTINE argMainLoop( i, j, k ) BIND(C, NAME = 'argMainLoop')
 USE, INTRINSIC :: ISO_C_BINDING 
 INTEGER, VALUE :: i
 TYPE(C_FUNPTR), VALUE :: j, k
 END
!
 INTEGER(C_INT) FUNCTION arVideoOpen(text) BIND(C, NAME = 'arVideoOpen')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(C_CHAR) :: text(*)
 END
!
 INTEGER(C_INT) FUNCTION arVideoInqSize(ix, iy) BIND(C, NAME = 'arVideoInqSize')
 USE, INTRINSIC :: ISO_C_BINDING
 INTEGER(C_INT) :: ix, iy
 END
! 
 INTEGER(C_INT) FUNCTION arParamLoad(text, i, xparam) BIND(C, NAME = 'arParamLoad')
 USE, INTRINSIC :: ISO_C_BINDING
 IMPORT :: t_ARParam
 CHARACTER(C_CHAR) :: text(*)
 TYPE(t_ARParam) :: xparam
 INTEGER(C_INT), VALUE :: i
 END
!
 SUBROUTINE arParamChangeSize(xparam, i, j, yparam) BIND(C, NAME = 'arParamChangeSize')
 USE, INTRINSIC :: ISO_C_BINDING
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam, yparam
 INTEGER(C_INT), VALUE :: i, j
 END
! 
 SUBROUTINE arInitCparam(xparam) BIND(C, NAME = 'arInitCparam' )
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
! 
 SUBROUTINE arParamDisp(xparam) BIND(C, NAME = 'arParamDisp')
 IMPORT :: t_ARParam
 TYPE(t_ARParam) :: xparam
 END
!
 INTEGER(C_INT) FUNCTION arLoadPatt(text) BIND(C, NAME = 'arLoadPatt')
 USE, INTRINSIC :: ISO_C_BINDING
 CHARACTER(C_CHAR) :: text(*)
 END
!
 SUBROUTINE argInit(iptr, d, i1, i2, i3, i4) BIND(C, NAME = 'argInit')
 USE, INTRINSIC :: ISO_C_BINDING
 TYPE(C_PTR), VALUE :: iptr
 REAL(C_DOUBLE)      , VALUE :: d
 INTEGER(C_INT)     , VALUE :: i1, i2, i3, i4
 END
!
 REAL(C_DOUBLE) FUNCTION arUtilTimer() BIND(C, NAME = 'arUtilTimer')
 USE, INTRINSIC :: ISO_C_BINDING
 END FUNCTION arUtilTimer
END INTERFACE
CONTAINS
!-------------------------------------------------
SUBROUTINE init() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
TYPE(t_ARParam) :: wparam
CHARACTER (LEN = 35) :: vconf = "Data/WDM_camera_flipV.xml" // c_null_char
CHARACTER (LEN = 35) :: cparam_name = "Data/camera_para.dat" // c_null_char
CHARACTER (LEN = 35) :: patt_name   = "Data/patt.hiro"       // c_null_char
!   /* open the video path */
IF( arVideoOpen( vconf ) < 0 ) STOP
!    /* find the size of the window */
IF( arVideoInqSize( ixsize, iysize ) < 0 ) STOP
PRINT *, "Image size (x,y) = ", ixsize, iysize
!    /* set the initial camera parameters */
IF ( arParamLoad( cparam_name, 1, wparam ) < 0 ) THEN
 PRINT *, "Camera parameter load error !!"
 STOP
END IF
CALL arParamChangeSize( wparam, ixsize, iysize, cparam )
CALL arInitCparam( cparam )
PRINT *, "*** Camera Parameter ***"
CALL arParamDisp( cparam )
ipatt_id = arLoadPatt( patt_name )
IF ( ipatt_id < 0 ) THEN
 PRINT *, "pattern load error !!"
 STOP
END IF
!    /* open the graphics window */
CALL argInit( C_LOC(cparam), 1.0_c_double, 0, 0, 0, 0 )
RETURN
END SUBROUTINE init
!-------------------------------------------------
SUBROUTINE mainLoop() BIND(C, NAME = 'mainLoop')
IMPLICIT NONE
INTEGER(C_INTPTR_T) :: dataPtr
TYPE(t_ARMarkerInfo), POINTER:: marker_info(:)
INTEGER(C_INT) :: marker_num
INTEGER(C_INT) :: j, k
LOGICAL, SAVE :: qfirst = .TRUE.
IF (qfirst) THEN 
 ALLOCATE(marker_info(50)) ! ?
 qfirst = .FALSE.
END IF
dataPtr = arVideoGetImage()
IF (dataPtr == 0) THEN
 CALL arUtilSleep(2)
 RETURN
END IF
IF (icount == 0) CALL arUtilTimerReset()
icount = icount + 1
CALL argDrawMode2D()
CALL argDispImage( dataPtr, 0, 0 )
IF ( arDetectMarker( dataPtr, ithresh, marker_info, marker_num ) < 0 ) THEN
 CALL cleanup()
 RETURN
END IF
CALL arVideoCapNext()
k = -1
DO j = 1, marker_num
 IF (ipatt_id == marker_info(j)%id) THEN
  IF (k == -1) THEN 
   k = j
  ELSE 
   IF (marker_info(k)%cf < marker_info(j)%cf ) k = j
  END IF
 END IF
END DO
IF ( k == -1 ) THEN 
 CALL argSwapBuffers()
 RETURN
END IF
! draw-main 
CALL arGetTransMat( C_LOC(marker_info(k)), patt_center, patt_width, patt_trans )
CALL draw()
CALL argSwapBuffers()
RETURN
END SUBROUTINE mainLoop
!-------------------------------------------------
SUBROUTINE cleanup() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE 
CALL arVideoCapStop()
CALL arVideoClose()
CALL argCleanup()
RETURN
END SUBROUTINE cleanup
!-------------------------------------------------
SUBROUTINE keyEvent( key, ix, iy ) BIND(C, NAME = 'keyEvent')
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
INTEGER(C_SIGNED_CHAR), VALUE :: key
INTEGER(C_INT)        , VALUE :: ix, iy
IF (key == 27) THEN
  PRINT *, "*** ", icount / arUtilTimer(), "(frame/sec)"
  CALL cleanup()
  STOP
END IF
RETURN
END SUBROUTINE keyEvent
!---------------------------------------------------------------
SUBROUTINE draw() BIND(C)
USE, INTRINSIC :: ISO_C_BINDING 
USE opengl_gl
USE opengl_glu
USE opengl_glut
IMPLICIT NONE
REAL(gldouble) :: gl_para(4, 4)
REAL(glfloat ) :: mat_ambient(4) = [ 0.0_glfloat, 0.0_glfloat, 1.0_glfloat, 1.0_glfloat ]
REAL(glfloat ) :: mat_flash(4)   = [ 0.0_glfloat, 0.0_glfloat, 1.0_glfloat, 1.0_glfloat ]
REAL(glfloat ) :: mat_flash_shiny(1) = [ 50.0 ]
REAL(glfloat ) :: ambi(4)            = [ 0.1_glfloat, 0.1_glfloat, 0.1_glfloat, 0.1_glfloat ]
REAL(glfloat ) :: lightZeroColor(4)  = [ 1.0_glfloat, 0.0_glfloat, 0.0_glfloat, 0.1_glfloat ]
REAL(glfloat ), PARAMETER :: a = 10.0_glfloat, b = 28.0_glfloat, &
                              c =  8.0_glfloat / 3.0_glfloat, d = 0.01_glfloat
! shadow  
REAL(gldouble) :: shadowMatrix(4, 4), plane(4)
REAL(glfloat ), SAVE :: xl = 50.0_glfloat, yl = 50.0_glfloat, zl = 200.0_glfloat, wl = 0.0_glfloat, theta = 0.0
!
! Lorenz
REAL(glfloat ) :: dx, dy, dz
REAL(glfloat ), SAVE :: x(3000), y(3000), z(3000)
INTEGER :: i
LOGICAL, SAVE :: first = .TRUE.


xl = 100.0 * COS(theta)
yl = 100.0 * SIN(theta)
theta = theta + 0.05
!
! Lorenz attractor
IF (first) THEN
 first = .FALSE.
 x(1) = 1.0_glfloat
 y(1) = 1.0_glfloat
 z(1) = 1.0_glfloat
 DO i = 1, 3000 - 1
  dx = d * ( a    * (y(i) - x(i))     )
  dy = d * ( x(i) * (b - z(i)) - y(i) )
  dz = d * ( x(i) * y(i) - c * z(i)   )
  x(i + 1) = x(i) + dx
  y(i + 1) = y(i) + dy
  z(i + 1) = z(i) + dz
 END DO
END IF
! 
CALL argDrawMode3D()
CALL argDraw3dCamera( 0, 0 )
CALL glClearDepth( 1.0_gldouble )
CALL glClear(GL_DEPTH_BUFFER_BIT)
CALL glEnable(GL_DEPTH_TEST)
CALL glDepthFunc(GL_LEQUAL)
CALL argConvGlpara(patt_trans, gl_para) ! /* load the camera transformation matrix */
CALL glMatrixMode(GL_MODELVIEW)
CALL glLoadMatrixd( gl_para )
CALL glLightfv(GL_LIGHT0, GL_POSITION, [xl, yl, zl, wl]) ! light_position + light kind
CALL glLightfv(GL_LIGHT0, GL_AMBIENT , ambi)
CALL glLightfv(GL_LIGHT0, GL_DIFFUSE , lightZeroColor)
CALL glMaterialfv(GL_FRONT, GL_SPECULAR , mat_flash)
CALL glMaterialfv(GL_FRONT, GL_SHININESS, mat_flash_shiny)	
CALL glMaterialfv(GL_FRONT, GL_AMBIENT  , mat_ambient)
CALL glMatrixMode(GL_MODELVIEW)
!
! light   position = [xl, yl, zl, wl]
!
CALL glPushMatrix()
CALL glPushAttrib(GL_LIGHTING_BIT)
CALL glColor3f(1.0_glfloat, 0.0_glfloat, 0.1_glfloat)
CALL glTranslatef( xl, yl, zl )
CALL glutSolidSphere(5.0_gldouble, 10, 10)
CALL glPopAttrib()
CALL glPopMatrix()
!
! Shadow
!
plane = [ 0.0_gldouble,   0.0_gldouble,  1.0_gldouble,  0.0_gldouble]     ! plane normal; 4th element = intercept
CALL glPushMatrix()
CALL projectShadowMatrix(shadowMatrix, plane, [xl, yl, zl, wl]) ! light vector; 4th element = scaling factor
CALL glMultMatrixd(shadowMatrix)
CALL glColor3f(0.1_glfloat, 0.1_glfloat, 0.1_glfloat)
!
CALL glLineWidth(2.0_glfloat)
CALL glBegin(GL_LINE_STRIP)
CALL glBegin(GL_POLYGON)
DO i = 1, 3000
 IF (i > 100) THEN
  CALL glVertex3f(2 * x(i), 2 * y(i), 2 * z(i))
 END IF 
END DO
CALL glEnd()
CALL glPopMatrix()
!
! Object
!
CALL glPushMatrix()
CALL glEnable(GL_LIGHTING)
CALL glEnable(GL_LIGHT0)
!
CALL glLineWidth(2.0_glfloat)
CALL glBegin(GL_LINE_STRIP)
CALL glBegin(GL_POLYGON)
DO i = 1, 3000
 IF (i > 100) THEN
  CALL glVertex3f(2 * x(i), 2 * y(i), 2 * z(i))
 END IF 
END DO
CALL glEnd()
CALL glPopMatrix()
!
CALL glDisable( GL_LIGHTING )
CALL glDisable( GL_DEPTH_TEST )
!
RETURN
END SUBROUTINE draw
!-------------------------------------------------------------
SUBROUTINE projectShadowMatrix(sm, plane, light)
REAL(gldouble), INTENT(OUT) :: sm(:, :)
REAL(gldouble), INTENT(IN ) :: plane(:)
REAL(glfloat ), INTENT(IN ) :: light(:)
REAL(gldouble) :: dot
INTEGER :: i
sm = 0.0
dot = DOT_PRODUCT(plane, REAL(light, gldouble))     ! <plane|light> I  dot_product * diagonal Matrix
DO i = 1, 4
 sm(:, i) = -REAL(light(:), gldouble) * plane(i)    ! |light><plane|   4x4 Matrix
END DO
DO i = 1, 4
 sm(i, i) = dot + sm(i, i)                          
END DO
RETURN
END SUBROUTINE projectShadowMatrix
!-------------------------------------------------------------
END MODULE m_ART
!=====================================================
PROGRAM shadow
USE m_art
IMPLICIT NONE
CALL glutInit()
CALL init()
CALL arVideoCapStart()
CALL argMainLoop( 0, C_FUNLOC(keyEvent), C_FUNLOC(mainLoop) )
STOP
END PROGRAM shadow