math_linear_algebra.scad

Go to the documentation of this file.

//! Linear algebra mathematical functions.
/***************************************************************************//**
  \file   math_linear_algebra.scad
  \author Roy Allen Sutton
  \date   2015-2017

  \copyright

    This file is part of [omdl] (https://github.com/royasutton/omdl),
    an OpenSCAD mechanical design library.

    The \em omdl is free software; you can redistribute it and/or modify
    it under the terms of the [GNU Lesser General Public License]
    (http://www.gnu.org/licenses/lgpl.html) as published by the Free
    Software Foundation; either version 2.1 of the License, or (at
    your option) any later version.

    The \em omdl is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
    Lesser General Public License for more details.

    You should have received a copy of the GNU Lesser General Public
    License along with the \em omdl; if not, write to the Free Software
    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
    02110-1301, USA; or see <http://www.gnu.org/licenses/>.

*******************************************************************************/

//----------------------------------------------------------------------------//
/***************************************************************************//**
  \addtogroup math
  @{

  \defgroup math_linear_algebra Linear Algebra
  \brief    Linear algebra transformations on Euclidean coordinates.
  @{
*******************************************************************************/
//----------------------------------------------------------------------------//

//! Multiply all coordinates by a 4x4 3d-transformation matrix.
/***************************************************************************//**
  \param    c <coords-3d> A list of 3d coordinate points.
  \param    m <matrix-4x4> A 4x4 transformation matrix (decimal-list-4-list4).

  \returns  <coords-3d> A list of 3d coordinate points multiplied by the
            transformation matrix.

  \details

    See [Wikipedia] and [multmatrix] for more information.

    [Wikipedia]: https://en.wikipedia.org/wiki/Transformation_matrix
    [multmatrix]: https://en.wikibooks.org/wiki/OpenSCAD_User_Manual/Transformations#multmatrix
*******************************************************************************/
function multmatrix_lp
(
  c,
  m
) =
  let
  (
    m11=m[0][0], m12=m[0][1], m13=m[0][2], m14=m[0][3],
    m21=m[1][0], m22=m[1][1], m23=m[1][2], m24=m[1][3],
    m31=m[2][0], m32=m[2][1], m33=m[2][2], m34=m[2][3]
  )
  [
    for (ci=c)
    let
    (
      x = ci[0], y = ci[1], z = ci[2]
    )
      [m11*x+m12*y+m13*z+m14, m21*x+m22*y+m23*z+m24, m31*x+m32*y+m33*z+m34]
  ];

//! Translate all coordinates dimensions.
/***************************************************************************//**
  \param    c <coords-nd> A list of nd coordinate points.
  \param    v <decimal-list-n> A list of translations for each dimension.

  \returns  <coords-nd> A list of translated coordinate points.

  \details

    See [Wikipedia] for more information and [transformation matrix].

    [Wikipedia]: https://en.wikipedia.org/wiki/Translation_(geometry)
    [transformation matrix]: https://en.wikipedia.org/wiki/Transformation_matrix
*******************************************************************************/
function translate_lp
(
  c,
  v
) = not_defined(v) ? c
  : let
    (
      d = len(first(c)),
      u = is_scalar(v) ? v : 0,
      w = [for (i=[0 : d-1]) edefined_or(v, i, u)]
    )
    [for (ci=c) [for (di=[0 : d-1]) ci[di] + w[di]]];

//! Rotate all coordinates about one or more axes in Euclidean 2d or 3d space.
/***************************************************************************//**
  \param    c <coords-3d|coords-2d> A list of 3d or 2d coordinate points.
  \param    a <decimal-list-3|decimal> The axis rotation angle.
            A list [ax, ay, az] or a single decimal to specify az only.
  \param    v <vector-3d> An arbitrary axis for the rotation. When
            specified, the rotation angle will be \p a or az about the
            line \p v that passes through point \p o.
  \param    o <point-3d> A 3d point origin for the rotation.
            Ignored when \p v is not specified.

  \returns  <coords-3d|coords-2d> A list of 3d or 2d rotated coordinates.
            Rotation order is rz, ry, rx.

  \details

    See [Wikipedia] for more information on [transformation matrix]
    and [axis rotation].

    [Wikipedia]: https://en.wikipedia.org/wiki/Rotation_matrix
    [transformation matrix]: https://en.wikipedia.org/wiki/Transformation_matrix
    [axis rotation]: http://inside.mines.edu/fs_home/gmurray/ArbitraryAxisRotation
*******************************************************************************/
function rotate_lp
(
  c,
  a,
  v,
  o = origin3d
) = not_defined(a) ? c
  : let
    (
      d  = len(first(c)),
      az = edefined_or(a, 2, is_scalar(a) ? a : 0),

      cg = cos(az), sg = sin(az),

      rc = (d == 2) ? [for (ci=c) [cg*ci[0]-sg*ci[1], sg*ci[0]+cg*ci[1]]]
         : (d != 3) ? c
         : (not_defined(v) || is_list(a)) ?
            let
            (
              ax = edefined_or(a, 0, 0),
              ay = edefined_or(a, 1, 0),

              ca = cos(ax), cb = cos(ay),
              sa = sin(ax), sb = sin(ay),

              m11 = cb*cg,
              m12 = cg*sa*sb-ca*sg,
              m13 = ca*cg*sb+sa*sg,

              m21 = cb*sg,
              m22 = ca*cg+sa*sb*sg,
              m23 = -cg*sa+ca*sb*sg,

              m31 = -sb,
              m32 = cb*sa,
              m33 = ca*cb
            )
            multmatrix_lp(c, [[m11,m12,m13,0], [m21,m22,m23,0], [m31,m32,m33,0]])
         :  let
            (
              vx  = v[0],  vy  = v[1],  vz  = v[2],
              vx2 = vx*vx, vy2 = vy*vy, vz2 = vz*vz,
              l2  = vx2 + vy2 + vz2
            )
            (l2 == 0) ? c
         :  let
            (
              ox = o[0], oy = o[1], oz = o[2],
              ll = sqrt(l2),
              oc = 1 - cg,

              m11 = vx2+(vy2+vz2)*cg,
              m12 = vx*vy*oc-vz*ll*sg,
              m13 = vx*vz*oc+vy*ll*sg,
              m14 = (ox*(vy2+vz2)-vx*(oy*vy+oz*vz))*oc+(oy*vz-oz*vy)*ll*sg,

              m21 = vx*vy*oc+vz*ll*sg,
              m22 = vy2+(vx2+vz2)*cg,
              m23 = vy*vz*oc-vx*ll*sg,
              m24 = (oy*(vx2+vz2)-vy*(ox*vx+oz*vz))*oc+(oz*vx-ox*vz)*ll*sg,

              m31 = vx*vz*oc-vy*ll*sg,
              m32 = vy*vz*oc+vx*ll*sg,
              m33 = vz2+(vx2+vy2)*cg,
              m34 = (oz*(vx2+vy2)-vz*(ox*vx+oy*vy))*oc+(ox*vy-oy*vx)*ll*sg
            )
            multmatrix_lp(c, [[m11,m12,m13,m14], [m21,m22,m23,m24], [m31,m32,m33,m34]])/l2
    )
    rc;

//! Scale all coordinates dimensions.
/***************************************************************************//**
  \param    c <coords-nd> A list of nd coordinate points.
  \param    v <decimal-list-n> A list of scalers for each dimension.

  \returns  <coords-nd> A list of scaled coordinate points.
*******************************************************************************/
function scale_lp
(
  c,
  v
) = not_defined(v) ? c
  : let
    (
      d = len(first(c)),
      u = is_scalar(v) ? v : 1,
      w = [for (i=[0 : d-1]) edefined_or(v, i, u)]
    )
    [for (ci=c) [for (di=[0 : d-1]) ci[di] * w[di]]];

//! Scale all coordinates dimensions proportionately to fit inside a region.
/***************************************************************************//**
  \param    c <coords-nd> A list of nd coordinate points.
  \param    v <decimal-list-n> A list of bounds for each dimension.

  \returns  <coords-nd> A list of proportionately scaled coordinate
            points which exactly fit the region bounds \p v.
*******************************************************************************/
function resize_lp
(
  c,
  v
) = not_defined(v) ? c
  : let
    (
      d = len(first(c)),
      u = is_scalar(v) ? v : 1,
      w = [for (i=[0 : d-1]) edefined_or(v, i, u)],
      m = [for (i=[0 : d-1]) let (cv = [for (ci=c) (ci[i])]) [min(cv), max(cv)]],
      s = [for (i=[0 : d-1]) second(m[i]) - first(m[i])]
    )
    [for (ci=c) [for (di=[0 : d-1]) ci[di]/s[di] * w[di]]];

//! @}
//! @}

//----------------------------------------------------------------------------//
// end of file
//----------------------------------------------------------------------------//