message: ENH: add simulations of Forbild cones
diff: @@ -77,7 +77,7 @@ ForbildPhantomFileReader ::GenerateOutputInformation()
  CreateForbildElliptCyl(line, fig);
  else if (fig.substr(0, 9) == "Ellipsoid")
  CreateForbildEllipsoid(line, fig);
  - else if (fig.substr(0, 9) == "Cone")
  + else if (fig.substr(0, 4) == "Cone")
  CreateForbildCone(line, fig);
  // Density
  @@ -296,9 +296,83 @@ ForbildPhantomFileReader ::CreateForbildEllipsoid(const std::string & s, const s
  }
  void
  -ForbildPhantomFileReader ::CreateForbildCone(const std::string & /*s*/, const std::string & /*fig*/)
  +ForbildPhantomFileReader::CreateForbildCone(const std::string & s, const std::string & fig)
  {
  - itkExceptionMacro(<< "Cones have not been implemented (yet).");
  + ScalarType l = 0.;
  + if (!FindParameterInString("l", s, l))
  + itkExceptionMacro(<< "Could not find l (length) in " << s);
  + size_t found = 0;
  + ScalarType r1 = 0.;
  + if (!FindParameterInString("r1", s, r1))
  + {
  + itkExceptionMacro(<< "Missing radius r1 in " << fig << ", " << found << " found in " << s);
  + }
  + ScalarType r2 = 0.;
  + if (!FindParameterInString("r2", s, r2))
  + {
  + itkExceptionMacro(<< "Missing radius r2 in " << fig << ", " << found << " found in " << s);
  + }
  +
  + VectorType axes;
  + axes.Fill(r1);
  + VectorType planeDir;
  + planeDir[0] = 0.;
  + planeDir[1] = 0.;
  + planeDir[2] = 1.;
  +
  + QuadricShape::Pointer q = QuadricShape::New();
  + PointType spatialOffset;
  + spatialOffset[0] = 0.;
  + spatialOffset[1] = 0.;
  + if (r2 > r1)
  + {
  + axes[2] = -l * r1 / (r2 - r1);
  + q->AddClipPlane(-1. * planeDir, axes[2]);
  + q->AddClipPlane(planeDir, l - axes[2]);
  + spatialOffset[2] = axes[2] - 0.5 * l;
  + }
  + else
  + {
  + axes[2] = -l * r1 / (r1 - r2);
  + q->AddClipPlane(-1. * planeDir, -axes[2]);
  + q->AddClipPlane(planeDir, l + axes[2]);
  + spatialOffset[2] = -axes[2] - 0.5 * l;
  + }
  + q->SetEllipsoid(VectorType(0.), axes);
  + q->SetJ(0.);
  +
  + RotationMatrixType rot;
  + rot.Fill(0.);
  + if (fig == "Cone_x")
  + {
  + rot[0][2] = 1.;
  + rot[1][0] = 1.;
  + rot[2][1] = 1.;
  + }
  + else if (fig == "Cone_y")
  + {
  + rot[0][1] = 1.;
  + rot[1][2] = 1.;
  + rot[2][0] = 1.;
  + }
  + else if (fig == "Cone_z")
  + {
  + rot.SetIdentity();
  + }
  + else if (fig == "Cone")
  + {
  + VectorType dir;
  + if (!FindVectorInString("axis", s, dir))
  + itkExceptionMacro(<< "Could not find axis in " << s);
  + rot = ComputeRotationMatrixBetweenVectors(planeDir, dir);
  + }
  + else
  + {
  + itkExceptionMacro(<< "Unknown figure: " << fig);
  + }
  + q->Rotate(rot);
  + q->Translate(m_Center + rot * spatialOffset);
  + m_ConvexShape = q.GetPointer();
  }
  void
  