title: Centripetal Catmull-Rom spline interpolation alpha parameter
body: <p>After searching almost every topic on Catmull-Rom splines and finally implementing it successfully I'm now stuck at a point where I don't know if I made a logical mistake or if my code is simply wrong.</p>
  
  <p>My problem is, that my code won't accept the alpha parameter for the centripetal catmull-rom spline interpolation. The results with different alpha values are always the same. </p>
  
  <p>My main goal is to interpolate 4 points (6 when I double the first and last points) with the centripetal catmull-rom algorithm. 
  I implemented basically the C# version which can be found here (<a href="https://en.wikipedia.org/wiki/Centripetal_Catmull%E2%80%93Rom_spline" rel="nofollow noreferrer" title="wiki">Centripetal Catmull-Rom Spline</a>) or also in a couple of other threads throughout SO (<a href="https://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections" title="Thread 1">Catmull-rom curve with no cusps and no self-intersections</a>, <a href="https://stackoverflow.com/questions/30748316/catmull-rom-interpolation-on-svg-paths" title="Thread 2">Catmull-Rom interpolation on SVG Paths</a>). </p>
  
  <p>My implementation:</p>
  
  <p>First I have an 3D image with 4 different voxels set to one. I try to identify those voxels and save the coordinates. For example z=250, y=100, x=323. I save those values and cast them as double to use them as control points for my CR-algorithm.   </p>
  
  <pre><code>for(int nx=0; nx&lt;Nx; ++nx)
  for(int ny=0; ny&lt;Ny; ++ny)
  for(int nz=0; nz&lt;Nz; ++nz)
  {
  if(TempArray[nz*Ny*Nx+ny*Nx+nx] &gt; 0)
  {
  Coordinates[counter*3+0] = nz;  // +0 = z0
  Coordinates[counter*3+1] = ny;  // +1 = y0
  Coordinates[counter*3+2] = nx;  // +2 = x0
  ++counter;
  }
  }
  </code></pre>
  
  <p>Just to clearify. The arrays have the following layout:
  Coordinates[z0,y0,x0,z1,y1,x1,...,zn,yn,xn].</p>
  
  <p>My implementation of the catamull-rom interpolation is as follows:</p>
  
  <pre><code>double P0z, P0y, P0x, P1z, P1y, P1x, P2z, P2y, P2x, P3z, P3y, P3x;
  P0z = Coordinates[n*3+0];   
  P1z = Coordinates[n*3+3];   
  P2z = Coordinates[n*3+6];   
  P3z = Coordinates[n*3+9];
  P0y = Coordinates[n*3+1];   
  P1y = Coordinates[n*3+4];   
  P2y = Coordinates[n*3+7];   
  P3y = Coordinates[n*3+10];
  P0x = Coordinates[n*3+2];   
  P1x = Coordinates[n*3+5];   
  P2x = Coordinates[n*3+8];   
  P3x = Coordinates[n*3+11];
  
  double disAxBx, disAyBy, disAzBz; disAzBz=disAyBy=disAxBx=0.0;
  
  disAzBz = pow( P1z - P0z, 2 );
  disAyBy = pow( P1y - P0y, 2 );
  disAxBx = pow( P1x - P0x, 2 );
  
  double distanceVec = pow( disAzBz + disAyBy + disAxBx, 0.5 );
  
  t0 = 0.0;
  t1 = pow( distanceVec, alpha) + t0;  // This is the part where the alpha comes in
  t2 = pow( distanceVec, alpha) + t1;
  t3 = pow( distanceVec, alpha) + t2;
  
  for ( double t=t1; t&lt;t2; t+= (t2-t1)/(NumberOfPoints) )  // t starts at t1 and runs till t2 (like in the examples)
  {
  z=y=x=0;
  
  A1[0] = ( ( (t1-t)/(t1-t0) ) * P0z ) + ( ( (t-t0)/(t1-t0) ) * P1z );
  A1[1] = ( ( (t1-t)/(t1-t0) ) * P0y ) + ( ( (t-t0)/(t1-t0) ) * P1y );
  A1[2] = ( ( (t1-t)/(t1-t0) ) * P0x ) + ( ( (t-t0)/(t1-t0) ) * P1x );
  
  A2[0] = ( ( (t2-t)/(t2-t1) ) * P1z ) + ( ( (t-t1)/(t2-t1) ) * P2z );
  A2[1] = ( ( (t2-t)/(t2-t1) ) * P1y ) + ( ( (t-t1)/(t2-t1) ) * P2y );
  A2[2] = ( ( (t2-t)/(t2-t1) ) * P1x ) + ( ( (t-t1)/(t2-t1) ) * P2x );
  
  A3[0] = ( ( (t3-t)/(t3-t2) ) * P2z ) + ( ( (t-t2)/(t3-t2) ) * P3z );
  A3[1] = ( ( (t3-t)/(t3-t2) ) * P2y ) + ( ( (t-t2)/(t3-t2) ) * P3y );
  A3[2] = ( ( (t3-t)/(t3-t2) ) * P2x ) + ( ( (t-t2)/(t3-t2) ) * P3x );
  
  B1[0] = ( ( (t2-t)/(t2-t0) ) * A1[0] ) + ( ( (t-t0)/(t2-t0) ) * A2[0] );
  B1[1] = ( ( (t2-t)/(t2-t0) ) * A1[1] ) + ( ( (t-t0)/(t2-t0) ) * A2[1] );
  B1[2] = ( ( (t2-t)/(t2-t0) ) * A1[2] ) + ( ( (t-t0)/(t2-t0) ) * A2[2] );
  
  B2[0] = ( ( (t3-t)/(t3-t1) ) * A2[0] ) + ( ( (t-t1)/(t3-t1) ) * A3[0] );
  B2[1] = ( ( (t3-t)/(t3-t1) ) * A2[1] ) + ( ( (t-t1)/(t3-t1) ) * A3[1] );
  B2[2] = ( ( (t3-t)/(t3-t1) ) * A2[2] ) + ( ( (t-t1)/(t3-t1) ) * A3[2] );
  
  C[0] =int ( ( ( (t2-t)/(t2-t1) ) * B1[0] ) + ( ( (t-t1)/(t2-t1) ) * B2[0] ) +0.5 );
  C[1] =int ( ( ( (t2-t)/(t2-t1) ) * B1[1] ) + ( ( (t-t1)/(t2-t1) ) * B2[1] ) +0.5 );
  C[2] =int ( ( ( (t2-t)/(t2-t1) ) * B1[2] ) + ( ( (t-t1)/(t2-t1) ) * B2[2] ) +0.5 );
  
  z=C[0]; y=C[1]; x=C[2];
  
  OutputArray[z*Nx*Ny+y*Nx+x] = 1.f;          
  }
  } // Loop over 4 consecutive points
  </code></pre>
  
  <p>Now what I try to do is, to set every voxel (or pixel in 2D) to 1, which is calculated by the CR-interpolation method. I basically want to calculate the coordinates with the CR algorithm. If I set alpha=0, the results are as expected (I used similar points to the wikipedia example). I get a nice self-intersection at the top. But if I change the value to 0.5 or 1, I still get the same results. </p>
  
  <p>Now I suspect that there is something wrong with the types I use. Probably it is not wise to cast the integer coordinates as double or to cast them back to integer (with +0.5). But that would not explain the self-intersection I get. 
  I actually did not provide an image, because it is very similar to the images we have in <a href="https://en.wikipedia.org/wiki/Centripetal_Catmull%E2%80%93Rom_spline" rel="nofollow noreferrer" title="wiki">1</a> and <a href="https://stackoverflow.com/questions/9489736/catmull-rom-curve-with-no-cusps-and-no-self-intersections" title="Thread 1">2</a>. 
  <strong>Thanks to everyone who even considers to read this.</strong></p>