Shuenhoy/vtk-m
1
0
Fork 0
mirror of https://gitlab.kitware.com/vtk/vtk-m synced 2024-09-16 17:22:55 +00:00
Code Issues 2 Actions Packages Projects Releases Wiki Activity
330bdb9ff3
vtk-m/vtkm/exec/CellInterpolate.h
Kenneth Moreland cc497e6a1b Remove cont/Assert.h and exec/Assert.h 
These asserts are consolidated into the unified Assert.h. Also made some
minor edits to add asserts where appropriate and a little bit of
reconfiguring as found.
2016-04-20 15:41:14 -06:00

481 lines
18 KiB
C++
Raw Blame History

//============================================================================
// Copyright (c) Kitware, Inc.
// All rights reserved.
// See LICENSE.txt for details.
// This software is distributed WITHOUT ANY WARRANTY; without even
// the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
// PURPOSE. See the above copyright notice for more information.
//
// Copyright 2015 Sandia Corporation.
// Copyright 2015 UT-Battelle, LLC.
// Copyright 2015 Los Alamos National Security.
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Under the terms of Contract DE-AC52-06NA25396 with Los Alamos National
// Laboratory (LANL), the U.S. Government retains certain rights in
// this software.
//============================================================================
#ifndef vtk_m_exec_Interpolate_h
#define vtk_m_exec_Interpolate_h
#include <vtkm/Assert.h>
#include <vtkm/CellShape.h>
#include <vtkm/Math.h>
#include <vtkm/VecRectilinearPointCoordinates.h>
#include <vtkm/VectorAnalysis.h>
#include <vtkm/exec/FunctorBase.h>
namespace vtkm {
namespace exec {
namespace internal {
// This is really the WorldCoorindatesToParametericCoordinates function, but
// moved to this header file because it is required to interpolate in a
// polygon, which is divided into triangles.
template<typename WorldCoordVector>
VTKM_EXEC_EXPORT
typename WorldCoordVector::ComponentType
ReverseInterpolateTriangle(
const WorldCoordVector &pointWCoords,
const typename WorldCoordVector::ComponentType &wcoords)
{
VTKM_ASSERT(pointWCoords.GetNumberOfComponents() == 3);
// We will solve the world to parametric coordinates problem geometrically.
// Consider the parallelogram formed by wcoords and p0 of the triangle and
// the two adjacent edges. This parallelogram is equivalent to the
// axis-aligned rectangle anchored at the origin of parametric space.
//
// p2 |\ (1,0) |\ //
// | \ | \ //
// | \ | \ //
// | \ | \ //
// | \ | \ //
// | \ | (u,v) \ //
// | --- \ |-------* \ //
// | ---*wcoords | | \ //
// | | \ | | \ //
// p0 *--- | \ (0,0) *------------------\ (1,0) //
// ---| \ //
// x-- \ //
// --- \ //
// ---\ p1 //
//
// In this diagram, the distance between p0 and the point marked x divided by
// the length of the edge it is on is equal, by proportionality, to the u
// parametric coordiante. (The v coordinate follows the other edge
// accordingly.) Thus, if we can find the interesection at x (or more
// specifically the distance between p0 and x), then we can find that
// parametric coordinate.
//
// Because the triangle is in 3-space, we are actually going to intersect the
// edge with a plane that is parallel to the opposite edge of p0 and
// perpendicular to the triangle. This is partially because it is easy to
// find the intersection between a plane and a line and partially because the
// computation will work for points not on the plane. (The result is
// equivalent to a point projected on the plane.)
//
// First, we define an implicit plane as:
//
// dot((p - wcoords), planeNormal) = 0
//
// where planeNormal is the normal to the plane (easily computed from the
// triangle), and p is any point in the plane. Next, we define the parametric
// form of the line:
//
// p(d) = (p1 - p0)d + p0
//
// Where d is the fraction of distance from p0 toward p1. Note that d is
// actually equal to the parametric coordinate we are trying to find. Once we
// compute it, we are done. We can skip the part about finding the actual
// coordinates of the intersection.
//
// Solving for the interesection is as simple as substituting the line's
// definition of p(d) into p for the plane equation. With some basic algebra
// you get:
//
// d = dot((wcoords - p0), planeNormal)/dot((p1-p0), planeNormal)
//
// From here, the u coordiante is simply d. The v coordinate follows
// similarly.
//
typedef typename WorldCoordVector::ComponentType Vector3;
typedef typename Vector3::ComponentType T;
Vector3 pcoords(T(0));
Vector3 triangleNormal =
vtkm::TriangleNormal(pointWCoords[0], pointWCoords[1], pointWCoords[2]);
for (vtkm::IdComponent dimension = 0; dimension < 2; dimension++)
{
Vector3 p0 = pointWCoords[0];
Vector3 p1 = pointWCoords[dimension+1];
Vector3 p2 = pointWCoords[2-dimension];
Vector3 planeNormal = vtkm::Cross(triangleNormal, p2-p0);
T d = vtkm::dot(wcoords - p0, planeNormal)/vtkm::dot(p1 - p0, planeNormal);
pcoords[dimension] = d;
}
return pcoords;
}
}
/// \brief Interpolate a point field in a cell.
///
/// Given the point field values for each node and the parametric coordinates
/// of a point within the cell, interpolates the field to that point.
///
template<typename FieldVecType,
typename ParametricCoordType>
VTKM_EXEC_EXPORT
typename FieldVecType::ComponentType
CellInterpolate(const FieldVecType &pointFieldValues,
const vtkm::Vec<ParametricCoordType,3> &parametricCoords,
vtkm::CellShapeTagGeneric shape,
const vtkm::exec::FunctorBase &worklet)
{
typename FieldVecType::ComponentType result;
switch (shape.Id)
{
vtkmGenericCellShapeMacro(
result = CellInterpolate(pointFieldValues,
parametricCoords,
CellShapeTag(),
worklet));
default:
worklet.RaiseError("Unknown cell shape sent to interpolate.");
return typename FieldVecType::ComponentType();
}
return result;
}
//-----------------------------------------------------------------------------
template<typename FieldVecType,
typename ParametricCoordType>
VTKM_EXEC_EXPORT
typename FieldVecType::ComponentType
CellInterpolate(const FieldVecType &,
const vtkm::Vec<ParametricCoordType,3> &,
vtkm::CellShapeTagEmpty,
const vtkm::exec::FunctorBase &worklet)
{
worklet.RaiseError("Attempted to interpolate an empty cell.");
return typename FieldVecType::ComponentType();
}
//-----------------------------------------------------------------------------
template<typename FieldVecType,
typename ParametricCoordType>
VTKM_EXEC_EXPORT
typename FieldVecType::ComponentType
CellInterpolate(const FieldVecType &pointFieldValues,
const vtkm::Vec<ParametricCoordType,3>,
vtkm::CellShapeTagVertex,
const vtkm::exec::FunctorBase &vtkmNotUsed(worklet))
{
VTKM_ASSERT(pointFieldValues.GetNumberOfComponents() == 1);
return pointFieldValues[0];
}
//-----------------------------------------------------------------------------
template<typename FieldVecType,
typename ParametricCoordType>
VTKM_EXEC_EXPORT
typename FieldVecType::ComponentType
CellInterpolate(const FieldVecType &pointFieldValues,
const vtkm::Vec<ParametricCoordType,3> &parametricCoords,
vtkm::CellShapeTagLine,
const vtkm::exec::FunctorBase &vtkmNotUsed(worklet))
{
VTKM_ASSERT(pointFieldValues.GetNumberOfComponents() == 2);
return vtkm::Lerp(pointFieldValues[0],
pointFieldValues[1],
parametricCoords[0]);
}
template<typename ParametricCoordType>
VTKM_EXEC_EXPORT
vtkm::Vec<vtkm::FloatDefault,3>
CellInterpolate(const vtkm::VecRectilinearPointCoordinates<1> &field,
const vtkm::Vec<ParametricCoordType,3> &pcoords,
vtkm::CellShapeTagLine,
const vtkm::exec::FunctorBase &)
{
typedef vtkm::Vec<vtkm::FloatDefault,3> T;
const T &origin = field.GetOrigin();
const T &spacing = field.GetSpacing();
return T(origin[0] + static_cast<vtkm::FloatDefault>(pcoords[0])*spacing[0],
origin[1],
origin[2]);
}
//-----------------------------------------------------------------------------
template<typename FieldVecType,
typename ParametricCoordType>
VTKM_EXEC_EXPORT
typename FieldVecType::ComponentType
CellInterpolate(const FieldVecType &field,
const vtkm::Vec<ParametricCoordType,3> &pcoords,
vtkm::CellShapeTagTriangle,
const vtkm::exec::FunctorBase &vtkmNotUsed(worklet))
{
VTKM_ASSERT(field.GetNumberOfComponents() == 3);
typedef typename FieldVecType::ComponentType T;
return static_cast<T>( (field[0] * (1 - pcoords[0] - pcoords[1]))
+ (field[1] * pcoords[0])
+ (field[2] * pcoords[1]));
}
//-----------------------------------------------------------------------------
template<typename FieldVecType,
typename ParametricCoordType>
VTKM_EXEC_EXPORT
typename FieldVecType::ComponentType
CellInterpolate(const FieldVecType &field,
const vtkm::Vec<ParametricCoordType,3> &pcoords,
vtkm::CellShapeTagPolygon,
const vtkm::exec::FunctorBase &worklet)
{
const vtkm::IdComponent numPoints = field.GetNumberOfComponents();
VTKM_ASSERT(numPoints > 0);
switch (numPoints)
{
case 1:
return CellInterpolate(field,pcoords,vtkm::CellShapeTagVertex(),worklet);
case 2:
return CellInterpolate(field,pcoords,vtkm::CellShapeTagLine(),worklet);
case 3:
return CellInterpolate(field,pcoords,vtkm::CellShapeTagTriangle(),worklet);
case 4:
return CellInterpolate(field,pcoords,vtkm::CellShapeTagQuad(),worklet);
}
// If we are here, then there are 5 or more points on this polygon.
// Arrange the points such that they are on the circle circumscribed in the
// unit square from 0 to 1. That is, the point are on the circle centered at
// coordinate 0.5,0.5 with radius 0.5. The polygon is divided into regions
// defined by they triangle fan formed by the points around the center. This
// is C0 continuous but not necessarily C1 continuous. It is also possible to
// have a non 1 to 1 mapping between parametric coordinates world coordinates
// if the polygon is not planar or convex.
typedef typename FieldVecType::ComponentType FieldType;
// Find the interpolation for the center point.
FieldType fieldCenter = field[0];
for (vtkm::IdComponent pointIndex = 1; pointIndex < numPoints; pointIndex++)
{
fieldCenter = fieldCenter + field[pointIndex];
}
fieldCenter = fieldCenter*FieldType(1.0f/static_cast<float>(numPoints));
if ((vtkm::Abs(pcoords[0]-0.5f) < 4*vtkm::Epsilon<ParametricCoordType>()) &&
(vtkm::Abs(pcoords[1]-0.5f) < 4*vtkm::Epsilon<ParametricCoordType>()))
{
return fieldCenter;
}
ParametricCoordType angle = vtkm::ATan2(pcoords[1]-0.5f, pcoords[0]-0.5f);
if (angle < 0)
{
angle += static_cast<ParametricCoordType>(2*vtkm::Pi());
}
const ParametricCoordType deltaAngle =
static_cast<ParametricCoordType>(2*vtkm::Pi()/numPoints);
vtkm::IdComponent firstPointIndex =
static_cast<vtkm::IdComponent>(vtkm::Floor(angle/deltaAngle));
vtkm::IdComponent secondPointIndex = firstPointIndex + 1;
if (secondPointIndex == numPoints)
{
secondPointIndex = 0;
}
// Transform pcoords for polygon into pcoords for triangle.
vtkm::Vec<vtkm::Vec<ParametricCoordType,3>,3> polygonCoords;
polygonCoords[0][0] = 0.5f;
polygonCoords[0][1] = 0.5f;
polygonCoords[0][2] = 0;
polygonCoords[1][0] = 0.5f*(vtkm::Cos(deltaAngle*static_cast<ParametricCoordType>(firstPointIndex))+1);
polygonCoords[1][1] = 0.5f*(vtkm::Sin(deltaAngle*static_cast<ParametricCoordType>(firstPointIndex))+1);
polygonCoords[1][2] = 0.0f;
polygonCoords[2][0] = 0.5f*(vtkm::Cos(deltaAngle*static_cast<ParametricCoordType>(secondPointIndex))+1);
polygonCoords[2][1] = 0.5f*(vtkm::Sin(deltaAngle*static_cast<ParametricCoordType>(secondPointIndex))+1);
polygonCoords[2][2] = 0.0f;
vtkm::Vec<ParametricCoordType,3> trianglePCoords =
vtkm::exec::internal::ReverseInterpolateTriangle(polygonCoords, pcoords);
// Set up parameters for triangle that pcoords is in
vtkm::Vec<FieldType,3> triangleField;
triangleField[0] = fieldCenter;
triangleField[1] = field[firstPointIndex];
triangleField[2] = field[secondPointIndex];
// Now use the triangle interpolate
return vtkm::exec::CellInterpolate(triangleField,
trianglePCoords,
vtkm::CellShapeTagTriangle(),
worklet);
}
//-----------------------------------------------------------------------------
template<typename FieldVecType,
typename ParametricCoordType>
VTKM_EXEC_EXPORT
typename FieldVecType::ComponentType
CellInterpolate(const FieldVecType &field,
const vtkm::Vec<ParametricCoordType,3> &pcoords,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase &vtkmNotUsed(worklet))
{
VTKM_ASSERT(field.GetNumberOfComponents() == 4);
typedef typename FieldVecType::ComponentType T;
T bottomInterp = vtkm::Lerp(field[0], field[1], pcoords[0]);
T topInterp = vtkm::Lerp(field[3], field[2], pcoords[0]);
return vtkm::Lerp(bottomInterp, topInterp, pcoords[1]);
}
template<typename ParametricCoordType>
VTKM_EXEC_EXPORT
vtkm::Vec<vtkm::FloatDefault,3>
CellInterpolate(const vtkm::VecRectilinearPointCoordinates<2> &field,
const vtkm::Vec<ParametricCoordType,3> &pcoords,
vtkm::CellShapeTagQuad,
const vtkm::exec::FunctorBase &)
{
typedef vtkm::Vec<vtkm::FloatDefault,3> T;
const T &origin = field.GetOrigin();
const T &spacing = field.GetSpacing();
return T(origin[0] + static_cast<vtkm::FloatDefault>(pcoords[0])*spacing[0],
origin[1] + static_cast<vtkm::FloatDefault>(pcoords[1])*spacing[1],
origin[2]);
}
//-----------------------------------------------------------------------------
template<typename FieldVecType,
typename ParametricCoordType>
VTKM_EXEC_EXPORT
typename FieldVecType::ComponentType
CellInterpolate(const FieldVecType &field,
const vtkm::Vec<ParametricCoordType,3> &pcoords,
vtkm::CellShapeTagTetra,
const vtkm::exec::FunctorBase &vtkmNotUsed(worklet))
{
VTKM_ASSERT(field.GetNumberOfComponents() == 4);
typedef typename FieldVecType::ComponentType T;
return static_cast<T>( (field[0] * (1-pcoords[0]-pcoords[1]-pcoords[2]))
+ (field[1] * pcoords[0])
+ (field[2] * pcoords[1])
+ (field[3] * pcoords[2]));
}
//-----------------------------------------------------------------------------
template<typename FieldVecType,
typename ParametricCoordType>
VTKM_EXEC_EXPORT
typename FieldVecType::ComponentType
CellInterpolate(const FieldVecType &field,
const vtkm::Vec<ParametricCoordType,3> &pcoords,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase &vtkmNotUsed(worklet))
{
VTKM_ASSERT(field.GetNumberOfComponents() == 8);
typedef typename FieldVecType::ComponentType T;
T bottomFrontInterp = vtkm::Lerp(field[0], field[1], pcoords[0]);
T bottomBackInterp = vtkm::Lerp(field[3], field[2], pcoords[0]);
T topFrontInterp = vtkm::Lerp(field[4], field[5], pcoords[0]);
T topBackInterp = vtkm::Lerp(field[7], field[6], pcoords[0]);
T bottomInterp = vtkm::Lerp(bottomFrontInterp, bottomBackInterp, pcoords[1]);
T topInterp = vtkm::Lerp(topFrontInterp, topBackInterp, pcoords[1]);
return vtkm::Lerp(bottomInterp, topInterp, pcoords[2]);
}
template<typename ParametricCoordType>
VTKM_EXEC_EXPORT
vtkm::Vec<vtkm::FloatDefault,3>
CellInterpolate(const vtkm::VecRectilinearPointCoordinates<3> &field,
const vtkm::Vec<ParametricCoordType,3> &pcoords,
vtkm::CellShapeTagHexahedron,
const vtkm::exec::FunctorBase &)
{
vtkm::Vec<vtkm::FloatDefault,3> pcoordsCast(
static_cast<vtkm::FloatDefault>(pcoords[0]),
static_cast<vtkm::FloatDefault>(pcoords[1]),
static_cast<vtkm::FloatDefault>(pcoords[2]));
return field.GetOrigin() + pcoordsCast*field.GetSpacing();
}
//-----------------------------------------------------------------------------
template<typename FieldVecType,
typename ParametricCoordType>
VTKM_EXEC_EXPORT
typename FieldVecType::ComponentType
CellInterpolate(const FieldVecType &field,
const vtkm::Vec<ParametricCoordType,3> &pcoords,
vtkm::CellShapeTagWedge,
const vtkm::exec::FunctorBase &vtkmNotUsed(worklet))
{
VTKM_ASSERT(field.GetNumberOfComponents() == 6);
typedef typename FieldVecType::ComponentType T;
T bottomInterp = static_cast<T>( (field[0] * (1 - pcoords[0] - pcoords[1]))
+ (field[1] * pcoords[1])
+ (field[2] * pcoords[0]));
T topInterp = static_cast<T>( (field[3] * (1 - pcoords[0] - pcoords[1]))
+ (field[4] * pcoords[1])
+ (field[5] * pcoords[0]));
return vtkm::Lerp(bottomInterp, topInterp, pcoords[2]);
}
//-----------------------------------------------------------------------------
template<typename FieldVecType,
typename ParametricCoordType>
VTKM_EXEC_EXPORT
typename FieldVecType::ComponentType
CellInterpolate(const FieldVecType &field,
const vtkm::Vec<ParametricCoordType,3> &pcoords,
vtkm::CellShapeTagPyramid,
const vtkm::exec::FunctorBase &vtkmNotUsed(worklet))
{
VTKM_ASSERT(field.GetNumberOfComponents() == 5);
typedef typename FieldVecType::ComponentType T;
T frontInterp = vtkm::Lerp(field[0], field[1], pcoords[0]);
T backInterp = vtkm::Lerp(field[3], field[2], pcoords[0]);
T baseInterp = vtkm::Lerp(frontInterp, backInterp, pcoords[1]);
return vtkm::Lerp(baseInterp, field[4], pcoords[2]);
}
}
} // namespace vtkm::exec
#endif //vtk_m_exec_Interpolate_h
Reference in New Issue View Git Blame Copy Permalink