/* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% %% %% Wroc{\l}aw, September 2022 %% %% %% %% %% %% %% %% LINEAR-TIME ALGORITHM FOR COMPUTING %% %% THE BERNSTEIN-B\'{E}ZIER COEFFICIENTS %% %% OF B-SPLINE BASIS FUNCTIONS %% %% (test program for surfaces) %% %% %% %% GNU C Compiler 11.2.0 %% %% %% %% Filip Chudy, Pawe\{l} Wo\'{z}ny %% %% %% %% %% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% */ #include #include #include #include #define NMULT 50 // N = 50 * n float knotsX[610]; // n + 2m float knotsY[610]; // n + 2m float samplesX[2510]; // N float samplesY[2510]; // N float points[70][70][3]; // (n + m) * (n + m) * d float coeffsX[510][90][20]; // n * (n + 2m) * (m + 1) float coeffsY[510][90][20]; // n * (n + 2m) * (m + 1) float points_work[20][70][3]; // m * (n + m) * d float points_inner_dBC[70][3]; // m * (n + m) * d float h[50]; // m float bsplineX[2510][90]; // N * (n + 2m) float bsplineY[2510][90]; // N * (n + 2m) float post_bezier[90]; // n + 2m float splines[20][90]; // m * (n + 2m) float result_points[2510][2510][3]; // N * N * d float result_points_new[2510][2510][3]; // N * N * d float precision[2510][2510]; float get_digits(float a, float b, float max) { if( (a - b) / b < 1e-8) return max; float result = -log10((a - b) / b); if(result > max) return max; } float float_rand( float min, float max ) { float scale = rand() / (float) RAND_MAX; // [0, 1.0] return min + scale * (max - min); // [min, max] } void generate_knots_and_samples(int n1, int m1, int n2, int m2) { for(int i = 0; i <= m1; i++) knotsX[i] = 1.0; samplesX[0] = knotsX[m1]; int s = 1; for(int i = 1; i <= n1; i++) { knotsX[m1 + i] = knotsX[m1 + i - 1] + float_rand(0.02, 1.0); float step = knotsX[m1 + i] - knotsX[m1 + i - 1]; for(int q = 1; q < NMULT; q++) { samplesX[s] = knotsX[m1 + i - 1] + step * q / NMULT; s++; } samplesX[s] = knotsX[m1 + i]; s++; } for(int i = 1; i <= m1; i++) knotsX[m1 + n1 + i] = knotsX[m1 + n1]; for(int i = 0; i <= m2; i++) knotsY[i] = 1.0; samplesY[0] = knotsY[m2]; s = 1; for(int i = 1; i <= n2; i++) { knotsY[m2 + i] = knotsY[m2 + i - 1] + float_rand(0.02, 1.0); float step = knotsY[m2 + i] - knotsY[m2 + i - 1]; for(int q = 1; q < NMULT; q++) { samplesY[s] = knotsY[m2 + i - 1] + step * q / NMULT; s++; } samplesY[s] = knotsY[m2 + i]; s++; } for(int i = 1; i <= m2; i++) knotsY[m2 + n2 + i] = knotsY[m2 + n2]; } void generate_points(int n1, int m1, int n2, int m2, int d) { for(int i = 0; i < n1 + m1 + 1; i++) { for (int j = 0; j < n2 + m2 + 1; j++) { for(int dim_c = 0; dim_c < d; dim_c++) { points[i][j][dim_c] = float_rand(-1.0, 1.0); } } } } void zero_coeffs(int n1, int m1, int n2, int m2) { for(int j = 0; j < n1; j++) { for(int i = 0; i < n1 + m1; i++) { for(int k = 0; k <= m1; k++) coeffsX[j][i][k] = 0.0; } } for(int j = 0; j < n2; j++) { for(int i = 0; i < n2 + m2; i++) { for(int k = 0; k <= m2; k++) coeffsY[j][i][k] = 0.0; } } } // ************* // NEW ALGORITHM // ************* void compute_new_coeffs(int n1, int m1, int n2, int m2) { float v, c1, c2, c3; for(int j = 0; j < n1; j++) { float numerator = pow(knotsX[j + 1 + m1] - knotsX[j + m1], m1 - 1); coeffsX[j][j + m1][m1] = numerator; for(int k = 2; k <= m1; k++) coeffsX[j][j + m1][m1] /= knotsX[j + k + m1] - knotsX[j + m1]; coeffsX[j][j - m1 + m1][0] = numerator; for(int k = 2; k <= m1; k++) coeffsX[j][j - m1 + m1][0] /= knotsX[j + 1 + m1] - knotsX[j + 1 - k + m1]; } for(int i = n1 - 2; i >= n1 - m1; i--) { c1 = (knotsX[n1 - 1 + m1] - knotsX[i + m1]) / (knotsX[n1 + m1] - knotsX[i + m1]); c2 = (knotsX[n1 + m1] - knotsX[n1 - 1 + m1]) / (knotsX[n1 + m1] - knotsX[i + 1 + m1]); for(int k = m1 - 1; k >= 0; k--) coeffsX[n1 - 1][i + m1][k] = c1 * coeffsX[n1 - 1][i + m1][k + 1] + c2 * coeffsX[n1 - 1][i + 1 + m1][k + 1]; } for(int j = n1 - 2; j >= 0; j--) { for(int i = j - 1; i >= j - m1 + 1; i--) { coeffsX[j][i + m1][m1] = coeffsX[j + 1][i + m1][0]; v = (knotsX[m1 + i + 1 + m1] - knotsX[i + m1]) / (knotsX[m1 + i + 2 + m1] - knotsX[i + 1 + m1]); float denom = (knotsX[j + 1 + m1] - knotsX[i + m1]); c1 = (knotsX[j + m1] - knotsX[i + m1]) / denom; c2 = (knotsX[j + 1 + m1] - knotsX[m1 + i + 2 + m1]) / denom; c3 = (knotsX[m1 + i + 2 + m1] - knotsX[j + m1]) / denom; for(int k = m1 - 1; k >= 0; k--) { coeffsX[j][i + m1][k] = c1 * coeffsX[j][i + m1][k + 1] + v * (c2 * coeffsX[j][i + 1 + m1][k] + c3 * coeffsX[j][i + 1 + m1][k + 1]); } } } for(int j = 0; j < n2; j++) { float numerator = pow(knotsY[j + 1 + m2] - knotsY[j + m2], m2 - 1); coeffsY[j][j + m2][m2] = numerator; for(int k = 2; k <= m2; k++) coeffsY[j][j + m2][m2] /= knotsY[j + k + m2] - knotsY[j + m2]; coeffsY[j][j - m2 + m2][0] = numerator; for(int k = 2; k <= m2; k++) coeffsY[j][j - m2 + m2][0] /= knotsY[j + 1 + m2] - knotsY[j + 1 - k + m2]; } for(int i = n2 - 2; i >= n2 - m2; i--) { c1 = (knotsY[n2 - 1 + m2] - knotsY[i + m2]) / (knotsY[n2 + m2] - knotsY[i + m2]); c2 = (knotsY[n2 + m2] - knotsY[n2 - 1 + m2]) / (knotsY[n2 + m2] - knotsY[i + 1 + m2]); for(int k = m2 - 1; k >= 0; k--) coeffsY[n2 - 1][i + m2][k] = c1 * coeffsY[n2 - 1][i + m2][k + 1] + c2 * coeffsY[n2 - 1][i + 1 + m2][k + 1]; } for(int j = n2 - 2; j >= 0; j--) { for(int i = j - 1; i >= j - m2 + 1; i--) { coeffsY[j][i + m2][m2] = coeffsY[j + 1][i + m2][0]; v = (knotsY[m2 + i + 1 + m2] - knotsY[i + m2]) / (knotsY[m2 + i + 2 + m2] - knotsY[i + 1 + m2]); float denom = (knotsY[j + 1 + m2] - knotsY[i + m2]); c1 = (knotsY[j + m2] - knotsY[i + m2]) / denom; c2 = (knotsY[j + 1 + m2] - knotsY[m2 + i + 2 + m2]) / denom; c3 = (knotsY[m2 + i + 2 + m2] - knotsY[j + m2]) / denom; for(int k = m2 - 1; k >= 0; k--) { coeffsY[j][i + m2][k] = c1 * coeffsY[j][i + m2][k + 1] + v * (c2 * coeffsY[j][i + 1 + m2][k] + c3 * coeffsY[j][i + 1 + m2][k + 1]); } } } } void compute_h(unsigned int n, float t) { float u = 1.0 - t; unsigned int n1 = n + 1; h[0] = 1.0; if(t <= 0.5) { u = t / u; for(int k = 1; k <= n; k++) { h[k] = h[k - 1] * u * (n1 - k); h[k] = h[k] / (k + h[k]); } } else { u = u / t; for(int k = 1; k <= n; k++) { h[k] = h[k - 1] * (n1 - k); h[k] = h[k] / (k * u + h[k]); } } return; } float eval_bezier_1dX(unsigned int n, float t, int j, int i) { float result; result = coeffsX[j][i + n][0]; if(t <= 0.5) { for(int k = 1; k <= n; k++) { float h1 = 1 - h[k]; result = h1 * result + h[k] * coeffsX[j][i + n][k]; } } else { for(int k = 1; k <= n; k++) { float h1 = 1 - h[k]; result = h1 * result + h[k] * coeffsX[j][i + n][k]; } } return result; } float eval_bezier_1dY(unsigned int n, float t, int j, int i) { float result; result = coeffsY[j][i + n][0]; if(t <= 0.5) { for(int k = 1; k <= n; k++) { float h1 = 1 - h[k]; result = h1 * result + h[k] * coeffsY[j][i + n][k]; } } else { for(int k = 1; k <= n; k++) { float h1 = 1 - h[k]; result = h1 * result + h[k] * coeffsY[j][i + n][k]; } } return result; } void eval_bspline_new(int n1, int m1, int n2, int m2, int d, int samples_noX, int samples_noY) { compute_new_coeffs(n1, m1, n2, m2); int j = 0; for(int s = 0; s < samples_noX; s++) { if(samplesX[s] >= knotsX[j + 1 + m1] && j < n1 - 1) j++; float t = (samplesX[s] - knotsX[j + m1]) / (knotsX[j + 1 + m1] - knotsX[j + m1]); compute_h(m1, t); for(int i = j - m1; i <= j; i++) { bsplineX[s][i + m1] = eval_bezier_1dX(m1, t, j, i); } } j = 0; for(int s = 0; s < samples_noY; s++) { if(samplesY[s] >= knotsY[j + 1 + m2] && j < n2 - 1) j++; float t = (samplesY[s] - knotsY[j + m2]) / (knotsY[j + 1 + m2] - knotsY[j + m2]); compute_h(m2, t); for(int i = j - m2; i <= j; i++) { bsplineY[s][i + m2] = eval_bezier_1dY(m2, t, j, i); } } int j1 = 0, j2 = 0; for(int s1 = 0; s1 < samples_noX; s1++) { if(samplesX[s1] >= knotsX[j1 + 1 + m1] && j1 < n1 - 1) j1++; j2 = 0; for(int s2 = 0; s2 < samples_noY; s2++) { if(samplesY[s2] >= knotsY[j2 + 1 + m1] && j2 < n2 - 1) j2++; for(int dim_c = 0; dim_c < d; dim_c++) { result_points_new[s1][s2][dim_c] = 0.; } for(int i1 = j1 - m1; i1 <= j1; i1++) { for(int i2 = j2 - m2; i2 <= j2; i2++) { for(int dim_c = 0; dim_c < d; dim_c++) { result_points_new[s1][s2][dim_c] += bsplineX[s1][i1 + m1] * bsplineY[s2][i2 + m2] * points[i1 + m1][i2 + m2][dim_c]; } } } } } } // ***************** // BLIND DE BOOR-COX // ***************** void deBoorCoxInner(int n, int m, int i1, int j, float u, int d) { for(int i = j - m; i <= j; i++) { for(int dim_c = 0; dim_c < d; dim_c++) { points_work[0][i + m][dim_c] = points[i1][i + m][dim_c]; } } for(int k = 1; k <= m; k++) { for(int i = j - m + k; i <= j; i++) { float a = (u - knotsY[i + m]) / (knotsY[i + 1 + m - k + m] - knotsY[i + m]); for(int dim_c = 0; dim_c < d; dim_c++) { points_work[k][i + m][dim_c] = a * points_work[k - 1][i + m][dim_c] + (1.0 - a) * points_work[k - 1][i - 1 + m][dim_c]; } } } // result stored in points_work[m][j + m][0..d-1] } void deBoorCoxOuter(int n, int m, int j, float u, int d) { for(int i = j - m; i <= j; i++) { for(int dim_c = 0; dim_c < d; dim_c++) { points_work[0][i + m][dim_c] = points_inner_dBC[i + m][dim_c]; } } for(int k = 1; k <= m; k++) { for(int i = j - m + k; i <= j; i++) { float a = (u - knotsX[i + m]) / (knotsX[i + 1 + m - k + m] - knotsX[i + m]); for(int dim_c = 0; dim_c < d; dim_c++) { points_work[k][i + m][dim_c] = a * points_work[k - 1][i + m][dim_c] + (1.0 - a) * points_work[k - 1][i - 1 + m][dim_c]; } } } // result stored in points_work[m][j + m][0..d-1] } void many_deBoorCox(int n1, int m1, int n2, int m2, int d, int samples_noX, int samples_noY) { int j1 = 0, j2 = 0; for(int s1 = 0; s1 < samples_noX; s1++) { if(samplesX[s1] >= knotsX[j1 + 1 + m1] && j1 < n1 - 1) j1++; j2 = 0; for(int s2 = 0; s2 < samples_noY; s2++) { if(samplesY[s2] >= knotsY[j2 + 1 + m2] && j2 < n2 - 1) j2++; for(int i1 = j1 - m1; i1 <= j1; i1++) { deBoorCoxInner(n2, m2, i1 + m1, j2, samplesY[s2], d); for(int dim_c = 0; dim_c < d; dim_c++) { points_inner_dBC[i1 + m1][dim_c] = points_work[m2][j2 + m2][dim_c]; } } deBoorCoxOuter(n1, m1, j1, samplesX[s1], d); for(int dim_c = 0; dim_c < d; dim_c++) { result_points[s1][s2][dim_c] = points_work[m1][j1 + m1][dim_c]; } } } // result stored in result_points[no. of sets][curve sampling][point dim] } // ************************* void toy_case() { // TOY CASE TO CHECK CORRECTNESS int n1 = 1, n2 = 1, m1 = 2, m2 = 2; knotsX[0] = 0.0; knotsX[1] = 0.0; knotsX[2] = 0.0; knotsX[3] = 1.0; knotsX[4] = 1.0; knotsX[5] = 1.0; samplesX[0] = 0.0; samplesX[1] = 0.5; samplesX[2] = 0.625; samplesX[3] = 0.874; samplesX[4] = 1.0; knotsY[0] = 0.0; knotsY[1] = 0.0; knotsY[2] = 0.0; knotsY[3] = 1.0; knotsY[4] = 1.0; knotsY[5] = 1.0; samplesY[0] = 0.0; samplesY[1] = 0.5; samplesY[2] = 0.625; samplesY[3] = 0.874; samplesY[4] = 1.0; int samples_noX = 5, samples_noY = 5; int d = 3; /* points[0][0][0] = 0.0; points[0][0][1] = 0.0; points[0][0][2] = 0.0; points[0][1][0] = 0.0; points[0][1][1] = 1.0; points[0][1][2] = -0.3; points[1][0][0] = 1.0; points[1][0][1] = 0.0; points[1][0][2] = 1.0; points[1][1][0] = 1.0; points[1][1][1] = 1.0; points[1][1][2] = 0.0; */ points[0][0][0] = 0.0; points[0][0][1] = 0.0; points[0][0][2] = 0.0; points[0][1][0] = 0.0; points[0][1][1] = 1.0; points[0][1][2] = 2.0; points[0][2][0] = 0.0; points[0][2][1] = 2.0; points[0][2][2] = 0.0; points[1][0][0] = 1.0; points[1][0][0] = 0.0; points[1][0][1] = -1.0; points[1][1][0] = 1.4; points[1][1][1] = 0.9; points[1][1][2] = 0.0; points[1][2][0] = 1.0; points[1][2][1] = 2.0; points[1][2][2] = 0.0; points[2][0][0] = 2.0; points[2][0][0] = 0.0; points[2][0][1] = 0.0; points[2][1][0] = 2.0; points[2][1][1] = 1.0; points[2][1][2] = 0.0; points[2][2][0] = 2.0; points[2][2][1] = 2.0; points[2][2][2] = 0.0; eval_bspline_new(n1, m1, n2, m2, d, samples_noX, samples_noY); many_deBoorCox(n1, m1, n2, m2, d, samples_noX, samples_noY); /* printf("BEZIER COEFFS\n"); for(int j = 0; j < n; j++) { for(int i = j - m; i <= j; i++) { printf("j = %d, i = %d:\n", j, i); printf("%1.3f %1.3f %1.3f\n", coeffs[j][i + m][0], coeffs[j][i + m][1], coeffs[j][i + m][2]); } } */ printf("LEFT = dBC; RIGHT = new\n"); for(int i = 0; i < samples_noX; i++) { for(int j = 0; j < samples_noY; j++) { printf("s(%f, %f) = [%1.3f, %1.3f, %1.3f] = [%1.3f, %1.3f, %1.3f]\n", samplesX[i], samplesY[j], result_points[i][j][0], result_points[i][j][1], result_points[i][j][2], result_points_new[i][j][0], result_points_new[i][j][1], result_points_new[i][j][2]); /* printf("%1.3f %1.3f %1.3f\n", (1. - samplesX[i]) * (1. - samplesY[j]) * points[0][0][0] + (1. - samplesX[i]) * samplesY[j] * points[0][1][0] + samplesX[i] * (1. - samplesY[j]) * points[1][0][0] + samplesX[i] * samplesY[j] * points[1][1][0], (1. - samplesX[i]) * (1. - samplesY[j]) * points[0][0][1] + (1. - samplesX[i]) * samplesY[j] * points[0][1][1] + samplesX[i] * (1. - samplesY[j]) * points[1][0][1] + samplesX[i] * samplesY[j] * points[1][1][1], (1. - samplesX[i]) * (1. - samplesY[j]) * points[0][0][2] + (1. - samplesX[i]) * samplesY[j] * points[0][1][2] + samplesX[i] * (1. - samplesY[j]) * points[1][0][2] + samplesX[i] * samplesY[j] * points[1][1][2]); */ } } } void do_tests(int max_m1, int max_n1, int step_m, int step_n) { FILE *out_table; FILE *out_prec; char filename[35]; out_table = fopen("result_surf.csv", "w"); int samples_noX, samples_noY; time_t randtime; srand((unsigned) time(&randtime)); clock_t start_t, end_t; clock_t total_start_t, total_end_t; double elapsed_total = 0.0; // total_start_t = clock(); int d = 3; fprintf(out_table, "m1,m2,n1,n2,samples_noX, samples_noY,deBoor,new\n"); for(int m1 = 3; m1 <= max_m1; m1 += step_m) { for(int m2 = m1 - 2; m2 <= m1; m2 += step_m) { for(int n1 = 10; n1 <= max_n1; n1 += step_n) { for(int n2 = n1; n2 <= n1; n2 += step_n) { double elapsed_new = 0.0, elapsed_dBC = 0.0; samples_noX = n1 * NMULT + 1; samples_noY = n2 * NMULT + 1; for(int x = 0; x < samples_noX; x++) { for(int y = 0; y < samples_noY; y++) { precision[x][y] = 0.0; } } for(int q = 0; q < 10; q++) { generate_knots_and_samples(n1, m1, n2, m2); generate_points(n1, m1, n2, m2, d); zero_coeffs(n1, m1, n2, m2); start_t = clock(); eval_bspline_new(n1, m1, n2, m2, d, samples_noX, samples_noY); end_t = clock(); elapsed_new += (double)(end_t - start_t) / CLOCKS_PER_SEC; start_t = clock(); many_deBoorCox(n1, m1, n2, m2, d, samples_noX, samples_noY); end_t = clock(); elapsed_dBC += (double)(end_t - start_t) / CLOCKS_PER_SEC; for(int x = 0; x < samples_noX; x++) { for(int y = 0; y < samples_noY; y++) { float prec = 0.0; for(int dim_c = 0; dim_c < d; dim_c++) { prec += get_digits(result_points_new[x][y][dim_c], result_points[x][y][dim_c], 8.0); } prec /= d; precision[x][y] = (prec + q * precision[x][y]) / (q + 1); } } } sprintf(filename, "prec_surf_%d_%d_%d_%d.csv", n1, n2, m1, m2); out_prec = fopen(filename, "w"); for(int x = 0; x < samples_noX; x++) { for(int y = 0; y < samples_noY; y++) { fprintf(out_prec, "%1.3f,", precision[x][y]); } fprintf(out_prec, "\n"); } fclose(out_prec); fprintf(out_table, "%d,%d,%d,%d,%d,%d,%lf,%lf\n", m1, m2, n1, n2, samples_noX, samples_noY, elapsed_dBC, elapsed_new); } } } } fclose(out_table); // total_end_t = clock(); // elapsed_total = (double)(total_end_t - total_start_t) / CLOCKS_PER_SEC; // printf("Total time: %lf\n", elapsed_total); } int main() { // toy_case(); // make_latex_table(); do_tests(9, 50, 2, 20); return 0; }