Альфа-параметр центростремительной сплайновой интерполяции Катмулла-Рома

После поиска почти каждой темы по сплайнам Catmull-Rom и, наконец, ее успешной реализации, я застрял в моменте, когда я не знаю, допустил ли я логическую ошибку или мой код просто неверен.

Моя проблема в том, что мой код не примет альфа-параметр для центростремительной интерполяции катмулла. Результаты с разными значениями альфа всегда одинаковы.

Моя главная цель — интерполировать 4 точки (6, когда я удваиваю первую и последнюю точку) с помощью алгоритма центростремительного катетера.
Я реализовал в основном версию C #, которую можно найти здесь (Центростремительный сплайм Catmull-Rom) или также в нескольких других потоках по всему SO (Кривая Catmull-rom без острых выступов и самопересечений, Интерполяция Катмулла-Рома на путях SVG).

Моя реализация:

Сначала у меня есть 3D-изображение с 4 различными вокселями, установленными на один. Я пытаюсь определить эти воксели и сохранить координаты. Например, z = 250, y = 100, x = 323. Я сохраняю эти значения и делаю их двойными, чтобы использовать их в качестве контрольных точек для моего CR-алгоритма.

for(int nx=0; nx<Nx; ++nx)
for(int ny=0; ny<Ny; ++ny)
for(int nz=0; nz<Nz; ++nz)
{
if(TempArray[nz*Ny*Nx+ny*Nx+nx] > 0)
{
Coordinates[counter*3+0] = nz;  // +0 = z0
Coordinates[counter*3+1] = ny;  // +1 = y0
Coordinates[counter*3+2] = nx;  // +2 = x0
++counter;
}
}

Просто чтобы прояснить. Массивы имеют следующую компоновку:
Координаты [z0, y0, x0, z1, y1, x1, …, гп, уп, хп].

Моя реализация интерполяции Catamull-rom выглядит следующим образом:

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<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

Теперь я пытаюсь установить для каждого вокселя (или пикселя в 2D) значение 1, которое рассчитывается методом интерполяции CR. Я в основном хочу рассчитать координаты с помощью алгоритма CR. Если я установлю альфа = 0, результаты будут такими, как ожидалось (я использовал точки, аналогичные примеру в Википедии). Я получаю хорошее самопересечение наверху. Но если я изменю значение на 0,5 или 1, я получу те же результаты.

Теперь я подозреваю, что с типами, которые я использую, что-то не так. Возможно, нецелесообразно приводить целочисленные координаты к двойным или возвращать их к целому числу (с +0,5). Но это не объясняет самопересечение, которое я получаю.
Я на самом деле не предоставил изображение, потому что оно очень похоже на изображения, которые мы имеем в 1 а также 2.
Спасибо всем, кто даже подумывает прочитать это.