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