Commit c2fff570 authored by Marcus D. Hanwell's avatar Marcus D. Hanwell Committed by Code Review
Browse files

Merge topic 'improved-ribbon-filter' into master

50eab6c3 Use vtkColor in the protein ribbon filter
166d57ce Improve the PDBReader to support secondary structures
parents 5abadb7a 50eab6c3
......@@ -25,6 +25,7 @@
#include "vtkRenderer.h"
#include "vtkProteinRibbonFilter.h"
#include "vtkPDBReader.h"
#include "vtkInteractorStyleSwitch.h"
int TestProteinRibbon(int argc, char *argv[])
{
......@@ -50,7 +51,6 @@ int TestProteinRibbon(int argc, char *argv[])
// setup actor
vtkNew<vtkActor> actor;
actor->SetMapper(polyDataMapper.GetPointer());
actor->GetProperty()->SetColor(0.0f, 1.0f, 0.0f);
// setup render window
vtkNew<vtkRenderer> ren;
......@@ -58,12 +58,18 @@ int TestProteinRibbon(int argc, char *argv[])
win->AddRenderer(ren.GetPointer());
vtkNew<vtkRenderWindowInteractor> iren;
iren->SetRenderWindow(win.GetPointer());
vtkInteractorStyleSwitch* is = vtkInteractorStyleSwitch::SafeDownCast(iren->GetInteractorStyle());
if (is)
{
is->SetCurrentStyleToTrackballCamera();
}
ren->AddActor(actor.GetPointer());
ren->SetBackground(0.0, 0.0, 0.0);
ren->SetBackground(0,0,0);
win->SetSize(450, 450);
win->Render();
ren->GetActiveCamera()->Zoom(2.2);
ren->ResetCamera();
ren->GetActiveCamera()->Zoom(1.5);
ren->ResetCameraClippingRange();
win->Render();
// Finally render the scene and compare the image to a reference image
win->SetMultiSamples(0);
......
......@@ -24,11 +24,59 @@
#include "vtkStringArray.h"
#include "vtkCellArray.h"
#include "vtkObjectFactory.h"
#include "vtkVector.h"
#include "vtkMath.h"
#include "vtkCellData.h"
#include "vtkPolyDataNormals.h"
#include "vtkSphereSource.h"
#include "vtkPeriodicTable.h"
#include "vtkMolecule.h"
#include "vtkVectorOperators.h"
#include <map>
vtkStandardNewMacro(vtkProteinRibbonFilter)
namespace
{
// Shamelessly copied from vtkColorSeries vtkColor3ubFromHex3, should be added
// to the vtkColor code.
vtkColor3ub ToColor3ubFromHex3(vtkTypeUInt32 hex)
{
int b = hex & 0xff;
hex >>= 8;
int g = hex & 0xff;
hex >>= 8;
int r = hex & 0xff;
return vtkColor3ub(r, g, b);
}
// Again, should be added to the vtkColor code.
vtkColor3ub ToColor3ubFromColor3f(const vtkColor3f& color)
{
return vtkColor3ub(static_cast<unsigned char>(color[0] * 255.0f),
static_cast<unsigned char>(color[1] * 255.0f),
static_cast<unsigned char>(color[2] * 255.0f));
}
} // End of anonymous namespace.
vtkProteinRibbonFilter::vtkProteinRibbonFilter()
{
this->CoilWidth = 0.3f;
this->HelixWidth = 1.3f;
this->SphereResolution = 60;
this->SubdivideFactor = 20;
this->ElementColors["H"] = ToColor3ubFromHex3(0xCCCCCC);
this->ElementColors["C"] = ToColor3ubFromHex3(0xAAAAAA);
this->ElementColors["O"] = ToColor3ubFromHex3(0xCC0000);
this->ElementColors["N"] = ToColor3ubFromHex3(0x0000CC);
this->ElementColors["S"] = ToColor3ubFromHex3(0xCCCC00);
this->ElementColors["P"] = ToColor3ubFromHex3(0x6622CC);
this->ElementColors["F"] = ToColor3ubFromHex3(0x00CC00);
this->ElementColors["CL"] = ToColor3ubFromHex3(0x00CC00);
this->ElementColors["BR"] = ToColor3ubFromHex3(0x882200);
this->ElementColors["I"] = ToColor3ubFromHex3(0x6600AA);
this->ElementColors["FE"] = ToColor3ubFromHex3(0xCC6600);
this->ElementColors["CA"] = ToColor3ubFromHex3(0xDDDDDD);
}
vtkProteinRibbonFilter::~vtkProteinRibbonFilter()
......@@ -51,56 +99,339 @@ int vtkProteinRibbonFilter::RequestData(vtkInformation *,
vtkInformationVector **inputVector,
vtkInformationVector *outputVector)
{
vtkPolyData *input = vtkPolyData::GetData(inputVector[0]);
vtkPolyData *output = vtkPolyData::GetData(outputVector);
// extract alpha-carbon backbone from input poly data
vtkStringArray *atomTypes =
vtkStringArray::SafeDownCast(
input->GetPointData()->GetAbstractArray("atom_types"));
if(!atomTypes)
vtkPolyData *input = vtkPolyData::GetData(inputVector[0]);
vtkPolyData *output = vtkPolyData::GetData(outputVector);
vtkPointData* pointData = input->GetPointData();
// Extract alpha-carbon backbone from input poly data
vtkStringArray *atomTypes =
vtkStringArray::SafeDownCast(pointData->GetAbstractArray("atom_types"));
if (!atomTypes)
{
vtkErrorMacro(<< "Atom Type String Array Required");
return 0;
}
// Extract secondary structures information from input poly data
vtkIdTypeArray *resi =
vtkIdTypeArray::SafeDownCast(pointData->GetAbstractArray("residue"));
vtkUnsignedCharArray *chain =
vtkUnsignedCharArray::SafeDownCast(pointData->GetAbstractArray("chain"));
vtkUnsignedCharArray *atom_ss =
vtkUnsignedCharArray::SafeDownCast(pointData->GetAbstractArray("secondary_structures"));
vtkUnsignedCharArray *atom_ssbegin =
vtkUnsignedCharArray::SafeDownCast(pointData->GetAbstractArray("secondary_structures_begin"));
vtkUnsignedCharArray *atom_ssend =
vtkUnsignedCharArray::SafeDownCast(pointData->GetAbstractArray("secondary_structures_end"));
vtkUnsignedCharArray *ishetatm =
vtkUnsignedCharArray::SafeDownCast(pointData->GetAbstractArray("ishetatm"));
if (!resi || !chain || !atom_ss || !atom_ssbegin || !atom_ssend || !ishetatm)
{
vtkErrorMacro(<< "Atom Secondary Structures Arrays Required");
return 0;
}
char currentChain = 0;
char ss = 0;
vtkIdType currentResi = 0;
vtkVector3f currentCA;
vtkVector3f prevCO;
bool hasPrevCO = false;
vtkNew<vtkPoints> strandPoints;
vtkNew<vtkPolyData> strand;
strand->Allocate();
strand->SetPoints(strandPoints.GetPointer());
vtkNew<vtkUnsignedCharArray> faceColors;
faceColors->SetName("RGB");
faceColors->SetNumberOfComponents(3);
// Initialize colors per point/atom
std::vector<vtkColor3ub> atomsColors;
this->SetColorByAtom(atomsColors, atomTypes);
this->SetColorByStructure(atomsColors, atomTypes, atom_ss,
ToColor3ubFromHex3(0xFF0080),
ToColor3ubFromHex3(0xFFC800));
std::vector<vtkColor3ub> colors;
std::vector<std::pair<vtkVector3f, bool> > borderPoints[2];
// Need this for radius / color lookups
vtkNew<vtkPeriodicTable> pTab;
for (int i = 0; i < input->GetNumberOfPoints(); i++)
{
vtkStdString type = atomTypes->GetValue(i);
if (ishetatm->GetValue(i))
{
vtkErrorMacro(<< "Atom Type String Array Required");
return 0;
if (type == "O")
{
continue;
}
unsigned short atomicNum = pTab->GetAtomicNumber(type);
vtkColor3f color = pTab->GetDefaultRGBTuple(atomicNum);
CreateAtomAsSphere(strand.GetPointer(), faceColors.GetPointer(),
input->GetPoint(i), ToColor3ubFromColor3f(color),
pTab->GetVDWRadius(atomicNum), 1.f);
}
else if (type == "CA")
{
// Create a ribbon between 2 CA atoms passing through each O atoms found in-between
double *xyz = input->GetPoint(i);
char atomChain = chain->GetValue(i);
vtkIdType atomResi = resi->GetValue(i);
vtkNew<vtkPoints> points;
for(int i = 0; i < input->GetNumberOfPoints(); i++)
if (currentChain != atomChain || currentResi + 1 != atomResi)
{
this->CreateThinStrip(strand.GetPointer(), faceColors.GetPointer(),
strandPoints.GetPointer(), borderPoints[0], borderPoints[1], colors);
borderPoints[0].clear();
borderPoints[1].clear();
colors.clear();
ss = 0;
hasPrevCO = false;
}
currentCA.Set(xyz[0], xyz[1], xyz[2]);
currentChain = atomChain;
currentResi = atomResi;
ss = atom_ss->GetValue(i);
colors.push_back(atomsColors[i]);
}
else if (type == "O")
{
vtkStdString type = atomTypes->GetValue(i);
if(type == "CA")
// Insert a new step in the next ribbon
double *xyz = input->GetPoint(i);
vtkVector3f p(xyz[0], xyz[1], xyz[2]);
p = (p - currentCA).Normalized() * ((ss == 'c') ? this->CoilWidth : this->HelixWidth);
if (hasPrevCO && p.Dot(prevCO) < 0)
{
points->InsertNextPoint(input->GetPoint(i));
p = p * -1.f;
}
hasPrevCO = true;
prevCO = p;
bool isSheet = (ss == 's');
borderPoints[0].push_back(
std::pair<vtkVector3f, bool>(currentCA - prevCO, isSheet));
borderPoints[1].push_back(
std::pair<vtkVector3f, bool>(currentCA + prevCO, isSheet));
}
}
// Create the last ribbon strip if needed
this->CreateThinStrip(strand.GetPointer(), faceColors.GetPointer(),
strandPoints.GetPointer(), borderPoints[0],
borderPoints[1], colors);
strand->GetCellData()->SetScalars(faceColors.GetPointer());
vtkNew<vtkCellArray> lines;
lines->InsertNextCell(points->GetNumberOfPoints());
for(int i = 0; i < points->GetNumberOfPoints(); i++)
// Compute the model normals
vtkNew<vtkPolyDataNormals> pdnormals;
pdnormals->SetInputData(strand.GetPointer());
pdnormals->SetFeatureAngle(60.0);
pdnormals->Update();
output->ShallowCopy(pdnormals->GetOutput());
return 1;
}
void vtkProteinRibbonFilter::SetColorByAtom(std::vector<vtkColor3ub>& colors,
vtkStringArray* atomTypes)
{
vtkNew<vtkPeriodicTable> pTab;
unsigned int len = atomTypes->GetNumberOfValues();
colors.resize(len);
for (unsigned int i = 0; i < len; i++)
{
colors[i] = this->ElementColors[atomTypes->GetValue(i)];
}
}
void vtkProteinRibbonFilter::SetColorByStructure(std::vector<vtkColor3ub>& colors,
vtkStringArray* atomTypes,
vtkUnsignedCharArray* ss,
const vtkColor3ub& helixColor,
const vtkColor3ub& sheetColor)
{
unsigned int len = atomTypes->GetNumberOfValues();
colors.resize(len);
for (unsigned int i = 0; i < len; i++)
{
if (ss->GetValue(i) == 's')
{
lines->InsertCellPoint(i);
colors[i] = sheetColor;
}
else if (ss->GetValue(i) == 'h')
{
colors[i] = helixColor;
}
}
}
vtkNew<vtkPolyData> backbone;
backbone->SetPoints(points.GetPointer());
backbone->SetLines(lines.GetPointer());
vtkNew<vtkSplineFilter> splineFilter;
splineFilter->SetNumberOfSubdivisions(750);
splineFilter->SetInputData(backbone.GetPointer());
void vtkProteinRibbonFilter::CreateAtomAsSphere(vtkPolyData* poly,
vtkUnsignedCharArray *faceColors,
double *pos,
const vtkColor3ub& color,
float radius, float scale)
{
// Create the sphere source at the atom size & position
vtkNew<vtkSphereSource> sphereSource;
sphereSource->SetPhiResolution(this->SphereResolution);
sphereSource->SetThetaResolution(this->SphereResolution);
sphereSource->SetCenter(pos);
sphereSource->SetRadius(radius * scale);
sphereSource->Update();
vtkNew<vtkTubeFilter> tubeFilter;
tubeFilter->SetInputConnection(splineFilter->GetOutputPort());
tubeFilter->SetRadius(0.4f);
tubeFilter->SetNumberOfSides(20);
tubeFilter->SetCapping(true);
tubeFilter->Update();
// Extract polydata from sphere
vtkPolyData *sphere = sphereSource->GetOutput();
vtkPoints *spherePoints = sphere->GetPoints();
vtkCellArray *spherePolys = sphere->GetPolys();
output->ShallowCopy(tubeFilter->GetOutput());
vtkPoints *points = poly->GetPoints();
// Get offset for the new point IDs that will be added to points
vtkIdType pointOffset = points->GetNumberOfPoints();
// Total number of new points
vtkIdType numPoints = spherePoints->GetNumberOfPoints();
vtkIdType numCellPoints, *cellPoints;
// Add new points
for (vtkIdType i = 0; i < numPoints; ++i)
{
points->InsertNextPoint(spherePoints->GetPoint(i));
}
return 1;
// Add new cells (polygons) that represent the sphere
spherePolys->InitTraversal();
while (spherePolys->GetNextCell(numCellPoints, cellPoints) != 0)
{
vtkIdType *newCellPoints = new vtkIdType[numCellPoints];
for (vtkIdType i = 0; i < numCellPoints; ++i)
{
// The new point ids should be offset by the pointOffset above
newCellPoints[i] = cellPoints[i] + pointOffset;
}
poly->InsertNextCell(VTK_TRIANGLE_STRIP, numCellPoints, newCellPoints);
for (int i = 0; i < 3; ++i)
{
faceColors->InsertNextValue(color[i]);
}
delete [] newCellPoints;
}
}
void vtkProteinRibbonFilter::CreateThinStrip(vtkPolyData* poly,
vtkUnsignedCharArray *faceColors,
vtkPoints* p,
std::vector<std::pair<vtkVector3f, bool> >& p1,
std::vector<std::pair<vtkVector3f, bool> >& p2,
std::vector<vtkColor3ub> &colors)
{
if (p1.size() < 2 || p2.size() < 2)
{
return;
}
// Get offset for the new point IDs that will be added to points
vtkIdType pointOffset = p->GetNumberOfPoints();
// Subdivide (smooth) the 2 ribbon borders
std::vector<vtkVector3f>* points1 = Subdivide(p1, this->SubdivideFactor);
std::vector<vtkVector3f>* points2 = Subdivide(p2, this->SubdivideFactor);
const int len = points1->size();
// Insert smoothed ribbon borders points into the polydata
for (int i = 0; i < len; i++)
{
p->InsertNextPoint((*points1)[i].GetData());
p->InsertNextPoint((*points2)[i].GetData());
}
delete points1;
delete points2;
// Fill in between the 2 ribbons borders with a single triangle strip
vtkIdType* connectivity = new vtkIdType[4 * (len - 1)];
for (int i = 0, j = 0; i < len - 1; i++, j += 4)
{
connectivity[j + 0] = pointOffset + 2 * i;
connectivity[j + 1] = pointOffset + 2 * i + 1;
connectivity[j + 2] = pointOffset + 2 * i + 2;
connectivity[j + 3] = pointOffset + 2 * i + 3;
vtkColor3ub color = colors[floor(0.5f + i /
static_cast<float>(this->SubdivideFactor))];
for (int k = 0; k < 4; ++k)
{
for (int ci = 0; ci < 3; ++ci)
{
faceColors->InsertNextValue(color[ci]);
}
}
}
poly->InsertNextCell(VTK_TRIANGLE_STRIP, 4 * (len - 1), connectivity);
}
std::vector<vtkVector3f>* vtkProteinRibbonFilter::Subdivide(std::vector<std::pair<vtkVector3f, bool> >& p,
int div)
{
std::vector<vtkVector3f>* ret = new std::vector<vtkVector3f>;
std::vector<vtkVector3f> points;
// Smoothing test
points.push_back(p[0].first);
for (int i = 1, lim = p.size() - 1; i < lim; i++)
{
vtkVector3f& p1 = p[i].first;
vtkVector3f& p2 = p[i+1].first;
if (p[i].second)
{
points.push_back((p1 + p2) * 0.5f);
}
else
{
points.push_back(p1);
}
}
points.push_back(p[p.size() - 1].first);
// Catmull-Rom subdivision
for (int i = -1, size = points.size(); i <= size - 3; i++)
{
vtkVector3f& p0 = points[(i == -1) ? 0 : i];
vtkVector3f& p1 = points[i+1];
vtkVector3f& p2 = points[i+2];
vtkVector3f& p3 = points[(i == size - 3) ? size - 1 : i + 3];
vtkVector3f v0 = (p2 - p0) * 0.5f;
vtkVector3f v1 = (p3 - p1) * 0.5f;
for (int j = 0; j < div; j++)
{
double t = 1.0 / div * j;
double t2 = t * t;
double x = p1.GetX() + t * v0.GetX()
+ t2 * (-3 * p1.GetX() + 3 * p2.GetX() - 2 * v0.GetX() - v1.GetX())
+ t2 * t * (2 * p1.GetX() - 2 * p2.GetX() + v0.GetX() + v1.GetX());
double y = p1.GetY() + t * v0.GetY()
+ t2 * (-3 * p1.GetY() + 3 * p2.GetY() - 2 * v0.GetY() - v1.GetY())
+ t2 * t * (2 * p1.GetY() - 2 * p2.GetY() + v0.GetY() + v1.GetY());
double z = p1.GetZ() + t * v0.GetZ()
+ t2 * (-3 * p1.GetZ() + 3 * p2.GetZ() - 2 * v0.GetZ() - v1.GetZ())
+ t2 * t * (2 * p1.GetZ() - 2 * p2.GetZ() + v0.GetZ() + v1.GetZ());
ret->push_back(vtkVector3f(x, y, z));
}
}
ret->push_back(points[points.size() - 1]);
return ret;
}
void vtkProteinRibbonFilter::PrintSelf(ostream& os, vtkIndent indent)
{
this->Superclass::PrintSelf(os, indent);
......
......@@ -18,12 +18,18 @@
// .NAME vtkProteinRibbonFilter - generates protein ribbons
// .SECTION Description
// vtkProteinRibbonFilter is an poly data algorithm which generates
// protein ribbons.
// vtkProteinRibbonFilter is a polydata algorithm that generates protein
// ribbons.
#include "vtkDomainsChemistryModule.h" // for export macro
#include "vtkPolyDataAlgorithm.h"
#include "vtkColor.h" // For vtkColor3ub.
#include <map> // For element to color map.
class vtkVector3f;
class vtkStringArray;
class VTKDOMAINSCHEMISTRY_EXPORT vtkProteinRibbonFilter
: public vtkPolyDataAlgorithm
{
......@@ -33,6 +39,18 @@ public:
static vtkProteinRibbonFilter* New();
vtkGetMacro(CoilWidth, float);
vtkSetMacro(CoilWidth, float);
vtkGetMacro(HelixWidth, float);
vtkSetMacro(HelixWidth, float);
vtkGetMacro(SphereResolution, int);
vtkSetMacro(SphereResolution, int);
vtkGetMacro(SubdivideFactor, int);
vtkSetMacro(SubdivideFactor, int);
protected:
vtkProteinRibbonFilter();
~vtkProteinRibbonFilter();
......@@ -43,6 +61,32 @@ protected:
vtkInformationVector **,
vtkInformationVector *);
void CreateThinStrip(vtkPolyData* poly, vtkUnsignedCharArray *faceColors,
vtkPoints* p, std::vector<std::pair<vtkVector3f, bool> >& p1,
std::vector<std::pair<vtkVector3f, bool> >& p2,
std::vector<vtkColor3ub> &colors);
void CreateAtomAsSphere(vtkPolyData* poly, vtkUnsignedCharArray *faceColors,
double *pos, const vtkColor3ub& color, float radius,
float scale);
static std::vector<vtkVector3f>* Subdivide(std::vector<std::pair<vtkVector3f, bool> >& p,
int div);
void SetColorByAtom( std::vector<vtkColor3ub>& colors, vtkStringArray* atomTypes);
void SetColorByStructure(std::vector<vtkColor3ub>& colors,
vtkStringArray* atomTypes, vtkUnsignedCharArray* ss,
const vtkColor3ub& helixColor,
const vtkColor3ub& sheetColor);
std::map<std::string, vtkColor3ub> ElementColors;
float CoilWidth;
float HelixWidth;
int SphereResolution;
int SubdivideFactor;
private:
vtkProteinRibbonFilter(const vtkProteinRibbonFilter&); // Not implemented.
void operator=(const vtkProteinRibbonFilter&); // Not implemented.
......
......@@ -128,6 +128,12 @@ vtkMoleculeReaderBase::vtkMoleculeReaderBase()
this->Points = NULL;
this->RGB = NULL;
this->Radii = NULL;
this->Chain = NULL;
this->Residue = NULL;
this->SecondaryStructures = NULL;
this->SecondaryStructuresBegin = NULL;
this->SecondaryStructuresEnd = NULL;
this->IsHetatm = NULL;
this->NumberOfAtoms = 0;
this->SetNumberOfInputPorts(0);
......@@ -159,6 +165,30 @@ vtkMoleculeReaderBase::~vtkMoleculeReaderBase()
{
this->Radii->Delete();
}
if(this->Chain)
{
this->Chain->Delete();
}
if(this->Residue)
{
this->Residue->Delete();
}
if(this->SecondaryStructures)
{
this->SecondaryStructures->Delete();
}
if(this->SecondaryStructuresBegin)
{
this->SecondaryStructuresBegin->Delete();
}
if(this->SecondaryStructuresEnd)
{
this->SecondaryStructuresEnd->Delete();
}
if(this->IsHetatm)
{
this->IsHetatm->Delete();
}
}
int vtkMoleculeReaderBase::RequestData(
......@@ -221,6 +251,72 @@ int vtkMoleculeReaderBase::ReadMolecule(FILE *fp, vtkPolyData *output)