Ignore:
Timestamp:
Jun 12, 2008, 7:14:46 PM (17 years ago)
Author:
Frederik Heber <heber@…>
Children:
7fcea6
Parents:
9185c8
Message:

VolumeOfConvexEnvelope: Works!

VolumeOfConvexEnvelope has been analysed into various smaller functions and approach is working.
two new files: boundary.?pp
various new functions:
class Tesselation with AddPoint(), TesselateOnBoundary() and GuessStartingTriangle() does the actual tesselation
CreateClustersinWater() will create the repetition of the cluster with correct spacing (unfinished).
GetDiametersOfCluster() calculate the greatest diameter in projection per axis
GetBoundaryPoints() gets the boundary on the convex envelope by projection for a molecular cluster
GetCommonEndpoint() finds the endpoint two lines are sharing

File:
1 edited

Legend:

Unmodified
Added
Removed

  • molecuilder/src/molecules.cpp

    r9185c8 re292c10  
    13701370  return No;
    13711371};
    1372 
    1373 /** Determines the volume of a cluster.
    1374  * Determines first the convex envelope, then tesselates it and calculates its volume.
    1375  * \param *out output stream for debugging
    1376  * \param IsAngstroem for the units on output
    1377  */
    1378 double molecule::VolumeOfConvexEnvelope(ofstream *out, bool IsAngstroem)
    1379 {
    1380         atom *Walker = NULL;
    1381        
    1382         // 1. calculate center of gravity
    1383         *out << endl;
    1384         vector *CenterOfGravity = DetermineCenterOfGravity(out);
    1385        
    1386         // 2. translate all points into CoG
    1387         *out << Verbose(1) << "Translating system to Center of Gravity." << endl;
    1388         Walker=start;
    1389         while (Walker->next != end) {
    1390           Walker = Walker->next;
    1391           Walker->x.Translate(CenterOfGravity);
    1392         }
    1393         // 2. calculate distance of each atom and sort into hash table
    1394 //      map<double, int> DistanceSet;
    1395 //     
    1396 //      Walker = start;
    1397 //      while (Walker->next != end) {
    1398 //              Walker = Walker->next;
    1399 //              DistanceSet.insert( pair<double, int>(ptr->x.distance(CenterOfGravity), ptr->nr) );
    1400 //      }
    1401        
    1402         // 3. Find all points on the boundary
    1403   *out << Verbose(1) << "Finding all boundary points." << endl;
    1404         Boundaries BoundaryPoints[NDIM];        // first is alpha, second is (r, nr)
    1405         BoundariesTestPair BoundaryTestPair;
    1406         vector AxisVector, AngleReferenceVector, AngleReferenceNormalVector;
    1407         double radius, angle;
    1408         // 3a. Go through every axis
    1409         for (int axis=0; axis<NDIM; axis++)  {
    1410                 AxisVector.Zero();
    1411     AngleReferenceVector.Zero();
    1412     AngleReferenceNormalVector.Zero();
    1413                 AxisVector.x[axis] = 1.;
    1414                 AngleReferenceVector.x[(axis+1)%NDIM] = 1.;
    1415                 AngleReferenceNormalVector.x[(axis+2)%NDIM] = 1.;
    1416 //              *out << Verbose(1) << "Axisvector is ";
    1417 //              AxisVector.Output(out);
    1418 //              *out << " and AngleReferenceVector is ";
    1419 //              AngleReferenceVector.Output(out);
    1420 //              *out << "." << endl;
    1421 //    *out << " and AngleReferenceNormalVector is ";
    1422 //    AngleReferenceNormalVector.Output(out);
    1423 //    *out << "." << endl;
    1424                 // 3b. construct set of all points, transformed into cylindrical system and with left and right neighbours
    1425                 Walker = start;
    1426                 while (Walker->next != end) {
    1427                         Walker = Walker->next;
    1428                         vector ProjectedVector;
    1429                         ProjectedVector.CopyVector(&Walker->x);
    1430                         ProjectedVector.ProjectOntoPlane(&AxisVector);
    1431                   // correct for negative side
    1432                   //if (Projection(y) < 0)
    1433                     //angle = 2.*M_PI - angle;
    1434                         radius = ProjectedVector.Norm();
    1435                         if (fabs(radius) > MYEPSILON)
    1436                           angle = ProjectedVector.Angle(&AngleReferenceVector);
    1437                         else
    1438                           angle = 0.;  // otherwise it's a vector in Axis Direction and unimportant for boundary issues
    1439                          
    1440                         //*out << "Checking sign in quadrant : " << ProjectedVector.Projection(&AngleReferenceNormalVector) << "." << endl;
    1441                         if (ProjectedVector.Projection(&AngleReferenceNormalVector) > 0) {
    1442                           angle = 2.*M_PI - angle;
    1443                         }
    1444                         *out << Verbose(2) << "Inserting " << *Walker << ": (r, alpha) = (" << radius << "," << angle << "): ";
    1445                         ProjectedVector.Output(out);
    1446                         *out << endl;
    1447                         BoundaryTestPair = BoundaryPoints[axis].insert( BoundariesPair (angle, DistanceNrPair (radius, Walker) ) );
    1448                         if (BoundaryTestPair.second) { // successfully inserted
    1449                         } else { // same point exists, check first r, then distance of original vectors to center of gravity
    1450                                 *out << Verbose(2) << "Encountered two vectors whose projection onto axis " << axis << " is equal: " << endl;
    1451                                 *out << Verbose(2) << "Present vector: ";
    1452                                 BoundaryTestPair.first->second.second->x.Output(out);
    1453                                 *out << endl;
    1454                                 *out << Verbose(2) << "New vector: ";
    1455                                 Walker->x.Output(out);
    1456                                 *out << endl;
    1457                                 double tmp = ProjectedVector.Norm();
    1458                                 if (tmp > BoundaryTestPair.first->second.first) {
    1459                                         BoundaryTestPair.first->second.first = tmp;
    1460                                         BoundaryTestPair.first->second.second = Walker;
    1461                                         *out << Verbose(2) << "Keeping new vector." << endl;
    1462                                 } else if (tmp == BoundaryTestPair.first->second.first) {
    1463                                         if (BoundaryTestPair.first->second.second->x.ScalarProduct(&BoundaryTestPair.first->second.second->x) < Walker->x.ScalarProduct(&Walker->x)) { // Norm() does a sqrt, which makes it a lot slower
    1464                                                 BoundaryTestPair.first->second.second = Walker;
    1465                                                 *out << Verbose(2) << "Keeping new vector." << endl;
    1466                                         } else {
    1467                                                 *out << Verbose(2) << "Keeping present vector." << endl;
    1468                                         }
    1469                                 } else {
    1470                                                 *out << Verbose(2) << "Keeping present vector." << endl;
    1471                                 }
    1472                         }
    1473                 }
    1474                 // printing all inserted for debugging
    1475                 {
    1476       *out << Verbose(2) << "Printing list of candidates for axis " << axis << " which we have inserted so far." << endl;
    1477       int i=0;
    1478       for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    1479         if (runner != BoundaryPoints[axis].begin())
    1480           *out << ", " << i << ": " << *runner->second.second;
    1481         else
    1482           *out << i << ": " << *runner->second.second;
    1483         i++;
    1484       }
    1485       *out << endl;
    1486                 }
    1487     // 3c. throw out points whose distance is less than the mean of left and right neighbours
    1488                 bool flag = false;
    1489                 do { // do as long as we still throw one out per round
    1490                   *out << Verbose(1) << "Looking for candidates to kick out by convex condition ... " << endl;
    1491                         flag = false;
    1492                         Boundaries::iterator left = BoundaryPoints[axis].end();
    1493                         Boundaries::iterator right = BoundaryPoints[axis].end();
    1494                         for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    1495                                 // set neighbours correctly
    1496                                 if (runner == BoundaryPoints[axis].begin()) {
    1497                                         left = BoundaryPoints[axis].end();
    1498                                 }       else {
    1499                                         left = runner;
    1500                                 }
    1501         left--;
    1502         right = runner;
    1503         right++;
    1504                                 if (right == BoundaryPoints[axis].end()) {
    1505                                         right = BoundaryPoints[axis].begin();
    1506                                 }
    1507                                 // check distance
    1508                                
    1509                                 // construct the vector of each side of the triangle on the projected plane (defined by normal vector AxisVector)
    1510                                 {
    1511           vector SideA, SideB, SideC, SideH;
    1512           SideA.CopyVector(&left->second.second->x);
    1513           SideA.ProjectOntoPlane(&AxisVector);
    1514 //          *out << "SideA: ";
    1515 //          SideA.Output(out);
    1516 //          *out << endl;
    1517          
    1518           SideB.CopyVector(&right->second.second->x);
    1519           SideB.ProjectOntoPlane(&AxisVector);
    1520 //          *out << "SideB: ";
    1521 //          SideB.Output(out);
    1522 //          *out << endl;
    1523          
    1524           SideC.CopyVector(&left->second.second->x);
    1525           SideC.SubtractVector(&right->second.second->x);
    1526           SideC.ProjectOntoPlane(&AxisVector);
    1527 //          *out << "SideC: ";
    1528 //          SideC.Output(out);
    1529 //          *out << endl;
    1530 
    1531           SideH.CopyVector(&runner->second.second->x);
    1532           SideH.ProjectOntoPlane(&AxisVector);
    1533 //          *out << "SideH: ";
    1534 //          SideH.Output(out);
    1535 //          *out << endl;
    1536          
    1537           // calculate each length
    1538           double a = SideA.Norm();
    1539           //double b = SideB.Norm();
    1540           //double c = SideC.Norm();
    1541           double h = SideH.Norm();
    1542           // calculate the angles
    1543           double alpha = SideA.Angle(&SideH);
    1544                                 double beta = SideA.Angle(&SideC);
    1545                                 double gamma = SideB.Angle(&SideH);
    1546                                 double delta = SideC.Angle(&SideH);
    1547                                 double MinDistance = a * sin(beta)/(sin(delta)) * (((alpha < M_PI/2.) || (gamma < M_PI/2.)) ? 1. : -1.);
    1548 //                              *out << Verbose(2) << " I calculated: a = " << a << ", h = " << h << ", beta(" << left->second.second->Name << "," << left->second.second->Name << "-" << right->second.second->Name << ") = " << beta << ", delta(" << left->second.second->Name << "," << runner->second.second->Name << ") = " << delta << ", Min = " << MinDistance << "." << endl;
    1549                                 *out << Verbose(1) << "Checking CoG distance of runner " << *runner->second.second << " " << h << " against triangle's side length spanned by (" << *left->second.second << "," << *right->second.second << ") of " << MinDistance << "." << endl;
    1550                                 if ((fabs(h/fabs(h) - MinDistance/fabs(MinDistance)) < MYEPSILON) && (h <  MinDistance)) {
    1551                                         // throw out point
    1552                                   *out << Verbose(1) << "Throwing out " << *runner->second.second << "." << endl;
    1553                                         BoundaryPoints[axis].erase(runner);
    1554                                         flag = true;
    1555                                 }
    1556                                 }
    1557                         }
    1558                 } while (flag);
    1559         }       
    1560         // 3d. put into boundary set only those points appearing in each of the NDIM sets
    1561         int *AtomList = new int[AtomCount];
    1562         for (int j=0; j<AtomCount; j++) // reset list
    1563                 AtomList[j] = 0;
    1564         for (int axis=0; axis<NDIM; axis++)  {  // increase whether it's present
    1565                 for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    1566                         AtomList[runner->second.second->nr]++;
    1567                 }
    1568         }
    1569 
    1570         // 3e. Points no more in the total list, have to be thrown out of each axis lists too!
    1571         for (int axis=0; axis<NDIM; axis++)  {
    1572                 for(Boundaries::iterator runner = BoundaryPoints[axis].begin(); runner != BoundaryPoints[axis].end(); runner++) {
    1573                         if (AtomList[runner->second.second->nr] < 1) {
    1574                           *out << Verbose(1) << "Throwing out " << *runner->second.second << " in axial projection of axis " << axis << "." << endl;
    1575                                 BoundaryPoints[axis].erase(runner);
    1576                         }
    1577                 }
    1578         }
    1579        
    1580         // listing boundary points
    1581         *out << Verbose(1) << "The following atoms are on the boundary:";
    1582   int BoundaryPointCount = 0;
    1583   Walker = start;
    1584   while (Walker->next != end) {
    1585     Walker = Walker->next;
    1586     if (AtomList[Walker->nr] > 0) {
    1587       BoundaryPointCount ++;
    1588       *out << " " << Walker->Name;
    1589     }
    1590   }
    1591   *out << endl;
    1592 
    1593         *out << Verbose(2) << "I found " << BoundaryPointCount << " points on the convex boundary." << endl;
    1594         // now we have the whole set of edge points in the BoundaryList
    1595        
    1596         // 4. Create each possible line, number them and put them into a list (3 * AtomCount possible)
    1597         LinesMap AllTriangleLines;
    1598        
    1599         // then create each line
    1600   *out << Verbose(1) << "Creating lines between adjacent boundary points." << endl;
    1601   atom *Runner = NULL;
    1602   int LineCount = 0;
    1603   Walker = start;
    1604   while (Walker->next != end) {
    1605     Walker = Walker->next;
    1606     if (AtomList[Walker->nr] > 0) { // boundary point!
    1607       *out << Verbose(1) << "Current Walker is " << *Walker << "." << endl;
    1608       // make a list of all neighbours
    1609       DistanceMap Distances;
    1610       double distance = 0;
    1611       Runner = start;
    1612       while (Runner->next != end) {
    1613         Runner = Runner->next;
    1614         if ((AtomList[Runner->nr] > 0) && (Runner != Walker)) { // don't mess with ourselves!
    1615           distance = Walker->x.Distance(&Runner->x);    // note this is squared distance
    1616           Distances.insert ( DistanceNrPair( distance, Runner) );
    1617           *out << Verbose(2) << "Inserted distance " << distance << " to " << *Runner << " successfully." << endl;
    1618         }
    1619       } // end of distance filling while loop
    1620       // take 3NNs to draw a line in between
    1621       int Counter = 0;  // counts up to three inserted lines
    1622       LinesTestPair LinesTest;
    1623       for (DistanceMap::iterator spanner = Distances.begin(); spanner != Distances.end(); spanner++) {
    1624         *out << Verbose(1) << "Current distance to insert is " << spanner->first << " to " << spanner->second+1 << "." << endl;
    1625         LinesTest = AllTriangleLines.insert ( LinesPair( (Walker->nr*AtomCount)+spanner->second->nr, LineCount) );
    1626         if (!LinesTest.second) {
    1627           *out << Verbose(2) << "Insertion of line between " << Walker->nr+1 << " and " << spanner->second->nr+1 << " failed, already present line nr. " << LinesTest.first->second << "." << endl;
    1628         } else {
    1629           LineCount++;
    1630           Counter++;
    1631           *out << Verbose(2) << "Insertion of line between " << (LinesTest.first->first/AtomCount)+1 << " and " << (LinesTest.first->first%AtomCount)+1 << " successful, inserting mirrored line and increased counter to " << Counter << "." << endl;
    1632           LinesTest = AllTriangleLines.insert ( LinesPair( Walker->nr+(spanner->second->nr*AtomCount), LineCount) );
    1633           if (!LinesTest.second) {
    1634             *out << Verbose(2) << "Insertion of mirror line between " << Walker->nr << " and " << spanner->second->nr << " failed, already present line nr. " << LinesTest.first->second << "." << endl;
    1635             *out << Verbose(0) << "SERIOUS ERROR!!!" << endl;
    1636           }
    1637         }
    1638         if (Counter >= 3) { // stop after three lines (if not run out of possible NNs already)
    1639           *out << Verbose(2) << "Ok, three lines inserted from current walker " << *Walker << ", stopping." << endl;
    1640           break;
    1641         }
    1642       }
    1643     } // end of if boundary pointy
    1644   } // end of loop through all atoms
    1645 
    1646         // listing created lines
    1647         *out << Verbose(2) << "The following lines were created:";
    1648         for (LinesMap::iterator liner = AllTriangleLines.begin(); liner != AllTriangleLines.end(); liner++)
    1649           if ((liner->first/AtomCount) < (liner->first%AtomCount))
    1650             *out << " " << (liner->first/AtomCount)+1 << "<->" << (liner->first%AtomCount)+1;
    1651         *out << endl;
    1652         *out << Verbose(2) << "I created " << LineCount << " edges." << endl;
    1653        
    1654         // 5. Create every possible triangle from the lines (numbering of each line must be i < j < k)
    1655         struct triangles {
    1656                 atom *x[3];
    1657                 int nr;
    1658         } TriangleList[2 * 3 * AtomCount];      // each line can be part in at most 2 triangles
    1659 
    1660   *out << Verbose(1) << "Creating triangles out of these lines." << endl;
    1661         int TriangleCount = 0;
    1662   LinetoAtomPair InsertionPair = LinetoAtomPair(-1, NULL);
    1663   int endpoints[3];
    1664   for (LinesMap::iterator liner = AllTriangleLines.begin(); liner != AllTriangleLines.end(); liner++) { // go through every line
    1665                 LinetoAtomMap LinetoAtom;
    1666                 // we have two points, check the lines at either end, whether they share a common endpoint
    1667     *out << Verbose(1) << "Examining line between " << (liner->first/AtomCount)+1 << " and " << (liner->first%AtomCount)+1 << "." << endl;
    1668     int enden[3][2];
    1669     enden[0][0] = (liner->first/AtomCount);
    1670     enden[0][1] = (liner->first%AtomCount);
    1671     if (enden[0][0] < enden[0][1]) {
    1672                 for (int endpoint=0;endpoint<2;endpoint++) {
    1673                         endpoints[0] = enden[0][endpoint];
    1674                         *out << Verbose(2) << "Current first endpoint is " << endpoints[0]+1 << "." << endl;
    1675                        
    1676                         LinesMap::iterator LineRangeStart = AllTriangleLines.lower_bound(endpoints[0]*AtomCount);
    1677                         for (LinesMap::iterator coach = LineRangeStart; (coach->first/AtomCount) == endpoints[0]; coach++) { // look at all line's other endpoint
    1678                     enden[1][0] = (coach->first/AtomCount);
    1679                     enden[1][1] = (coach->first%AtomCount);
    1680           //*out << Verbose(3) << "Line is " << coach->first << " with nr. " << coach->second << ": " << enden[1][0+1] << "<->" << enden[1][1]+1 << "." << endl;
    1681                           endpoints[1] = (enden[1][0] != endpoints[0] ) ? enden[1][0] : enden[1][1];
    1682                           if ((endpoints[1] > endpoints[0]) && (endpoints[1] != enden[0][(endpoint+1)%2])) {  // don't go back the very same line!
    1683                           *out << Verbose(3) << "Current second endpoint is " << endpoints[1]+1 << "." << endl; 
    1684                          
    1685               LinesMap::iterator LineRangeStart2 = AllTriangleLines.lower_bound(endpoints[1]*AtomCount);
    1686               for (LinesMap::iterator coacher = LineRangeStart2; (coacher->first/AtomCount) == endpoints[1]; coacher++) { // look at all line's other endpoint
    1687                 enden[2][0] = (coacher->first/AtomCount);
    1688                 enden[2][1] = (coacher->first%AtomCount);
    1689               //*out << Verbose(3) << "Line is " << coacher->first << " with nr. " << coacher->second << ": " << enden[2][0]+1 << "<->" << enden[2][1]+1 << "." << endl;
    1690                 endpoints[2] = (enden[2][0] != endpoints[1]) ? enden[2][0] : enden[2][1];
    1691                 if ((endpoints[2] > endpoints[1]) && (endpoints[2] != endpoints[0])) {  // don't go back the very same line!
    1692                 *out << Verbose(4) << "Current third endpoint is " << endpoints[2]+1 << "." << endl;
    1693                
    1694                 if (endpoints[2] == enden[0][(endpoint+1)%2]) { // jo, closing triangle
    1695                   *out << Verbose(3) << "MATCH: Triangle is ";
    1696                   for (int k=0;k<3;k++) {
    1697                     *out << endpoints[k]+1 << "\t";
    1698                     TriangleList[TriangleCount].x[k] = FindAtom(endpoints[k]);
    1699                   }
    1700                   *out << endl;
    1701                   TriangleList[TriangleCount].nr = TriangleCount;
    1702                   TriangleCount++;
    1703                 } else {
    1704                   *out << Verbose(3) << "No match!" << endl;
    1705                 }
    1706               }
    1707               }
    1708                         }
    1709                         }
    1710                 }
    1711     }
    1712         }
    1713   // listing created lines
    1714   *out << Verbose(2) << "The following triangles were created:";
    1715   for (int i=0;i<TriangleCount;i++) {
    1716     *out << " " << TriangleList[i].x[0]->Name << "<->" << TriangleList[i].x[1]->Name << "<->" << TriangleList[i].x[2]->Name;
    1717   }
    1718   *out << endl;
    1719         *out << Verbose(2) << "I created " << TriangleCount << " triangles." << endl;
    1720 
    1721         // 6. Every triangle forms a pyramid with the center of gravity as its peak, sum up the volumes
    1722   *out << Verbose(1) << "Calculating the volume of the pyramids formed out of triangles and center of gravity." << endl;
    1723         double volume = 0.;
    1724         double PyramidVolume = 0.;
    1725         double G,h;
    1726         vector x,y;
    1727         double a,b,c;
    1728         for (int i=0; i<TriangleCount; i++) { // go through every triangle, calculate volume of its pyramid with CoG as peak
    1729                 x.CopyVector(&TriangleList[i].x[0]->x);
    1730                 x.SubtractVector(&TriangleList[i].x[1]->x);
    1731                 y.CopyVector(&TriangleList[i].x[0]->x);
    1732                 y.SubtractVector(&TriangleList[i].x[2]->x);
    1733                 a = sqrt(TriangleList[i].x[0]->x.Distance(&TriangleList[i].x[1]->x));
    1734                 b = sqrt(TriangleList[i].x[0]->x.Distance(&TriangleList[i].x[2]->x));
    1735                 c = sqrt(TriangleList[i].x[2]->x.Distance(&TriangleList[i].x[1]->x));
    1736                 G =  sqrt( ( (a*a+b*b+c*c)*(a*a+b*b+c*c) - 2*(a*a*a*a + b*b*b*b + c*c*c*c) )/16.); // area of tesselated triangle
    1737                 x.MakeNormalVector(&TriangleList[i].x[0]->x, &TriangleList[i].x[1]->x, &TriangleList[i].x[2]->x);
    1738                 x.Scale(TriangleList[i].x[1]->x.Projection(&x));
    1739                 h = x.Norm(); // distance of CoG to triangle
    1740                 PyramidVolume = (1./3.) * G * h;                // this formula holds for _all_ pyramids (independent of n-edge base or (not) centered peak)
    1741                 *out << Verbose(2) << "Area of triangle is " << G << "A^2, height is " << h << " and the volume is " << PyramidVolume << "." << endl;
    1742                 volume += PyramidVolume;
    1743         }
    1744         *out << Verbose(1) << "The summed volume is " << volume << " " << (IsAngstroem ? "A^3" : "a_0^3") << "." << endl;
    1745 
    1746   // 7. translate all points back from CoG
    1747   *out << Verbose(1) << "Translating system back from Center of Gravity." << endl;
    1748         CenterOfGravity->Scale(-1);
    1749   Walker=start;
    1750   while (Walker->next != end) {
    1751     Walker = Walker->next;
    1752     Walker->x.Translate(CenterOfGravity);
    1753   }
    1754 
    1755   // free reference lists
    1756         delete[](AtomList);
    1757 
    1758         return volume;
    1759 };
    1760 
    17611372/** Returns Shading as a char string.
    17621373 * \param color the Shading
Note: See TracChangeset for help on using the changeset viewer.