Changeset 3839

Timestamp:

Jul 25, 2013 2:07:37 PM (10 years ago)

Author:

ldelgass

Message:

Add automatic bond generation to VTK molecules: if no bonds are specified (as
lines in VTK file), computes simple single bonds based only on covalent radii.
DataSets? must include the element array for bonds to be computed.

Location:

trunk/packages/vizservers/vtkvis

Files:

• Molecule.cpp (modified) (2 diffs)
• Molecule.h (modified) (1 diff)

trunk/packages/vizservers/vtkvis/Molecule.cpp

@@ -183 +183 @@
     addLabelArray(ds);
     addRadiusArray(ds, _atomScaling, _radiusScale);
+    computeBonds(ds);
 
     vtkPolyData *pd = vtkPolyData::SafeDownCast(ds);
@@ -1059 +1060 @@
     return elementCmap;
 }
+
+int Molecule::computeBonds(vtkDataSet *dataSet, float tolerance)
+{
+    vtkPolyData *pd = vtkPolyData::SafeDownCast(dataSet);
+    if (pd == NULL) {
+        ERROR("DataSet not a PolyData");
+        return -1;
+    }
+
+    if (pd->GetNumberOfPoints() < 2 ||
+        pd->GetNumberOfLines() > 0) {
+        return 0;
+    }
+
+    vtkDataArray *elements = NULL;
+    if (pd->GetPointData() != NULL &&
+        pd->GetPointData()->GetScalars() != NULL &&
+        strcmp(pd->GetPointData()->GetScalars()->GetName(), "element") == 0) {
+        elements = pd->GetPointData()->GetScalars();
+    } else {
+        TRACE("Can't compute bonds without an element array");
+        return -1;
+    }
+    int numAtoms = pd->GetNumberOfPoints();
+    vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
+    for (int i = 0; i < numAtoms; i++) {
+        double pt1[3];
+        int elem1 = (int)elements->GetComponent(i, 0);
+        float r1 = g_covalentRadii[elem1];
+        pd->GetPoint(i, pt1);
+        for (int j = i+1; j < numAtoms; j++) {
+            float r2, dist;
+            double pt2[3];
+            int elem2 = (int)elements->GetComponent(j, 0);
+
+            if (elem1 == 1 && elem2 == 1)
+                continue;
+
+            r2 = g_covalentRadii[elem2];
+            pd->GetPoint(j, pt2);
+
+            double x = pt1[0] - pt2[0];
+            double y = pt1[1] - pt2[1];
+            double z = pt1[2] - pt2[2];
+            dist = (float)sqrt(x*x + y*y + z*z);
+            //TRACE("i=%d,j=%d elem1=%d,elem2=%d r1=%g,r2=%g, dist=%g", i, j, elem1, elem2, r1, r2, dist); 
+            if (dist < (r1 + r2 + tolerance)) {
+                vtkIdType pts[2];
+                pts[0] = i;
+                pts[1] = j;
+                cells->InsertNextCell(2, pts);
+            }
+        }
+    }
+
+    if (cells->GetNumberOfCells() > 0) {
+        pd->SetLines(cells);
+    }
+
+    TRACE("Generated %d bonds", cells->GetNumberOfCells());
+    return cells->GetNumberOfCells();
+}

trunk/packages/vizservers/vtkvis/Molecule.h

@@ -190 +190 @@
     static void addRadiusArray(vtkDataSet *dataSet, AtomScaling scaling, double scaleFactor);
 
+    static int computeBonds(vtkDataSet *dataSet, float tolerance = 0.45f);
+
     float _bondColor[3];
     double _radiusScale;