Merge pull request #667 from moritzdannhauer/origin/CleaverFix

Origin/cleaver fix

Author: dcwhite

src/Core/Algorithms/Field/InterfaceWithCleaverAlgorithm.cc

@@ -31,7 +31,10 @@ ToDo: Padding is always enabled because of exit() in cleaver lib
 */
 
 ///TODO: fix include path to remove Externals/ part
+
 #include <Externals/cleaver/lib/FloatField.h>
+#include <Externals/cleaver/lib/vec3.h>
+#include <Externals/cleaver/lib/BoundingBox.h>
 #include <Externals/cleaver/lib/Cleaver.h>
 #include <Externals/cleaver/lib/InverseField.h>
 #include <Externals/cleaver/lib/PaddedVolume.h>
@@ -49,6 +52,7 @@ ToDo: Padding is always enabled because of exit() in cleaver lib
 #include <Core/Utils/StringUtil.h>
 #include <boost/scoped_ptr.hpp>
 #include <Core/Logging/Log.h>
+#include <Core/Math/MiscMath.h>
 
 using namespace SCIRun::Core::Algorithms;
 using namespace SCIRun::Core::Algorithms::Fields;
@@ -74,38 +78,41 @@ InterfaceWithCleaverAlgorithm::InterfaceWithCleaverAlgorithm()
   addParameter(VolumeScalingZ,1.0);  
 }
 
-boost::shared_ptr<Cleaver::ScalarField> InterfaceWithCleaverAlgorithm::makeCleaverFieldFromLatVol(FieldHandle field)
+boost::shared_ptr<Cleaver::ScalarField> InterfaceWithCleaverAlgorithm::makeCleaverFieldFromLatVol(FieldHandle field )
 {
-  //TODO: this function assumes input is completely valid, may want to move various checks from run() function here.
-
   VMesh*  vmesh   = field->vmesh();
   VField* vfield = field->vfield();
   VMesh::dimension_type dims;
   vmesh->get_dimensions( dims ); 
+  
   float* ptr = static_cast<float*>(vfield->fdata_pointer());
-  auto cleaverField = boost::make_shared<Cleaver::FloatField>(dims[0], dims[1], dims[2], ptr);
-
-  //0 = constant, 1 = linear
-  if (0 == vfield->basis_order())
-    cleaverField->setCenter(Cleaver::FloatField::CellCentered);
-  else if (1 == vfield->basis_order())
-    cleaverField->setCenter(Cleaver::FloatField::NodeCentered);
-  //else: TODO? handle other bases, probably by throwing exception.
-
-  //TODO: Cleaver FloatField needs more setup (constructor is not sufficient)
-  // 1.
-  // setBounds(BoundingBox). Convert vmesh->get_bounding_box() to Cleaver BoundingBox.
-  // 2.
-  // setScale(vec3). Need to figure out which vmesh function to call, and convert to Cleaver::vec3.
-  // 3. unit test heavily.
-
+  
+  auto cleaverField = boost::make_shared<Cleaver::FloatField>(dims[0], dims[1], dims[2], ptr);  
+    
+  BBox bbox=vmesh->get_bounding_box();
+  Point bmin, bmax;
+  if (bbox.valid())
+  {
+   bmin = bbox.min();
+   bmax = bbox.max();
+  }
+  Cleaver::BoundingBox  bb = Cleaver::BoundingBox(Cleaver::vec3::zero, Cleaver::vec3(dims[0],dims[1],dims[2]));
+  cleaverField->setBounds(bb);
+  const Transform &transform = vmesh->get_transform();
+  
+  int x_spacing=fabs(transform.get_mat_val(0,0)), y_spacing=fabs(transform.get_mat_val(1,1)), z_spacing=fabs(transform.get_mat_val(2,2));
+  if (IsNan(x_spacing)) x_spacing=1;
+  if (IsNan(y_spacing)) y_spacing=1;
+  if (IsNan(z_spacing)) z_spacing=1;
+  
+  cleaverField->setScale(Cleaver::vec3(x_spacing,y_spacing,z_spacing));
+  
   return cleaverField;
 }
 
 FieldHandle InterfaceWithCleaverAlgorithm::run(const std::vector<FieldHandle>& input) const
 {
   FieldHandle output;      
-
   std::vector<FieldHandle> inputs;
   std::copy_if(input.begin(), input.end(), std::back_inserter(inputs), [](FieldHandle f) { return f; });
 
@@ -119,10 +126,11 @@ FieldHandle InterfaceWithCleaverAlgorithm::run(const std::vector<FieldHandle>& i
     THROW_ALGORITHM_INPUT_ERROR(" At least 2 indicator functions stored as float values are needed to run cleaver! " );
     return FieldHandle();
   }
-
+  
   std::vector<boost::shared_ptr<Cleaver::ScalarField>> fields;  
+
   VMesh::dimension_type dims; int x=0,y=0,z=0; 
-  for (size_t p=1; p<inputs.size(); p++)
+  for (size_t p=0; p<inputs.size(); p++)
   {
     FieldHandle input = inputs[p];
     VMesh*  imesh1   = input->vmesh();
@@ -141,7 +149,7 @@ FieldHandle InterfaceWithCleaverAlgorithm::run(const std::vector<FieldHandle>& i
       }
 
       imesh1->get_dimensions( dims ); 
-      if (p==1)
+      if (p==0)
       {
         x=dims[0]; y=dims[1]; z=dims[2];
         if (x<1 || y<1 || z<1)
@@ -163,18 +171,27 @@ FieldHandle InterfaceWithCleaverAlgorithm::run(const std::vector<FieldHandle>& i
         return FieldHandle();
       }    
 
+      //0 = constant, 1 = linear
+       if (1 != vfield1->basis_order())
+       {
+        THROW_ALGORITHM_INPUT_ERROR("Input data need to be defined on input mesh nodes.");
+       } 
+      
       if (vfield1->is_float())
-      {
+      {    
         float* ptr = static_cast<float*>(vfield1->fdata_pointer());
-        if (ptr)
-        {
+	if (ptr)
+        {	
           fields.push_back(makeCleaverFieldFromLatVol(input));
         }
         else
         {
           THROW_ALGORITHM_INPUT_ERROR(" float field is NULL pointer");
           return FieldHandle();
         }
+      } else
+      {
+       THROW_ALGORITHM_INPUT_ERROR(" Input field needs to be a structured mesh (best would be a LatVol) with float values defnied on the elements. "); 
       }
 
     }
@@ -199,35 +216,43 @@ FieldHandle InterfaceWithCleaverAlgorithm::run(const std::vector<FieldHandle>& i
       double newX = xScale*volume->size().x;
       double newY = yScale*volume->size().y;
       double newZ = zScale*volume->size().z;
-      std::cout << "Cleaver setting volume size to " << newX << "," << newY << "," << newZ << std::endl;
       volume->setSize(newX, newY, newZ);
     }
     else // None
+    {
       volume->setSize(dims[0],dims[1],dims[2]);
+    }
   }
   else
   {
     THROW_ALGORITHM_INPUT_ERROR(" Invalid Scaling. Use Input sizes.");
   }
-
+    
   /// Padding is now optional! 
   boost::shared_ptr<Cleaver::AbstractVolume> paddedVolume(volume);
-
   const bool verbose = get(Verbose).toBool();
   const bool pad = get(Padding).toBool();
   if (pad)
     paddedVolume.reset(new Cleaver::PaddedVolume(volume.get()));
+    
+  std::cout << "Creating Mesh with Volume Size " << paddedVolume->size().toString() << std::endl;
+    
   boost::scoped_ptr<Cleaver::TetMesh> mesh(Cleaver::createMeshFromVolume(paddedVolume.get(), verbose));
 
-  FieldInformation fi("TetVolMesh",0,"double");   ///create output field
+  auto nr_of_tets  = mesh->tets.size();
+  auto nr_of_verts = mesh->verts.size();
+
+  if (nr_of_tets==0 || nr_of_verts==0)
+  {
+    THROW_ALGORITHM_INPUT_ERROR(" Number of resulting tetrahedral nodes or elements is 0. If you disabled padding enable it and execute again. ");
+  }
 
+  FieldInformation fi("TetVolMesh",0,"double");   ///create output field
+  
   output = CreateField(fi);
   auto omesh = output->vmesh();
   auto ofield = output->vfield();
 
-  auto nr_of_tets  = mesh->tets.size();
-  auto nr_of_verts = mesh->verts.size();
-
   omesh->node_reserve(nr_of_verts);
   omesh->elem_reserve(nr_of_tets);
 
@@ -247,16 +272,15 @@ FieldHandle InterfaceWithCleaverAlgorithm::run(const std::vector<FieldHandle>& i
     vdata[2]=mesh->tets[i]->verts[2]->tm_v_index;
     vdata[3]=mesh->tets[i]->verts[3]->tm_v_index;
     omesh->add_elem(vdata);
-    auto mat_label = mesh->tets[i]->mat_label+1;
+    auto mat_label = mesh->tets[i]->mat_label;
     values[i]=mat_label;
   }
   ofield->resize_values();
   ofield->set_values(values);
   mesh->computeAngles();
   std::ostringstream ostr1;
-  ostr1 << "Number of tetrahedral elements:" << ofield->vmesh()->num_elems() <<  std::endl;
-  ostr1 << "Number of tetrahedral nodes:" << ofield->vmesh()->num_nodes() << std::endl;
-  ostr1 << "Number of tetrahedral nodes:" << ofield->vmesh()->num_nodes() << std::endl;
+  ostr1 << "Number of tetrahedral elements:" << nr_of_tets <<  std::endl;
+  ostr1 << "Number of tetrahedral nodes:" << nr_of_verts << std::endl;
   ostr1 << "Worst Angle (min):" <<  mesh->min_angle << std::endl;
   ostr1 << "Worst Angle (max):" <<  mesh->max_angle << std::endl;
   ostr1 << "Volume:" << volume->size().toString() << std::endl;

src/Core/Algorithms/Field/Tests/InterfaceWithCleaverAlgorithmTests.cc

@@ -42,46 +42,48 @@ using namespace SCIRun::Core::Algorithms::Fields;
 using namespace SCIRun::TestUtils;
 
 
-  FieldHandle SphereSignedDistanceField(bool negative)
+  FieldHandle BoxSignedDistanceField(bool negative)
   {  
 
    FieldInformation f("LatVolMesh", 1, "float");
-   size_type size = 5;
-   Point minb(-2.0, -2.0, -2.0);
-   Point maxb(2.0, 2.0, 2.0);
+   size_type size = 3;
+   Point minb(-1.0, -1.0, -1.0);
+   Point maxb(1.0, 1.0, 1.0);
    MeshHandle omesh = CreateMesh(f, size, size, size, minb, maxb);
    FieldHandle ofield = CreateField(f,omesh);
-   std::vector<double> values(125);
-   values[0]=0.9641;values[1]=0.5;values[2]=0.32843;values[3]=0.5;values[4]=0.9641;values[5]=0.5;
-   values[6]=-0.05051;values[7]=-0.26393;values[8]=-0.05051;values[9]=0.5;values[10]=0.32843;
-   values[11]=-0.26393;values[12]=-0.5;values[13]=-0.26393;values[14]=0.32843;values[15]=0.5;
-   values[16]=-0.05051;values[17]=-0.26393;values[18]=-0.05051;values[19]=0.5;values[20]=0.9641;
-   values[21]=0.5;values[22]=0.32843;values[23]=0.5;values[24]=0.9641;values[25]=0.5;
-   values[26]=-0.05051;values[27]=-0.26393;values[28]=-0.05051;values[29]=0.5;
-   values[30]=-0.05051;values[31]=-0.76795;values[32]=-1.0858;values[33]=-0.76795;
-   values[34]=-0.05051;values[35]=-0.26393;values[36]=-1.0858;values[37]=-1.5;
-   values[38]=-1.0858;values[39]=-0.26393;values[40]=-0.05051;values[41]=-0.76795;
-   values[42]=-1.0858;values[43]=-0.76795;values[44]=-0.05051;values[45]=0.5;values[46]=-0.05051;
-   values[47]=-0.26393;values[48]=-0.05051;values[49]=0.5;values[50]=0.32843;values[51]=-0.26393;
-   values[52]=-0.5;values[53]=-0.26393;values[54]=0.32843;values[55]=-0.26393;values[56]=-1.0858;
-   values[57]=-1.5;values[58]=-1.0858;values[59]=-0.26393;values[60]=-0.5;values[61]=-1.5;
-   values[62]=-2.5;values[63]=-1.5;values[64]=-0.5;values[65]=-0.26393;values[66]=-1.0858;
-   values[67]=-1.5;values[68]=-1.0858;values[69]=-0.26393;values[70]=0.32843;values[71]=-0.26393;
-   values[72]=-0.5;values[73]=-0.26393;values[74]=0.32843;values[75]=0.5;values[76]=-0.05051;
-   values[77]=-0.26393;values[78]=-0.05051;values[79]=0.5;values[80]=-0.05051;values[81]=-0.76795;
-   values[82]=-1.0858;values[83]=-0.76795;values[84]=-0.05051;values[85]=-0.26393;values[86]=-1.0858;
-   values[87]=-1.5;values[88]=-1.0858;values[89]=-0.26393;values[90]=-0.05051;values[91]=-0.76795;
-   values[92]=-1.0858;values[93]=-0.76795;values[94]=-0.05051;values[95]=0.5;values[96]=-0.05051;
-   values[97]=-0.26393;values[98]=-0.05051;values[99]=0.5;values[100]=0.9641;values[101]=0.5;
-   values[102]=0.32843;values[103]=0.5;values[104]=0.9641;values[105]=0.5;values[106]=-0.05051;
-   values[107]=-0.26393;values[108]=-0.05051;values[109]=0.5;values[110]=0.32843;values[111]=-0.26393;
-   values[112]=-0.5;values[113]=-0.26393;values[114]=0.32843;values[115]=0.5;values[116]=-0.05051;
-   values[117]=-0.26393;values[118]=-0.05051;values[119]=0.5;values[120]=0.9641;values[121]=0.5;
-   values[122]=0.32843;values[123]=0.5;values[124]=0.9641;  
+   std::vector<double> values(27);
+   
+   values[0]=-1.73205077648162841796875;
+   values[1]=-1.41421353816986083984375;
+   values[2]=-1.73205077648162841796875;
+   values[3]=-1.41421353816986083984375;
+   values[4]=-1;
+   values[5]=-1.41421353816986083984375;
+   values[6]=-1.73205077648162841796875;
+   values[7]=-1.41421353816986083984375;
+   values[8]=-1.73205077648162841796875;
+   values[9]=-1.41421353816986083984375;
+   values[10]=-1;
+   values[11]=-1.41421353816986083984375;
+   values[12]=-1;
+   values[13]=0;
+   values[14]=-1;
+   values[15]=-1.41421353816986083984375;
+   values[16]=-1;
+   values[17]=-1.41421353816986083984375;
+   values[18]=-1.73205077648162841796875;
+   values[19]=-1.41421353816986083984375;
+   values[20]=-1.73205077648162841796875;
+   values[21]=-1.41421353816986083984375;
+   values[22]=-1;
+   values[23]=-1.41421353816986083984375;
+   values[24]=-1.73205077648162841796875;
+   values[25]=-1.41421353816986083984375;
+   values[26]=-1.73205077648162841796875;
 
    if (negative)
    {
-    for (int i=0;i<125;i++)
+    for (int i=0;i<27;i++)
      {
        values[i]=-values[i];
      }
@@ -95,16 +97,19 @@ using namespace SCIRun::TestUtils;
 TEST(CleaverInterfaceTest, SphereSignedDistanceFieldMatrix1)
 {
   InterfaceWithCleaverAlgorithm algo;
-  FAIL() << "todo";
 
-/* auto info = algo.run(SphereSignedDistanceField(false),SphereSignedDistanceField(true));
+  std::vector<FieldHandle> inputs;
+  inputs.push_back(BoxSignedDistanceField(true));
+  inputs.push_back(BoxSignedDistanceField(false));
+  
+  auto info = algo.run(inputs);
 
- std::cout << "Number of mesh elements: " <<  info->vmesh()->num_elems() << std::endl;
- std::cout << "Number of mesh nodes: " <<  info->vmesh()->num_nodes() << std::endl;
- std::cout << "Number of mesh values: " <<  info->vfield()->num_values() << std::endl;
+  std::cout << "Number of mesh elements: " <<  info->vmesh()->num_elems() << std::endl;
+  std::cout << "Number of mesh nodes: " <<  info->vmesh()->num_nodes() << std::endl;
+  std::cout << "Number of mesh values: " <<  info->vfield()->num_values() << std::endl;
 
- ASSERT_TRUE(info->vmesh()->num_elems() != 98650);
- ASSERT_TRUE(info->vmesh()->num_nodes() != 18367);
- ASSERT_TRUE(info->vfield()->num_values() != 98650);
- */
+  ASSERT_TRUE(info->vmesh()->num_elems() == 4326);
+  ASSERT_TRUE(info->vmesh()->num_nodes() == 908);
+  ASSERT_TRUE(info->vfield()->num_values() == 4326);
+ 
 }