Skip to content

Commit

Permalink
refs #7261. Fix edge case.
Browse files Browse the repository at this point in the history
  • Loading branch information
OwenArnold committed Jun 7, 2013
1 parent 758b31e commit adfc199
Show file tree
Hide file tree
Showing 2 changed files with 41 additions and 14 deletions.
53 changes: 40 additions & 13 deletions Code/Mantid/Framework/Crystal/src/PeaksInRegion.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -120,7 +120,7 @@ namespace Crystal
setPropertySettings("PeakRadius", new EnabledWhenProperty("CheckPeakExtents", IS_NOT_DEFAULT) );
}

void validateExtentsInput(std::vector<double>& extents)
void validateExtentsInput(const std::vector<double>& extents)
{
if(extents.size() != 6)
{
Expand All @@ -140,6 +140,28 @@ namespace Crystal
}
}

bool pointOutsideAnyExtents(const V3D& testPoint, const std::vector<double>& extents)
{
return testPoint[0] < extents[0] || testPoint[0] > extents[1]
|| testPoint[1] < extents[2] || testPoint[1] > extents[3]
|| testPoint[2] < extents[4] || testPoint[2] > extents[5];
}

bool pointInsideAllExtents(const V3D& testPoint, const std::vector<double>& extents)
{
return testPoint[0] >= extents[0] && testPoint[0] <= extents[1]
&& testPoint[1] >= extents[2] && testPoint[1] <= extents[3]
&& testPoint[2]>= extents[4] && testPoint[2] <= extents[5];
}

void checkTouchPoint(const V3D& touchPoint,const V3D& normal,const V3D& peakCenter,const V3D& faceVertex, const double distance)
{
if( normal.scalar_prod(touchPoint - faceVertex) != 0)
{
throw std::runtime_error("Debugging. Calculation is wrong. touch point should always be on the plane!"); // Remove this line later. Check that geometry is setup properly.
}
}

//----------------------------------------------------------------------------------------------
/** Execute the algorithm.
*/
Expand Down Expand Up @@ -217,8 +239,8 @@ namespace Crystal
normals[i].normalize();
}


Progress prog(this, 0, nPeaks, nPeaks);
const int reportEveryNumber = 100;
Progress prog(this, 0, nPeaks/reportEveryNumber, nPeaks/reportEveryNumber);

PARALLEL_FOR2(ws, outputWorkspace)
for(int i = 0; i < nPeaks; ++i)
Expand All @@ -227,15 +249,12 @@ namespace Crystal
IPeak* peak = ws->getPeakPtr(i);
V3D peakCenter = coordFrameFunc(peak);

if(i%100 == 0)
prog.doReport();
if(i%reportEveryNumber == 0)
prog.report();

bool doesIntersect = true;
if (peakCenter[0] < extents[0] || peakCenter[0] >= extents[1]
|| peakCenter[1] < extents[2] || peakCenter[1] >= extents[3]
|| peakCenter[2] < extents[4] || peakCenter[2] >= extents[5])
if (pointOutsideAnyExtents(peakCenter, extents))
{

// Out of bounds.
doesIntersect = false;

Expand All @@ -244,11 +263,19 @@ namespace Crystal
// Take account of radius spherical extents.
for(int i = 0; i < 6; ++i)
{
double distance = normals[i].scalar_prod(peakCenter - faces[i][0]); // Distance between plane and peak center.
if(peakRadius >= std::abs(distance)) // Sphere passes through one of the faces, so intersects the box.
double distance = normals[i].scalar_prod(faces[i][0] - peakCenter); // Distance between plane and peak center.
if(peakRadius >= std::abs(distance)) // Sphere passes through one of the PLANES defined by the box faces.
{
doesIntersect = true;
break;
// Check that it is actually within the face boundaries.
V3D touchPoint = (normals[i] * distance) + peakCenter; // Vector equation of line give touch point on plane.

//checkTouchPoint(touchPoint, normals[i], peakCenter, faces[i][0], distance); // Debugging line.

if(pointInsideAllExtents(touchPoint, extents))
{
doesIntersect = true;
break;
}
}
}
}
Expand Down
2 changes: 1 addition & 1 deletion Code/Mantid/Framework/Crystal/test/PeaksInRegionTest.h
Original file line number Diff line number Diff line change
Expand Up @@ -340,7 +340,7 @@ class PeaksInRegionTest : public CxxTest::TestSuite
PeaksWorkspace_sptr ws = WorkspaceCreationHelper::createPeaksWorkspace(1);
auto detectorIds = ws->getInstrument()->getDetectorIDs();
Peak& peak = ws->getPeak(0);
peak.setHKL(Mantid::Kernel::V3D(0, 0, 2)); // This point is actually on the y = 0 plane, i.e. satisfies the plane equation. aX + bY + cZ = 0, but is outside the box.
peak.setHKL(Mantid::Kernel::V3D(2, 0, 0)); // This point is actually on the y = 0 plane, i.e. satisfies the plane equation. aX + bY + cZ = 0, but is outside the box.

const std::string outName = "OutWS";

Expand Down

0 comments on commit adfc199

Please sign in to comment.