Index: /trunk/psLib/src/astro/psCoord.c
===================================================================
--- /trunk/psLib/src/astro/psCoord.c	(revision 6408)
+++ /trunk/psLib/src/astro/psCoord.c	(revision 6409)
@@ -10,6 +10,6 @@
 *  @author GLG, MHPCC
 *
-*  @version $Revision: 1.115 $ $Name: not supported by cvs2svn $
-*  @date $Date: 2006-02-09 01:57:37 $
+*  @version $Revision: 1.116 $ $Name: not supported by cvs2svn $
+*  @date $Date: 2006-02-10 00:51:00 $
 *
 *  Copyright 2004-2005 Maui High Performance Computing Center, University of Hawaii
@@ -306,18 +306,6 @@
         out = psPlaneAlloc();
     }
-    out->x = psPolynomial4DEval(
-                 distort->x,
-                 coords->x,
-                 coords->y,
-                 mag,
-                 color
-             );
-    out->y = psPolynomial4DEval(
-                 distort->y,
-                 coords->x,
-                 coords->y,
-                 mag,
-                 color
-             );
+    out->x = psPolynomial4DEval(distort->x, coords->x, coords->y, mag, color);
+    out->y = psPolynomial4DEval(distort->y, coords->x, coords->y, mag, color);
     return (out);
 }
@@ -547,6 +535,4 @@
 multiplies them.  Basically, for each non-zero coeff in the trans1 coeff[][]
 array, you must multiply by all non-zero coeffs in trans2.
- 
-XXX: Inefficient in that the out polynomial is allocated every time.
  *****************************************************************************/
 static psPolynomial2D *multiplyDPoly2D(
@@ -579,9 +565,187 @@
 
 
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
 /*****************************************************************************
 psPlaneTransformCombine(out, trans1, trans2)
- 
-XXX: Much room for optimization.  Currently, we call the polyMultiply
-routine far too many times.
  *****************************************************************************/
 psPlaneTransform *psPlaneTransformCombine(
@@ -656,43 +820,74 @@
     // and its coefficients are added into the myPT coeff matrix.
     //
-    // XXX: This is horribly inefficient in that the trans1 polys are repeatedly
-    // multiplied against themselves.  This can easily be improved.
-    //
-
-    //
-    // Determine the new x-polynomial
-    //
+    // trans1XPolys[i]: contains a polynomial corresponding to trans1->x raised to the i-th power.
+    //
+
+    psS32 order = PS_MAX(PS_MAX(PS_MAX(trans2->x->nX, trans2->x->nY), trans2->y->nX), trans2->y->nY);
+    psPolynomial2D **trans1XPolys = (psPolynomial2D **) psAlloc((order + 1) * sizeof(psPolynomial2D *));
+    psPolynomial2D **trans1YPolys = (psPolynomial2D **) psAlloc((order + 1) * sizeof(psPolynomial2D *));
+
+    //
+    // Raise the trans1 polynomials to whatever power is need in the trans2 polynomials.
+    //
+    trans1XPolys[0] = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 0, 0);
+    trans1XPolys[0]->coeff[0][0] = 1.0;
+    trans1YPolys[0] = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 0, 0);
+    trans1YPolys[0]->coeff[0][0] = 1.0;
+
+    for (psS32 c = 1 ; c < (order + 1) ; c++) {
+        trans1XPolys[c] = multiplyDPoly2D(trans1XPolys[c-1], trans1->x);
+        trans1YPolys[c] = multiplyDPoly2D(trans1YPolys[c-1], trans1->y);
+    }
+
     psTrace(__func__, 5, "Determine the new x-polynomial\n");
-    for (psS32 t2x = 0 ; t2x < (1 + trans2->x->nX) ; t2x++) {
-        for (psS32 t2y = 0 ; t2y < (1 + trans2->x->nY) ; t2y++) {
+    for (psS32 t2x = 0 ; t2x < (trans2->x->nX + 1) ; t2x++) {
+        for (psS32 t2y = 0 ; t2y < (trans2->x->nY + 1) ; t2y++) {
             psTrace(__func__, 6, "X: -------------------- (t2x, t2y) (%d, %d) --------------------\n", t2x, t2y);
             if (trans2->x->mask[t2x][t2y] == 0) {
                 psTrace(__func__, 6, "In this iteration, we raise trans1->x to the %d power and trans1->y to the %d-power.\n", t2x, t2y);
-                // Create the constant "f(x, y) = 1" polynomial.
-                psPolynomial2D *currPoly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 0, 0);
-                currPoly->coeff[0][0] = 1.0;
-
-                psPolynomial2D *newPoly = NULL;
-                // Must raise trans1->x to the t2x-power.
-                for (psS32 c = 0 ; c < t2x; c++) {
-                    newPoly = multiplyDPoly2D(currPoly, trans1->x);
-                    psFree(currPoly);
-                    currPoly = newPoly;
+                psPolynomial2D *newPoly = multiplyDPoly2D(trans1XPolys[t2x], trans1YPolys[t2y]);
+
+                if (psTraceGetLevel(__func__) >= 6) {
+                    PS_POLY_PRINT_2D(newPoly);
                 }
 
-                // Must raise trans1->y to the t2y-power.
-                for (psS32 c = 0 ; c < t2y; c++) {
-                    newPoly = multiplyDPoly2D(currPoly, trans1->y);
-                    psFree(currPoly);
-                    currPoly = newPoly;
+                // Set the appropriate coeffs in myPT->x
+                for (psS32 i = 0 ; i < (1 + newPoly->nX) ; i++) {
+                    for (psS32 j = 0 ; j < (1 + newPoly->nY) ; j++) {
+                        myPT->x->coeff[i][j]+= newPoly->coeff[i][j] * trans2->x->coeff[t2x][t2y];
+                    }
                 }
+
                 if (psTraceGetLevel(__func__) >= 6) {
-                    PS_POLY_PRINT_2D(currPoly);
+                    PS_POLY_PRINT_2D(myPT->x);
                 }
+                psFree(newPoly);
+            }
+        }
+    }
+    if (psTraceGetLevel(__func__) >= 6) {
+        psTrace(__func__, 6, "The final x-polynomial\n");
+        PS_POLY_PRINT_2D(myPT->x);
+    }
+
+    //
+    // Determine the new y-polynomial
+    //
+    psTrace(__func__, 5, "Determine the new y-polynomial\n");
+    for (psS32 t2x = 0 ; t2x < (trans2->y->nX + 1) ; t2x++) {
+        for (psS32 t2y = 0 ; t2y < (trans2->y->nY + 1) ; t2y++) {
+            psTrace(__func__, 5, "Y: -------------------- (t2x, t2y) (%d, %d) --------------------\n", t2x, t2y);
+            if (trans2->y->mask[t2x][t2y] == 0) {
+                psTrace(__func__, 5, "In this iteration, we raise trans1->x to the %d power and trans1->y to the %d-power.\n", t2x, t2y);
+                psPolynomial2D *newPoly = multiplyDPoly2D(trans1XPolys[t2x], trans1YPolys[t2y]);
+
+                if (psTraceGetLevel(__func__) >= 6) {
+                    PS_POLY_PRINT_2D(newPoly);
+                }
 
                 // Set the appropriate coeffs in myPT->x
-                for (psS32 i = 0 ; i < (1 + currPoly->nX) ; i++) {
-                    for (psS32 j = 0 ; j < (1 + currPoly->nY) ; j++) {
-                        myPT->x->coeff[i][j]+= currPoly->coeff[i][j] * trans2->x->coeff[t2x][t2y];
+                for (psS32 i = 0 ; i < (1 + newPoly->nX) ; i++) {
+                    for (psS32 j = 0 ; j < (1 + newPoly->nY) ; j++) {
+                        myPT->y->coeff[i][j]+= newPoly->coeff[i][j] * trans2->y->coeff[t2x][t2y];
                     }
                 }
@@ -700,54 +895,5 @@
                     PS_POLY_PRINT_2D(myPT->x);
                 }
-                psFree(currPoly);
-            }
-        }
-    }
-    if (psTraceGetLevel(__func__) >= 6) {
-        psTrace(__func__, 6, "The final x-polynomial\n");
-        PS_POLY_PRINT_2D(myPT->x);
-    }
-
-    //
-    // Determine the new y-polynomial
-    //
-    psTrace(__func__, 5, "Determine the new y-polynomial\n");
-    for (psS32 t2x = 0 ; t2x < (1 + trans2->y->nX) ; t2x++) {
-        for (psS32 t2y = 0 ; t2y < (1 + trans2->y->nY) ; t2y++) {
-            psTrace(__func__, 5, "Y: -------------------- (t2x, t2y) (%d, %d) --------------------\n", t2x, t2y);
-            if (trans2->y->mask[t2x][t2y] == 0) {
-                psTrace(__func__, 5, "In this iteration, we raise trans1->x to the %d power and trans1->y to the %d-power.\n", t2x, t2y);
-                // Create the constant "f(x, y) = 1" polynomial.
-                psPolynomial2D *currPoly = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 0, 0);
-                currPoly->coeff[0][0] = 1.0;
-
-                psPolynomial2D *newPoly = NULL;
-                // Must raise trans1->x to the t2x-power.
-                for (psS32 c = 0 ; c < t2x; c++) {
-                    newPoly = multiplyDPoly2D(currPoly, trans1->x);
-                    psFree(currPoly);
-                    currPoly = newPoly;
-                }
-
-                // Must raise trans1->y to the t2y-power.
-                for (psS32 c = 0 ; c < t2y; c++) {
-                    newPoly = multiplyDPoly2D(currPoly, trans1->y);
-                    psFree(currPoly);
-                    currPoly = newPoly;
-                }
-                if (psTraceGetLevel(__func__) >= 6) {
-                    PS_POLY_PRINT_2D(currPoly);
-                }
-
-                // Set the appropriate coeffs in myPT->x
-                for (psS32 i = 0 ; i < (1 + currPoly->nX) ; i++) {
-                    for (psS32 j = 0 ; j < (1 + currPoly->nY) ; j++) {
-                        myPT->y->coeff[i][j]+= currPoly->coeff[i][j] * trans2->y->coeff[t2x][t2y];
-                    }
-                }
-                if (psTraceGetLevel(__func__) >= 6) {
-                    PS_POLY_PRINT_2D(myPT->x);
-                }
-                psFree(currPoly);
+                psFree(newPoly);
             }
         }
@@ -757,4 +903,11 @@
         PS_POLY_PRINT_2D(myPT->y);
     }
+
+    for (psS32 c = 0 ; c < (order + 1) ; c++) {
+        psFree(trans1XPolys[c]);
+        psFree(trans1YPolys[c]);
+    }
+    psFree(trans1XPolys);
+    psFree(trans1YPolys);
 
     psTrace(__func__, 3, "---- %s() end ----\n", __func__);