barycentric.cpp

/*
Dependencies: LAPACK

Ubuntu: sudo apt-get install liblapack3 liblapack-dev

Compile: g++ -o barycentric barycentric.cpp -O3 -Wall -Wextra -pedantic -std=gnu++11 -llapack
*/

#include <iostream>
#include <cassert>
#include <algorithm>
#include <vector>

// do not use lapacke
//#include <lapacke.h>

// just describe here dgesv_ function
extern "C"
{
  void dgesv_(int * n, int * nrhs, double * a, int * lda,
    int * ipiv, double * b, int * ldb, int * info);
}

typedef std::vector<double> point_t;

point_t find_barycentric_coords(point_t const & point, std::vector<point_t> const & simplex)
{
  size_t size = simplex[0].size();

  // at least 2-dimensional
  assert(size >= 2);

  // simplex consists of size+1 vertecies
  assert(simplex.size() == size + 1);

  // given point and vertecies have right size
  assert(point.size() == size);
  for (auto const & vertex : simplex)
  {
    assert(vertex.size() == size);
  }

  std::vector<double> a(size * size);
  for (size_t i = 0; i < size; ++i)
  {
    for (size_t j = 0; j < size; ++j)
    {
      a[i*size+j] = simplex[i][j] - simplex[size][j];
    }
  }

  std::vector<double> b(size + 1);
  for (size_t i = 0; i < size; ++i)
  {
    b[i] = point[i] - simplex[size][i];
  }

  // additional variables needed for LAPACK
  std::vector<int> ipiv(size);
  int info;
  int size_int = static_cast<int>(size);
  int b_size = 1;

  // solve a*x = b
  // result will be in b
  dgesv_(&size_int, &b_size, &a[0], &size_int, &ipiv[0],
    &b[0], &size_int, &info);

  // calculate last coordinate
  b[size] = 1 - std::accumulate(b.begin(), b.end() - 1, 0.);

  return b;
}

int main()
{
  std::vector<point_t> simplex = {
    {1, 2},
    {17, 3},
    {9, 12}
  };
  point_t point = {10, 15};
  point_t results = find_barycentric_coords(point, simplex);
  for (auto r : results)
  {
    std::cout << r << std::endl;
  }
  return 0;
}