Changeset 3839


Ignore:
Timestamp:
Jul 25, 2013 2:07:37 PM (7 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:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/packages/vizservers/vtkvis/Molecule.cpp

    r3795 r3839  
    183183    addLabelArray(ds);
    184184    addRadiusArray(ds, _atomScaling, _radiusScale);
     185    computeBonds(ds);
    185186
    186187    vtkPolyData *pd = vtkPolyData::SafeDownCast(ds);
     
    10591060    return elementCmap;
    10601061}
     1062
     1063int Molecule::computeBonds(vtkDataSet *dataSet, float tolerance)
     1064{
     1065    vtkPolyData *pd = vtkPolyData::SafeDownCast(dataSet);
     1066    if (pd == NULL) {
     1067        ERROR("DataSet not a PolyData");
     1068        return -1;
     1069    }
     1070
     1071    if (pd->GetNumberOfPoints() < 2 ||
     1072        pd->GetNumberOfLines() > 0) {
     1073        return 0;
     1074    }
     1075
     1076    vtkDataArray *elements = NULL;
     1077    if (pd->GetPointData() != NULL &&
     1078        pd->GetPointData()->GetScalars() != NULL &&
     1079        strcmp(pd->GetPointData()->GetScalars()->GetName(), "element") == 0) {
     1080        elements = pd->GetPointData()->GetScalars();
     1081    } else {
     1082        TRACE("Can't compute bonds without an element array");
     1083        return -1;
     1084    }
     1085    int numAtoms = pd->GetNumberOfPoints();
     1086    vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
     1087    for (int i = 0; i < numAtoms; i++) {
     1088        double pt1[3];
     1089        int elem1 = (int)elements->GetComponent(i, 0);
     1090        float r1 = g_covalentRadii[elem1];
     1091        pd->GetPoint(i, pt1);
     1092        for (int j = i+1; j < numAtoms; j++) {
     1093            float r2, dist;
     1094            double pt2[3];
     1095            int elem2 = (int)elements->GetComponent(j, 0);
     1096
     1097            if (elem1 == 1 && elem2 == 1)
     1098                continue;
     1099
     1100            r2 = g_covalentRadii[elem2];
     1101            pd->GetPoint(j, pt2);
     1102
     1103            double x = pt1[0] - pt2[0];
     1104            double y = pt1[1] - pt2[1];
     1105            double z = pt1[2] - pt2[2];
     1106            dist = (float)sqrt(x*x + y*y + z*z);
     1107            //TRACE("i=%d,j=%d elem1=%d,elem2=%d r1=%g,r2=%g, dist=%g", i, j, elem1, elem2, r1, r2, dist);
     1108            if (dist < (r1 + r2 + tolerance)) {
     1109                vtkIdType pts[2];
     1110                pts[0] = i;
     1111                pts[1] = j;
     1112                cells->InsertNextCell(2, pts);
     1113            }
     1114        }
     1115    }
     1116
     1117    if (cells->GetNumberOfCells() > 0) {
     1118        pd->SetLines(cells);
     1119    }
     1120
     1121    TRACE("Generated %d bonds", cells->GetNumberOfCells());
     1122    return cells->GetNumberOfCells();
     1123}
  • trunk/packages/vizservers/vtkvis/Molecule.h

    r3795 r3839  
    190190    static void addRadiusArray(vtkDataSet *dataSet, AtomScaling scaling, double scaleFactor);
    191191
     192    static int computeBonds(vtkDataSet *dataSet, float tolerance = 0.45f);
     193
    192194    float _bondColor[3];
    193195    double _radiusScale;
Note: See TracChangeset for help on using the changeset viewer.