17_test_ann_vs_geom_trimesh.cpp

This is an example of how to use the netdem library.

#include "contact_solver_factory.hpp"
#include "distribution_uniform.hpp"
#include "general_net.hpp"
#include "igl_wrapper.hpp"
#include "mlpack_utils.hpp"
#include "model_volume_based.hpp"
#include "particle.hpp"
#include "regression_net.hpp"
#include "shape_trimesh.hpp"
#include "solver_boolean_pp.hpp"
#include "solver_boolean_pw.hpp"
#include "spherical_voronoi.hpp"
#include "utils_math.hpp"
#include <fstream>
#include <iostream>
#include <random>
#include <sstream>
#include <string>
 
using namespace netdem;
using namespace std;
 
void TestANNvsGeometricTrimesh_v2() {
  // load particle
  TriMesh tri_mesh_1;
  tri_mesh_1.InitFromSTL("data/particle_template.stl");
  tri_mesh_1.Decimate(200);
  tri_mesh_1.AlignAxes();
  tri_mesh_1.SetSize(1.0);
 
  Particle obj_p1 = Particle(&tri_mesh_1);
  obj_p1.need_update_stl_model = true;
 
  Particle obj_p2 = Particle(&tri_mesh_1);
  obj_p2.need_update_stl_model = true;
  cout << "particle created ... " << endl;
 
  obj_p2.SetRodrigues(1.2, 0, 1, 0);
  obj_p2.SetPosition(0, 0, 0.94);
 
  SolverBooleanPP cnt_solver;
  VolumeBased cnt_model;
 
  string root_dir = "local/examples/netdem/ann_models/trimesh_trimesh/";
  GeneralNet classifier;
  classifier.Load(root_dir + "ann_classifier.xml", "detection");
  RegressionNet regressor;
  regressor.Load(root_dir + "ann_regressor.xml", "resolution");
 
  // random generator
  UniformDistribution uniform_dist(0.0, 1.0);
 
  // random cases
  for (int trial = 0; trial < 100; trial++) {
    Vec3d pos;
    pos[0] = obj_p2.pos[0] + (uniform_dist.Get() - 0.5) * 0.001;
    pos[1] = obj_p2.pos[1] + (uniform_dist.Get() - 0.5) * 0.001;
    pos[2] = obj_p2.pos[2] + (uniform_dist.Get() - 0.5) * 0.001;
 
    double rot_angle{1.0e-4};
    Vec3d rot_axis;
    rot_axis[0] = uniform_dist.Get() * 2.0 - 1.0;
    rot_axis[1] = uniform_dist.Get() * 2.0 - 1.0;
    rot_axis[2] = uniform_dist.Get() * 2.0 - 1.0;
 
    auto dquat = Math::Quaternion::FromRodrigues(rot_angle, rot_axis);
    auto quat = Math::Quaternion::Multiply(dquat, obj_p2.quaternion);
    Math::Quaternion::Normalize(&quat);
 
    // apply to particle
    obj_p2.SetPosition(pos[0], pos[1], pos[2]);
    obj_p2.SetQuaternion(quat[0], quat[1], quat[2], quat[3]);
 
    // contact detection and resolution
    arma::mat input(7, 1, arma::fill::zeros);
    input(0, 0) = obj_p2.pos[0];
    input(1, 0) = obj_p2.pos[1];
    input(2, 0) = obj_p2.pos[2];
    input(3, 0) = obj_p2.quaternion[0] * Math::Sign(obj_p2.quaternion[0]);
    input(4, 0) = obj_p2.quaternion[1] * Math::Sign(obj_p2.quaternion[0]);
    input(5, 0) = obj_p2.quaternion[2] * Math::Sign(obj_p2.quaternion[0]);
    input(6, 0) = obj_p2.quaternion[3] * Math::Sign(obj_p2.quaternion[0]);
    auto output = classifier.Classify(input);
 
    cnt_solver.Init(&obj_p1, &obj_p2);
    auto cnt_flag_geo = cnt_solver.Detect();
    cout << "ann vs geometric: " << output(0) << ", " << cnt_flag_geo << endl;
 
    if (cnt_flag_geo) {
      auto cnt = ContactPP(&obj_p1, &obj_p2);
      cnt.SetCollisionModel(&cnt_model);
      cnt_solver.ResolveInit(&cnt, 1.0e-4);
      auto &cnt_geoms = cnt.collision_entries[0].cnt_geoms;
 
      // skip the contact if overlap is too large
      if (cnt_geoms.vol * cnt_geoms.sn > 6.0e-4)
        continue;
 
      auto output = regressor.Predict(input);
      cout << ">> ann: " << pow(output(0, 0) / 40.0, 2.0) << ", "
           << output(1, 0) / 40.0 << ", " << output(2, 0) << ", "
           << output(3, 0) << ", " << output(4, 0) << ", " << output(5, 0)
           << ", " << output(6, 0) << ", " << output(7, 0) << endl;
 
      cout << ">> geo: " << cnt_geoms.vol << ", " << cnt_geoms.sn << ", "
           << cnt_geoms.dir_n[0] << ", " << cnt_geoms.dir_n[1] << ", "
           << cnt_geoms.dir_n[2] << ", " << cnt_geoms.pos[0] << ", "
           << cnt_geoms.pos[1] << ", " << cnt_geoms.pos[2] << endl;
    }
  }
}