normgrid.cpp

/* Mondschein Engine is a fast, powerful, and easy-to-use 3D realtime rendering engine.
 *
 * Copyright (C) 2009-2013 by Andreas Amereller
 * E-Mail: andreas.amereller.dev@googlemail.com
 *
 * This program is free software; you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the Free
 * Software Foundation; either version 3 of the License, or (at your option)
 * any later version.
 *
 * This program is distributed in the hope that it will be useful, but WITHOUT
 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 * FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for
 * more details.
 *
 * You should have received a copy of the GNU General Public License along
 * with this program; if not, write to the Free Software Foundation, Inc.,
 * 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
 */

#include "core/normgrid.h"
#include "core/aitken_neville.h"

using namespace mondschein;
using namespace math;

std::array<Eigen::Vector3d,2> make_triangle(const Eigen::Vector3d &_a,const Eigen::Vector3d &_b,const Eigen::Vector3d &_c);
Eigen::Vector3d calc_normal(const std::array<Eigen::Vector3d,2> &_t);
Eigen::Vector3d interpolate_vector(const std::array<Eigen::Vector3d,4> &_v);
void expand(std::vector<std::vector<Eigen::Vector3d> > &_grid);

std::array<Eigen::Vector3d,2> make_triangle(const Eigen::Vector3d &_a,const Eigen::Vector3d &_b,const Eigen::Vector3d &_c)
{
  std::array<Eigen::Vector3d,2> t;
  t.at(0)=(_b-_a);
  t.at(1)=(_c-_a);
  return t;
}

Eigen::Vector3d calc_normal(const std::array<Eigen::Vector3d,2> &_t)
{
  return _t.at(0).cross(_t.at(1)).normalized();
}

Eigen::Vector3d interpolate_vector(const std::array<Eigen::Vector3d,4> &_v)
{
  Eigen::Vector3d v(0.0,0.0,0.0);
  for (uint32 i=0; i<3; ++i)
  {
    std::map<float64,float64> p;
    for (uint32 j=0; j<3; j++) p.insert(std::make_pair(j,_v.at(j)(i)));
    Aitken_Neville_p ip(new Aitken_Neville(p));
    v(i)=ip->interpolate(-1.0);
  }
  return v;
}

void expand(std::vector<std::vector<Eigen::Vector3d> > &_grid)
{
  std::vector<Eigen::Vector3d> b,e;
  std::array<Eigen::Vector3d,4> a;
  uint32 s1=_grid.size();
  uint32 s2=_grid.at(0).size();
  for (uint32 j=0; j<s2; ++j)
  {
    for (uint32 k=0; k<4; ++k) a.at(k)=_grid.at(k).at(j);
    b.push_back(interpolate_vector(a));
    for (uint32 k=0; k<4; ++k) a.at(k)=_grid.at(s1-1-k).at(j);
    e.push_back(interpolate_vector(a));
  }
  _grid.insert(_grid.begin(),b);
  _grid.push_back(e);
  s1=_grid.size();
  s2=_grid.at(0).size();
  for (uint32 i=0; i<s1; ++i)
  {
    for (uint32 k=0; k<4; ++k) a.at(k)=_grid.at(i).at(k);
    _grid.at(i).insert(_grid.at(i).begin(),interpolate_vector(a));
    for (uint32 k=0; k<4; ++k) a.at(k)=_grid.at(i).at(s2-2-k);
    _grid.at(i).push_back(interpolate_vector(a));
  }
  return;
}

Normal_Grid::Normal_Grid()
{
  return;
}

Normal_Grid::Normal_Grid(Normal_Grid_c _nc)
{
  return;
}

Normal_Grid::~Normal_Grid()
{
  return;
}

Normal_Grid &Normal_Grid::operator=(Normal_Grid_c _nc)
{
  return *this;
}

std::vector<Eigen::Vector3d> Normal_Grid::flat(const std::vector<std::vector<Eigen::Vector3d> > &_grid)
{
  std::vector<Eigen::Vector3d> n;
  for (uint32 i=0; i<_grid.size()-1; ++i)
  {
    for (uint32 j=0; j<_grid.at(i).size()-1; ++j)
    {
      n.push_back(calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i+1).at(j),_grid.at(i).at(j+1))));
      n.push_back(calc_normal(make_triangle(_grid.at(i+1).at(j+1),_grid.at(i).at(j+1),_grid.at(i+1).at(j))));
    }
  }
  return n;
}

std::vector<Eigen::Vector3d> Normal_Grid::flat_inv(const std::vector<std::vector<Eigen::Vector3d> > &_grid)
{
  std::vector<Eigen::Vector3d> n;
  for (uint32 i=0; i<_grid.size()-1; ++i)
  {
    for (uint32 j=0; j<_grid.at(i).size()-1; ++j)
    {
      n.push_back(calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i+1).at(j),_grid.at(i).at(j+1)))*(-1));
      n.push_back(calc_normal(make_triangle(_grid.at(i+1).at(j+1),_grid.at(i).at(j+1),_grid.at(i+1).at(j)))*(-1));
    }
  }
  return n;
}

std::vector<Eigen::Vector3d> Normal_Grid::intermediate(std::vector<std::vector<Eigen::Vector3d> > _grid)
{
  expand(_grid);
  std::vector<Eigen::Vector3d> n;
  Eigen::Vector3d tn;
  for (uint32 i=1; i<_grid.size()-1; ++i)
  {
    for (uint32 j=1; j<_grid.at(i).size()-1; ++j)
    {
      tn=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i+1).at(j),_grid.at(i+1).at(j+1)));
      tn+=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i+1).at(j+1),_grid.at(i).at(j+1)));
      tn+=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i).at(j+1),_grid.at(i-1).at(j+1)));
      tn+=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i-1).at(j+1),_grid.at(i-1).at(j)));
      tn+=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i-1).at(j),_grid.at(i-1).at(j-1)));
      tn+=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i-1).at(j-1),_grid.at(i).at(j-1)));
      tn+=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i).at(j-1),_grid.at(i).at(j+1)));
      tn+=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i).at(j+1),_grid.at(i+1).at(j))).normalized();
      n.push_back(tn);
    }
  }
  return n;
}

std::vector<Eigen::Vector3d> Normal_Grid::intermediate_inv(std::vector<std::vector<Eigen::Vector3d> > _grid)
{
  expand(_grid);
  std::vector<Eigen::Vector3d> n;
  Eigen::Vector3d tn;
  for (uint32 i=0; i<_grid.size()-1; ++i)
  {
    for (uint32 j=0; j<_grid.at(i).size()-1; ++j)
    {
      tn=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i+1).at(j),_grid.at(i+1).at(j+1)));
      tn+=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i+1).at(j+1),_grid.at(i).at(j+1)));
      tn+=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i).at(j+1),_grid.at(i-1).at(j+1)));
      tn+=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i-1).at(j+1),_grid.at(i-1).at(j)));
      tn+=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i-1).at(j),_grid.at(i-1).at(j-1)));
      tn+=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i-1).at(j-1),_grid.at(i).at(j-1)));
      tn+=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i).at(j-1),_grid.at(i).at(j+1)));
      tn+=calc_normal(make_triangle(_grid.at(i).at(j),_grid.at(i).at(j+1),_grid.at(i+1).at(j))).normalized();
      n.push_back(tn);
    }
  }
  return n;
}