Commit b7908438 authored by Gigg, Martyn Anthony's avatar Gigg, Martyn Anthony
Browse files

Slight fix to cylinder triangulation and added cone triangulation as well. Updated

Object test also. Closes #693
parent 2fd3d78e
......@@ -104,6 +104,10 @@ class DLLExport Cone : public Quadratic
///This will get the bounding box for the cone
void getBoundingBox(double& xmax,double &ymax,double &zmax,double &xmin,double &ymin,double &zmin);
/// The number of slices to approximate a cone
static int g_nslices;
/// The number of stacks to approximate a cone
static int g_nstacks;
};
} // NAMESPACE MonteCarlo
......
......@@ -88,6 +88,9 @@ class DLLExport Object
double CylinderSolidAngle(const V3D & observer, const Mantid::Geometry::V3D & centre,
const Mantid::Geometry::V3D & axis,
const double radius, const double height) const;
double ConeSolidAngle(const V3D & observer, const Mantid::Geometry::V3D & centre,
const Mantid::Geometry::V3D & axis,
const double radius, const double height) const;
/// Geometry Handle for rendering
boost::shared_ptr<GeometryHandler> handle;
......
......@@ -84,7 +84,12 @@ namespace Mantid
void setBaseEqn();
///Writes the sphere equatation in MCNP format
void write(std::ostream&) const;
void getBoundingBox(double &xmax,double &ymax,double &zmax,double &xmin,double &ymin,double &zmin);
void getBoundingBox(double &xmax,double &ymax,double &zmax,double &xmin,double &ymin,double &zmin);
/// The number of slices to approximate a sphere
static int g_nslices;
/// The number of stacks to approximate a sphere
static int g_nstacks;
};
} // NAMESPACE Geometry
......
......@@ -28,8 +28,16 @@ namespace Mantid
namespace Geometry
{
Kernel::Logger& Cone::PLog(Kernel::Logger::get("Cone"));
// The number of slices to use to approximate a cylinder
int Cone::g_nslices = 10;
// The number of slices to use to approximate a cylinder
int Cone::g_nstacks = 1;
/// Floating point tolerance
//const double CTolerance(1e-6);
......
......@@ -2,6 +2,8 @@
#include "MantidGeometry/ObjComponent.h"
#include "MantidGeometry/Quat.h"
#include "MantidGeometry/Cylinder.h"
#include "MantidGeometry/Cone.h"
#include "MantidGeometry/Sphere.h"
#include <climits>
#ifdef _WIN32
#include "windows.h"
......@@ -142,7 +144,7 @@ namespace Mantid
gluQuadricNormals(qobj,GL_SMOOTH);
glPushMatrix();
glTranslated(center[0],center[1],center[2]);
gluSphere(qobj,radius,5,5);
gluSphere(qobj,radius, Geometry::Sphere::g_nslices, Geometry::Sphere::g_nstacks);
glPopMatrix();
gluDeleteQuadric(qobj);
}
......@@ -211,9 +213,9 @@ namespace Mantid
Quat rot(unit,axis);
rot.GLMatrix(mat);
glMultMatrixd(mat);
gluCylinder(qobj,0,radius,height,10,5);
gluCylinder(qobj,0,radius,height,Geometry::Cone::g_nslices,Geometry::Cone::g_nstacks);
glTranslated(0.0,0.0,height);
gluDisk(qobj,0,radius,10,1);
gluDisk(qobj,0,radius,Geometry::Cone::g_nslices,1);
glPopMatrix();
}
......
......@@ -31,6 +31,7 @@
#include "MantidGeometry/vtkGeometryCacheReader.h"
#include "MantidGeometry/vtkGeometryCacheWriter.h"
#include "MantidGeometry/Cylinder.h"
#include "MantidGeometry/Cone.h"
namespace Mantid
{
......@@ -933,7 +934,6 @@ namespace Mantid
xmin=ymin=zmin=-big;
xmax=ymax=zmax=big;
getBoundingBox(xmax,ymax,zmax,xmin,ymin,zmin);
std::cerr << "Ray trace bb " << xmin << " " << xmax << " " << ymin << " " << ymax << " " << zmin << " " << zmax << "\n";
// Is the bounding box a reasonable one?
if( xmax<big && ymax<big && zmax<big && xmin>-big && ymin>-big && zmin>-big )
{
......@@ -1099,6 +1099,8 @@ namespace Mantid
return SphereSolidAngle(observer, geometry_vectors, radius);
else if(type==3)
return CylinderSolidAngle(observer, geometry_vectors[0], geometry_vectors[1], radius, height);
else if(type==4)
return ConeSolidAngle(observer, geometry_vectors[0], geometry_vectors[1], radius, height);
return rayTraceSolidAngle(observer);
}
......@@ -1266,7 +1268,7 @@ namespace Mantid
}
/**
* Calculate the solid angle for a sphere using triangulation.
* Calculate the solid angle for a cylinder using triangulation.
* @param observer The observer's point
* @param centre The centre vector
* @param axis The axis vector
......@@ -1278,11 +1280,10 @@ namespace Mantid
const Mantid::Geometry::V3D & axis,
const double radius, const double height) const
{
// The cone is broken down into three pieces and then in turn broken down into triangles. Any triangle
// The cylinder is broken down into three pieces and then in turn broken down into triangles. Any triangle
// that has a normal facing away from the observer gives a negative solid angle and is excluded
// For simplicity the triangulation points are constructed such that the cone axis points up the +Z axis
// and then rotated into their final position
Geometry::V3D axis_direction = axis;
axis_direction.normalize();
// Required rotation
......@@ -1301,17 +1302,25 @@ namespace Mantid
for( int sl = 0; sl < nslices; ++sl )
{
int vertex = sl;
cos_table[vertex] = radius*cos(angle_step*vertex);
sin_table[vertex] = radius*sin(angle_step*vertex);
Geometry::V3D pt2 = Geometry::V3D(cos_table[vertex], sin_table[vertex], 0.0) + centre;
if( sl < nslices - 1 ) vertex = sl + 1;
cos_table[vertex] = radius*std::cos(angle_step*vertex);
sin_table[vertex] = radius*std::sin(angle_step*vertex);
Geometry::V3D pt2 = Geometry::V3D(cos_table[vertex], sin_table[vertex], 0.0);
if( sl < nslices - 1 )
{
vertex = sl + 1;
cos_table[vertex] = radius*std::cos(angle_step*vertex);
sin_table[vertex] = radius*std::sin(angle_step*vertex);
}
else vertex = 0;
Geometry::V3D pt3 = Geometry::V3D(cos_table[vertex], sin_table[vertex], 0.0) + centre;
Geometry::V3D pt3 = Geometry::V3D(cos_table[vertex], sin_table[vertex], 0.0);
transform.rotate(pt2);
transform.rotate(pt3);
double sa = getTriangleSolidAngle(centre, pt3, pt2, observer);
pt2 += centre;
pt3 += centre;
double sa = getTriangleSolidAngle(centre, pt2, pt3, observer);
if( sa > 0.0 )
{
solid_angle += sa;
......@@ -1319,22 +1328,25 @@ namespace Mantid
}
// Second the top cap
Geometry::V3D top_centre = Geometry::V3D(0.0, 0.0, height) + centre;
Geometry::V3D top_centre = Geometry::V3D(0.0, 0.0, height);
transform.rotate(top_centre);
top_centre += centre;
for( int sl = 0; sl < nslices; ++sl )
{
int vertex = sl;
Geometry::V3D pt2 = Geometry::V3D(cos_table[vertex], sin_table[vertex], height) + centre;
Geometry::V3D pt2 = Geometry::V3D(cos_table[vertex], sin_table[vertex], height);
if( sl < nslices - 1 ) vertex = sl + 1;
else vertex = 0;
Geometry::V3D pt3 = Geometry::V3D(cos_table[vertex], sin_table[vertex], height) + centre;
Geometry::V3D pt3 = Geometry::V3D(cos_table[vertex], sin_table[vertex], height);
// Rotate them to the correct axis orientation
transform.rotate(pt2);
transform.rotate(pt3);
double sa = getTriangleSolidAngle(top_centre, pt2, pt3, observer);
transform.rotate(pt3);
pt2 += centre;
pt3 += centre;
double sa = getTriangleSolidAngle(top_centre, pt3, pt2, observer);
if( sa > 0.0 )
{
solid_angle += sa;
......@@ -1353,28 +1365,33 @@ namespace Mantid
for( int sl = 0; sl < nslices; ++sl )
{
int vertex = sl;
Geometry::V3D pt1 = Geometry::V3D(cos_table[vertex], sin_table[vertex], z0) + centre;
Geometry::V3D pt1 = Geometry::V3D(cos_table[vertex], sin_table[vertex], z0);
if( sl < nslices - 1 ) vertex = sl + 1;
else vertex = 0;
Geometry::V3D pt3 = Geometry::V3D(cos_table[vertex], sin_table[vertex], z0) + centre;
Geometry::V3D pt3 = Geometry::V3D(cos_table[vertex], sin_table[vertex], z0);
vertex = sl;
Geometry::V3D pt2 = Geometry::V3D(cos_table[vertex], sin_table[vertex], z1) + centre;
Geometry::V3D pt2 = Geometry::V3D(cos_table[vertex], sin_table[vertex], z1);
if( sl < nslices - 1 ) vertex = sl + 1;
else vertex = 0;
Geometry::V3D pt4 = Geometry::V3D(cos_table[vertex], sin_table[vertex], z1) + centre;
Geometry::V3D pt4 = Geometry::V3D(cos_table[vertex], sin_table[vertex], z1);
// Rotations
transform.rotate(pt1);
transform.rotate(pt3);
transform.rotate(pt2);
transform.rotate(pt4);
double sa = getTriangleSolidAngle(pt1, pt3, pt4, observer);
pt1 += centre;
pt2 += centre;
pt3 += centre;
pt4 += centre;
double sa = getTriangleSolidAngle(pt1, pt4, pt3, observer);
if( sa > 0.0 )
{
solid_angle += sa;
}
sa = getTriangleSolidAngle(pt1, pt4, pt2, observer);
sa = getTriangleSolidAngle(pt1, pt2, pt4, observer);
if( sa > 0.0 )
{
solid_angle += sa;
......@@ -1385,6 +1402,154 @@ namespace Mantid
z1 += z_step;
}
delete [] cos_table;
delete [] sin_table;
return solid_angle;
}
/**
* Calculate the solid angle for a cone using triangulation.
* @param observer The observer's point
* @param centre The centre vector
* @param axis The axis vector
* @param radius The radius
* @param height The height
* @returns The solid angle value
*/
double Object::ConeSolidAngle(const V3D & observer, const Mantid::Geometry::V3D & centre,
const Mantid::Geometry::V3D & axis,
const double radius, const double height) const
{
// The cone is broken down into three pieces and then in turn broken down into triangles. Any triangle
// that has a normal facing away from the observer gives a negative solid angle and is excluded
// For simplicity the triangulation points are constructed such that the cone axis points up the +Z axis
// and then rotated into their final position
Geometry::V3D axis_direction = axis;
axis_direction.normalize();
// Required rotation
Geometry::V3D initial_axis = Geometry::V3D(0., 0., 1.0);
Geometry::V3D final_axis = axis_direction;
Geometry::Quat transform(initial_axis, final_axis);
// Do the base cap which is a point at the centre and nslices points around it
const int nslices(Mantid::Geometry::Cone::g_nslices);
const double angle_step = 2*M_PI/(double)nslices;
// Store the (x,y) points as they are used quite frequently
double *cos_table = new double[nslices];
double *sin_table = new double[nslices];
double solid_angle(0.0);
for( int sl = 0; sl < nslices; ++sl )
{
int vertex = sl;
cos_table[vertex] = std::cos(angle_step*vertex);
sin_table[vertex] = std::sin(angle_step*vertex);
Geometry::V3D pt2 = Geometry::V3D(radius*cos_table[vertex], radius*sin_table[vertex], 0.0);
if( sl < nslices - 1 )
{
vertex = sl + 1;
cos_table[vertex] = std::cos(angle_step*vertex);
sin_table[vertex] = std::sin(angle_step*vertex);
}
else vertex = 0;
Geometry::V3D pt3 = Geometry::V3D(radius*cos_table[vertex], radius*sin_table[vertex], 0.0);
transform.rotate(pt2);
transform.rotate(pt3);
pt2 += centre;
pt3 += centre;
double sa = getTriangleSolidAngle(centre, pt2, pt3, observer);
if( sa > 0.0 )
{
solid_angle += sa;
}
}
// Now the main section
const int nstacks(Mantid::Geometry::Cone::g_nstacks);
const double z_step = height / nstacks;
const double r_step = height / nstacks;
double z0(0.0), z1(z_step);
double r0(radius), r1(r0 - r_step);
for( int st = 1; st < nstacks; ++st )
{
if( st == nstacks ) z1 = height;
for( int sl = 0; sl < nslices; ++sl )
{
int vertex = sl;
Geometry::V3D pt1 = Geometry::V3D(r0*cos_table[vertex], r0*sin_table[vertex], z0);
if( sl < nslices - 1 ) vertex = sl + 1;
else vertex = 0;
Geometry::V3D pt3 = Geometry::V3D(r0*cos_table[vertex], r0*sin_table[vertex], z0);
vertex = sl;
Geometry::V3D pt2 = Geometry::V3D(r1*cos_table[vertex], r1*sin_table[vertex], z1);
if( sl < nslices - 1 ) vertex = sl + 1;
else vertex = 0;
Geometry::V3D pt4 = Geometry::V3D(r1*cos_table[vertex], r1*sin_table[vertex], z1);
// Rotations
transform.rotate(pt1);
transform.rotate(pt3);
transform.rotate(pt2);
transform.rotate(pt4);
pt1 += centre;
pt2 += centre;
pt3 += centre;
pt4 += centre;
double sa = getTriangleSolidAngle(pt1, pt4, pt3, observer);
if( sa > 0.0 )
{
solid_angle += sa;
}
sa = getTriangleSolidAngle(pt1, pt2, pt4, observer);
if( sa > 0.0 )
{
solid_angle += sa;
}
}
z0 = z1;
r0 = r1;
z1 += z_step;
r1 -= r_step;
}
// Top section
Geometry::V3D top_centre = Geometry::V3D(0.0, 0.0, height) + centre;
transform.rotate(top_centre);
top_centre += centre;
for( int sl = 0; sl < nslices; ++sl )
{
int vertex = sl;
Geometry::V3D pt2 = Geometry::V3D(r0*cos_table[vertex], r0*sin_table[vertex], height);
if( sl < nslices - 1 ) vertex = sl + 1;
else vertex = 0;
Geometry::V3D pt3 = Geometry::V3D(r0*cos_table[vertex], r0*sin_table[vertex], height);
// Rotate them to the correct axis orientation
transform.rotate(pt2);
transform.rotate(pt3);
pt2 += centre;
pt3 += centre;
double sa = getTriangleSolidAngle(top_centre, pt3, pt2, observer);
if( sa > 0.0 )
{
solid_angle += sa;
}
}
delete [] cos_table;
delete [] sin_table;
......
......@@ -31,6 +31,13 @@ Kernel::Logger& Sphere::PLog(Kernel::Logger::get("Sphere"));
//const double STolerance(1e-6); ///< Tolerance (should be replaced with boost::)
// The number of slices to use to approximate a sphere
int Sphere::g_nslices = 5;
// The number of slices to use to approximate a sphere
int Sphere::g_nstacks = 5;
Sphere::Sphere() : Quadratic(),
Centre(0,0,0),Radius(0.0)
/*!
......
......@@ -17,6 +17,8 @@
#include "MantidGeometry/SurfaceFactory.h"
#include "MantidGeometry/Track.h"
#include "MantidGeometry/GluGeometryHandler.h"
using namespace Mantid;
using namespace Geometry;
......@@ -564,29 +566,33 @@ public:
Test solid angle calculation for a capped cylinder
*/
{
Object A = createCappedCylinder();
double satol=2e-2; // tolerance for solid angle
Object A = createSmallCappedCylinder();
// Want to test triangulation so setup a geometry handler
boost::shared_ptr<GluGeometryHandler> h = boost::shared_ptr<GluGeometryHandler>(new GluGeometryHandler(&A));
h->setCylinder(V3D(-1.0,0.0,0.0), V3D(1., 0.0, 0.0), 0.005, 0.003);
A.setGeometryHandler(h);
// solid angle at distance 4 from capped cyl -3.2 1.2 in x, rad 3
double satol(1e-4); // tolerance for solid angle
// solid angle at point -0.5 from capped cyl -1.0 -0.997 in x, rad 0.005 - approx WISH cylinder
//
// soild angle of circle radius 3, distance 3 is 2pi(1-cos(t)) where
// t is atan(3/3), should be 1.840302,
// Work round for reverse track intercept has been made in Object.cpp
TS_ASSERT_DELTA(A.rayTraceSolidAngle(V3D(4.2,0,0)),1.840302,satol);
TS_ASSERT_DELTA(A.rayTraceSolidAngle(V3D(-7.2,0,0)),1.25663708,satol);
// t is atan(3/3), should be 0.000317939
TS_ASSERT_DELTA(A.triangleSolidAngle(V3D(-0.5, 0.0, 0.0)), 0.000317939, satol);
// Other end
TS_ASSERT_DELTA(A.triangleSolidAngle(V3D(-1.497, 0.0, 0.0)), 0.000317939, satol);
// No analytic value for side on SA, using hi-res value
TS_ASSERT_DELTA(A.rayTraceSolidAngle(V3D(0,0,7)),0.7531,satol);
TS_ASSERT_DELTA(A.rayTraceSolidAngle(V3D(0,7,0)),0.7531,satol);
TS_ASSERT_DELTA(A.triangleSolidAngle(V3D(0, 0, 0.1)), 8.03225e-05, satol);
TS_ASSERT_DELTA(A.triangleSolidAngle(V3D(0, 0.1, 0)), 8.03225e-05, satol);
// internal point (should be 4pi)
TS_ASSERT_DELTA(A.rayTraceSolidAngle(V3D(0,0,0)),4*M_PI,satol);
// surface point
TS_ASSERT_DELTA(A.rayTraceSolidAngle(V3D(1.2,0,0)),2*M_PI,satol);
TS_ASSERT_DELTA(A.rayTraceSolidAngle(V3D(-3.2,0,0)),2*M_PI,satol);
TS_ASSERT_DELTA(A.rayTraceSolidAngle(V3D(0,3,0)),2*M_PI,satol);
// distant points
TS_ASSERT_DELTA(A.rayTraceSolidAngle(V3D(20,0,0)),0.07850147,satol);
TS_ASSERT_DELTA(A.rayTraceSolidAngle(V3D(200,0,0)),0.000715295,satol);
TS_ASSERT_DELTA(A.rayTraceSolidAngle(V3D(2000,0,0)),7.08131e-6,satol);
TS_ASSERT_DELTA(A.triangleSolidAngle(V3D(-0.999, 0.0, 0.0)),4*M_PI,satol);
// surface points
TS_ASSERT_DELTA(A.triangleSolidAngle(V3D(-1.0, 0.0, 0.0)),2*M_PI,satol);
TS_ASSERT_DELTA(A.triangleSolidAngle(V3D(-0.997, 0.0, 0.0)),2*M_PI,satol);
}
void testSolidAngleCubeTriangles()
......@@ -813,6 +819,38 @@ private:
return retVal;
}
// This creates a cylinder to test the solid angle that is more realistic in size
// for a detector cylinder
Object createSmallCappedCylinder()
{
std::string C31="cx 0.005"; // cylinder x-axis radius 0.005 and height 0.003
std::string C32="px -0.997";
std::string C33="px -1.0";
// First create some surfaces
std::map<int,Surface*> CylSurMap;
CylSurMap[31]=new Cylinder();
CylSurMap[32]=new Plane();
CylSurMap[33]=new Plane();
CylSurMap[31]->setSurface(C31);
CylSurMap[32]->setSurface(C32);
CylSurMap[33]->setSurface(C33);
CylSurMap[31]->setName(31);
CylSurMap[32]->setName(32);
CylSurMap[33]->setName(33);
// Capped cylinder (id 21)
// using surface ids: 31 (cylinder) 32 (plane (top) ) and 33 (plane (base))
std::string ObjCapCylinder="-31 -32 33";
Object retVal;
retVal.setObject(21,ObjCapCylinder);
retVal.populate(CylSurMap);
return retVal;
}
Object createSphere()
{
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment