18 #ifndef PIC_ALGORITHMS_POISSON_SOLVER_HPP 19 #define PIC_ALGORITHMS_POISSON_SOLVER_HPP 21 #include "../base.hpp" 23 #include "../image.hpp" 25 #ifndef PIC_DISABLE_EIGEN 27 #ifndef PIC_EIGEN_NOT_BUNDLED 28 #include "../externals/Eigen/Sparse" 29 #include "../externals/Eigen/src/SparseCore/SparseMatrix.h" 31 #include <Eigen/Sparse> 32 #include <Eigen/src/SparseCore/SparseMatrix.h> 39 #ifndef PIC_DISABLE_EIGEN 60 int tot = height * width;
63 printf(
"Init matrix...");
66 std::vector< Eigen::Triplet< double > > tL;
68 for(
int i = 0; i < height; i++) {
71 for(
int j = 0; j < width; j++) {
74 tL.push_back(Eigen::Triplet< double > (indI, indI, 4.0f));
76 if(((indI + 1) < tot) &&
77 ((indI % width) != (width - 1))) {
79 tL.push_back(Eigen::Triplet< double > (indI, indI + 1, -1.0f));
80 tL.push_back(Eigen::Triplet< double > (indI + 1, indI, -1.0f));
85 for(
int i = 0; i < (tot - width); i++) {
86 tL.push_back(Eigen::Triplet< double > (i + width, i , -1.0f));
87 tL.push_back(Eigen::Triplet< double > (i , i + width, -1.0f));
95 Eigen::SparseMatrix<double> A = Eigen::SparseMatrix<double>(tot, tot);
96 A.setFromTriplets(tL.begin(), tL.end());
97 Eigen::SimplicialCholesky<Eigen::SparseMatrix<double> > solver(A);
99 for(
int k = 0; k < f->
channels; k++) {
101 Eigen::VectorXd b, x;
102 b = Eigen::VectorXd::Zero(tot);
105 for(
int i = 0; i < height; i++) {
106 int tmpI = i * width;
107 for(
int j = 0; j < width; j++) {
108 int indI = (tmpI + j);
115 if(solver.info() != Eigen::Success) {
117 printf(
"SOLVER FAILED!\n");
124 printf(
"SOLVER SUCCESS!\n");
127 for(
int i = 0; i < height; i++) {
128 int tmpI = i * width;
130 for(
int j = 0; j < width; j++) {
131 (*ret)(j, i)[k] =
float(x(tmpI + j));
150 std::vector<int> coords,
154 printf(
"Iterative Poisson solver... ");
162 Image *tmpSwap = NULL;
166 for(
int i = 0; i < maxSteps; i++) {
167 for(
unsigned int j = 0; j < coords.size(); j++) {
168 int coord = coords[j];
171 float workValue = -laplacian->
data[coord];
174 workValue += img->
data[c];
177 workValue += img->
data[c];
180 workValue += img->
data[c];
183 workValue += img->
data[c];
185 tmpImg->
data[coord] = workValue / 4.0f;
void reverseAddress(int ind, int &x, int &y)
reverseAddress computes (x, y) given a memory address
Definition: image.hpp:680
float * data
data is the main buffer where pixel values are stored.
Definition: image.hpp:91
int channels
Definition: image.hpp:80
int getAddress(int x, int y)
getAddress calculates a memory address from (x, y)
Definition: image.hpp:650
PIC_INLINE Image * computePoissonSolverIterative(Image *img, Image *laplacian, std::vector< int > coords, int maxSteps=100)
computePoissonSolverIterative
Definition: poisson_solver.hpp:149
PIC_INLINE Image * computePoissonSolver(Image *f, Image *ret=NULL)
computePoissonSolver
Definition: poisson_solver.hpp:47
#define PIC_INLINE
Definition: base.hpp:33
The Image class stores an image as buffer of float.
Definition: image.hpp:60
Image * clone() const
Clone creates a deep copy of the calling instance.
Image * allocateSimilarOne()
allocateSimilarOne creates an Image with similar size of the calling instance.
Definition: bilateral_separation.hpp:25
int width
Definition: image.hpp:80
int height
Definition: image.hpp:80