~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina6/tina-libs/tina/math/mathGeom_mat2.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**********
  2  *
  3  * Copyright (c) 2003, Division of Imaging Science and Biomedical Engineering,
  4  * University of Manchester, UK.  All rights reserved.
  5  * 
  6  * Redistribution and use in source and binary forms, with or without modification, 
  7  * are permitted provided that the following conditions are met:
  8  * 
  9  *   . Redistributions of source code must retain the above copyright notice, 
 10  *     this list of conditions and the following disclaimer.
 11  *    
 12  *   . Redistributions in binary form must reproduce the above copyright notice,
 13  *     this list of conditions and the following disclaimer in the documentation 
 14  *     and/or other materials provided with the distribution.
 15  * 
 16  *   . Neither the name of the University of Manchester nor the names of its
 17  *     contributors may be used to endorse or promote products derived from this 
 18  *     software without specific prior written permission.
 19  * 
 20  * 
 21  * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" 
 22  * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE 
 23  * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE 
 24  * ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE 
 25  * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR 
 26  * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF 
 27  * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
 28  * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN 
 29  * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) 
 30  * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE 
 31  * POSSIBILITY OF SUCH DAMAGE.
 32  *
 33  **********
 34  *
 35  * Program :    TINA
 36  * File    :  $Source: /home/tina/cvs/tina-libs/tina/math/mathGeom_mat2.c,v $
 37  * Date    :  $Date: 2003/09/23 11:32:20 $
 38  * Version :  $Revision: 1.5 $
 39  * CVS Id  :  $Id: mathGeom_mat2.c,v 1.5 2003/09/23 11:32:20 matts Exp $
 40  *
 41  * Author  :  Legacy TINA
 42  *
 43  * Notes : Mat2 handling
 44  *
 45  *********
 46 */
 47 /** 
 48  *  @file
 49  *  @brief  Functions for handling 2 by 2 matrices (Mat2).      
 50  *
 51  *  Four main types of functions;
 52  *
 53  *  - Mat2 allocation and initialisation routines.
 54  *  - Functions for accessing and reordering the values in a Mat2.
 55  *  - Mat2 manipulations (eg addition of 2 Mat2's, transpose of a Mat2 etc).
 56  *  - Mat2 file/print/format routines.
 57  * 
 58 */
 59 
 60 #include "mathGeom_mat2.h"
 61 
 62 #if HAVE_CONFIG_H
 63 #include <config.h>
 64 #endif
 65 
 66 
 67 #include <math.h>
 68 #include <stdio.h>
 69 #include <string.h>
 70 #include <tina/sys/sysDef.h>
 71 #include <tina/sys/sysPro.h>
 72 #include <tina/math/math_GeomDef.h>
 73 #include <tina/math/mathGeom_vec2.h>
 74 
 75 
 76 
 77 void           *mat2_alloc(void)
 78 {
 79         Mat2           *m = ts_ralloc(Mat2);
 80 
 81         m->el[0][0] = 0.0;
 82         m->el[0][1] = 0.0;
 83         m->el[1][0] = 0.0;
 84         m->el[1][1] = 0.0;
 85         return ((void *) m);
 86 }
 87 
 88 void           *mat2_make(Mat2 n)
 89 {
 90         Mat2           *m = ts_ralloc(Mat2);
 91 
 92         *m = n;
 93         return ((void *) m);
 94 }
 95 
 96 void            mat2_free(void *m)
 97 {
 98         rfree((void *) m);
 99 }
100 
101 Mat2            mat2(double mxx, double mxy, double myx, double myy)
102 {
103         Mat2            m = {Mat2_id};
104 
105         m.el[0][0] = mxx;
106         m.el[0][1] = mxy;
107         m.el[1][0] = myx;
108         m.el[1][1] = myy;
109         return (m);
110 }
111 
112 Mat2            mat2_unit(void)
113 {
114         return (mat2(1.0, 0.0,
115                      0.0, 1.0));
116 }
117 
118 Mat2            mat2_zero(void)
119 {
120         static Mat2     mat2_0 = {Mat2_id};             /* static data! */
121 
122         return (mat2_0);
123 }
124 
125 void            mat2_comps(Mat2 m, float *mxx, float *mxy, float *myx, float *myy)
126 {
127         *mxx = m.el[0][0];
128         *mxy = m.el[0][1];
129         *myx = m.el[1][0];
130         *myy = m.el[1][1];
131 }
132 
133 Vec2            mat2_rowx(Mat2 m)
134 {
135         return (vec2(m.el[0][0], m.el[0][1]));
136 }
137 
138 Vec2            mat2_rowy(Mat2 m)
139 {
140         return (vec2(m.el[1][0], m.el[1][1]));
141 }
142 
143 Vec2            mat2_colx(Mat2 m)
144 {
145         return (vec2(m.el[0][0], m.el[1][0]));
146 }
147 
148 Vec2            mat2_coly(Mat2 m)
149 {
150         return (vec2(m.el[0][1], m.el[1][1]));
151 }
152 
153 Mat2            mat2_of_rows(Vec2 rx, Vec2 ry)
154 {
155         Mat2            m = {Mat2_id};
156 
157         m.el[0][0] = rx.el[0];
158         m.el[1][0] = rx.el[1];
159         m.el[0][1] = ry.el[0];
160         m.el[1][1] = ry.el[1];
161         return (m);
162 }
163 
164 Mat2            mat2_of_cols(Vec2 cx, Vec2 cy)
165 {
166         Mat2            m = {Mat2_id};
167 
168         m.el[0][0] = cx.el[0];
169         m.el[0][1] = cx.el[1];
170         m.el[1][0] = cy.el[0];
171         m.el[1][1] = cy.el[1];
172         return (m);
173 }
174 
175 Mat2            mat2_sum(Mat2 m, Mat2 n)
176 {
177         Mat2            sum = {Mat2_id};
178         int             i, j;
179 
180         for (i = 0; i < 2; ++i)
181                 for (j = 0; j < 2; ++j)
182                         sum.el[i][j] = m.el[i][j] + n.el[i][j];
183         return (sum);
184 }
185 
186 Mat2            mat2_diff(Mat2 m, Mat2 n)
187 {
188         Mat2            diff = {Mat2_id};
189         int             i, j;
190 
191         for (i = 0; i < 2; ++i)
192                 for (j = 0; j < 2; ++j)
193                         diff.el[i][j] = m.el[i][j] - n.el[i][j];
194         return (diff);
195 }
196 
197 Mat2            mat2_minus(Mat2 m)
198 {
199         Mat2            minus = {Mat2_id};
200         int             i, j;
201 
202         for (i = 0; i < 2; ++i)
203                 for (j = 0; j < 2; ++j)
204                         minus.el[i][j] = -m.el[i][j];
205         return (minus);
206 }
207 
208 Mat2            mat2_times(double k, Mat2 m)
209 {
210         Mat2            prod = {Mat2_id};
211         int             i, j;
212 
213         for (i = 0; i < 2; ++i)
214                 for (j = 0; j < 2; ++j)
215                         prod.el[i][j] = k * m.el[i][j];
216         return (prod);
217 }
218 
219 Mat2            mat2_prod(Mat2 m, Mat2 n)
220 {
221         Mat2            prod = {Mat2_id};
222         int             i, j, k;
223         double          sum;
224 
225         for (i = 0; i < 2; ++i)
226                 for (j = 0; j < 2; ++j)
227                 {
228                         sum = 0.0;
229                         for (k = 0; k < 2; ++k)
230                                 sum += m.el[i][k] * n.el[k][j];
231                         prod.el[i][j] = sum;
232                 }
233         return (prod);
234 }
235 
236 Mat2            mat2_inverse(Mat2 m)
237 {
238         double          det, k;
239         float           mxx, mxy, myx, myy;
240 
241         mat2_comps(m, &mxx, &mxy, &myx, &myy);
242         det = mxx * myy - mxy * myx;
243         if (det == 0.0)
244                 return (mat2_unit());
245         k = 1.0 / det;
246         return (mat2(k * myy, -k * mxy,
247                      -k * myx, k * mxx));
248 }
249 
250 Mat2            mat2_transpose(Mat2 m)
251 {
252         Mat2            trans = {Mat2_id};
253         int             i, j;
254 
255         for (i = 0; i < 2; ++i)
256                 for (j = 0; j < 2; ++j)
257                         trans.el[i][j] = m.el[j][i];
258         return (trans);
259 }
260 
261 Mat2            mat2_sym(Mat2 m)
262 {
263         double          mxx = mat2_xx(m);
264         double          mxy = mat2_xy(m);
265         double          myx = mat2_yx(m);
266         double          myy = mat2_yy(m);
267 
268         mxy = 0.5 * (mxy + myx);
269         return (mat2(mxx, mxy, mxy, myy));
270 }
271 
272 double          mat2_trace(Mat2 m)
273 {
274         return (m.el[0][0] + m.el[1][1]);
275 }
276 
277 double          mat2_det(Mat2 m)
278 {
279         return (mat2_xx(m) * mat2_yy(m) - mat2_xy(m) * mat2_yx(m));
280 }
281 
282 Bool            mat2_posdef(Mat2 m)
283 {
284         if (mat2_xx(m) <= 0.0)
285                 return (false);
286         else if (mat2_det(m) <= 0.0)
287                 return (false);
288         else
289                 return (true);
290 }
291 
292 Vec2            mat2_vprod(Mat2 m, Vec2 v)
293 {
294         Vec2            prod = {Vec2_id};
295         int             i, j;
296         double          sum;
297 
298         for (i = 0; i < 2; ++i)
299         {
300                 sum = 0.0;
301                 for (j = 0; j < 2; ++j)
302                         sum += m.el[i][j] * v.el[j];
303                 prod.el[i] = sum;
304         }
305         return (prod);
306 }
307 
308 double          mat2_sprod(Vec2 v, Mat2 m, Vec2 w)
309 {
310         double          sum = 0.0;
311         int             i, j;
312 
313         for (i = 0; i < 2; ++i)
314                 for (j = 0; j < 2; ++j)
315                         sum += v.el[i] * m.el[i][j] * w.el[j];
316         return (sum);
317 }
318 
319 Mat2            mat2_tensor(Vec2 v, Vec2 w)
320 {
321         Mat2            prod = {Mat2_id};
322         int             i, j;
323 
324         for (i = 0; i < 2; ++i)
325                 for (j = 0; j < 2; ++j)
326                         prod.el[i][j] = v.el[i] * w.el[j];
327         return (prod);
328 }
329 
330 Mat2            mat2_read(FILE * fp)
331 {
332         Mat2            temp = {Mat2_id};
333         int             i;
334 
335         for (i = 0; i < 2; ++i)
336                 (void) fscanf(fp, "%f %f", &temp.el[i][0], &temp.el[i][1]);
337         return (temp);
338 }
339 
340 void            mat2_print(FILE * fp, Mat2 m)
341 {
342         int             i;
343 
344         for (i = 0; i < 2; ++i)
345                 (void) fprintf(fp, "%f %f\n", m.el[i][0], m.el[i][1]);
346 }
347 
348 void            mat2_pprint(FILE * fp, char *msg, Mat2 m)
349 {
350         int             i, n;
351 
352         n = strlen(msg);
353         (void) fprintf(fp, "%s", msg);
354         (void) fprintf(fp, "|%15.6f%15.6f|\n", m.el[0][0], m.el[0][1]);
355         for (i = 0; i < n; ++i)
356                 (void) fputc(' ', stderr);
357         (void) fprintf(fp, "|%15.6f%15.6f|\n", m.el[1][0], m.el[1][1]);
358 }
359 
360 void            mat2_format(Mat2 m)
361 {
362         int             i;
363 
364         for (i = 0; i < 2; ++i)
365                 format("%f %f \n", m.el[i][0], m.el[i][1]);
366 }
367 
368 /**
369 eigendirection and eigenvalues of symmetric 2by2 matrix
370 |a b|
371 |b c|
372 **/
373 
374 void            mat2_sym_eigen(Mat2 m, double *theta, double *lambda1, double *lambda2)
375 {
376         double          a, b, c, p, q, r;
377 
378         a = mat2_xx(m);
379         b = 0.5 * (mat2_xy(m) + mat2_yx(m));
380         c = mat2_yy(m);
381         p = (a + c) / 2.0;
382         q = (a - c) / 2.0;
383         r = sqrt(q * q + b * b);
384         if (b == 0.0 && q == 0.0)
385                 *theta = 0.0;
386         else
387                 *theta = 0.5 * atan2(b, q);
388         *lambda1 = p + r;
389         *lambda2 = p - r;
390 }
391 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.