diff --git a/include/fcl/BV/BV_node.h b/include/fcl/BV/BV_node.h
index 6dd62d743..97ed69e10 100644
--- a/include/fcl/BV/BV_node.h
+++ b/include/fcl/BV/BV_node.h
@@ -156,7 +156,7 @@ Vector3<typename BVNode<BV>::Scalar> BVNode<BV>::getCenter() const
 template <typename Scalar, typename BV>
 struct GetOrientationImpl
 {
-  Matrix3<Scalar> operator()(/*const BVNode<BV>& node*/)
+  static Matrix3<Scalar> run(const BVNode<BV>& /*node*/)
   {
     return Matrix3<Scalar>::Identity();
   }
@@ -166,15 +166,14 @@ struct GetOrientationImpl
 template <typename BV>
 Matrix3<typename BV::Scalar> BVNode<BV>::getOrientation() const
 {
-  GetOrientationImpl<typename BV::Scalar, BV> getOrientationImpl;
-  return getOrientationImpl(bv);
+  return GetOrientationImpl<typename BV::Scalar, BV>::run(bv);
 }
 
 //==============================================================================
 template <typename Scalar>
 struct GetOrientationImpl<Scalar, OBB<Scalar>>
 {
-  Matrix3<Scalar> operator()(const OBB<Scalar>& bv)
+  static Matrix3<Scalar> run(const OBB<Scalar>& bv)
   {
     return bv.frame.linear();
   }
@@ -184,7 +183,7 @@ struct GetOrientationImpl<Scalar, OBB<Scalar>>
 template <typename Scalar>
 struct GetOrientationImpl<Scalar, RSS<Scalar>>
 {
-  Matrix3<Scalar> operator()(const RSS<Scalar>& bv)
+  static Matrix3<Scalar> run(const RSS<Scalar>& bv)
   {
     return bv.frame.linear();
   }
@@ -194,7 +193,7 @@ struct GetOrientationImpl<Scalar, RSS<Scalar>>
 template <typename Scalar>
 struct GetOrientationImpl<Scalar, OBBRSS<Scalar>>
 {
-  Matrix3<Scalar> operator()(const OBBRSS<Scalar>& bv)
+  static Matrix3<Scalar> run(const OBBRSS<Scalar>& bv)
   {
     return bv.obb.frame.linear();
   }
diff --git a/include/fcl/BV/kDOP.h b/include/fcl/BV/kDOP.h
index db3556e3b..fd8139039 100644
--- a/include/fcl/BV/kDOP.h
+++ b/include/fcl/BV/kDOP.h
@@ -412,7 +412,7 @@ void minmax(Scalar p, Scalar& minv, Scalar& maxv)
 template <typename Scalar, std::size_t N>
 struct GetDistancesImpl
 {
-  void operator()(const Vector3<Scalar>& p, Scalar* d)
+  static void run(const Vector3<Scalar>& /*p*/, Scalar* /*d*/)
   {
     // Do nothing
   }
@@ -422,15 +422,14 @@ struct GetDistancesImpl
 template <typename Scalar, std::size_t N>
 void getDistances(const Vector3<Scalar>& p, Scalar* d)
 {
-  GetDistancesImpl<Scalar, N> getDistancesImpl;
-  getDistancesImpl(p, d);
+  GetDistancesImpl<Scalar, N>::run(p, d);
 }
 
 //==============================================================================
 template <typename Scalar>
 struct GetDistancesImpl<Scalar, 5>
 {
-  void operator()(const Vector3<Scalar>& p, Scalar* d)
+  static void run(const Vector3<Scalar>& p, Scalar* d)
   {
     d[0] = p[0] + p[1];
     d[1] = p[0] + p[2];
@@ -444,7 +443,7 @@ struct GetDistancesImpl<Scalar, 5>
 template <typename Scalar>
 struct GetDistancesImpl<Scalar, 6>
 {
-  void operator()(const Vector3<Scalar>& p, Scalar* d)
+  static void run(const Vector3<Scalar>& p, Scalar* d)
   {
     d[0] = p[0] + p[1];
     d[1] = p[0] + p[2];
@@ -459,7 +458,7 @@ struct GetDistancesImpl<Scalar, 6>
 template <typename Scalar>
 struct GetDistancesImpl<Scalar, 9>
 {
-  void operator()(const Vector3<Scalar>& p, Scalar* d)
+  static void run(const Vector3<Scalar>& p, Scalar* d)
   {
     d[0] = p[0] + p[1];
     d[1] = p[0] + p[2];
diff --git a/include/fcl/BVH/BVH_model.h b/include/fcl/BVH/BVH_model.h
index 30f214805..4876ff2f6 100644
--- a/include/fcl/BVH/BVH_model.h
+++ b/include/fcl/BVH/BVH_model.h
@@ -378,7 +378,7 @@ OBJECT_TYPE BVHModel<BV>::getObjectType() const
 template <typename BV>
 struct GetNodeTypeImpl
 {
-  NODE_TYPE operator()()
+  static NODE_TYPE run()
   {
     return BV_UNKNOWN;
   }
@@ -388,8 +388,7 @@ struct GetNodeTypeImpl
 template <typename BV>
 NODE_TYPE BVHModel<BV>::getNodeType() const
 {
-  GetNodeTypeImpl<BV> getNodeTypeImpl;
-  return getNodeTypeImpl();
+  return GetNodeTypeImpl<BV>::run();
 }
 
 //==============================================================================
@@ -1205,7 +1204,7 @@ int BVHModel<BV>::recursiveRefitTree_bottomup(int bv_id)
 template <typename Scalar, typename BV>
 struct MakeParentRelativeRecurseImpl
 {
-  void operator()(BVHModel<BV>& model,
+  static void run(BVHModel<BV>& model,
                   int bv_id,
                   const Matrix3<Scalar>& parent_axis,
                   const Vector3<Scalar>& parent_c)
@@ -1230,8 +1229,8 @@ void BVHModel<BV>::makeParentRelativeRecurse(
     const Matrix3<Scalar>& parent_axis,
     const Vector3<Scalar>& parent_c)
 {
-  MakeParentRelativeRecurseImpl<typename BV::Scalar, BV> tmp;
-  tmp(*this, bv_id, parent_axis, parent_c);
+  MakeParentRelativeRecurseImpl<typename BV::Scalar, BV>::run(
+        *this, bv_id, parent_axis, parent_c);
 }
 
 //==============================================================================
@@ -1278,7 +1277,7 @@ void BVHModel<BV>::computeLocalAABB()
 template <typename Scalar>
 struct MakeParentRelativeRecurseImpl<Scalar, OBB<Scalar>>
 {
-  void operator()(BVHModel<OBB<Scalar>>& model,
+  static void run(BVHModel<OBB<Scalar>>& model,
                   int bv_id,
                   const Matrix3<Scalar>& parent_axis,
                   const Vector3<Scalar>& parent_c)
@@ -1303,7 +1302,7 @@ struct MakeParentRelativeRecurseImpl<Scalar, OBB<Scalar>>
 template <typename Scalar>
 struct MakeParentRelativeRecurseImpl<Scalar, RSS<Scalar>>
 {
-  void operator()(BVHModel<RSS<Scalar>>& model,
+  static void run(BVHModel<RSS<Scalar>>& model,
                   int bv_id,
                   const Matrix3<Scalar>& parent_axis,
                   const Vector3<Scalar>& parent_c)
@@ -1328,7 +1327,7 @@ struct MakeParentRelativeRecurseImpl<Scalar, RSS<Scalar>>
 template <typename Scalar>
 struct MakeParentRelativeRecurseImpl<Scalar, OBBRSS<Scalar>>
 {
-  void operator()(BVHModel<OBBRSS<Scalar>>& model,
+  static void run(BVHModel<OBBRSS<Scalar>>& model,
                   int bv_id,
                   const Matrix3<Scalar>& parent_axis,
                   const Vector3<Scalar>& parent_c)
@@ -1357,7 +1356,7 @@ struct MakeParentRelativeRecurseImpl<Scalar, OBBRSS<Scalar>>
 template <typename Scalar>
 struct GetNodeTypeImpl<AABB<Scalar>>
 {
-  NODE_TYPE operator()()
+  static NODE_TYPE run()
   {
     return BV_AABB;
   }
@@ -1367,7 +1366,7 @@ struct GetNodeTypeImpl<AABB<Scalar>>
 template <typename Scalar>
 struct GetNodeTypeImpl<OBB<Scalar>>
 {
-  NODE_TYPE operator()()
+  static NODE_TYPE run()
   {
     return BV_OBB;
   }
@@ -1377,7 +1376,7 @@ struct GetNodeTypeImpl<OBB<Scalar>>
 template <typename Scalar>
 struct GetNodeTypeImpl<RSS<Scalar>>
 {
-  NODE_TYPE operator()()
+  static NODE_TYPE run()
   {
     return BV_RSS;
   }
@@ -1387,7 +1386,7 @@ struct GetNodeTypeImpl<RSS<Scalar>>
 template <typename Scalar>
 struct GetNodeTypeImpl<kIOS<Scalar>>
 {
-  NODE_TYPE operator()()
+  static NODE_TYPE run()
   {
     return BV_kIOS;
   }
@@ -1397,7 +1396,7 @@ struct GetNodeTypeImpl<kIOS<Scalar>>
 template <typename Scalar>
 struct GetNodeTypeImpl<OBBRSS<Scalar>>
 {
-  NODE_TYPE operator()()
+  static NODE_TYPE run()
   {
     return BV_OBBRSS;
   }
@@ -1407,7 +1406,7 @@ struct GetNodeTypeImpl<OBBRSS<Scalar>>
 template <typename Scalar>
 struct GetNodeTypeImpl<KDOP<Scalar, 16>>
 {
-  NODE_TYPE operator()()
+  static NODE_TYPE run()
   {
     return BV_KDOP16;
   }
@@ -1417,7 +1416,7 @@ struct GetNodeTypeImpl<KDOP<Scalar, 16>>
 template <typename Scalar>
 struct GetNodeTypeImpl<KDOP<Scalar, 18>>
 {
-  NODE_TYPE operator()()
+  static NODE_TYPE run()
   {
     return BV_KDOP18;
   }
@@ -1427,7 +1426,7 @@ struct GetNodeTypeImpl<KDOP<Scalar, 18>>
 template <typename Scalar>
 struct GetNodeTypeImpl<KDOP<Scalar, 24>>
 {
-  NODE_TYPE operator()()
+  static NODE_TYPE run()
   {
     return BV_KDOP24;
   }
diff --git a/include/fcl/BVH/BV_splitter.h b/include/fcl/BVH/BV_splitter.h
index 21256f561..994af8f44 100644
--- a/include/fcl/BVH/BV_splitter.h
+++ b/include/fcl/BVH/BV_splitter.h
@@ -246,7 +246,7 @@ void BVSplitter<BV>::computeRule(
 template <typename Scalar, typename BV>
 struct ApplyImpl
 {
-  bool operator()(
+  static bool run(
       const BVSplitter<BV>& splitter, const Vector3<Scalar>& q)
   {
     return q[splitter.split_axis] > splitter.split_value;
@@ -257,15 +257,14 @@ struct ApplyImpl
 template <typename BV>
 bool BVSplitter<BV>::apply(const Vector3<Scalar>& q) const
 {
-  ApplyImpl<Scalar, BV> applyImpl;
-  return applyImpl(*this, q);
+  return ApplyImpl<Scalar, BV>::run(*this, q);
 }
 
 //==============================================================================
 template <typename Scalar, typename BV>
 struct ComputeRuleCenterImpl
 {
-  void operator()(
+  static void run(
       BVSplitter<BV>& splitter,
       const BV& bv,
       unsigned int* /*primitive_indices*/,
@@ -289,15 +288,15 @@ template <typename BV>
 void BVSplitter<BV>::computeRule_bvcenter(
     const BV& bv, unsigned int* primitive_indices, int num_primitives)
 {
-  ComputeRuleCenterImpl<Scalar, BV> computeRuleCenterImpl;
-  computeRuleCenterImpl(*this, bv, primitive_indices, num_primitives);
+  ComputeRuleCenterImpl<Scalar, BV>::run(
+        *this, bv, primitive_indices, num_primitives);
 }
 
 //==============================================================================
 template <typename Scalar, typename BV>
 struct ComputeRuleMeanImpl
 {
-  void operator()(
+  static void run(
       BVSplitter<BV>& splitter,
       const BV& bv,
       unsigned int* primitive_indices,
@@ -342,15 +341,15 @@ template <typename BV>
 void BVSplitter<BV>::computeRule_mean(
     const BV& bv, unsigned int* primitive_indices, int num_primitives)
 {
-  ComputeRuleMeanImpl<Scalar, BV> computeRuleMeanImpl;
-  computeRuleMeanImpl(*this, bv, primitive_indices, num_primitives);
+  ComputeRuleMeanImpl<Scalar, BV>::run(
+        *this, bv, primitive_indices, num_primitives);
 }
 
 //==============================================================================
 template <typename Scalar, typename BV>
 struct ComputeRuleMedianImpl
 {
-  void operator()(
+  static void run(
       BVSplitter<BV>& splitter,
       const BV& bv,
       unsigned int* primitive_indices,
@@ -400,15 +399,15 @@ template <typename BV>
 void BVSplitter<BV>::computeRule_median(
     const BV& bv, unsigned int* primitive_indices, int num_primitives)
 {
-  ComputeRuleMedianImpl<Scalar, BV> computeRuleMedianImpl;
-  computeRuleMedianImpl(*this, bv, primitive_indices, num_primitives);
+  ComputeRuleMedianImpl<Scalar, BV>::run(
+        *this, bv, primitive_indices, num_primitives);
 }
 
 //==============================================================================
 template <typename Scalar>
 struct ComputeRuleCenterImpl<Scalar, OBB<Scalar>>
 {
-  void operator()(
+  static void run(
       BVSplitter<OBB<Scalar>>& splitter,
       const OBB<Scalar>& bv,
       unsigned int* /*primitive_indices*/,
@@ -423,7 +422,7 @@ struct ComputeRuleCenterImpl<Scalar, OBB<Scalar>>
 template <typename Scalar>
 struct ComputeRuleMeanImpl<Scalar, OBB<Scalar>>
 {
-  void operator()(
+  static void run(
       BVSplitter<OBB<Scalar>>& splitter,
       const OBB<Scalar>& bv,
       unsigned int* primitive_indices,
@@ -440,7 +439,7 @@ struct ComputeRuleMeanImpl<Scalar, OBB<Scalar>>
 template <typename Scalar>
 struct ComputeRuleMedianImpl<Scalar, OBB<Scalar>>
 {
-  void operator()(
+  static void run(
       BVSplitter<OBB<Scalar>>& splitter,
       const OBB<Scalar>& bv,
       unsigned int* primitive_indices,
@@ -457,7 +456,7 @@ struct ComputeRuleMedianImpl<Scalar, OBB<Scalar>>
 template <typename Scalar>
 struct ComputeRuleCenterImpl<Scalar, RSS<Scalar>>
 {
-  void operator()(
+  static void run(
       BVSplitter<RSS<Scalar>>& splitter,
       const RSS<Scalar>& bv,
       unsigned int* /*primitive_indices*/,
@@ -472,7 +471,7 @@ struct ComputeRuleCenterImpl<Scalar, RSS<Scalar>>
 template <typename Scalar>
 struct ComputeRuleMeanImpl<Scalar, RSS<Scalar>>
 {
-  void operator()(
+  static void run(
       BVSplitter<RSS<Scalar>>& splitter,
       const RSS<Scalar>& bv,
       unsigned int* primitive_indices,
@@ -489,7 +488,7 @@ struct ComputeRuleMeanImpl<Scalar, RSS<Scalar>>
 template <typename Scalar>
 struct ComputeRuleMedianImpl<Scalar, RSS<Scalar>>
 {
-  void operator()(
+  static void run(
       BVSplitter<RSS<Scalar>>& splitter,
       const RSS<Scalar>& bv,
       unsigned int* primitive_indices,
@@ -506,7 +505,7 @@ struct ComputeRuleMedianImpl<Scalar, RSS<Scalar>>
 template <typename Scalar>
 struct ComputeRuleCenterImpl<Scalar, kIOS<Scalar>>
 {
-  void operator()(
+  static void run(
       BVSplitter<kIOS<Scalar>>& splitter,
       const kIOS<Scalar>& bv,
       unsigned int* /*primitive_indices*/,
@@ -521,7 +520,7 @@ struct ComputeRuleCenterImpl<Scalar, kIOS<Scalar>>
 template <typename Scalar>
 struct ComputeRuleMeanImpl<Scalar, kIOS<Scalar>>
 {
-  void operator()(
+  static void run(
       BVSplitter<kIOS<Scalar>>& splitter,
       const kIOS<Scalar>& bv,
       unsigned int* primitive_indices,
@@ -538,7 +537,7 @@ struct ComputeRuleMeanImpl<Scalar, kIOS<Scalar>>
 template <typename Scalar>
 struct ComputeRuleMedianImpl<Scalar, kIOS<Scalar>>
 {
-  void operator()(
+  static void run(
       BVSplitter<kIOS<Scalar>>& splitter,
       const kIOS<Scalar>& bv,
       unsigned int* primitive_indices,
@@ -555,7 +554,7 @@ struct ComputeRuleMedianImpl<Scalar, kIOS<Scalar>>
 template <typename Scalar>
 struct ComputeRuleCenterImpl<Scalar, OBBRSS<Scalar>>
 {
-  void operator()(
+  static void run(
       BVSplitter<OBBRSS<Scalar>>& splitter,
       const OBBRSS<Scalar>& bv,
       unsigned int* /*primitive_indices*/,
@@ -570,7 +569,7 @@ struct ComputeRuleCenterImpl<Scalar, OBBRSS<Scalar>>
 template <typename Scalar>
 struct ComputeRuleMeanImpl<Scalar, OBBRSS<Scalar>>
 {
-  void operator()(
+  static void run(
       BVSplitter<OBBRSS<Scalar>>& splitter,
       const OBBRSS<Scalar>& bv,
       unsigned int* primitive_indices,
@@ -587,7 +586,7 @@ struct ComputeRuleMeanImpl<Scalar, OBBRSS<Scalar>>
 template <typename Scalar>
 struct ComputeRuleMedianImpl<Scalar, OBBRSS<Scalar>>
 {
-  void operator()(
+  static void run(
       BVSplitter<OBBRSS<Scalar>>& splitter,
       const OBBRSS<Scalar>& bv,
       unsigned int* primitive_indices,
@@ -604,7 +603,7 @@ struct ComputeRuleMedianImpl<Scalar, OBBRSS<Scalar>>
 template <typename Scalar>
 struct ApplyImpl<Scalar, OBB<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const BVSplitter<OBB<Scalar>>& splitter,
       const Vector3<Scalar>& q)
   {
@@ -616,7 +615,7 @@ struct ApplyImpl<Scalar, OBB<Scalar>>
 template <typename Scalar>
 struct ApplyImpl<Scalar, RSS<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const BVSplitter<RSS<Scalar>>& splitter,
       const Vector3<Scalar>& q)
   {
@@ -628,7 +627,7 @@ struct ApplyImpl<Scalar, RSS<Scalar>>
 template <typename Scalar>
 struct ApplyImpl<Scalar, kIOS<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const BVSplitter<kIOS<Scalar>>& splitter,
       const Vector3<Scalar>& q)
   {
@@ -640,7 +639,7 @@ struct ApplyImpl<Scalar, kIOS<Scalar>>
 template <typename Scalar>
 struct ApplyImpl<Scalar, OBBRSS<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const BVSplitter<OBBRSS<Scalar>>& splitter,
       const Vector3<Scalar>& q)
   {
@@ -661,7 +660,7 @@ void BVSplitter<BV>::clear()
 template <typename Scalar, typename BV>
 struct ComputeSplitVectorImpl
 {
-  void operator()(const BV& bv, Vector3<Scalar>& split_vector)
+  static void run(const BV& bv, Vector3<Scalar>& split_vector)
   {
     split_vector = bv.frame.linear().col(0);
   }
@@ -671,15 +670,14 @@ struct ComputeSplitVectorImpl
 template <typename Scalar, typename BV>
 void computeSplitVector(const BV& bv, Vector3<Scalar>& split_vector)
 {
-  ComputeSplitVectorImpl<Scalar, BV> computeSplitVectorImpl;
-  computeSplitVectorImpl(bv, split_vector);
+  ComputeSplitVectorImpl<Scalar, BV>::run(bv, split_vector);
 }
 
 //==============================================================================
 template <typename Scalar>
 struct ComputeSplitVectorImpl<Scalar, kIOS<Scalar>>
 {
-  void operator()(const kIOS<Scalar>& bv, Vector3<Scalar>& split_vector)
+  static void run(const kIOS<Scalar>& bv, Vector3<Scalar>& split_vector)
   {
     /*
       switch(bv.num_spheres)
@@ -717,7 +715,7 @@ struct ComputeSplitVectorImpl<Scalar, kIOS<Scalar>>
 template <typename Scalar>
 struct ComputeSplitVectorImpl<Scalar, OBBRSS<Scalar>>
 {
-  void operator()(const OBBRSS<Scalar>& bv, Vector3<Scalar>& split_vector)
+  static void run(const OBBRSS<Scalar>& bv, Vector3<Scalar>& split_vector)
   {
     split_vector = bv.obb.frame.linear().col(0);
   }
diff --git a/include/fcl/broadphase/hierarchy_tree.h b/include/fcl/broadphase/hierarchy_tree.h
index df8af5720..7021ad7d4 100644
--- a/include/fcl/broadphase/hierarchy_tree.h
+++ b/include/fcl/broadphase/hierarchy_tree.h
@@ -94,15 +94,23 @@ bool nodeBaseLess(NodeBase<BV>* a, NodeBase<BV>* b, int d)
 
 //==============================================================================
 template <typename Scalar, typename BV>
-struct SelectImpl1
+struct SelectImpl
 {
-  std::size_t operator()(
+  static std::size_t run(
       const NodeBase<BV>& /*query*/,
       const NodeBase<BV>& /*node1*/,
       const NodeBase<BV>& /*node2*/)
   {
     return 0;
   }
+
+  static std::size_t run(
+      const BV& /*query*/,
+      const NodeBase<BV>& /*node1*/,
+      const NodeBase<BV>& /*node2*/)
+  {
+    return 0;
+  }
 };
 
 /// @brief select from node1 and node2 which is close to a given query. 0 for node1 and 1 for node2
@@ -112,23 +120,9 @@ size_t select(
     const NodeBase<BV>& node1,
     const NodeBase<BV>& node2)
 {
-  SelectImpl1<typename BV::Scalar, BV> selectImpl;
-  return selectImpl(query, node1, node2);
+  return SelectImpl<typename BV::Scalar, BV>::run(query, node1, node2);
 }
 
-//==============================================================================
-template <typename Scalar, typename BV>
-struct SelectImpl2
-{
-  std::size_t operator()(
-      const BV& /*query*/,
-      const NodeBase<BV>& /*node1*/,
-      const NodeBase<BV>& /*node2*/)
-  {
-    return 0;
-  }
-};
-
 /// @brief select from node1 and node2 which is close to a given query bounding volume. 0 for node1 and 1 for node2
 template<typename BV>
 size_t select(
@@ -136,30 +130,18 @@ size_t select(
     const NodeBase<BV>& node1,
     const NodeBase<BV>& node2)
 {
-  SelectImpl2<typename BV::Scalar, BV> selectImpl;
-  return selectImpl(query, node1, node2);
+  return SelectImpl<typename BV::Scalar, BV>::run(query, node1, node2);
 }
 
 /// @brief Class for hierarchy tree structure
 template<typename BV>
 class HierarchyTree
 {
+public:
 
   using Scalar = typename BV::Scalar;
 
   typedef NodeBase<BV> NodeType;
-  typedef typename std::vector<NodeBase<BV>* >::iterator NodeVecIterator;
-  typedef typename std::vector<NodeBase<BV>* >::const_iterator NodeVecConstIterator;
-
-  struct SortByMorton
-  {
-    bool operator() (const NodeType* a, const NodeType* b) const
-    {
-      return a->code < b->code;
-    }
-  };
-
-public:
 
   /// @brief Create hierarchy tree with suitable setting.
   /// bu_threshold decides the height of tree node to start bottom-up construction / optimization;
@@ -230,6 +212,17 @@ class HierarchyTree
 
 private:
 
+  typedef typename std::vector<NodeBase<BV>* >::iterator NodeVecIterator;
+  typedef typename std::vector<NodeBase<BV>* >::const_iterator NodeVecConstIterator;
+
+  struct SortByMorton
+  {
+    bool operator() (const NodeType* a, const NodeType* b) const
+    {
+      return a->code < b->code;
+    }
+  };
+
   /// @brief construct a tree for a set of leaves from bottom -- very heavy way 
   void bottomup(const NodeVecIterator lbeg, const NodeVecIterator lend);
 
@@ -388,9 +381,14 @@ struct nodeBaseLess
 
 //==============================================================================
 template <typename Scalar, typename BV>
-struct SelectImpl3
+struct SelectImpl
 {
-  bool operator()(size_t query, size_t node1, size_t node2, NodeBase<BV>* nodes)
+  static bool run(size_t query, size_t node1, size_t node2, NodeBase<BV>* nodes)
+  {
+    return 0;
+  }
+
+  static bool run(const BV& query, size_t node1, size_t node2, NodeBase<BV>* nodes)
   {
     return 0;
   }
@@ -400,26 +398,14 @@ struct SelectImpl3
 template<typename BV> 
 size_t select(size_t query, size_t node1, size_t node2, NodeBase<BV>* nodes)
 {
-  SelectImpl3<typename BV::Scalar, BV> selectImpl;
-  return selectImpl(query, node1, node2, nodes);
+  return SelectImpl<typename BV::Scalar, BV>::run(query, node1, node2, nodes);
 }
 
-//==============================================================================
-template <typename Scalar, typename BV>
-struct SelectImpl4
-{
-  bool operator()(const BV& query, size_t node1, size_t node2, NodeBase<BV>* nodes)
-  {
-    return 0;
-  }
-};
-
 /// @brief select the node from node1 and node2 which is close to the query AABB<Scalar>. 0 for node1 and 1 for node2.
 template<typename BV>
 size_t select(const BV& query, size_t node1, size_t node2, NodeBase<BV>* nodes)
 {
-  SelectImpl4<typename BV::Scalar, BV> selectImpl;
-  return selectImpl(query, node1, node2, nodes);
+  return SelectImpl<typename BV::Scalar, BV>::run(query, node1, node2, nodes);
 }
 
 /// @brief Class for hierarchy tree structure
@@ -732,9 +718,9 @@ bool HierarchyTree<BV>::update(NodeType* leaf, const BV& bv)
 
 //==============================================================================
 template <typename Scalar, typename BV>
-struct UpdateImpl1
+struct UpdateImpl
 {
-  bool operator()(
+  static bool run(
       const HierarchyTree<BV>& tree,
       typename HierarchyTree<BV>::NodeType* leaf,
       const BV& bv,
@@ -745,36 +731,29 @@ struct UpdateImpl1
     tree.update_(leaf, bv);
     return true;
   }
-};
-
-template<typename BV>
-bool HierarchyTree<BV>::update(NodeType* leaf, const BV& bv, const Vector3<Scalar>& vel, Scalar margin)
-{
-  UpdateImpl1<typename BV::Scalar, BV> updateImpl;
-  return updateImpl(leaf, bv, vel, margin);
-}
 
-//==============================================================================
-template <typename Scalar, typename BV>
-struct UpdateImpl2
-{
-  bool operator()(
+  static bool run(
       const HierarchyTree<BV>& tree,
       typename HierarchyTree<BV>::NodeType* leaf,
       const BV& bv,
-      const Vector3<Scalar>& vel)
+      const Vector3<Scalar>& /*vel*/)
   {
     if(leaf->bv.contain(bv)) return false;
-    update_(leaf, bv);
+    tree.update_(leaf, bv);
     return true;
   }
 };
 
+template<typename BV>
+bool HierarchyTree<BV>::update(NodeType* leaf, const BV& bv, const Vector3<Scalar>& vel, Scalar margin)
+{
+  return UpdateImpl<typename BV::Scalar, BV>::run(*this, leaf, bv, vel, margin);
+}
+
 template<typename BV>
 bool HierarchyTree<BV>::update(NodeType* leaf, const BV& bv, const Vector3<Scalar>& vel)
 {
-  UpdateImpl2<typename BV::Scalar, BV> updateImpl;
-  return updateImpl(leaf, bv, vel);
+  return UpdateImpl<typename BV::Scalar, BV>::run(*this, leaf, bv, vel);
 }
 
 template<typename BV>
@@ -2504,9 +2483,9 @@ void HierarchyTree<BV>::fetchLeaves(size_t root, NodeType*& leaves, int depth)
 
 //==============================================================================
 template <typename Scalar>
-struct SelectImpl1<Scalar, AABB<Scalar>>
+struct SelectImpl<Scalar, AABB<Scalar>>
 {
-  std::size_t operator()(
+  static std::size_t run(
       const NodeBase<AABB<Scalar>>& node,
       const NodeBase<AABB<Scalar>>& node1,
       const NodeBase<AABB<Scalar>>& node2)
@@ -2521,13 +2500,8 @@ struct SelectImpl1<Scalar, AABB<Scalar>>
     Scalar d2 = fabs(v2[0]) + fabs(v2[1]) + fabs(v2[2]);
     return (d1 < d2) ? 0 : 1;
   }
-};
 
-//==============================================================================
-template <typename Scalar>
-struct SelectImpl2<Scalar, AABB<Scalar>>
-{
-  std::size_t operator()(
+  static std::size_t run(
       const AABB<Scalar>& query,
       const NodeBase<AABB<Scalar>>& node1,
       const NodeBase<AABB<Scalar>>& node2)
@@ -2542,13 +2516,8 @@ struct SelectImpl2<Scalar, AABB<Scalar>>
     Scalar d2 = fabs(v2[0]) + fabs(v2[1]) + fabs(v2[2]);
     return (d1 < d2) ? 0 : 1;
   }
-};
 
-//==============================================================================
-template <typename Scalar>
-struct UpdateImpl1<Scalar, AABB<Scalar>>
-{
-  bool operator()(
+  static bool run(
       const HierarchyTree<AABB<Scalar>>& tree,
       typename HierarchyTree<AABB<Scalar>>::NodeType* leaf,
       const AABB<Scalar>& bv_,
@@ -2569,13 +2538,8 @@ struct UpdateImpl1<Scalar, AABB<Scalar>>
     tree.update(leaf, bv);
     return true;
   }
-};
 
-//==============================================================================
-template <typename Scalar>
-struct UpdateImpl2<Scalar, AABB<Scalar>>
-{
-  bool operator()(
+  static bool run(
       const HierarchyTree<AABB<Scalar>>& tree,
       typename HierarchyTree<AABB<Scalar>>::NodeType* leaf,
       const AABB<Scalar>& bv_,
@@ -2599,9 +2563,9 @@ namespace implementation_array
 
 //==============================================================================
 template <typename Scalar>
-struct SelectImpl3<Scalar, AABB<Scalar>>
+struct SelectImpl<Scalar, AABB<Scalar>>
 {
-  bool operator()(size_t query, size_t node1, size_t node2, NodeBase<AABB<Scalar>>* nodes)
+  static bool run(size_t query, size_t node1, size_t node2, NodeBase<AABB<Scalar>>* nodes)
   {
     const AABB<Scalar>& bv = nodes[query].bv;
     const AABB<Scalar>& bv1 = nodes[node1].bv;
@@ -2613,13 +2577,8 @@ struct SelectImpl3<Scalar, AABB<Scalar>>
     Scalar d2 = fabs(v2[0]) + fabs(v2[1]) + fabs(v2[2]);
     return (d1 < d2) ? 0 : 1;
   }
-};
 
-//==============================================================================
-template <typename Scalar>
-struct SelectImpl4<Scalar, AABB<Scalar>>
-{
-  bool operator()(const AABB<Scalar>& query, size_t node1, size_t node2, NodeBase<AABB<Scalar>>* nodes)
+  static bool run(const AABB<Scalar>& query, size_t node1, size_t node2, NodeBase<AABB<Scalar>>* nodes)
   {
     const AABB<Scalar>& bv = query;
     const AABB<Scalar>& bv1 = nodes[node1].bv;
diff --git a/include/fcl/ccd/conservative_advancement.h b/include/fcl/ccd/conservative_advancement.h
index 2a909780a..dd4001269 100644
--- a/include/fcl/ccd/conservative_advancement.h
+++ b/include/fcl/ccd/conservative_advancement.h
@@ -610,9 +610,9 @@ bool conservativeAdvancementShapeMeshOriented(const S& o1,
 
 //==============================================================================
 template <typename Scalar, typename S, typename NarrowPhaseSolver>
-struct ConservativeAdvancementImpl1
+struct ConservativeAdvancementImpl
 {
-  bool operator()(
+  static bool run(
       const S& o1,
       const MotionBase<Scalar>* motion1,
       const BVHModel<RSS<Scalar>>& o2,
@@ -620,17 +620,12 @@ struct ConservativeAdvancementImpl1
       const NarrowPhaseSolver* nsolver,
       const CollisionRequest<Scalar>& request,
       CollisionResult<Scalar>& result,
-      typename NarrowPhaseSolver::Scalar& toc)
+      Scalar& toc)
   {
     return details::conservativeAdvancementShapeMeshOriented<S, RSS<Scalar>, NarrowPhaseSolver, ShapeMeshConservativeAdvancementTraversalNodeRSS<S, NarrowPhaseSolver> >(o1, motion1, o2, motion2, nsolver, request, result, toc);
   }
-};
 
-//==============================================================================
-template <typename Scalar, typename S, typename NarrowPhaseSolver>
-struct ConservativeAdvancementImpl2
-{
-  bool operator()(
+  static bool run(
       const S& o1,
       const MotionBase<Scalar>* motion1,
       const BVHModel<OBBRSS<Scalar>>& o2,
@@ -654,8 +649,9 @@ bool conservativeAdvancement(const S& o1,
                              CollisionResult<typename NarrowPhaseSolver::Scalar>& result,
                              typename NarrowPhaseSolver::Scalar& toc)
 {
-  ConservativeAdvancementImpl1<typename NarrowPhaseSolver::Scalar, S, NarrowPhaseSolver> conservativeAdvancementImpl;
-  return conservativeAdvancementImpl(o1, motion1, o2, motion2, nsolver, request, result, toc);
+  return ConservativeAdvancementImpl<
+      typename NarrowPhaseSolver::Scalar, S, NarrowPhaseSolver>::run(
+        o1, motion1, o2, motion2, nsolver, request, result, toc);
 }
 
 template<typename S, typename NarrowPhaseSolver>
@@ -668,15 +664,16 @@ bool conservativeAdvancement(const S& o1,
                              CollisionResult<typename NarrowPhaseSolver::Scalar>& result,
                              typename NarrowPhaseSolver::Scalar& toc)
 {
-  ConservativeAdvancementImpl2<typename NarrowPhaseSolver::Scalar, S, NarrowPhaseSolver> conservativeAdvancementImpl;
-  return conservativeAdvancementImpl(o1, motion1, o2, motion2, nsolver, request, result, toc);
+  return ConservativeAdvancementImpl<
+      typename NarrowPhaseSolver::Scalar, S, NarrowPhaseSolver>::run(
+        o1, motion1, o2, motion2, nsolver, request, result, toc);
 }
 
 //==============================================================================
 template <typename Scalar, typename NarrowPhaseSolver>
-struct ConservativeAdvancementImpl1<Scalar, BVHModel<RSS<Scalar>>, NarrowPhaseSolver>
+struct ConservativeAdvancementImpl<Scalar, BVHModel<RSS<Scalar>>, NarrowPhaseSolver>
 {
-  bool operator()(
+  static bool run(
       const BVHModel<RSS<Scalar>>& o1,
       const MotionBase<Scalar>* motion1,
       const BVHModel<RSS<Scalar>>& o2,
@@ -692,9 +689,9 @@ struct ConservativeAdvancementImpl1<Scalar, BVHModel<RSS<Scalar>>, NarrowPhaseSo
 
 //==============================================================================
 template <typename Scalar, typename NarrowPhaseSolver>
-struct ConservativeAdvancementImpl2<Scalar, BVHModel<OBBRSS<Scalar>>, NarrowPhaseSolver>
+struct ConservativeAdvancementImpl<Scalar, BVHModel<OBBRSS<Scalar>>, NarrowPhaseSolver>
 {
-  bool operator()(
+  static bool run(
       const BVHModel<OBBRSS<Scalar>>& o1,
       const MotionBase<Scalar>* motion1,
       const BVHModel<OBBRSS<Scalar>>& o2,
diff --git a/include/fcl/ccd/motion_base.h b/include/fcl/ccd/motion_base.h
index 7a016f2f2..ae2460bdb 100644
--- a/include/fcl/ccd/motion_base.h
+++ b/include/fcl/ccd/motion_base.h
@@ -203,7 +203,7 @@ using MotionBasePtr = std::shared_ptr<MotionBase<Scalar>>;
 template <typename Scalar, typename BV, typename MotionT>
 struct TBVMotionBoundVisitorVisitImpl
 {
-  Scalar operator()(
+  static Scalar run(
       const TBVMotionBoundVisitor<BV>& /*visitor*/,
       const MotionT& /*motion*/)
   {
@@ -224,8 +224,9 @@ typename BV::Scalar TBVMotionBoundVisitor<BV>::visit(
     const SplineMotion<Scalar>& motion) const
 {
   using Scalar = typename BV::Scalar;
-  TBVMotionBoundVisitorVisitImpl<Scalar, BV, SplineMotion<Scalar>> visitMotionBaseImpl;
-  return visitMotionBaseImpl(*this, motion);
+
+  return TBVMotionBoundVisitorVisitImpl<
+      Scalar, BV, SplineMotion<Scalar>>::run(*this, motion);
 }
 
 //==============================================================================
@@ -234,8 +235,9 @@ typename BV::Scalar TBVMotionBoundVisitor<BV>::visit(
     const ScrewMotion<Scalar>& motion) const
 {
   using Scalar = typename BV::Scalar;
-  TBVMotionBoundVisitorVisitImpl<Scalar, BV, ScrewMotion<Scalar>> visitMotionBaseImpl;
-  return visitMotionBaseImpl(*this, motion);
+
+  return TBVMotionBoundVisitorVisitImpl<
+      Scalar, BV, ScrewMotion<Scalar>>::run(*this, motion);
 }
 
 //==============================================================================
@@ -244,8 +246,9 @@ typename BV::Scalar TBVMotionBoundVisitor<BV>::visit(
     const InterpMotion<Scalar>& motion) const
 {
   using Scalar = typename BV::Scalar;
-  TBVMotionBoundVisitorVisitImpl<Scalar, BV, InterpMotion<Scalar>> visitMotionBaseImpl;
-  return visitMotionBaseImpl(*this, motion);
+
+  return TBVMotionBoundVisitorVisitImpl<
+      Scalar, BV, InterpMotion<Scalar>>::run(*this, motion);
 }
 
 //==============================================================================
@@ -254,15 +257,16 @@ typename BV::Scalar TBVMotionBoundVisitor<BV>::visit(
     const TranslationMotion<Scalar>& motion) const
 {
   using Scalar = typename BV::Scalar;
-  TBVMotionBoundVisitorVisitImpl<Scalar, BV, TranslationMotion<Scalar>> visitMotionBaseImpl;
-  return visitMotionBaseImpl(*this, motion);
+
+  return TBVMotionBoundVisitorVisitImpl<
+      Scalar, BV, TranslationMotion<Scalar>>::run(*this, motion);
 }
 
 //==============================================================================
 template <typename Scalar>
 struct TBVMotionBoundVisitorVisitImpl<Scalar, RSS<Scalar>, SplineMotion<Scalar>>
 {
-  Scalar operator()(
+  static Scalar run(
       const TBVMotionBoundVisitor<RSS<Scalar>>& visitor,
       const SplineMotion<Scalar>& motion)
   {
@@ -324,7 +328,7 @@ struct TBVMotionBoundVisitorVisitImpl<Scalar, RSS<Scalar>, SplineMotion<Scalar>>
 template <typename Scalar>
 struct TBVMotionBoundVisitorVisitImpl<Scalar, RSS<Scalar>, ScrewMotion<Scalar>>
 {
-  Scalar operator()(
+  static Scalar run(
       const TBVMotionBoundVisitor<RSS<Scalar>>& visitor,
       const ScrewMotion<Scalar>& motion)
   {
@@ -365,7 +369,7 @@ struct TBVMotionBoundVisitorVisitImpl<Scalar, RSS<Scalar>, ScrewMotion<Scalar>>
 template <typename Scalar>
 struct TBVMotionBoundVisitorVisitImpl<Scalar, RSS<Scalar>, InterpMotion<Scalar>>
 {
-  Scalar operator()(
+  static Scalar run(
       const TBVMotionBoundVisitor<RSS<Scalar>>& visitor,
       const InterpMotion<Scalar>& motion)
   {
@@ -401,7 +405,7 @@ struct TBVMotionBoundVisitorVisitImpl<Scalar, RSS<Scalar>, InterpMotion<Scalar>>
 template <typename Scalar>
 struct TBVMotionBoundVisitorVisitImpl<Scalar, RSS<Scalar>, TranslationMotion<Scalar>>
 {
-  Scalar operator()(
+  static Scalar run(
       const TBVMotionBoundVisitor<RSS<Scalar>>& visitor,
       const TranslationMotion<Scalar>& motion)
   {
@@ -413,7 +417,7 @@ struct TBVMotionBoundVisitorVisitImpl<Scalar, RSS<Scalar>, TranslationMotion<Sca
 template <typename Scalar, typename MotionT>
 struct TriangleMotionBoundVisitorVisitImpl
 {
-  Scalar operator()(
+  static Scalar run(
       const TriangleMotionBoundVisitor<Scalar>& /*visitor*/,
       const MotionT& /*motion*/)
   {
@@ -426,8 +430,8 @@ template<typename Scalar>
 Scalar TriangleMotionBoundVisitor<Scalar>::visit(
     const SplineMotion<Scalar>& motion) const
 {
-  TriangleMotionBoundVisitorVisitImpl<Scalar, SplineMotion<Scalar>> visitMotionBaseImpl;
-  return visitMotionBaseImpl(*this, motion);
+  return TriangleMotionBoundVisitorVisitImpl<
+      Scalar, SplineMotion<Scalar>>::run(*this, motion);
 }
 
 //==============================================================================
@@ -435,8 +439,8 @@ template<typename Scalar>
 Scalar TriangleMotionBoundVisitor<Scalar>::visit(
     const ScrewMotion<Scalar>& motion) const
 {
-  TriangleMotionBoundVisitorVisitImpl<Scalar, ScrewMotion<Scalar>> visitMotionBaseImpl;
-  return visitMotionBaseImpl(*this, motion);
+  return TriangleMotionBoundVisitorVisitImpl<
+      Scalar, ScrewMotion<Scalar>>::run(*this, motion);
 }
 
 //==============================================================================
@@ -444,8 +448,8 @@ template<typename Scalar>
 Scalar TriangleMotionBoundVisitor<Scalar>::visit(
     const InterpMotion<Scalar>& motion) const
 {
-  TriangleMotionBoundVisitorVisitImpl<Scalar, InterpMotion<Scalar>> visitMotionBaseImpl;
-  return visitMotionBaseImpl(*this, motion);
+  return TriangleMotionBoundVisitorVisitImpl<
+      Scalar, InterpMotion<Scalar>>::run(*this, motion);
 }
 
 //==============================================================================
@@ -453,8 +457,8 @@ template<typename Scalar>
 Scalar TriangleMotionBoundVisitor<Scalar>::visit(
     const TranslationMotion<Scalar>& motion) const
 {
-  TriangleMotionBoundVisitorVisitImpl<Scalar, TranslationMotion<Scalar>> visitMotionBaseImpl;
-  return visitMotionBaseImpl(*this, motion);
+  return TriangleMotionBoundVisitorVisitImpl<
+      Scalar, TranslationMotion<Scalar>>::run(*this, motion);
 }
 
 //==============================================================================
@@ -465,7 +469,7 @@ Scalar TriangleMotionBoundVisitor<Scalar>::visit(
 template <typename Scalar>
 struct TriangleMotionBoundVisitorVisitImpl<Scalar, ScrewMotion<Scalar>>
 {
-  Scalar operator()(
+  static Scalar run(
       const TriangleMotionBoundVisitor<Scalar>& visitor,
       const ScrewMotion<Scalar>& motion)
   {
@@ -502,7 +506,7 @@ struct TriangleMotionBoundVisitorVisitImpl<Scalar, ScrewMotion<Scalar>>
 template <typename Scalar>
 struct TriangleMotionBoundVisitorVisitImpl<Scalar, InterpMotion<Scalar>>
 {
-  Scalar operator()(
+  static Scalar run(
       const TriangleMotionBoundVisitor<Scalar>& visitor,
       const InterpMotion<Scalar>& motion)
   {
@@ -535,7 +539,7 @@ struct TriangleMotionBoundVisitorVisitImpl<Scalar, InterpMotion<Scalar>>
 template <typename Scalar>
 struct TriangleMotionBoundVisitorVisitImpl<Scalar, SplineMotion<Scalar>>
 {
-  Scalar operator()(
+  static Scalar run(
       const TriangleMotionBoundVisitor<Scalar>& visitor,
       const SplineMotion<Scalar>& motion)
   {
@@ -564,7 +568,7 @@ struct TriangleMotionBoundVisitorVisitImpl<Scalar, SplineMotion<Scalar>>
 template <typename Scalar>
 struct TriangleMotionBoundVisitorVisitImpl<Scalar, TranslationMotion<Scalar>>
 {
-  Scalar operator()(
+  static Scalar run(
       const TriangleMotionBoundVisitor<Scalar>& visitor,
       const TranslationMotion<Scalar>& motion)
   {
diff --git a/include/fcl/collision_func_matrix.h b/include/fcl/collision_func_matrix.h
index 88dbc4fb5..c10fd77fb 100644
--- a/include/fcl/collision_func_matrix.h
+++ b/include/fcl/collision_func_matrix.h
@@ -529,7 +529,7 @@ struct BVHShapeCollider<OBBRSS<typename NarrowPhaseSolver::Scalar>, T_SH, Narrow
 template <typename Scalar, typename T_BVH>
 struct BVHCollideImpl
 {
-  std::size_t operator()(
+  static std::size_t run(
       const CollisionGeometry<Scalar>* o1,
       const Transform3<Scalar>& tf1,
       const CollisionGeometry<Scalar>* o2,
@@ -567,8 +567,8 @@ std::size_t BVHCollide(
     const CollisionRequest<typename T_BVH::Scalar>& request,
     CollisionResult<typename T_BVH::Scalar>& result)
 {
-  BVHCollideImpl<typename T_BVH::Scalar, T_BVH> tmp;
-  return tmp(o1, tf1, o2, tf2, request, result);
+  return BVHCollideImpl<typename T_BVH::Scalar, T_BVH>::run(
+        o1, tf1, o2, tf2, request, result);
 }
 
 namespace details
@@ -602,7 +602,7 @@ std::size_t orientedMeshCollide(
 template <typename Scalar>
 struct BVHCollideImpl<Scalar, OBB<Scalar>>
 {
-  std::size_t operator()(
+  static std::size_t run(
       const CollisionGeometry<Scalar>* o1,
       const Transform3<Scalar>& tf1,
       const CollisionGeometry<Scalar>* o2,
@@ -620,7 +620,7 @@ struct BVHCollideImpl<Scalar, OBB<Scalar>>
 template <typename Scalar>
 struct BVHCollideImpl<Scalar, OBBRSS<Scalar>>
 {
-  std::size_t operator()(
+  static std::size_t run(
       const CollisionGeometry<Scalar>* o1,
       const Transform3<Scalar>& tf1,
       const CollisionGeometry<Scalar>* o2,
@@ -638,7 +638,7 @@ struct BVHCollideImpl<Scalar, OBBRSS<Scalar>>
 template <typename Scalar>
 struct BVHCollideImpl<Scalar, kIOS<Scalar>>
 {
-  std::size_t operator()(
+  static std::size_t run(
       const CollisionGeometry<Scalar>* o1,
       const Transform3<Scalar>& tf1,
       const CollisionGeometry<Scalar>* o2,
diff --git a/include/fcl/distance_func_matrix.h b/include/fcl/distance_func_matrix.h
index 644a517ac..636645b32 100644
--- a/include/fcl/distance_func_matrix.h
+++ b/include/fcl/distance_func_matrix.h
@@ -347,7 +347,7 @@ struct BVHShapeDistancer<OBBRSS<typename T_SH::Scalar>, T_SH, NarrowPhaseSolver>
 template <typename Scalar, typename T_BVH>
 struct BVHDistanceImpl
 {
-  Scalar operator()(
+  static Scalar run(
       const CollisionGeometry<Scalar>* o1,
       const Transform3<Scalar>& tf1,
       const CollisionGeometry<Scalar>* o2,
@@ -383,8 +383,8 @@ typename T_BVH::Scalar BVHDistance(
     const DistanceRequest<typename T_BVH::Scalar>& request,
     DistanceResult<typename T_BVH::Scalar>& result)
 {
-  BVHDistanceImpl<typename T_BVH::Scalar, T_BVH> tmp;
-  return tmp(o1, tf1, o2, tf2, request, result);
+  return BVHDistanceImpl<typename T_BVH::Scalar, T_BVH>::run(
+        o1, tf1, o2, tf2, request, result);
 }
 
 namespace details
@@ -416,7 +416,7 @@ typename T_BVH::Scalar orientedMeshDistance(
 template <typename Scalar>
 struct BVHDistanceImpl<Scalar, RSS<Scalar>>
 {
-  Scalar operator()(
+  static Scalar run(
       const CollisionGeometry<Scalar>* o1,
       const Transform3<Scalar>& tf1,
       const CollisionGeometry<Scalar>* o2,
@@ -434,7 +434,7 @@ struct BVHDistanceImpl<Scalar, RSS<Scalar>>
 template <typename Scalar>
 struct BVHDistanceImpl<Scalar, kIOS<Scalar>>
 {
-  Scalar operator()(
+  static Scalar run(
       const CollisionGeometry<Scalar>* o1,
       const Transform3<Scalar>& tf1,
       const CollisionGeometry<Scalar>* o2,
@@ -452,7 +452,7 @@ struct BVHDistanceImpl<Scalar, kIOS<Scalar>>
 template <typename Scalar>
 struct BVHDistanceImpl<Scalar, OBBRSS<Scalar>>
 {
-  Scalar operator()(
+  static Scalar run(
       const CollisionGeometry<Scalar>* o1,
       const Transform3<Scalar>& tf1,
       const CollisionGeometry<Scalar>* o2,
diff --git a/include/fcl/narrowphase/gjk_solver_indep.h b/include/fcl/narrowphase/gjk_solver_indep.h
index 5ee05895c..50ce95764 100755
--- a/include/fcl/narrowphase/gjk_solver_indep.h
+++ b/include/fcl/narrowphase/gjk_solver_indep.h
@@ -221,7 +221,7 @@ bool GJKSolver_indep<Scalar>::shapeIntersect(const S1& s1, const Transform3<Scal
 template<typename Scalar, typename S1, typename S2>
 struct ShapeIntersectIndepImpl
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& gjkSolver,
       const S1& s1,
       const Transform3<Scalar>& tf1,
@@ -285,8 +285,8 @@ bool GJKSolver_indep<Scalar>::shapeIntersect(
     const Transform3<Scalar>& tf2,
     std::vector<ContactPoint<Scalar>>* contacts) const
 {
-  ShapeIntersectIndepImpl<Scalar, S1, S2> shapeIntersectImpl;
-  return shapeIntersectImpl(*this, s1, tf1, s2, tf2, contacts);
+  return ShapeIntersectIndepImpl<Scalar, S1, S2>::run(
+        *this, s1, tf1, s2, tf2, contacts);
 }
 
 // Shape intersect algorithms not using built-in GJK algorithm
@@ -317,7 +317,7 @@ bool GJKSolver_indep<Scalar>::shapeIntersect(
   template <typename Scalar>\
   struct ShapeIntersectIndepImpl<Scalar, SHAPE1<Scalar>, SHAPE2<Scalar>>\
   {\
-    bool operator()(\
+    static bool run(\
         const GJKSolver_indep<Scalar>& /*gjkSolver*/,\
         const SHAPE1<Scalar>& s1,\
         const Transform3<Scalar>& tf1,\
@@ -333,7 +333,7 @@ bool GJKSolver_indep<Scalar>::shapeIntersect(
   template <typename Scalar>\
   struct ShapeIntersectIndepImpl<Scalar, SHAPE2<Scalar>, SHAPE1<Scalar>>\
   {\
-    bool operator()(\
+    static bool run(\
         const GJKSolver_indep<Scalar>& /*gjkSolver*/,\
         const SHAPE2<Scalar>& s1,\
         const Transform3<Scalar>& tf1,\
@@ -376,7 +376,7 @@ FCL_GJK_INDEP_SHAPE_SHAPE_INTERSECT(Cone, Plane, details::conePlaneIntersect)
 template <typename Scalar>
 struct ShapeIntersectIndepImpl<Scalar, Halfspace<Scalar>, Halfspace<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& /*gjkSolver*/,
       const Halfspace<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -395,7 +395,7 @@ struct ShapeIntersectIndepImpl<Scalar, Halfspace<Scalar>, Halfspace<Scalar>>
 template <typename Scalar>
 struct ShapeIntersectIndepImpl<Scalar, Plane<Scalar>, Plane<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& /*gjkSolver*/,
       const Plane<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -410,7 +410,7 @@ struct ShapeIntersectIndepImpl<Scalar, Plane<Scalar>, Plane<Scalar>>
 template <typename Scalar>
 struct ShapeIntersectIndepImpl<Scalar, Plane<Scalar>, Halfspace<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& /*gjkSolver*/,
       const Plane<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -429,7 +429,7 @@ struct ShapeIntersectIndepImpl<Scalar, Plane<Scalar>, Halfspace<Scalar>>
 template <typename Scalar>
 struct ShapeIntersectIndepImpl<Scalar, Halfspace<Scalar>, Plane<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& /*gjkSolver*/,
       const Halfspace<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -449,7 +449,7 @@ struct ShapeIntersectIndepImpl<Scalar, Halfspace<Scalar>, Plane<Scalar>>
 template<typename Scalar, typename S>
 struct ShapeTriangleIntersectIndepImpl
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& gjkSolver,
       const S& s,
       const Transform3<Scalar>& tf,
@@ -516,8 +516,7 @@ bool GJKSolver_indep<Scalar>::shapeTriangleIntersect(
     Scalar* penetration_depth,
     Vector3<Scalar>* normal) const
 {
-  ShapeTriangleIntersectIndepImpl<Scalar, S> shapeTriangleIntersectImpl;
-  return shapeTriangleIntersectImpl(
+  return ShapeTriangleIntersectIndepImpl<Scalar, S>::run(
         *this, s, tf, P1, P2, P3, contact_points, penetration_depth, normal);
 }
 
@@ -525,7 +524,7 @@ bool GJKSolver_indep<Scalar>::shapeTriangleIntersect(
 template<typename Scalar>
 struct ShapeTriangleIntersectIndepImpl<Scalar, Sphere<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& /*gjkSolver*/,
       const Sphere<Scalar>& s,
       const Transform3<Scalar>& tf,
@@ -546,7 +545,7 @@ struct ShapeTriangleIntersectIndepImpl<Scalar, Sphere<Scalar>>
 template<typename Scalar, typename S>
 struct ShapeTransformedTriangleIntersectIndepImpl
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& gjkSolver,
       const S& s,
       const Transform3<Scalar>& tf1,
@@ -615,8 +614,7 @@ bool GJKSolver_indep<Scalar>::shapeTriangleIntersect(
     Scalar* penetration_depth,
     Vector3<Scalar>* normal) const
 {
-  ShapeTransformedTriangleIntersectIndepImpl<Scalar, S> shapeTriangleIntersectImpl;
-  return shapeTriangleIntersectImpl(
+  return ShapeTransformedTriangleIntersectIndepImpl<Scalar, S>::run(
         *this, s, tf1, P1, P2, P3, tf2,
         contact_points, penetration_depth, normal);
 }
@@ -625,7 +623,7 @@ bool GJKSolver_indep<Scalar>::shapeTriangleIntersect(
 template<typename Scalar>
 struct ShapeTransformedTriangleIntersectIndepImpl<Scalar, Sphere<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& /*gjkSolver*/,
       const Sphere<Scalar>& s,
       const Transform3<Scalar>& tf1,
@@ -647,7 +645,7 @@ struct ShapeTransformedTriangleIntersectIndepImpl<Scalar, Sphere<Scalar>>
 template<typename Scalar>
 struct ShapeTransformedTriangleIntersectIndepImpl<Scalar, Halfspace<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& /*gjkSolver*/,
       const Halfspace<Scalar>& s,
       const Transform3<Scalar>& tf1,
@@ -669,7 +667,7 @@ struct ShapeTransformedTriangleIntersectIndepImpl<Scalar, Halfspace<Scalar>>
 template<typename Scalar>
 struct ShapeTransformedTriangleIntersectIndepImpl<Scalar, Plane<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& /*gjkSolver*/,
       const Plane<Scalar>& s,
       const Transform3<Scalar>& tf1,
@@ -691,7 +689,7 @@ struct ShapeTransformedTriangleIntersectIndepImpl<Scalar, Plane<Scalar>>
 template<typename Scalar, typename S1, typename S2>
 struct ShapeDistanceIndepImpl
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& gjkSolver,
       const S1& s1,
       const Transform3<Scalar>& tf1,
@@ -751,8 +749,8 @@ bool GJKSolver_indep<Scalar>::shapeDistance(
     Vector3<Scalar>* p1,
     Vector3<Scalar>* p2) const
 {
-  ShapeDistanceIndepImpl<Scalar, S1, S2> shapeDistanceImpl;
-  return shapeDistanceImpl(*this, s1, tf1, s2, tf2, dist, p1, p2);
+  return ShapeDistanceIndepImpl<Scalar, S1, S2>::run(
+        *this, s1, tf1, s2, tf2, dist, p1, p2);
 }
 
 // Shape distance algorithms not using built-in GJK algorithm
@@ -783,7 +781,7 @@ bool GJKSolver_indep<Scalar>::shapeDistance(
 template<typename Scalar>
 struct ShapeDistanceIndepImpl<Scalar, Sphere<Scalar>, Capsule<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& /*gjkSolver*/,
       const Sphere<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -801,7 +799,7 @@ struct ShapeDistanceIndepImpl<Scalar, Sphere<Scalar>, Capsule<Scalar>>
 template<typename Scalar>
 struct ShapeDistanceIndepImpl<Scalar, Capsule<Scalar>, Sphere<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& /*gjkSolver*/,
       const Capsule<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -819,7 +817,7 @@ struct ShapeDistanceIndepImpl<Scalar, Capsule<Scalar>, Sphere<Scalar>>
 template<typename Scalar>
 struct ShapeDistanceIndepImpl<Scalar, Sphere<Scalar>, Sphere<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& /*gjkSolver*/,
       const Sphere<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -837,7 +835,7 @@ struct ShapeDistanceIndepImpl<Scalar, Sphere<Scalar>, Sphere<Scalar>>
 template<typename Scalar>
 struct ShapeDistanceIndepImpl<Scalar, Capsule<Scalar>, Capsule<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& /*gjkSolver*/,
       const Capsule<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -855,7 +853,7 @@ struct ShapeDistanceIndepImpl<Scalar, Capsule<Scalar>, Capsule<Scalar>>
 template<typename Scalar, typename S>
 struct ShapeTriangleDistanceIndepImpl
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& gjkSolver,
       const S& s,
       const Transform3<Scalar>& tf,
@@ -917,8 +915,7 @@ bool GJKSolver_indep<Scalar>::shapeTriangleDistance(
     Vector3<Scalar>* p1,
     Vector3<Scalar>* p2) const
 {
-  ShapeTriangleDistanceIndepImpl<Scalar, S> shapeTriangleDistanceImpl;
-  return shapeTriangleDistanceImpl(
+  return ShapeTriangleDistanceIndepImpl<Scalar, S>::run(
         *this, s, tf, P1, P2, P3, dist, p1, p2);
 }
 
@@ -926,7 +923,7 @@ bool GJKSolver_indep<Scalar>::shapeTriangleDistance(
 template<typename Scalar>
 struct ShapeTriangleDistanceIndepImpl<Scalar, Sphere<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& /*gjkSolver*/,
       const Sphere<Scalar>& s,
       const Transform3<Scalar>& tf,
@@ -945,7 +942,7 @@ struct ShapeTriangleDistanceIndepImpl<Scalar, Sphere<Scalar>>
 template<typename Scalar, typename S>
 struct ShapeTransformedTriangleDistanceIndepImpl
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& gjkSolver,
       const S& s,
       const Transform3<Scalar>& tf1,
@@ -1009,8 +1006,7 @@ bool GJKSolver_indep<Scalar>::shapeTriangleDistance(
     Vector3<Scalar>* p1,
     Vector3<Scalar>* p2) const
 {
-  ShapeTransformedTriangleDistanceIndepImpl<Scalar, S> shapeTriangleDistanceImpl;
-  return shapeTriangleDistanceImpl(
+  return ShapeTransformedTriangleDistanceIndepImpl<Scalar, S>::run(
         *this, s, tf1, P1, P2, P3, tf2, dist, p1, p2);
 }
 
@@ -1018,7 +1014,7 @@ bool GJKSolver_indep<Scalar>::shapeTriangleDistance(
 template<typename Scalar>
 struct ShapeTransformedTriangleDistanceIndepImpl<Scalar, Sphere<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_indep<Scalar>& /*gjkSolver*/,
       const Sphere<Scalar>& s,
       const Transform3<Scalar>& tf1,
diff --git a/include/fcl/narrowphase/gjk_solver_libccd.h b/include/fcl/narrowphase/gjk_solver_libccd.h
index 7c7f9447b..4bbb2ba47 100755
--- a/include/fcl/narrowphase/gjk_solver_libccd.h
+++ b/include/fcl/narrowphase/gjk_solver_libccd.h
@@ -213,7 +213,7 @@ bool GJKSolver_libccd<Scalar>::shapeIntersect(
 template<typename Scalar, typename S1, typename S2>
 struct ShapeIntersectLibccdImpl
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& gjkSolver,
       const S1& s1, const Transform3<Scalar>& tf1,
       const S2& s2, const Transform3<Scalar>& tf2,
@@ -273,8 +273,8 @@ bool GJKSolver_libccd<Scalar>::shapeIntersect(
     const S2& s2, const Transform3<Scalar>& tf2,
     std::vector<ContactPoint<Scalar>>* contacts) const
 {
-  ShapeIntersectLibccdImpl<Scalar, S1, S2> shapeIntersectImpl;
-  return shapeIntersectImpl(*this, s1, tf1, s2, tf2, contacts);
+  return ShapeIntersectLibccdImpl<Scalar, S1, S2>::run(
+        *this, s1, tf1, s2, tf2, contacts);
 }
 
 // Shape intersect algorithms not using libccd
@@ -305,7 +305,7 @@ bool GJKSolver_libccd<Scalar>::shapeIntersect(
   template <typename Scalar>\
   struct ShapeIntersectLibccdImpl<Scalar, SHAPE1<Scalar>, SHAPE2<Scalar>>\
   {\
-    bool operator()(\
+    static bool run(\
         const GJKSolver_libccd<Scalar>& /*gjkSolver*/,\
         const SHAPE1<Scalar>& s1,\
         const Transform3<Scalar>& tf1,\
@@ -321,7 +321,7 @@ bool GJKSolver_libccd<Scalar>::shapeIntersect(
   template <typename Scalar>\
   struct ShapeIntersectLibccdImpl<Scalar, SHAPE2<Scalar>, SHAPE1<Scalar>>\
   {\
-    bool operator()(\
+    static bool run(\
         const GJKSolver_libccd<Scalar>& /*gjkSolver*/,\
         const SHAPE2<Scalar>& s1,\
         const Transform3<Scalar>& tf1,\
@@ -364,7 +364,7 @@ FCL_GJK_LIBCCD_SHAPE_SHAPE_INTERSECT(Cone, Plane, details::conePlaneIntersect)
 template <typename Scalar>
 struct ShapeIntersectLibccdImpl<Scalar, Halfspace<Scalar>, Halfspace<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& /*gjkSolver*/,
       const Halfspace<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -383,7 +383,7 @@ struct ShapeIntersectLibccdImpl<Scalar, Halfspace<Scalar>, Halfspace<Scalar>>
 template <typename Scalar>
 struct ShapeIntersectLibccdImpl<Scalar, Plane<Scalar>, Plane<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& /*gjkSolver*/,
       const Plane<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -398,7 +398,7 @@ struct ShapeIntersectLibccdImpl<Scalar, Plane<Scalar>, Plane<Scalar>>
 template <typename Scalar>
 struct ShapeIntersectLibccdImpl<Scalar, Plane<Scalar>, Halfspace<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& /*gjkSolver*/,
       const Plane<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -417,7 +417,7 @@ struct ShapeIntersectLibccdImpl<Scalar, Plane<Scalar>, Halfspace<Scalar>>
 template <typename Scalar>
 struct ShapeIntersectLibccdImpl<Scalar, Halfspace<Scalar>, Plane<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& /*gjkSolver*/,
       const Halfspace<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -437,7 +437,7 @@ struct ShapeIntersectLibccdImpl<Scalar, Halfspace<Scalar>, Plane<Scalar>>
 template<typename Scalar, typename S>
 struct ShapeTriangleIntersectLibccdImpl
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& gjkSolver,
       const S& s,
       const Transform3<Scalar>& tf,
@@ -483,8 +483,7 @@ bool GJKSolver_libccd<Scalar>::shapeTriangleIntersect(
     Scalar* penetration_depth,
     Vector3<Scalar>* normal) const
 {
-  ShapeTriangleIntersectLibccdImpl<Scalar, S> shapeTriangleIntersectImpl;
-  return shapeTriangleIntersectImpl(
+  return ShapeTriangleIntersectLibccdImpl<Scalar, S>::run(
         *this, s, tf, P1, P2, P3, contact_points, penetration_depth, normal);
 }
 
@@ -492,7 +491,7 @@ bool GJKSolver_libccd<Scalar>::shapeTriangleIntersect(
 template<typename Scalar>
 struct ShapeTriangleIntersectLibccdImpl<Scalar, Sphere<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& /*gjkSolver*/,
       const Sphere<Scalar>& s,
       const Transform3<Scalar>& tf,
@@ -512,7 +511,7 @@ struct ShapeTriangleIntersectLibccdImpl<Scalar, Sphere<Scalar>>
 template<typename Scalar, typename S>
 struct ShapeTransformedTriangleIntersectLibccdImpl
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& gjkSolver,
       const S& s,
       const Transform3<Scalar>& tf1,
@@ -560,8 +559,7 @@ bool GJKSolver_libccd<Scalar>::shapeTriangleIntersect(
     Scalar* penetration_depth,
     Vector3<Scalar>* normal) const
 {
-  ShapeTransformedTriangleIntersectLibccdImpl<Scalar, S> shapeTriangleIntersectImpl;
-  return shapeTriangleIntersectImpl(
+  return ShapeTransformedTriangleIntersectLibccdImpl<Scalar, S>::run(
         *this, s, tf1, P1, P2, P3, tf2,
         contact_points, penetration_depth, normal);
 }
@@ -570,7 +568,7 @@ bool GJKSolver_libccd<Scalar>::shapeTriangleIntersect(
 template<typename Scalar>
 struct ShapeTransformedTriangleIntersectLibccdImpl<Scalar, Sphere<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& /*gjkSolver*/,
       const Sphere<Scalar>& s,
       const Transform3<Scalar>& tf1,
@@ -592,7 +590,7 @@ struct ShapeTransformedTriangleIntersectLibccdImpl<Scalar, Sphere<Scalar>>
 template<typename Scalar>
 struct ShapeTransformedTriangleIntersectLibccdImpl<Scalar, Halfspace<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& /*gjkSolver*/,
       const Halfspace<Scalar>& s,
       const Transform3<Scalar>& tf1,
@@ -614,7 +612,7 @@ struct ShapeTransformedTriangleIntersectLibccdImpl<Scalar, Halfspace<Scalar>>
 template<typename Scalar>
 struct ShapeTransformedTriangleIntersectLibccdImpl<Scalar, Plane<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& /*gjkSolver*/,
       const Plane<Scalar>& s,
       const Transform3<Scalar>& tf1,
@@ -636,7 +634,7 @@ struct ShapeTransformedTriangleIntersectLibccdImpl<Scalar, Plane<Scalar>>
 template<typename Scalar, typename S1, typename S2>
 struct ShapeDistanceLibccdImpl
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& gjkSolver,
       const S1& s1,
       const Transform3<Scalar>& tf1,
@@ -684,8 +682,8 @@ bool GJKSolver_libccd<Scalar>::shapeDistance(
     Vector3<Scalar>* p1,
     Vector3<Scalar>* p2) const
 {
-  ShapeDistanceLibccdImpl<Scalar, S1, S2> shapeDistanceImpl;
-  return shapeDistanceImpl(*this, s1, tf1, s2, tf2, dist, p1, p2);
+  return ShapeDistanceLibccdImpl<Scalar, S1, S2>::run(
+        *this, s1, tf1, s2, tf2, dist, p1, p2);
 }
 
 // Shape distance algorithms not using libccd
@@ -716,7 +714,7 @@ bool GJKSolver_libccd<Scalar>::shapeDistance(
 template<typename Scalar>
 struct ShapeDistanceLibccdImpl<Scalar, Sphere<Scalar>, Capsule<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& /*gjkSolver*/,
       const Sphere<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -734,7 +732,7 @@ struct ShapeDistanceLibccdImpl<Scalar, Sphere<Scalar>, Capsule<Scalar>>
 template<typename Scalar>
 struct ShapeDistanceLibccdImpl<Scalar, Capsule<Scalar>, Sphere<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& /*gjkSolver*/,
       const Capsule<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -752,7 +750,7 @@ struct ShapeDistanceLibccdImpl<Scalar, Capsule<Scalar>, Sphere<Scalar>>
 template<typename Scalar>
 struct ShapeDistanceLibccdImpl<Scalar, Sphere<Scalar>, Sphere<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& /*gjkSolver*/,
       const Sphere<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -770,7 +768,7 @@ struct ShapeDistanceLibccdImpl<Scalar, Sphere<Scalar>, Sphere<Scalar>>
 template<typename Scalar>
 struct ShapeDistanceLibccdImpl<Scalar, Capsule<Scalar>, Capsule<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& /*gjkSolver*/,
       const Capsule<Scalar>& s1,
       const Transform3<Scalar>& tf1,
@@ -788,7 +786,7 @@ struct ShapeDistanceLibccdImpl<Scalar, Capsule<Scalar>, Capsule<Scalar>>
 template<typename Scalar, typename S>
 struct ShapeTriangleDistanceLibccdImpl
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& gjkSolver,
       const S& s,
       const Transform3<Scalar>& tf,
@@ -834,8 +832,7 @@ bool GJKSolver_libccd<Scalar>::shapeTriangleDistance(
     Vector3<Scalar>* p1,
     Vector3<Scalar>* p2) const
 {
-  ShapeTriangleDistanceLibccdImpl<Scalar, S> shapeTriangleDistanceImpl;
-  return shapeTriangleDistanceImpl(
+  return ShapeTriangleDistanceLibccdImpl<Scalar, S>::run(
         *this, s, tf, P1, P2, P3, dist, p1, p2);
 }
 
@@ -843,7 +840,7 @@ bool GJKSolver_libccd<Scalar>::shapeTriangleDistance(
 template<typename Scalar>
 struct ShapeTriangleDistanceLibccdImpl<Scalar, Sphere<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& /*gjkSolver*/,
       const Sphere<Scalar>& s,
       const Transform3<Scalar>& tf,
@@ -862,7 +859,7 @@ struct ShapeTriangleDistanceLibccdImpl<Scalar, Sphere<Scalar>>
 template<typename Scalar, typename S>
 struct ShapeTransformedTriangleDistanceLibccdImpl
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& gjkSolver,
       const S& s,
       const Transform3<Scalar>& tf1,
@@ -911,8 +908,7 @@ bool GJKSolver_libccd<Scalar>::shapeTriangleDistance(
     Vector3<Scalar>* p1,
     Vector3<Scalar>* p2) const
 {
-  ShapeTransformedTriangleDistanceLibccdImpl<Scalar, S> shapeTriangleDistanceImpl;
-  return shapeTriangleDistanceImpl(
+  return ShapeTransformedTriangleDistanceLibccdImpl<Scalar, S>::run(
         *this, s, tf1, P1, P2, P3, tf2, dist, p1, p2);
 }
 
@@ -920,7 +916,7 @@ bool GJKSolver_libccd<Scalar>::shapeTriangleDistance(
 template<typename Scalar>
 struct ShapeTransformedTriangleDistanceLibccdImpl<Scalar, Sphere<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const GJKSolver_libccd<Scalar>& /*gjkSolver*/,
       const Sphere<Scalar>& s,
       const Transform3<Scalar>& tf1,
diff --git a/include/fcl/traversal/distance/mesh_conservative_advancement_traversal_node.h b/include/fcl/traversal/distance/mesh_conservative_advancement_traversal_node.h
index ccfebe510..a29f6c242 100644
--- a/include/fcl/traversal/distance/mesh_conservative_advancement_traversal_node.h
+++ b/include/fcl/traversal/distance/mesh_conservative_advancement_traversal_node.h
@@ -323,7 +323,7 @@ void MeshConservativeAdvancementTraversalNode<BV>::leafTesting(int b1, int b2) c
 template <typename Scalar, typename BV>
 struct CanStopImpl
 {
-  bool operator()(
+  static bool run(
       const MeshConservativeAdvancementTraversalNode<BV>& node, Scalar c)
   {
     if((c >= node.w * (node.min_distance - node.abs_err))
@@ -389,15 +389,14 @@ template <typename BV>
 bool MeshConservativeAdvancementTraversalNode<BV>::canStop(
     MeshConservativeAdvancementTraversalNode<BV>::Scalar c) const
 {
-  CanStopImpl<typename BV::Scalar, BV> canStopImpl;
-  return canStopImpl(*this, c);
+  return CanStopImpl<typename BV::Scalar, BV>::run(*this, c);
 }
 
 //==============================================================================
 template <typename Scalar>
 struct CanStopImpl<Scalar, OBB<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const MeshConservativeAdvancementTraversalNode<OBB<Scalar>>& node,
       Scalar c)
   {
@@ -420,7 +419,7 @@ struct CanStopImpl<Scalar, OBB<Scalar>>
 template <typename Scalar>
 struct CanStopImpl<Scalar, RSS<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const MeshConservativeAdvancementTraversalNode<RSS<Scalar>>& node,
       Scalar c)
   {
@@ -443,7 +442,7 @@ struct CanStopImpl<Scalar, RSS<Scalar>>
 template <typename Scalar>
 struct CanStopImpl<Scalar, OBBRSS<Scalar>>
 {
-  bool operator()(
+  static bool run(
       const MeshConservativeAdvancementTraversalNode<OBBRSS<Scalar>>& node,
       Scalar c)
   {