Home Explore Help
Sign In
3DFaces
/
libigl_Martin
3
0
Fork 0
Files Issues 0 Pull Requests 0 Wiki
Browse Source

added svd using blas

Former-commit-id: 0e3e52132593a8ac700c86c3e4f18f2d2edf7ed3
Alec Jacobson (jalec 12 years ago
parent
43cd86b468
commit
85581d0f96
9 changed files with 297 additions and 0 deletions
Split View Show Diff Stats
  1. 4 0
      examples/bbw/README
  2. 24 0
      examples/svd/Makefile
  3. 96 0
      examples/svd/example.cpp
  4. 25 0
      file-formats/bf.html
  5. 39 0
      file-formats/dmat.html
  6. 17 0
      file-formats/index.html
  7. 32 0
      file-formats/tgf.html
  8. 38 0
      include/igl/svd.cpp
  9. 22 0
      include/igl/svd.h

+ 4 - 0
examples/bbw/README
View File 

@@ -23,6 +23,10 @@ This will produce at least 2 files:
 
   examples/brick-volume.mesh  Tetrahedral mesh over which weights were computed
 
   examples/brick-volume.dmat  Coefficients of weights defined over .mesh
 
 
+= File Formats =
 
+
 
+See file-formats/index.html
 
+
 
 The dependencies are:
 
   igl_lib
 
     eigen3

+ 24 - 0
examples/svd/Makefile
View File 

@@ -0,0 +1,24 @@
 
+
 
+.PHONY: all
 
+
 
+# Shared flags etc.
 
+include ../../Makefile.conf
 
+
 
+all: example
 
+
 
+.PHONY: example
 
+
 
+igl_lib=../../
 
+
 
+CFLAGS=-g
 
+inc=-I$(igl_lib)/include
 
+lib=-L$(igl_lib)/lib -ligl -framework Accelerate
 
+
 
+example: example.o
 
+	g++ $(CFLAGS) -o example example.o $(lib)
 
+
 
+example.o: example.cpp
 
+	g++ $(CFLAGS) -c example.cpp -o example.o $(inc)
 
+clean:
 
+	rm -f example.o
 
+	rm -f example

+ 96 - 0
examples/svd/example.cpp
View File 

@@ -0,0 +1,96 @@
 
+#include <igl/svd.h>
 
+#include <cstdlib>
 
+#include <Accelerate/Accelerate.h>
 
+#include <cstdio>
 
+
 
+/* Auxiliary routines prototypes */
 
+extern void print_matrix( char* desc, int m, int n, double* a, int lda );
 
+
 
+/* Parameters */
 
+
 
+void print3x3(const char * s, double * a)
 
+{
 
+  printf("%s =\n",s);
 
+  for(int i = 0;i<3;i++)
 
+  {
 
+    for(int j = 0;j<3;j++)
 
+    {
 
+      printf("%g ",a[j*3+i]);
 
+    }
 
+    printf("\n");
 
+  }
 
+  printf("\n");
 
+}
 
+
 
+int main(int argc, char * argv[])
 
+{
 
+  //// List of rest positions
 
+  ////        (0,1)
 
+  ////         / \
 
+  ////        /   \
 
+  ////       /     \
 
+  ////      /       \
 
+  ////  (-1,0)-----(1,0)
 
+  ////
 
+  //double rest[3][3] = {
 
+  //  {-1,0,0},
 
+  //  {1,0,0},
 
+  //  {0,1,0}};
 
+  //// List of pose positions
 
+  //// 
+  //// (0,1)
 
+  ////  |   \
 
+  ////  |    \
 
+  ////  |     (1,0)
 
+  ////  |    /
 
+  ////  |   /
 
+  //// (0,-1)
 
+  //double pose[3][3] = {
 
+  //  {0,1,0},
 
+  //  {0,-1,0},
 
+  //  {1,0,0}};
 
+  //// Compute covariance matrix C
 
+  //double C[3*3];
 
+  //// Initialize to zero
 
+  //for(int i = 0;i<3*3;i++)
 
+  //{
 
+  //  C[i] = 0;
 
+  //}
 
+
 
+  //// Loop over vertices
 
+  //for(int i = 0;i<3;i++)
 
+  //{
 
+  //  // Compute outer product rest[i] * pose[i]
 
+  //  // Loop over coordinates
 
+  //  for(int j = 0;j<3;j++)
 
+  //  {
 
+  //    // Loop over coordinates
 
+  //    for(int k = 0;k<3;k++)
 
+  //    {
 
+  //      C[k*3+j] = rest[i][j] * pose[i][k];
 
+  //    }
 
+  //  }
 
+  //}
 
+  //print3x3("C",C);
 
+
 
+
 
+  //
 
+  //double C[3*3] = {8,3,4,1,5,9,6,7,2};
 
+  double C[3*3] = {5242.55,3364,-0,-8170.15,-5242.56,0,-0,-0,0};
 
+  double u[3*3],s[3],vt[3*3];
 
+  print3x3("C",C);
 
+  // Compute SVD of C
 
+  igl::svd3x3(C,u,s,vt);
 
+  print3x3("u",u);
 
+  print3x3("vt",vt);
 
+
 
+  // Compute R = u*vt
 
+  double R[3*3];
 
+  const double _3 = 3;
 
+  const double _1 = 1;
 
+  cblas_dgemm(CblasColMajor, CblasNoTrans,CblasNoTrans,3,3,3,1,u,3,vt,3,1,R,3);
 
+  print3x3("RT (transposed to be row-major)",R);
 
+
 
+
 
+  return 0;
 
+}

+ 25 - 0
file-formats/bf.html
View File 

@@ -0,0 +1,25 @@
 
+<!DOCTYPE HTML>
 
+<html>
 
+  <head>
 
+    <link rel='stylesheet' type='text/css' href='../style.css' >
 
+    <title>IGL_LIB file formats | .bf</title>
 
+  </head>
 
+  <body>
 
+    <h1>.bf - bone forests</h1>
 
+    <hr>
 
+    <p>
 
+A .bf file contains a bone forest. That is possibly multiple bone trees, or
 
+skeletons. Singleton trees will be interpreted as point handles. Otherwise
 
+each edge in a tree will be interpreted as a bone segment.
 
+    </p>
 
+    <p>
 
+Each line contains data about a vertex (joint) of a bone tree:
 
+    </p>
 
+    <pre><code>[weight index] [parent index] [x] [y] [z] [undocument optional data]</code></pre>
 
+    <p>
 
+Indices begin with 0.  The weight index is -1 if the bone does not have an
 
+associated weight. The parent index is -1 for root nodes.
 
+    </p>
 
+  </body>
 
+</html>
 
+

+ 39 - 0
file-formats/dmat.html
View File 

@@ -0,0 +1,39 @@
 
+<!DOCTYPE HTML>
 
+<html>
 
+  <head>
 
+    <link rel='stylesheet' type='text/css' href='../style.css' >
 
+    <title>IGL_LIB file formats | .dmat</title>
 
+  </head>
 
+  <body>
 
+    <h1>.dmat - dense matrices</h1>
 
+    <hr>
 
+    <p>
 
+A .dmat file contains a dense matrix in column major order. It can contain
 
+ASCII or binary data. Note that it is uncompressed so binary only reduces the
 
+file size by 50%. But writing and reading binary is usualy faster. In MATLAB,
 
+binary is almost 100x faster.
 
+    </p>
 
+    <h2>ASCII</h2>
 
+    <p>
 
+The first line is a header containing:
 
+    </p>
 
+    <pre><code>[#cols] [#rows]</code></pre>
 
+    <p>
 
+Then the coefficients are printed in column-major order separated by spaces. 
+    </p>
 
+    <h2>Binary</h2>
 
+    <p>
 
+Binary files will also contain the ascii header, but it should read:
 
+    </p>
 
+    <pre><code>0 0</code></pre>
 
+    <p>
 
+Then there should be another header containing the size of the binary part:
 
+    </p>
 
+    <pre><code>[#cols] [#rows]</code></pre>
 
+    <p>
 
+Then coefficents are written in column-major order in Little-endian 4-byte
 
+double precision IEEE floating point format.
 
+    </p>
 
+  </body>
 
+</html>
 
+

+ 17 - 0
file-formats/index.html
View File 

@@ -0,0 +1,17 @@
 
+<!DOCTYPE HTML>
 
+<html>
 
+  <head>
 
+    <link rel='stylesheet' type='text/css' href='../style.css' >
 
+    <title>IGL_LIB file formats</title>
 
+  </head>
 
+  <body>
 
+  <h1>IGL_LIB file formats</h1>
 
+  <ul>
 
+    <li><a href="./dmat.html">.dmat</a> uncompressed ASCII/binary files for dense matrices</li>
 
+    <li><a href="./bf.html">.bf</a> ASCII files for representing skeletal bone "forests"</li>
 
+    <li><a href="./tgf.html">.tgf</a> ASCII files for representing control handle graphs</li>
 
+    <li><a href="http://tetgen.berlios.de/fformats.off.html">.off</a> Geomview's polyhedral file format</li>
 
+    <li><a href="http://www.ann.jussieu.fr/frey/publications/RT-0253.pdf#page=33">.mesh</a> Medit's triangle surface mesh + tetrahedral volume mesh file format, see page 33, section 7.2.1</li>
 
+  </ul>
 
+  </body>
 
+</html>

+ 32 - 0
file-formats/tgf.html
View File 

@@ -0,0 +1,32 @@
 
+<!DOCTYPE HTML>
 
+<html>
 
+  <head>
 
+    <link rel='stylesheet' type='text/css' href='../style.css' >
 
+    <title>IGL_LIB file formats | .tgf</title>
 
+  </head>
 
+  <body>
 
+    <h1>.tgf - control handle graphs</h1>
 
+    <hr>
 
+    <p>
 
+A .tgf file contains a graph of describing a set of control handles/structures:
 
+point controls, bones of a skelton and cages made of "cage edges".
 
+    </p>
 
+    <p>
 
+The first part of the file consists of lines regarding each vertex of the graph. Each line reads:
 
+    </p>
 
+    <pre><code>[index] [x] [y] [z] [undocument optional data]</code></pre>
 
+    <p>
 
+Indices begin with 1 and should proceed in order.
 
+Then there should be a line with a sole:
 
+    </p>
 
+    <pre><code>#</code></pre>
 
+    <p>
 
+The next section concerns the edges of the graph. Each line corresponds to an edge:
 
+    </p>
 
+    <pre><code>[source index] [dest index] [is bone] [is pseudo-edge] [is cage edge] [undocument other data]</code></pre>
 
+    <p>
 
+Bone edges trump pseudo and cage edges. 
+    </p>
 
+  </body>
 
+</html>
 
+

+ 38 - 0
include/igl/svd.cpp
View File 

@@ -0,0 +1,38 @@
 
+#include "svd.h"
 
+
 
+#include <cstdlib>
 
+#include <Accelerate/Accelerate.h>
 
+#include <cstdio>
 
+
 
+bool igl::svd3x3(double * a, double * u, double * s, double * vt)
 
+{
 
+  /* Locals */
 
+  int m = 3, n = 3, lda = 3, ldu = 3, ldvt = 3, info, lwork;
 
+  double wkopt;
 
+  double* work;
 
+  /* Local arrays */
 
+  /* iwork dimension should be at least 8*min(m,n) */
 
+  int iwork[8*3];
 
+  //double s[3], u[3*3], vt[3*3];
 
+  //double a[3*3] = {8,3,4,1,5,9,6,7,2};
 
+  /* Query and allocate the optimal workspace */
 
+  lwork = -1;
 
+  dgesdd_( 
+    "Singular vectors",
 
+    &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &wkopt, &lwork, iwork, &info);
 
+  lwork = (int)wkopt;
 
+  work = (double*)malloc( lwork*sizeof(double) );
 
+  /* Compute SVD */
 
+  dgesdd_(
 
+    "Singular vectors",
 
+    &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, iwork, &info );
 
+  /* Check for convergence */
 
+  if( info > 0 )
 
+  {
 
+    printf("The algorithm computing SVD failed to converge.\n" );
 
+    return false;
 
+  }
 
+  /* Free workspace */
 
+  free( (void*)work );
 
+  return true;
 
+}

+ 22 - 0
include/igl/svd.h
View File 

@@ -0,0 +1,22 @@
 
+#ifndef IGL_SVD_H
 
+#define IGL_SVD_H
 
+namespace igl
 
+{
 
+  // Compute 3x3 SVD using lapack's dgesdd_ function
 
+  //
 
+  // Input:
 
+  //   a  pointer to 3x3 matrix in COLUMN MAJOR order
 
+  // Outputs:
 
+  //   u  pointer to 3x3 matrix in COLUMN MAJOR order
 
+  //   s  pointer to 3 vector 
+  //   vt  pointer to 3x3 matrix in COLUMN MAJOR order, or think of this as v in
 
+  //     row-major order
 
+  // Returns true on success, false on failure
 
+  bool svd3x3(double * a, double * u, double * s, double * vt);
 
+};
 
+
 
+#ifdef IGL_HEADER_ONLY
 
+#include "svd.cpp"
 
+#endif
 
+
 
+#endif