libk3dmesh::detail::position_t libk3dmesh::detail::normal(const position_t& p1, const position_t& p2, const position_t& p3, const position_t& p4)
{
 // normal calculation taken from k3d::normal()
 position_t n(0, 0, 0);

 std::cerr << debug << "p1: " << p1[0] << ", " << p1[1] << ", " << p1[2] << std::endl;
 std::cerr << debug << "p2: " << p2[0] << ", " << p2[1] << ", " << p2[2] << std::endl;
 std::cerr << debug << "p3: " << p3[0] << ", " << p3[1] << ", " << p3[2] << std::endl;
 std::cerr << debug << "p4: " << p4[0] << ", " << p4[1] << ", " << p4[2] << std::endl;
 
 position_t Points[4] = {p1,p2,p3,p4};
 
 for(unsigned long k = 0; k < 4; ++k)
 {
  const position_t& i = *Points[k];
  const position_t& j = *Points[(k+1)%4];

  std::cerr << debug << "i: " << i[0] << ", " << i[1] << ", " << i[2] << std::endl;
  std::cerr << debug << "j: " << j[0] << ", " << j[1] << ", " << j[2] << std::endl;
  
  n[0] += (i[1] + j[1]) * (j[2] - i[2]);
  n[1] += (i[2] + j[2]) * (j[0] - i[0]);
  n[2] += (i[0] + j[0]) * (j[1] - i[1]);
  std::cerr << debug << "n: " << n[0] << ", " << n[1] << ", " << n[2] << std::endl;

 }
 return 0.5 * n;
}