home *** CD-ROM | disk | FTP | other *** search
/ 3D GFX / 3D GFX.iso / amiutils / i_l / irit5 / triv_lib / triv_der.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-12-31  |  10.5 KB  |  280 lines
  1. /******************************************************************************
  2. * Triv_Der.c - Compute derivatives of tri-variates.                  *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Sep. 94.                          *
  5. ******************************************************************************/
  6.  
  7. #include <ctype.h>
  8. #include <stdio.h>
  9. #include <string.h>
  10. #include "triv_loc.h"
  11.  
  12. /* Define some marcos to make some of the routines below look better. They  */
  13. /* calculate the index of the U, V, W point of the control mesh in Points.  */
  14. #define DERIVED_TV(U, V, W)    TRIV_MESH_UVW(DerivedTV, (U), (V), (W))
  15. #define TV(U, V, W)        TRIV_MESH_UVW(TV, (U), (V), (W))
  16.  
  17.  
  18. /*****************************************************************************
  19. * DESCRIPTION:                                                               M
  20. * Given a trivariate, computes its partial derivative trivariate in          M
  21. * direction Dir.                                 M
  22. *                                                                            *
  23. * PARAMETERS:                                                                M
  24. *   TV:       Trivariate to differentiate.                                   M
  25. *   Dir:      Direction of differentiation. Either U or V or W.              M
  26. *                                                                            *
  27. * RETURN VALUE:                                                              M
  28. *   TrivTVStruct *:   Differentiated trivariate in direction Dir.            M
  29. *                                                                            *
  30. * KEYWORDS:                                                                  M
  31. *   TrivTVDerive, trivariates                                                M
  32. *****************************************************************************/
  33. TrivTVStruct *TrivTVDerive(TrivTVStruct *TV, TrivTVDirType Dir)
  34. {
  35.     switch (TV -> GType) {
  36.     case TRIV_TVBEZIER_TYPE:
  37.         return TrivBzrTVDerive(TV, Dir);
  38.     case TRIV_TVBSPLINE_TYPE:
  39.         return TrivBspTVDerive(TV, Dir);
  40.     default:
  41.         TRIV_FATAL_ERROR(TRIV_ERR_UNDEF_GEOM);
  42.         return NULL;
  43.     }
  44. }
  45.  
  46. /*****************************************************************************
  47. * DESCRIPTION:                                                               M
  48. * Given a Bezier trivariate, computes its partial derivative trivariate in   M
  49. * direction Dir.                                 M
  50. * Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then:   M
  51. * Q(i) = (k - 1) * (P(i+1) - P(i)), i = 0 to k-2.                 V
  52. *                                                                            *
  53. * PARAMETERS:                                                                M
  54. *   TV:       Trivariate to differentiate.                                   M
  55. *   Dir:      Direction of differentiation. Either U or V or W.              M
  56. *                                                                            *
  57. * RETURN VALUE:                                                              M
  58. *   TrivTVStruct *:   Differentiated trivariate in direction Dir. A Bezier   M
  59. *              trivariate.                         M
  60. *                                                                            *
  61. * KEYWORDS:                                                                  M
  62. *   TrivBzrTVDerive, trivariates                                             M
  63. *****************************************************************************/
  64. TrivTVStruct *TrivBzrTVDerive(TrivTVStruct *TV, TrivTVDirType Dir)
  65. {
  66.     CagdBType
  67.     IsNotRational = !TRIV_IS_RATIONAL_TV(TV);
  68.     int i, j, k, l,
  69.     ULength = TV -> ULength,
  70.     VLength = TV -> VLength,
  71.     WLength = TV -> WLength,
  72.     MaxCoord = CAGD_NUM_OF_PT_COORD(TV -> PType);
  73.     TrivTVStruct
  74.         *DerivedTV = NULL;
  75.  
  76.     if (!IsNotRational) {
  77.     TRIV_FATAL_ERROR(TRIV_ERR_RATIONAL_NO_SUPPORT);
  78.     return NULL;
  79.     }
  80.  
  81.     switch (Dir) {
  82.     case TRIV_CONST_U_DIR:
  83.         DerivedTV = TrivBzrTVNew(MAX(ULength - 1, 1), VLength, WLength,
  84.                      TV -> PType);
  85.  
  86.         for (i = 0; i < MAX(ULength - 1, 1); i++)
  87.         for (j = 0; j < VLength; j++)
  88.             for (k = 0; k < WLength; k++) {
  89.             for (l = IsNotRational; l <= MaxCoord; l++)
  90.                 DerivedTV -> Points[l][DERIVED_TV(i, j, k)] =
  91.                     ULength < 2 ? 0.0
  92.                         : (ULength - 1) *
  93.                        (TV -> Points[l][TV(i + 1, j, k)] -
  94.                         TV -> Points[l][TV(i, j, k)]);
  95.             }
  96.         break;
  97.     case TRIV_CONST_V_DIR:
  98.         DerivedTV = TrivBzrTVNew(ULength, MAX(VLength - 1, 1), WLength,
  99.                      TV -> PType);
  100.  
  101.         for (i = 0; i < ULength; i++)
  102.         for (j = 0; j < MAX(VLength - 1, 1); j++)
  103.             for (k = 0; k < WLength; k++) {
  104.             for (l = IsNotRational; l <= MaxCoord; l++)
  105.                 DerivedTV -> Points[l][DERIVED_TV(i, j, k)] =
  106.                 VLength < 2 ? 0.0
  107.                         : (VLength - 1) *
  108.                        (TV -> Points[l][TV(i, j + 1, k)] -
  109.                         TV -> Points[l][TV(i, j, k)]);
  110.             }
  111.         break;
  112.     case TRIV_CONST_W_DIR:
  113.         DerivedTV = TrivBzrTVNew(ULength, VLength, MAX(WLength - 1, 1),
  114.                      TV -> PType);
  115.  
  116.         for (i = 0; i < ULength; i++)
  117.         for (j = 0; j < VLength; j++)
  118.             for (k = 0; k < MAX(WLength - 1, 1); k++) {
  119.             for (l = IsNotRational; l <= MaxCoord; l++)
  120.                 DerivedTV -> Points[l][DERIVED_TV(i, j, k)] =
  121.                 WLength < 2 ? 0.0
  122.                         : (WLength - 1) *
  123.                        (TV -> Points[l][TV(i, j, k + 1)] -
  124.                         TV -> Points[l][TV(i, j, k)]);
  125.             }
  126.         break;
  127.     default:
  128.         TRIV_FATAL_ERROR(TRIV_ERR_DIR_NOT_VALID);
  129.         break;
  130.     }
  131.  
  132.     return DerivedTV;
  133. }
  134. /*****************************************************************************
  135. * DESCRIPTION:                                                               M
  136. * Given a Bspline trivariate, computes its partial derivative trivariate in  M
  137. * direction Dir.                                 M
  138. *   Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then: M
  139. * Q(i) = (k - 1) * (P(i+1) - P(i)) / (Kv(i + k) - Kv(i + 1)), i = 0 to k-2.  V
  140. *                                                                            *
  141. * PARAMETERS:                                                                M
  142. *   TV:       Trivariate to differentiate.                                   M
  143. *   Dir:      Direction of differentiation. Either U or V or W.              M
  144. *                                                                            *
  145. * RETURN VALUE:                                                              M
  146. *   TrivTVStruct *:   Differentiated trivariate in direction Dir. A Bspline  M
  147. *              trivariate.                         M
  148. *                                                                            *
  149. * KEYWORDS:                                                                  M
  150. *   TrivBspTVDerive, trivariates                                             M
  151. *****************************************************************************/
  152. TrivTVStruct *TrivBspTVDerive(TrivTVStruct *TV, TrivTVDirType Dir)
  153. {
  154.     CagdBType
  155.     IsNotRational = !TRIV_IS_RATIONAL_TV(TV);
  156.     int i, j, k, l, NewUOrder, NewVOrder, NewWOrder,
  157.         NewULength, NewVLength, NewWLength,
  158.     ULength = TV -> ULength,
  159.     VLength = TV -> VLength,
  160.     WLength = TV -> WLength,
  161.     UOrder = TV -> UOrder,
  162.     VOrder = TV -> VOrder,
  163.     WOrder = TV -> WOrder,
  164.     MaxCoord = CAGD_NUM_OF_PT_COORD(TV -> PType);
  165.     CagdRType **DPoints,
  166.     *UKv = TV -> UKnotVector,
  167.     *VKv = TV -> VKnotVector,
  168.     *WKv = TV -> WKnotVector,
  169.     **Points = TV -> Points;
  170.     TrivTVStruct
  171.         *DerivedTV = NULL;
  172.  
  173.     if (!IsNotRational) {
  174.     TRIV_FATAL_ERROR(TRIV_ERR_RATIONAL_NO_SUPPORT);
  175.     return NULL;
  176.     }
  177.  
  178.     switch (Dir) {
  179.     case TRIV_CONST_U_DIR:
  180.         NewULength = UOrder < 2 ? ULength : ULength - 1;
  181.         NewUOrder = MAX(UOrder - 1, 1);
  182.         DerivedTV = TrivBspTVNew(NewULength, VLength, WLength,
  183.                      NewUOrder, VOrder, WOrder, TV -> PType);
  184.         CAGD_GEN_COPY(DerivedTV -> UKnotVector, &UKv[UOrder < 2 ? 0 : 1],
  185.               sizeof(CagdRType) * (NewULength + NewUOrder));
  186.         CAGD_GEN_COPY(DerivedTV -> VKnotVector, VKv,
  187.               sizeof(CagdRType) * (VLength + VOrder));
  188.         CAGD_GEN_COPY(DerivedTV -> WKnotVector, WKv,
  189.               sizeof(CagdRType) * (WLength + WOrder));
  190.         DPoints = DerivedTV -> Points;
  191.  
  192.         for (i = 0; i < NewULength; i++)
  193.         for (j = 0; j < VLength; j++)
  194.             for (k = 0; k < WLength; k++) {
  195.             CagdRType
  196.                 Denom = UKv[i + UOrder] - UKv[i + 1];
  197.  
  198.             for (l = IsNotRational; l <= MaxCoord; l++) {
  199.                 if (APX_EQ(Denom, 0.0))
  200.                     Denom = INFINITY;
  201.  
  202.                 DPoints[l][DERIVED_TV(i, j, k)] =
  203.                 UOrder < 2 ? 0.0
  204.                        : (UOrder - 1) *
  205.                           (Points[l][TV(i + 1, j, k)] -
  206.                            Points[l][TV(i, j, k)]) / Denom;
  207.             }
  208.             }
  209.         break;
  210.     case TRIV_CONST_V_DIR:
  211.         NewVLength = VOrder < 2 ? VLength : VLength - 1;
  212.         NewVOrder = MAX(VOrder - 1, 1);
  213.         DerivedTV = TrivBspTVNew(ULength, NewVLength, WLength,
  214.                      UOrder, NewVOrder, WOrder, TV -> PType);
  215.         CAGD_GEN_COPY(DerivedTV -> UKnotVector, UKv,
  216.               sizeof(CagdRType) * (ULength + UOrder));
  217.         CAGD_GEN_COPY(DerivedTV -> VKnotVector, &VKv[VOrder < 2 ? 0 : 1],
  218.               sizeof(CagdRType) * (NewVLength + NewVOrder));
  219.         CAGD_GEN_COPY(DerivedTV -> WKnotVector, WKv,
  220.               sizeof(CagdRType) * (WLength + WOrder));
  221.         DPoints = DerivedTV -> Points;
  222.  
  223.         for (i = 0; i < ULength; i++)
  224.         for (j = 0; j < NewVLength; j++)
  225.             for (k = 0; k < WLength; k++) {
  226.             CagdRType
  227.                 Denom = VKv[j + VOrder] - VKv[j + 1];
  228.  
  229.             for (l = IsNotRational; l <= MaxCoord; l++) {
  230.                 if (APX_EQ(Denom, 0.0))
  231.                 Denom = INFINITY;
  232.  
  233.                 DPoints[l][DERIVED_TV(i, j, k)] =
  234.                 VOrder < 2 ? 0.0
  235.                        : (VOrder - 1) *
  236.                           (Points[l][TV(i, j + 1, k)] -
  237.                            Points[l][TV(i, j, k)]) / Denom;
  238.             }
  239.         }
  240.         break;
  241.     case TRIV_CONST_W_DIR:
  242.         NewWLength = WOrder < 2 ? WLength : WLength - 1;
  243.         NewWOrder = MAX(WOrder - 1, 1);
  244.         DerivedTV = TrivBspTVNew(ULength, VLength, NewWLength,
  245.                      UOrder, VOrder, NewWOrder, TV -> PType);
  246.         CAGD_GEN_COPY(DerivedTV -> UKnotVector, UKv,
  247.               sizeof(CagdRType) * (ULength + UOrder));
  248.         CAGD_GEN_COPY(DerivedTV -> VKnotVector, VKv,
  249.               sizeof(CagdRType) * (VLength + VOrder));
  250.         CAGD_GEN_COPY(DerivedTV -> WKnotVector, &WKv[WOrder < 2 ? 0 : 1],
  251.               sizeof(CagdRType) * (NewWLength + NewWOrder));
  252.         DPoints = DerivedTV -> Points;
  253.  
  254.         for (i = 0; i < ULength; i++)
  255.         for (j = 0; j < VLength; j++)
  256.             for (k = 0; k < NewWLength; k++) {
  257.             CagdRType
  258.                 Denom = WKv[k + WOrder] - WKv[k + 1];
  259.  
  260.             for (l = IsNotRational; l <= MaxCoord; l++) {
  261.                 if (APX_EQ(Denom, 0.0))
  262.                 Denom = INFINITY;
  263.  
  264.                 DPoints[l][DERIVED_TV(i, j, k)] =
  265.                 WOrder < 2 ? 0.0
  266.                        : (WOrder - 1) *
  267.                           (Points[l][TV(i, j, k + 1)] -
  268.                            Points[l][TV(i, j, k)]) / Denom;
  269.             }
  270.         }
  271.         break;
  272.     default:
  273.         TRIV_FATAL_ERROR(TRIV_ERR_DIR_NOT_VALID);
  274.         break;
  275.     }
  276.  
  277.     return DerivedTV;
  278. }
  279.  
  280.