Merge branch 'mask'

Seems to work

CMakeLists.txt

@@ -9,6 +9,7 @@ project (seam-carving)
 
 # Add an executable with the above sources
 add_executable(seam-carving src/seam-carving.cpp)
+add_library(utils src/utils.cpp)
 
 # Set the directories that should be included in the build command for this target
 # when running g++ these will be included as -I/directory/path/
@@ -16,3 +17,9 @@ target_include_directories(seam-carving
     PRIVATE 
         ${PROJECT_SOURCE_DIR}/deps
 )
+target_link_libraries(seam-carving utils)
+
+target_include_directories(utils
+    PRIVATE 
+        ${PROJECT_SOURCE_DIR}/deps
+)

src/image.h

@@ -0,0 +1,19 @@
+#ifndef STBIDEF
+#ifdef STB_IMAGE_WRITE_STATIC
+#define STBIDEF  static
+#else
+#ifdef __cplusplus
+#define STBIDEF  extern "C"
+#else
+#define STBIDEF  extern
+#endif
+#endif
+#endif
+
+typedef unsigned char stbi_uc;
+
+STBIDEF int stbi_write_png(char const *filename, int w, int h, int comp, const void  *data, int stride_in_bytes);
+
+STBIDEF stbi_uc *stbi_load(char const *filename, int *x, int *y, int *comp, int req_comp);
+
+STBIDEF void     stbi_image_free      (void *retval_from_stbi_load);

src/seam-carving.cpp

@@ -1,4 +1,5 @@
 #include <iostream>
+#include <cassert>
 #include <random>
 #include <string>
 #include <vector>
@@ -6,26 +7,22 @@
 #include <CLI11.hpp>
 
 // Image filtering and I/O
-#define STB_IMAGE_IMPLEMENTATION
-#include <stb_image.h>
-#define STB_IMAGE_WRITE_IMPLEMENTATION
 #include <SimpleProgressBar.hpp>
-#include <stb_image_write.h>
+#include "utils.hpp"
+#include "image.h"
+
+#define DEFAULT_SEAMS 1
 
 // Global flags
 bool silent = false;
 bool test_energy = false;
-bool energy_recompute_all = false;
+bool until_mask_removal = false;
 int max_step = 1;
 
 // Get index for any table indexed by [width*(i : height) + (j : width)], but a
 // : dim_long, b : dim_large
 #define im_index(a, b) (vertical ? (width * a + b) : (width * b + a))
 
-bool nearly_equal(float a, float b) {
-  return std::nextafter(a, std::numeric_limits<float>::lowest()) <= 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) {
@@ -42,8 +39,8 @@ void export_image(const char *filename, const void *data, int width, int height,
 std::string function = "grad";
 
 #define compute_energy_for_pixel(source, width, height, i, j, nbChannels,      \
-                                 nbColorChannels, energy)                      \
   /* computes the energy at pixel i,j, i.e. energy[width*j+i]*/                \
+                                 nbColorChannels, dest)                        \
   auto indexPixel = (nbChannels) * (width * (j) + (i));                        \
   auto indexPixel_up =                                                         \
       ((j) - 1 > 0) ? (nbChannels) * (width * ((j) - 1) + (i)) : indexPixel;   \
@@ -55,10 +52,10 @@ std::string function = "grad";
   auto indexPixel_right = ((i) + 1 < width)                                    \
                               ? (nbChannels) * (width * (j) + ((i) + 1))       \
                               : indexPixel;                                    \
-  energy[width * j + i] = 0;                                                   \
+  dest =0;								\
   if (function == "gradnorm") {                                                \
     for (auto ch = 0; ch < (nbColorChannels); ch++) {                          \
-      energy[(width) * (j) + (i)] +=                                           \
+      dest+=								\
           (std::pow(fabs((float)source[indexPixel_up + ch] -                   \
                          source[indexPixel + ch]),                             \
                     2) +                                                       \
@@ -76,7 +73,7 @@ std::string function = "grad";
     }                                                                          \
   } else if (function == "grad") {                                             \
     for (auto ch = 0; ch < (nbColorChannels); ch++) {                          \
-      energy[(width) * (j) + (i)] +=                                           \
+      dest +=                                           \
           ((fabs((float)source[indexPixel_up + ch] -                           \
                  source[indexPixel + ch])) +                                   \
                                                                                \
@@ -90,7 +87,7 @@ std::string function = "grad";
     }                                                                          \
   } else if (function == "gradnorminf") {                                      \
     for (auto ch = 0; ch < (nbColorChannels); ch++) {                          \
-      energy[(width) * (j) + (i)] +=                                           \
+      dest+=                                           \
           std::max(std::max((fabs((float)source[indexPixel_up + ch] -          \
                                   source[indexPixel + ch])),                   \
                                                                                \
@@ -107,7 +104,6 @@ std::string function = "grad";
     std::cerr << "function " << function << " not available" << std::endl;     \
     exit(1);                                                                   \
   };
-
 // 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
 // slower...
@@ -122,20 +118,21 @@ std::vector<float> energy_e1(std::vector<unsigned char> source, int width,
 
   for (auto i = 0; i < width; i++) {
     for (auto j = 0; j < height; j++) {
-      compute_energy_for_pixel(source, width, height, i, j, nbChannels,
+      compute_energy_for_pixel(
-                               nbColorChannels, energy);
+        source, width, height,
-      energy[width * j + i] = (energy[width * j + i]);
+        i, j, nbChannels, nbColorChannels, energy[width*j+i]
+      );
     }
   }
   return energy;
 }
 
-std::vector<int> optimal_seam(std::vector<float> energy, int width, int height,
-                              bool vertical) {
 /** Given the energy value, returns the optimal seam */
+template <typename T> std::vector<int> optimal_seam(std::vector<T> energy, int width, int height,
+                              bool vertical) {
 
   // dyn_energy is indexed by [dim_large*(i : dim_long) + (j : dim_large)]
-  std::vector<float> dyn_energy(width * height);
+  std::vector<T> dyn_energy(width * height);
 
   int dim_large = vertical ? width : height;
   int dim_long = vertical ? height : width; // Number of elements in the seam
@@ -147,14 +144,15 @@ std::vector<int> optimal_seam(std::vector<float> energy, int width, int height,
 
   for (auto i = 1; i < dim_long; i++) { // Propagate dyn_energy
     for (auto j = 0; j < dim_large; j++) {
-      dyn_energy[dim_large * i + j] = __FLT_MAX__;
+      dyn_energy[dim_large * i + j] = limits::max_energy();
       int lower_bound = std::max(j - max_step, 0);
       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
-        dyn_energy[dim_large * i + j] = std::min(
+          if (dyn_energy[dim_large*(i-1)+k] < dyn_energy[dim_large*i+j]) {
-            dyn_energy[dim_large * i + j], 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)];
     }
@@ -163,7 +161,8 @@ std::vector<int> optimal_seam(std::vector<float> energy, int width, int height,
   std::vector<int> result(dim_long);
   // Find the seam end
   int min_idx = -1;
-  float min_val = __FLT_MAX__;
+  T min_val = limits::max_energy();
+
   for (auto j = 0; j < dim_large; j++) {
     if (min_val > dyn_energy[dim_large * (dim_long - 1) + j]) {
       min_idx = j;
@@ -184,7 +183,7 @@ std::vector<int> optimal_seam(std::vector<float> energy, int width, int height,
       (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
+    // This is not the nicest way to do this but we want to check in priority
     // at the center to have straight seams
     bool found = false;
     if (is_next_idx(min_idx)) {
@@ -215,8 +214,8 @@ std::vector<int> optimal_seam(std::vector<float> energy, int width, int height,
 
   return result;
 }
-/** Carves an image by one seam. Returns the optimal seam used */
 
+/** Carves an image by one seam. Returns the optimal seam used */
 template <typename T>
 void remove_seam(const std::vector<T> source, std::vector<T> &output, int width,
                  int height, int nbChannels, bool vertical,
@@ -239,75 +238,117 @@ void remove_seam(const std::vector<T> source, std::vector<T> &output, int width,
     }
   }
 }
-void update_energy_along_seam(std::vector<int> opt_seam,
+
-                              std::vector<unsigned char> output_img,
+
-                              std::vector<float> &output_energy, int height,
+// It would be preferable to use templates only for the value assignation but this is in fact far less efficient
-                              int width, int nbChannels, bool vertical,
+void recompute_energy_along_seam(
-                              int nbColorChannels) {
+    std::vector<unsigned char> carved_img, std::vector<float> &output_energy, std::vector<int> opt_seam,
-  // Recompute the energy along the seam
+    int width, int height,  int nbChannels, int nbColorChannels, bool vertical
-  if (energy_recompute_all) {
+  ) {
-    std::vector<float> energy =
+  int dim_large = vertical ? width : height;
-        energy_e1(output_img, vertical ? width - 1 : width,
+  int dim_long = vertical ? height : width;
-                  vertical ? height : height - 1, nbChannels);
+
-    std::copy(energy.begin(), energy.end(), output_energy.begin());
+  int newWidth  = vertical ? width-1 : width;
-  } else {
+  int newHeight = vertical ? height : height-1;
-    if (vertical) { // Vertical seam
+
-      for (auto j0 = 0; j0 < height; j0++) {
+  for (auto j = 0; j < dim_long; j++) {
-        auto i0 = opt_seam[j0];
+    for (auto i = -1; i < 2; i++) {
-        // recompute on a square around i0,j0
+      int x = vertical ? (opt_seam[j] + i) : j;
-        for (auto i_offset = -max_step - 2; i_offset <= max_step + 2;
+      int y = vertical ? j : (opt_seam[j] + i);
-             i_offset++) {
+      if ((0 < (opt_seam[j] + i)) && ((opt_seam[j] + i) < dim_large - 1)) {
-          for (auto j_offset = -max_step - 2; j_offset <= max_step + 2;
-               j_offset++) {
-            if (((0 < (i0 + i_offset)) && ((i0 + i_offset) < width - 1)) &&
-                ((0 < j0 + j_offset) && (j0 + j_offset < (height)))) {
         compute_energy_for_pixel(
-                  output_img, (width - 1), height, (i0 + i_offset),
+          carved_img, newWidth, newHeight,
-                  (j0 + j_offset), nbChannels, nbColorChannels, output_energy);
+          x, y, nbChannels,
+          nbColorChannels, output_energy[width*y+x]
+        );
       }
     }
   }
 }
-    } else {
+
-      for (auto i0 = 0; i0 < width; i0++) {
+void recompute_energy_along_seam(
-        auto j0 = opt_seam[i0];
+  std::vector<unsigned char> carved_img, std::vector<std::pair<int, float>> &output_energy, std::vector<int> opt_seam,
-        for (auto j_offset = -max_step - 2; j_offset <= max_step + 2;
+  int width, int height,  int nbChannels, int nbColorChannels, bool vertical
-             j_offset++) {
+  ) {
-          for (auto i_offset = -max_step - 2; i_offset <= max_step + 2;
+  int dim_large = vertical ? width : height;
-               i_offset++) {
+  int dim_long = vertical ? height : width;
-            if (((0 < (j0 + j_offset)) && ((j0 + j_offset) < height - 1)) &&
+
-                ((0 < (i0 + i_offset)) && (i0 + i_offset) < width)) {
+  int newWidth  = vertical ? width-1 : width;
+  int newHeight = vertical ? height : height-1;
+
+  for (auto j = 0; j < dim_long; j++) {
+    for (auto i = -1; i < 2; i++) {
+      int x = vertical ? (opt_seam[j] + i) : j;
+      int y = vertical ? j : (opt_seam[j] + i);
+      if ((0 < (opt_seam[j] + i)) && ((opt_seam[j] + i) < dim_large - 1)) {
         compute_energy_for_pixel(
-                  output_img, width, height - 1, (i0 + i_offset),
+          carved_img, newWidth, newHeight,
-                  (j0 + j_offset), nbChannels, nbColorChannels, output_energy);
+          x, y, nbChannels,
+          nbColorChannels, output_energy[width*y+x].first
+        );
       }
     }
   }
 }
-    }
+
-  }
+
-}
-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 */
+template <typename T>
+
+std::vector<int> carving_step(const std::vector<unsigned char> source_img,
+                              std::vector<T> source_energy,
+                              std::vector<unsigned char> &output_img,
+                              std::vector<T> &output_energy, int width,
+                              int height, int nbChannels, int nbColorChannels,
+                              bool vertical) {
+
   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);
-  update_energy_along_seam(opt_seam, output_img, output_energy, height, width,
+  // Recompute the energy along the seam, we need a separate function for templating
-                           nbChannels, vertical, nbColorChannels);
+  recompute_energy_along_seam(
+    source_img, output_energy, opt_seam,
+    width, height, nbChannels, nbColorChannels,
+    vertical
+  );
+
+
   return opt_seam;
 }
 
+
+
+
+std::vector<std::pair<int, float>> mask_energy(std::vector<float> energy,
+                                              int width, int height, unsigned char* mask) {
+  std::vector<std::pair<int, float>> output(width*height);
+
+  for (auto i=0; i < width*height; i++) {
+    output[i] = {mask[i], energy[i]};
+  }
+
+  return output;
+}
+
+bool does_seam_remove_mask(unsigned char* mask, int width, int height, int nbChannels,
+                          std::vector<int> opt_seam, bool vertical)
+{
+  int dim_large = vertical ? width : height;
+  int dim_long = vertical ? height : width;
+
+  for (int i=0; i < dim_long; i++) {
+    if (mask[im_index(i, opt_seam[i])] == 0) return true;
+  }
+  return false;
+}
+
 auto  seam_carving(unsigned char *source, int width, int height, int nbChannels,
                   int nbSeams, bool vertical, unsigned char* mask=nullptr) {
+
   int nbColorChannels = nbChannels > 3 ? 3 : nbChannels;
   int curWidth = width;
   int curHeight = height;
@@ -320,6 +361,9 @@ auto seam_carving(unsigned char *source, int width, int height, int nbChannels,
   // Contains at each step the carved image
   std::vector<float> source_energy(width * height);
   // Contains at each step the carved energy
+  std::vector<std::pair<int, float>> masked_energy;
+  std::vector<std::pair<int, float>> output_masked_energy(width*height);
+  // Source energy with (-1, 0, 1) on first element according to mask value
   std::vector<unsigned char> output_img(width * height * nbChannels);
   // Receives at each step the newly carved image
   std::vector<float> output_energy(width * height);
@@ -336,6 +380,8 @@ auto seam_carving(unsigned char *source, int width, int height, int nbChannels,
   }
 
   source_energy = energy_e1(source_img, width, height, nbChannels);
+  if (mask)
+    masked_energy = mask_energy(source_energy, width, height, mask);
 
   if (test_energy) {
     ini_energy = energy_e1(source_img, width, height, nbChannels);
@@ -364,18 +410,46 @@ auto seam_carving(unsigned char *source, int width, int height, int nbChannels,
       }
       if (nbChannels == 4)
         test_energy_output[nbChannels * k + 3] = source_img[nbChannels * k + 3];
+
+      if (mask) {
+        if (mask[k] == 2) // Green
+          test_energy_output[nbChannels*k + 1] = 125;
+        else if (mask[k] == 0) {// Red
+          test_energy_output[nbChannels * k] = 125;
+        }
+      }
     }
   }
 
-  SimpleProgressBar::ProgressBar bar(nbSeams);
+  SimpleProgressBar::ProgressBar bar(until_mask_removal ? dim_large : nbSeams);
+  bar.print();
 
-  for (auto seam_index = 0; seam_index < nbSeams; seam_index++) {
+  auto seam_index = 0;
-    std::vector<int> opt_seam = carving_step(
+  while (seam_index++ < nbSeams || until_mask_removal) {
-        source_img, source_energy, output_img, output_energy, curWidth,
+    std::vector<int> opt_seam;
-        curHeight, nbChannels, nbColorChannels, vertical);
+    if (mask) {
+      opt_seam = carving_step(
+        source_img, masked_energy, output_img, output_masked_energy,
+        curWidth, curHeight, nbChannels, nbColorChannels, vertical
+      );
+
+      if (
+        until_mask_removal &&
+        !does_seam_remove_mask(mask, width, height, nbChannels, opt_seam, vertical)
+      ) break;
+
+    } else {
+      opt_seam = carving_step(
+        source_img, source_energy, output_img, output_energy,
+        curWidth, 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());
+    if (mask)
+      std::copy(output_masked_energy.begin(), output_masked_energy.end(),
+        masked_energy.begin());
 
     vertical ? curWidth-- : curHeight--; // We just reduced the dimension
 
@@ -400,6 +474,7 @@ auto seam_carving(unsigned char *source, int width, int height, int nbChannels,
     bar.increment();
     bar.print();
   }
+  std::cout << std::endl; // Add newline after ProgressBar
 
   if (test_energy) {
 
@@ -415,7 +490,7 @@ int main(int argc, char **argv) {
   std::string maskImage;
   std::string sourceImage;
   std::string outputImage = "output.png";
-  int nbSeams = 1;
+  int nbSeams = DEFAULT_SEAMS;
   bool vertical = false;
 
   app.add_option("-s,--source", sourceImage, "Source image")
@@ -435,10 +510,10 @@ int main(int argc, char **argv) {
   app.add_flag("--silent", silent, "No verbose messages");
   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");
   app.add_option("-f,--function", function,
                  "The function to apply to compute the energy");
+  app.add_flag("-u,--until-mask-removal", until_mask_removal,
+               "Carve the image until there are no more red pixels in the mask");
   CLI11_PARSE(app, argc, argv);
 
   // Image loading
@@ -465,14 +540,27 @@ int main(int argc, char **argv) {
       unsigned char r, g, b;
       for (auto i=0; i < width*height; i++) {
         r = mask[maskChannels*i];
-      g = mask[maskChannels * i];
+        g = mask[maskChannels*i+1];
-      b = mask[maskChannels * i];
+        b = mask[maskChannels*i+2];
-      mask[2 * i] = (r == 0) && (g == 255) && (b == 0);
+        bool positive = (g > r && g > b && g > 100); // Mask images are not always the cleanest
-      mask[2 * i + 1] = (r == 255) && (g == 0) && (b == 0);
+        bool negative = (r > g && r > b && r > 100);
+
+        mask[i] = positive ? 2 : (negative ? 0 : 1);
       }
-    //* From now on, mask has the same dimensions as source and exactly 2
+      //* From now on, mask has the same dimensions as source and one single channel
-    // channels
+      //* The values are:
-    //* The first channel is positive, the second one negative.
+      //* . (2) we want to keep the pixel
+      //* . (1) nothing in particular
+      //* . (0) we want to remove the pixel
+  } 
+  if (until_mask_removal && maskImage.empty()) {
+    std::cerr << "Flag --until-mask-removal but no mask provided." << std::endl;
+    until_mask_removal = false;
+  }
+  if (until_mask_removal && nbSeams != DEFAULT_SEAMS) {
+    std::cerr << "Flag --nb-seams specified but --until-mask-removal provided." << std::endl;
+    nbSeams = DEFAULT_SEAMS;
+
   }
 
   nbSeams = std::min(

src/utils.cpp

@@ -0,0 +1,26 @@
+#define STB_IMAGE_IMPLEMENTATION
+#include <stb_image.h>
+#define STB_IMAGE_WRITE_IMPLEMENTATION
+#include <stb_image_write.h>
+
+#include <iostream>
+#include <random>
+#include <vector>
+
+bool nearly_equal(float a, float b) {
+  return std::nextafter(a, std::numeric_limits<float>::lowest()) <= b &&
+    std::nextafter(a, std::numeric_limits<float>::max()) >= b;
+}
+
+std::pair<int, float> operator+(std::pair<int, float>& p1, std::pair<int, float>& p2) {
+  return {
+    // If one of the two pixels is "protected" (ie 2), we want to prevent this line removing
+    // Else, we want to keep the information of "is there a pixel to remove in this seam" (ie 0)
+    (p1.first==2) || (p2.first==2) ? 2 : std::min(p1.first, p2.first),
+    p1.second+p2.second
+  };
+}
+
+void operator+=(std::pair<int, float>& p1, std::pair<int, float>& p2) {
+  p1 = p1+p2;
+}

src/utils.hpp

@@ -0,0 +1,19 @@
+#include <random>
+
+
+bool nearly_equal(float a, float b);
+
+std::pair<int, float> operator+(std::pair<int, float>& p1, std::pair<int, float>& p2);
+
+void operator+=(std::pair<int, float>& p1, std::pair<int, float>& p2);
+
+namespace limits {
+  struct max_energy {
+    template<class T> operator T() {
+      return std::numeric_limits<T>::max();
+    }
+    operator std::pair<int, float>() {
+      return {3, std::numeric_limits<float>::max()};
+    }
+  };
+}