|
LR-splines
0.5
|
00001 00002 #include <stdio.h> 00003 #include <stdlib.h> 00004 #include <iostream> 00005 #include <string.h> 00006 #include <sstream> 00007 #include <GoTools/geometry/SplineSurface.h> 00008 #include <GoTools/trivariate/SplineVolume.h> 00009 #include "LRSpline/LRSplineSurface.h" 00010 #include "LRSpline/LRSplineVolume.h" 00011 00012 using namespace Go; 00013 using namespace LR; 00014 using namespace std; 00015 00016 int main(int argc, char **argv) { 00017 00018 // set default parameter values 00019 const double TOL = 1e-6; 00020 int p1 = 3; 00021 int p2 = 4; 00022 int p3 = 5; 00023 int n1 = 8; 00024 int n2 = 5; 00025 int n3 = 7; 00026 int dim = 4; 00027 bool rat = false; 00028 bool vol = false; 00029 stringstream parameters; 00030 parameters << " parameters: \n" \ 00031 " -p1 <n> polynomial ORDER (degree+1) in first parametric direction\n" \ 00032 " -p2 <n> polynomial order in second parametric direction\n" \ 00033 " -p3 <n> polynomial order in third parametric direction (volume only)\n" \ 00034 " -n1 <n> number of basis functions in first parametric direction\n" \ 00035 " -n2 <n> number of basis functions in second parametric direction\n" \ 00036 " -n2 <n> number of basis functions in third parametric direction (volume only)\n" \ 00037 " -dim <n> dimension of the controlpoints\n" \ 00038 " -vol test trivariate volumes instead\n" \ 00039 " -help display (this) help information\n" ; 00040 parameters << " default values\n"; 00041 parameters << " -p = { " << p1 << ", " << p2 << ", " << p3 << " }\n"; 00042 parameters << " -n = { " << n1 << ", " << n2 << ", " << n3 << " }\n"; 00043 parameters << " -vol = { " << ((vol)?"true":"false") << " }\n"; 00044 00045 // read input 00046 for(int i=1; i<argc; i++) { 00047 if(strcmp(argv[i], "-p1") == 0) 00048 p1 = atoi(argv[++i]); 00049 else if(strcmp(argv[i], "-p2") == 0) 00050 p2 = atoi(argv[++i]); 00051 else if(strcmp(argv[i], "-p3") == 0) { 00052 p3 = atoi(argv[++i]); 00053 vol = true; 00054 } else if(strcmp(argv[i], "-n1") == 0) 00055 n1 = atoi(argv[++i]); 00056 else if(strcmp(argv[i], "-n2") == 0) 00057 n2 = atoi(argv[++i]); 00058 else if(strcmp(argv[i], "-n3") == 0) { 00059 n3 = atoi(argv[++i]); 00060 vol = true; 00061 } else if(strcmp(argv[i], "-dim") == 0) 00062 dim = atoi(argv[++i]); 00063 else if(strcmp(argv[i], "-vol") == 0) 00064 vol = true; 00065 else if(strcmp(argv[i], "-help") == 0) { 00066 cout << "usage: " << argv[0] << endl << parameters.str(); 00067 exit(0); 00068 } else { 00069 cout << "usage: " << argv[0] << endl << parameters.str(); 00070 exit(1); 00071 } 00072 } 00073 00074 // do some error testing on input 00075 if(n1 < p1) { 00076 cerr << "ERROR: n1 must be greater or equal to p1\n"; 00077 exit(2); 00078 } else if(n2 < p2) { 00079 cerr << "ERROR: n2 must be greater or equal to p2\n"; 00080 exit(2); 00081 } 00082 00083 00084 // make a uniform integer knot vector 00085 double knot_u[n1+p1]; 00086 double knot_v[n2+p2]; 00087 double knot_w[n3+p3]; 00088 for(int i=0; i<p1+n1; i++) 00089 knot_u[i] = (i<p1) ? 0 : (i>n1) ? n1-p1+1 : i-p1+1; 00090 for(int i=0; i<p2+n2; i++) 00091 knot_v[i] = (i<p2) ? 0 : (i>n2) ? n2-p2+1 : i-p2+1; 00092 for(int i=0; i<p3+n3; i++) 00093 knot_w[i] = (i<p3) ? 0 : (i>n3) ? n3-p3+1 : i-p3+1; 00094 00095 // create a list of "random" control points (all between 0.1 and 1.1) 00096 int nCP = (vol) ? n1*n2*n3 : n1*n2; 00097 nCP *= (dim+rat); 00098 double cp[nCP]; 00099 int k=0; 00100 for(int i=0; i<nCP; i++) // 839 as a generator over Z_853 gives a period of 425. Should suffice 00101 cp[k++] = (i*839 % 853) / 853.0 + 0.1; // rational weights also random and thus we need >0 00102 00103 00104 bool oneFail = false; 00105 if(!vol) { 00106 // make two identical surfaces 00107 SplineSurface ss(n1, n2, p1, p2, knot_u, knot_v, cp, dim, rat); 00108 LRSplineSurface lr(n1, n2, p1, p2, knot_u, knot_v, cp, dim, rat); 00109 00110 // compare function values on edges, knots and in between the knots 00111 // as well as all derivatives (up to first derivatives) 00112 vector<Point> lr_pts(6), ss_pts(6); 00113 vector<bool> assertion_passed; 00114 vector<double> par_u_values; 00115 vector<double> par_v_values; 00116 for(double u=0; u<=n1-p1+1; u+=0.5) { 00117 for(double v=0; v<=n2-p2+1; v+=0.5) { 00118 lr.point(lr_pts, u,v, 2); 00119 ss.point(ss_pts, u,v, 2); 00120 00121 bool correct = true; 00122 for(int i=0; i<6; i++) 00123 for(int d=0; d<dim; d++) 00124 if( fabs(lr_pts[i][d]-ss_pts[i][d]) > TOL ) // absolute error 00125 correct = false; 00126 cout << " LR(" << u << ", " << v << ") = " << lr_pts[0] << endl; 00127 cout << " SS(" << u << ", " << v << ") = " << ss_pts[0] << endl; 00128 cout << "dx LR(" << u << ", " << v << ") = " << lr_pts[1] << endl; 00129 cout << "dx SS(" << u << ", " << v << ") = " << ss_pts[1] << endl; 00130 cout << "dy LR(" << u << ", " << v << ") = " << lr_pts[2] << endl; 00131 cout << "dy SS(" << u << ", " << v << ") = " << ss_pts[2] << endl; 00132 cout << "d2x LR(" << u << ", " << v << ") = " << lr_pts[3] << endl; 00133 cout << "d2x SS(" << u << ", " << v << ") = " << ss_pts[3] << endl; 00134 cout << "dxdy LR(" << u << ", " << v << ") = " << lr_pts[4] << endl; 00135 cout << "dxdy SS(" << u << ", " << v << ") = " << ss_pts[4] << endl; 00136 cout << "d2y LR(" << u << ", " << v << ") = " << lr_pts[5] << endl; 00137 cout << "d2y SS(" << u << ", " << v << ") = " << ss_pts[5] << endl; 00138 00139 // collect results for summary at the end 00140 assertion_passed.push_back(correct); 00141 par_u_values.push_back(u); 00142 par_v_values.push_back(v); 00143 00144 if(correct) 00145 cout << "Parameter (" << u << ", " << v << ") evaluated OK\n"; 00146 else { 00147 cout << "ASSERTION FAILED at (" << u << ", " << v << ")\n"; 00148 } 00149 cout << endl; 00150 } 00151 } 00152 00153 // display results 00154 cout << " ======= RESULT (SURFACE) SUMMARY ======== \n\n"; 00155 cout << "_____u____________v_____________ASSERT__\n"; 00156 for(uint i=0; i<assertion_passed.size(); i++) { 00157 printf("%10.4g %10.4g ", par_u_values[i], par_v_values[i]); 00158 if(assertion_passed[i]) 00159 cout << "OK\n"; 00160 else 00161 cout << "FAIL\n"; 00162 if(!assertion_passed[i]) 00163 oneFail = true; 00164 } 00165 00166 } else { 00167 // make two identical volumes 00168 SplineVolume sv(n1, n2, n3, p1, p2, p3, knot_u, knot_v, knot_w, cp, dim, rat); 00169 LRSplineVolume lr(n1, n2, n3, p1, p2, p3, knot_u, knot_v, knot_w, cp, dim, rat); 00170 00171 // compare function values on edges, knots and in between the knots 00172 // as well as all derivatives (up to first derivatives) 00173 vector<Point> lr_pts(10), sv_pts(10); 00174 vector<bool> assertion_passed; 00175 vector<double> par_u_values; 00176 vector<double> par_v_values; 00177 vector<double> par_w_values; 00178 for(double u=0; u<=n1-p1+1; u+=0.5) { 00179 for(double v=0; v<=n2-p2+1; v+=0.5) { 00180 for(double w=0; w<=n3-p3+1; w+=0.5) { 00181 lr.point(lr_pts, u,v,w,2); 00182 sv.point(sv_pts, u,v,w,2); 00183 00184 bool correct = true; 00185 for(int i=0; i<10; i++) 00186 for(int d=0; d<dim; d++) 00187 if( fabs(lr_pts[i][d]-sv_pts[i][d]) > TOL ) // absolute error 00188 correct = false; 00189 cout << " LR(" << u << ", " << v << ") = " << lr_pts[0] << endl; 00190 cout << " SS(" << u << ", " << v << ") = " << sv_pts[0] << endl; 00191 cout << "dx LR(" << u << ", " << v << ") = " << lr_pts[1] << endl; 00192 cout << "dx SS(" << u << ", " << v << ") = " << sv_pts[1] << endl; 00193 cout << "dy LR(" << u << ", " << v << ") = " << lr_pts[2] << endl; 00194 cout << "dy SS(" << u << ", " << v << ") = " << sv_pts[2] << endl; 00195 cout << "dy LR(" << u << ", " << v << ") = " << lr_pts[3] << endl; 00196 cout << "dy SS(" << u << ", " << v << ") = " << sv_pts[3] << endl; 00197 cout << "d2x LR(" << u << ", " << v << ") = " << lr_pts[4] << endl; 00198 cout << "d2x SS(" << u << ", " << v << ") = " << sv_pts[4] << endl; 00199 cout << "dxdy LR(" << u << ", " << v << ") = " << lr_pts[5] << endl; 00200 cout << "dxdy SS(" << u << ", " << v << ") = " << sv_pts[5] << endl; 00201 cout << "dxdz LR(" << u << ", " << v << ") = " << lr_pts[6] << endl; 00202 cout << "dxdz SS(" << u << ", " << v << ") = " << sv_pts[6] << endl; 00203 cout << "d2y LR(" << u << ", " << v << ") = " << lr_pts[7] << endl; 00204 cout << "d2y SS(" << u << ", " << v << ") = " << sv_pts[7] << endl; 00205 cout << "dydz LR(" << u << ", " << v << ") = " << lr_pts[8] << endl; 00206 cout << "dydz SS(" << u << ", " << v << ") = " << sv_pts[8] << endl; 00207 cout << "d2z LR(" << u << ", " << v << ") = " << lr_pts[9] << endl; 00208 cout << "d2z SS(" << u << ", " << v << ") = " << sv_pts[9] << endl; 00209 00210 // collect results for summary at the end 00211 assertion_passed.push_back(correct); 00212 par_u_values.push_back(u); 00213 par_v_values.push_back(v); 00214 par_w_values.push_back(w); 00215 00216 if(correct) 00217 cout << "Parameter (" << u << ", " << v << ") evaluated OK\n"; 00218 else { 00219 cout << "ASSERTION FAILED at (" << u << ", " << v << ")\n"; 00220 } 00221 cout << endl; 00222 } 00223 } 00224 } 00225 00226 // display results 00227 cout << " ======= RESULT (VOLUME) SUMMARY ======== \n\n"; 00228 cout << "_____u____________v____________w_____________ASSERT__\n"; 00229 for(uint i=0; i<assertion_passed.size(); i++) { 00230 printf("%10.4g %10.4g %10.4g ", par_u_values[i], par_v_values[i], par_w_values[i]); 00231 if(assertion_passed[i]) 00232 cout << "OK\n"; 00233 else 00234 cout << "FAIL\n"; 00235 if(!assertion_passed[i]) 00236 oneFail = true; 00237 } 00238 00239 } 00240 cout << "----------------------------------------\n"; 00241 if(oneFail) 00242 cout << " test FAILED\n"; 00243 else 00244 cout << " all assertions passed\n"; 00245 cout << "========================================\n"; 00246 00247 }
1.7.6.1