SOR spline type

Splines in POV-Ray are splitted into several things. When You want to use some object placed "evenly" on surface of sor object you can't do it with spline. But sor is listed in documentation in spline based objects list! So let's add spline base on sor equations. Below description can be also considered as tutorial for introducing new types of splines. Building it I have used equations available in POV-Ray documentation evaluated by my in some old post. This patch introduces new keyword sor_spline. To move data from sor object to sor_spline order of coordinates has to be changed. Old y coordinate is now clock in spline. Advantage is that one splane can hold data from five old sors. That's becouse it can operate up to five dimensions along clock value.

#version 3.5;

#declare S=spline{
   sor_spline
  -1.000000,0.000000*x
   0.000000,0.118143*x
   0.540084,0.620253*x
   0.827004,0.210970*x
   0.962025,0.194093*x
   1.000000,0.286920*x
   1.033755,0.468354*x
};

union{
  sor{
    7
    <0.000000, -1.000000>
    <0.118143,  0.000000>
    <0.620253,  0.540084>
    <0.210970,  0.827004>
    <0.194093,  0.962025>
    <0.286920,  1.000000>
    <0.468354,  1.033755>
  }
  #local C=0;
  #while(C<=1)
    #local R=0;
    #while(R<360)
      sphere{<0,C,-S(C).x> , .01 rotate y*R}
      #local R=R+30;
    #end
    #local C=C+.01;
  #end
  translate 1.5*z-y/2
  pigment{color rgb 1}
}

light_source{-9 1}
background{1}
camera{right x up y location y/9}

How it is implemented

Implementation is a litle complicated but not the worst one :-) First we have to add ne keyword to the system. To call it SOR_SPLINE_TOKEN we have to open parse.h and find enum TOKEN_IDS. Some entries in this list describe tokens for types of splines so let's add SOR_SPLINE_TOKEN there:

enum TOKEN_IDS
{
  ABS_TOKEN = 0,
  :
  :
  LINEAR_SPLINE_TOKEN,
  QUADRATIC_SPLINE_TOKEN,
  CUBIC_SPLINE_TOKEN,
  BEZIER_SPLINE_TOKEN,
  SOR_SPLINE_TOKEN,         // add this line to enable token recognition
  SOR_SPLINE_TOKEN,
  POLYGON_TOKEN,
  PRISM_TOKEN,
  :
  :
};

Now we have to connect SOR_SPLINE_TOKEN with string in source file for proper parsing. This connection is realised via Reserved_Words array in tokenize.cpp. Let's add necessary changes there:

RESERVED_WORD Reserved_Words[LAST_TOKEN] = {
  {AA_THRESHOLD_TOKEN, "aa_threshold"},
  :
  :
  {SOLID_TOKEN, "solid" },
  {SOR_TOKEN, "sor"},
  {SOR_SPLINE_TOKEN, "sor_spline"},     // add this line to connect token with string
  {SPECULAR_TOKEN, "specular"},
  {SPHERE_SWEEP_TOKEN, "sphere_sweep"},
  :
  :
};

Now we have to define this new type of spline for spline code. This list is defined in splines.h. Added here SOR_SPLINE should have value different than any other spline type:

:
:
#define LINEAR_SPLINE         1
#define QUADRATIC_SPLINE      2
#define NATURAL_SPLINE        3
#define CATMULL_ROM_SPLINE    4
#define SOR_SPLINE            5    // add this definition with value different than any other
:
:

Having keyword, token and internal spline type we can force parser to parse it correctly. Parse_Spline function is defined in express.cpp file. Let's open it. In the second expectation block in Parse_Spline there is recognition of type. Similiar part appear in Parse_Spline_Call. Let's add there recognition for our sor spline type:

void Parse_Spline_Call(EXPRESS Express, int *Terms)
{
        :
        :
	if(Token.Token_Id == COMMA_TOKEN)
	{
                :
                :
		switch(Token.Token_Id)
		{
			case SOR_SPLINE_TOKEN:
				spline->Type = SOR_SPLINE;
				break;
			default:
				Error("linear_spline, quadratic_spline, natural_spline, sor_spline or cubic_spline expected.");
				break;
		}
	        :
	        :
        }
        :
        :
}

SPLINE * Parse_Spline()
{
   :
   :
   /* Determine kind of spline */
   EXPECT
     :
     :
     CASE(SOR_SPLINE_TOKEN)           //
       Type = SOR_SPLINE;             //  Add this block to connect scene keyword with spline type
     END_CASE                         // 
     :
     :
   END_EXPECT
   :
   :
}

Since sor calculations are complicated we have to optimize it adding precomputed coefficients. For this purpose we have to change spline structure definition. It is written in frame.h. Let's add array for coefficients and field to determine if coefficients are already precalculated:

struct Spline_Entry
{
  :
  :
  DBL sor_coeff[5][4];      // add this field for four sor coefficents per each dimmension
};

struct Spline_Struct {
  :
  :
  int SOR_Coeffs_Computed;  // add this field to determine if sor coefficients are already precalculated
};

Now rest of changes will be focused on the main spline.cpp file. Now we have to point to the engine when coefficents should be (re)calculated. There are two situations. After spline copy or creation and after every new entry added to list of control points.

SPLINE * Create_Spline(int Type)
{
  :
  :
  New->Coeffs_Computed = false;
  New->SOR_Coeffs_Computed = false;  // add this line to force first precomputation
  :
  :
}

SPLINE * Copy_Spline(SPLINE * Old)
{
  :
  :
  New->Coeffs_Computed = Old->Coeffs_Computed;
  New->SOR_Coeffs_Computed = Old->SOR_Coeffs_Computed;  // add this line to inherit field from base spline
  :
  :
}

void Insert_Spline_Entry(SPLINE * sp, DBL p, EXPRESS v)
{
  :
  :
  sp->Coeffs_Computed = false;
  sp->SOR_Coeffs_Computed = false;      // add this line for sor_spline type
  :
  :
}

Now we have to modify function Get_Spline_Val. There should be block to call (if necessary) precomputation function at begining and block to call function with sor interpolation.

DBL Get_Spline_Val(SPLINE *sp, DBL p, EXPRESS v, int *Terms)
{
    :
    :
    // add below block at begining of function
    if(!sp->SOR_Coeffs_Computed)
    {
        switch(sp->Type)
        {
            case SOR_SPLINE:
                Precompute_SOR_Coeffs(sp);
                break;
            default:
                break;
        }
    }
    :
    :
    switch(sp->Type)
    {
        :
        :
        case SOR_SPLINE:
            for(k=0; k<5; k++)
            {
                /* If outside the spline range, return the first or last point */
                if(i <= 1)
                    v[k] = se[0].vec[k];
                else if(i >= last)
                    v[k] = se[last-1].vec[k];
                else if( last <= 2 )
                    v[k] = se[max(0,last-1)].vec[k];
                /* Else, normal case. */
                else
                    v[k] = SOR_interpolate(se, i-1, k, p);
            }
            break;
        :
        :
    }
    :
    :
}

Finnaly we have to write both functions mentioned above. One for precomputation and one for interpolation.

/*****************************************************************************
*
* FUNCTION
*
*       Precompute_SOR_Coeffs
*
* INPUT
*
*       sp : a pointer to the spline to compute interpolation coefficients for
*
* OUTPUT
*
* RETURNS
*
* AUTHOR
*
*   ABX (abx@abx.art.pl)
*
* DESCRIPTION
*
*       Computes the coefficients used in sor_interpolate.
*
* CHANGES
*
*       2002.08.09 - Initial version by ABX
*
******************************************************************************/

void Precompute_SOR_Coeffs(SPLINE *sp)
{
    int i, k;
    // temporary variables for operations
    // I hope that compilers remove part of it
    DBL b0,b1,b2,b3;
    DBL M00,M01,M02;
    DBL M10,M11,M12;
    DBL M20,M21,M30,M31;
    DBL M2131,M1101,M1202;
    DBL M1101_1,M1101_2;
    DBL A,B,C,D;
    
    for(k = 0; k < 5; k++)
    {
        for(i = 2; i < sp->Number_Of_Entries - 1; i++)
        {
            b0=pow(sp->SplineEntries[i-1].vec[k],2);
            b1=pow(sp->SplineEntries[i].vec[k],2);
            b2=2*sp->SplineEntries[i-1].vec[k]*(sp->SplineEntries[i].vec[k]-sp->SplineEntries[i-2].vec[k])/
               (sp->SplineEntries[i].par-sp->SplineEntries[i-2].par);
            b3=2*sp->SplineEntries[i].vec[k]*(sp->SplineEntries[i+1].vec[k]-sp->SplineEntries[i-1].vec[k])/
               (sp->SplineEntries[i+1].par-sp->SplineEntries[i-1].par);
            M00=pow(sp->SplineEntries[i-1].par,3);
            M01=pow(sp->SplineEntries[i-1].par,2);
            M02=sp->SplineEntries[i-1].par;
            M10=pow(sp->SplineEntries[i].par,3);
            M11=pow(sp->SplineEntries[i].par,2);
            M12=sp->SplineEntries[i].par;
            M20=3*pow(sp->SplineEntries[i-1].par,2);
            M21=2*sp->SplineEntries[i-1].par;
            M30=3*pow(sp->SplineEntries[i].par,2);
            M31=2*sp->SplineEntries[i].par;
            M2131=M21-M31;
            M1101=M11-M01;
            M1202=M12-M02;
            M1101_1=M1101-M31*M1202;
            M1101_2=M21*M1202-M1101;
            A=((b0-b1)*M2131+b2*M1101_1+b3*M1101_2)/
              ((M00-M10)*M2131+M20*M1101_1+M30*M1101_2);
            B=(b2-b3+A*(M30-M20))/M2131;
            C=b3-M30*A-M31*B;
            D=b1-M10*A-M11*B-M12*C;
            
            sp->SplineEntries[i].SOR_coeff[k][0]=A;
            sp->SplineEntries[i].SOR_coeff[k][1]=B;
            sp->SplineEntries[i].SOR_coeff[k][2]=C;
            sp->SplineEntries[i].SOR_coeff[k][3]=D;
        }
    }
    sp->SOR_Coeffs_Computed = true;
}

/*****************************************************************************
*
* FUNCTION
*
*       SOR_interpolate
*
* INPUT
*
*       se : a pointer to the entries in the spline
*       i  : the first point in the interpolation interval
*       k  : which dimension of the 5D vector to interpolate in
*       p  : the parameter to interpolate the value for
*
* OUTPUT
*
* RETURNS
*
*       The value of the kth dimension of the SOR interpolation of the
*        vector at p
*
* AUTHOR
*
*       ABX (abx@abx.art.pl)
*
* DESCRIPTION
*
* CHANGES
*
*       2002.08.09 initial version
*
******************************************************************************/

DBL SOR_interpolate(SPLINE_ENTRY * se, int i, int k, DBL p)
{
    DBL Result=
      sqrt(
        se[i+1].SOR_coeff[k][0]*pow(p,3)
       +se[i+1].SOR_coeff[k][1]*pow(p,2)
       +se[i+1].SOR_coeff[k][2]*p
       +se[i+1].SOR_coeff[k][3]
      );
    return Result;
}

Notes

I think that Coeffs_Computed and SOR_Coeffs_Computed could be connected as one field but at this moment seems it will require recalculations after every type change. Perhaps somebody will imagine better solution.

I wonder if sor_coeff and coeff fields could be moved from Spline_Entry to Spline_Struct and could be made dynamically allocated when necessary. It could save some memory when type of spline not require coefficients.

SOR spline type

Why it is introduced

How it is implemented

Notes

Contact