#include "CImg.h"
#include <iostream>
#include <vector>
#include <algorithm>
#include "poissonimed.h"
#include "armadillo"
using namespace arma;
CImg < double >poissonmerge(const CImg < double >source,
const CImg < double >dest,
const vector < int >sourcex,
const vector < int >sourcey, const int destx,
const int desty)
{
int sourceminx = *min_element(sourcex.begin(), sourcex.end());
int sourcemaxx = *max_element(sourcex.begin(), sourcex.end());
int sourceminy = *min_element(sourcey.begin(), sourcey.end());
int sourcemaxy = *max_element(sourcey.begin(), sourcey.end());
unsigned int numrows = sourcemaxy - sourceminy + 1;
unsigned int numcols = sourcemaxx - sourceminx + 1;
int width = dest.width();
int height = dest.height();
int midx, midy;
int pix_idx, nbhd_idx, sourceimgx, sourceimgy, destimgx, destimgy;
bool pntstatus, nbhd_status;
unsigned char red[] = { 255, 0, 0 };
vector < int >destvertexx, destvertexy;
vector < int >regionvertexx, regionvertexy;
CImg < double >destcopy = dest;
CImgList < double >gradients;
CImg < double >vx, vy;
mat sysmat;
colvec rhs(numcols * numrows);
colvec newdestchunk;
ofstream outfile("poissonimed.log");
polyCentroid(sourcex.size(), sourcex, sourcey, &midx, &midy);
for (int idx = 0; idx < sourcex.size(); idx++) {
destvertexx.push_back(sourcex[idx] - midx + destx);
destvertexy.push_back(sourcey[idx] - midy + desty);
regionvertexx.push_back(sourcex[idx] - sourceminx);
regionvertexy.push_back(sourcey[idx] - sourceminy);
}
/* for(int idx=0; idx < sourcex.size()-1; idx++)
destcopy.draw_line(destvertexx[idx], destvertexy[idx], destvertexx[idx+1], destvertexy[idx+1], red);
destcopy.draw_line(destvertexx[0], destvertexy[0], destvertexx[sourcex.size()-1], destvertexy[sourcex.size()-1], red);
destcopy.draw_circle(destx, desty, 2, red);
*/
for (int channel = 0; channel < dest.spectrum(); channel++) {
sysmat = zeros < mat > (numcols * numrows, numcols * numrows);
rhs = zeros < colvec > (numcols * numrows);
gradients =
(CImg <
double >(source.get_channel(channel), false)).crop(sourceminx,
sourceminy,
sourcemaxx,
sourcemaxy).
get_gradient();
vx = gradients(0);
vy = gradients(1);
for (int xstep = 0; xstep < numcols; xstep++) {
for (int ystep = 0; ystep < numrows; ystep++) {
sourceimgx = sourceminx + xstep;
sourceimgy = sourceminy + ystep;
destimgx = sourceimgx - midx + destx;
destimgy = sourceimgy - midy + desty;
pix_idx = ystep * numcols + xstep;
pntstatus =
pntpoly(sourcex.size(), regionvertexx, regionvertexy,
xstep, ystep);
// if not inside of the region, then set to the destination image
if (!pntstatus) {
rhs(pix_idx) = dest(destimgx, destimgy, channel);
sysmat(pix_idx, pix_idx) = 1;
} else {
// if inside the region, check each neighboring pixel
// left neighbor
if (destimgx - 1 >= 0) {
nbhd_idx = ystep * numcols + xstep - 1;
sysmat(pix_idx, pix_idx)++;
if (pntpoly
(sourcex.size(), regionvertexx, regionvertexy,
xstep - 1, ystep)) {
sysmat(pix_idx, nbhd_idx) = -1;
rhs(pix_idx) +=
projInterpVect(vx(xstep, ystep),
vy(xstep, ystep),
vx(xstep - 1, ystep),
vy(xstep - 1, ystep), 1, 0);
} else
rhs(pix_idx) +=
dest(destimgx - 1, destimgy, channel);
}
// right neighbor
if (destimgx + 1 < width) {
nbhd_idx = ystep * numcols + xstep + 1;
sysmat(pix_idx, pix_idx)++;
if (pntpoly
(sourcex.size(), regionvertexx, regionvertexy,
xstep + 1, ystep)) {
sysmat(pix_idx, nbhd_idx) = -1;
rhs(pix_idx) +=
projInterpVect(vx(xstep, ystep),
vy(xstep, ystep),
vx(xstep + 1, ystep),
vy(xstep + 1, ystep), -1,
0);
} else
rhs(pix_idx) +=
dest(destimgx + 1, destimgy, channel);
}
// top neighbor
if (destimgy - 1 >= 0) {
nbhd_idx = (ystep - 1) * numcols + xstep;
sysmat(pix_idx, pix_idx)++;
if (pntpoly
(sourcex.size(), regionvertexx, regionvertexy,
xstep, ystep - 1)) {
sysmat(pix_idx, nbhd_idx) = -1;
rhs(pix_idx) +=
projInterpVect(vx(xstep, ystep),
vy(xstep, ystep), vx(xstep,
ystep -
1),
vy(xstep, ystep - 1), 0, 1);
} else
rhs(pix_idx) +=
dest(destimgx, destimgy - 1, channel);
}
// bottom neighbor
if (destimgx + 1 < height) {
nbhd_idx = (ystep + 1) * numcols + xstep;
sysmat(pix_idx, pix_idx)++;
if (pntpoly
(sourcex.size(), regionvertexx, regionvertexy,
xstep, ystep + 1)) {
sysmat(pix_idx, nbhd_idx) = -1;
rhs(pix_idx) +=
projInterpVect(vx(xstep, ystep),
vy(xstep, ystep), vx(xstep,
ystep +
1),
vy(xstep, ystep + 1), 0,
-1);
} else
rhs(pix_idx) +=
dest(destimgx, destimgy + 1, channel);
}
}
}
}
newdestchunk = solve(sysmat, rhs);
for (int xstep = 0; xstep < numcols; xstep++) {
for (int ystep = 0; ystep < numrows; ystep++) {
sourceimgx = sourceminx + xstep;
sourceimgy = sourceminy + ystep;
destimgx = sourceimgx - midx + destx;
destimgy = sourceimgy - midy + desty;
pix_idx = ystep * numcols + xstep;
destcopy(destimgx, destimgy, channel) =
newdestchunk(pix_idx);
}
}
}
return destcopy;
}
int main(int argc, char **argv)
{
CImg < double >sourceimg("source.jpg"), destimg("dest.jpg");
CImg < double >sourcedispimg = sourceimg, destdispimg = destimg;
CImgDisplay sourcedisp(sourcedispimg, "Select a polygonal region"),
destdisp(destdispimg, "Select an insertion point");
vector < int >sourcex, sourcey;
vector < int >::iterator xpointiter, ypointiter;
int destx, desty;
bool destset = false;
bool sourcepolyset = false;
const unsigned char red[] = { 255, 0, 0 };
while (!sourcedisp.is_closed() && !destdisp.is_closed()) {
CImgDisplay::wait(sourcedisp, destdisp);
if (sourcedisp.is_event())
switch (sourcedisp.button()) {
case CIMG_LEFT_MOUSEBUTTON:
sourcex.push_back(sourcedisp.mouse_x());
sourcey.push_back(sourcedisp.mouse_y());
sourcedispimg = sourceimg;
cout << "Current source polygon vertices: ";
xpointiter = sourcex.begin();
ypointiter = sourcey.begin();
while (xpointiter != sourcex.end()) {
if (xpointiter != sourcex.begin())
sourcedispimg.draw_line(*(xpointiter - 1),
*(ypointiter - 1),
*xpointiter, *ypointiter,
red);
cout << "(" << *(xpointiter++) << "," <<
*(ypointiter++) << ") ";
}
cout << endl;
cout <<
"(Use the right mouse button to close the polygon or middle mouse button to clear it)"
<< endl;
sourcedisp.display(sourcedispimg);
break;
case CIMG_RIGHT_MOUSEBUTTON:
// TODO: check the polygon is nice and has at least three vertices
cout << "Closed source polygon" << endl;
sourcepolyset = true;
sourcedispimg = sourceimg;
xpointiter = sourcex.begin();
ypointiter = sourcey.begin();
while (xpointiter != sourcex.end()) {
if (xpointiter != sourcex.begin())
sourcedispimg.draw_line(*(xpointiter - 1),
*(ypointiter - 1),
*xpointiter, *ypointiter,
red);
xpointiter++;
ypointiter++;
}
sourcedispimg.draw_line(sourcex[sourcex.size() - 1],
sourcey[sourcey.size() - 1],
sourcex[0], sourcey[0], red);
int midx, midy;
polyCentroid(sourcex.size(), sourcex, sourcey, &midx,
&midy);
sourcedispimg.draw_circle(midx, midy, 2, red);
sourcedisp.display(sourcedispimg);
sourcedisp.paint();
if (destset) {
cout << "Merging the images... " << endl;
destdispimg =
poissonmerge(sourceimg, destimg, sourcex, sourcey,
destx, desty);
destdisp.display(destdispimg).paint();
// TODO: reset the context so can