summaryrefslogtreecommitdiffstats
path: root/src/libnr/nr-matrix.cpp
diff options
context:
space:
mode:
Diffstat (limited to 'src/libnr/nr-matrix.cpp')
-rw-r--r--src/libnr/nr-matrix.cpp607
1 files changed, 607 insertions, 0 deletions
diff --git a/src/libnr/nr-matrix.cpp b/src/libnr/nr-matrix.cpp
new file mode 100644
index 000000000..24fa2e206
--- /dev/null
+++ b/src/libnr/nr-matrix.cpp
@@ -0,0 +1,607 @@
+#define __NR_MATRIX_C__
+
+/** \file
+ * Various matrix routines. Currently includes some NR::rotate etc. routines too.
+ */
+
+/*
+ * Authors:
+ * Lauris Kaplinski <lauris@kaplinski.com>
+ *
+ * This code is in public domain
+ */
+
+#include "nr-matrix.h"
+
+
+
+/**
+ * Multiply two NRMatrices together, storing the result in d.
+ */
+NRMatrix *
+nr_matrix_multiply(NRMatrix *d, NRMatrix const *m0, NRMatrix const *m1)
+{
+ if (m0) {
+ if (m1) {
+ NR::Coord d0 = m0->c[0] * m1->c[0] + m0->c[1] * m1->c[2];
+ NR::Coord d1 = m0->c[0] * m1->c[1] + m0->c[1] * m1->c[3];
+ NR::Coord d2 = m0->c[2] * m1->c[0] + m0->c[3] * m1->c[2];
+ NR::Coord d3 = m0->c[2] * m1->c[1] + m0->c[3] * m1->c[3];
+ NR::Coord d4 = m0->c[4] * m1->c[0] + m0->c[5] * m1->c[2] + m1->c[4];
+ NR::Coord d5 = m0->c[4] * m1->c[1] + m0->c[5] * m1->c[3] + m1->c[5];
+
+ NR::Coord *dest = d->c;
+ *dest++ = d0;
+ *dest++ = d1;
+ *dest++ = d2;
+ *dest++ = d3;
+ *dest++ = d4;
+ *dest = d5;
+ } else {
+ *d = *m0;
+ }
+ } else {
+ if (m1) {
+ *d = *m1;
+ } else {
+ nr_matrix_set_identity(d);
+ }
+ }
+
+ return d;
+}
+
+
+
+
+/**
+ * Store the inverted value of Matrix m in d
+ */
+NRMatrix *
+nr_matrix_invert(NRMatrix *d, NRMatrix const *m)
+{
+ if (m) {
+ NR::Coord const det = m->c[0] * m->c[3] - m->c[1] * m->c[2];
+ if (!NR_DF_TEST_CLOSE(det, 0.0, NR_EPSILON)) {
+
+ NR::Coord const idet = 1.0 / det;
+ NR::Coord *dest = d->c;
+
+ /* Cache m->c[0] and m->c[4] in case d == m. */
+ NR::Coord const m_c0(m->c[0]);
+ NR::Coord const m_c4(m->c[4]);
+
+ /*0*/ *dest++ = m->c[3] * idet;
+ /*1*/ *dest++ = -m->c[1] * idet;
+ /*2*/ *dest++ = -m->c[2] * idet;
+ /*3*/ *dest++ = m_c0 * idet;
+ /*4*/ *dest++ = -m_c4 * d->c[0] - m->c[5] * d->c[2];
+ /*5*/ *dest = -m_c4 * d->c[1] - m->c[5] * d->c[3];
+
+ } else {
+ nr_matrix_set_identity(d);
+ }
+ } else {
+ nr_matrix_set_identity(d);
+ }
+
+ return d;
+}
+
+
+
+
+
+/**
+ * Set this matrix to a translation of x and y
+ */
+NRMatrix *
+nr_matrix_set_translate(NRMatrix *m, NR::Coord const x, NR::Coord const y)
+{
+ NR::Coord *dest = m->c;
+
+ *dest++ = 1.0; //0
+ *dest++ = 0.0; //1
+ *dest++ = 0.0; //2
+ *dest++ = 1.0; //3
+ *dest++ = x; //4
+ *dest = y; //5
+
+ return m;
+}
+
+
+
+
+
+/**
+ * Set this matrix to a scaling transform in sx and sy
+ */
+NRMatrix *
+nr_matrix_set_scale(NRMatrix *m, NR::Coord const sx, NR::Coord const sy)
+{
+ NR::Coord *dest = m->c;
+
+ *dest++ = sx; //0
+ *dest++ = 0.0; //1
+ *dest++ = 0.0; //2
+ *dest++ = sy; //3
+ *dest++ = 0.0; //4
+ *dest = 0.0; //5
+
+ return m;
+}
+
+
+
+
+
+/**
+ * Set this matrix to a rotating transform of angle 'theta' radians
+ */
+NRMatrix *
+nr_matrix_set_rotate(NRMatrix *m, NR::Coord const theta)
+{
+ NR::Coord const s = sin(theta);
+ NR::Coord const c = cos(theta);
+
+ NR::Coord *dest = m->c;
+
+ *dest++ = c; //0
+ *dest++ = s; //1
+ *dest++ = -s; //2
+ *dest++ = c; //3
+ *dest++ = 0.0; //4
+ *dest = 0.0; //5
+
+ return m;
+}
+
+
+
+
+
+
+
+
+
+/**
+ * Implement NR functions and methods
+ */
+namespace NR {
+
+
+
+
+
+/**
+ * Constructor. Assign to nr if not null, else identity
+ */
+Matrix::Matrix(NRMatrix const *nr)
+{
+ if (nr) {
+ assign(nr->c);
+ } else {
+ set_identity();
+ }
+}
+
+
+
+
+
+/**
+ * Multiply two matrices together
+ */
+Matrix operator*(Matrix const &m0, Matrix const &m1)
+{
+ NR::Coord const d0 = m0[0] * m1[0] + m0[1] * m1[2];
+ NR::Coord const d1 = m0[0] * m1[1] + m0[1] * m1[3];
+ NR::Coord const d2 = m0[2] * m1[0] + m0[3] * m1[2];
+ NR::Coord const d3 = m0[2] * m1[1] + m0[3] * m1[3];
+ NR::Coord const d4 = m0[4] * m1[0] + m0[5] * m1[2] + m1[4];
+ NR::Coord const d5 = m0[4] * m1[1] + m0[5] * m1[3] + m1[5];
+
+ Matrix ret( d0, d1, d2, d3, d4, d5 );
+
+ return ret;
+}
+
+
+
+
+
+/**
+ * Multiply a matrix by another
+ */
+Matrix &Matrix::operator*=(Matrix const &o)
+{
+ *this = *this * o;
+ return *this;
+}
+
+
+
+
+
+/**
+ * Multiply by a scaling matrix
+ */
+Matrix &Matrix::operator*=(scale const &other)
+{
+ /* This loop is massive overkill. Let's unroll.
+ * o _c[] goes from 0..5
+ * o other[] alternates between 0 and 1
+ */
+ /*
+ * for (unsigned i = 0; i < 3; ++i) {
+ * for (unsigned j = 0; j < 2; ++j) {
+ * this->_c[i * 2 + j] *= other[j];
+ * }
+ * }
+ */
+
+ NR::Coord const xscale = other[0];
+ NR::Coord const yscale = other[1];
+ NR::Coord *dest = _c;
+
+ /*i=0 j=0*/ *dest++ *= xscale;
+ /*i=0 j=1*/ *dest++ *= yscale;
+ /*i=1 j=0*/ *dest++ *= xscale;
+ /*i=1 j=1*/ *dest++ *= yscale;
+ /*i=2 j=0*/ *dest++ *= xscale;
+ /*i=2 j=1*/ *dest *= yscale;
+
+ return *this;
+}
+
+
+
+
+
+/**
+ * Return the inverse of this matrix. If an inverse is not defined,
+ * then return the identity matrix.
+ */
+Matrix Matrix::inverse() const
+{
+ Matrix d(0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
+
+ NR::Coord const det = _c[0] * _c[3] - _c[1] * _c[2];
+ if (!NR_DF_TEST_CLOSE(det, 0.0, NR_EPSILON)) {
+
+ NR::Coord const idet = 1.0 / det;
+ NR::Coord *dest = d._c;
+
+ /*0*/ *dest++ = _c[3] * idet;
+ /*1*/ *dest++ = -_c[1] * idet;
+ /*2*/ *dest++ = -_c[2] * idet;
+ /*3*/ *dest++ = _c[0] * idet;
+ /*4*/ *dest++ = -_c[4] * d._c[0] - _c[5] * d._c[2];
+ /*5*/ *dest = -_c[4] * d._c[1] - _c[5] * d._c[3];
+
+ } else {
+ d.set_identity();
+ }
+
+ return d;
+}
+
+
+
+
+
+/**
+ * Set this matrix to Identity
+ */
+void Matrix::set_identity()
+{
+ NR::Coord *dest = _c;
+
+ *dest++ = 1.0; //0
+ *dest++ = 0.0; //1
+ *dest++ = 0.0; //2
+ *dest++ = 1.0; //3
+ // translation
+ *dest++ = 0.0; //4
+ *dest = 0.0; //5
+}
+
+
+
+
+
+/**
+ * return an Identity matrix
+ */
+Matrix identity()
+{
+ Matrix ret(1.0, 0.0,
+ 0.0, 1.0,
+ 0.0, 0.0);
+ return ret;
+}
+
+
+
+
+
+/**
+ *
+ */
+Matrix from_basis(Point const x_basis, Point const y_basis, Point const offset)
+{
+ Matrix const ret(x_basis[X], y_basis[X],
+ x_basis[Y], y_basis[Y],
+ offset[X], offset[Y]);
+ return ret;
+}
+
+
+
+
+/**
+ * Returns a rotation matrix corresponding by the specified angle (in radians) about the origin.
+ *
+ * \see NR::rotate_degrees
+ *
+ * Angle direction in Inkscape code: If you use the traditional mathematics convention that y
+ * increases upwards, then positive angles are anticlockwise as per the mathematics convention. If
+ * you take the common non-mathematical convention that y increases downwards, then positive angles
+ * are clockwise, as is common outside of mathematics.
+ */
+rotate::rotate(NR::Coord const theta) :
+ vec(cos(theta),
+ sin(theta))
+{
+}
+
+
+
+
+
+/**
+ * Return the determinant of the Matrix
+ */
+NR::Coord Matrix::det() const
+{
+ return _c[0] * _c[3] - _c[1] * _c[2];
+}
+
+
+
+
+
+/**
+ * Return the scalar of the descriminant of the Matrix
+ */
+NR::Coord Matrix::descrim2() const
+{
+ return fabs(det());
+}
+
+
+
+
+
+/**
+ * Return the descriminant of the Matrix
+ */
+NR::Coord Matrix::descrim() const
+{
+ return sqrt(descrim2());
+}
+
+
+
+
+
+/**
+ * Assign a matrix to a given coordinate array
+ */
+Matrix &Matrix::assign(Coord const *array)
+{
+ assert(array != NULL);
+
+ Coord const *src = array;
+ Coord *dest = _c;
+
+ *dest++ = *src++; //0
+ *dest++ = *src++; //1
+ *dest++ = *src++; //2
+ *dest++ = *src++; //3
+ *dest++ = *src++; //4
+ *dest = *src ; //5
+
+ return *this;
+}
+
+
+
+
+
+/**
+ * Copy this matrix's value to a NRMatrix
+ */
+NRMatrix *Matrix::copyto(NRMatrix *nrm) const {
+
+ assert(nrm != NULL);
+
+ Coord const *src = _c;
+ Coord *dest = nrm->c;
+
+ *dest++ = *src++; //0
+ *dest++ = *src++; //1
+ *dest++ = *src++; //2
+ *dest++ = *src++; //3
+ *dest++ = *src++; //4
+ *dest = *src ; //5
+
+ return nrm;
+}
+
+
+
+
+/**
+ * Copy this matrix's values to an array
+ */
+NR::Coord *Matrix::copyto(NR::Coord *array) const {
+
+ assert(array != NULL);
+
+ Coord const *src = _c;
+ Coord *dest = array;
+
+ *dest++ = *src++; //0
+ *dest++ = *src++; //1
+ *dest++ = *src++; //2
+ *dest++ = *src++; //3
+ *dest++ = *src++; //4
+ *dest = *src ; //5
+
+ return array;
+}
+
+
+
+
+
+/**
+ *
+ */
+double expansion(Matrix const &m) {
+ return sqrt(fabs(m.det()));
+}
+
+
+
+
+
+/**
+ *
+ */
+double Matrix::expansion() const {
+ return sqrt(fabs(det()));
+}
+
+
+
+
+
+/**
+ *
+ */
+double Matrix::expansionX() const {
+ return sqrt(_c[0] * _c[0] + _c[1] * _c[1]);
+}
+
+
+
+
+
+/**
+ *
+ */
+double Matrix::expansionY() const {
+ return sqrt(_c[2] * _c[2] + _c[3] * _c[3]);
+}
+
+
+
+
+
+/**
+ *
+ */
+bool Matrix::is_translation(Coord const eps) const {
+ return ( fabs(_c[0] - 1.0) < eps &&
+ fabs(_c[3] - 1.0) < eps &&
+ fabs(_c[1]) < eps &&
+ fabs(_c[2]) < eps );
+}
+
+
+
+
+
+/**
+ *
+ */
+bool Matrix::test_identity() const {
+ return NR_MATRIX_DF_TEST_CLOSE(this, &NR_MATRIX_IDENTITY, NR_EPSILON);
+}
+
+
+
+
+
+/**
+ *
+ */
+bool transform_equalp(Matrix const &m0, Matrix const &m1, NR::Coord const epsilon) {
+ return NR_MATRIX_DF_TEST_TRANSFORM_CLOSE(&m0, &m1, epsilon);
+}
+
+
+
+
+
+/**
+ *
+ */
+bool translate_equalp(Matrix const &m0, Matrix const &m1, NR::Coord const epsilon) {
+ return NR_MATRIX_DF_TEST_TRANSLATE_CLOSE(&m0, &m1, epsilon);
+}
+
+
+
+
+
+/**
+ *
+ */
+bool matrix_equalp(Matrix const &m0, Matrix const &m1, NR::Coord const epsilon)
+{
+ return ( NR_MATRIX_DF_TEST_TRANSFORM_CLOSE(&m0, &m1, epsilon) &&
+ NR_MATRIX_DF_TEST_TRANSLATE_CLOSE(&m0, &m1, epsilon) );
+}
+
+
+
+
+
+/**
+ * A home-made assertion. Stop if the two matrixes are not 'close' to
+ * each other.
+ */
+void assert_close(Matrix const &a, Matrix const &b)
+{
+ if (!matrix_equalp(a, b, 1e-3)) {
+ fprintf(stderr,
+ "a = | %g %g |,\tb = | %g %g |\n"
+ " | %g %g | \t | %g %g |\n"
+ " | %g %g | \t | %g %g |\n",
+ a[0], a[1], b[0], b[1],
+ a[2], a[3], b[2], b[3],
+ a[4], a[5], b[4], b[5]);
+ abort();
+ }
+}
+
+
+
+} //namespace NR
+
+
+
+/*
+ Local Variables:
+ mode:c++
+ c-file-style:"stroustrup"
+ c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
+ indent-tabs-mode:nil
+ fill-column:99
+ End:
+*/
+// vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :
generated by cgit v1.2.3 (git 2.25.1) at 2026-05-19 01:57:05 +0000