1 #define SEEN_TRANSF_MAT_3x4_C
3 /*
4 * 3x4 transformation matrix to map points from projective 3-space into the projective plane
5 *
6 * Authors:
7 * Maximilian Albert <Anhalter42@gmx.de>
8 *
9 * Copyright (C) 2007 Authors
10 *
11 * Released under GNU GPL, read the file 'COPYING' for more information
12 */
14 #include "transf_mat_3x4.h"
15 #include <gtk/gtk.h>
16 #include <2geom/matrix.h>
17 #include "svg/stringstream.h"
18 #include "syseq.h"
19 #include "document.h"
20 #include "inkscape.h"
22 namespace Proj {
24 TransfMat3x4::TransfMat3x4 () {
25 for (unsigned int i = 0; i < 3; ++i) {
26 for (unsigned int j = 0; j < 4; ++j) {
27 tmat[i][j] = (i == j ? 1 : 0); // or should we initialize all values with 0? does it matter at all?
28 }
29 }
30 }
32 TransfMat3x4::TransfMat3x4 (Proj::Pt2 vp_x, Proj::Pt2 vp_y, Proj::Pt2 vp_z, Proj::Pt2 origin) {
33 for (unsigned int i = 0; i < 3; ++i) {
34 tmat[i][0] = vp_x[i];
35 tmat[i][1] = vp_y[i];
36 tmat[i][2] = vp_z[i];
37 tmat[i][3] = origin[i];
38 }
39 }
41 TransfMat3x4::TransfMat3x4(TransfMat3x4 const &rhs) {
42 for (unsigned int i = 0; i < 3; ++i) {
43 for (unsigned int j = 0; j < 4; ++j) {
44 tmat[i][j] = rhs.tmat[i][j];
45 }
46 }
47 }
49 Pt2
50 TransfMat3x4::column (Proj::Axis axis) const {
51 return Proj::Pt2 (tmat[0][axis], tmat[1][axis], tmat[2][axis]);
52 }
54 Pt2
55 TransfMat3x4::image (Pt3 const &point) {
56 double x = tmat[0][0] * point[0] + tmat[0][1] * point[1] + tmat[0][2] * point[2] + tmat[0][3] * point[3];
57 double y = tmat[1][0] * point[0] + tmat[1][1] * point[1] + tmat[1][2] * point[2] + tmat[1][3] * point[3];
58 double w = tmat[2][0] * point[0] + tmat[2][1] * point[1] + tmat[2][2] * point[2] + tmat[2][3] * point[3];
60 return Pt2 (x, y, w);
61 }
63 Pt3
64 TransfMat3x4::preimage (Geom::Point const &pt, double coord, Proj::Axis axis) {
65 double x[4];
66 double v[3];
67 v[0] = pt[Geom::X];
68 v[1] = pt[Geom::Y];
69 v[2] = 1.0;
70 int index = (int) axis;
72 SysEq::SolutionKind sol = SysEq::gaussjord_solve<3,4>(tmat, x, v, index, coord, true);
74 if (sol != SysEq::unique) {
75 if (sol == SysEq::no_solution) {
76 g_print ("No solution. Please investigate.\n");
77 } else {
78 g_print ("Infinitely many solutions. Please investigate.\n");
79 }
80 }
81 return Pt3(x[0], x[1], x[2], x[3]);
82 }
84 void
85 TransfMat3x4::set_image_pt (Proj::Axis axis, Proj::Pt2 const &pt) {
86 // FIXME: Do we need to adapt the coordinates in any way or can we just use them as they are?
87 for (int i = 0; i < 3; ++i) {
88 tmat[i][axis] = pt[i];
89 }
90 }
92 void
93 TransfMat3x4::toggle_finite (Proj::Axis axis) {
94 g_return_if_fail (axis != Proj::W);
95 if (has_finite_image(axis)) {
96 Geom::Point dir (column(axis).affine());
97 Geom::Point origin (column(Proj::W).affine());
98 dir -= origin;
99 set_column (axis, Proj::Pt2(dir[Geom::X], dir[Geom::Y], 0));
100 } else {
101 Proj::Pt2 dir (column(axis));
102 Proj::Pt2 origin (column(Proj::W).affine());
103 dir = dir + origin;
104 dir[2] = 1.0;
105 set_column (axis, dir);
106 }
107 }
109 gchar *
110 TransfMat3x4::pt_to_str (Proj::Axis axis) {
111 Inkscape::SVGOStringStream os;
112 os << tmat[0][axis] << " : "
113 << tmat[1][axis] << " : "
114 << tmat[2][axis];
115 return g_strdup(os.str().c_str());
116 }
118 /* Check for equality (with a small tolerance epsilon) */
119 bool
120 TransfMat3x4::operator==(const TransfMat3x4 &rhs) const
121 {
122 // Should we allow a certain tolerance or "normalize" the matrices first?
123 for (int i = 0; i < 3; ++i) {
124 Proj::Pt2 pt1 = column(Proj::axes[i]);
125 Proj::Pt2 pt2 = rhs.column(Proj::axes[i]);
126 if (pt1 != pt2) {
127 return false;
128 }
129 }
130 return true;
131 }
133 /* Multiply a projective matrix by an affine matrix (by only multiplying the 'affine part' of the
134 * projective matrix) */
135 TransfMat3x4
136 TransfMat3x4::operator*(Geom::Matrix const &A) const {
137 TransfMat3x4 ret;
139 for (int j = 0; j < 4; ++j) {
140 ret.tmat[0][j] = A[0]*tmat[0][j] + A[2]*tmat[1][j] + A[4]*tmat[2][j];
141 ret.tmat[1][j] = A[1]*tmat[0][j] + A[3]*tmat[1][j] + A[5]*tmat[2][j];
142 ret.tmat[2][j] = tmat[2][j];
143 }
145 return ret;
146 }
148 // FIXME: Shouldn't rather operator* call operator*= for efficiency? (Because in operator*=
149 // there is in principle no need to create a temporary object, which happens in the assignment)
150 TransfMat3x4 &
151 TransfMat3x4::operator*=(Geom::Matrix const &A) {
152 *this = *this * A;
153 return *this;
154 }
156 void
157 TransfMat3x4::copy_tmat(double rhs[3][4]) {
158 for (int i = 0; i < 3; ++i) {
159 for (int j = 0; j < 4; ++j) {
160 rhs[i][j] = tmat[i][j];
161 }
162 }
163 }
165 void
166 TransfMat3x4::print () const {
167 g_print ("Transformation matrix:\n");
168 for (int i = 0; i < 3; ++i) {
169 g_print (" ");
170 for (int j = 0; j < 4; ++j) {
171 g_print ("%8.2f ", tmat[i][j]);
172 }
173 g_print ("\n");
174 }
175 }
177 void
178 TransfMat3x4::normalize_column (Proj::Axis axis) {
179 Proj::Pt2 new_col(column(axis));
180 new_col.normalize();
181 set_image_pt(axis, new_col);
182 }
185 } // namespace Proj
187 /*
188 Local Variables:
189 mode:c++
190 c-file-style:"stroustrup"
191 c-file-offsets:((innamespace . 0)(inline-open . 0)(case-label . +))
192 indent-tabs-mode:nil
193 fill-column:99
194 End:
195 */
196 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=99 :