Home My Page Projects Code Snippets Project Openings 3D graphics for Standard ML
 Summary Activity SCM

# SCM Repository

[sml3d] View of /src/common/matrix4d.sml
 [sml3d] / src / common / matrix4d.sml

# View of /src/common/matrix4d.sml

Tue Jun 17 03:29:03 2008 UTC (11 years, 3 months ago) by jhr
File size: 10520 byte(s)
```  Added rotations.
```
```(* matrix4d.sml
*
* COPYRIGHT (c) 2008 John Reppy (http://cs.uchicago.edu/~jhr)
*
* Double-precision 4x4 matrices
*)

structure Matrix4d : sig

type mat4d

type double = SML3dTypes.double
type vec3d = SML3dTypes.vec3d
type vec4d = SML3dTypes.vec4d

val mat : {
m11 : double, m12 : double, m13 : double, m14 : double,
m21 : double, m22 : double, m23 : double, m24 : double,
m31 : double, m32 : double, m33 : double, m34 : double,
m41 : double, m42 : double, m43 : double, m44 : double
} -> mat4d

(* standard matrices *)
val identity : mat4d			(* identity matrix *)
val isoscale : double -> mat4d		(* isotropic scaling *)
val scale    : vec3d -> mat4d		(* anisotropic scaling *)
val rotateX  : double -> mat4d		(* rotation around X axis (in degrees) *)
val rotateY  : double -> mat4d		(* rotation around Y axis (in degrees) *)
val rotateZ  : double -> mat4d		(* rotation around Z axis (in degrees) *)
val rotate   : double * vec3d -> mat4d	(* rotation around an axis (in degrees) *)

val perspective : {fov : double, aspect : double, near : double, far : double} -> mat4d

val mxv : mat4d * vec4d -> vec4d
val vxm : vec4d * mat4d -> vec4d
val mxm : mat4d * mat4d -> mat4d

val transpose : mat4d -> mat4d

val inverse : mat4d -> mat4d option

(* returns matrix as flat vector in column-major order *)
val toVector : mat4d -> double Vector.vector

end = struct

type double = SML3dTypes.double
type vec3d = SML3dTypes.vec3d
type vec4d = SML3dTypes.vec4d

val epsilond : double = 1.0e~7

(* represented in column-major order, since that is what OpenGL LoadMatrix
* expects.
*)
datatype mat4d = M of double Vector.vector

fun ! (a : double Vector.vector, (r, c)) = Unsafe.Vector.sub(a, 4*c + r - 5)
infix !

fun new {m11, m12, m13, m14, m21, m22, m23, m24, m31, m32, m33, m34, m41, m42, m43, m44} =
M(Vector.fromList [
m11, m21, m31, m41, m12, m22, m32, m42, m13, m23, m33, m43, m14, m24, m34, m44
])

val mat = new

val identity = new{
m11 = 1.0, m12 = 0.0, m13 = 0.0, m14 = 0.0,
m21 = 0.0, m22 = 1.0, m23 = 0.0, m24 = 0.0,
m31 = 0.0, m32 = 0.0, m33 = 1.0, m34 = 0.0,
m41 = 0.0, m42 = 0.0, m43 = 0.0, m44 = 1.0
}

fun isoscale s = new{
m11 = s,   m12 = 0.0, m13 = 0.0, m14 = 0.0,
m21 = 0.0, m22 = s,   m23 = 0.0, m24 = 0.0,
m31 = 0.0, m32 = 0.0, m33 = s,   m34 = 0.0,
m41 = 0.0, m42 = 0.0, m43 = 0.0, m44 = 1.0
}

fun scale {x, y, z} = new{
m11 = x,   m12 = 0.0, m13 = 0.0, m14 = 0.0,
m21 = 0.0, m22 = y,   m23 = 0.0, m24 = 0.0,
m31 = 0.0, m32 = 0.0, m33 = z,   m34 = 0.0,
m41 = 0.0, m42 = 0.0, m43 = 0.0, m44 = 1.0
}

fun translate {x, y, z} = new{
m11 = 1.0, m12 = 0.0, m13 = 0.0, m14 = x,
m21 = 0.0, m22 = 1.0, m23 = 0.0, m24 = y,
m31 = 0.0, m32 = 0.0, m33 = 1.0, m34 = z,
m41 = 0.0, m42 = 0.0, m43 = 0.0, m44 = 1.0
}

(* rotation around X axis (in degrees) *)
fun rotateX angle = let
val s = Double.sin angle
val c = Double.cos angle
in
new {
m11 = 1.0, m12 = 0.0, m13 = 0.0, m14 = 0.0,
m21 = 0.0, m22 = c,   m23 = ~s,	 m24 = 0.0,
m31 = 0.0, m32 = s,   m33 = c,   m34 = 0.0,
m41 = 0.0, m42 = 0.0, m43 = 0.0, m44 = 1.0
}
end

(* rotation around Y axis (in degrees) *)
fun rotateY angle = let
val s = Double.sin angle
val c = Double.cos angle
in
new {
m11 = c,   m12 = 0.0, m13 = s,   m14 = 0.0,
m21 = 0.0, m22 = 1.0, m23 = 0.0, m24 = 0.0,
m31 = ~s,  m32 = 0.0, m33 = c,   m34 = 0.0,
m41 = 0.0, m42 = 0.0, m43 = 0.0, m44 = 1.0
}
end

(* rotation around Z axis (in degrees) *)
fun rotateZ angle = let
val s = Double.sin angle
val c = Double.cos angle
in
new {
m11 = c,   m12 = ~s,  m13 = 0.0, m14 = 0.0,
m21 = s,   m22 = c,   m23 = 0.0, m24 = 0.0,
m31 = 0.0, m32 = 0.0, m33 = 1.0, m34 = 0.0,
m41 = 0.0, m42 = 0.0, m43 = 0.0, m44 = 1.0
}
end

(* rotation around an axis (in degrees) *)
fun rotate (angle, {x, y, z}) = let
val s = Double.sin angle
val c = Double.cos angle
val c' = 1.0 - c
in
new {
m11 = x*x*c'+c,   m12 = x*y*c'-z*s, m13 = x*z*c'+y*s, m14 = 0.0,
m21 = y*x*c'+z*s, m22 = y*y*c'+c,   m23 = y*z*c'-x*s, m24 = 0.0,
m31 = x*z*c'-y*s, m32 = y*z*c'+x*s, m33 = z*z*c'+c,   m34 = 0.0,
m41 = 0.0,        m42 = 0.0,        m43 = 0.0,        m44 = 1.0
}
end

(* construct the OpenGL perspective projection matrix *)
fun perspective {fov, aspect, near, far} = let
val f = 1.0 / Double.tan ((fov * Double.M_PI) / 360.0)
val d' = 1.0 / (near - far)
in
new {
m11 = f / aspect, m12 = 0.0,	m13 = 0.0,	     m14 = 0.0,
m21 = 0.0,	  m22 = f,	m23 = 0.0,	     m24 = 0.0,
m31 = 0.0,	  m32 = 0.0,	m33 = d'*(near+far), m34 = 2.0*d'*near*far,
m41 = 0.0,	  m42 = 0.0,	m43 = ~1.0,	     m44 = 0.0
}
end

(* multiply a matrix times a column vector *)
fun mxv (M mat, {x, y, z, w}) = let
fun m arg = mat ! arg
in {
x = m(1,1) * x + m(1,2) * y + m(1,3) * z + m(1,4) * w,
y = m(2,1) * x + m(2,2) * y + m(2,3) * z + m(2,4) * w,
z = m(3,1) * x + m(3,2) * y + m(3,3) * z + m(3,4) * w,
w = m(4,1) * x + m(4,2) * y + m(4,3) * z + m(4,4) * w
} end

(* multiply a row vector times a matrix *)
fun vxm ({x, y, z, w}, M mat) = let
fun m arg = mat ! arg
in {
x = m(1,1) * x + m(2,1) * y + m(3,1) * z + m(4,1) * w,
y = m(1,2) * x + m(2,2) * y + m(3,2) * z + m(4,2) * w,
z = m(1,3) * x + m(2,3) * y + m(3,3) * z + m(4,3) * w,
w = m(1,4) * x + m(2,4) * y + m(3,4) * z + m(4,4) * w
} end

(* matrix transpose *)
fun transpose (M mat) = let
fun m arg = mat ! arg
in
new {
m11 = m(1,1), m12 = m(2,1), m13 = m(3,1), m14 = m(4,1),
m21 = m(1,2), m22 = m(2,2), m23 = m(3,2), m24 = m(4,2),
m31 = m(1,3), m32 = m(2,3), m33 = m(3,3), m34 = m(4,3),
m41 = m(1,4), m42 = m(2,4), m43 = m(3,4), m44 = m(4,4)
}
end

(* matrix multiplication *)
fun mxm (M aMat, M bMat) = let
fun a arg = aMat ! arg
fun b arg = bMat ! arg
in
new {
m11 = a(1,1) * b(1,1) + a(1,2) * b(2,1) + a(1,3) * b(3,1) + a(1,4) * b(4,1),
m12 = a(1,1) * b(1,2) + a(1,2) * b(2,2) + a(1,3) * b(3,2) + a(1,4) * b(4,2),
m13 = a(1,1) * b(1,3) + a(1,2) * b(2,3) + a(1,3) * b(3,3) + a(1,4) * b(4,3),
m14 = a(1,1) * b(1,4) + a(1,2) * b(2,4) + a(1,3) * b(3,4) + a(1,4) * b(4,4),
m21 = a(2,1) * b(1,1) + a(2,2) * b(2,1) + a(2,3) * b(3,1) + a(2,4) * b(4,1),
m22 = a(2,1) * b(1,2) + a(2,2) * b(2,2) + a(2,3) * b(3,2) + a(2,4) * b(4,2),
m23 = a(2,1) * b(1,3) + a(2,2) * b(2,3) + a(2,3) * b(3,3) + a(2,4) * b(4,3),
m24 = a(2,1) * b(1,4) + a(2,2) * b(2,4) + a(2,3) * b(3,4) + a(2,4) * b(4,4),
m31 = a(3,1) * b(1,1) + a(3,2) * b(2,1) + a(3,3) * b(3,1) + a(3,4) * b(4,1),
m32 = a(3,1) * b(1,2) + a(3,2) * b(2,2) + a(3,3) * b(3,2) + a(3,4) * b(4,2),
m33 = a(3,1) * b(1,3) + a(3,2) * b(2,3) + a(3,3) * b(3,3) + a(3,4) * b(4,3),
m34 = a(3,1) * b(1,4) + a(3,2) * b(2,4) + a(3,3) * b(3,4) + a(3,4) * b(4,4),
m41 = a(4,1) * b(1,1) + a(4,2) * b(2,1) + a(4,3) * b(3,1) + a(4,4) * b(4,1),
m42 = a(4,1) * b(1,2) + a(4,2) * b(2,2) + a(4,3) * b(3,2) + a(4,4) * b(4,2),
m43 = a(4,1) * b(1,3) + a(4,2) * b(2,3) + a(4,3) * b(3,3) + a(4,4) * b(4,3),
m44 = a(4,1) * b(1,4) + a(4,2) * b(2,4) + a(4,3) * b(3,4) + a(4,4) * b(4,4)
}
end

fun inverse (M mat) = let
fun m arg = (mat ! arg)
(* 2x2 subdeterminants *)
val tbt34C12 = m(3,1)*m(4,2) - m(3,2)*m(4,1)
val tbt34C13 = m(3,1)*m(4,3) - m(3,3)*m(4,1)
val tbt34C14 = m(3,1)*m(4,4) - m(3,4)*m(4,1)
val tbt34C23 = m(3,2)*m(4,3) - m(3,3)*m(4,2)
val tbt34C24 = m(3,2)*m(4,4) - m(3,4)*m(4,2)
val tbt34C34 = m(3,3)*m(4,4) - m(3,4)*m(4,3)
(* 2x2 subdeterminants *)
val tbt24C12 = m(2,1)*m(4,2) - m(2,2)*m(4,1)
val tbt24C13 = m(2,1)*m(4,3) - m(2,3)*m(4,1)
val tbt24C14 = m(2,1)*m(4,4) - m(2,4)*m(4,1)
val tbt24C23 = m(2,2)*m(4,3) - m(2,3)*m(4,2)
val tbt24C24 = m(2,2)*m(4,4) - m(2,4)*m(4,2)
val tbt24C34 = m(2,3)*m(4,4) - m(2,4)*m(4,3)
(* 2x2 subdeterminants *)
val tbt23C12 = m(2,1)*m(3,2) - m(2,2)*m(3,1)
val tbt23C13 = m(2,1)*m(3,3) - m(2,3)*m(3,1)
val tbt23C14 = m(2,1)*m(3,4) - m(2,4)*m(3,1)
val tbt23C23 = m(2,2)*m(3,3) - m(2,3)*m(3,2)
val tbt23C24 = m(2,2)*m(3,4) - m(2,4)*m(3,2)
val tbt23C34 = m(2,3)*m(3,4) - m(2,4)*m(3,3)
(* 3x3 subdeterminants *)
val sd11 = m(2,2)*tbt34C34 - m(2,3)*tbt34C24 + m(2,4)*tbt34C23
val sd12 = m(2,1)*tbt34C34 - m(2,3)*tbt34C14 + m(2,4)*tbt34C13
val sd13 = m(2,1)*tbt34C24 - m(2,2)*tbt34C14 + m(2,4)*tbt34C12
val sd14 = m(2,1)*tbt34C23 - m(2,2)*tbt34C13 + m(2,3)*tbt34C12
(* 3x3 subdeterminants *)
val sd21 = m(1,2)*tbt34C34 - m(1,3)*tbt34C24 + m(1,4)*tbt34C23
val sd22 = m(1,1)*tbt34C34 - m(1,3)*tbt34C14 + m(1,4)*tbt34C13
val sd23 = m(1,1)*tbt34C24 - m(1,2)*tbt34C14 + m(1,4)*tbt34C12
val sd24 = m(1,1)*tbt34C23 - m(1,2)*tbt34C13 + m(1,3)*tbt34C12
(* 3x3 subdeterminants *)
val sd31 = m(1,2)*tbt24C34 - m(1,3)*tbt24C24 + m(1,4)*tbt24C23
val sd32 = m(1,1)*tbt24C34 - m(1,3)*tbt24C14 + m(1,4)*tbt24C13
val sd33 = m(1,1)*tbt24C24 - m(1,2)*tbt24C14 + m(1,4)*tbt24C12
val sd34 = m(1,1)*tbt24C23 - m(1,2)*tbt24C13 + m(1,3)*tbt24C12
(* 3x3 subdeterminants *)
val sd41 = m(1,2)*tbt23C34 - m(1,3)*tbt23C24 + m(1,4)*tbt23C23
val sd42 = m(1,1)*tbt23C34 - m(1,3)*tbt23C14 + m(1,4)*tbt23C13
val sd43 = m(1,1)*tbt23C24 - m(1,2)*tbt23C14 + m(1,4)*tbt23C12
val sd44 = m(1,1)*tbt23C23 - m(1,2)*tbt23C13 + m(1,3)*tbt23C12
(* determinant *)
val det = m(1,1)*sd11 - m(1,2)*sd12 + m(1,3)*sd13 - m(1,4)*sd14
in
if ((det > ~epsilond) andalso (det < epsilond))
then NONE
else let
val detInv = 1.0 / det
in
SOME(new {
m11 = sd11*detInv,
m21 = ~sd12*detInv,
m31 = sd13*detInv,
m41 = ~sd14*detInv,
m12 = ~sd21*detInv,
m22 = sd22*detInv,
m32 = ~sd23*detInv,
m42 = sd24*detInv,
m13 = sd31*detInv,
m23 = ~sd32*detInv,
m33 = sd33*detInv,
m43 = ~sd34*detInv,
m14 = ~sd41*detInv,
m24 = sd42*detInv,
m34 = ~sd43*detInv,
m44 = sd44*detInv
})
end
end

fun toVector (M v) = v

end
```