source: vtkvis/branches/1.7/Molecule.cpp @ 4778

Last change on this file since 4778 was 3841, checked in by ldelgass, 11 years ago

More trace verbosity fixes

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