1 /* This file is part of the GNU plotutils package. Copyright (C) 1995,
2 1996, 1997, 1998, 1999, 2000, 2005, 2008, Free Software Foundation, Inc.
3
4 The GNU plotutils package is free software. You may redistribute it
5 and/or modify it under the terms of the GNU General Public License as
6 published by the Free Software foundation; either version 2, or (at your
7 option) any later version.
8
9 The GNU plotutils package is distributed in the hope that it will be
10 useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 General Public License for more details.
13
14 You should have received a copy of the GNU General Public License along
15 with the GNU plotutils package; see the file COPYING. If not, write to
16 the Free Software Foundation, Inc., 51 Franklin St., Fifth Floor,
17 Boston, MA 02110-1301, USA. */
18
19 #include "sys-defines.h"
20 #include "extern.h"
21
22 /* Computes the product of two PS-style transformation matrices
23 (i.e. matrix representations of affine transformations). */
24
25 void
_matrix_product(const double m[6],const double n[6],double product[6])26 _matrix_product (const double m[6], const double n[6], double product[6])
27 {
28 double local_product[6];
29
30 local_product[0] = m[0] * n[0] + m[1] * n[2];
31 local_product[1] = m[0] * n[1] + m[1] * n[3];
32
33 local_product[2] = m[2] * n[0] + m[3] * n[2];
34 local_product[3] = m[2] * n[1] + m[3] * n[3];
35
36 local_product[4] = m[4] * n[0] + m[5] * n[2] + n[4];
37 local_product[5] = m[4] * n[1] + m[5] * n[3] + n[5];
38
39 memcpy (product, local_product, 6 * sizeof (double));
40
41 return;
42 }
43
44 /* Computes the inverse of a PS-style transformation matrix, which should
45 be nonsingular for correct results. */
46
47 void
_matrix_inverse(const double m[6],double inverse[6])48 _matrix_inverse (const double m[6], double inverse[6])
49 {
50 double det = m[0] * m[3] - m[1] * m[2];
51
52 if (det == 0.0)
53 /* bogus */
54 inverse[0] = inverse[1] = inverse[2] = inverse[3]
55 = inverse[4] = inverse[5] = 0.0;
56 else
57 {
58 double invdet = 1.0 / det;
59
60 inverse[0] = invdet * m[3];
61 inverse[1] = - invdet * m[1];
62 inverse[2] = - invdet * m[2];
63 inverse[3] = invdet * m[0];
64 inverse[4] = invdet * (m[2] * m[5] - m[3] * m[4]);
65 inverse[5] = invdet * (m[1] * m[4] - m[0] * m[5]);
66 }
67 }
68
69 /* _matrix_norm computes the matrix norm (in the l^2 sense) of the linear
70 transformation part of a PS-style transformation matrix. Actually we
71 compute instead the geometric mean of the l^1 and l^infinity norms. By
72 Hadamard's 3-line theorem, this geometric mean is an upper bound on the
73 true l^2 norm.
74
75 This function is called only to select appropriate line widths and font
76 sizes. For the purposes of those functions, the above approximation
77 should suffice. */
78
79 double
_matrix_norm(const double m[6])80 _matrix_norm (const double m[6])
81 {
82 double mt[4], pm[4];
83 double norm1, norm2;
84 double a[4];
85 int i;
86
87 mt[0] = m[0]; /* transpose of m */
88 mt[1] = m[2];
89 mt[2] = m[1];
90 mt[3] = m[3];
91
92 pm[0] = m[0] * mt[0] + m[1] * mt[2]; /* pm = m * mt */
93 pm[1] = m[0] * mt[1] + m[1] * mt[3];
94 pm[2] = m[2] * mt[0] + m[3] * mt[2];
95 pm[3] = m[2] * mt[1] + m[3] * mt[3];
96
97 for (i = 0; i < 4; i++)
98 a[i] = fabs(pm[i]);
99
100 /* compute l^1 and l^infinity norms of m * mt */
101 norm1 = DMAX(a[0]+a[1], a[2]+a[3]);
102 norm2 = DMAX(a[0]+a[2], a[1]+a[3]);
103
104 /* l^2 norm of m is sqrt of l^2 norm of m * mt */
105 return sqrt(sqrt(norm1 * norm2));
106 }
107
108 /* Compute the minimum and maximum singular values of a 2-by-2 matrix M.
109 The singular values are the square roots of the eigenvalues of M times
110 its transpose. */
111
112 void
_matrix_sing_vals(const double m[6],double * min_sing_val,double * max_sing_val)113 _matrix_sing_vals (const double m[6], double *min_sing_val, double *max_sing_val)
114 {
115 double mt[4], pm[4];
116 double trace, det, disc, sqrtdisc;
117 double s1, s2;
118
119 mt[0] = m[0]; /* transpose of m */
120 mt[1] = m[2];
121 mt[2] = m[1];
122 mt[3] = m[3];
123
124 pm[0] = m[0] * mt[0] + m[1] * mt[2]; /* pm = m * mt */
125 pm[1] = m[0] * mt[1] + m[1] * mt[3];
126 pm[2] = m[2] * mt[0] + m[3] * mt[2];
127 pm[3] = m[2] * mt[1] + m[3] * mt[3];
128
129 trace = pm[0] + pm[3];
130 det = pm[0] * pm[3] - pm[1] * pm[2];
131 /* s^2 + b s + c = 0, where b = -trace, c = det */
132 disc = trace * trace - 4.0 * det;
133 disc = DMAX(0.0, disc); /* paranoia */
134 sqrtdisc = sqrt (disc);
135 s1 = 0.5 * (trace - sqrtdisc);
136 s2 = 0.5 * (trace + sqrtdisc);
137 s1 = DMAX(0.0, s1); /* paranoia */
138 s2 = DMAX(0.0, s2); /* paranoia */
139
140 *min_sing_val = sqrt(s1);
141 *max_sing_val = sqrt(s2);
142 }
143