source: vtkvis/branches/1.8/Molecule.cpp @ 6178

Last change on this file since 6178 was 5836, checked in by ldelgass, 9 years ago

merge r5835 from vtkvis trunk (remove support for VTK 5)

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