source: vtkvis/trunk/Molecule.cpp @ 5073

Last change on this file since 5073 was 5073, checked in by ldelgass, 5 years ago

merge r5072 from release branch

  • Property svn:eol-style set to native
File size: 34.5 KB
Line 
1/* -*- mode: c++; c-basic-offset: 4; indent-tabs-mode: nil -*- */
2/*
3 * Copyright (C) 2004-2012  HUBzero Foundation, LLC
4 *
5 * Author: Leif Delgass <ldelgass@purdue.edu>
6 */
7
8#include <cstdio>
9#include <cfloat>
10#include <cassert>
11
12#include <vtkDataSet.h>
13#include <vtkCellArray.h>
14#include <vtkPointData.h>
15#include <vtkFloatArray.h>
16#include <vtkDoubleArray.h>
17#include <vtkIntArray.h>
18#include <vtkStringArray.h>
19#include <vtkPolyData.h>
20#include <vtkPolyDataMapper.h>
21#include <vtkActor.h>
22#include <vtkProperty.h>
23#include <vtkSphereSource.h>
24#include <vtkCylinderSource.h>
25#include <vtkTransform.h>
26#include <vtkTransformPolyDataFilter.h>
27#include <vtkGlyph3DMapper.h>
28#include <vtkPointSetToLabelHierarchy.h>
29#include <vtkLabelPlacementMapper.h>
30#include <vtkTextProperty.h>
31#include <vtkTransformPolyDataFilter.h>
32
33#include "Molecule.h"
34#include "MoleculeData.h"
35#include "Renderer.h"
36#include "Trace.h"
37
38using namespace VtkVis;
39
40Molecule::Molecule() :
41    GraphicsObject(),
42    _radiusScale(0.3),
43    _atomScaling(COVALENT_RADIUS),
44    _labelsOn(false),
45    _colorMap(NULL),
46    _colorMode(COLOR_BY_ELEMENTS),
47    _colorFieldType(DataSet::POINT_DATA)
48{
49    _bondColor[0] = _bondColor[1] = _bondColor[2] = 1.0f;
50    _colorFieldRange[0] = DBL_MAX;
51    _colorFieldRange[1] = -DBL_MAX;
52    _faceCulling = true;
53}
54
55Molecule::~Molecule()
56{
57#ifdef WANT_TRACE
58    if (_dataSet != NULL)
59        TRACE("Deleting Molecule for %s", _dataSet->getName().c_str());
60    else
61        TRACE("Deleting Molecule with NULL DataSet");
62#endif
63}
64
65void Molecule::setDataSet(DataSet *dataSet,
66                          Renderer *renderer)
67{
68    if (_dataSet != dataSet) {
69        _dataSet = dataSet;
70
71        _renderer = renderer;
72
73        if (renderer->getUseCumulativeRange()) {
74            renderer->getCumulativeDataRange(_dataRange,
75                                             _dataSet->getActiveScalarsName(),
76                                             1);
77            renderer->getCumulativeDataRange(_vectorMagnitudeRange,
78                                             _dataSet->getActiveVectorsName(),
79                                             3);
80            for (int i = 0; i < 3; i++) {
81                renderer->getCumulativeDataRange(_vectorComponentRange[i],
82                                                 _dataSet->getActiveVectorsName(),
83                                                 3, i);
84            }
85        } else {
86            _dataSet->getScalarRange(_dataRange);
87            _dataSet->getVectorRange(_vectorMagnitudeRange);
88            for (int i = 0; i < 3; i++) {
89                _dataSet->getVectorRange(_vectorComponentRange[i], i);
90            }
91        }
92
93        update();
94    }
95}
96
97/**
98 * \brief Create and initialize VTK Props to render a Molecule
99 */
100void Molecule::initProp()
101{
102    if (_atomProp == NULL) {
103        _atomProp = vtkSmartPointer<vtkActor>::New();
104        if (_faceCulling && _opacity == 1.0)
105            setCulling(_atomProp->GetProperty(), true);
106        _atomProp->GetProperty()->EdgeVisibilityOff();
107        _atomProp->GetProperty()->SetEdgeColor(_edgeColor[0], _edgeColor[1], _edgeColor[2]);
108        _atomProp->GetProperty()->SetLineWidth(_edgeWidth);
109        _atomProp->GetProperty()->SetOpacity(_opacity);
110        _atomProp->GetProperty()->SetAmbient(.2);
111        _atomProp->GetProperty()->SetSpecular(.2);
112        _atomProp->GetProperty()->SetSpecularPower(80.0);
113        if (!_lighting)
114            _atomProp->GetProperty()->LightingOff();
115    }
116    if (_bondProp == NULL) {
117        _bondProp = vtkSmartPointer<vtkActor>::New();
118        if (_faceCulling && _opacity == 1.0)
119            setCulling(_bondProp->GetProperty(), true);
120        _bondProp->GetProperty()->EdgeVisibilityOff();
121        _bondProp->GetProperty()->SetColor(_bondColor[0], _bondColor[1], _bondColor[2]);
122        _bondProp->GetProperty()->SetEdgeColor(_edgeColor[0], _edgeColor[1], _edgeColor[2]);
123        _bondProp->GetProperty()->SetLineWidth(_edgeWidth);
124        _bondProp->GetProperty()->SetOpacity(_opacity);
125        _bondProp->GetProperty()->SetAmbient(.2);
126        _bondProp->GetProperty()->SetSpecular(.2);
127        _bondProp->GetProperty()->SetSpecularPower(80.0);
128        if (!_lighting)
129            _bondProp->GetProperty()->LightingOff();
130    }
131    if (_labelProp == NULL) {
132        _labelProp = vtkSmartPointer<vtkActor2D>::New();
133    }
134    if (_prop == NULL) {
135        _prop = vtkSmartPointer<vtkAssembly>::New();
136    }
137}
138
139/**
140 * \brief Internal method to set up pipeline after a state change
141 */
142void Molecule::update()
143{
144    if (_dataSet == NULL) {
145        return;
146    }
147    vtkDataSet *ds = _dataSet->getVtkDataSet();
148
149    if (_atomMapper == NULL) {
150        _atomMapper = vtkSmartPointer<vtkGlyph3DMapper>::New();
151        _atomMapper->SetResolveCoincidentTopologyToPolygonOffset();
152        _atomMapper->ScalarVisibilityOn();
153    }
154    if (_bondMapper == NULL) {
155        _bondMapper = vtkSmartPointer<vtkGlyph3DMapper>::New();
156        _bondMapper->SetResolveCoincidentTopologyToPolygonOffset();
157        _bondMapper->ScalarVisibilityOn();
158    }
159    if (_labelMapper == NULL) {
160        _labelMapper = vtkSmartPointer<vtkLabelPlacementMapper>::New();
161        _labelMapper->SetShapeToRoundedRect();
162        _labelMapper->SetBackgroundColor(1.0, 1.0, 0.7);
163        _labelMapper->SetBackgroundOpacity(0.8);
164        _labelMapper->SetMargin(3);
165    }
166
167    if (_lut == NULL) {
168        if (ds->GetPointData() == NULL ||
169            ds->GetPointData()->GetScalars() == NULL ||
170            strcmp(ds->GetPointData()->GetScalars()->GetName(), "element") != 0) {
171            WARN("No element array in dataset %s", _dataSet->getName().c_str());
172            setColorMap(ColorMap::getDefault());
173            if (_atomScaling != NO_ATOM_SCALING)
174                _atomScaling = NO_ATOM_SCALING;
175        } else {
176            TRACE("Using element default colormap");
177            setColorMap(ColorMap::getElementDefault());
178        }
179    }
180
181    initProp();
182
183    addLabelArray(ds);
184    addRadiusArray(ds, _atomScaling, _radiusScale);
185    computeBonds(ds);
186
187    vtkPolyData *pd = vtkPolyData::SafeDownCast(ds);
188    if (pd) {
189        TRACE("Points: %d Verts: %d Lines: %d Polys: %d Strips: %d",
190              pd->GetNumberOfPoints(),
191              pd->GetNumberOfVerts(),
192              pd->GetNumberOfLines(),
193              pd->GetNumberOfPolys(),
194              pd->GetNumberOfStrips());
195        // DataSet is a vtkPolyData
196        if (pd->GetNumberOfLines() > 0) {
197            // Bonds
198            setupBondPolyData();
199
200            if (_cylinderSource == NULL) {
201                _cylinderSource = vtkSmartPointer<vtkCylinderSource>::New();
202                _cylinderSource->SetRadius(0.075);
203                _cylinderSource->SetHeight(1.0);
204                _cylinderSource->SetResolution(12);
205                _cylinderSource->CappingOff();
206                _cylinderTrans = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
207                _cylinderTrans->SetInputConnection(_cylinderSource->GetOutputPort());
208                vtkSmartPointer<vtkTransform> trans = vtkSmartPointer<vtkTransform>::New();
209                trans->RotateZ(-90.0);
210                _cylinderTrans->SetTransform(trans);
211            }
212            if (_lineSource == NULL) {
213                _lineSource = vtkSmartPointer<vtkLineSource>::New();
214                _lineSource->SetPoint1(-0.5, 0, 0);
215                _lineSource->SetPoint2(0.5, 0, 0);
216            }
217
218            _bondMapper->SetSourceConnection(_cylinderTrans->GetOutputPort());
219#ifdef USE_VTK6
220            _bondMapper->SetInputData(_bondPD);
221#else
222            _bondMapper->SetInputConnection(_bondPD->GetProducerPort());
223#endif
224            _bondMapper->SetOrientationArray("bond_orientations");
225            _bondMapper->SetOrientationModeToDirection();
226            _bondMapper->OrientOn();
227            _bondMapper->SetScaleArray("bond_scales");
228            _bondMapper->SetScaleModeToScaleByVectorComponents();
229            _bondMapper->ScalingOn();
230            _bondMapper->ClampingOff();
231
232            _bondProp->SetMapper(_bondMapper);
233            getAssembly()->AddPart(_bondProp);
234        }
235        if (pd->GetNumberOfPoints() > 0) {
236            if (_labelHierarchy == NULL) {
237                _labelHierarchy = vtkSmartPointer<vtkPointSetToLabelHierarchy>::New();
238            }
239            if (_labelTransform == NULL) {
240                _labelTransform = vtkSmartPointer<vtkTransform>::New();
241            }
242            vtkSmartPointer<vtkTransformPolyDataFilter> transformFilter = vtkSmartPointer<vtkTransformPolyDataFilter>::New();
243#ifdef USE_VTK6             
244            transformFilter->SetInputData(pd);
245#else
246            transformFilter->SetInput(pd);
247#endif
248            transformFilter->SetTransform(_labelTransform);
249            _labelHierarchy->SetInputConnection(transformFilter->GetOutputPort());
250            _labelHierarchy->SetLabelArrayName("_atom_labels");
251            _labelHierarchy->GetTextProperty()->SetColor(0, 0, 0);
252            _labelMapper->SetInputConnection(_labelHierarchy->GetOutputPort());
253            _labelProp->SetMapper(_labelMapper);
254
255            // Atoms
256            if (_sphereSource == NULL) {
257                _sphereSource = vtkSmartPointer<vtkSphereSource>::New();
258                _sphereSource->SetRadius(1.0);
259                _sphereSource->SetThetaResolution(14);
260                _sphereSource->SetPhiResolution(14);
261            }
262
263            _atomMapper->SetSourceConnection(_sphereSource->GetOutputPort());
264#ifdef USE_VTK6
265            _atomMapper->SetInputData(pd);
266#else
267            _atomMapper->SetInputConnection(pd->GetProducerPort());
268#endif
269            _atomMapper->SetScaleArray("_radii");
270            _atomMapper->SetScaleModeToScaleByMagnitude();
271            _atomMapper->ScalingOn();
272            _atomMapper->OrientOff();
273
274            _atomProp->SetMapper(_atomMapper);
275            getAssembly()->AddPart(_atomProp);
276        }
277    } else {
278        // DataSet is NOT a vtkPolyData
279        ERROR("DataSet is not a PolyData");
280        return;
281    }
282
283    setAtomLabelVisibility(_labelsOn);
284
285    if (pd->GetNumberOfPoints() > 0) {
286        _atomMapper->Update();
287        _labelMapper->Update();
288    }
289    if (pd->GetNumberOfLines() > 0) {
290        _bondMapper->Update();
291    }
292}
293
294void Molecule::updateLabelTransform()
295{
296    if (_labelTransform != NULL && getAssembly() != NULL) {
297        _labelTransform->SetMatrix(getAssembly()->GetMatrix());
298    }
299}
300
301void Molecule::updateRanges(Renderer *renderer)
302{
303    if (_dataSet == NULL) {
304        ERROR("called before setDataSet");
305        return;
306    }
307
308    if (renderer->getUseCumulativeRange()) {
309        renderer->getCumulativeDataRange(_dataRange,
310                                         _dataSet->getActiveScalarsName(),
311                                         1);
312        renderer->getCumulativeDataRange(_vectorMagnitudeRange,
313                                         _dataSet->getActiveVectorsName(),
314                                         3);
315        for (int i = 0; i < 3; i++) {
316            renderer->getCumulativeDataRange(_vectorComponentRange[i],
317                                             _dataSet->getActiveVectorsName(),
318                                             3, i);
319        }
320    } else {
321        _dataSet->getScalarRange(_dataRange);
322        _dataSet->getVectorRange(_vectorMagnitudeRange);
323        for (int i = 0; i < 3; i++) {
324            _dataSet->getVectorRange(_vectorComponentRange[i], i);
325        }
326    }
327
328    // Need to update color map ranges
329    double *rangePtr = _colorFieldRange;
330    if (_colorFieldRange[0] > _colorFieldRange[1]) {
331        rangePtr = NULL;
332    }
333    setColorMode(_colorMode, _colorFieldType, _colorFieldName.c_str(), rangePtr);
334}
335
336void Molecule::setColorMode(ColorMode mode)
337{
338    _colorMode = mode;
339    if (_dataSet == NULL)
340        return;
341
342    switch (mode) {
343    case COLOR_BY_ELEMENTS: // Assume default scalar is "element" array
344        setColorMode(mode,
345                     DataSet::POINT_DATA,
346                     "element");
347        break;
348    case COLOR_BY_SCALAR:
349        setColorMode(mode,
350                     _dataSet->getActiveScalarsType(),
351                     _dataSet->getActiveScalarsName());
352        break;
353    case COLOR_BY_VECTOR_MAGNITUDE:
354        setColorMode(mode,
355                     _dataSet->getActiveVectorsType(),
356                     _dataSet->getActiveVectorsName());
357        break;
358    case COLOR_BY_VECTOR_X:
359        setColorMode(mode,
360                     _dataSet->getActiveVectorsType(),
361                     _dataSet->getActiveVectorsName());
362        break;
363    case COLOR_BY_VECTOR_Y:
364        setColorMode(mode,
365                     _dataSet->getActiveVectorsType(),
366                     _dataSet->getActiveVectorsName());
367        break;
368    case COLOR_BY_VECTOR_Z:
369        setColorMode(mode,
370                     _dataSet->getActiveVectorsType(),
371                     _dataSet->getActiveVectorsName());
372        break;
373    case COLOR_CONSTANT:
374    default:
375        setColorMode(mode, DataSet::POINT_DATA, NULL, NULL);
376        break;
377    }
378}
379
380void Molecule::setColorMode(ColorMode mode,
381                            const char *name, double range[2])
382{
383    if (_dataSet == NULL)
384        return;
385    DataSet::DataAttributeType type = DataSet::POINT_DATA;
386    int numComponents = 1;
387    if (name != NULL && strlen(name) > 0 &&
388        !_dataSet->getFieldInfo(name, &type, &numComponents)) {
389        ERROR("Field not found: %s", name);
390        return;
391    }
392    setColorMode(mode, type, name, range);
393}
394
395void Molecule::setColorMode(ColorMode mode, DataSet::DataAttributeType type,
396                            const char *name, double range[2])
397{
398    _colorMode = mode;
399    _colorFieldType = type;
400    if (name == NULL)
401        _colorFieldName.clear();
402    else
403        _colorFieldName = name;
404    if (range == NULL) {
405        _colorFieldRange[0] = DBL_MAX;
406        _colorFieldRange[1] = -DBL_MAX;
407    } else {
408        memcpy(_colorFieldRange, range, sizeof(double)*2);
409    }
410
411    if (_dataSet == NULL || (_atomMapper == NULL && _bondMapper == NULL))
412        return;
413
414    switch (type) {
415    case DataSet::POINT_DATA:
416        if (_atomMapper != NULL)
417            _atomMapper->SetScalarModeToUsePointFieldData();
418        if (_bondMapper != NULL)
419            _bondMapper->SetScalarModeToUsePointFieldData();
420        break;
421    case DataSet::CELL_DATA:
422        if (_atomMapper != NULL)
423            _atomMapper->SetScalarModeToUseCellFieldData();
424        if (_bondMapper != NULL)
425            _bondMapper->SetScalarModeToUseCellFieldData();
426        break;
427    default:
428        ERROR("Unsupported DataAttributeType: %d", type);
429        return;
430    }
431
432    if (name != NULL && strlen(name) > 0) {
433        if (_atomMapper != NULL)
434            _atomMapper->SelectColorArray(name);
435        if (_bondMapper != NULL)
436            _bondMapper->SelectColorArray(name);
437    } else {
438        if (_atomMapper != NULL)
439            _atomMapper->SetScalarModeToDefault();
440        if (_bondMapper != NULL)
441            _bondMapper->SetScalarModeToDefault();
442    }
443
444    if (_lut != NULL) {
445        if (range != NULL) {
446            _lut->SetRange(range);
447        } else if (mode == COLOR_BY_ELEMENTS) {
448            double range[2];
449            range[0] = 0;
450            range[1] = NUM_ELEMENTS;
451            _lut->SetRange(range);
452        } else if (name != NULL && strlen(name) > 0) {
453            double r[2];
454            int comp = -1;
455            if (mode == COLOR_BY_VECTOR_X)
456                comp = 0;
457            else if (mode == COLOR_BY_VECTOR_Y)
458                comp = 1;
459            else if (mode == COLOR_BY_VECTOR_Z)
460                comp = 2;
461
462            if (_renderer->getUseCumulativeRange()) {
463                int numComponents;
464                if  (!_dataSet->getFieldInfo(name, type, &numComponents)) {
465                    ERROR("Field not found: %s, type: %d", name, type);
466                    return;
467                } else if (numComponents < comp+1) {
468                    ERROR("Request for component %d in field with %d components",
469                          comp, numComponents);
470                    return;
471                }
472                _renderer->getCumulativeDataRange(r, name, type, numComponents, comp);
473            } else {
474                _dataSet->getDataRange(r, name, type, comp);
475            }
476            _lut->SetRange(r);
477        }
478    }
479
480    switch (mode) {
481    case COLOR_BY_ELEMENTS:
482        if (_atomMapper != NULL)
483            _atomMapper->ScalarVisibilityOn();
484        if (_bondMapper != NULL)
485            _bondMapper->ScalarVisibilityOn();
486        break;
487    case COLOR_BY_SCALAR:
488        if (_atomMapper != NULL)
489            _atomMapper->ScalarVisibilityOn();
490        if (_bondMapper != NULL)
491            _bondMapper->ScalarVisibilityOn();
492        break;
493    case COLOR_BY_VECTOR_MAGNITUDE:
494        if (_atomMapper != NULL)
495            _atomMapper->ScalarVisibilityOn();
496        if (_bondMapper != NULL)
497            _bondMapper->ScalarVisibilityOn();
498        if (_lut != NULL) {
499            _lut->SetVectorModeToMagnitude();
500        }
501        break;
502    case COLOR_BY_VECTOR_X:
503        if (_atomMapper != NULL)
504            _atomMapper->ScalarVisibilityOn();
505        if (_bondMapper != NULL)
506            _bondMapper->ScalarVisibilityOn();
507        if (_lut != NULL) {
508            _lut->SetVectorModeToComponent();
509            _lut->SetVectorComponent(0);
510        }
511        break;
512    case COLOR_BY_VECTOR_Y:
513        if (_atomMapper != NULL)
514            _atomMapper->ScalarVisibilityOn();
515        if (_bondMapper != NULL)
516            _bondMapper->ScalarVisibilityOn();
517        if (_lut != NULL) {
518            _lut->SetVectorModeToComponent();
519            _lut->SetVectorComponent(1);
520        }
521        break;
522    case COLOR_BY_VECTOR_Z:
523        if (_atomMapper != NULL)
524            _atomMapper->ScalarVisibilityOn();
525        if (_bondMapper != NULL)
526            _bondMapper->ScalarVisibilityOn();
527        if (_lut != NULL) {
528            _lut->SetVectorModeToComponent();
529            _lut->SetVectorComponent(2);
530        }
531        break;
532    case COLOR_CONSTANT:
533    default:
534        if (_atomMapper != NULL)
535            _atomMapper->ScalarVisibilityOff();
536        if (_bondMapper != NULL)
537            _bondMapper->ScalarVisibilityOff();
538        break;
539    }
540}
541
542/**
543 * \brief Called when the color map has been edited
544 */
545void Molecule::updateColorMap()
546{
547    setColorMap(_colorMap);
548}
549
550/**
551 * \brief Associate a colormap lookup table with the DataSet
552 */
553void Molecule::setColorMap(ColorMap *cmap)
554{
555    if (cmap == NULL)
556        return;
557
558    _colorMap = cmap;
559 
560    if (_lut == NULL) {
561        _lut = vtkSmartPointer<vtkLookupTable>::New();
562        if (_atomMapper != NULL) {
563            _atomMapper->UseLookupTableScalarRangeOn();
564            _atomMapper->SetLookupTable(_lut);
565        }
566        if (_bondMapper != NULL) {
567            _bondMapper->UseLookupTableScalarRangeOn();
568            _bondMapper->SetLookupTable(_lut);
569        }
570        _lut->DeepCopy(cmap->getLookupTable());
571        switch (_colorMode) {
572        case COLOR_BY_ELEMENTS: {
573            double range[2];
574            range[0] = 0;
575            range[1] = NUM_ELEMENTS;
576            _lut->SetRange(range);
577        }
578            break;
579        case COLOR_CONSTANT:
580        case COLOR_BY_SCALAR:
581            _lut->SetRange(_dataRange);
582            break;
583        case COLOR_BY_VECTOR_MAGNITUDE:
584            _lut->SetRange(_vectorMagnitudeRange);
585            break;
586        case COLOR_BY_VECTOR_X:
587            _lut->SetRange(_vectorComponentRange[0]);
588            break;
589        case COLOR_BY_VECTOR_Y:
590            _lut->SetRange(_vectorComponentRange[1]);
591            break;
592        case COLOR_BY_VECTOR_Z:
593            _lut->SetRange(_vectorComponentRange[2]);
594            break;
595        default:
596            break;
597        }
598    } else {
599        double range[2];
600        _lut->GetTableRange(range);
601        _lut->DeepCopy(cmap->getLookupTable());
602        _lut->SetRange(range);
603        _lut->Modified();
604    }
605}
606
607void Molecule::setAtomLabelField(const char *fieldName)
608{
609    if (_labelHierarchy != NULL) {
610        if (strcmp(fieldName, "default") == 0) {
611            _labelHierarchy->SetLabelArrayName("_atom_labels");
612        } else {
613            _labelHierarchy->SetLabelArrayName(fieldName);
614        }
615    }
616}
617
618/**
619 * \brief Turn on/off rendering of atom labels
620 */
621void Molecule::setAtomLabelVisibility(bool state)
622{
623    _labelsOn = state;
624    if (_labelProp != NULL) {
625        _labelProp->SetVisibility((state ? 1 : 0));
626    }
627}
628
629/**
630 * \brief Turn on/off rendering of the atoms
631 */
632void Molecule::setAtomVisibility(bool state)
633{
634    if (_atomProp != NULL) {
635        _atomProp->SetVisibility((state ? 1 : 0));
636    }
637}
638
639/**
640 * \brief Turn on/off rendering of the bonds
641 */
642void Molecule::setBondVisibility(bool state)
643{
644    if (_bondProp != NULL) {
645        _bondProp->SetVisibility((state ? 1 : 0));
646    }
647}
648
649/**
650 * \brief Toggle visibility of the prop
651 */
652void Molecule::setVisibility(bool state)
653{
654    GraphicsObject::setVisibility(state);
655    if (_labelProp != NULL) {
656        if (!state)
657            _labelProp->SetVisibility(0);
658        else
659            setAtomLabelVisibility(_labelsOn);
660    }
661}
662
663/**
664 * \brief Set opacity of molecule
665 */
666void Molecule::setOpacity(double opacity)
667{
668    GraphicsObject::setOpacity(opacity);
669    if (_labelMapper != NULL) {
670        _labelMapper->SetBackgroundOpacity(opacity);
671    }
672}
673
674void Molecule::setAtomQuality(double quality)
675{
676    if (_sphereSource == NULL)
677        return;
678
679    if (quality > 10.0)
680        quality = 10.0;
681
682    int thetaRes = (int)(quality * 14.0);
683    int phiRes = (int)(quality * 14.0);
684    if (thetaRes < 4) thetaRes = 4;
685    if (phiRes < 3) phiRes = 3;
686
687    _sphereSource->SetThetaResolution(thetaRes);
688    _sphereSource->SetPhiResolution(phiRes);
689
690    if (_atomMapper != NULL) {
691        _atomMapper->Modified();
692        _atomMapper->Update();
693    }
694}
695
696void Molecule::setBondQuality(double quality)
697{
698    if (_cylinderSource == NULL)
699        return;
700
701    if (quality > 10.0)
702        quality = 10.0;
703
704    int res = (int)(quality * 12.0);
705    if (res < 3) res = 3;
706
707    _cylinderSource->SetResolution(res);
708
709    if (_bondMapper != NULL) {
710        _bondMapper->Modified();
711        _bondMapper->Update();
712    }
713}
714
715void Molecule::setBondStyle(BondStyle style)
716{
717    switch (style) {
718    case BOND_STYLE_CYLINDER:
719        if (_bondProp != NULL) {
720            _bondProp->GetProperty()->SetLineWidth(_edgeWidth);
721            _bondProp->GetProperty()->SetLighting(_lighting ? 1 : 0);
722        }
723        if (_bondMapper != NULL && _cylinderTrans != NULL &&
724            _bondMapper->GetInputConnection(1, 0) != _cylinderTrans->GetOutputPort()) {
725            _cylinderTrans->Modified();
726            _bondMapper->SetSourceConnection(_cylinderTrans->GetOutputPort());
727            _bondMapper->Modified();
728        }
729        break;
730    case BOND_STYLE_LINE:
731        if (_bondProp != NULL) {
732            _bondProp->GetProperty()->LightingOff();
733        }
734        if (_bondMapper != NULL && _lineSource != NULL &&
735            _bondMapper->GetInputConnection(1, 0) != _lineSource->GetOutputPort()) {
736            _lineSource->Modified();
737            _bondMapper->SetSourceConnection(_lineSource->GetOutputPort());
738            _bondMapper->Modified();
739        }
740         break;
741    default:
742        WARN("Unknown bond style");   
743    }
744}
745
746/**
747 * \brief Set constant bond color
748 */
749void Molecule::setBondColor(float color[3])
750{
751    _bondColor[0] = color[0];
752    _bondColor[1] = color[1];
753    _bondColor[2] = color[2];
754    if (_bondProp != NULL) {
755        _bondProp->GetProperty()->SetColor(_bondColor[0], _bondColor[1], _bondColor[2]);
756    }
757}
758
759/**
760 * \brief Set mode to determine how bonds are colored
761 */
762void Molecule::setBondColorMode(BondColorMode mode)
763{
764    if (_bondMapper == NULL) return;
765
766    switch (mode) {
767    case BOND_COLOR_BY_ELEMENTS:
768        _bondMapper->ScalarVisibilityOn();
769        break;
770    case BOND_COLOR_CONSTANT:
771        _bondMapper->ScalarVisibilityOff();
772        break;
773    default:
774        WARN("Unknown bond color mode");
775    }
776}
777
778/**
779 * \brief Set a group of world coordinate planes to clip rendering
780 *
781 * Passing NULL for planes will remove all cliping planes
782 */
783void Molecule::setClippingPlanes(vtkPlaneCollection *planes)
784{
785    if (_atomMapper != NULL) {
786        _atomMapper->SetClippingPlanes(planes);
787    }
788    if (_bondMapper != NULL) {
789        _bondMapper->SetClippingPlanes(planes);
790    }
791}
792
793/**
794 * \brief Set the radius type used for scaling atoms
795 */
796void Molecule::setAtomScaling(AtomScaling state)
797{
798    if (state == _atomScaling)
799        return;
800
801    _atomScaling = state;
802    if (_dataSet != NULL) {
803        vtkDataSet *ds = _dataSet->getVtkDataSet();
804        addRadiusArray(ds, _atomScaling, _radiusScale);
805        if (_atomMapper != NULL) {
806             assert(ds->GetPointData() != NULL &&
807                    ds->GetPointData()->GetArray("_radii") != NULL);
808            _atomMapper->SetScaleModeToScaleByMagnitude();
809            _atomMapper->SetScaleArray("_radii");
810            _atomMapper->ScalingOn();
811        }
812    }
813}
814
815/**
816 * \brief Set the constant radius scaling factor for atoms.  This
817 * can be used to convert from Angstroms to atom coordinates units.
818 */
819void Molecule::setAtomRadiusScale(double scale)
820{
821    if (scale == _radiusScale)
822        return;
823
824    _radiusScale = scale;
825    if (_dataSet != NULL) {
826        vtkDataSet *ds = _dataSet->getVtkDataSet();
827        addRadiusArray(ds, _atomScaling, _radiusScale);
828    }
829}
830
831/**
832 * \brief Set the constant radius scaling factor for bonds.
833 */
834void Molecule::setBondRadiusScale(double scale)
835{
836    if (_cylinderSource != NULL &&
837        _cylinderSource->GetRadius() != scale) {
838        _cylinderSource->SetRadius(scale);
839        // Workaround bug with source modification not causing
840        // mapper to be updated
841        if (_bondMapper != NULL) {
842            _bondMapper->Modified();
843        }
844    }
845}
846
847void Molecule::setupBondPolyData()
848{
849    if (_dataSet == NULL)
850        return;
851    if (_bondPD == NULL) {
852        _bondPD = vtkSmartPointer<vtkPolyData>::New();
853    } else {
854        _bondPD->Initialize();
855    }
856    vtkPolyData *pd = vtkPolyData::SafeDownCast(_dataSet->getVtkDataSet());
857    if (pd == NULL)
858        return;
859    vtkCellArray *lines = pd->GetLines();
860    lines->InitTraversal();
861    vtkIdType npts, *pts;
862    vtkSmartPointer<vtkPoints> bondPoints = vtkSmartPointer<vtkPoints>::New();
863    vtkSmartPointer<vtkIntArray> bondElements = vtkSmartPointer<vtkIntArray>::New();
864    vtkSmartPointer<vtkDoubleArray> bondVectors = vtkSmartPointer<vtkDoubleArray>::New();
865    vtkSmartPointer<vtkDoubleArray> bondScales = vtkSmartPointer<vtkDoubleArray>::New();
866    bondElements->SetName("element");
867    bondVectors->SetName("bond_orientations");
868    bondVectors->SetNumberOfComponents(3);
869    bondScales->SetName("bond_scales");
870    bondScales->SetNumberOfComponents(3);
871    vtkDataArray *elements = NULL;
872    if (pd->GetPointData() != NULL &&
873        pd->GetPointData()->GetScalars() != NULL &&
874        strcmp(pd->GetPointData()->GetScalars()->GetName(), "element") == 0) {
875        elements = pd->GetPointData()->GetScalars();
876    }
877    for (int i = 0; lines->GetNextCell(npts, pts); i++) {
878        assert(npts == 2);
879        double pt0[3], pt1[3];
880        double newPt0[3], newPt1[3];
881        //double *pt0 = pd->GetPoint(pts[0]);
882        //double *pt1 = pd->GetPoint(pts[1]);
883        pd->GetPoint(pts[0], pt0);
884        pd->GetPoint(pts[1], pt1);
885        double center[3];
886
887        for (int j = 0; j < 3; j++)
888            center[j] = pt0[j] + (pt1[j] - pt0[j]) * 0.5;
889        for (int j = 0; j < 3; j++)
890            newPt0[j] = pt0[j] + (pt1[j] - pt0[j]) * 0.25;
891        for (int j = 0; j < 3; j++)
892            newPt1[j] = pt0[j] + (pt1[j] - pt0[j]) * 0.75;
893
894        bondPoints->InsertNextPoint(newPt0);
895        bondPoints->InsertNextPoint(newPt1);
896
897#ifdef DEBUG
898        TRACE("Bond %d: (%g,%g,%g)-(%g,%g,%g)-(%g,%g,%g)", i,
899              pt0[0], pt0[1], pt0[2],
900              center[0], center[1], center[2],
901              pt1[0], pt1[1], pt1[2]);
902#endif
903
904        double vec[3];
905        for (int j = 0; j < 3; j++)
906            vec[j] = center[j] - pt0[j];
907
908        bondVectors->InsertNextTupleValue(vec);
909        bondVectors->InsertNextTupleValue(vec);
910#ifdef DEBUG
911        TRACE("Bond %d, vec: %g,%g,%g", i, vec[0], vec[1], vec[2]);
912#endif
913
914        double scale[3];
915        scale[0] = sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]);
916        scale[1] = 1.0;
917        scale[2] = 1.0;
918
919        bondScales->InsertNextTupleValue(scale);
920        bondScales->InsertNextTupleValue(scale);
921
922        if (elements != NULL) {
923            int element = (int)elements->GetComponent(pts[0], 0);
924#ifdef DEBUG
925            TRACE("Bond %d, elt 0: %d", i, element);
926#endif
927            bondElements->InsertNextValue(element);
928            element = (int)elements->GetComponent(pts[1], 0);
929#ifdef DEBUG
930            TRACE("Bond %d, elt 1: %d", i, element);
931#endif
932            bondElements->InsertNextValue(element);
933        }
934    }
935    _bondPD->SetPoints(bondPoints);
936    if (elements != NULL)
937        _bondPD->GetPointData()->SetScalars(bondElements);
938    _bondPD->GetPointData()->AddArray(bondVectors);
939    _bondPD->GetPointData()->AddArray(bondScales);
940}
941
942void Molecule::addLabelArray(vtkDataSet *dataSet)
943{
944    vtkDataArray *elements = NULL;
945    if (dataSet->GetPointData() != NULL &&
946        dataSet->GetPointData()->GetScalars() != NULL &&
947        strcmp(dataSet->GetPointData()->GetScalars()->GetName(), "element") == 0) {
948        elements = dataSet->GetPointData()->GetScalars();
949    } else {
950        WARN("Can't label atoms without an element array");
951    }
952
953    vtkSmartPointer<vtkStringArray> labelArray = vtkSmartPointer<vtkStringArray>::New();
954    labelArray->SetName("_atom_labels");
955    vtkPolyData *pd = vtkPolyData::SafeDownCast(dataSet);
956    if (pd == NULL) {
957        ERROR("DataSet not a PolyData");
958        return;
959    }
960    for (int i = 0; i < pd->GetNumberOfPoints(); i++) {
961        char buf[32];
962        if (elements != NULL) {
963            int elt = (int)elements->GetComponent(i, 0);
964            sprintf(buf, "%s%d", g_elementNames[elt], i);
965        } else {
966            sprintf(buf, "%d", i);
967        }
968        labelArray->InsertNextValue(buf);
969    }
970    dataSet->GetPointData()->AddArray(labelArray);
971}
972
973/**
974 * \brief Add a scalar array to dataSet with sizes for the elements
975 * specified in the "element" scalar array
976 */
977void Molecule::addRadiusArray(vtkDataSet *dataSet, AtomScaling scaling, double scaleFactor)
978{
979    vtkDataArray *elements = NULL;
980    if (dataSet->GetPointData() != NULL &&
981        dataSet->GetPointData()->GetScalars() != NULL &&
982        strcmp(dataSet->GetPointData()->GetScalars()->GetName(), "element") == 0) {
983        elements = dataSet->GetPointData()->GetScalars();
984    } else if (scaling != NO_ATOM_SCALING) {
985        WARN("Can't use non-constant scaling without an element array");
986    }
987    const float *radiusSource = NULL;
988    switch (scaling) {
989    case VAN_DER_WAALS_RADIUS:
990        radiusSource = g_vdwRadii;
991        break;
992    case COVALENT_RADIUS:
993        radiusSource = g_covalentRadii;
994        break;
995    case ATOMIC_RADIUS:
996        radiusSource = g_atomicRadii;
997        break;
998    case NO_ATOM_SCALING:
999    default:
1000        ;
1001    }
1002    vtkPolyData *pd = vtkPolyData::SafeDownCast(dataSet);
1003    if (pd == NULL) {
1004        ERROR("DataSet not a PolyData");
1005        return;
1006    }
1007    vtkDataArray *array = dataSet->GetPointData()->GetArray("_radii");
1008    vtkSmartPointer<vtkFloatArray> radii;
1009    if (array == NULL) {
1010        radii = vtkSmartPointer<vtkFloatArray>::New();
1011        radii->SetName("_radii");
1012        radii->SetNumberOfComponents(1);
1013        radii->SetNumberOfValues(pd->GetNumberOfPoints());
1014    } else {
1015        radii = vtkFloatArray::SafeDownCast(array);
1016        assert(radii != NULL);
1017    }
1018    for (int i = 0; i < pd->GetNumberOfPoints(); i++) {
1019        float value;
1020        if (elements != NULL && radiusSource != NULL) {
1021            int elt = (int)elements->GetComponent(i, 0);
1022            value = radiusSource[elt] * scaleFactor;
1023        } else {
1024            value = scaleFactor;
1025        }
1026        radii->SetValue(i, value);
1027    }
1028    radii->Modified();
1029    if (array == NULL) {
1030        dataSet->GetPointData()->AddArray(radii);
1031    }
1032}
1033
1034/**
1035 * \brief Create a color map to map atomic numbers to element colors
1036 */
1037ColorMap *Molecule::createElementColorMap()
1038{
1039    ColorMap *elementCmap = new ColorMap("elementDefault");
1040    ColorMap::ControlPoint cp[NUM_ELEMENTS+1];
1041
1042    elementCmap->setNumberOfTableEntries(NUM_ELEMENTS+1);
1043    for (int i = 0; i <= NUM_ELEMENTS; i++) {
1044        cp[i].value = i/((double)NUM_ELEMENTS);
1045        for (int c = 0; c < 3; c++) {
1046            cp[i].color[c] = ((double)g_elementColors[i][c])/255.;
1047        }
1048        elementCmap->addControlPoint(cp[i]);
1049    }
1050    ColorMap::OpacityControlPoint ocp[2];
1051    ocp[0].value = 0;
1052    ocp[0].alpha = 1.0;
1053    ocp[1].value = 1.0;
1054    ocp[1].alpha = 1.0;
1055    elementCmap->addOpacityControlPoint(ocp[0]);
1056    elementCmap->addOpacityControlPoint(ocp[1]);
1057    elementCmap->build();
1058    double range[2];
1059    range[0] = 0;
1060    range[1] = NUM_ELEMENTS;
1061    elementCmap->getLookupTable()->SetRange(range);
1062
1063    return elementCmap;
1064}
1065
1066int Molecule::computeBonds(vtkDataSet *dataSet, float tolerance)
1067{
1068    vtkPolyData *pd = vtkPolyData::SafeDownCast(dataSet);
1069    if (pd == NULL) {
1070        ERROR("DataSet not a PolyData");
1071        return -1;
1072    }
1073
1074    if (pd->GetNumberOfPoints() < 2 ||
1075        pd->GetNumberOfLines() > 0) {
1076        return 0;
1077    }
1078
1079    vtkDataArray *elements = NULL;
1080    if (pd->GetPointData() != NULL &&
1081        pd->GetPointData()->GetScalars() != NULL &&
1082        strcmp(pd->GetPointData()->GetScalars()->GetName(), "element") == 0) {
1083        elements = pd->GetPointData()->GetScalars();
1084    } else {
1085        TRACE("Can't compute bonds without an element array");
1086        return -1;
1087    }
1088    int numAtoms = pd->GetNumberOfPoints();
1089    vtkSmartPointer<vtkCellArray> cells = vtkSmartPointer<vtkCellArray>::New();
1090    for (int i = 0; i < numAtoms; i++) {
1091        double pt1[3];
1092        int elem1 = (int)elements->GetComponent(i, 0);
1093        float r1 = g_covalentRadii[elem1];
1094        pd->GetPoint(i, pt1);
1095        for (int j = i+1; j < numAtoms; j++) {
1096            float r2, dist;
1097            double pt2[3];
1098            int elem2 = (int)elements->GetComponent(j, 0);
1099
1100            if (elem1 == 1 && elem2 == 1)
1101                continue;
1102
1103            r2 = g_covalentRadii[elem2];
1104            pd->GetPoint(j, pt2);
1105
1106            double x = pt1[0] - pt2[0];
1107            double y = pt1[1] - pt2[1];
1108            double z = pt1[2] - pt2[2];
1109            dist = (float)sqrt(x*x + y*y + z*z);
1110            //TRACE("i=%d,j=%d elem1=%d,elem2=%d r1=%g,r2=%g, dist=%g", i, j, elem1, elem2, r1, r2, dist);
1111            if (dist < (r1 + r2 + tolerance)) {
1112                vtkIdType pts[2];
1113                pts[0] = i;
1114                pts[1] = j;
1115                cells->InsertNextCell(2, pts);
1116            }
1117        }
1118    }
1119
1120    if (cells->GetNumberOfCells() > 0) {
1121        pd->SetLines(cells);
1122    }
1123
1124    TRACE("Generated %d bonds", cells->GetNumberOfCells());
1125    return cells->GetNumberOfCells();
1126}
Note: See TracBrowser for help on using the repository browser.