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

SCM Repository

[sml3d] View of /trunk/sml3d/src/sml3d/base/matrix4f.sml
ViewVC logotype

View of /trunk/sml3d/src/sml3d/base/matrix4f.sml

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1261 - (download) (annotate)
Wed Jan 18 02:21:03 2012 UTC (6 years, 9 months ago) by jhr
File size: 15401 byte(s)
  Working on new SML3d core libraries
(* matrix4f.sml
 *
 * COPYRIGHT (c) 2012 The SML3d Project (http://sml3d.cs.uchicago.edu)
 * All rights reserved.
 *
 * Single-precision 4x4 matrices
 *)

structure Matrix4f :> MATRIX4 where type flt = SML3dTypes.float
  = struct

    type flt = SML3dTypes.float
    type vec3 = SML3dTypes.vec3f
    type vec4 = SML3dTypes.vec4f

    val epsilond = Double.epsilon

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

  (* index = 4*(c-1) + r - 1 = 4*c + r - 5 *)
    val usub = Unsafe.Vector.sub
    fun ! (a : flt Vector.vector, (r, c)) = usub(a, 4*c + r - 5)
    infix 8 !

    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

  (* convert 16-element vector to matrix in column-major order *)
    fun fromVector v = if (Vector.length v = 16) then M v else raise Size
  (* convert 16-element vector to matrix in row-major order *)
    fun fromVectorT v = if (Vector.length v = 16)
	  then M(Vector.fromList [
	      v!(1,1), v!(1,2), v!(1,3), v!(1,4),	(* column 1 *)
	      v!(2,1), v!(2,2), v!(2,3), v!(2,4),	(* column 2 *)
	      v!(3,1), v!(3,2), v!(3,3), v!(3,4),	(* column 3 *)
	      v!(4,1), v!(4,2), v!(4,3), v!(4,4)	(* column 4 *)
	    ])
	  else raise Size

  (* returns matrix as flat vector in column-major order *)
    fun toVector (M v) = v
  (* returns matrix as flat vector in row-major order *)
    fun toVectorT (M v) = Vector.fromList [
	    v!(1,1), v!(1,2), v!(1,3), v!(1,4),	(* column 1 *)
	    v!(2,1), v!(2,2), v!(2,3), v!(2,4),	(* column 2 *)
	    v!(3,1), v!(3,2), v!(3,3), v!(3,4),	(* column 3 *)
	    v!(4,1), v!(4,2), v!(4,3), v!(4,4)	(* column 4 *)
	  ]

  (* create a matrix from three column vectors *)
    fun fromCols (c1 : vec4, c2 : vec4, c3 : vec4, c4 : vec4) = new {
	    m11 = #1 c1, m12 = #1 c2, m13 = #1 c3, m14 = #1 c4,
	    m21 = #2 c1, m22 = #2 c2, m23 = #2 c3, m24 = #2 c4,
	    m31 = #3 c1, m32 = #3 c2, m33 = #3 c3, m34 = #3 c4,
	    m41 = #3 c1, m42 = #3 c2, m43 = #3 c3, m44 = #3 c4
	  }
  (* create a matrix from three row vectors *)
    fun fromRows (r1 : vec4, r2 : vec4, r3 : vec4, r4 : vec4) = new {
	    m11 = #1 r1, m12 = #2 r1, m13 = #3 r1, m14 = #4 r1,
	    m21 = #1 r2, m22 = #2 r2, m23 = #3 r2, m24 = #4 r2,
	    m31 = #1 r3, m32 = #2 r3, m33 = #3 r3, m34 = #4 r3,
	    m41 = #1 r4, m42 = #2 r4, m43 = #3 r4, m44 = #4 r4
	  }

  (* return the columns of the matrix *)
    fun toCols (M v) = (
	    (v!(1,1), v!(2,1), v!(3,1), v!(4,1)),	(* column 1 *)
	    (v!(1,2), v!(2,2), v!(3,2), v!(4,2)),	(* column 2 *)
	    (v!(1,3), v!(2,3), v!(3,3), v!(4,3)),	(* column 3 *)
	    (v!(1,4), v!(2,4), v!(3,4), v!(4,4))	(* column 4 *)
	  )
  (* return the rows of the matrix *)
    fun toRows (M v) = (
	    (v!(1,1), v!(1,2), v!(1,3), v!(1,4)),	(* row 1 *)
	    (v!(2,1), v!(2,2), v!(2,3), v!(2,4)),	(* row 2 *)
	    (v!(3,1), v!(3,2), v!(3,3), v!(3,4)),	(* row 3 *)
	    (v!(4,1), v!(4,2), v!(4,3), v!(4,4))	(* row 4 *)
	  )

  (* project specific columns/rows *)
    fun col1 (M v) = (v!(1,1), v!(2,1), v!(3,1), v!(4,1))
    fun col2 (M v) = (v!(1,2), v!(2,2), v!(3,2), v!(4,2))
    fun col3 (M v) = (v!(1,3), v!(2,3), v!(3,3), v!(4,3))
    fun col4 (M v) = (v!(1,4), v!(2,4), v!(3,4), v!(4,4))
    fun row1 (M v) = (v!(1,1), v!(1,2), v!(1,3), v!(1,4))
    fun row2 (M v) = (v!(2,1), v!(2,2), v!(2,3), v!(2,4))
    fun row3 (M v) = (v!(3,1), v!(3,2), v!(3,3), v!(3,4))
    fun row4 (M v) = (v!(4,1), v!(4,2), v!(4,3), v!(4,4))

    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 angle = Float.toRadians angle
	  val s = Float.sin angle
	  val c = Float.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 angle = Float.toRadians angle
	  val s = Float.sin angle
	  val c = Float.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 angle = Float.toRadians angle
	  val s = Float.sin angle
	  val c = Float.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 angle = Float.toRadians angle
	  val s = Float.sin angle
	  val c = Float.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 / Float.tan ((fov * Float.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

  (* arithmetic *)
    fun add (M v1, M v2) =
	  M(Vector.tabulate(16, fn i => usub(v1, i) + usub(v2, i)))
    fun adds (M v1, s, M v2) =
	  M(Vector.tabulate(16, fn i => usub(v1, i) + s*usub(v2, i)))
    fun sub (M v1, M v2) =
	  M(Vector.tabulate(16, fn i => usub(v1, i) - usub(v2, i)))
    fun negate (M v) = M(Vector.map ~ v)

  (* multiply a scalar times a matrix *)
    fun sxm (s, M v) = M(Vector.map (fn x => s*x) v)

  (* multiply a matrix times a column vector *)
    fun mxv (M mat, (x, y, z, w)) = let
	  fun m arg = mat ! arg
	  in (
	    m(1,1) * x + m(1,2) * y + m(1,3) * z + m(1,4) * w,
	    m(2,1) * x + m(2,2) * y + m(2,3) * z + m(2,4) * w,
	    m(3,1) * x + m(3,2) * y + m(3,3) * z + m(3,4) * 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 (
	    m(1,1) * x + m(2,1) * y + m(3,1) * z + m(4,1) * w,
	    m(1,2) * x + m(2,2) * y + m(3,2) * z + m(4,2) * w,
	    m(1,3) * x + m(2,3) * y + m(3,3) * z + m(4,3) * w,
	    m(1,4) * x + m(2,4) * y + m(3,4) * z + m(4,4) * w
	  ) 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

  (* matrix times point (i.e., w = 0) *)
    fun mxpt (M mat, (x, y, z)) = let
	  fun m arg = mat ! arg
	  in (
	    m(1,1) * x + m(1,2) * y + m(1,3) * z,
	    m(2,1) * x + m(2,2) * y + m(2,3) * z,
	    m(3,1) * x + m(3,2) * y + m(3,3) * z
	  ) end

  (* matrix times point (i.e., w = 1) *)
    fun mxvec (M mat, (x, y, z)) = let
	  fun m arg = mat ! arg
	  in (
	    m(1,1) * x + m(1,2) * y + m(1,3) * z + m(4,1),
	    m(2,1) * x + m(2,2) * y + m(2,3) * z + m(4,2),
	    m(3,1) * x + m(3,2) * y + m(3,3) * z + m(4,3)
	  ) 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

    fun trace (M mat) = mat!(1,1) + mat!(2,2) + mat!(3,3) + mat!(4,4)

    fun det (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
	  in
	    m(1,1)*sd11 - m(1,2)*sd12 + m(1,3)*sd13 - m(1,4)*sd14
	  end

    fun inverse (M mat) = let
	  fun m arg = FP.ftod(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 = FP.dtof(sd11*detInv),
		      m21 = FP.dtof(~sd12*detInv),
		      m31 = FP.dtof(sd13*detInv),
		      m41 = FP.dtof(~sd14*detInv),
		      m12 = FP.dtof(~sd21*detInv),
		      m22 = FP.dtof(sd22*detInv),
		      m32 = FP.dtof(~sd23*detInv),
		      m42 = FP.dtof(sd24*detInv),
		      m13 = FP.dtof(sd31*detInv),
		      m23 = FP.dtof(~sd32*detInv),
		      m33 = FP.dtof(sd33*detInv),
		      m43 = FP.dtof(~sd34*detInv),
		      m14 = FP.dtof(~sd41*detInv),
		      m24 = FP.dtof(sd42*detInv),
		      m34 = FP.dtof(~sd43*detInv),
		      m44 = FP.dtof(sd44*detInv)
		    })
		end
	  end

  end

root@smlnj-gforge.cs.uchicago.edu
ViewVC Help
Powered by ViewVC 1.0.0