1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkMoleculeAppend.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 "vtkMoleculeAppend.h"
16
17 #include "vtkAlgorithmOutput.h"
18 #include "vtkDataArray.h"
19 #include "vtkDataSetAttributes.h"
20 #include "vtkExecutive.h"
21 #include "vtkInformation.h"
22 #include "vtkInformationVector.h"
23 #include "vtkMergePoints.h"
24 #include "vtkMolecule.h"
25 #include "vtkNew.h"
26 #include "vtkObjectFactory.h"
27 #include "vtkPoints.h"
28 #include "vtkUnsignedCharArray.h"
29
30 #include <set>
31 #include <utility>
32
33 vtkStandardNewMacro(vtkMoleculeAppend);
34
35 //------------------------------------------------------------------------------
PrintSelf(ostream & os,vtkIndent indent)36 void vtkMoleculeAppend::PrintSelf(ostream& os, vtkIndent indent)
37 {
38 this->Superclass::PrintSelf(os, indent);
39 os << indent << "MergeCoincidentAtoms: " << this->MergeCoincidentAtoms << endl;
40 }
41
42 //------------------------------------------------------------------------------
vtkMoleculeAppend()43 vtkMoleculeAppend::vtkMoleculeAppend()
44 : MergeCoincidentAtoms(true)
45 {
46 }
47
48 //------------------------------------------------------------------------------
GetInput(int idx)49 vtkDataObject* vtkMoleculeAppend::GetInput(int idx)
50 {
51 if (this->GetNumberOfInputConnections(0) <= idx)
52 {
53 return nullptr;
54 }
55 return vtkMolecule::SafeDownCast(this->GetExecutive()->GetInputData(0, idx));
56 }
57
58 //------------------------------------------------------------------------------
RequestData(vtkInformation *,vtkInformationVector ** inputVector,vtkInformationVector * outputVector)59 int vtkMoleculeAppend::RequestData(
60 vtkInformation*, vtkInformationVector** inputVector, vtkInformationVector* outputVector)
61 {
62 vtkMolecule* output = vtkMolecule::GetData(outputVector, 0);
63 vtkDataSetAttributes* outputAtomData = output->GetAtomData();
64 vtkDataSetAttributes* outputBondData = output->GetBondData();
65
66 // ********************
67 // Create output data arrays following first input arrays.
68 vtkMolecule* mol0 = vtkMolecule::SafeDownCast(this->GetInput());
69 outputAtomData->CopyStructure(mol0->GetAtomData());
70 outputBondData->CopyStructure(mol0->GetBondData());
71 output->SetAtomicNumberArrayName(mol0->GetAtomicNumberArrayName());
72 output->SetBondOrdersArrayName(mol0->GetBondOrdersArrayName());
73 vtkUnsignedCharArray* outputGhostAtoms = output->GetAtomGhostArray();
74 vtkUnsignedCharArray* outputGhostBonds = output->GetBondGhostArray();
75
76 // ********************
77 // Initialize unique atoms/bonds containers
78 vtkNew<vtkMergePoints> uniquePoints;
79 vtkNew<vtkPoints> uniquePointsList;
80 double bounds[6] = { 0., 0., 0., 0., 0., 0. };
81 uniquePoints->InitPointInsertion(uniquePointsList, bounds, 0);
82 std::set<std::pair<vtkIdType, vtkIdType>> uniqueBonds;
83
84 // ********************
85 // Process each input
86 for (int idx = 0; idx < this->GetNumberOfInputConnections(0); ++idx)
87 {
88 vtkMolecule* input = vtkMolecule::GetData(inputVector[0], idx);
89
90 // --------------------
91 // Sanity check on input
92 int inputNbAtomArrays = input->GetAtomData()->GetNumberOfArrays();
93 if (inputNbAtomArrays != outputAtomData->GetNumberOfArrays())
94 {
95 vtkErrorMacro(<< "Input " << idx << ": Wrong number of atom array. Has " << inputNbAtomArrays
96 << " instead of " << outputAtomData->GetNumberOfArrays());
97 return 0;
98 }
99
100 int inputNbBondArrays = input->GetBondData()->GetNumberOfArrays();
101 if (input->GetNumberOfBonds() > 0 && inputNbBondArrays != outputBondData->GetNumberOfArrays())
102 {
103 vtkErrorMacro(<< "Input " << idx << ": Wrong number of bond array. Has " << inputNbBondArrays
104 << " instead of " << outputBondData->GetNumberOfArrays());
105 return 0;
106 }
107
108 for (vtkIdType ai = 0; ai < inputNbAtomArrays; ai++)
109 {
110 vtkAbstractArray* inArray = input->GetAtomData()->GetAbstractArray(ai);
111 if (!this->CheckArrays(inArray, outputAtomData->GetAbstractArray(inArray->GetName())))
112 {
113 vtkErrorMacro(<< "Input " << idx << ": atoms arrays do not match with output");
114 return 0;
115 }
116 }
117
118 for (vtkIdType ai = 0; ai < inputNbBondArrays; ai++)
119 {
120 vtkAbstractArray* inArray = input->GetBondData()->GetAbstractArray(ai);
121 if (!this->CheckArrays(inArray, outputBondData->GetAbstractArray(inArray->GetName())))
122 {
123 vtkErrorMacro(<< "Input " << idx << ": bonds arrays do not match with output");
124 return 0;
125 }
126 }
127
128 // --------------------
129 // add atoms and bonds without duplication
130
131 // map from 'input molecule atom ids' to 'output molecule atom ids'
132 std::vector<vtkIdType> atomIdMap(input->GetNumberOfAtoms(), -1);
133
134 vtkIdType previousNbOfAtoms = output->GetNumberOfAtoms();
135 int nbOfAtoms = 0;
136 for (vtkIdType i = 0; i < input->GetNumberOfAtoms(); i++)
137 {
138 double pt[3];
139 input->GetAtomicPositionArray()->GetPoint(i, pt);
140 bool addAtom = true;
141 if (this->MergeCoincidentAtoms)
142 {
143 addAtom = uniquePoints->InsertUniquePoint(pt, atomIdMap[i]) == 1;
144 }
145 else
146 {
147 atomIdMap[i] = previousNbOfAtoms + nbOfAtoms;
148 }
149
150 if (addAtom)
151 {
152 nbOfAtoms++;
153 vtkAtom atom = input->GetAtom(i);
154 output->AppendAtom(atom.GetAtomicNumber(), atom.GetPosition()).GetId();
155 if (outputGhostAtoms)
156 {
157 outputGhostAtoms->InsertValue(atomIdMap[i], 255);
158 }
159 }
160 }
161 vtkIdType previousNbOfBonds = output->GetNumberOfBonds();
162 int nbOfBonds = 0;
163 for (vtkIdType i = 0; i < input->GetNumberOfBonds(); i++)
164 {
165 vtkBond bond = input->GetBond(i);
166 // as bonds are undirected, put min atom number at first to avoid duplication.
167 vtkIdType atom1 = atomIdMap[bond.GetBeginAtomId()];
168 vtkIdType atom2 = atomIdMap[bond.GetEndAtomId()];
169 auto result =
170 uniqueBonds.insert(std::make_pair(std::min(atom1, atom2), std::max(atom1, atom2)));
171 if (result.second)
172 {
173 nbOfBonds++;
174 output->AppendBond(atom1, atom2, bond.GetOrder());
175 }
176 }
177
178 // --------------------
179 // Reset arrays size (and allocation if needed)
180 for (vtkIdType ai = 0; ai < input->GetAtomData()->GetNumberOfArrays(); ai++)
181 {
182 vtkAbstractArray* inArray = input->GetAtomData()->GetAbstractArray(ai);
183 vtkAbstractArray* outArray = output->GetAtomData()->GetAbstractArray(inArray->GetName());
184 outArray->Resize(previousNbOfAtoms + nbOfAtoms);
185 }
186
187 for (vtkIdType ai = 0; ai < input->GetBondData()->GetNumberOfArrays(); ai++)
188 {
189 // skip bond orders array as it is auto-filled by AppendBond method
190 vtkAbstractArray* inArray = input->GetBondData()->GetAbstractArray(ai);
191 if (!strcmp(inArray->GetName(), input->GetBondOrdersArrayName()))
192 {
193 continue;
194 }
195 vtkAbstractArray* outArray = output->GetBondData()->GetAbstractArray(inArray->GetName());
196 outArray->Resize(previousNbOfBonds + nbOfBonds);
197 }
198
199 // --------------------
200 // Fill DataArrays
201 for (vtkIdType i = 0; i < input->GetNumberOfAtoms(); i++)
202 {
203 for (vtkIdType ai = 0; ai < input->GetAtomData()->GetNumberOfArrays(); ai++)
204 {
205 vtkAbstractArray* inArray = input->GetAtomData()->GetAbstractArray(ai);
206 vtkAbstractArray* outArray = output->GetAtomData()->GetAbstractArray(inArray->GetName());
207 // Use Value of non-ghost atom.
208 if (outputGhostAtoms && outputGhostAtoms->GetValue(atomIdMap[i]) == 0)
209 {
210 continue;
211 }
212 outArray->InsertTuple(atomIdMap[i], i, inArray);
213 }
214 }
215 for (vtkIdType i = 0; i < input->GetNumberOfBonds(); i++)
216 {
217 vtkBond bond = input->GetBond(i);
218 vtkIdType outputBondId =
219 output->GetBondId(atomIdMap[bond.GetBeginAtomId()], atomIdMap[bond.GetEndAtomId()]);
220
221 for (vtkIdType ai = 0; ai < input->GetBondData()->GetNumberOfArrays(); ai++)
222 {
223 // skip bond orders array as it is auto-filled by AppendBond method
224 vtkAbstractArray* inArray = input->GetBondData()->GetAbstractArray(ai);
225 if (!strcmp(inArray->GetName(), input->GetBondOrdersArrayName()))
226 {
227 continue;
228 }
229 vtkAbstractArray* outArray = output->GetBondData()->GetAbstractArray(inArray->GetName());
230 outArray->InsertTuple(outputBondId, i, inArray);
231 }
232 }
233 }
234
235 if (outputGhostBonds)
236 {
237 outputGhostBonds->SetNumberOfTuples(output->GetNumberOfBonds());
238 outputGhostBonds->Fill(0);
239 for (vtkIdType bondId = 0; bondId < output->GetNumberOfBonds(); bondId++)
240 {
241 vtkIdType atom1 = output->GetBond(bondId).GetBeginAtomId();
242 vtkIdType atom2 = output->GetBond(bondId).GetEndAtomId();
243 if (outputGhostAtoms->GetValue(atom1) == 1 || outputGhostAtoms->GetValue(atom2) == 1)
244 {
245 outputGhostBonds->SetValue(bondId, 1);
246 }
247 }
248 }
249
250 return 1;
251 }
252
253 //------------------------------------------------------------------------------
FillInputPortInformation(int i,vtkInformation * info)254 int vtkMoleculeAppend::FillInputPortInformation(int i, vtkInformation* info)
255 {
256 info->Set(vtkAlgorithm::INPUT_IS_REPEATABLE(), 1);
257 return this->Superclass::FillInputPortInformation(i, info);
258 }
259
260 //------------------------------------------------------------------------------
CheckArrays(vtkAbstractArray * array1,vtkAbstractArray * array2)261 bool vtkMoleculeAppend::CheckArrays(vtkAbstractArray* array1, vtkAbstractArray* array2)
262 {
263 if (strcmp(array1->GetName(), array2->GetName()) != 0)
264 {
265 vtkErrorMacro(<< "Execute: input name (" << array1->GetName() << "), must match output name ("
266 << array2->GetName() << ")");
267 return false;
268 }
269
270 if (array1->GetDataType() != array2->GetDataType())
271 {
272 vtkErrorMacro(<< "Execute: input ScalarType (" << array1->GetDataType()
273 << "), must match output ScalarType (" << array2->GetDataType() << ")");
274 return false;
275 }
276
277 if (array1->GetNumberOfComponents() != array2->GetNumberOfComponents())
278 {
279 vtkErrorMacro("Components of the inputs do not match");
280 return false;
281 }
282
283 return true;
284 }
285