// **************************************************
//
//  Computes the limit point and limit normal vector for
//  a given one-ring of a Loop subdivision control mesh.
//
//  Computes the normal height of the bounding prism
//  over a Loop patch's chord triangle.
//
//
//     Leif P. Kobbelt   18/03/98
//
//     kobbelt@informatik.uni-erlangen.de
//
//     University of Erlangen
//     IMMD 9 (Computer Graphics Group)
//     Am Weichselgarten 9
//     91058 Erlangen
//     Germany
//
// **************************************************

#define  MIN_VALENCE   3
#define  MAX_VALENCE  12

// **************************************************

float LimitMask_3[] = {
 .4000, .2000, .2000, .2000,
 .0000, .4082,-.3809,-.0273,
 .0000, .0000,-.3953, .3953 };

float LimitMask_4[] = {
 .4364, .1409, .1409, .1409, .1409,
 .0000, .3162,-.3776,-.3162, .3776,
 .0000, .0000,-.4925, .0000, .4925 };

float LimitMask_5[] = {
 .4714, .1057, .1057, .1057, .1057, .1057,
 .0000, .0000, .2575, .1592,-.1592,-.2575,
 .0000, .2582, .0022,-.2568,-.1610, .1573 };

float LimitMask_6[] = {
 .5000, .0833, .0833, .0833, .0833, .0833, .0833,
 .0000, .0334, .2046, .1712,-.0334,-.2046,-.1712,
 .0000, .2192, .1009,-.1183,-.2192,-.1009, .1183 };

float LimitMask_7[] = {
 .5222, .0683, .0683, .0683, .0683, .0683, .0683, .0683,
 .0000,-.1644,-.0293, .1278, .1887, .1075,-.0547,-.1756,
 .0000,-.0999,-.1878,-.1343, .0203, .1597, .1788, .0632 };

float LimitMask_8[] = {
 .5391, .0576, .0576, .0576, .0576, .0576, .0576, .0576, .0576,
 .0000, .1601, .0691,-.0624,-.1573,-.1601,-.0691, .0624, .1573,
 .0000,-.0217,-.1358,-.1704,-.1052, .0217, .1358, .1704, .1052 };

float LimitMask_9[] = {
 .5522, .0498, .0498, .0498, .0498, .0498, .0498, .0498, .0498, .0498,
 .0000, .1437, .1357, .0642,-.0373,-.1214,-.1486,-.1063,-.0143, .0845,
 .0000, .0370,-.0645,-.1358,-.1436,-.0842, .0146, .1066, .1487, .1212 };

float LimitMask_10[] = {
 .5624, .0438, .0438, .0438, .0438, .0438, .0438, .0438, .0438, .0438, .0438,
 .0000,-.1188,-.0577, .0255, .0989, .1346, .1188, .0577,-.0255,-.0989,-.1346,
 .0000, .0523, .1159, .1352, .1028, .0312,-.0523,-.1159,-.1352,-.1028,-.0312 };

float LimitMask_11[] = {
 .5704, .0391, .0391, .0391, .0391, .0391, .0391, .0391, .0391, .0391, .0391, .0391,
 .0000,-.0109, .0571, .1070, .1229, .0998, .0450,-.0241,-.0855,-.1198,-.1161,-.0755,
 .0000, .1227, .1085, .0598,-.0079,-.0731,-.1150,-.1205,-.0877,-.0270, .0422, .0980 };

float LimitMask_12[] = {
 .5768, .0353, .0353, .0353, .0353, .0353, .0353, .0353, .0353, .0353, .0353, .0353, .0353,
 .0000,-.0354, .0231, .0755, .1076, .1109, .0845, .0354,-.0231,-.0755,-.1076,-.1109,-.0845,
 .0000,-.1065,-.1115,-.0866,-.0385, .0199, .0730, .1065, .1115, .0866, .0385,-.0199,-.0730 };


float *LimitMasks[] = { LimitMask3, LimitMask4, LimitMask5,  LimitMask6,  LimitMask7,
                        LimitMask8, LimitMask9, LimitMask10, LimitMask11, LimitMask12 };

// **************************************************

float RangeMin_3[]  = {
 .0834, .0834, .0834, .0000, .0000, .0000, .0000, .0000, .0313 };
float RangeMax_3[]  = {
 .4000, .5000, .5000, .0834, .0834, .1355, .0834, .0834, .2000 };

float RangeMin_4[]  = {
 .0834, .0834, .0834, .0000, .0000, .0000, .0000, .0000, .0000, .0000 };
float RangeMax_4[]  = {
 .4364, .5000, .5000, .0834, .0834, .1355, .0834, .0834, .1436, .1434 };

float RangeMin_5[]  = {
 .0834, .0834, .0834, .0000, .0000, .0000, .0000, .0000, .0000, .0000, .0000 };
float RangeMax_5[]  = {
 .4715, .5000, .5000, .0834, .0834, .1355, .0834, .0834, .1375, .1058, .1375 };

float RangeMin_6[]  = {
 .0834, .0834, .0834, .0000, .0000, .0000, .0000, .0000, .0000, .0000, .0000,
 .0000 };
float RangeMax_6[]  = {
 .5000, .5000, .5000, .0834, .0834, .1355, .0834, .0834, .1355, .0834, .0834,
 .1355 };

float RangeMin_7[]  = {
 .0834, .0683, .0683, .0000, .0000, .0000, .0000, .0000, .0000, .0000, .0000,
 .0000, .0000 };
float RangeMax_7[]  = {
 .5223, .5000, .5000, .0834, .0834, .1355, .0834, .0834, .1344, .0683, .0683,
 .0683, .1344 };

float RangeMin_8[]  = {
 .0834, .0576, .0576, .0000, .0000, .0000, .0000, .0000, .0000, .0000, .0000,
 .0000, .0000, .0000 };
float RangeMax_8[]  = {
 .5393, .5000, .5000, .0834, .0834, .1355, .0834, .0834, .1338, .0576, .0576,
 .0576, .0576, .1338 };

float RangeMin_9[]  = {
 .0834, .0498, .0498, .0000, .0000, .0000, .0000, .0000, .0000, .0000, .0000,
 .0000, .0000, .0000, .0000 };
float RangeMax_9[]  = {
 .5524, .5000, .5000, .0834, .0834, .1355, .0834, .0834, .1334, .0501, .0498,
 .0498, .0498, .0501, .1334 };

float RangeMin_10[] = {
 .0834, .0438, .0438, .0000, .0000, .0000, .0000, .0000, .0000, .0000, .0000,
 .0000, .0000, .0000, .0000, .0000 };
float RangeMax_10[] = {
 .5627, .5000, .5000, .0834, .0834, .1355, .0834, .0834, .1331, .0449, .0438,
 .0438, .0438, .0438, .0449, .1331 };

float RangeMin_11[] = {
 .0834, .0391, .0391, .0000, .0000, .0000, .0000, .0000, .0000, .0000, .0000,
 .0000, .0000, .0000, .0000, .0000, .0000 };
float RangeMax_11[] = {
 .5707, .5000, .5000, .0834, .0834, .1355, .0834, .0834, .1329, .0411, .0391,
 .0391, .0391, .0391, .0391, .0411, .1329 };

float RangeMin_12[] = {
 .0834, .0353, .0353, .0000, .0000, .0000, .0000, .0000, .0000, .0000, .0000,
.0000, .0000, .0000, .0000, .0000, .0000, .0000 };
float RangeMax_12[] = {
 .5771, .5000, .5000, .0834, .0834, .1355, .0834, .0834, .1327, .0383, .0353,
 .0353, .0353, .0353, .0353, .0353, .0383, .1327 };


float *RangeMin[] = { RangeMin_3, RangeMin_4, RangeMin_5,  RangeMin_6,  RangeMin_7,
                      RangeMin_8, RangeMin_9, RangeMin_10, RangeMin_11, RangeMin_12 };

float *RangeMax[] = { RangeMax_3, RangeMax_4, RangeMax_5,  RangeMax_6,  RangeMax_7,
                      RangeMax_8, RangeMax_9, RangeMax_10, RangeMax_11, RangeMax_12 };

// **************************************************

void GetLimitPoint(int n, Vec *Q,   // <- one-ring sub mesh
                          Vec *P,   // -> limit point
                          Vec *N)   // -> limit normal
{
  assert((n >= MIN_VALENCE) && (n <= MAX_VALENCE));

  Vec Du(0.0,0.0,0.0);
  Vec Dv(0.0,0.0,0.0);

  *P = Q[0] * LimitMasks[n-3][0];

  for (int i=1 ; i<=n ; i++)
    {
      *P += Q[i] * LimitMasks[n-3][i        ];
      Du += Q[i] * LimitMasks[n-3][i+  (n+1)];
      Dv += Q[i] * LimitMasks[n-3][i+2*(n+1)];
    }

  *N = Du % Dv;     // check orientation !!!

  N->normalize();
}

// **************************************************

//  Vertex order in Q is:          ... ...
//
//                               n+5  0   8
//
//                              3   1   2   7
//
//                                4   5   6
//
//  Vertex Q0 has valence n, Vertices Q1 and Q2 are regular.

void GetPrisma(int n, Vec *Q,             // <- one-ring neighborhood of triangle Q0,Q1,Q2
               Vec P0, Vec P1, Vec P2,    // <- limit points for Q0,Q1,Q2
               float &min,&max);          // -> prisma height in normal direction
{
  assert((n >= MIN_VALENCE) && (n <= MAX_VALENCE));

  Vec N = (P1 - P0) % (P2 - P0);     // check orientation !!!

  N.normalize();

  float p = P0 * N;

  min = 0.0;
  max = 0.0;

  for (int i=0 ; i<=n+5 ; i++)
    {
      float d = Q[i] * N - p;

      if (d > 0.0)
        {
          max += d * RangeMax[n-3][i];
          min += d * RangeMin[n-3][i];
        }
      else
        {
          max += d * RangeMin[n-3][i];
          min += d * RangeMax[n-3][i];
        }
    }
}

// **************************************************