Adaptation des ent�tes.

1 files changed, 684 insertions, 0 deletions

diff --git a/generic/Transfo.c b/generic/Transfo.c
new file mode 100644
index 0000000..aa50303
--- /dev/null
+++ b/generic/Transfo.c
@@ -0,0 +1,684 @@
+/*
+ * Transfo.c -- Implementation of transformation routines.
+ *
+ * Authors		: Patrick Lecoanet.
+ * Creation date	: 
+ *
+ * $Id$
+ */
+
+/*
+ *  Copyright (c) 1993 - 1999 CENA, Patrick Lecoanet --
+ *
+ * This code is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Library General Public
+ * License as published by the Free Software Foundation; either
+ * version 2 of the License, or (at your option) any later version.
+ *
+ * This code 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
+ * Library General Public License for more details.
+ *
+ * You should have received a copy of the GNU Library General Public
+ * License along with this code; if not, write to the Free
+ * Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+ *
+ */
+
+/*
+ * This package deals with *AFFINE* 3x3 matrices.
+ * This means that you should not try to feed it with matrices
+ * containing perspective changes. It is assumed that the third
+ * column is always [0 0 1] as this is the case for affine matrices.
+ * Furthermore affine matrices are known to be invertible (non singular).
+ * Despite this, various tests are done to test the invertibility because
+ * of numerical precision or limit.
+ * Any of the operations in this module yield an affine matrix. Composition
+ * of two affine matrices and inversion of an affine matrix result in an
+ * affine matrix (Affine matrices Group property). Rotation, translation
+ * anamorphic scaling, xy shear and yx shear also preserve the property.
+ *
+ */
+
+
+#include "Item.h"
+#include "Geo.h"
+#include "Transfo.h"
+#include "Types.h"
+
+#include <stdlib.h>
+
+
+static const char rcsid[] = "$Imagine: Transfo.c,v 1.7 1997/01/24 14:33:37 lecoanet Exp $";
+static const char compile_id[]="$Compile: " __FILE__ " " __DATE__ " " __TIME__ " $";
+
+
+
+/*
+ *************************************************************************
+ *
+ * The transformation primitives are based on affines matrices retricted
+ * to the following pattern:
+ *
+ *	x x 0
+ *	x x 0
+ *	x x 1
+ *
+ * It is necessary to only feed such matrices to these primitives as they
+ * do optimizations based on their properties. Furthermore the package
+ * stores only the first two columns, the third is constant. There is no
+ * way to describe perspective transformation with these transformation
+ * matrices.
+ *
+ *************************************************************************
+ */
+
+/*
+ *************************************************************************
+ *
+ * RadarTransfoNew --
+ *	Create a new transformation and return it initialized to
+ *	identity.
+ * 
+ *************************************************************************
+ */
+RadarTransfo *
+RadarTransfoNew()
+{
+  RadarTransfo	*t;
+
+  t = (RadarTransfo *) malloc(sizeof(RadarTransfo));
+  RadarTransfoSetIdentity(t);
+  
+  return t;
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarTransfoDuplicate --
+ *	Create a new transformation identical to the model t.
+ * 
+ *************************************************************************
+ */
+RadarTransfo *
+RadarTransfoDuplicate(RadarTransfo *t)
+{
+  RadarTransfo	*nt;
+
+  nt = (RadarTransfo *) RadarMalloc(sizeof(RadarTransfo));
+  if (t) {
+    *nt = *t;
+  }
+  else {
+    RadarTransfoSetIdentity(nt);
+  }
+  
+  return nt;
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarTransfoFree --
+ *	Delete a transformation and free its memory.
+ * 
+ *************************************************************************
+ */
+void
+RadarTransfoFree(RadarTransfo	*t)
+{
+  RadarFree(t);
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarPrintTransfo --
+ *	Print the transfo matrix on stdout.
+ * 
+ *************************************************************************
+ */
+void
+RadarPrintTransfo(RadarTransfo	*t)
+{
+  printf("(%5g %5g\n %5g %5g\n %5g %5g)\n",
+	 t->_[0][0], t->_[0][1],
+	 t->_[1][0], t->_[1][1],
+	 t->_[2][0], t->_[2][1]);
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarTransfoIsIdentity --
+ *	Tell if the given transfo is (close to) identity.
+ * 
+ *************************************************************************
+ */
+RadarBool
+RadarTransfoIsIdentity(RadarTransfo	*t)
+{
+  RadarReal	tmp;
+  RadarBool	res = False;
+
+  tmp = t->_[0][0] - 1.0;
+  res = res & ((tmp < PRECISION_LIMIT) && (tmp > -PRECISION_LIMIT));
+  tmp = t->_[1][1] - 1.0;
+  res = res & ((tmp < PRECISION_LIMIT) && (tmp > -PRECISION_LIMIT));
+  tmp = t->_[0][1];
+  res = res & ((tmp < PRECISION_LIMIT) && (tmp > -PRECISION_LIMIT));
+  tmp = t->_[1][0];
+  res = res & ((tmp < PRECISION_LIMIT) && (tmp > -PRECISION_LIMIT));
+  tmp = t->_[2][0];
+  res = res & ((tmp < PRECISION_LIMIT) && (tmp > -PRECISION_LIMIT));
+  tmp = t->_[2][1];
+  res = res & ((tmp < PRECISION_LIMIT) && (tmp > -PRECISION_LIMIT));
+  return res;
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarTransfoSetIdentity --
+ *	Initialize the given transfo to identity.
+ * 
+ *************************************************************************
+ */
+void
+RadarTransfoSetIdentity(RadarTransfo	*t)
+{
+  *t = ((RadarTransfo) {{{1, 0}, {0, 1}, {0, 0}}});
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarTransfoCompose --
+ *	Combine two transformations t1 and t2 by post-concatenation.
+ *	Returns the resulting transformation.
+ *	t2 can be NULL, meaning identity transform. This is used in
+ *	the toolkit to optimize some cases.
+ *
+ *	All the parameters must be distincts transforms.
+ *
+ *************************************************************************
+ */
+RadarTransfo *
+RadarTransfoCompose(RadarTransfo	*res,
+		    RadarTransfo	*t1,
+		    RadarTransfo	*t2)
+{
+  if ((t1 != NULL) && (t2 != NULL)) {
+    register RadarReal	tmp;
+
+    tmp = t1->_[0][0];
+    res->_[0][0] = tmp*t2->_[0][0] + t1->_[0][1]*t2->_[1][0];
+    res->_[0][1] = tmp*t2->_[0][1] + t1->_[0][1]*t2->_[1][1];
+    tmp = t1->_[1][0];
+    res->_[1][0] = tmp*t2->_[0][0] + t1->_[1][1]*t2->_[1][0];
+    res->_[1][1] = tmp*t2->_[0][1] + t1->_[1][1]*t2->_[1][1];
+    tmp = t1->_[2][0];
+    res->_[2][0] = tmp*t2->_[0][0] + t1->_[2][1]*t2->_[1][0] + t2->_[2][0];
+    res->_[2][1] = tmp*t2->_[0][1] + t1->_[2][1]*t2->_[1][1] + t2->_[2][1];
+  }
+  else if (t1 == NULL) {
+    if (res != t2) {
+      *res = *t2;
+    }
+  }
+  else if (t2 == NULL) {
+    if (res != t2) {
+      *res = *t1;
+    }
+  }
+  else {
+    RadarTransfoSetIdentity(res);
+  }
+  
+  return res;
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarTransfoInvert --
+ *	Compute the inverse of the given matrix and return it. This
+ *	function makes the assumption that the matrix is affine to
+ *	optimize the job. Do not give it a general matrix, this will
+ *	fail. This code is from Graphics Gems II. Anyway an affine
+ *	matrix is always invertible for affine matrices form a sub
+ *	group of the non-singular matrices.
+ *
+ *************************************************************************
+ */
+RadarTransfo *
+RadarTransfoInvert(RadarTransfo	*t,
+		   RadarTransfo	*inv)
+{
+  RadarReal	pos, neg, temp, det_l;
+
+  if (t == NULL) {
+    return NULL;
+  }
+
+  /*
+   * Compute the determinant of the upper left 2x2 sub matrix to see
+   * if it is singular.
+   */
+  pos = neg = 0.0;
+  temp = t->_[0][0] * t->_[1][1];
+  if (temp >= 0.0) {
+    pos += temp;
+  }
+  else {
+    neg += temp;
+  }
+  temp = - t->_[0][1] * t->_[1][0];
+  if (temp >= 0.0) {
+    pos += temp;
+  }
+  else {
+    neg += temp;
+  }
+  det_l = pos + neg;
+  temp = det_l / (pos - neg); /* Why divide by (pos - neg) ?? */
+  
+  if (ABS(temp) < PRECISION_LIMIT) {
+    RadarWarning("RadarTransfoInvert : singular matrix");
+    return NULL;
+  }
+  
+  det_l = 1.0/ det_l;
+  inv->_[0][0] = t->_[1][1] * det_l;
+  inv->_[0][1] = - t->_[0][1] * det_l;
+  inv->_[1][0] = - t->_[1][0] * det_l;
+  inv->_[1][1] = t->_[0][0] * det_l;
+  /*
+   * The code below is equivalent to:
+   *   inv->_[2][0] = (t->_[1][0] * t->_[2][1] - t->_[1][1] * t->_[2][0]) * det_l;
+   *   inv->_[2][1] = - (t->_[0][0] * t->_[2][1] - t->_[0][1] * t->_[2][0]) * det_l;
+   *
+   * with some operations factored (already computed) to increase speed.
+   */
+  inv->_[2][0] = - (inv->_[0][0] * t->_[2][0] + inv->_[1][0] * t->_[2][1]);
+  inv->_[2][1] = - (inv->_[0][1] * t->_[2][0] + inv->_[1][1] * t->_[2][1]);
+
+  return inv;
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarTransfoDecompose --
+ *	Decompose an affine matrix into translation, scale, shear and
+ *	rotation. The different values are stored in the locations
+ *	pointed to by the pointer parameters. If some values are not
+ *	needed a NULL pointer can be given instead. The resulting shear
+ *	shears x coordinate when y change.
+ *	This code is taken from Graphics Gems II.
+ *
+ *************************************************************************
+ */
+void
+RadarTransfoDecompose(RadarTransfo	*t,
+		      RadarPoint	*scale,
+		      RadarPoint	*trans,
+		      RadarReal		*rotation,
+		      RadarReal		*shearxy)
+{
+  RadarTransfo	local;
+  RadarReal	shear, len, rot, det;
+  
+  if (t == NULL) {
+    /* Identity transform */
+    if (scale) {
+      scale->x = 1.0;
+      scale->y = 1.0;
+    }
+    if (trans) {
+      trans->x = 0.0;
+      trans->y = 0.0;
+    }
+    if (rotation) {
+      *rotation = 0.0;
+    }
+    if (shearxy) {
+      *shearxy = 0.0;
+    }
+    //printf("Transfo is identity\n");
+    return;
+  }
+
+  det = (t->_[0][0]*t->_[1][1] - t->_[0][1]*t->_[1][0]);
+  if (ABS(det) < PRECISION_LIMIT) {
+    RadarWarning("RadarTransfoDecompose : singular matrix");
+    return;
+  }
+  
+  local = *t;
+  //RadarPrintTransfo(&local);
+  /* Get translation part if needed */
+  if (trans) {
+    trans->x = ABS(local._[2][0]) < PRECISION_LIMIT ? 0 : local._[2][0];
+    trans->y = ABS(local._[2][1]) < PRECISION_LIMIT ? 0 : local._[2][1];
+  }
+  if (!scale && !shearxy && !rotation) {
+    return;
+  }
+
+  /* Get scale and shear */
+  len = sqrt(local._[0][0]*local._[0][0] +
+	     local._[0][1]*local._[0][1]); /* Get x scale from 1st row */
+  if (scale) {
+    scale->x = len < PRECISION_LIMIT ? 0.0 : len;
+  }
+  local._[0][0] /= len;			 /* Normalize 1st row */
+  local._[0][1] /= len;
+  shear = (local._[0][0]*local._[1][0] +
+	   local._[0][1]*local._[1][1]); /* Shear is dot product of 1st row & 2nd row */
+  /* Make the 2nd row orthogonal to the 1st row
+   * by linear combinaison:
+   * row1.x = row1.x + row0.x*-shear &
+   * row1.y = row1.y + row0.y*-shear
+   */
+  local._[1][0] -= local._[0][0]*shear;
+  local._[1][1] -= local._[0][1]*shear;
+  len = sqrt(local._[1][0]*local._[1][0] +
+	     local._[1][1]*local._[1][1]); /* Get y scale from 2nd row */
+  if (scale)
+    scale->y = len < PRECISION_LIMIT ? 0.0 : len;
+  local._[1][0] /= len;			 /* Normalize 2nd row */
+  local._[1][1] /= len;
+  shear /= len;
+  if (shearxy) {
+    *shearxy = ABS(shear) < PRECISION_LIMIT ? 0.0 : shear;
+    //printf("shear %f\n", *shearxy);
+  }
+  //printf("Matrix after scale & shear extracted\n");
+  //RadarPrintTransfo(&local);
+  /* Get rotation */
+  /* Check for a coordinate system flip. If det of upper-left 2x2
+   * is -1, there is a reflection. If the rotation is < 180° negate
+   * the y scale. If the rotation is > 180° then negate the x scale
+   * and report a rotation between 0 and 180°. This dissymetry is
+   * the result of computing (z) rotation from the first row (x component
+   * of the axis system basis).
+   */
+  det = (local._[0][0]*local._[1][1]- local._[0][1]*local._[1][0]);
+  
+  rot = atan2(local._[0][1], local._[0][0]);
+  if (rot < 0) {
+    rot = 2*M_PI+rot;
+  }
+  rot = rot < PRECISION_LIMIT ? 0.0 : rot;
+  if (rot >= M_PI) {
+    /*rot -= M_PI;  Why that, I'll have to check Graphic Gems ??? */
+    if (scale && det < 0) {
+      scale->x *= -1;
+    }
+  }
+  else if (scale && det < 0) {
+    scale->y *= -1;
+  }
+  
+  //printf("scalex %f\n", scale->x);
+  //printf("scaley %f\n", scale->y);
+  //printf("rotation %f\n", rot*180.0/3.1415);
+
+  if (rotation) {
+    *rotation = rot;
+  }
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarTransfoEqual --
+ *	Return True if t1 and t2 are equal (i.e they have the same
+ *	rotation, shear scales and translations). If include_translation
+ *	is True the translations are considered in the test.
+ *
+ *************************************************************************
+ */
+RadarBool
+RadarTransfoEqual(RadarTransfo	*t1,
+		  RadarTransfo	*t2,
+		  RadarBool	include_translation)
+{
+  if (include_translation) {
+    return (t1->_[0][0] == t2->_[0][0] &&
+	    t1->_[0][1] == t2->_[0][1] &&
+	    t1->_[1][0] == t2->_[1][0] &&
+	    t1->_[1][1] == t2->_[1][1] &&
+	    t1->_[2][0] == t2->_[2][0] &&
+	    t1->_[2][1] == t2->_[2][1]);
+  }
+  else {
+    return (t1->_[0][0] == t2->_[0][0] &&
+	    t1->_[0][1] == t2->_[0][1] &&
+	    t1->_[1][0] == t2->_[1][0] &&
+	    t1->_[1][1] == t2->_[1][1]);
+  }
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarTransfoHasShear --
+ *	Return True if t has a shear factor in x or y or describe a
+ *	rotation or both.
+ *
+ *************************************************************************
+ */
+RadarBool
+RadarTransfoHasShear(RadarTransfo	*t)
+{
+  return t->_[0][1] != 0.0 || t->_[1][0] != 0.0;
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarTransfoIsTranslation --
+ *	Return True if t is a pure translation.
+ *
+ *************************************************************************
+ */
+RadarBool
+RadarTransfoIsTranslation(RadarTransfo	*t)
+{
+  return (t->_[0][0] == 0.0 &&
+	  t->_[0][1] == 0.0 &&
+	  t->_[1][0] == 0.0 &&
+	  t->_[1][1] == 0.0);
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarTransformPoint --
+ *	Apply the transformation to the point. The point is
+ *	modified and returned as the value of the function.
+ *	A NULL transformation means identity. This is only used
+ *	in the toolkit to optimize some cases. It should never
+ *	happen in user code.
+ *
+ *************************************************************************
+ */
+RadarPoint *
+RadarTransformPoint(RadarTransfo	*t,
+		    register RadarPoint	*p,
+		    RadarPoint		*xp)
+{
+  if (t == NULL) {
+    xp->x = p->x;
+    xp->y = p->y;
+  }
+  else {
+    xp->x = t->_[0][0]*p->x + t->_[1][0]*p->y + t->_[2][0];
+    xp->y = t->_[0][1]*p->x + t->_[1][1]*p->y + t->_[2][1];
+  }
+  return xp;
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarTransformPoints --
+ *	Apply the transformation to the points in p returning points in xp.
+ *	The number of points is in num.
+ *	A NULL transformation means identity. This is only used
+ *	in the toolkit to optimize some cases. It should never
+ *	happen in user code.
+ *
+ *************************************************************************
+ */
+void
+RadarTransformPoints(RadarTransfo	*t,
+		     RadarPoint		*p,
+		     RadarPoint		*xp,
+		     int		num)
+{
+  if (t == NULL) {
+    memcpy(xp, p, sizeof(RadarPoint)*num);
+  }
+  else {
+    int i;
+    
+    for (i = 0; i < num; i++) {
+      xp[i].x = t->_[0][0]*p[i].x + t->_[1][0]*p[i].y + t->_[2][0];
+      xp[i].y = t->_[0][1]*p[i].x + t->_[1][1]*p[i].y + t->_[2][1];
+    }
+  }
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarTranslate --
+ *	Translate the given transformation by delta_x, delta_y. Returns
+ *	the resulting transformation.
+ *
+ * RadarSetTranslation --
+ *	Set the translation instead of combining it into the
+ *	transformation.
+ *
+ *************************************************************************
+ */
+RadarTransfo *
+RadarTranslate(RadarTransfo	*t,
+	       RadarReal	delta_x,
+	       RadarReal	delta_y)
+{
+  t->_[2][0] = t->_[2][0] + delta_x;
+  t->_[2][1] = t->_[2][1] + delta_y;
+
+  return t;
+}
+
+RadarTransfo *
+RadarSetTranslation(RadarTransfo	*t,
+		    RadarReal		delta_x,
+		    RadarReal		delta_y)
+{
+  t->_[2][0] = delta_x;
+  t->_[2][1] = delta_y;
+
+  return t;
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarScale --
+ *	Scale the given transformation by scale_x, scale_y. Returns the
+ *	resulting transformation.
+ *
+ *************************************************************************
+ */
+RadarTransfo *
+RadarScale(RadarTransfo		*t,
+	   register RadarReal	scale_x,
+	   register RadarReal	scale_y)
+{
+  t->_[0][0] = t->_[0][0]*scale_x;
+  t->_[0][1] = t->_[0][1]*scale_y;
+  t->_[1][0] = t->_[1][0]*scale_x;
+  t->_[1][1] = t->_[1][1]*scale_y;
+  t->_[2][0] = t->_[2][0]*scale_x;
+  t->_[2][1] = t->_[2][1]*scale_y;
+  
+  return t;
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarRotateRad --
+ *	Rotate the given transformation by angle radians
+ *	counter-clockwise around the origin. Returns the resulting
+ *	transformation.
+ *
+ *************************************************************************
+ */
+RadarTransfo *
+RadarRotateRad(RadarTransfo	*t,
+	       RadarReal	angle)
+{
+  register RadarReal	c = cos(angle);
+  register RadarReal	s = sin(angle);
+  register RadarReal	tmp;
+
+  tmp = t->_[0][0];
+  t->_[0][0] = tmp*c - t->_[0][1]*s;
+  t->_[0][1] = tmp*s + t->_[0][1]*c;
+  tmp = t->_[1][0];
+  t->_[1][0] = tmp*c - t->_[1][1]*s;
+  t->_[1][1] = tmp*s + t->_[1][1]*c;
+  tmp = t->_[2][0];
+  t->_[2][0] = tmp*c - t->_[2][1]*s;
+  t->_[2][1] = tmp*s + t->_[2][1]*c;
+  
+  return t;
+}
+
+
+/*
+ *************************************************************************
+ *
+ * RadarRotateDeg --
+ *	Rotate the given transformation by angle degrees
+ *	counter-clockwise around the origin. Returns the resulting
+ *	transformation.
+ *
+ *************************************************************************
+ */
+RadarTransfo *
+RadarRotateDeg(RadarTransfo	*t,
+	       RadarReal	angle)
+{
+  return RadarRotateRad(t, DegreesToRadian(angle));
+}
+
+  
+#undef PRECISION_LIMIT
+
+