#include "DGtal/io/readers/VolReader.h"
#include "DGtal/io/viewers/Viewer3D.h"
#include "DGtal/topology/SurfelAdjacency.h"
#include "DGtal/topology/LightImplicitDigitalSurface.h"
#include "DGtal/topology/DigitalSurface.h"
#include "DGtal/math/linalg/EigenSupport.h"
#include "DGtal/dec/DiscreteExteriorCalculus.h"
#include "DGtal/dec/DiscreteExteriorCalculusSolver.h"
#include "DGtal/dec/DiscreteExteriorCalculusFactory.h"
#include "ConfigExamples.h"
template <typename Predicate, typename Domain>
struct FalseOutsideDomain
{
typedef typename Predicate::Point Point;
predicate(&predicate), domain(&domain)
{
}
bool
operator()(const Point& point) const
{
if (!domain->
isInside(point))
return false;
return (*predicate)(point);
}
const Predicate* predicate;
const Domain* domain;
};
void
alcapone_3d()
{
using std::endl;
const std::string filename = examplesPath + "samples/Al.100.vol";
typedef FalseOutsideDomain<Image, DGtal::Z3i::Domain> ImageExtended;
const ImageExtended image_extended(image, domain);
const SurfelAdjacency surfel_adjacency(true);
const DigitalSurfaceContainer digital_surface_container(kspace, image_extended, surfel_adjacency, cell_bel);
const DigitalSurface digital_surface(digital_surface_container);
trace.
info() <<
"digital_surface_size=" << digital_surface.size() << endl;
{
viewer->setWindowTitle("alcapone surface");
for (DigitalSurface::ConstIterator si=digital_surface.begin(), se=digital_surface.end(); si!=se; si++)
{
(*viewer) << cell;
}
}
const Calculus calculus = CalculusFactory::createFromNSCells<2>(digital_surface.begin(), digital_surface.end(), false);
trace.
info() <<
"calculus=" << calculus << endl;
{
viewer->setWindowTitle("alcapone calculus");
}
Calculus::DualForm0 rho(calculus);
for (int index=0; index<rho.length(); index++)
{
const Calculus::SCell cell = rho.getSCell(index);
const Calculus::Point point = kspace.
sKCoords(cell);
if (point[2]>1) continue;
rho.myContainer(index) = point[1]>100 ? 1 : -1;
}
trace.
info() <<
"rho_range=" << rho.myContainer.minCoeff() <<
"/" << rho.myContainer.maxCoeff() << endl;
{
viewer->setWindowTitle("alcapone poisson rho");
}
const Calculus::DualIdentity0 laplace = calculus.laplace<
DUAL>();
PoissonSolver solver;
const Calculus::DualForm0 phi = solver.
solve(rho);
trace.
info() <<
"phi_range=" << phi.myContainer.minCoeff() <<
"/" << phi.myContainer.maxCoeff() << endl;
{
viewer->setWindowTitle("alcapone poisson phi");
}
}
void
pyramid_3d()
{
using std::endl;
trace.
info() <<
"input_domain=" << input_domain << endl;
kspace.
init(input_domain.lowerBound(), input_domain.upperBound(),
true);
trace.
info() <<
"input_set_size=" << input_set.size() << endl;
{
viewer->setWindowTitle("input set");
(*viewer) << input_set;
}
const SurfelAdjacency surfel_adjacency(true);
const DigitalSurfaceContainer digital_surface_container(kspace, input_set, surfel_adjacency, cell_bel);
const DigitalSurface digital_surface(digital_surface_container);
for (int ii=0; ii<3; ii++)
{
std::stringstream ss;
for (int jj=0; jj<3; jj++)
ss << surfel_adjacency.getAdjacency(ii, jj) << " ";
}
trace.
info() <<
"digital_surface_container=" << digital_surface_container << endl;
trace.
info() <<
"digital_surface_size=" << digital_surface.size() << endl;
{
viewer->setWindowTitle("digital surface");
for (DigitalSurface::ConstIterator si=digital_surface.begin(), se=digital_surface.end(); si!=se; si++)
{
(*viewer) << cell;
}
}
const Calculus calculus = CalculusFactory::createFromNSCells<2>(digital_surface.begin(), digital_surface.end());
trace.
info() <<
"calculus=" << calculus << endl;
{
viewer->setWindowTitle("discrete exterior calculus");
}
const Calculus::PrimalForm2 primal_surfel_area = calculus.hodge<0,
DUAL>() *
Calculus::DualForm0(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(0, DUAL)));
const Calculus::DualForm2 dual_surfel_area = calculus.hodge<0,
PRIMAL>() *
Calculus::PrimalForm0(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(0, PRIMAL)));
const double area_th = calculus.kFormLength(2, PRIMAL);
const double area_primal = dual_surfel_area.myContainer.sum();
const double area_dual = primal_surfel_area.myContainer.sum();
trace.
info() <<
"area_th=" << area_th << endl;
trace.
info() <<
"area_primal=" << area_primal << endl;
trace.
info() <<
"area_dual=" << area_dual << endl;
FATAL_ERROR( area_th == area_primal );
FATAL_ERROR( area_th == area_dual );
const Calculus::DualForm1 dual_edge_length = calculus.hodge<1,
PRIMAL>() *
Calculus::PrimalForm1(calculus, Eigen::VectorXd::Ones(calculus.kFormLength(1, PRIMAL)));
const bool dual_edge_length_match = dual_edge_length.myContainer == Eigen::VectorXd::Ones(calculus.kFormLength(1, DUAL)) ;
trace.
info() <<
"dual_edge_length_match=" << dual_edge_length_match << endl;
FATAL_ERROR( dual_edge_length_match );
for (Calculus::Index index=0; index<dual_surfel_area.length(); index++)
{
const Calculus::SCell cell = dual_surfel_area.getSCell(index);
const Calculus::Scalar area = dual_surfel_area.myContainer(index);
trace.
info() << index <<
" " << cell <<
" " << area << endl;
}
}
int main(int argc, char* argv[])
{
QApplication app(argc, argv);
alcapone_3d();
return app.exec();
}