Added lim tutorial

Former-commit-id: e63f9f14e9bec7dae3918649b6773a460ac106f7

include/igl/lim/compute_lim.cpp

@@ -0,0 +1,79 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Christian Schüller <schuellchr@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#include "compute_lim.h"
+#include <LIMSolverInterface.h>
+
+IGL_INLINE int igl::compute_lim(
+  Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
+  const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
+  const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
+  const std::vector<int>& borderVertices,
+  const Eigen::Matrix<double,Eigen::Dynamic,1>& gradients,
+  const Eigen::SparseMatrix<double>& constraintMatrix,
+  const Eigen::Matrix<double,Eigen::Dynamic,1>& constraintTargets,
+  int energyType,
+  double tolerance,
+  int maxIteration,
+  bool findLocalMinima)
+{
+  return ComputeLIM(
+    vertices,
+  initialVertices,
+  elements,
+  borderVertices,
+  gradients,
+  constraintMatrix,
+  constraintTargets,
+  energyType,
+  tolerance,
+  maxIteration,
+  findLocalMinima,
+  true,
+  true,
+  -1,
+  -1);
+}
+
+IGL_INLINE int igl::compute_lim(
+  Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
+  const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
+  const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
+  const std::vector<int>& borderVertices,
+  const Eigen::Matrix<double,Eigen::Dynamic,1>& gradients,
+  const Eigen::SparseMatrix<double>& constraintMatrix,
+  const Eigen::Matrix<double,Eigen::Dynamic,1>& constraintTargets,
+  int energyType,
+  double tolerance,
+  int maxIteration,
+  bool findLocalMinima,
+  bool enableOuput,
+  bool enableAlphaUpdate,
+  double beta,
+  double eps)
+{
+  return ComputeLIM(
+    vertices,
+  initialVertices,
+  elements,
+  borderVertices,
+  gradients,
+  constraintMatrix,
+  constraintTargets,
+  energyType,
+  tolerance,
+  maxIteration,
+  findLocalMinima,
+  enableOuput,
+  enableAlphaUpdate,
+  beta,
+  eps);
+}
+
+#ifdef IGL_STATIC_LIBRARY
+// Explicit template instanciation
+#endif

include/igl/lim/compute_lim.h

@@ -0,0 +1,84 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+//
+// Copyright (C) 2014 Christian Schüller <schuellchr@gmail.com>
+//
+// This Source Code Form is subject to the terms of the Mozilla Public License
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_COMPUTE_LIM_H
+#define IGL_COMPUTE_LIM_H
+#include <igl/igl_inline.h>
+#include <Eigen/Core>
+
+namespace igl
+{
+  // Computes a locally injective mapping of a triangle or tet-mesh based on a deformation energy
+  // subject to some provided linear positional constraints Cv-d.
+  //
+  // Inputs:
+  //   vertices          vx3 matrix containing vertex position of the mesh
+  //   initialVertices   vx3 matrix containing vertex position of initial rest pose mesh
+  //   elements          exd matrix containing vertex indices of all elements
+  //   borderVertices    (only needed for 2D LSCM) vector containing indices of border vertices
+  //   gradients         (only needed for 2D Poisson) vector containing partial derivatives of target element gradients (structure is: [xx_0, xy_0, xx_1, xy_1, ..., xx_v, xy_v, yx_0, yy_0, yx_1, yy_1, ..., yx_v, yy_v]')
+  //   constraintMatrix  C: (c)x(v*(d-1)) sparse linear positional constraint matrix. X an Y-coordinates are alternatingly stacked per row (structure for triangles: [x_1, y_1, x_2, y_2, ..., x_v,y_v])
+  //   constraintTargets d: c vector target positions
+  //   energyType        type of used energy: 0=Dirichlet,1=Laplacian,2=Green,3=ARAP,4=LSCM
+  //   tolerance         max squared positional constraints error
+  //   maxIteration      max number of iterations
+  //   findLocalMinima   iterating until a local minima is found. If not enabled only tolerance must be fulfilled.
+  //   enableOutput      (optional) enables the output (#itaration / hessian correction / step size / positional constraints / barrier constraints / deformation energy) (default : true)
+  //   enableAlphaUpdate (optional) enables dynamic alpha weight adjustment (default = true)
+  //   beta              (optional) steepness factor of barrier slopes (default: ARAP/LSCM = 0.01, Green = 1)
+  //   eps               (optional) smallest valid triangle area (default: 1e-5 * smallest triangle)
+  //   
+  // where:
+  //   v : # vertices
+  //   c : # linear constraints
+  //   e : # elements of mesh
+  //   d : # vertices per element (triangle = 3, tet = 4)
+  //----------------------------------------------------------------------------------------
+  // Output:
+  // vertices          vx3 matrix containing resulting vertex position of the mesh
+  //----------------------------------------------------------------------------------------
+  // Return values:
+  //  1 : Successful optimization with fulfilled tolerance
+  // -1 : Max iteration reached before tolerance was fulfilled
+  // -2 : not feasible -> has inverted elements (may want to decrease eps?)
+
+  int compute_lim(
+    Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
+    const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
+    const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
+    const std::vector<int>& borderVertices,
+    const Eigen::Matrix<double,Eigen::Dynamic,1>& gradients,
+    const Eigen::SparseMatrix<double>& constraintMatrix,
+    const Eigen::Matrix<double,Eigen::Dynamic,1>& constraintTargets,
+    int energyType,
+    double tolerance,
+    int maxIteration,
+    bool findLocalMinima);
+
+  int compute_lim(
+    Eigen::Matrix<double,Eigen::Dynamic,3>& vertices,
+    const Eigen::Matrix<double,Eigen::Dynamic,3>& initialVertices,
+    const Eigen::Matrix<int,Eigen::Dynamic,Eigen::Dynamic>& elements,
+    const std::vector<int>& borderVertices,
+    const Eigen::Matrix<double,Eigen::Dynamic,1>& gradients,
+    const Eigen::SparseMatrix<double>& constraintMatrix,
+    const Eigen::Matrix<double,Eigen::Dynamic,1>& constraintTargets,
+    int energyType,
+    double tolerance,
+    int maxIteration,
+    bool findLocalMinima,
+    bool enableOuput,
+    bool enableAlphaUpdate,
+    double beta,
+    double eps);
+}
+
+#ifndef IGL_STATIC_LIBRARY
+#  include "compute_lim.cpp"
+#endif
+
+#endif

include/igl/linsolve_with_fixed.h

@@ -0,0 +1,109 @@
+// This file is part of libigl, a simple c++ geometry processing library.
+// 
+// Copyright (C) 2013 Alec Jacobson <alecjacobson@gmail.com>
+// 
+// This Source Code Form is subject to the terms of the Mozilla Public License 
+// v. 2.0. If a copy of the MPL was not distributed with this file, You can 
+// obtain one at http://mozilla.org/MPL/2.0/.
+#ifndef IGL_LINSOLVE_WITH_FIXED_H
+#define IGL_LINSOLVE_WITH_FIXED_H
+#include "igl_inline.h"
+
+#define EIGEN_YES_I_KNOW_SPARSE_MODULE_IS_NOT_STABLE_YET
+#include <Eigen/Core>
+#include <Eigen/Dense>
+#include <Eigen/Sparse>
+#include <Eigen/SparseExtra>
+// Bug in unsupported/Eigen/SparseExtra needs iostream first
+#include <iostream>
+#include <unsupported/Eigen/SparseExtra>
+
+namespace igl
+{
+  // LINSOLVE_WITH_FIXED Solves a sparse linear system Ax=b with
+  // fixed known values in x
+
+  // Templates:
+  //   T  should be a eigen matrix primitive type like float or double
+  // Inputs:
+  //   A  m by n sparse matrix
+  //   b  n by k dense matrix, k is bigger than 1 for multiple right handsides
+  //   known list of indices to known rows in x (must be sorted in ascending order)
+  //   Y  list of fixed values corresponding to known rows in Z
+  // Outputs:
+  //   x  n by k dense matrix containing the solution
+
+  template<typename T, typename Derived>
+  inline void SolveLinearSystemWithBoundaryConditions(const SparseMatrix<T>& A, const MatrixBase<Derived>& b, const vector<int>& known, MatrixBase<Derived>& x)
+  {
+  const int rows = x.rows();
+  const int cols = x.cols();
+  SparseMatrix<T> AP(rows,rows);
+
+  const int knownCount = known.size();
+  const int unknownCount = rows - knownCount;
+  int knownIndex = 0;
+  int unknownIndex = 0;
+
+  if(knownCount >= rows)
+  {
+  std::cerr << "No unknowns to solve for!\n";
+  return;
+  }
+
+  // Permutate to sort by known boundary conditions -> [A1,A2] * [x1;x2]
+  VectorXi transpositions;
+  transpositions.resize(rows);
+  for(int n=0;n<rows;n++)
+  {
+    if(knownIndex < known.size() && n == known[knownIndex])
+    {
+      transpositions[unknownCount + knownIndex] = n;
+      knownIndex++;
+    }
+    else
+    {
+      transpositions[unknownIndex] = n;
+      unknownIndex++;
+    }
+  }
+  PermutationMatrix<Dynamic,Dynamic,int> P(transpositions);
+  PermutationMatrix<Dynamic,Dynamic,int> PInv = P.inverse();
+  A.twistedBy(PInv).evalTo(AP);
+
+  // Split in kown and unknown parts A1,A2,x2,b1
+  SparseMatrix<T> A1, A2;
+  internal::plain_matrix_type<Derived>::type x1, x2, b1;
+  x2 = (PInv*x).block(unknownCount,0,knownCount,cols);
+  b1 = (PInv*b).block(0,0,unknownCount,cols);
+
+  SparseMatrix<T> A1Temp(rows,unknownCount);
+  SparseMatrix<T> A2Temp(rows,knownCount);
+  A1Temp = AP.innerVectors(0,unknownCount).transpose();
+  A2Temp = AP.innerVectors(unknownCount,knownCount).transpose();
+  A1 = A1Temp.innerVectors(0,unknownCount).transpose();
+  A2 = A2Temp.innerVectors(0,unknownCount).transpose();
+
+  // Solve A1*x1 = b - A2*x2
+  SparseLU<SparseMatrix<T>> solver;
+  solver.compute(A1);
+  if (solver.info() != Eigen::Success)
+  {
+  std::cerr << "Matrix can not be decomposed.\n";
+  return;
+  }
+
+  internal::plain_matrix_type<Derived>::type t = b1-A2*x2;
+  x1 = solver.solve(t);
+  if (solver.info() != Eigen::Success)
+  {
+  std::cerr << "Solving failed.\n";
+  return;
+  }
+
+  // Compose resulting x
+  x.block(0,0,unknownCount,cols) = x1;
+  x.block(unknownCount,0,knownCount,cols) = x2;
+  x = P * x;
+  };
+}

tutorial/607_LIM/CMakeLists.txt

@@ -0,0 +1,14 @@
+cmake_minimum_required(VERSION 2.6)
+project(607_LIM)
+
+include("../CMakeLists.shared")
+find_package(LIM REQUIRED)
+
+include_directories(${LIM_INCLUDE_DIR})
+
+set(SOURCES
+${PROJECT_SOURCE_DIR}/main.cpp
+)
+
+add_executable(${PROJECT_NAME} ${SOURCES} ${SHARED_SOURCES} ${LIM_SOURCES})
+target_link_libraries(${PROJECT_NAME} ${SHARED_LIBRARIES})

tutorial/607_LIM/main.cpp

@@ -0,0 +1,97 @@
+#include <igl/readOFF.h>
+#include <igl/viewer/Viewer.h>
+#include <igl/per_vertex_normals.h>
+#include <igl/avg_edge_length.h>
+#include <iostream>
+
+#include "igl/lim/compute_lim.h"
+
+using namespace Eigen;
+using namespace std;
+
+// Mesh
+Eigen::MatrixX3d V0;
+Eigen::MatrixX3d V1;
+Eigen::MatrixXi F;
+
+SparseMatrix<double> C;
+Eigen::VectorXd b;
+vector<int> bv; // border vertices (don't need to be defined exept for LSCM)
+Eigen::VectorXd g; // gradients (don't need to be defined execpt for Poisson)
+
+// This function is called every time a keyboard button is pressed
+// keys: 0:Original Mesh / 1:Harmonic / 2:Biharmonic / 3:Green / 4:ARAP
+bool key_down(igl::Viewer& viewer,unsigned char key,int modifier)
+{
+  using namespace std;
+  using namespace Eigen;
+
+  if(key >= '0' && key <= '4')
+  {
+    // compute locally injective map
+    int energy = key - 49;
+
+    V1 = V0;
+    if(energy >= 0)
+      igl::compute_lim(V1,V0,F,bv,g,C,b,energy,1e-8,100,true);
+
+    // set mesh
+    viewer.set_mesh(V1,F);
+  }
+
+
+  return false;
+}
+
+int main(int argc, char *argv[])
+{
+  using namespace std;
+  using namespace Eigen;
+
+  // load a mesh in OFF format
+  igl::readOFF("../shared/grid.off",V0,F);
+  V1 = V0;
+
+  // find bottom and left boundary vertices
+  vector<int> fixedVertices;
+  for(int i=0;i<V0.rows();i++)
+  {
+    if(V0.row(i)[0] == 0 || V0.row(i)[1] == 0)
+      fixedVertices.push_back(i);
+  }
+
+  // fix boundaries
+  C.resize(2*fixedVertices.size()+2,V0.rows()*2);
+  for(int i=0;i<fixedVertices.size();i++)
+  {
+    C.insert(2*i,2*fixedVertices[i]) = 1;
+    C.insert(2*i+1,2*fixedVertices[i]+1) = 1;
+  }
+
+  // constraint targets
+  b.resize(2*fixedVertices.size()+2);
+ 
+  for(int i=0;i<fixedVertices.size();i++)
+  {
+    b(2*i) = V0.row(fixedVertices[i])[0];
+    b(2*i+1) = V0.row(fixedVertices[i])[1];
+  }
+
+  // drag top left corner vertex to the center
+  int id = 2;
+  C.insert(2*fixedVertices.size(),2*id) = 1;
+  C.insert(2*fixedVertices.size()+1,2*id+1) = 1;
+  b(2*fixedVertices.size()) = 0.2;
+  b(2*fixedVertices.size()+1) = 0.2;
+
+  // compute locally injective map
+  igl::compute_lim(V1,V0,F,bv,g,C,b,0,1e-8,100,true);
+
+  // Show mesh
+  igl::Viewer viewer;
+  viewer.callback_key_down = &key_down;
+  viewer.set_mesh(V1, F);
+  viewer.options.show_lines = true;
+  viewer.options.lighting_factor = 0.0f;
+  viewer.launch();
+}

tutorial/CMakeLists.txt

@@ -37,3 +37,5 @@ add_subdirectory("601_Serialization")
 add_subdirectory("604_Triangle")
 add_subdirectory("605_Tetgen")
 add_subdirectory("606_AmbientOcclusion")
+
+add_subdirectory("607_LIM")

tutorial/cmake/FindLIM.cmake

@@ -0,0 +1,47 @@
+# - Try to find the LIM library
+# Once done this will define
+#
+#  LIM_FOUND - system has LIM
+#  LIM_INCLUDE_DIR - the LIM include directory
+#  LIM_SOURCES - the LIM source files
+
+FIND_PATH(LIM_INCLUDE_DIR LIMSolverInterface.h
+   /usr/include
+   /usr/local/include
+   ${PROJECT_SOURCE_DIR}/../libigl/external/lim/
+   ${PROJECT_SOURCE_DIR}/../../external/lim/
+   NO_DEFAULT_PATH
+)
+
+set(
+  LIM_SOURCES
+  ${LIM_INCLUDE_DIR}/NMSolver.cpp
+  ${LIM_INCLUDE_DIR}/LIMSolver.cpp
+  ${LIM_INCLUDE_DIR}/LIMSolver2D.cpp
+  ${LIM_INCLUDE_DIR}/LIMSolver3D.cpp
+  ${LIM_INCLUDE_DIR}/TriangleMesh.cpp
+  ${LIM_INCLUDE_DIR}/TetrahedronMesh.cpp
+  ${LIM_INCLUDE_DIR}/Dirichlet_LIMSolver2D.cpp
+  ${LIM_INCLUDE_DIR}/Dirichlet_LIMSolver3D.cpp
+  ${LIM_INCLUDE_DIR}/UniformLaplacian_LIMSolver2D.cpp
+  ${LIM_INCLUDE_DIR}/UniformLaplacian_LIMSolver3D.cpp
+  ${LIM_INCLUDE_DIR}/Laplacian_LIMSolver2D.cpp
+  ${LIM_INCLUDE_DIR}/Laplacian_LIMSolver3D.cpp
+  ${LIM_INCLUDE_DIR}/LGARAP_LIMSolver2D.cpp
+  ${LIM_INCLUDE_DIR}/LGARAP_LIMSolver3D.cpp
+  ${LIM_INCLUDE_DIR}/GreenStrain_LIMSolver2D.cpp
+  ${LIM_INCLUDE_DIR}/GreenStrain_LIMSolver3D.cpp
+  ${LIM_INCLUDE_DIR}/LSConformal_LIMSolver2D.cpp
+  ${LIM_INCLUDE_DIR}/Poisson_LIMSolver2D.cpp
+  )
+
+if(LIM_INCLUDE_DIR)
+   message(STATUS "Found LIM: ${LIM_INCLUDE_DIR}")
+else(LIM_INCLUDE_DIR)
+   message(FATAL_ERROR "could NOT find LIM")
+endif(LIM_INCLUDE_DIR)
+
+set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DANSI_DECLARATORS")
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DANSI_DECLARATORS")
+
+MARK_AS_ADVANCED(LIM_INCLUDE_DIR LIM_LIBRARIES LIM_SOURCES)

tutorial/shared/grid.off

@@ -0,0 +1,403 @@
+OFF
+145 256 0
+0 0 0 
+1 0 0 
+1 1 0 
+0 1 0 
+0.5 0.5 0 
+0.5 0 0 
+0.75 0.25 0 
+0.25 0.25 0 
+1 0.5 0 
+0.75 0.75 0 
+0.5 1 0 
+0.25 0.75 0 
+0 0.5 0 
+0.625 0.125 0 
+0.5 0.25 0 
+0.375 0.125 0 
+0.875 0.625 0 
+0.75 0.5 0 
+0.875 0.375 0 
+0.375 0.875 0 
+0.5 0.75 0 
+0.625 0.875 0 
+0.125 0.375 0 
+0.25 0.5 0 
+0.125 0.625 0 
+0.25 0 0 
+0.125 0.125 0 
+0.75 0 0 
+0.875 0.125 0 
+0.625 0.375 0 
+0.375 0.375 0 
+1 0.25 0 
+1 0.75 0 
+0.875 0.875 0 
+0.625 0.625 0 
+0.75 1 0 
+0.25 1 0 
+0.125 0.875 0 
+0.375 0.625 0 
+0 0.75 0 
+0 0.25 0 
+0.5625 0.1875 0 
+0.4375 0.1875 0 
+0.5 0.125 0 
+0.8125 0.5625 0 
+0.8125 0.4375 0 
+0.875 0.5 0 
+0.4375 0.8125 0 
+0.5625 0.8125 0 
+0.5 0.875 0 
+0.1875 0.4375 0 
+0.1875 0.5625 0 
+0.125 0.5 0 
+0.3125 0.0625 0 
+0.25 0.125 0 
+0.1875 0.0625 0 
+0.8125 0.0625 0 
+0.75 0.125 0 
+0.6875 0.0625 0 
+0.5625 0.3125 0 
+0.5 0.375 0 
+0.4375 0.3125 0 
+0.9375 0.3125 0 
+0.875 0.25 0 
+0.9375 0.1875 0 
+0.9375 0.8125 0 
+0.875 0.75 0 
+0.9375 0.6875 0 
+0.6875 0.5625 0 
+0.625 0.5 0 
+0.6875 0.4375 0 
+0.6875 0.9375 0 
+0.75 0.875 0 
+0.8125 0.9375 0 
+0.1875 0.9375 0 
+0.25 0.875 0 
+0.3125 0.9375 0 
+0.4375 0.6875 0 
+0.5 0.625 0 
+0.5625 0.6875 0 
+0.0625 0.6875 0 
+0.125 0.75 0 
+0.0625 0.8125 0 
+0.0625 0.1875 0 
+0.125 0.25 0 
+0.0625 0.3125 0 
+0.3125 0.4375 0 
+0.375 0.5 0 
+0.3125 0.5625 0 
+0.5625 0.0625 0 
+0.4375 0.0625 0 
+0.6875 0.1875 0 
+0.625 0.25 0 
+0.375 0.25 0 
+0.3125 0.1875 0 
+0.9375 0.5625 0 
+0.9375 0.4375 0 
+0.8125 0.6875 0 
+0.75 0.625 0 
+0.75 0.375 0 
+0.8125 0.3125 0 
+0.4375 0.9375 0 
+0.5625 0.9375 0 
+0.3125 0.8125 0 
+0.375 0.75 0 
+0.625 0.75 0 
+0.6875 0.8125 0 
+0.0625 0.4375 0 
+0.0625 0.5625 0 
+0.1875 0.3125 0 
+0.25 0.375 0 
+0.25 0.625 0 
+0.1875 0.6875 0 
+0.125 0 0 
+0.0625 0.0625 0 
+0.375 0 0 
+0.1875 0.1875 0 
+0.625 0 0 
+0.875 0 0 
+0.9375 0.0625 0 
+0.8125 0.1875 0 
+0.3125 0.3125 0 
+0.6875 0.3125 0 
+0.5625 0.4375 0 
+0.4375 0.4375 0 
+1 0.125 0 
+1 0.375 0 
+1 0.625 0 
+1 0.875 0 
+0.9375 0.9375 0 
+0.8125 0.8125 0 
+0.6875 0.6875 0 
+0.5625 0.5625 0 
+0.875 1 0 
+0.625 1 0 
+0.375 1 0 
+0.125 1 0 
+0.0625 0.9375 0 
+0.1875 0.8125 0 
+0.3125 0.6875 0 
+0.4375 0.5625 0 
+0 0.875 0 
+0 0.625 0 
+0 0.375 0 
+0 0.125 0 
+3 41 42 43
+3 44 45 46
+3 47 48 49
+3 50 51 52
+3 53 54 55
+3 56 57 58
+3 59 60 61
+3 62 63 64
+3 65 66 67
+3 68 69 70
+3 71 72 73
+3 74 75 76
+3 77 78 79
+3 80 81 82
+3 83 84 85
+3 86 87 88
+3 89 43 90
+3 91 92 41
+3 42 93 94
+3 95 46 96
+3 97 98 44
+3 45 99 100
+3 101 49 102
+3 103 104 47
+3 48 105 106
+3 107 52 108
+3 109 110 50
+3 51 111 112
+3 113 55 114
+3 115 90 53
+3 54 94 116
+3 117 58 89
+3 118 119 56
+3 57 120 91
+3 93 61 121
+3 92 122 59
+3 60 123 124
+3 125 64 119
+3 126 96 62
+3 63 100 120
+3 127 67 95
+3 128 129 65
+3 66 130 97
+3 99 70 122
+3 98 131 68
+3 69 132 123
+3 133 73 129
+3 134 102 71
+3 72 106 130
+3 135 76 101
+3 136 137 74
+3 75 138 103
+3 105 79 131
+3 104 139 77
+3 78 140 132
+3 141 82 137
+3 142 108 80
+3 81 112 138
+3 143 85 107
+3 144 114 83
+3 84 116 109
+3 111 88 139
+3 110 121 86
+3 87 124 140
+3 13 41 43
+3 41 14 42
+3 43 42 15
+3 16 44 46
+3 44 17 45
+3 46 45 18
+3 19 47 49
+3 47 20 48
+3 49 48 21
+3 22 50 52
+3 50 23 51
+3 52 51 24
+3 25 53 55
+3 53 15 54
+3 55 54 26
+3 27 56 58
+3 56 28 57
+3 58 57 13
+3 14 59 61
+3 59 29 60
+3 61 60 30
+3 31 62 64
+3 62 18 63
+3 64 63 28
+3 32 65 67
+3 65 33 66
+3 67 66 16
+3 17 68 70
+3 68 34 69
+3 70 69 29
+3 35 71 73
+3 71 21 72
+3 73 72 33
+3 36 74 76
+3 74 37 75
+3 76 75 19
+3 20 77 79
+3 77 38 78
+3 79 78 34
+3 39 80 82
+3 80 24 81
+3 82 81 37
+3 40 83 85
+3 83 26 84
+3 85 84 22
+3 23 86 88
+3 86 30 87
+3 88 87 38
+3 5 89 90
+3 89 13 43
+3 90 43 15
+3 13 91 41
+3 91 6 92
+3 41 92 14
+3 15 42 94
+3 42 14 93
+3 94 93 7
+3 8 95 96
+3 95 16 46
+3 96 46 18
+3 16 97 44
+3 97 9 98
+3 44 98 17
+3 18 45 100
+3 45 17 99
+3 100 99 6
+3 10 101 102
+3 101 19 49
+3 102 49 21
+3 19 103 47
+3 103 11 104
+3 47 104 20
+3 21 48 106
+3 48 20 105
+3 106 105 9
+3 12 107 108
+3 107 22 52
+3 108 52 24
+3 22 109 50
+3 109 7 110
+3 50 110 23
+3 24 51 112
+3 51 23 111
+3 112 111 11
+3 0 113 114
+3 113 25 55
+3 114 55 26
+3 25 115 53
+3 115 5 90
+3 53 90 15
+3 26 54 116
+3 54 15 94
+3 116 94 7
+3 5 117 89
+3 117 27 58
+3 89 58 13
+3 27 118 56
+3 118 1 119
+3 56 119 28
+3 13 57 91
+3 57 28 120
+3 91 120 6
+3 7 93 121
+3 93 14 61
+3 121 61 30
+3 14 92 59
+3 92 6 122
+3 59 122 29
+3 30 60 124
+3 60 29 123
+3 124 123 4
+3 1 125 119
+3 125 31 64
+3 119 64 28
+3 31 126 62
+3 126 8 96
+3 62 96 18
+3 28 63 120
+3 63 18 100
+3 120 100 6
+3 8 127 95
+3 127 32 67
+3 95 67 16
+3 32 128 65
+3 128 2 129
+3 65 129 33
+3 16 66 97
+3 66 33 130
+3 97 130 9
+3 6 99 122
+3 99 17 70
+3 122 70 29
+3 17 98 68
+3 98 9 131
+3 68 131 34
+3 29 69 123
+3 69 34 132
+3 123 132 4
+3 2 133 129
+3 133 35 73
+3 129 73 33
+3 35 134 71
+3 134 10 102
+3 71 102 21
+3 33 72 130
+3 72 21 106
+3 130 106 9
+3 10 135 101
+3 135 36 76
+3 101 76 19
+3 36 136 74
+3 136 3 137
+3 74 137 37
+3 19 75 103
+3 75 37 138
+3 103 138 11
+3 9 105 131
+3 105 20 79
+3 131 79 34
+3 20 104 77
+3 104 11 139
+3 77 139 38
+3 34 78 132
+3 78 38 140
+3 132 140 4
+3 3 141 137
+3 141 39 82
+3 137 82 37
+3 39 142 80
+3 142 12 108
+3 80 108 24
+3 37 81 138
+3 81 24 112
+3 138 112 11
+3 12 143 107
+3 143 40 85
+3 107 85 22
+3 40 144 83
+3 144 0 114
+3 83 114 26
+3 22 84 109
+3 84 26 116
+3 109 116 7
+3 11 111 139
+3 111 23 88
+3 139 88 38
+3 23 110 86
+3 110 7 121
+3 86 121 30
+3 38 87 140
+3 87 30 124
+3 140 124 4