1 /* 2 * Copyright (C) 1998, 2000-2007, 2010, 2011, 2012, 2013 SINTEF ICT, 3 * Applied Mathematics, Norway. 4 * 5 * Contact information: E-mail: tor.dokken@sintef.no 6 * SINTEF ICT, Department of Applied Mathematics, 7 * P.O. Box 124 Blindern, 8 * 0314 Oslo, Norway. 9 * 10 * This file is part of SISL. 11 * 12 * SISL is free software: you can redistribute it and/or modify 13 * it under the terms of the GNU Affero General Public License as 14 * published by the Free Software Foundation, either version 3 of the 15 * License, or (at your option) any later version. 16 * 17 * SISL is distributed in the hope that it will be useful, 18 * but WITHOUT ANY WARRANTY; without even the implied warranty of 19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 20 * GNU Affero General Public License for more details. 21 * 22 * You should have received a copy of the GNU Affero General Public 23 * License along with SISL. If not, see 24 * <http://www.gnu.org/licenses/>. 25 * 26 * In accordance with Section 7(b) of the GNU Affero General Public 27 * License, a covered work must retain the producer line in every data 28 * file that is created or manipulated using SISL. 29 * 30 * Other Usage 31 * You can be released from the requirements of the license by purchasing 32 * a commercial license. Buying such a license is mandatory as soon as you 33 * develop commercial activities involving the SISL library without 34 * disclosing the source code of your own applications. 35 * 36 * This file may be used in accordance with the terms contained in a 37 * written agreement between you and SINTEF ICT. 38 */ 39 40 #include "sisl-copyright.h" 41 42 /* 43 * 44 * $Id: s1990.c,v 1.5 2005-02-28 09:04:49 afr Exp $ 45 * 46 */ 47 48 49 #define S1990 50 51 #include "sislP.h" 52 /* 53 * Forward declarations. 54 * --------------------- 55 */ 56 57 #if defined(SISLNEEDPROTOTYPES) 58 static void s1990_s9edg(double [],double [],double [],double,double *, 59 int,int *); 60 /* 61 static void s1990_s9smooth(double [],int,int,int,double,double [],int *); 62 */ 63 #else 64 static void s1990_s9edg(); 65 /* static void s1990_s9smooth(); */ 66 #endif 67 68 #if defined(SISLNEEDPROTOTYPES) 69 void s1990(SISLSurf * ps,double aepsge,int * jstat)70 s1990(SISLSurf *ps,double aepsge,int *jstat) 71 #else 72 void s1990(ps,aepsge,jstat) 73 SISLSurf *ps; 74 double aepsge; 75 int *jstat; 76 #endif 77 /* 78 ********************************************************************* 79 * 80 ********************************************************************* 81 * 82 * PURPOSE : To make the orientation surface on the unit sphere to 83 * a b-spline surface, the surface is representated with 84 * a surrounding cone piced from the unit sphere. 85 * 86 * 87 * 88 * INPUT : ps - The original B-spline surface. 89 * aepsge - Geometry resolution. 90 * 91 * 92 * OUTPUT : jstat - status messages 93 * > 0 : warning 94 * = 0 : ok 95 * < 0 : error 96 * 97 * 98 * METHOD : We are making a cone surrounding the orientating surface 99 * on the unit sphere. The cone is representated with senter 100 * coordinates and an angle. The orientation is computed 101 * from aproximation of the normal to the surface. 102 * 103 * 104 * REFERENCES : 105 * 106 *- 107 * CALLS : 108 * 109 * WRITTEN BY : Arne Laksaa, SI, 89-01. 110 * UJK, Changed to accept equality between vertices. 111 * REWISED BY : Vibeke Skytt, SI, 91-02. 112 * UJK, SI, 91-10.Accepting surf's degenarated to a curve lying in 113 * a plane. Necessary for 2D ! 114 * Revised by : Christophe Rene Birkeland, SINTEF OSLO, June 1993. 115 * 116 ********************************************************************* 117 */ 118 { 119 int kpos = 0; /* Position of the error. */ 120 int kstat = 0; /* Local status variable. */ 121 int kfirst = 1; /* Flag to mark if the first patch is treating. */ 122 int kcount; /* Counts number of vanishing normals. */ 123 int kn1; /* Number of vertices of surface in 1. par. direction.*/ 124 int kn2; /* Number of vertices of surface in 2. par. direction.*/ 125 int kdim; /* Dimension of the space in which the objects lie. */ 126 int kdim4; /* Help variable to contain 4*kdim. */ 127 int kver,khor; /* The index to the vertice in the upper left corner 128 to the patch to treat. */ 129 int k1,k2,k3,k4; /* Control variables in loop. */ 130 int ki; /* Control variable in loop. */ 131 int lcone[4]; /* Flag telling if the cone has been generated. */ 132 double *t=SISL_NULL; /* Allocating t[5][kdim]. Five tangents around the 133 patch, the first and the last is the same. */ 134 double *tn; /* Allocating tn[4][kdim]. Four normals in the corner 135 of the patch. */ 136 double *tsen; /* Allocating tsen[4][kdim] for senter in edge cones. */ 137 double *ttan; /* Allocating ttan[kdim] for tangent on edges. */ 138 double tmax,tmin; /* Maximum and minimum coordinates to the narmals in 139 the first patch. */ 140 double tlen; /* The length of a vector. */ 141 double tnlen; /* The length of a normal vector. */ 142 double tang; /* An angle between two vectors. */ 143 double t1,t2; /* Help variables. */ 144 double sang[4]; /* Angel to the cones to edges. */ 145 double svec1[3]; /* Vectors used to determin degeneration. */ 146 double svec2[3]; /* Vectors used to determin degeneration. */ 147 double *scoef; /* Pointer to smoothed coefficient vector. */ 148 double slen[5]; /* Distances between coefficients. */ 149 double scorn[4]; /* Angle between derivatives in corner of patch. */ 150 double aepsge2 = min(aepsge, 1.0e-6); 151 152 /* Initiate output status */ 153 154 *jstat = 0; 155 156 /* Test if the surfaces already have been treated. */ 157 158 if (ps->pdir != SISL_NULL) goto out; 159 160 /* Initialate dimentions. */ 161 162 kdim = ps -> idim; 163 kn1 = ps -> in1; 164 kn2 = ps -> in2; 165 kdim4 = 4*kdim; 166 167 lcone[0] = 1; 168 lcone[1] = 1; 169 lcone[2] = 1; 170 lcone[3] = 1; 171 172 /*Make a new direction cone. */ 173 174 if ((ps->pdir = newdir(kdim)) == SISL_NULL) goto err101; 175 176 ps->pdir->aang = DZERO; 177 for (k1=0;k1<kdim;k1++) ps->pdir->ecoef[k1] = DZERO; 178 179 /* Allocate scratch for smoothed coefficients. */ 180 181 if ((ps->pdir->esmooth = newarray(kn1*kn2*kdim,DOUBLE)) == SISL_NULL) goto err101; 182 scoef = ps->pdir->esmooth; 183 184 /* Compute coefficients of smoothed curve. */ 185 186 /* s1990_s9smooth(ps->ecoef,kn1,kn2,kdim,aepsge,scoef,&kstat); 187 if (kstat < 0) goto error; */ 188 189 memcopy(scoef,ps->ecoef,kn1*kn2*kdim,DOUBLE); 190 191 /* Allocate local used matrices, t[5][kdim] and tn[4][kdim]. */ 192 193 if ((t = newarray(14*kdim,double)) == SISL_NULL) goto err101; 194 tn = t + 5*kdim; 195 tsen = tn + 4*kdim; 196 ttan = tsen + 4*kdim; 197 198 /* Here we are treating each patch in the control polygon separately.*/ 199 200 for (kver=0; kver < (kn2-1); kver++) 201 for (khor=0; khor < (kn1-1); khor++) 202 { 203 slen[0] = slen[1] = slen[2] = slen[3] = DZERO; 204 scorn[0] = scorn[1] = scorn[2] = scorn[3] = DZERO; 205 206 /* Here we make the tangents in each corner of the patch, 207 and in direction with the clock. The first and the last 208 vector contains both the first tangent. */ 209 210 k2 = (kver*kn1+khor)*kdim; 211 212 for (k1=0; k1 < kdim; k1++,k2++) 213 { 214 t[kdim+k1] = scoef[k2+kdim] - scoef[k2]; 215 t[2*kdim+k1] = scoef[k2+(kn1+1)*kdim]-scoef[k2+kdim]; 216 t[3*kdim+k1] = scoef[k2+kn1*kdim]-scoef[k2+(kn1+1)*kdim]; 217 t[kdim4+k1] = t[k1] = scoef[k2]-scoef[k2+kn1*kdim]; 218 219 slen[0] += t[k1]*t[k1]; 220 slen[1] += t[k1+kdim]*t[k1+kdim]; 221 slen[2] += t[k1+2*kdim]*t[k1+2*kdim]; 222 slen[3] += t[k1+3*kdim]*t[k1+3*kdim]; 223 } 224 slen[4] = slen[0] = sqrt(slen[0]); 225 slen[1] = sqrt(slen[1]); 226 slen[2] = sqrt(slen[2]); 227 slen[3] = sqrt(slen[3]); 228 229 scorn[0] = s6ang(t,t+kdim,kdim); 230 scorn[1] = s6ang(t+kdim,t+2*kdim,kdim); 231 scorn[2] = s6ang(t+2*kdim,t+3*kdim,kdim); 232 scorn[3] = s6ang(t+3*kdim,t,kdim); 233 234 /* If problems on edges is found we jump to the surface. */ 235 236 if (ps->pdir->igtpi > 0) goto next; 237 238 /* Computing cones of edges in ends of parameter two. */ 239 240 if (kver == 0) 241 { 242 if (lcone[0]) 243 { 244 /* First time to generate cone. */ 245 246 memcopy(tsen,t+kdim,kdim,DOUBLE); 247 tlen = slen[1]; 248 249 if (tlen > aepsge2) 250 { 251 for (k1=0; k1 < kdim; k1++) tsen[k1] /= tlen; 252 lcone[0] = 0; 253 sang[0] = (double)0; 254 } 255 } 256 else 257 { 258 /* Modify existing cone. */ 259 s1990_s9edg(t+(kdim),ttan,tsen,aepsge2,sang,kdim,&kstat); 260 261 if (kstat) ps->pdir->igtpi = 10; 262 } 263 } 264 if (kver == kn2-2) 265 { 266 if (lcone[1]) 267 { 268 /* First time to generate cone. */ 269 270 memcopy(tsen+kdim,t+3*kdim,kdim,DOUBLE); 271 tlen = slen[3]; 272 273 if (tlen > aepsge2) 274 { 275 for (k1=0; k1 < kdim; k1++) tsen[kdim+k1] /= tlen; 276 lcone[1] = 0; 277 sang[1] = (double)0; 278 } 279 } 280 else 281 { 282 s1990_s9edg(t+(3*kdim),ttan,tsen+kdim,aepsge2,sang+1,kdim,&kstat); 283 if (kstat) ps->pdir->igtpi = 10; 284 } 285 } 286 287 /* Computing cones of edges in ends of parameter one. */ 288 289 if (khor == 0) 290 { 291 if (lcone[2]) 292 /* First time to generate cone. */ 293 { 294 memcopy(tsen+2*kdim,t,kdim,DOUBLE); 295 tlen = slen[0]; 296 297 if (tlen > aepsge2) 298 { 299 for (k1=0; k1 < kdim; k1++) tsen[2*kdim+k1] /= tlen; 300 lcone[2] = 0; 301 sang[2] = (double)0; 302 } 303 } 304 else 305 { 306 s1990_s9edg(t,ttan,tsen+(2*kdim),aepsge2,sang+2,kdim,&kstat); 307 if (kstat) ps->pdir->igtpi = 10; 308 } 309 } 310 if (khor == kn1-2) 311 { 312 if (lcone[3]) 313 { 314 memcopy(tsen+3*kdim,t+2*kdim,kdim,DOUBLE); 315 tlen = slen[2]; 316 317 if (tlen > aepsge2) 318 { 319 for (k1=0; k1 < kdim; k1++) tsen[3*kdim+k1] /= tlen; 320 lcone[3] = 0; 321 sang[3] = (double)0; 322 } 323 } 324 else 325 { 326 s1990_s9edg(t+(2*kdim),ttan,tsen+(3*kdim),aepsge2,sang+3,kdim,&kstat); 327 if (kstat) ps->pdir->igtpi = 10; 328 } 329 } 330 331 next: 332 333 /* Here we makes the normales in each corner of the patch. 334 We are using a cross product between two tangents. 335 The normals is also normalized by deviding with its 336 own length. */ 337 338 for (kcount=0, ki=0, k1=0; k1 < kdim4; k1+=kdim, ki++) 339 { 340 for (tlen=DZERO,k2=0,k3=1,k4=2; k2 < kdim; k2++,k3++,k4++) 341 { 342 if(k3 == kdim) k3 = 0; 343 if(k4 == kdim) k4 = 0; 344 tn[k1+k2] = t[k1+k3]*t[k1+kdim+k4]-t[k1+k4]*t[k1+kdim+k3]; 345 346 tlen += tn[k1+k2]*tn[k1+k2]; 347 } 348 tlen = sqrt(tlen); 349 /* KYS 070494 : multiplied ANGULAR_TOLERANCE by 1.0e-2 */ 350 if (slen[ki]>aepsge2 && slen[ki+1]>aepsge2 && 351 scorn[ki] > 1.0e-2*ANGULAR_TOLERANCE) 352 for (k2=0; k2 < kdim; k2++) tn[k1+k2] /= tlen; 353 else 354 { 355 for (k2=0; k2 < kdim; k2++) tn[k1+k2] = ps->pdir->ecoef[k2]; 356 kcount++; 357 } 358 } 359 360 if (kcount == 4) continue; /* Degenerate control polygon patch */ 361 362 /* We are treating the first patch. */ 363 364 if (kfirst) 365 { 366 /* Computing the center coordinates of the cone.*/ 367 368 for (tlen=DZERO,k1=0; k1 < kdim; k1++) 369 { 370 tmin = (double)1.0; 371 tmax = - tmin; 372 for (k2=0; k2 < kdim4; k2+=kdim) 373 { 374 tmax = max(tn[k2+k1],tmax); 375 tmin = min(tn[k2+k1],tmin); 376 } 377 ps->pdir->ecoef[k1]=(tmax+tmin)/(double)2; 378 379 tlen += ps->pdir->ecoef[k1]*ps->pdir->ecoef[k1]; 380 } 381 tlen = sqrt(tlen); 382 if (tlen > DZERO) 383 for (k1=0; k1 < kdim; k1++) ps->pdir->ecoef[k1] /= tlen; 384 else 385 /* KYS 070494 : 'continue' replaced by the following block {} */ 386 /* There are nonzero normals pointing in 387 opposite directions, i.e. not simple case */ 388 { 389 if (khor <= kver) 390 ps->pdir->igtpi = 2; 391 else 392 ps->pdir->igtpi = 1; 393 ps->pdir->aang = PI; 394 goto out; 395 } 396 397 /* Computing the angle of the cone. */ 398 399 for (ps->pdir->aang=DZERO,k1=0; k1<kdim4; k1+=kdim) 400 { 401 for (tnlen=DZERO,tlen=DZERO,k2=0;k2<kdim;k2++) 402 { 403 tlen += ps->pdir->ecoef[k2]*tn[k1+k2]; 404 tnlen += tn[k1+k2]*tn[k1+k2]; 405 } 406 407 if (tlen >= DZERO) tlen = min((double)1.0,tlen); 408 else tlen = max((double)-1.0,tlen); 409 410 tlen = acos(tlen); 411 if (sqrt(tnlen) < aepsge2) tlen = DZERO; 412 413 ps->pdir->aang = max(ps->pdir->aang,tlen); 414 } 415 416 kfirst = 0; /* The first patch have been treated.*/ 417 } 418 else 419 for (k1=0; k1<kdim4; k1+=kdim) 420 { 421 /* Computing the angle beetween the senter of the cone 422 and the normal. */ 423 424 for (tnlen=DZERO,tang=DZERO,k2=0;k2<kdim;k2++) 425 { 426 tang += ps->pdir->ecoef[k2]*tn[k1+k2]; 427 tnlen += tn[k1+k2]*tn[k1+k2]; 428 } 429 430 if (tang >= DZERO) tang = MIN((double)1.0,tang); 431 else tang = MAX((double)-1.0,tang); 432 433 tang = acos(tang); 434 if (sqrt(tnlen) < aepsge2) tang = DZERO; 435 436 if (tang + ps->pdir->aang >= PI) 437 { 438 /* The angle is to great, give a meesage 439 how to subdivied and exit this function. */ 440 441 if (khor <= kver) 442 ps->pdir->igtpi = 2; 443 else 444 ps->pdir->igtpi = 1; 445 goto out; 446 } 447 else if (tang > ps->pdir->aang) 448 { 449 /* The normal is not inside the cone, than we 450 have to compute a new cone. */ 451 452 /* Computing the center coordinates.*/ 453 454 double sin_tang = sin(tang); /*@ hke */ 455 double delta = (tang - ps->pdir->aang)/2.0; /*@ hke */ 456 457 t1 = sin(delta)/sin_tang; /*@ hke */ 458 t2 = sin(tang - delta)/sin_tang; /*@ hke */ 459 460 /* 461 t1 = (tang - ps->pdir->aang)/((double)2*tang); 462 t2 = (double)1 - t1; 463 */ 464 465 for (tlen=DZERO,k2=0; k2<kdim; k2++) 466 { 467 ps->pdir->ecoef[k2] = 468 ps->pdir->ecoef[k2]*t2 + tn[k1+k2]*t1; 469 tlen += ps->pdir->ecoef[k2]*ps->pdir->ecoef[k2]; 470 } 471 tlen = sqrt(tlen); 472 473 for (k2=0; k2 < kdim; k2++) ps->pdir->ecoef[k2] /= tlen; 474 475 /* Computing the angle of the cone. */ 476 477 ps->pdir->aang = (tang + ps->pdir->aang)/(double)2; 478 } 479 } 480 481 if (ps->pdir->aang >= SIMPLECASE) 482 { 483 /* The angle is to great, give a meesage 484 how to subdivied and exit this function. */ 485 486 if (khor <= kver) 487 ps->pdir->igtpi = 20; 488 else 489 ps->pdir->igtpi = 10; 490 } 491 } 492 493 /* A final check if we have made a cone. */ 494 /* UJK, SI, 91-10, when 2D, return values from edge case */ 495 if (kfirst && kdim != 2) 496 { 497 /* No cone has been generated. We must examin if the surface is 498 degenerated to a point or line. */ 499 for (k1 = 1; k1 < kn1*kn2; k1++) 500 if (s6dist(scoef,scoef + (k1*kdim),kdim) >aepsge2) break; 501 502 if (k1 == kn1*kn2) 503 { 504 /* Degenerated to a point. */ 505 ps->pdir->igtpi = 0; 506 ps->pdir->aang = DZERO; 507 ps->pdir->ecoef[0] = (double) 1.0; 508 for (k1 = 1; k1 < kdim; k1++) ps->pdir->ecoef[k1] = DZERO; 509 } 510 else 511 { 512 s6diff(scoef,scoef + (k1*kdim),kdim,svec1); 513 514 for (k2 = k1 + 1; k2 < kn1*kn2; k2++) 515 if (s6dist(scoef,scoef + (k2*kdim),kdim) >aepsge2) 516 { 517 s6diff(scoef,scoef + (k2*kdim),kdim,svec2); 518 if (s6ang(svec1,svec2,kdim) > 1.0e-2*ANGULAR_TOLERANCE) break; 519 } 520 521 if (k2 == kn1*kn2) 522 { 523 /* Degenerated to a line. */ 524 ps->pdir->igtpi = 0; 525 ps->pdir->aang = DZERO; 526 ps->pdir->ecoef[0] = (double) 1.0; 527 for (k1 = 1; k1 < kdim; k1++) ps->pdir->ecoef[k1] = DZERO; 528 } 529 else 530 { 531 /* Three points describing a plane found, continue subdividing. */ 532 if (ps->et1[kn1] - ps->et1[ps->ik1-1] >= 533 ps->et2[kn2] - ps->et2[ps->ik2-1]) 534 ps->pdir->igtpi = 1; 535 else 536 ps->pdir->igtpi = 2; 537 } 538 } 539 } 540 541 /* success */ 542 543 goto out; 544 545 /* Error in space allacation. */ 546 547 err101: 548 *jstat = -101; 549 s6err("s1990",*jstat,kpos); 550 goto out; 551 552 /* Error in lower level routine. */ 553 554 /* error : 555 *jstat = kstat; 556 goto out; 557 */ 558 559 /* Free local used memory. */ 560 561 out: 562 if (t != SISL_NULL) freearray(t); 563 } 564 565 #if defined(SISLNEEDPROTOTYPES) 566 static void s1990_s9edg(double et[],double etan[],double esen[],double aepsge,double * cang,int idim,int * jstat)567 s1990_s9edg(double et[],double etan[],double esen[],double aepsge, 568 double *cang,int idim,int *jstat) 569 #else 570 static void s1990_s9edg(et,etan,esen,aepsge,cang,idim,jstat) 571 double et[]; 572 double etan[]; 573 double esen[]; 574 double aepsge; 575 double *cang; 576 int idim; 577 int *jstat; 578 #endif 579 /* 580 ********************************************************************* 581 * 582 ********************************************************************* 583 * 584 * PURPOSE : To make the orientation surface on the unit sphere to 585 * the edges of a b-spline surface, the surface is 586 * representated with a surrounding cone piced from 587 * the unit sphere. 588 * 589 * 590 * 591 * INPUT : et[] - The tangent vector to the edge. 592 * etan[] - The normalized tangent vector to the edge. 593 * idim - The dimention of the surface. 594 * 595 * 596 * INPUT/OUTPUT:esen[] - The senter vector of the cone. 597 * cang - The angel of the cone.. 598 * 599 * 600 * 601 * OUTPUT : jstat - status messages 602 * > 0 : warning 603 * = 0 : ok 604 * < 0 : error 605 * 606 * 607 * METHOD : 608 * 609 * 610 * REFERENCES : 611 * 612 *- 613 * CALLS : 614 * 615 * WRITTEN BY : Arne Laksaa, SI, 89-05. 616 * 617 ********************************************************************* 618 */ 619 { 620 int ki; 621 double tlen; 622 double tang; 623 double t1,t2; 624 625 626 /* Normalizing the tangent. */ 627 628 for (tlen = DZERO,ki=0; ki < idim; ki++) 629 { 630 etan[ki] = et[ki]; 631 tlen += etan[ki]*etan[ki]; 632 } 633 tlen = sqrt(tlen); 634 635 if (tlen > aepsge) 636 for (ki=0; ki < idim; ki++) etan[ki] /= tlen; 637 else 638 { 639 *jstat = 0; 640 goto out; 641 } 642 643 644 /* Computing the angle beetween the senter of the cone 645 and the tangent. */ 646 647 for (tang=DZERO,ki=0;ki<idim;ki++) 648 tang += esen[ki]*etan[ki]; 649 650 if (tang >= DZERO) tang = min((double)1.0,tang); 651 else tang = max((double)-1.0,tang); 652 653 tang = acos(tang); 654 655 656 if (tang + *cang >= PI) 657 { 658 /* The angle is to great, give a meesage 659 to subdivied and exit this function. */ 660 661 *jstat = 1; 662 goto out; 663 } 664 else if (tang > *cang) 665 { 666 /* The tangent is not inside the cone, and we 667 have to compute a new cone. */ 668 669 /* Computing the center coordinates.*/ 670 671 t1 = (tang - *cang)/((double)2*tang); 672 t2 = (double)1 - t1; 673 674 for (tlen=DZERO,ki=0; ki<idim; ki++) 675 { 676 esen[ki] = esen[ki]*t2 + etan[ki]*t1; 677 tlen += esen[ki]*esen[ki]; 678 } 679 tlen = sqrt(tlen); 680 681 if (tlen > DZERO) 682 for (ki=0; ki < idim; ki++) esen[ki] /= tlen; 683 else 684 { 685 /* Vi have to be aware of colapsed polygon. */ 686 687 *jstat = 1; 688 goto out; 689 } 690 691 /* Computing the angle of the cone. */ 692 693 *cang = (tang + *cang)/(double)2; 694 } 695 696 697 if (*cang >= SIMPLECASE) 698 { 699 /* The angle is to large, give a meesage 700 to subdivied and exit this function. */ 701 702 *jstat = 1; 703 goto out; 704 } 705 706 707 *jstat = 0; 708 709 out: ; 710 } 711 712 #if 0 713 #if defined(SISLNEEDPROTOTYPES) 714 static void 715 s1990_s9smooth(double ecoef1[],int in1,int in2,int idim, 716 double aepsge,double ecoef2[],int *jstat) 717 #else 718 static void s1990_s9smooth(ecoef1,in1,in2,idim,aepsge,ecoef2,jstat) 719 double ecoef1[]; 720 int in1; 721 int in2; 722 int idim; 723 double aepsge; 724 double ecoef2[]; 725 int *jstat; 726 #endif 727 /* 728 ********************************************************************* 729 * 730 ********************************************************************* 731 * 732 * PURPOSE : Perform noise filthering at the corners of the control 733 * polygon of a B-spline surface. 734 * 735 * 736 * 737 * INPUT : ecoef1 - Original coefficients of surface 738 * in1 - Number of coefficients in 1. par dir. 739 * in2 - Number of coefficients in 2. par dir. 740 * idim - Dimension of geometry space. 741 * aepsge - Geometry resolution. 742 * 743 * 744 * 745 * OUTPUT : ecoef2 - New coefficients after smoothing. 746 * jstat - status messages 747 * > 0 : warning 748 * = 0 : ok 749 * < 0 : error 750 * 751 * 752 * METHOD : 753 * 754 * 755 * 756 * REFERENCES : 757 * 758 *- 759 * CALLS : s6dplane - Distance to given plane. 760 * s6dline - Distance to given line. 761 * s6dist - Distance between two points. 762 * 763 * WRITTEN BY : Vibeke Skytt, SI, 91-02. 764 * 765 ********************************************************************* 766 */ 767 { 768 int kstat = 0; /* Local status variable. */ 769 int kn = MIN(in1/2,in2/2)+1; /* Maximum numbers of 770 coefficients to smooth. */ 771 int ki,kj,kh,kl; /* Counters. */ 772 int kc; /* Index of current corner. */ 773 int k1; /* Sign of change in 1. par dir */ 774 int k2; /* Sign of change in 2. par dir */ 775 int lcorn[4]; /* Indexes of corners. */ 776 int lsgn1[4]; /* Sign of changes in 1. par dir */ 777 int lsgn2[4]; /* Sign of changes in 2. par dir */ 778 double tdist; /* Distance to closest point in plane. */ 779 780 /* Set contents of arrays. */ 781 782 lcorn[0] = 0; 783 lcorn[1] = (in1-1)*idim; 784 lcorn[2] = (in1*in2-1)*idim; 785 lcorn[3] = in1*(in2-1)*idim; 786 787 lsgn1[0] = 1; 788 lsgn1[1] = -1; 789 lsgn1[2] = -1; 790 lsgn1[3] = 1; 791 792 lsgn2[0] = 1; 793 lsgn2[1] = 1; 794 lsgn2[2] = -1; 795 lsgn2[3] = -1; 796 797 /* Copy coefficients to output array. */ 798 799 memcopy(ecoef2,ecoef1,in1*in2*idim,DOUBLE); 800 801 /* For each corner, try to smooth the coefficients in the 802 neighbourhood of the corner. */ 803 804 for (ki=0; ki<4; ki++) 805 { 806 kc = lcorn[ki]; /* Index of current corner. */ 807 k1 = lsgn1[ki]; /* Sign change in 1. par dir. */ 808 k2 = lsgn2[ki]; /* Sign change in 2. par dir. */ 809 810 /* Try to smooth coefficients on center line. */ 811 812 for (kj=2; kj<kn; kj++) 813 { 814 if (s6dist(ecoef2+kc,ecoef2+kc+(k2*kj*in1+k1*kj)*idim, 815 idim) < aepsge) continue; 816 817 for (kh=1; kh<kj; kh++) 818 { 819 tdist = s6dline(ecoef2+kc,ecoef2+kc+(k2*kj*in1+k1*kj)*idim, 820 ecoef2+kc+(k2*kh*in1+k1*kh)*idim,idim,&kstat); 821 if (kstat < 0) goto error; 822 if (kstat || tdist >= aepsge) break; 823 } 824 if (kh < kj) break; 825 } 826 827 /* Perform smoothing. */ 828 829 kj--; 830 for (kh=1; kh<kj; kh++) 831 memcopy(ecoef2+kc+(k2*kh*in1+k1*kh)*idim,ecoef2+kc, 832 idim,DOUBLE); 833 834 /* Try to smooth coefficients on lower triangle. */ 835 836 for (kj=2; kj<kn; kj++) 837 { 838 for (kh=1; kh<kj; kh++) 839 { 840 for (kl=0; kl<kh; kl++) 841 { 842 tdist = s6dplane(ecoef2+kc,ecoef2+kc+k1*kj*idim, 843 ecoef2+kc+(k2*kj*in1+k1*kj)*idim, 844 ecoef2+kc+(k2*kl*in1+k1*kh)*idim, 845 idim,&kstat); 846 if (tdist >= aepsge) break; 847 } 848 if (tdist >= aepsge) break; 849 } 850 if (kh < kj) break; 851 } 852 853 /* Perform smoothing. */ 854 855 kj--; 856 for (kh=1; kh<kj; kh++) 857 for (kl=0; kl<kh; kl++) 858 memcopy(ecoef2+kc+(k2*kl*in1+k1*kh)*idim,ecoef2+kc, 859 idim,DOUBLE); 860 861 /* Try to smooth coefficients on upper triangle. */ 862 863 for (kj=2; kj<kn; kj++) 864 { 865 for (kh=0; kh<kj; kh++) 866 { 867 for (kl=kh+1; kl<kj; kl++) 868 { 869 tdist = s6dplane(ecoef2+kc,ecoef2+kc+k2*kj*in1*idim, 870 ecoef2+kc+(k2*kj*in1+k1*kj)*idim, 871 ecoef2+kc+(k2*kl*in1+k1*kh)*idim, 872 idim,&kstat); 873 if (tdist >= aepsge) break; 874 } 875 if (tdist >= aepsge) break; 876 } 877 if (kh < kj) break; 878 } 879 880 /* Perform smoothing. */ 881 882 kj--; 883 for (kh=0; kh<kj; kh++) 884 for (kl=kh+1; kl<kj; kl++) 885 memcopy(ecoef2+kc+(k2*kl*in1+k1*kh)*idim,ecoef2+kc, 886 idim,DOUBLE); 887 } 888 889 /* Smoothing performed. */ 890 *jstat = 0; 891 goto out; 892 893 /* Error in lower level routine. */ 894 895 error : *jstat = kstat; 896 goto out; 897 898 out : 899 900 return; 901 } 902 903 #endif /* if 0 */ 904