ATTEMPT TO Merge branch 'do-not-recompute-energy'

This commit is contained in:
François Colin de Verdière 2025-04-01 17:22:43 +02:00
commit 1918b7bea4

View File

@ -1,103 +1,112 @@
#include <iostream> #include <iostream>
#include <string>
#include <random> #include <random>
#include <string>
#include <vector> #include <vector>
//Command-line parsing // Command-line parsing
#include <CLI11.hpp> #include <CLI11.hpp>
//Image filtering and I/O // Image filtering and I/O
#define STB_IMAGE_IMPLEMENTATION #define STB_IMAGE_IMPLEMENTATION
#include <stb_image.h> #include <stb_image.h>
#define STB_IMAGE_WRITE_IMPLEMENTATION #define STB_IMAGE_WRITE_IMPLEMENTATION
#include <stb_image_write.h>
#include <SimpleProgressBar.hpp> #include <SimpleProgressBar.hpp>
#include <stb_image_write.h>
//Global flag to silent verbose messages // Global flag to silent verbose messages
bool silent; bool silent;
bool test_energy; bool test_energy;
bool energy_recompute_all = false;
// Get index for any table indexed by [width*(i : height) + (j : width)], but a
// Get index for any table indexed by [width*(i : height) + (j : width)], but a : dim_long, b : dim_large // : dim_long, b : dim_large
#define im_index(a, b) \ #define im_index(a, b) (vertical ? (width * a + b) : (width * b + a))
(vertical ? (width*a + b) : (width*b + a))
bool nearly_equal(float a, float b) { bool nearly_equal(float a, float b) {
return std::nextafter(a, std::numeric_limits<float>::lowest()) <= b return std::nextafter(a, std::numeric_limits<float>::lowest()) <= b &&
&& std::nextafter(a, std::numeric_limits<float>::max()) >= b; std::nextafter(a, std::numeric_limits<float>::max()) >= b;
} }
void export_image(const char* filename, const void* data, int width, int height, int nbChannels) { void export_image(const char *filename, const void *data, int width, int height,
if (!silent) std::cout << "Exporting to \"" << filename << "\".." << std::endl; int nbChannels) {
int errcode = stbi_write_png(filename, width, height, nbChannels, data, nbChannels*width); if (!silent)
std::cout << "Exporting to \"" << filename << "\".." << std::endl;
int errcode = stbi_write_png(filename, width, height, nbChannels, data,
nbChannels * width);
if (!errcode) { if (!errcode) {
std::cerr << "Error while exporting the resulting image." << std::endl; std::cerr << "Error while exporting the resulting image." << std::endl;
exit(errcode); exit(errcode);
} }
} }
#define compute_energy_for_pixel(source, width, height, i, j, nbChannels, \
nbColorChannels, energy) \
auto indexPixel = (nbChannels) * (width * (j) + (i)); \
auto indexPixel_up = \
((j) - 1 > 0) ? (nbChannels) * (width * ((j) - 1) + (i)) : indexPixel; \
auto indexPixel_down = ((j) + 1 < height) \
? (nbChannels) * (width * ((j) + 1) + (i)) \
: indexPixel; \
auto indexPixel_left = \
((i) - 1 > 0) ? (nbChannels) * (width * (j) + ((i) - 1)) : indexPixel; \
auto indexPixel_right = ((i) + 1 < width) \
? (nbChannels) * (width * (j) + ((i) + 1)) \
: indexPixel; \
energy[width * j + i] = 0; \
for (auto ch = 0; ch < (nbColorChannels); ch++) { \
energy[(width) * (j) + (i)] += \
(fabs((float)source[indexPixel_up + ch] - source[indexPixel + ch]) + \
fabs((float)source[indexPixel_down + ch] - source[indexPixel + ch]) + \
fabs((float)source[indexPixel_left + ch] - source[indexPixel + ch]) + \
fabs((float)source[indexPixel_right + ch] - \
source[indexPixel + ch])); \
} \
// Le alpha n'est pas pris en compte dans l'énergie
// Here, we use this /ugly/ macro to avoid defining a function that would be way
// longer...
/** e_1 energy, energy is always normalized between 0 and 1 */ /** e_1 energy, energy is always normalized between 0 and 1 */
std::vector<float> energy_e1(std::vector<unsigned char> source, int width, int height, int nbChannels) { std::vector<float> energy_e1(std::vector<unsigned char> source, int width,
int nbColorChannels = nbChannels > 3 ? 3 : nbChannels; // nombre de canaux, excepté le alpha int height, int nbChannels) {
std::vector<float> energy(width*height); int nbColorChannels =
nbChannels > 3 ? 3 : nbChannels; // nombre de canaux, excepté le alpha
std::vector<float> energy(width * height);
float max_energy = 0; float max_energy = 0;
for(auto i=0 ; i < width; ++i) { for (auto i = 0; i < width; i++) {
for(auto j=0; j < height; ++j) { for (auto j = 0; j < height; j++) {
auto indexPixel = nbChannels*(width*j+i); compute_energy_for_pixel(source, width, height, i, j, nbChannels,
auto indexPixel_up = (j-1 > 0) ? nbChannels*(width*(j-1)+i) : indexPixel; nbColorChannels, energy);
auto indexPixel_down = (j+1 < height) ? nbChannels*(width*(j+1)+i) : indexPixel;
auto indexPixel_left = (i-1 > 0) ? nbChannels*(width*j+(i-1)) : indexPixel;
auto indexPixel_right = (i+1 < width) ? nbChannels*(width*j+(i+1)) : indexPixel;
energy[width*j+i] = 0;
for (auto ch=0; ch < nbColorChannels; ch++) { // Le alpha n'est pas pris en compte dans l'énergie
energy[width*j+i] += (
fabs((float)source[indexPixel_up+ch] - source[indexPixel+ch])
+fabs((float)source[indexPixel_down+ch] - source[indexPixel+ch])
+fabs((float)source[indexPixel_left+ch] - source[indexPixel+ch])
+fabs((float)source[indexPixel_right+ch] - source[indexPixel+ch])
);
}
max_energy = std::max(max_energy, energy[width*j+i]);
} }
} }
if (max_energy == 0) { return energy; }
for (auto k=0; k < width*height; k++) {
energy[k] = energy[k]/max_energy;
}
return energy; return energy;
} }
/** Given the energy value, returns the optimal seam */ std::vector<int> optimal_seam(std::vector<float> energy, int width, int height,
std::vector<int> optimal_seam(std::vector<float> energy, int width, int height, bool vertical, int max_step=1) { bool vertical, int max_step = 1) {
/** Given the energy value, returns the optimal seam */
// dyn_energy is indexed by [dim_large*(i : dim_long) + (j : dim_large)] // dyn_energy is indexed by [dim_large*(i : dim_long) + (j : dim_large)]
std::vector<float> dyn_energy(width*height); std::vector<float> dyn_energy(width * height);
int dim_large = vertical ? width : height; int dim_large = vertical ? width : height;
int dim_long = vertical ? height : width; // Number of elements in the seam int dim_long = vertical ? height : width; // Number of elements in the seam
//* Find an end of the minimal connected vertical/horizontal seam //* Find an end of the minimal connected vertical/horizontal seam
for (auto i=0; i < dim_large; i++) { for (auto i = 0; i < dim_large; i++) {
dyn_energy[i] = energy[i]; dyn_energy[i] = energy[i];
} }
for (auto i=1; i < dim_long; i++) { // Propagate dyn_energy for (auto i = 1; i < dim_long; i++) { // Propagate dyn_energy
for (auto j=0; j < dim_large; j++) { for (auto j = 0; j < dim_large; j++) {
dyn_energy[dim_large*i+j] = __FLT_MAX__; dyn_energy[dim_large * i + j] = __FLT_MAX__;
int lower_bound = std::max(j - max_step, 0); int lower_bound = std::max(j - max_step, 0);
int upper_bound = std::min(j+max_step, dim_large-1); int upper_bound = std::min(j + max_step, dim_large - 1);
for (auto k=lower_bound; k <= upper_bound; k++) { // Compute energy based on predecessors for (auto k = lower_bound; k <= upper_bound;
dyn_energy[dim_large*i+j] = std::min( k++) { // Compute energy based on predecessors
dyn_energy[dim_large*i+j], dyn_energy[dim_large * i + j] = std::min(
dyn_energy[dim_large*(i-1)+k] dyn_energy[dim_large * i + j], dyn_energy[dim_large * (i - 1) + k]);
);
} }
dyn_energy[dim_large*i+j] += energy[im_index(i, j)]; dyn_energy[dim_large * i + j] += energy[im_index(i, j)];
} }
} }
@ -105,179 +114,270 @@ std::vector<int> optimal_seam(std::vector<float> energy, int width, int height,
// Find the seam end // Find the seam end
int min_idx = -1; int min_idx = -1;
float min_val = __FLT_MAX__; float min_val = __FLT_MAX__;
for (auto j=0; j < dim_large; j++) { for (auto j = 0; j < dim_large; j++) {
if (min_val > dyn_energy[dim_large*(dim_long-1) + j]) { if (min_val > dyn_energy[dim_large * (dim_long - 1) + j]) {
min_idx = j; min_idx = j;
min_val = dyn_energy[dim_large*(dim_long-1) + j]; min_val = dyn_energy[dim_large * (dim_long - 1) + j];
} }
} }
result[dim_long-1] = min_idx; result[dim_long - 1] = min_idx;
//* Backtracking to find the path //* Backtracking to find the path
for (auto i=dim_long-1; i > 0; i--) { for (auto i = dim_long - 1; i > 0; i--) {
// We want to find either (bot_l, bot_c, bot_r) with dyn_energy[.] = min_val - energy[cur] // We want to find either (bot_l, bot_c, bot_r) with dyn_energy[.] = min_val -
// energy[cur]
//Idea : float next_energy = min_val - energy[width*i + min_idx]; // Idea : float next_energy = min_val - energy[width*i + min_idx];
//! With floats, we don't always have x + y - y == x, so we check is x+y == x+y //! With floats, we don't always have x + y - y == x, so we check is x+y == x+y
// This define is a bit ugly but 200x faster than using a lambda function // This define is a bit ugly but 200x faster than using a lambda function
#define is_next_idx(idx) \ #define is_next_idx(idx) \
(dyn_energy[(i-1)*dim_large + idx]+energy[im_index(i, min_idx)] == min_val) (dyn_energy[(i - 1) * dim_large + idx] + energy[im_index(i, min_idx)] == \
min_val)
// This is not the nicest way to do thiss but we want to check in priority at the center to have straight seams // This is not the nicest way to do thiss but we want to check in priority
// at the center to have straight seams
bool found = false; bool found = false;
if (is_next_idx(min_idx)) { if (is_next_idx(min_idx)) {
min_val = dyn_energy[(i-1)*dim_large+min_idx];
min_val = dyn_energy[(i - 1) * dim_large + min_idx];
found = true; found = true;
} }
for (auto k=1; !found && k <= max_step; k++) { for (auto k = 1; !found && k <= max_step; k++) {
if (min_idx+k < dim_large && is_next_idx(min_idx+k)) { if (min_idx + k < dim_large && is_next_idx(min_idx + k)) {
min_val = dyn_energy[(i-1)*dim_large+min_idx+k]; min_val = dyn_energy[(i - 1) * dim_large + min_idx + k];
min_idx = min_idx+k; min_idx = min_idx + k;
found = true; found = true;
} else if (min_idx-k >= 0 && is_next_idx(min_idx-k)) { } else if (min_idx - k >= 0 && is_next_idx(min_idx - k)) {
min_val = dyn_energy[(i-1)*dim_large+min_idx-k]; min_val = dyn_energy[(i - 1) * dim_large + min_idx - k];
min_idx = min_idx-k; min_idx = min_idx - k;
found = true; found = true;
} }
} }
if (!found) { if (!found) {
std::cerr << "Unable to backtrack path !" << std::endl; std::cerr << "Unable to backtrack path !" << std::endl;
exit(1); exit(1);
} }
result[i-1] = min_idx; result[i - 1] = min_idx;
} }
return result; return result;
} }
/** Carves an image by one seam. Returns the optimal seam used */ /** Carves an image by one seam. Returns the optimal seam used */
std::vector<int> carving_step(const std::vector<unsigned char> source, std::vector<unsigned char> &output, int width, int height, int nbChannels, bool vertical, int max_step=1) {
std::vector<float> energy = energy_e1(source, width, height, nbChannels);
std::vector<int> opt_seam = optimal_seam(energy, width, height, vertical, max_step=max_step);
std::vector<bool> blacklist(width*height);
template <typename T>
void remove_seam(const std::vector<T> source, std::vector<T> &output, int width,
int height, int nbChannels, bool vertical,
const std::vector<int> seam) {
// remove the given seam from the image, the result is in output
// the output must have at least the right size!
int dim_large = vertical ? width : height; int dim_large = vertical ? width : height;
int dim_long = vertical ? height : width; int dim_long = vertical ? height : width;
for (auto i = 0; i < dim_long; i++) {
for (auto k=0; k < width*height; k++) { blacklist[k] = false; }
for (auto i=0; i < dim_long; i++) {
int index = vertical ? opt_seam[i]+i*width : i+width*opt_seam[i];
blacklist[index] = true;
}
for (auto i=0; i < dim_long; i++) {
int cur_j = 0; int cur_j = 0;
for (auto j=0; cur_j < dim_large-1 && j < dim_large; j++) { for (auto j = 0; cur_j < dim_large - 1 && j < dim_large; j++) {
if (!blacklist[im_index(i, j)]) { if (seam[i] != j) {
int out_pixelIndex = nbChannels*(vertical ? ((width-1)*i + cur_j) : (width*cur_j + i)); int out_pixelIndex = nbChannels * (vertical ? ((width - 1) * i + cur_j)
int src_pixelIndex = nbChannels*im_index(i, j); : (width * cur_j + i));
int src_pixelIndex = nbChannels * im_index(i, j);
for (auto ch=0; ch < nbChannels; ch++) for (auto ch = 0; ch < nbChannels; ch++)
output[out_pixelIndex+ch] = source[src_pixelIndex+ch]; output[out_pixelIndex + ch] = source[src_pixelIndex + ch];
cur_j++; cur_j++;
} }
} }
} }
}
std::vector<int> carving_step(const std::vector<unsigned char> source_img,
std::vector<float> source_energy,
std::vector<unsigned char> &output_img,
std::vector<float> &output_energy, int width,
int height, int nbChannels, int nbColorChannels,
bool vertical) {
/** Carves an image and its energy by one seam, and recomputes the energy.
Returns the optimal seam used */
std::vector<int> opt_seam =
optimal_seam(source_energy, width, height, vertical);
remove_seam(source_img, output_img, width, height, nbChannels, vertical,
opt_seam);
remove_seam(source_energy, output_energy, width, height, 1, vertical,
opt_seam);
// Recompute the energy along the seam
if (energy_recompute_all) {
std::vector<float> energy =
energy_e1(output_img, vertical ? width - 1 : width,
vertical ? height : height - 1, nbChannels);
std::copy(energy.begin(), energy.end(), output_energy.begin());
} else {
// ASSUME WE ARE DOING A VERTICAL SEAM
if (vertical) {
for (auto j = 0; j < height; j++) {
for (auto i = -1; i < 2; i++) {
if ((0 < (opt_seam[j] + i)) && ((opt_seam[j] + i) < width - 1)) {
compute_energy_for_pixel(output_img, (width - 1), height,
(opt_seam[j] + i), j, nbChannels,
nbColorChannels, output_energy);
}
}
}
} else {
for (auto i = 0; i < width; i++) {
for (auto j = -1; j < 2; j++) {
if ((0 < (opt_seam[i] + j)) && ((opt_seam[i] + j) < height - 1)) {
compute_energy_for_pixel(output_img, width, height - 1, i,
(opt_seam[i] + j), nbChannels,
nbColorChannels, output_energy);
}
}
}
}
}
return opt_seam; return opt_seam;
} }
void seam_carving(unsigned char *source, int width, int height, int nbChannels,
void seam_carving(unsigned char* source, int width, int height, int nbChannels, const char* out_filename, int nbSeams, bool vertical, bool test_energy=false, int max_step=1) { const char *out_filename, int nbSeams, bool vertical,
bool test_energy = false, int max_step = 1) {
int nbColorChannels = nbChannels > 3 ? 3 : nbChannels; int nbColorChannels = nbChannels > 3 ? 3 : nbChannels;
int curWidth = width; int curWidth = width;
int curHeight = height; int curHeight = height;
int dim_large = vertical ? width : height; int dim_large = vertical ? width : height;
int dim_long = vertical ? height : width; int dim_long = vertical ? height : width;
// dim_long=longueur des seam
std::vector<unsigned char> carve_output(width*height*nbChannels); // Receives at each step the newly carved image std::vector<unsigned char> source_img(width * height * nbChannels);
std::vector<unsigned char> source_img(width*height*nbChannels); // Contains at each step the carved image // Contains at each step the carved image
std::vector<bool> complete_blacklist(width*height); // Contains all removed pixels, for "test_energy" std::vector<float> source_energy(width * height);
std::vector<float> ini_energy; // Contains the initial energy, only for "test_energy" // Contains at each step the carved energy
std::vector<unsigned char> test_energy_output(width*height*nbChannels); // Final output for "test_energy" std::vector<unsigned char> output_img(width * height * nbChannels);
// Receives at each step the newly carved image
std::vector<float> output_energy(width * height);
// Contains at each step the carved energy
std::vector<bool> complete_blacklist(width * height);
// Contains all removed pixels, for "test_energy"
std::vector<float> ini_energy;
// Contains the initial energy, only for "test_energy"
std::vector<unsigned char> test_energy_output(width * height * nbChannels);
// Final output for "test_energy"
for (auto i = 0; i < width * height * nbChannels; i++) {
source_img[i] = source[i];
}
for (auto i=0; i < width*height*nbChannels; i++) { source_img[i] = source[i]; } source_energy = energy_e1(source_img, width, height, nbChannels);
if (test_energy) { if (test_energy) {
ini_energy = energy_e1(source_img, width, height, nbChannels); ini_energy = energy_e1(source_img, width, height, nbChannels);
for (auto k=0; k < width*height; k++) { complete_blacklist[k] = false; } for (auto k = 0; k < width * height; k++) {
complete_blacklist[k] = false;
}
//* Prepare final output //* Prepare final output
for (auto k=0; k < width*height; k++) {
//for (auto i=0; i < nbColorChannels; i++) //* Uncomment if you prefer to see darkened source image
// output[nbChannels*k+i] = source_img[nbChannels*k+i]/nbChannels;
for (auto i=0; i < nbColorChannels; i++)
test_energy_output[nbChannels*k+i] = ini_energy[k]*255;
if (nbChannels==4) float max_energy = __FLT_MAX__;
test_energy_output[nbChannels*k+3] = source_img[nbChannels*k+3]; for (auto k = 0; k < width * height; k++) {
max_energy = fmax(max_energy, ini_energy[k]);
}
if (max_energy != 0) {
for (auto k = 0; k < width * height; k++) {
ini_energy[k] /= max_energy;
}
}
for (auto k = 0; k < width * height; k++) {
// for (auto i=0; i < nbColorChannels; i++) //* Uncomment if you prefer to
// see darkened source image
// output[nbChannels*k+i] = source_img[nbChannels*k+i]/nbChannels;
for (auto i = 0; i < nbColorChannels; i++)
test_energy_output[nbChannels * k + i] = ini_energy[k] * 255;
if (nbChannels == 4)
test_energy_output[nbChannels * k + 3] = source_img[nbChannels * k + 3];
} }
} }
SimpleProgressBar::ProgressBar bar(nbSeams); SimpleProgressBar::ProgressBar bar(nbSeams);
for (auto seam=0; seam < nbSeams; seam++) {
std::vector<int> opt_seam = carving_step(source_img, carve_output, curWidth, curHeight, nbChannels, vertical, max_step=max_step);
std::copy(carve_output.begin(), carve_output.end(), source_img.begin());
if (vertical) // We just reduced the dimension for (auto seam_index = 0; seam_index < nbSeams; seam_index++) {
curWidth--; std::vector<int> opt_seam = carving_step(
else source_img, source_energy, output_img, output_energy, curWidth,
curHeight--; curHeight, nbChannels, nbColorChannels, vertical);
std::copy(output_img.begin(), output_img.end(), source_img.begin());
std::copy(output_energy.begin(), output_energy.end(),
source_energy.begin());
vertical ? curWidth-- : curHeight--; // We just reduced the dimension
if (test_energy) { // Update blacklist if (test_energy) { // Update blacklist
for (auto i=0; i < dim_long; i++) { for (auto i = 0; i < dim_long; i++) {
int j, cur_j = 0; // cur_j is the index relative to the current carved image. j is absolute in the source image int j, cur_j = 0; // cur_j is the index relative to the current carved
for (j=0; j < dim_large && (cur_j < opt_seam[i] || complete_blacklist[im_index(i, j)]); j++) { // image. j is absolute in the source image
if (!complete_blacklist[im_index(i, j)]) { cur_j++; } for (j = 0; j < dim_large &&
(cur_j < opt_seam[i] || complete_blacklist[im_index(i, j)]);
j++) {
if (!complete_blacklist[im_index(i, j)]) {
cur_j++;
} }
assert(cur_j == opt_seam[i]); // Else, j == width and cur_j is not in the source image.. }
assert(cur_j == opt_seam[i]); // Else, j == width and cur_j is not in
// the source image..
complete_blacklist[im_index(i, j)] = true; complete_blacklist[im_index(i, j)] = true;
test_energy_output[nbChannels*im_index(i, j)] = 255; // Set carved pixel to red test_energy_output[nbChannels * im_index(i, j)] =
255; // Set carved pixel to red
} }
} }
bar.increment(); bar.increment();
bar.print(); bar.print();
} }
std::cout << std::endl;
if (test_energy) { if (test_energy) {
export_image(out_filename, test_energy_output.data(), width, height, nbChannels); export_image(out_filename, test_energy_output.data(), width, height,
nbChannels);
} else { } else {
export_image(out_filename, source_img.data(), curWidth, curHeight, nbChannels); export_image(out_filename, source_img.data(), curWidth, curHeight,
nbChannels);
} }
} }
int main(int argc, char **argv) { int main(int argc, char **argv) {
CLI::App app{"seam-carving"}; CLI::App app{"seam-carving"};
std::string sourceImage; std::string sourceImage;
app.add_option("-s,--source", sourceImage, "Source image")->required()->check(CLI::ExistingFile);; app.add_option("-s,--source", sourceImage, "Source image")
std::string outputImage= "output.png"; ->required()
->check(CLI::ExistingFile);
;
std::string outputImage = "output.png";
app.add_option("-o,--output", outputImage, "Output image")->required(); app.add_option("-o,--output", outputImage, "Output image")->required();
int nbSeams = 1; int nbSeams = 1;
app.add_option("-n,--nb-seams", nbSeams, "Number of seams"); app.add_option("-n,--nb-seams", nbSeams, "Number of seams");
int max_step = 1; int max_step = 1;
app.add_option("--max-step", max_step, "Max width of step to find a seam"); app.add_option("--max-step", max_step, "Max width of step to find a seam");
bool vertical = false; bool vertical = false;
app.add_flag("--vertical", vertical, "Vertical carving"); app.add_flag("--vertical", vertical,
"Vertical carving (remove vertical seams)");
silent = false; silent = false;
app.add_flag("--silent", silent, "No verbose messages"); app.add_flag("--silent", silent, "No verbose messages");
test_energy = false; test_energy = false;
app.add_flag("--test-energy", test_energy, "Don't resize image, just try the specified energy function"); app.add_flag("--test-energy", test_energy,
"Don't resize image, just try the specified energy function");
app.add_flag("--energy-recompute-all", energy_recompute_all,
"recompute the whole energy at each step");
CLI11_PARSE(app, argc, argv); CLI11_PARSE(app, argc, argv);
//Image loading // Image loading
int width, height, nbChannels; int width, height, nbChannels;
unsigned char *source = stbi_load(sourceImage.c_str(), &width, &height, &nbChannels, 0); unsigned char *source =
stbi_load(sourceImage.c_str(), &width, &height, &nbChannels, 0);
nbSeams = std::min(nbSeams, width); nbSeams = std::min(nbSeams, width);
//std::cout << "channels: " << nbChannels << std::endl; seam_carving(source, width, height, nbChannels, outputImage.c_str(), nbSeams,
vertical, test_energy = test_energy, max_step = max_step);
seam_carving(source, width, height, nbChannels, outputImage.c_str(), nbSeams, vertical, test_energy=test_energy, max_step=max_step);
stbi_image_free(source); stbi_image_free(source);
exit(0); exit(0);