1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkMoleculeToBondStickFilter.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 #include "vtkMoleculeToBondStickFilter.h"
16
17 #include "vtkCellArray.h"
18 #include "vtkCylinderSource.h"
19 #include "vtkObjectFactory.h"
20 #include "vtkMath.h"
21 #include "vtkMolecule.h"
22 #include "vtkPointData.h"
23 #include "vtkPoints.h"
24 #include "vtkPolyData.h"
25 #include "vtkTransform.h"
26 #include "vtkUnsignedShortArray.h"
27
28 vtkStandardNewMacro(vtkMoleculeToBondStickFilter);
29
30 //----------------------------------------------------------------------------
vtkMoleculeToBondStickFilter()31 vtkMoleculeToBondStickFilter::vtkMoleculeToBondStickFilter()
32 {
33 }
34
35 //----------------------------------------------------------------------------
~vtkMoleculeToBondStickFilter()36 vtkMoleculeToBondStickFilter::~vtkMoleculeToBondStickFilter()
37 {
38 }
39
40 //----------------------------------------------------------------------------
RequestData(vtkInformation *,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)41 int vtkMoleculeToBondStickFilter::RequestData(
42 vtkInformation *,
43 vtkInformationVector **inputVector,
44 vtkInformationVector *outputVector)
45 {
46 vtkMolecule *input = vtkMolecule::SafeDownCast
47 (vtkDataObject::GetData(inputVector[0]));
48 vtkPolyData *output = vtkPolyData::SafeDownCast
49 (vtkDataObject::GetData(outputVector));
50
51 // Extract data from input molecule
52 vtkIdType numBonds = input->GetNumberOfBonds();
53
54 // Prep the output
55 output->Initialize();
56 vtkCellArray *polys = vtkCellArray::New();
57 vtkPoints *points = vtkPoints::New();
58 vtkUnsignedShortArray *bondOrders = vtkUnsignedShortArray::New();
59
60 // Initialize a CylinderSource
61 vtkCylinderSource *cylSource = vtkCylinderSource::New();
62 cylSource->SetResolution(20);
63 cylSource->SetHeight(1.0);
64 cylSource->Update();
65
66 // Preallocate memory
67 points->Allocate(3 * numBonds * cylSource->GetOutput()->GetPoints()->
68 GetNumberOfPoints());
69 polys->Allocate(3 * numBonds * cylSource->GetOutput()->GetPolys()->
70 GetNumberOfCells());
71 bondOrders->Allocate(points->GetNumberOfPoints());
72
73 // Create a transform object to map the cylinder source to the bond
74 vtkTransform *xform = vtkTransform::New();
75 xform->PostMultiply();
76
77 // Declare some variables for later
78 vtkIdType numCellPoints, *cellPoints;
79 unsigned short bondOrder;
80 double bondLength;
81 double radius;
82 double delta[3];
83 double initialDisp[3];
84 double rotAngle;
85 double rotAxis[3];
86 double bondCenter[3];
87 double pos1[3], pos2[3];
88 // Normalized vector pointing along the cylinder (y axis);
89 const static double cylVec[3] = {0.0, 1.0, 0.0};
90 // Normalized vector pointing along bond
91 double bondVec[3];
92 // Unit z vector
93 const static double unitZ[3] = {0.0, 0.0, 1.0};
94
95 // Build a sphere for each atom and append it's data to the output
96 // arrays.
97 for (vtkIdType bondInd = 0; bondInd < numBonds; ++bondInd)
98 {
99 // Extract bond info
100 vtkBond bond = input->GetBond(bondInd);
101 bondOrder = bond.GetOrder();
102 bond.GetBeginAtom().GetPosition(pos1);
103 bond.GetEndAtom().GetPosition(pos2);
104
105 // Compute additional bond info
106 // - Normalized vector in direction of bond
107 vtkMath::Subtract(pos2, pos1, bondVec);
108 bondLength = vtkMath::Normalize(bondVec);
109 // - Axis for cylinder rotation, bondVec [cross] cylVec
110 vtkMath::Cross(bondVec, cylVec, rotAxis);
111 // Rotation angle
112 rotAngle = -vtkMath::DegreesFromRadians
113 (acos(vtkMath::Dot(bondVec, cylVec)));
114 // - Center of bond for translation
115 vtkMath::Add(pos2, pos1, bondCenter);
116 vtkMath::MultiplyScalar(bondCenter, 0.5);
117
118 // Set up delta step vector and bond radius from bond order:
119 switch (bondOrder)
120 {
121 case 1:
122 default:
123 radius = 0.1;
124 delta[0] = 0.0; delta[1] = 0.0; delta[2] = 0.0;
125 initialDisp[0] = 0.0; initialDisp[1] = 0.0; initialDisp[2] = 0.0;
126 break;
127 case 2:
128 radius = 0.1;
129 vtkMath::Cross(bondVec, unitZ, delta);
130 //vtkMath::Normalize(delta);
131 vtkMath::MultiplyScalar
132 (delta, radius + radius);
133 initialDisp[0] = delta[0] * -0.5;
134 initialDisp[1] = delta[1] * -0.5;
135 initialDisp[2] = delta[2] * -0.5;
136 break;
137 case 3:
138 radius = 0.1;
139 vtkMath::Cross(bondVec, unitZ, delta);
140 //vtkMath::Normalize(delta);
141 vtkMath::MultiplyScalar
142 (delta, radius + radius);
143 initialDisp[0] = -delta[0];
144 initialDisp[1] = -delta[1];
145 initialDisp[2] = -delta[2];
146 break;
147 }
148
149 // Construct transform
150 xform->Identity();
151 xform->Scale(radius, bondLength, radius);
152 xform->RotateWXYZ(rotAngle, rotAxis);
153 xform->Translate(bondCenter[0], bondCenter[1], bondCenter[2]);
154 xform->Translate(initialDisp);
155
156 // For each bond order, add a cylinder to output, translate by
157 // delta, and repeat.
158 for (unsigned short iter = 0; iter < bondOrder; ++iter)
159 {
160 vtkPolyData *cylinder = cylSource->GetOutput();
161 vtkPoints *cylPoints = cylinder->GetPoints();
162 vtkCellArray *cylPolys = cylinder->GetPolys();
163 xform->TransformPoints(cylPoints, points);
164
165 // Get offset for the new point IDs that will be added to points
166 vtkIdType pointOffset = points->GetNumberOfPoints();
167 // Total number of new points
168 vtkIdType numPoints = cylPoints->GetNumberOfPoints();
169
170 // Use bond order for point scalar data.
171 for (vtkIdType i = 0; i < numPoints; ++i)
172 {
173 bondOrders->InsertNextValue(bondOrder);
174 }
175
176 // Add new cells (polygons) that represent the cylinder
177 cylPolys->InitTraversal();
178 while (cylPolys->GetNextCell(numCellPoints, cellPoints) != 0)
179 {
180 vtkIdType *newCellPoints = new vtkIdType[numCellPoints];
181 for (vtkIdType i = 0; i < numCellPoints; ++i)
182 {
183 // The new point ids should be offset by the pointOffset above
184 newCellPoints[i] = cellPoints[i] + pointOffset;
185 }
186 polys->InsertNextCell(numCellPoints, newCellPoints);
187 delete [] newCellPoints;
188 }
189
190 // Setup for the next cylinder in a multi-bond
191 xform->Translate(delta);
192 }
193 }
194
195 // Release extra memory
196 points->Squeeze();
197 bondOrders->Squeeze();
198 polys->Squeeze();
199
200 // update output
201 output->SetPoints(points);
202 output->GetPointData()->SetScalars(bondOrders);
203 output->SetPolys(polys);
204
205 // Clean up
206 xform->Delete();
207 polys->Delete();
208 points->Delete();
209 bondOrders->Delete();
210 cylSource->Delete();
211
212 return 1;
213 }
214
215 //----------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)216 void vtkMoleculeToBondStickFilter::PrintSelf(ostream& os, vtkIndent indent)
217 {
218 this->Superclass::PrintSelf(os,indent);
219 }
220