1 #include <U2Core/disable-warnings.h> 2 U2_DISABLE_WARNINGS 3 4 // -*- Mode: C++; tab-width: 2; -*- 5 // vi: set ts=2: 6 // 7 8 //# include <BALL/KERNEL/atom.h> 9 //# include <BALL/KERNEL/molecule.h> 10 //# include <BALL/KERNEL/system.h> 11 //# include <BALL/KERNEL/PTE.h> 12 //# include <BALL/DATATYPE/string.h> 13 14 #include <BALL/STRUCTURE/solventAccessibleSurface.h> 15 #include <BALL/STRUCTURE/triangle.h> 16 #include <BALL/STRUCTURE/triangleEdge.h> 17 #include <BALL/STRUCTURE/trianglePoint.h> 18 #include <BALL/STRUCTURE/triangulatedSAS.h> 19 #include <BALL/STRUCTURE/triangulatedSurface.h> 20 #include <BALL/MATHS/angle.h> 21 #include <BALL/MATHS/matrix44.h> 22 #include <BALL/MATHS/plane3.h> 23 #include <BALL/MATHS/quaternion.h> 24 #include <BALL/MATHS/vector3.h> 25 #include <BALL/MATHS/vector4.h> 26 27 #include <list> 28 #include <vector> 29 30 namespace BALL 31 { 32 PartitionOfCircle(const TCircle3<double> & circle,std::list<Vector3> & partition)33 void PartitionOfCircle(const TCircle3<double>& circle, std::list<Vector3>& partition) 34 { 35 Size density = 64; 36 Vector3 center(circle.p.x,circle.p.y,circle.p.z); 37 Vector3 normal(circle.n.x,circle.n.y,circle.n.z); 38 Vector4 p; 39 p.set(normal.y,-normal.x,0.0,0.0); 40 if (p == Vector4::getZero()) 41 { 42 p.set(normal.z,0.0,-normal.x,0.0); 43 } 44 p.normalize(); 45 p *= circle.radius; 46 Quaternion rotate(normal,Angle(Constants::PI/density,true)); 47 Matrix4x4 rotation; 48 rotate.getRotationMatrix(rotation); 49 partition.push_back(Vector3(p.x,p.y,p.z)+center); 50 for (Size i = 0; i < 2*density+1; i++) 51 { 52 p = rotation*p; 53 partition.push_back(Vector3(p.x,p.y,p.z)+center); 54 } 55 } 56 TriangulatedSAS()57 TriangulatedSAS::TriangulatedSAS() 58 : TriangulatedSurface(), 59 sas_(NULL), 60 density_(4.5) 61 { 62 } 63 TriangulatedSAS(const TriangulatedSAS & surface,bool deep)64 TriangulatedSAS::TriangulatedSAS(const TriangulatedSAS& surface, bool deep) 65 : TriangulatedSurface(surface,deep), 66 sas_(surface.sas_), 67 density_(surface.density_) 68 { 69 } 70 TriangulatedSAS(SolventAccessibleSurface * sas,const double & density)71 TriangulatedSAS::TriangulatedSAS(SolventAccessibleSurface* sas, const double& density) 72 : TriangulatedSurface(), 73 sas_(sas), 74 density_(density) 75 { 76 } 77 ~TriangulatedSAS()78 TriangulatedSAS::~TriangulatedSAS() 79 { 80 } 81 set(const TriangulatedSAS & surface,bool)82 void TriangulatedSAS::set(const TriangulatedSAS& surface, bool) 83 { 84 if (this != &surface) 85 { 86 TriangulatedSurface::set(surface); 87 sas_ = surface.sas_; 88 density_ = surface.density_; 89 } 90 } 91 operator =(const TriangulatedSAS & surface)92 TriangulatedSAS& TriangulatedSAS::operator = (const TriangulatedSAS& surface) 93 { 94 if (this != &surface) 95 { 96 TriangulatedSurface::operator = (surface); 97 sas_ = surface.sas_; 98 density_ = surface.density_; 99 } 100 101 return *this; 102 } 103 setDensity(const double & density)104 void TriangulatedSAS::setDensity(const double& density) 105 { 106 density_ = density; 107 } 108 getDensity() const109 double TriangulatedSAS::getDensity() const 110 { 111 return density_; 112 } 113 compute()114 void TriangulatedSAS::compute() 115 { 116 SASTriangulator sast(this); 117 sast.run(); 118 } 119 120 ////////////////////////////////////////////////// 121 SASTriangulator()122 SASTriangulator::SASTriangulator() 123 : tsas_(NULL), 124 sqrt_density_(0.0), 125 edge_(), 126 template_spheres_() 127 { 128 } 129 SASTriangulator(TriangulatedSAS * tsas)130 SASTriangulator::SASTriangulator(TriangulatedSAS* tsas) 131 : tsas_(tsas), 132 sqrt_density_(sqrt(tsas->density_)), 133 edge_(tsas_->sas_->number_of_edges_), 134 template_spheres_() 135 { 136 } 137 ~SASTriangulator()138 SASTriangulator::~SASTriangulator() 139 { 140 } 141 run()142 void SASTriangulator::run() 143 { 144 // build template spheres with different densities and outside normal vectors 145 buildTemplateSpheres(); 146 147 // use these to triangulate the sas faces 148 for (Position i = 0; i < tsas_->sas_->number_of_faces_; i++) 149 { 150 triangulateFace(tsas_->sas_->faces_[i]); 151 } 152 } 153 triangulateFace(SASFace * face)154 void SASTriangulator::triangulateFace(SASFace* face) 155 { 156 // store the planes of the SAS edges 157 std::list< std::pair<TPlane3<double>, double> > planes; 158 createPlanes(face,planes); 159 160 // get template sphere 161 TSphere3<double>* sphere = &face->sphere_; 162 HashMap<Size,TriangulatedSurface>::ConstIterator s; 163 164 s = template_spheres_.find(numberOfRefinements(tsas_->density_, sphere->radius)); 165 166 TriangulatedSurface part = s->second; 167 part.blowUp(sphere->radius); 168 part.shift(sphere->p); 169 170 // tag inside points 171 tagPoints(part, planes); 172 173 // remove inside triangles and then isolated edges and points 174 removeInsideTriangles(part); 175 part.deleteIsolatedEdges(); 176 part.deleteIsolatedPoints(); 177 178 // create HashGrid of triangle points 179 // HashGrid3<TrianglePoint*> point_grid(createHashGrid(part)); 180 // create points on the cutting planes 181 // createPoints(part,planes,point_grid); 182 // create new triangles 183 // createNewTriangles(part,point_grid); 184 // remove isolated edges and points 185 // part.deleteIsolatedEdges(); 186 // part.deleteIsolatedPoints(); 187 // join the triangulation of this sas face with the sas 188 tsas_->join(part); 189 } 190 createPlanes(SASFace * face,std::list<std::pair<TPlane3<double>,double>> & planes)191 void SASTriangulator::createPlanes(SASFace* face, 192 std::list< std::pair<TPlane3<double>, double> >& planes) 193 { 194 std::pair<TPlane3<double>, double> plane; 195 std::list<SASEdge*>::iterator edge; 196 std::list<bool>::iterator o = face->orientation_.begin(); 197 for (edge = face->edge_.begin(); edge != face->edge_.end(); edge++) 198 { 199 plane.first.p = (*edge)->circle_.p; 200 if (*o) 201 { 202 plane.first.n = (*edge)->circle_.n; 203 } 204 else 205 { 206 plane.first.n = -(*edge)->circle_.n; 207 } 208 plane.second = plane.first.n*plane.first.p; 209 planes.push_back(plane); 210 o++; 211 } 212 } 213 tagPoints(TriangulatedSurface & part,const std::list<std::pair<TPlane3<double>,double>> & planes)214 void SASTriangulator::tagPoints(TriangulatedSurface& part, 215 const std::list< std::pair<TPlane3<double>, double> >& planes) 216 { 217 for (TriangulatedSurface::PointIterator p = part.beginPoint(); 218 p != part.endPoint(); 219 p++) 220 { 221 (*p)->index_ = 0; 222 std::list< std::pair<TPlane3<double>, double> >::const_iterator plane = planes.begin(); 223 for (; plane != planes.end(); ++plane) 224 { 225 if (Maths::isLessOrEqual(plane->first.n*(*p)->point_, plane->second)) 226 { 227 (*p)->index_ = 1; 228 break; 229 } 230 } 231 } 232 } 233 removeInsideTriangles(TriangulatedSurface & part)234 void SASTriangulator::removeInsideTriangles(TriangulatedSurface& part) 235 { 236 TriangulatedSurface::TriangleIterator t1, t2; 237 238 t1 = part.beginTriangle(); 239 while (t1 != part.endTriangle()) 240 { 241 if ((*t1)->vertex_[0]->index_ + 242 (*t1)->vertex_[1]->index_ + 243 (*t1)->vertex_[2]->index_ == 3) 244 { 245 t2 = t1; 246 t2++; 247 if (t2 == part.endTriangle()) 248 { 249 part.remove(t1); 250 t1 = part.endTriangle(); 251 } 252 else 253 { 254 part.remove(t1); 255 t1 = t2; 256 } 257 } 258 else 259 { 260 t1++; 261 } 262 } 263 } 264 createHashGrid(const TriangulatedSurface & part)265 HashGrid3<TrianglePoint*> SASTriangulator::createHashGrid(const TriangulatedSurface& part) 266 { 267 double x_min = (*(part.beginPoint()))->point_.x; 268 double y_min = (*(part.beginPoint()))->point_.y; 269 double z_min = (*(part.beginPoint()))->point_.z; 270 double x_max = (*(part.beginPoint()))->point_.x; 271 double y_max = (*(part.beginPoint()))->point_.y; 272 double z_max = (*(part.beginPoint()))->point_.z; 273 TriangulatedSurface::ConstPointIterator p = part.beginPoint()++; 274 while (p != part.endPoint()) 275 { 276 x_min = Maths::min(x_min,(*p)->point_.x); 277 y_min = Maths::min(y_min,(*p)->point_.y); 278 z_min = Maths::min(z_min,(*p)->point_.z); 279 x_max = Maths::max(x_min,(*p)->point_.x); 280 y_max = Maths::max(y_max,(*p)->point_.y); 281 z_max = Maths::max(z_max,(*p)->point_.z); 282 p++; 283 } 284 //double dist = ses_->reduced_surface_->r_max_/3; 285 double dist = 1.0; 286 Position nx = (Position)((x_max-x_min)/dist+5); 287 Position ny = (Position)((y_max-y_min)/dist+5); 288 Position nz = (Position)((z_max-z_min)/dist+5); 289 Vector3 origin(x_min-2*dist,y_min-2*dist,z_min-2*dist); 290 HashGrid3<TrianglePoint*> grid = HashGrid3<TrianglePoint*>(origin,nx,ny,nz,dist); 291 p = part.beginPoint()++; 292 while (p != part.endPoint()) 293 { 294 grid.insert(Vector3((*p)->point_.x,(*p)->point_.y,(*p)->point_.z),*p); 295 p++; 296 } 297 return grid; 298 } 299 createPoints(TriangulatedSurface & part,const std::list<std::pair<TPlane3<double>,double>> & planes,HashGrid3<TrianglePoint * > & grid)300 void SASTriangulator::createPoints(TriangulatedSurface& part, 301 const std::list< std::pair<TPlane3<double>,double> >& planes, 302 HashGrid3<TrianglePoint*>& grid) 303 { 304 TriangulatedSurface::EdgeIterator edge = part.beginEdge(); 305 std::list< std::pair<TPlane3<double>,double> >::const_iterator plane; 306 while (edge != part.endEdge()) 307 { 308 if ((*edge)->vertex_[0]->index_+(*edge)->vertex_[1]->index_ == 1) 309 { 310 // edge intersects one of the cutting planes 311 TrianglePoint* v1; 312 TrianglePoint* v2; 313 if ((*edge)->vertex_[0]->index_ == 0) 314 { 315 v1 = (*edge)->vertex_[0]; 316 v2 = (*edge)->vertex_[1]; 317 } 318 else 319 { 320 v1 = (*edge)->vertex_[1]; 321 v2 = (*edge)->vertex_[0]; 322 } 323 TVector3<double> diff(v2->point_-v1->point_); 324 TVector3<double> point(v2->point_); 325 double min = 1.0; 326 double div; 327 Size counter = 0; 328 for (plane = planes.begin(); plane != planes.end(); plane++) 329 { 330 // intersecting point = v1 + lambda * (v2-v1) 331 div = plane->first.n*diff; 332 if (Maths::isNotZero(div)) 333 { 334 double lambda = (plane->second-(plane->first.n*v1->point_))/div; 335 if (Maths::isGreaterOrEqual(lambda,0) && 336 (Maths::isLessOrEqual(lambda,min))) 337 { 338 min = lambda; 339 point.set(v1->point_+(lambda*diff)); 340 (*edge)->index_ = counter; 341 } 342 } 343 counter++; 344 } 345 v2->edges_.erase(*edge); 346 TrianglePoint* new_point = vertexExists(point,grid); 347 if (new_point == NULL) 348 { 349 new_point = new TrianglePoint; 350 new_point->point_ = point; 351 new_point->index_ = -1; 352 new_point->edges_.insert(*edge), 353 part.insert(new_point); 354 grid.insert(Vector3(point.x,point.y,point.z),new_point); 355 } 356 if ((*edge)->vertex_[0] == v1) 357 { 358 (*edge)->vertex_[1] = new_point; 359 } 360 else 361 { 362 (*edge)->vertex_[0] = new_point; 363 } 364 } 365 else 366 { 367 (*edge)->index_ = -1; 368 } 369 edge++; 370 } 371 } 372 createNewTriangles(TriangulatedSurface & part,HashGrid3<TrianglePoint * > & grid)373 void SASTriangulator::createNewTriangles(TriangulatedSurface& part, HashGrid3<TrianglePoint*>& grid) 374 { 375 TriangulatedSurface::TriangleIterator t = part.beginTriangle(); 376 while (t != part.endTriangle()) 377 { 378 Position type = 0; 379 if ((*t)->vertex_[0]->index_ == 1) 380 { 381 type++; 382 } 383 if ((*t)->vertex_[1]->index_ == 1) 384 { 385 type += 2; 386 } 387 if ((*t)->vertex_[2]->index_ == 1) 388 { 389 type += 4; 390 } 391 switch (type) 392 { 393 case 0 : break; // triangle not intersected 394 case 1 : onePointOutside(0,*t,part,grid); break; 395 case 2 : onePointOutside(1,*t,part,grid); break; 396 case 3 : twoPointsOutside(0,1,*t,part,grid); break; 397 case 4 : onePointOutside(2,*t,part,grid); break; 398 case 5 : twoPointsOutside(2,0,*t,part,grid); break; 399 case 6 : twoPointsOutside(1,2,*t,part,grid); break; 400 case 7 : break; // should never happen 401 } 402 t++; 403 } 404 } 405 onePointOutside(Index outside,Triangle * t,TriangulatedSurface & part,HashGrid3<TrianglePoint * > & grid)406 void SASTriangulator::onePointOutside(Index outside, Triangle* t, 407 TriangulatedSurface& part, HashGrid3<TrianglePoint*>& grid) 408 { 409 // get the relative indices of the intersected edges 410 Position edge[3]; 411 Position i = 0; 412 for (Position j = 0; j < 3; j++) 413 { 414 if (t->edge_[j]->index_ != -1) 415 { 416 edge[i] = j; 417 i++; 418 } 419 else 420 { 421 edge[2] = j; 422 } 423 } 424 // v1 = new vertex of the first intersected edge 425 // v2 = old vertex of the first intersected edge 426 // v3 = new vertex of the second intersected edge 427 // v4 = old vertex of the second intersected edge 428 TrianglePoint* v1; 429 TrianglePoint* v2; 430 TrianglePoint* v3; 431 TrianglePoint* v4; 432 Position test = ((t->edge_[edge[0]]->vertex_[0]->index_ == -1) ? 0 : 1); 433 v1 = t->edge_[edge[0]]->vertex_[test]; 434 v2 = t->edge_[edge[0]]->vertex_[1-test]; 435 test = ((t->edge_[edge[1]]->vertex_[0]->index_ == -1) ? 0 : 1); 436 v3 = t->edge_[edge[1]]->vertex_[test]; 437 v4 = t->edge_[edge[1]]->vertex_[1-test]; 438 439 // get the relative index of v4 which is needed 440 // to compute the correct orientation of the first new triangle 441 Index index = 0; 442 for (Position j = 0; j < 3; j++) 443 { 444 if (t->vertex_[j] == v4) 445 { 446 index = j; 447 } 448 } 449 // resize the original triangle 450 t->vertex_[outside]->faces_.erase(t); 451 t->vertex_[outside] = v1; 452 v1->faces_.insert(t); 453 // create a new triangle 454 Triangle* new_triangle = new Triangle; 455 new_triangle->vertex_[0] = v1; 456 if ((outside-index == 1) || (outside-index == -2)) 457 { 458 new_triangle->vertex_[1] = v4; 459 new_triangle->vertex_[2] = v3; 460 } 461 else 462 { 463 new_triangle->vertex_[1] = v3; 464 new_triangle->vertex_[2] = v4; 465 } 466 v1->faces_.insert(new_triangle); 467 v3->faces_.insert(new_triangle); 468 v4->faces_.insert(new_triangle); 469 part.insert(new_triangle); 470 // does the triangle intersect two different SAS edges? 471 if (t->edge_[edge[0]]->index_ != t->edge_[edge[1]]->index_) 472 { 473 // ACHTUNG: pos MUSS NOCH BERECHNET WERDEN !!! 474 TVector3<double> pos(v1->point_); 475 TrianglePoint* new_point = vertexExists(pos,grid); 476 if (new_point == NULL) 477 { 478 new_point = new TrianglePoint; 479 new_point->point_ = pos; 480 new_point->index_ = -1; 481 part.insert(new_point); 482 grid.insert(Vector3(pos.x,pos.y,pos.z),new_point); 483 } 484 Triangle* third_triangle = new Triangle; 485 third_triangle->vertex_[0] = v1; 486 if ((outside-index == 1) || (outside-index == -2)) 487 { 488 third_triangle->vertex_[1] = v3; 489 third_triangle->vertex_[2] = new_point; 490 } 491 else 492 { 493 third_triangle->vertex_[1] = new_point; 494 third_triangle->vertex_[2] = v3; 495 } 496 v1->faces_.insert(third_triangle); 497 v3->faces_.insert(third_triangle); 498 new_point->faces_.insert(third_triangle); 499 part.insert(third_triangle); 500 } 501 } 502 twoPointsOutside(Position outside1,Position outside2,Triangle * t,TriangulatedSurface & part,HashGrid3<TrianglePoint * > & grid)503 void SASTriangulator::twoPointsOutside(Position outside1, Position outside2, 504 Triangle* t, TriangulatedSurface& part, 505 HashGrid3<TrianglePoint*>& grid) 506 { 507 // get the relative indices of the intersected edges 508 Position edge[3]; 509 Position i = 0; 510 for (Position j = 0; j < 3; j++) 511 { 512 if (t->edge_[j]->index_ != -1) 513 { 514 edge[i] = j; 515 i++; 516 } 517 else 518 { 519 edge[2] = j; 520 } 521 } 522 // v1 = new vertex of the first intersected edge 523 // v2 = old vertex of the first intersected edge 524 // v3 = new vertex of the second intersected edge 525 TrianglePoint* v1; 526 TrianglePoint* v2; 527 TrianglePoint* v3; 528 Position test = ((t->edge_[edge[0]]->vertex_[0]->index_ == -1) ? 0 : 1); 529 v1 = t->edge_[edge[0]]->vertex_[test]; 530 v2 = t->edge_[edge[0]]->vertex_[1-test]; 531 test = ((t->edge_[edge[1]]->vertex_[0]->index_ == -1) ? 0 : 1); 532 v3 = t->edge_[edge[1]]->vertex_[test]; 533 // resize the original triangle 534 t->vertex_[outside1]->faces_.erase(t); 535 t->vertex_[outside2]->faces_.erase(t); 536 TLine3<double> line(v1->point_,v2->point_,TLine3<double>::FORM__TWO_POINTS); 537 bool one = line.has(t->vertex_[outside1]->point_); 538 if (one) 539 { 540 t->vertex_[outside1] = v1; 541 t->vertex_[outside2] = v3; 542 } 543 else 544 { 545 t->vertex_[outside1] = v3; 546 t->vertex_[outside2] = v1; 547 } 548 v1->faces_.insert(t); 549 v2->faces_.insert(t); 550 // does the triangle intersect two different SAS edges? 551 if (t->edge_[edge[0]]->index_ != t->edge_[edge[1]]->index_) 552 { 553 // ACHTUNG: pos MUSS NOCH BERECHNET WERDEN !!! 554 TVector3<double> pos(v1->point_); 555 TrianglePoint* new_point = vertexExists(pos,grid); 556 if (new_point == NULL) 557 { 558 new_point = new TrianglePoint; 559 new_point->point_ = pos; 560 new_point->index_ = -1; 561 part.insert(new_point); 562 grid.insert(Vector3(pos.x,pos.y,pos.z),new_point); 563 } 564 Triangle* new_triangle = new Triangle; 565 new_triangle->vertex_[0] = t->vertex_[outside2]; 566 new_triangle->vertex_[1] = t->vertex_[outside1]; 567 new_triangle->vertex_[2] = new_point; 568 v1->faces_.insert(new_triangle); 569 v3->faces_.insert(new_triangle); 570 new_point->faces_.insert(new_triangle); 571 part.insert(new_triangle); 572 } 573 } 574 vertexExists(const TVector3<double> & point,HashGrid3<TrianglePoint * > & grid)575 TrianglePoint* SASTriangulator::vertexExists(const TVector3<double>& point, 576 HashGrid3<TrianglePoint*>& grid) 577 { 578 double epsilon = Constants::EPSILON; 579 Constants::EPSILON = 0.001; 580 Vector3 p(point.x,point.y,point.z); 581 HashGridBox3<TrianglePoint*>* box = grid.getBox(p); 582 HashGridBox3<TrianglePoint*>::ConstBoxIterator b; 583 HashGridBox3<TrianglePoint*>::ConstDataIterator d; 584 if (box != NULL) 585 { 586 for (b = box->beginBox(); b != box->endBox(); b++) 587 { 588 for (d = b->beginData(); d != b->endData(); d++) 589 { 590 if ((*d)->point_ == point) 591 { 592 Constants::EPSILON = epsilon; 593 return *d; 594 } 595 } 596 } 597 } 598 Constants::EPSILON = epsilon; 599 return NULL; 600 } 601 numberOfRefinements(const double & density,const double & radius)602 Size SASTriangulator::numberOfRefinements(const double& density, const double& radius) 603 { 604 double test0 = (4.0*density*Constants::PI*radius*radius-12.0)/30.0; 605 Size n = 0; 606 if (Maths::isGreaterOrEqual(test0,0.0)) 607 { 608 double test1 = 1; 609 double test2 = 1; 610 while (Maths::isLess(test2,test0)) 611 { 612 test1 = test2; 613 test2 *= 4; 614 n++; 615 } 616 if (Maths::isLess(test2-test0,test0-test1)) 617 { 618 n++; 619 } 620 } 621 if (n > 4) 622 { 623 n = 4; 624 } 625 return n; 626 } 627 buildTemplateSpheres()628 void SASTriangulator::buildTemplateSpheres() 629 { 630 TriangulatedSphere sphere; 631 sphere.icosaeder(true); 632 sphere.setIndices(); 633 template_spheres_[0] = sphere; 634 sphere.refine(1,true); 635 sphere.setIndices(); 636 template_spheres_[1] = sphere; 637 sphere.refine(1,true); 638 sphere.setIndices(); 639 template_spheres_[2] = sphere; 640 sphere.refine(1,true); 641 sphere.setIndices(); 642 template_spheres_[3] = sphere; 643 sphere.refine(1,true); 644 sphere.setIndices(); 645 template_spheres_[4] = sphere; 646 } 647 648 } //namespace BALL 649