vtkMolecule.cxx 13.7 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
/*=========================================================================

Program:   Visualization Toolkit
Module:    vtkMolecule.cxx
Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
All rights reserved.
See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

This software is distributed WITHOUT ANY WARRANTY; without even
the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
PURPOSE.  See the above copyright notice for more information.

=========================================================================*/
#include "vtkMolecule.h"

16
#include "vtkAbstractElectronicData.h"
17
18
19
#include "vtkDataSetAttributes.h"
#include "vtkEdgeListIterator.h"
#include "vtkIdTypeArray.h"
20
#include "vtkNew.h"
21
22
23
24
#include "vtkObjectFactory.h"
#include "vtkPlane.h"
#include "vtkPoints.h"
#include "vtkUnsignedShortArray.h"
25
#include "vtkFloatArray.h"
26
27
28
#include "vtkVector.h"
#include "vtkVectorOperators.h"

29
#include <cassert>
30
31
32
33
34
35

//----------------------------------------------------------------------------
vtkStandardNewMacro(vtkMolecule);

//----------------------------------------------------------------------------
vtkMolecule::vtkMolecule()
36
  : ElectronicData(NULL)
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
{
  this->Initialize();
}

//----------------------------------------------------------------------------
void vtkMolecule::Initialize()
{
  // Reset underlying data structure
  this->Superclass::Initialize();

  // Setup vertex data
  vtkDataSetAttributes *vertData = this->GetVertexData();
  vertData->AllocateArrays(1); // atomic nums

  // Atomic numbers
52
  vtkNew<vtkUnsignedShortArray> atomicNums;
53
54
  atomicNums->SetNumberOfComponents(1);
  atomicNums->SetName("Atomic Numbers");
55
  vertData->SetScalars(atomicNums.GetPointer());
56
57

  // Nuclear coordinates
58
59
60
  vtkPoints *points = vtkPoints::New();
  this->SetPoints(points);
  points->Delete();
61
62
63
64
65

  // Setup edge data
  vtkDataSetAttributes *edgeData = this->GetEdgeData();
  edgeData->AllocateArrays(1); // Bond orders

66
  vtkNew<vtkUnsignedShortArray> bondOrders;
67
68
  bondOrders->SetNumberOfComponents(1);
  bondOrders->SetName("Bond Orders");
69
  edgeData->SetScalars(bondOrders.GetPointer());
70
71
72

  this->UpdateBondList();

73
  // Electronic data
74
  this->SetElectronicData(NULL);
75

76
77
78
79
80
81
  this->Modified();
}

//----------------------------------------------------------------------------
vtkMolecule::~vtkMolecule()
{
82
  this->SetElectronicData(NULL);
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
}

//----------------------------------------------------------------------------
void vtkMolecule::PrintSelf(ostream &os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  vtkIndent subIndent = indent.GetNextIndent();

  os << indent << "Atoms:\n";
  for (vtkIdType i = 0; i < this->GetNumberOfAtoms(); ++i)
    {
    this->GetAtom(i).PrintSelf(os, subIndent);
    }

  os << indent << "Bonds:\n";
  for (vtkIdType i = 0; i < this->GetNumberOfBonds(); ++i)
    {
    os << subIndent << "===== Bond " << i << ": =====\n";
    this->GetBond(i).PrintSelf(os, subIndent);
    }
104
105
106
107
108
109
110
111
112
113

  os << indent << "Electronic Data:\n";
  if (this->ElectronicData)
    {
    this->ElectronicData->PrintSelf(os, subIndent);
    }
  else
    {
    os << subIndent << "Not set.\n";
    }
114
115
116
}

//----------------------------------------------------------------------------
117
118
vtkAtom vtkMolecule::AppendAtom(unsigned short atomicNumber,
                                const vtkVector3f &pos)
119
{
120
  vtkUnsignedShortArray *atomicNums = vtkArrayDownCast<vtkUnsignedShortArray>
121
122
123
124
125
126
127
128
    (this->GetVertexData()->GetScalars());

  assert(atomicNums);

  vtkIdType id;
  this->AddVertexInternal(0, &id);

  atomicNums->InsertValue(id, atomicNumber);
129
  vtkIdType coordID = this->Points->InsertNextPoint(pos.GetData());
130
  (void)coordID;
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
  assert("point ids synced with vertex ids" && coordID == id);

  this->Modified();
  return vtkAtom(this, id);
}

//----------------------------------------------------------------------------
vtkAtom vtkMolecule::GetAtom(vtkIdType atomId)
{
  assert(atomId >= 0 && atomId < this->GetNumberOfAtoms());

  vtkAtom atom (this, atomId);
  return atom;
}

//----------------------------------------------------------------------------
unsigned short vtkMolecule::GetAtomAtomicNumber(vtkIdType id)
{
  assert(id >= 0 && id < this->GetNumberOfAtoms());

151
  vtkUnsignedShortArray *atomicNums = vtkArrayDownCast<vtkUnsignedShortArray>
152
153
154
155
156
157
158
159
160
161
162
163
    (this->GetVertexData()->GetScalars());

  assert(atomicNums);

  return atomicNums->GetValue(id);
}

//----------------------------------------------------------------------------
void vtkMolecule::SetAtomAtomicNumber(vtkIdType id, unsigned short atomicNum)
{
  assert(id >= 0 && id < this->GetNumberOfAtoms());

164
  vtkUnsignedShortArray *atomicNums = vtkArrayDownCast<vtkUnsignedShortArray>
165
166
167
168
169
170
171
172
173
    (this->GetVertexData()->GetScalars());

  assert(atomicNums);

  atomicNums->SetValue(id, atomicNum);
  this->Modified();
}

//----------------------------------------------------------------------------
174
void vtkMolecule::SetAtomPosition(vtkIdType id, const vtkVector3f &pos)
175
176
{
  assert(id >= 0 && id < this->GetNumberOfAtoms());
177
  this->Points->SetPoint(id, pos.GetData());
178
179
180
181
  this->Modified();
}

//----------------------------------------------------------------------------
182
void vtkMolecule::SetAtomPosition(vtkIdType id, double x, double y, double z)
183
184
{
  assert(id >= 0 && id < this->GetNumberOfAtoms());
185
186
  this->Points->SetPoint(id, x, y, z);
  this->Modified();
187
188
189
}

//----------------------------------------------------------------------------
190
vtkVector3f vtkMolecule::GetAtomPosition(vtkIdType id)
191
192
{
  assert(id >= 0 && id < this->GetNumberOfAtoms());
193
  vtkFloatArray *positions = vtkArrayDownCast<vtkFloatArray>(this->Points->GetData());
194
195
196
  assert(positions != NULL);
  float *data = static_cast<float *>(positions->GetVoidPointer(id * 3));
  return vtkVector3f(data);
197
198
199
}

//----------------------------------------------------------------------------
200
void vtkMolecule::GetAtomPosition(vtkIdType id, float pos[3])
201
{
202
  vtkVector3f position = this->GetAtomPosition(id);
203
204
205
  pos[0] = position.GetX();
  pos[1] = position.GetY();
  pos[2] = position.GetZ();
206
207
208
209
210
211
212
213
214
}

//----------------------------------------------------------------------------
vtkIdType vtkMolecule::GetNumberOfAtoms()
{
  return this->GetNumberOfVertices();
}

//----------------------------------------------------------------------------
215
vtkBond vtkMolecule::AppendBond(const vtkIdType atom1, const vtkIdType atom2,
216
217
                             const unsigned short order)
{
218
  vtkUnsignedShortArray *bondOrders = vtkArrayDownCast<vtkUnsignedShortArray>
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
    (this->GetEdgeData()->GetScalars());
  assert(bondOrders);

  vtkEdgeType edgeType;
  this->AddEdgeInternal(atom1, atom2, false, 0, &edgeType);
  this->SetBondListDirty();

  vtkIdType id = edgeType.Id;
  bondOrders->InsertValue(id, order);
  this->Modified();
  return vtkBond(this, id, atom1, atom2);
}

//----------------------------------------------------------------------------
vtkBond vtkMolecule::GetBond(vtkIdType bondId)
{
  assert(bondId >= 0 && bondId < this->GetNumberOfBonds());

237
  vtkIdTypeArray *bonds = this->GetBondList();
238
239
240
241
242
243
244
245
246
247
  // An array with two components holding the bonded atom's ids
  vtkIdType *ids = bonds->GetPointer(2 * bondId);
  return vtkBond (this, bondId, ids[0], ids[1]);
}

//----------------------------------------------------------------------------
void vtkMolecule::SetBondOrder(vtkIdType bondId, unsigned short order)
{
  assert(bondId >= 0 && bondId < this->GetNumberOfBonds());

248
  vtkUnsignedShortArray *bondOrders = vtkArrayDownCast<vtkUnsignedShortArray>
249
250
251
252
253
254
255
256
257
258
259
260
261
    (this->GetEdgeData()->GetScalars());

  assert(bondOrders);

  this->Modified();
  return bondOrders->SetValue(bondId, order);
}

//----------------------------------------------------------------------------
unsigned short vtkMolecule::GetBondOrder(vtkIdType bondId)
{
  assert(bondId >= 0 && bondId < this->GetNumberOfBonds());

262
  vtkUnsignedShortArray *bondOrders = vtkArrayDownCast<vtkUnsignedShortArray>
263
264
265
266
267
268
269
270
271
272
273
274
275
    (this->GetEdgeData()->GetScalars());

  assert(bondOrders);

  return bondOrders->GetValue(bondId);
}

//----------------------------------------------------------------------------
double vtkMolecule::GetBondLength(vtkIdType bondId)
{
  assert(bondId >= 0 && bondId < this->GetNumberOfBonds());

  // Get list of bonds
276
  vtkIdTypeArray *bonds = this->GetBondList();
277
278
279
280
  // An array of length two holding the bonded atom's ids
  vtkIdType *ids = bonds->GetPointer(bondId);

  // Get positions
281
282
  vtkVector3f pos1 = this->GetAtomPosition(ids[0]);
  vtkVector3f pos2 = this->GetAtomPosition(ids[1]);
283
284
285
286
287
288
289
290
291
292
293
294
295

  return (pos2 - pos1).Norm();
}

//----------------------------------------------------------------------------
vtkPoints * vtkMolecule::GetAtomicPositionArray()
{
  return this->Points;
}

//----------------------------------------------------------------------------
vtkUnsignedShortArray * vtkMolecule::GetAtomicNumberArray()
{
296
  vtkUnsignedShortArray *atomicNums = vtkArrayDownCast<vtkUnsignedShortArray>
297
298
299
300
301
302
303
304
305
306
307
308
309
    (this->GetVertexData()->GetScalars());

  assert(atomicNums);

  return atomicNums;
}

//----------------------------------------------------------------------------
vtkIdType vtkMolecule::GetNumberOfBonds()
{
  return this->GetNumberOfEdges();
}

310
311
312
//----------------------------------------------------------------------------
vtkCxxSetObjectMacro(vtkMolecule, ElectronicData, vtkAbstractElectronicData);

313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
//----------------------------------------------------------------------------
void vtkMolecule::ShallowCopy(vtkDataObject *obj)
{
  vtkMolecule *m = vtkMolecule::SafeDownCast(obj);
  if (!m)
    {
    vtkErrorMacro("Can only shallow copy from vtkMolecule or subclass.");
    return;
    }
  this->ShallowCopyStructure(m);
  this->ShallowCopyAttributes(m);
}

//----------------------------------------------------------------------------
void vtkMolecule::DeepCopy(vtkDataObject *obj)
{
  vtkMolecule *m = vtkMolecule::SafeDownCast(obj);
  if (!m)
    {
David C. Lonie's avatar
David C. Lonie committed
332
    vtkErrorMacro("Can only deep copy from vtkMolecule or subclass.");
333
334
335
336
337
338
    return;
    }
  this->DeepCopyStructure(m);
  this->DeepCopyAttributes(m);
}

339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
//----------------------------------------------------------------------------
bool vtkMolecule::CheckedShallowCopy(vtkGraph *g)
{
  bool result = this->Superclass::CheckedShallowCopy(g);
  this->BondListIsDirty = true;
  return result;
}

//----------------------------------------------------------------------------
bool vtkMolecule::CheckedDeepCopy(vtkGraph *g)
{
  bool result = this->Superclass::CheckedDeepCopy(g);
  this->BondListIsDirty = true;
  return result;
}
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390

//----------------------------------------------------------------------------
void vtkMolecule::ShallowCopyStructure(vtkMolecule *m)
{
  this->CopyStructureInternal(m, false);
}

//----------------------------------------------------------------------------
void vtkMolecule::DeepCopyStructure(vtkMolecule *m)
{
  this->CopyStructureInternal(m, true);
}

//----------------------------------------------------------------------------
void vtkMolecule::ShallowCopyAttributes(vtkMolecule *m)
{
  this->CopyAttributesInternal(m, false);
}

//----------------------------------------------------------------------------
void vtkMolecule::DeepCopyAttributes(vtkMolecule *m)
{
  this->CopyAttributesInternal(m, true);
}

//----------------------------------------------------------------------------
void vtkMolecule::CopyStructureInternal(vtkMolecule *m, bool deep)
{
  // Call superclass
  if (deep)
    {
    this->Superclass::DeepCopy(m);
    }
  else
    {
    this->Superclass::ShallowCopy(m);
    }
391
  this->BondListIsDirty = true;
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
  }

//----------------------------------------------------------------------------
void vtkMolecule::CopyAttributesInternal(vtkMolecule *m, bool deep)
{
  if (deep)
    {
    if (m->ElectronicData)
      this->ElectronicData->DeepCopy(m->ElectronicData);
    }
  else
    {
    this->SetElectronicData(m->ElectronicData);
    }
}

408
//----------------------------------------------------------------------------
409
410
411
412
413
414
void vtkMolecule::UpdateBondList()
{
  this->BuildEdgeList();
  this->BondListIsDirty = false;
}

415
416
417
418
419
420
421
422
423
424
425
426
427
428
//----------------------------------------------------------------------------
vtkIdTypeArray *vtkMolecule::GetBondList()
{
  // Create the edge list if it doesn't exist, or is marked as dirty.
  vtkIdTypeArray *edgeList = this->BondListIsDirty ? NULL : this->GetEdgeList();
  if (!edgeList)
    {
    this->UpdateBondList();
    edgeList = this->GetEdgeList();
    }
  assert(edgeList != NULL);
  return edgeList;
}

429
430
431
432
433
434
435
436
437
438
439
440
//----------------------------------------------------------------------------
bool vtkMolecule::GetPlaneFromBond(const vtkBond &bond,
                                   const vtkVector3f &normal,
                                   vtkPlane *plane)
{
  return vtkMolecule::GetPlaneFromBond(bond.GetBeginAtom(), bond.GetEndAtom(),
                                       normal, plane);
}

//----------------------------------------------------------------------------
bool vtkMolecule::GetPlaneFromBond(const vtkAtom &atom1, const vtkAtom &atom2,
                                   const vtkVector3f &normal, vtkPlane *plane)
441
442
{
  if (plane == NULL)
443
    {
444
    return false;
445
    }
446

447
  vtkVector3f v(atom1.GetPosition() - atom2.GetPosition());
448

449
  vtkVector3f n_i(normal);
450
  vtkVector3f unitV(v.Normalized());
451
452
453

  // Check if vectors are (nearly) parallel
  if (unitV.Compare(n_i.Normalized(), 1e-7))
454
    {
455
    return false;
456
    }
457
458
459
460
461

  // calculate projection of n_i onto v
  // TODO Remove or restore this when scalar mult. is supported again
  // vtkVector3d proj (unitV * n_i.Dot(unitV));
  double n_iDotUnitV = n_i.Dot(unitV);
462
  vtkVector3f proj (unitV[0] * n_iDotUnitV, unitV[1] * n_iDotUnitV,
463
464
465
466
                    unitV[2] * n_iDotUnitV);
  // end vtkVector reimplementation TODO

  // Calculate actual normal:
467
  vtkVector3f realNormal(n_i - proj);
468
469

  // Create plane:
470
471
472
  vtkVector3f pos(atom1.GetPosition());
  plane->SetOrigin(pos.Cast<double>().GetData());
  plane->SetNormal(realNormal.Cast<double>().GetData());
473
474
  return true;
}