Chaste  Build::
VesselSurfaceGenerator.cpp
1 /*
2 
3  Copyright (c) 2005-2015, University of Oxford.
4  All rights reserved.
5 
6  University of Oxford means the Chancellor, Masters and Scholars of the
7  University of Oxford, having an administrative office at Wellington
8  Square, Oxford OX1 2JD, UK.
9 
10  This file is part of Chaste.
11 
12  Redistribution and use in source and binary forms, with or without
13  modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26  ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29  GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30  HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31  LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32  OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34  */
35 
36 #define _BACKWARD_BACKWARD_WARNING_H 1 //Cut out the vtk deprecated warning for now (gcc4.3)
37 #include <vtkPolygon.h>
38 #include <stdlib.h>
39 #include "Exception.hpp"
40 #include "Vessel.hpp"
41 #include "VesselNode.hpp"
42 #include "VesselSegment.hpp"
43 #include "UblasCustomFunctions.hpp"
44 #include "UblasIncludes.hpp"
45 #include "DimensionalChastePoint.hpp"
46 #include "Part.hpp"
47 #include "BaseUnits.hpp"
48 
49 #include "VesselSurfaceGenerator.hpp"
50 
51 template<unsigned DIM>
53  mpVesselNetwork(pVesselNetwork),
54  mpSurface(vtkSmartPointer<vtkPolyData>::New()),
55  mReferenceLength(BaseUnits::Instance()->GetReferenceLengthScale())
56 {
57 }
58 
59 template<unsigned DIM>
61 {
62 }
63 
64 template<unsigned DIM>
65 std::vector<DimensionalChastePoint<DIM> > VesselSurfaceGenerator<DIM>::GetHoles()
66 {
67  std::vector<DimensionalChastePoint<DIM> > hole_locations;
68  std::vector<boost::shared_ptr<VesselSegment<DIM> > > segments = mpVesselNetwork->GetVesselSegments();
69  for (unsigned idx = 0; idx < segments.size(); idx++)
70  {
71  hole_locations.push_back(segments[idx]->GetMidPoint());
72  }
73  return hole_locations;
74 }
75 
76 template<unsigned DIM>
77 std::vector<std::vector<boost::shared_ptr<Polygon<DIM> > > > VesselSurfaceGenerator<DIM>::GetSurface()
78 {
79  if(DIM==2)
80  {
81  EXCEPTION("The surface generator currently only works in 3D");
82  }
83  // Define the global axes
84  c_vector<double, DIM> z_axis = unit_vector<double>(DIM,2);
85  c_vector<double, DIM> y_axis = unit_vector<double>(DIM,1);
86 
87  // Generate a surface for each segment
88  std::vector<std::vector<boost::shared_ptr<Polygon<DIM> > > > segment_polygons;
89  std::vector<boost::shared_ptr<VesselSegment<DIM> > > segments = mpVesselNetwork->GetVesselSegments();
90 
91  for (unsigned idx = 0; idx < segments.size(); idx++)
92  {
93  boost::shared_ptr<VesselNode<DIM> > p_start_node = segments[idx]->GetNode(0);
94  boost::shared_ptr<VesselNode<DIM> > p_end_node = segments[idx]->GetNode(1);
95  c_vector<double, DIM> segment_tangent = segments[idx]->GetUnitTangent();
96 
97  // Create the precursor points
98  std::vector<c_vector<double, DIM> > start_points = MakeCircle(p_start_node->GetRadius()/mReferenceLength);
99  std::vector<c_vector<double, DIM> > end_points = MakeCircle(p_end_node->GetRadius()/mReferenceLength);
100 
101  double angle = std::acos(inner_prod(z_axis, segment_tangent));
102  if (std::abs(inner_prod(z_axis, segment_tangent)) < 1.0 - 1.e-6)
103  {
104  c_vector<double, DIM> axis = VectorProduct(z_axis, segment_tangent);
105  RotateAboutAxis(start_points, axis, angle);
106  RotateAboutAxis(end_points, axis, angle);
107  }
108 
109  // Get the bisection planes at the nodes
110  std::vector<vtkSmartPointer<vtkPlane> > start_planes;
111  std::vector<vtkSmartPointer<vtkPlane> > end_planes;
112  c_vector<double, DIM> average_start_normal;
113  c_vector<double, DIM> average_end_normal;
114 
115  if (p_start_node->GetNumberOfSegments() == 1)
116  {
117  c_vector<double, DIM> node_location = p_start_node->rGetLocation().GetLocation(mReferenceLength);
118  vtkSmartPointer<vtkPlane> p_plane = vtkSmartPointer<vtkPlane>::New();
119  if(DIM==2)
120  {
121  p_plane->SetOrigin(node_location[0], node_location[1], node_location[2]);
122  p_plane->SetNormal(segment_tangent[0], segment_tangent[1], segment_tangent[2]);
123  }
124  else
125  {
126  p_plane->SetOrigin(node_location[0], node_location[1], 0.0);
127  p_plane->SetNormal(segment_tangent[0], segment_tangent[1], 0.0);
128  }
129  start_planes.push_back(p_plane);
130  }
131  else
132  {
133  std::vector<boost::shared_ptr<VesselSegment<DIM> > > node_segments = p_start_node->GetSegments();
134  for (unsigned jdx = 0; jdx < node_segments.size(); jdx++)
135  {
136  if (node_segments[jdx] != segments[idx])
137  {
138  c_vector<double, DIM> other_segment_tangent = node_segments[jdx]->GetUnitTangent();
139  if (node_segments[jdx]->GetNode(0) == p_start_node)
140  {
141  other_segment_tangent = -other_segment_tangent;
142  }
143 
144  average_start_normal += VectorProduct(segment_tangent, other_segment_tangent);
145  c_vector<double, DIM> node_location = p_start_node->rGetLocation().GetLocation(mReferenceLength);
146  vtkSmartPointer<vtkPlane> p_plane = vtkSmartPointer<vtkPlane>::New();
147  if(DIM==2)
148  {
149  p_plane->SetOrigin(node_location[0], node_location[1], 0.0);
150  }
151  else
152  {
153  p_plane->SetOrigin(node_location[0], node_location[1], node_location[2]);
154  }
155 
156  c_vector<double, DIM> bisection_vector = segment_tangent + other_segment_tangent;
157  bisection_vector /= norm_2(bisection_vector);
158  if(DIM==2)
159  {
160  p_plane->SetNormal(bisection_vector[0], bisection_vector[1], 0.0);
161  }
162  else
163  {
164  p_plane->SetNormal(bisection_vector[0], bisection_vector[1], bisection_vector[2]);
165  }
166  start_planes.push_back(p_plane);
167  }
168  }
169  }
170 
171  if (p_end_node->GetNumberOfSegments() == 1)
172  {
173  c_vector<double, DIM> node_location = p_end_node->rGetLocation().GetLocation(mReferenceLength);
174  vtkSmartPointer<vtkPlane> p_plane = vtkSmartPointer<vtkPlane>::New();
175  p_plane->SetOrigin(node_location[0], node_location[1], node_location[2]);
176  p_plane->SetNormal(segment_tangent[0], segment_tangent[1], segment_tangent[2]);
177  end_planes.push_back(p_plane);
178  }
179  else
180  {
181  std::vector<boost::shared_ptr<VesselSegment<DIM> > > node_segments = p_end_node->GetSegments();
182  for (unsigned jdx = 0; jdx < node_segments.size(); jdx++)
183  {
184  if (node_segments[jdx] != segments[idx])
185  {
186  c_vector<double, DIM> other_segment_tangent = node_segments[jdx]->GetUnitTangent();
187  if (node_segments[jdx]->GetNode(1) == p_end_node)
188  {
189  other_segment_tangent = -other_segment_tangent;
190  }
191  average_end_normal += VectorProduct(segment_tangent, other_segment_tangent);
192 
193  c_vector<double, DIM> node_location = p_end_node->rGetLocation().GetLocation(mReferenceLength);
194  vtkSmartPointer<vtkPlane> p_plane = vtkSmartPointer<vtkPlane>::New();
195  if(DIM==2)
196  {
197  p_plane->SetOrigin(node_location[0], node_location[1], 0.0);
198  }
199  else
200  {
201  p_plane->SetOrigin(node_location[0], node_location[1], node_location[2]);
202  }
203 
204  c_vector<double, DIM> bisection_vector = segment_tangent + other_segment_tangent;
205  bisection_vector /= norm_2(bisection_vector);
206  if(DIM==2)
207  {
208  p_plane->SetNormal(bisection_vector[0], bisection_vector[1], 0.0);
209  }
210  else
211  {
212  p_plane->SetNormal(bisection_vector[0], bisection_vector[1], bisection_vector[2]);
213  }
214  end_planes.push_back(p_plane);
215  }
216  }
217  }
218 
219  // Project the precursor points onto the first plane they hit
220  Translate(start_points, segments[idx]->GetMidPoint().GetLocation(mReferenceLength));
221  Translate(end_points, segments[idx]->GetMidPoint().GetLocation(mReferenceLength));
222 
223  std::vector<c_vector<double, DIM> > projected_start_points = start_points;
224  std::vector<c_vector<double, DIM> > projected_end_points = end_points;
225 
226  for (unsigned jdx = 0; jdx < start_planes.size(); jdx++)
227  {
228  if (jdx == 0)
229  {
230  ProjectOnPlane(projected_start_points, -segment_tangent, 2.0 * (segments[idx]->GetLength()/mReferenceLength), start_planes[jdx]);
231  }
232  else
233  {
234  std::vector<c_vector<double, DIM> > candidate_points = start_points;
235  ProjectOnPlane(candidate_points, -segment_tangent, 2.0 * (segments[idx]->GetLength()/mReferenceLength), start_planes[jdx]);
236  for (unsigned mdx = 0; mdx < projected_start_points.size(); mdx++)
237  {
238  if (norm_2(candidate_points[mdx] - start_points[mdx])
239  < norm_2(projected_start_points[mdx] - start_points[mdx]))
240  {
241  projected_start_points[mdx] = candidate_points[mdx];
242  }
243  }
244  }
245  }
246 
247  for (unsigned jdx = 0; jdx < end_planes.size(); jdx++)
248  {
249  if (jdx == 0)
250  {
251  ProjectOnPlane(projected_end_points, -segment_tangent, 2.0 * (segments[idx]->GetLength()/mReferenceLength),end_planes[jdx]);
252  }
253  else
254  {
255  std::vector<c_vector<double, DIM> > candidate_points = end_points;
256  ProjectOnPlane(candidate_points, -segment_tangent, 2.0 * (segments[idx]->GetLength()/mReferenceLength), end_planes[jdx]);
257  for (unsigned mdx = 0; mdx < projected_end_points.size(); mdx++)
258  {
259  if (norm_2(candidate_points[mdx] - end_points[mdx])
260  < norm_2(projected_end_points[mdx] - end_points[mdx]))
261  {
262  projected_end_points[mdx] = candidate_points[mdx];
263  }
264  }
265  }
266  }
267 
268  // Create the vertices and polygons
269  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > start_vertices;
270  std::vector<boost::shared_ptr<DimensionalChastePoint<DIM> > > end_vertices;
271 
272  for (unsigned jdx = 0; jdx < projected_start_points.size(); jdx++)
273  {
274  start_vertices.push_back(DimensionalChastePoint<DIM>::Create(projected_start_points[jdx], mReferenceLength));
275  }
276  for (unsigned jdx = 0; jdx < projected_end_points.size(); jdx++)
277  {
278  end_vertices.push_back(DimensionalChastePoint<DIM>::Create(projected_end_points[jdx], mReferenceLength));
279  }
280 
281  //
282  std::vector<boost::shared_ptr<Polygon<DIM> > > polygons;
283  for (unsigned jdx = 0; jdx < start_vertices.size(); jdx++)
284  {
285  unsigned index2;
286  if (jdx != start_vertices.size() - 1)
287  {
288  index2 = jdx + 1;
289  }
290  else
291  {
292  index2 = 0;
293  }
294  boost::shared_ptr<Polygon<DIM> > p_polygon = Polygon<DIM>::Create(start_vertices[jdx]);
295  p_polygon->AddVertex(start_vertices[index2]);
296  p_polygon->AddVertex(end_vertices[index2]);
297  p_polygon->AddVertex(end_vertices[jdx]);
298  polygons.push_back(p_polygon);
299  }
300  // If a node has connectivity one add an end-cap
301  if (segments[idx]->GetNode(0)->GetNumberOfSegments() == 1)
302  {
303  polygons.push_back(Polygon<DIM>::Create(start_vertices));
304  }
305 
306  if (segments[idx]->GetNode(1)->GetNumberOfSegments() == 1)
307  {
308  polygons.push_back(Polygon<DIM>::Create(end_vertices));
309  }
310  segment_polygons.push_back(polygons);
311  }
312  return segment_polygons;
313 }
314 
315 template<unsigned DIM>
316 std::vector<boost::shared_ptr<Polygon<DIM> > > VesselSurfaceGenerator<DIM>::GetSurfacePolygons()
317 {
318  std::vector<std::vector<boost::shared_ptr<Polygon<DIM> > > > segment_polygons = GetSurface();
319  std::vector<boost::shared_ptr<Polygon<DIM> > > polygons;
320  for (unsigned idx = 0; idx < segment_polygons.size(); idx++)
321  {
322  for (unsigned jdx = 0; jdx < segment_polygons[idx].size(); jdx++)
323  {
324  polygons.push_back(segment_polygons[idx][jdx]);
325  }
326  }
327  return polygons;
328 }
329 
330 template<unsigned DIM>
331 vtkSmartPointer<vtkPolyData> VesselSurfaceGenerator<DIM>::GetVtkSurface()
332 {
333  std::vector<std::vector<boost::shared_ptr<Polygon<DIM> > > > segment_polygons = GetSurface();
334 
335  // Add the polygons to a part
336  Part<DIM> part;
337  for (unsigned idx = 0; idx < segment_polygons.size(); idx++)
338  {
339  for (unsigned jdx = 0; jdx < segment_polygons[idx].size(); jdx++)
340  {
341  part.AddPolygon(segment_polygons[idx][jdx]->GetVertices(), true);
342  }
343  }
344 
345  mpSurface = part.GetVtk();
346  return mpSurface;
347 }
348 
349 template<unsigned DIM>
350 std::vector<c_vector<double, DIM> > VesselSurfaceGenerator<DIM>::MakeCircle(double radius, unsigned numberOfSegments)
351 {
352 
353  double increment = 2.0 * M_PI / double(numberOfSegments);
354  double angle = 0.0;
355  std::vector<c_vector<double, DIM> > points;
356 
357  for (unsigned idx = 0; idx < numberOfSegments; idx++)
358  {
359  c_vector<double, DIM> point;
360  point[0] = radius * std::cos(angle);
361  point[1] = radius * std::sin(angle);
362  point[2] = 0.0;
363  points.push_back(point);
364  angle += increment;
365  }
366  return points;
367 }
368 
369 template<unsigned DIM>
370 void VesselSurfaceGenerator<DIM>::ProjectOnPlane(std::vector<c_vector<double, DIM> >& rPoints,
371  c_vector<double, DIM> directionVector, double length,
372  vtkSmartPointer<vtkPlane> plane)
373 {
374  for (unsigned idx = 0; idx < rPoints.size(); idx++)
375  {
376  c_vector<double, DIM> point_on_line = rPoints[idx] + length * directionVector;
377  c_vector<double, DIM> projected_point;
378  double parametric_distance;
379  plane->IntersectWithLine(&rPoints[idx][0], &point_on_line[0], parametric_distance, &projected_point[0]);
380  rPoints[idx] = projected_point;
381  }
382 }
383 
384 template<unsigned DIM>
385 void VesselSurfaceGenerator<DIM>::RotateAboutAxis(std::vector<c_vector<double, DIM> >& rPoints, c_vector<double, DIM> axis,
386  double angle)
387 {
388  double sin_a = std::sin(angle);
389  double cos_a = std::cos(angle);
390  c_vector<double, DIM> unit_axis = axis / norm_2(axis);
391 
392  for (unsigned idx = 0; idx < rPoints.size(); idx++)
393  {
394  double dot_product = inner_prod(rPoints[idx], unit_axis);
395  c_vector<double, DIM> new_point;
396  new_point[0] = (unit_axis[0] * dot_product * (1.0 - cos_a) + rPoints[idx][0] * cos_a
397  + (-unit_axis[2] * rPoints[idx][1] + unit_axis[1] * rPoints[idx][2]) * sin_a);
398  new_point[1] = (unit_axis[1] * dot_product * (1.0 - cos_a) + rPoints[idx][1] * cos_a
399  + (unit_axis[2] * rPoints[idx][0] - unit_axis[0] * rPoints[idx][2]) * sin_a);
400  new_point[2] = (unit_axis[2] * dot_product * (1.0 - cos_a) + rPoints[idx][2] * cos_a
401  + (-unit_axis[1] * rPoints[idx][0] + unit_axis[0] * rPoints[idx][1]) * sin_a);
402 
403  rPoints[idx] = new_point;
404  }
405 }
406 
407 template<unsigned DIM>
408 void VesselSurfaceGenerator<DIM>::Translate(std::vector<c_vector<double, DIM> >& rPoints,
409  c_vector<double, DIM> translationVector)
410 {
411  for (unsigned idx = 0; idx < rPoints.size(); idx++)
412  {
413  rPoints[idx] += translationVector;
414  }
415 }
416 
417 template class VesselSurfaceGenerator<2> ;
418 template class VesselSurfaceGenerator<3> ;
vtkSmartPointer< vtkPolyData > GetVtk()
Definition: Part.cpp:511
std::vector< c_vector< double, DIM > > MakeCircle(double radius, unsigned numberOfSegments=16)
vtkSmartPointer< vtkPolyData > mpSurface
std::vector< std::vector< boost::shared_ptr< Polygon< DIM > > > > GetSurface()
vtkSmartPointer< vtkPolyData > GetVtkSurface()
VesselSurfaceGenerator(boost::shared_ptr< VesselNetwork< DIM > > pVesselNetwork)
void RotateAboutAxis(std::vector< c_vector< double, DIM > > &rPoints, c_vector< double, DIM > axis, double angle)
void Translate(std::vector< c_vector< double, DIM > > &rPoints, c_vector< double, DIM > translationVector)
units::quantity< unit::length > mReferenceLength
boost::shared_ptr< VesselNetwork< DIM > > mpVesselNetwork
static boost::shared_ptr< Polygon > Create(std::vector< boost::shared_ptr< DimensionalChastePoint< DIM > > > vertices)
Definition: Polygon.cpp:60
std::vector< boost::shared_ptr< Polygon< DIM > > > GetSurfacePolygons()
boost::shared_ptr< Polygon< DIM > > AddPolygon(std::vector< boost::shared_ptr< DimensionalChastePoint< DIM > > > vertices, bool newFacet=false, boost::shared_ptr< Facet< DIM > > pFacet=boost::shared_ptr< Facet< DIM > >())
Definition: Part.cpp:128
std::vector< DimensionalChastePoint< DIM > > GetHoles()
void ProjectOnPlane(std::vector< c_vector< double, DIM > > &rPoints, c_vector< double, DIM > directionVector, double length, vtkSmartPointer< vtkPlane > pPlane)
Definition: Part.hpp:60