home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 2 / Apprentice-Release2.iso / Source Code / C / Applications / RTrace 1.0 / source / patch.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-09-17  |  10.6 KB  |  355 lines  |  [TEXT/KAHL]
  1. /*
  2.  * Copyright (c) 1988, 1992 Antonio Costa, INESC-Norte.
  3.  * All rights reserved.
  4.  *
  5.  * This code received contributions from the following people:
  6.  *
  7.  *  Roman Kuchkuda      - basic ray tracer
  8.  *  Mark VandeWettering - MTV ray tracer
  9.  *  Augusto Sousa       - overall, shading model
  10.  *
  11.  * Redistribution and use in source and binary forms are permitted
  12.  * provided that the above copyright notice and this paragraph are
  13.  * duplicated in all such forms and that any documentation,
  14.  * advertising materials, and other materials related to such
  15.  * distribution and use acknowledge that the software was developed
  16.  * by Antonio Costa, at INESC-Norte. The name of the author and
  17.  * INESC-Norte may not be used to endorse or promote products derived
  18.  * from this software without specific prior written permission.
  19.  * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
  20.  * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
  21.  * WARRANTIES OF MERCHANTIBILITY AND FITNESS FOR A PARTICULAR PURPOSE.
  22.  */
  23. #include "defs.h"
  24. #include "extern.h"
  25.  
  26. /**********************************************************************
  27.  *    RAY TRACING - Patch (part1) - Version 7.3.2                     *
  28.  *                                                                    *
  29.  *    MADE BY    : Antonio Costa, INESC-Norte, October 1988           *
  30.  *    ADAPTED BY : Antonio Costa, INESC-Norte, June 1989              *
  31.  *    MODIFIED BY: Antonio Costa, INESC-Norte, August 1992            *
  32.  **********************************************************************/
  33.  
  34. /***** Patch *****/
  35. #include "patch.h"
  36.  
  37. static real     small_distance;
  38. static xyz_struct k0, k1, k2, k3, k4, k5, k6, k7, p[4];
  39. static patch_ptr patch;
  40.  
  41. void
  42. subpatch_enclose(p0, p1, p2, p3, pmin, pmax)
  43.   xyz_ptr         p0, p1, p2, p3, pmin, pmax;
  44. #define CORRECTION ((real) 0.25)/* 25% correction */
  45. {
  46.   xyz_struct      center;
  47.  
  48.   STRUCT_ASSIGN(*pmax, *p0);
  49.   STRUCT_ASSIGN(*pmin, *p0);
  50.   pmax->x = MAX(pmax->x, p1->x);
  51.   pmax->y = MAX(pmax->y, p1->y);
  52.   pmax->z = MAX(pmax->z, p1->z);
  53.   pmin->x = MIN(pmin->x, p1->x);
  54.   pmin->y = MIN(pmin->y, p1->y);
  55.   pmin->z = MIN(pmin->z, p1->z);
  56.   pmax->x = MAX(pmax->x, p2->x);
  57.   pmax->y = MAX(pmax->y, p2->y);
  58.   pmax->z = MAX(pmax->z, p2->z);
  59.   pmin->x = MIN(pmin->x, p2->x);
  60.   pmin->y = MIN(pmin->y, p2->y);
  61.   pmin->z = MIN(pmin->z, p2->z);
  62.   pmax->x = MAX(pmax->x, p3->x);
  63.   pmax->y = MAX(pmax->y, p3->y);
  64.   pmax->z = MAX(pmax->z, p3->z);
  65.   pmin->x = MIN(pmin->x, p3->x);
  66.   pmin->y = MIN(pmin->y, p3->y);
  67.   pmin->z = MIN(pmin->z, p3->z);
  68.   if (ABS(pmax->x - pmin->x) < threshold_distance)
  69.     pmax->x += threshold_distance;
  70.   if (ABS(pmax->y - pmin->y) < threshold_distance)
  71.     pmax->y += threshold_distance;
  72.   if (ABS(pmax->z - pmin->z) < threshold_distance)
  73.     pmax->z += threshold_distance;
  74.   center.x = (pmax->x + pmin->x) * 0.5;
  75.   center.y = (pmax->y + pmin->y) * 0.5;
  76.   center.z = (pmax->z + pmin->z) * 0.5;
  77.   pmax->x += (pmax->x - center.x) * CORRECTION;
  78.   pmax->y += (pmax->y - center.y) * CORRECTION;
  79.   pmax->z += (pmax->z - center.z) * CORRECTION;
  80.   pmin->x += (pmin->x - center.x) * CORRECTION;
  81.   pmin->y += (pmin->y - center.y) * CORRECTION;
  82.   pmin->z += (pmin->z - center.z) * CORRECTION;
  83. #undef CORRECTION
  84. }
  85. static boolean
  86. solve(a, b, c)
  87.   array2          a;            /* It is modified */
  88.   array1          b, c;
  89. {
  90.   REG int         row_pivot, i, j, k;
  91.   REG real        row_max, col_max, col_ratio, temp0;
  92.   array1          temp;
  93.   int             pivot[3];
  94.  
  95.   for (i = 0; i <= 2; POSINC(i))
  96.   {
  97.     pivot[i] = i;
  98.     row_max = 0.0;
  99.     for (j = 0; j <= 2; POSINC(j))
  100.       row_max = MAX(row_max, ABS(a[i][j]));
  101.     if (ABS(row_max) < ROUNDOFF)
  102.       return FALSE;
  103.     temp[i] = row_max;
  104.   }
  105.   for (i = 0; i <= 1; POSINC(i))
  106.   {
  107.     col_max = ABS(a[i][i]) / temp[i];
  108.     row_pivot = i;
  109.     for (j = SUCC(i); j <= 2; POSINC(j))
  110.     {
  111.       col_ratio = ABS(a[j][i]) / temp[j];
  112.       if (col_ratio > col_max)
  113.       {
  114.         col_max = col_ratio;
  115.         row_pivot = j;
  116.       }
  117.     }
  118.     if (ABS(col_max) < ROUNDOFF)
  119.       return FALSE;
  120.     if (row_pivot > i)
  121.     {
  122.       j = pivot[row_pivot];
  123.       pivot[row_pivot] = pivot[i];
  124.       pivot[i] = j;
  125.       temp0 = temp[row_pivot];
  126.       temp[row_pivot] = temp[i];
  127.       temp[i] = temp0;
  128.       for (j = 0; j <= 2; POSINC(j))
  129.       {
  130.         temp0 = a[row_pivot][j];
  131.         a[row_pivot][j] = a[i][j];
  132.         a[i][j] = temp0;
  133.       }
  134.     }
  135.     for (j = SUCC(i); j <= 2; POSINC(j))
  136.     {
  137.       a[j][i] /= a[i][i];
  138.       temp0 = a[j][i];
  139.       for (k = SUCC(i); k <= 2; POSINC(k))
  140.         a[j][k] -= temp0 * a[i][k];
  141.     }
  142.   }
  143.   if (ABS(a[2][2]) < ROUNDOFF)
  144.     return FALSE;
  145.   c[0] = b[pivot[0]];
  146.   for (i = 1; i <= 2; POSINC(i))
  147.   {
  148.     temp0 = 0.0;
  149.     for (j = 0; j <= PRED(i); POSINC(j))
  150.       temp0 += a[i][j] * c[j];
  151.     c[i] = b[pivot[i]] - temp0;
  152.   }
  153.   c[2] /= a[2][2];
  154.   for (i = 1; i >= 0; POSDEC(i))
  155.   {
  156.     temp0 = 0.0;
  157.     for (j = SUCC(i); j <= 2; POSINC(j))
  158.       temp0 += a[i][j] * c[j];
  159.     c[i] = (c[i] - temp0) / a[i][i];
  160.   }
  161.   return TRUE;
  162. }
  163.  
  164. #define ERROR_LIM ((real) 1.0e+5)
  165. #define F_ERROR(f)\
  166. ((ABS((f)[0]) > ERROR_LIM) OR(ABS((f)[1]) > ERROR_LIM)\
  167. OR(ABS((f)[2]) > ERROR_LIM) ?\
  168. INFINITY : SQR((f)[0]) + SQR((f)[1]) + SQR((f)[2]))
  169.  
  170. boolean
  171. patch_damped_nr(u, v, position, vector)
  172.   real            u, v;
  173.   xyz_ptr         position, vector;
  174. /***** Damped Newton-Raphson *****/
  175. #define DELTA_U_MAX ROUNDOFF
  176. #define DELTA_V_MAX DELTA_U_MAX
  177.  
  178. #define U_MAX ((real) 2.0)
  179. #define U_MIN ((real) -1.0)
  180. #define V_MAX U_MAX
  181. #define V_MIN U_MIN
  182.  
  183. #define ERROR_MAX   ROUNDOFF
  184. #define F_ERROR_MAX (SQR(ERROR_MAX) * 3.0)
  185. #define F_ERROR_LIM (SQR(ERROR_LIM) * 3.0)
  186.  
  187. #define ITERATION_MAX         (24)
  188. #define ITERATION_DAMPING_MAX (12)
  189.  
  190. {
  191.   boolean         intersect;
  192.   REG int         i, j;
  193.   REG real        b0, b1, b2, b3, u0, v0, distance, distance0, f_error;
  194.   REG real        damping, k;
  195.   array1          f, f0, delta;
  196.   array2          df;
  197.   xyz_struct      temp;
  198.  
  199.   /* Initial values */
  200.   b0 = BLEND0(v);
  201.   b1 = BLEND1(v);
  202.   CUBIC(patch, 1, v, k1);
  203.   CUBIC(patch, 3, v, k3);
  204.   MULTIPLY(p[0], b0, k4);
  205.   MULTIPLY(p[1], b0, k5);
  206.   MULTIPLY(p[2], b1, k6);
  207.   MULTIPLY(p[3], b1, k7);
  208.   SUBTRACT(k1, k4, k6, k1);
  209.   SUBTRACT(k3, k5, k7, k3);
  210.   CUBIC(patch, 0, u, k0);
  211.   CUBIC(patch, 2, u, k2);
  212.   b2 = BLEND0(u);
  213.   b3 = BLEND1(u);
  214.   temp.x = k0.x * b0 + k1.x * b2 + k2.x * b1 + k3.x * b3 - position->x;
  215.   temp.y = k0.y * b0 + k1.y * b2 + k2.y * b1 + k3.y * b3 - position->y;
  216.   temp.z = k0.z * b0 + k1.z * b2 + k2.z * b1 + k3.z * b3 - position->z;
  217.   distance = LENGTH(temp);
  218.   f[0] = temp.x - distance * vector->x;
  219.   f[1] = temp.y - distance * vector->y;
  220.   f[2] = temp.z - distance * vector->z;
  221.   i = 0;
  222.   if (F_ERROR(f) > F_ERROR_LIM)
  223.     return FALSE;
  224.   intersect = FALSE;
  225.   while ((i < ITERATION_MAX) AND NOT intersect)
  226.   {
  227.     POSINC(i);
  228.     /* 1st partial derivative (u) */
  229.     b0 = BLEND0(v);
  230.     b1 = BLEND1(v);
  231.     CUBIC(patch, 1, v, k1);
  232.     CUBIC(patch, 3, v, k3);
  233.     MULTIPLY(p[0], b0, k4);
  234.     MULTIPLY(p[1], b0, k5);
  235.     MULTIPLY(p[2], b1, k6);
  236.     MULTIPLY(p[3], b1, k7);
  237.     SUBTRACT(k1, k4, k6, k1);
  238.     SUBTRACT(k3, k5, k7, k3);
  239.     DERIV_CUBIC(patch, 0, u, k0);
  240.     DERIV_CUBIC(patch, 2, u, k2);
  241.     b2 = DERIV_BLEND0(u);
  242.     b3 = DERIV_BLEND1(u);
  243.     df[0][0] = k0.x * b0 + k1.x * b2 + k2.x * b1 + k3.x * b3;
  244.     df[1][0] = k0.y * b0 + k1.y * b2 + k2.y * b1 + k3.y * b3;
  245.     df[2][0] = k0.z * b0 + k1.z * b2 + k2.z * b1 + k3.z * b3;
  246.     /* 2nd partial derivative (v) */
  247.     b0 = BLEND0(u);
  248.     b1 = BLEND1(u);
  249.     CUBIC(patch, 0, u, k0);
  250.     CUBIC(patch, 2, u, k2);
  251.     MULTIPLY(p[0], b0, k4);
  252.     MULTIPLY(p[1], b1, k5);
  253.     MULTIPLY(p[2], b0, k6);
  254.     MULTIPLY(p[3], b1, k7);
  255.     SUBTRACT(k0, k4, k5, k0);
  256.     SUBTRACT(k2, k6, k7, k2);
  257.     DERIV_CUBIC(patch, 1, v, k1);
  258.     DERIV_CUBIC(patch, 3, v, k3);
  259.     b2 = DERIV_BLEND0(v);
  260.     b3 = DERIV_BLEND1(v);
  261.     df[0][1] = k0.x * b2 + k1.x * b0 + k2.x * b3 + k3.x * b1;
  262.     df[1][1] = k0.y * b2 + k1.y * b0 + k2.y * b3 + k3.y * b1;
  263.     df[2][1] = k0.z * b2 + k1.z * b0 + k2.z * b3 + k3.z * b1;
  264.     df[0][2] = -vector->x;
  265.     df[1][2] = -vector->y;
  266.     df[2][2] = -vector->z;
  267.     if (NOT solve(df, f, delta))
  268.       return FALSE;
  269.     /* Damping */
  270.     u0 = u - delta[0];
  271.     v0 = v - delta[1];
  272.     if ((u0 < U_MIN) OR(u0 > U_MAX) OR(v0 < V_MIN) OR(v0 > V_MAX))
  273.       return FALSE;
  274.     distance0 = distance - delta[2];
  275.     b0 = BLEND0(v0);
  276.     b1 = BLEND1(v0);
  277.     CUBIC(patch, 1, v0, k1);
  278.     CUBIC(patch, 3, v0, k3);
  279.     MULTIPLY(p[0], b0, k4);
  280.     MULTIPLY(p[1], b0, k5);
  281.     MULTIPLY(p[2], b1, k6);
  282.     MULTIPLY(p[3], b1, k7);
  283.     SUBTRACT(k1, k4, k6, k1);
  284.     SUBTRACT(k3, k5, k7, k3);
  285.     CUBIC(patch, 0, u0, k0);
  286.     CUBIC(patch, 2, u0, k2);
  287.     b2 = BLEND0(u0);
  288.     b3 = BLEND1(u0);
  289.     f[0] = k0.x * b0 + k1.x * b2 + k2.x * b1 + k3.x * b3 -
  290.       position->x - distance0 * vector->x;
  291.     f[1] = k0.y * b0 + k1.y * b2 + k2.y * b1 + k3.y * b3 -
  292.       position->y - distance0 * vector->y;
  293.     f[2] = k0.z * b0 + k1.z * b2 + k2.z * b1 + k3.z * b3 -
  294.       position->z - distance0 * vector->z;
  295.     damping = 1.0;
  296.     f_error = F_ERROR(f);
  297.     k = 1.0;
  298.     j = 0;
  299.     while (j < ITERATION_DAMPING_MAX)
  300.     {
  301.       POSINC(j);
  302.       k = k * 0.5;
  303.       u0 = u - delta[0] * k;
  304.       v0 = v - delta[1] * k;
  305.       distance0 = distance - delta[2] * k;
  306.       b0 = BLEND0(v0);
  307.       b1 = BLEND1(v0);
  308.       CUBIC(patch, 1, v0, k1);
  309.       CUBIC(patch, 3, v0, k3);
  310.       MULTIPLY(p[0], b0, k4);
  311.       MULTIPLY(p[1], b0, k5);
  312.       MULTIPLY(p[2], b1, k6);
  313.       MULTIPLY(p[3], b1, k7);
  314.       SUBTRACT(k1, k4, k6, k1);
  315.       SUBTRACT(k3, k5, k7, k3);
  316.       CUBIC(patch, 0, u0, k0);
  317.       CUBIC(patch, 2, u0, k2);
  318.       b2 = BLEND0(u0);
  319.       b3 = BLEND1(u0);
  320.       f0[0] = k0.x * b0 + k1.x * b2 + k2.x * b1 + k3.x * b3 -
  321.         position->x - distance0 * vector->x;
  322.       f0[1] = k0.y * b0 + k1.y * b2 + k2.y * b1 + k3.y * b3 -
  323.         position->y - distance0 * vector->y;
  324.       f0[2] = k0.z * b0 + k1.z * b2 + k2.z * b1 + k3.z * b3 -
  325.         position->z - distance0 * vector->z;
  326.       if (F_ERROR(f0) < f_error)
  327.       {
  328.         damping = k;
  329.         ARRAY_ASSIGN(f, f0, 1);
  330.         f_error = F_ERROR(f);
  331.       } else
  332.         break;
  333.     }
  334.     if (f_error > F_ERROR_LIM)
  335.       return FALSE;
  336.     /* New values */
  337.     u -= delta[0] * damping;
  338.     v -= delta[1] * damping;
  339.     distance -= delta[2] * damping;
  340.     intersect = (boolean) ((f_error <= F_ERROR_MAX)
  341.                            AND(ABS(delta[0]) <= DELTA_U_MAX)
  342.                            AND(ABS(delta[1]) <= DELTA_V_MAX));
  343.   }
  344.   if (intersect
  345.       AND(u >= -DELTA_U_MAX) AND(u <= 1.0 + DELTA_U_MAX)
  346.       AND(v >= -DELTA_V_MAX) AND(v <= 1.0 + DELTA_V_MAX)
  347.       AND(distance > threshold_distance) AND(distance < small_distance))
  348.   {
  349.     small_distance = distance;
  350.     patch->u_hit = MAX(0.0, MIN(1.0, u));
  351.     patch->v_hit = MAX(0.0, MIN(1.0, v));
  352.   }
  353.   return intersect;
  354. }
  355.