Add MCLaplace

deps/SimpleProgressBar.hpp

@@ -0,0 +1,256 @@
+/*
+ * MIT License
+ *
+ * Copyright (c) 2024 Robert Myers
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in all
+ * copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ * SOFTWARE.
+ */
+
+/**
+ * @file SimpleProgressBar.h
+ * 
+ * @brief A simple progress bar for console applications.
+ * 
+ * @author Robert Myers
+ * Contact: romyers@umich.edu
+ */
+
+ #pragma once
+
+ #include <iostream>
+ 
+ namespace SimpleProgressBar {
+ 
+     /**
+      * A simple progress bar for console applications.
+      */
+     class ProgressBar {
+ 
+     public:
+ 
+         /**
+          * Constructor.
+          * 
+          * @param totalSteps The total number of progress bar steps.
+          */
+         ProgressBar(unsigned int totalSteps = 100U);
+ 
+         /**
+          * Increments the progress displayed by the progress bar by a given
+          * number of steps (one by default).
+          * 
+          * If increment() brings the current progress above totalSteps,
+          * the current progress will be snapped to totalSteps.
+          * 
+          * @param steps The number of steps by which to increment the displayed
+          * progress.
+          */
+         void increment(unsigned int steps = 1U);
+ 
+         /**
+          * Retrieve the total number of steps in the progress bar.
+          * 
+          * @return The total number of steps
+          */
+         unsigned int getTotalSteps() const;
+ 
+         /**
+          * Prints the progress bar in its current state to a stream.
+          * Prints to std::cout by default.
+          * 
+          * @param out The stream to print to.
+          */
+         void print(std::ostream &out = std::cout) const;
+ 
+         /**
+          * Sets the width of the progress bar display in number of characters.
+          * The width of the progress bar is set to 80 characters by default.
+          * 
+          * REQUIRES: w is at least 2.
+          * 
+          * @param w The desired width of the progress bar, as a number of 
+          * characters.
+          */
+         void setWidth(unsigned int w);
+ 
+         /**
+          * Sets the symbol for the progress bar's left endcap.
+          * Set to '[' by default.
+          * 
+          * @param c The new symbol.
+          */
+         void setLeftEndcapSymbol(char c);
+ 
+         /**
+          * Sets the symbol for the progress bar's right endcap.
+          * Set to ']' by default.
+          * 
+          * @param c The new symbol.
+          */
+         void setRightEndcapSymbol(char c);
+ 
+         /**
+          * Sets the symbol indicating completed progress.
+          * Set to '=' by default.
+          * 
+          * @param c The new symbol.
+          */
+         void setDoneSymbol(char c);
+ 
+         /**
+          * Sets the symbol indicating progress that hasn't been completed yet.
+          * Set to ' ' by default.
+          */
+         void setTodoSymbol(char c);
+ 
+         /**
+          * Sets the progress bar to overwrite the current output line on each
+          * call to print(). This is equivalent to adding the '\r' character to
+          * the end of the printed output.
+          */
+         void enableOverwrite();
+ 
+         /**
+          * Sets the progress bar not to overwrite the current output line on
+          * each call to print(). Successive calls to print will be printed
+          * one after another.
+          */
+         void disableOverwrite();
+ 
+     private:
+ 
+         unsigned int progress;
+         unsigned int totalSteps;
+ 
+         char leftEndcapSymbol;
+         char rightEndcapSymbol;
+         char doneSymbol;
+         char todoSymbol;
+ 
+         unsigned int width;
+ 
+         bool overwrite;
+ 
+     };
+ 
+ }
+ 
+ ///////////////////////////////////////////////////////////////////////////////
+ /////////////////////////////// IMPLEMENTATION ////////////////////////////////
+ ///////////////////////////////////////////////////////////////////////////////
+ 
+ SimpleProgressBar::ProgressBar::ProgressBar(unsigned int totalSteps) : 
+     progress(0),
+     leftEndcapSymbol('['),
+     rightEndcapSymbol(']'),
+     doneSymbol('='),
+     todoSymbol(' '),
+     width(80U),
+     overwrite(true),
+     totalSteps(totalSteps) {}
+ 
+ unsigned int SimpleProgressBar::ProgressBar::getTotalSteps() const {
+ 
+     return totalSteps;
+ 
+ }
+ 
+ void SimpleProgressBar::ProgressBar::setLeftEndcapSymbol(char c) {
+ 
+     leftEndcapSymbol = c;
+ 
+ }
+ 
+ void SimpleProgressBar::ProgressBar::setRightEndcapSymbol(char c) {
+ 
+     rightEndcapSymbol = c;
+ 
+ }
+ 
+ void SimpleProgressBar::ProgressBar::setDoneSymbol(char c) {
+ 
+     doneSymbol = c;
+ 
+ }
+ 
+ void SimpleProgressBar::ProgressBar::setTodoSymbol(char c) {
+ 
+     todoSymbol = c;
+ 
+ }
+ 
+ void SimpleProgressBar::ProgressBar::setWidth(unsigned int w) {
+ 
+     width = w;
+ 
+ }
+ 
+ void SimpleProgressBar::ProgressBar::increment(unsigned int steps) {
+ 
+     progress += steps;
+ 
+     if(progress > totalSteps) {
+ 
+         progress = totalSteps;
+ 
+     }
+ 
+ }
+ 
+ void SimpleProgressBar::ProgressBar::print(std::ostream &out) const {
+ 
+     out << leftEndcapSymbol;
+ 
+     double proportionalProgress = static_cast<double>(progress) 
+                                 / static_cast<double>(totalSteps);
+ 
+     // For nonnegative values, casting to int is equivalent to flooring.
+     unsigned int numBarChars = static_cast<unsigned int>(
+         proportionalProgress * (width - 2)
+     );
+ 
+     for(unsigned int w = 0; w < numBarChars; ++w) {
+ 
+         out << doneSymbol;
+     
+     }
+ 
+     for(unsigned int w = numBarChars; w < width - 2; ++w) {
+ 
+         out << todoSymbol;
+ 
+     }
+ 
+     out << rightEndcapSymbol;
+     out << std::flush;
+     if(overwrite) out << '\r';
+ 
+ }
+ 
+ void SimpleProgressBar::ProgressBar::enableOverwrite() {
+ 
+     overwrite = true;
+ 
+ }
+ 
+ void SimpleProgressBar::ProgressBar::disableOverwrite() {
+ 
+     overwrite = false;
+ 
+ }

geometry-processing/7-MonteCarloGeometryProcessing/c++/MCLaplace.cpp

@@ -37,13 +37,62 @@
 #include "stb_image.h"
 #define STB_IMAGE_WRITE_IMPLEMENTATION
 #include "stb_image_write.h"
+#include "SimpleProgressBar.hpp"
 
 //Global flag to silent verbose messages
 bool silent;
 
 
-int main(int argc, char **argv)
+// Returns the distance and the position of the nearest non-white point
-{
+std::tuple<float, std::tuple<int, int>> nearest_neighbor(std::vector<bool> is_white, int width, int height, int px, int py) {
+  
+  float min_dist = INT_MAX;
+  std::tuple<int, int> nearest_index;
+  for (auto i=0; i < height; i++) {
+    for (auto j=0; j < width; j++) {
+      if (!is_white[i*width+j]) {
+        int dist = (i-px)*(i-px) + (j-py)*(j-py);
+        if (dist < min_dist) {
+          min_dist = dist;
+          nearest_index = {i, j};
+        }
+      }
+    }
+  }
+  return {min_dist, nearest_index};
+}
+
+std::tuple<float, std::tuple<int, int>> random_walk(std::vector<bool> is_white, int width, int height, int px, int py, int cur_depth, int max_depth) {
+  /*
+  * a. Set x0=x
+  * b. For a given point xi, compute the min distance d to ∂Ω
+  * c. Draw a random point xi+1 on ∂B(xi,d) (ball centered at xi with radius d)
+  * d. If xi+1 is close to the boundary (ϵ-close), retrieve the boundary conditions value g(xi+1)
+  * e. Otherwise go to step b.
+  */
+  auto nearest = nearest_neighbor(is_white, width, height, px, py);
+
+  float distance = std::get<0>(nearest);
+  if (distance <= 1) {
+    return nearest;
+  } else if (max_depth < cur_depth) {
+    return {__FLT_MAX__, {px, py}};
+  }
+
+  std::tuple<int, int> near_p = std::get<1>(nearest);
+
+  float r = distance;//*std::rand()/(float)RAND_MAX;
+  float alpha = M_2_PI*std::rand()/(float)RAND_MAX;
+
+  int next_px = px+r*cos(alpha);
+  int next_py = py+r*sin(alpha);
+
+  //std::cout << "("<<px<<", "<<py<<") -"<<distance<<"-> "<<"("<<r*cos(alpha)<<", "<<r*sin(alpha)<<")" << std::endl;
+  
+  return random_walk(is_white, width, height, next_px, next_py, cur_depth+1, max_depth);
+}
+
+int main(int argc, char **argv) {
   CLI::App app{"MCMonteCarlo"};
   std::string sourceImage;
   app.add_option("-i,--input", sourceImage, "Source image")->required()->check(CLI::ExistingFile);;
@@ -69,9 +118,55 @@ int main(int argc, char **argv)
    
   //Main computation
   std::vector<unsigned char> output(width*height*nbChannels);
-  
+  std::vector<bool> is_white(width*height);
-  for(auto i=0; i < width*height*nbChannels; ++i)
+
+  for(auto i=0; i < width*height*nbChannels; ++i) {
     output[i] = source[i];
+  }
+  for (auto i=0; i < height; i++) {
+    for (auto j=0; j < width; j++) {
+      is_white[i*width+j] = (
+        source[nbChannels*(i*width+j)] == 255
+        && source[nbChannels*(i*width+j)+1] == 255
+        && source[nbChannels*(i*width+j)+2] == 255
+      );
+    }
+  }
+  
+
+    
+  SimpleProgressBar::ProgressBar bar(height*width);
+  for (auto i=0; i < height; i++) {
+    for (auto j=0; j < width; j++) {
+      if (is_white[i*width+j]) {
+        /*
+        * a. Set x0=x
+        * b. For a given point xi, compute the min distance d to ∂Ω
+        * c. Draw a random point xi+1 on ∂B(xi,d) (ball centered at xi with radius d)
+        * d. If xi+1 is close to the boundary (ϵ-close), retrieve the boundary conditions value g(xi+1)
+        * e. Otherwise go to step b.
+        */
+       std::vector<float> rescol(nbChannels);
+       for (auto ch=0; ch < nbChannels; ch++)
+        rescol[ch] = 0.;
+       
+       for (auto k=0; k < nbSpp; k++) {
+          auto res = random_walk(is_white, width, height, i, j, 0, nbSteps);
+          float distance = std::get<0>(res);
+          int px = std::get<0>(std::get<1>(res));
+          int py = std::get<1>(std::get<1>(res));
+
+          for (auto ch=0; ch < nbChannels; ch++)
+            rescol[ch] += (distance == __FLT_MAX__ ? 255 : source[nbChannels*(px*width+py)+ch])/(float)nbSpp;
+        }
+        for (auto ch=0; ch < nbChannels; ch++)
+         output[nbChannels*(i*width+j)+ch] = rescol[ch];
+      }
+      bar.increment();
+      bar.print();
+    }
+  }
+  std::cout << std::endl;
   
   
   
@@ -80,8 +175,7 @@ int main(int argc, char **argv)
   //Final export
   if (!silent) std::cout<<"Exporting.."<<std::endl;
   int errcode = stbi_write_png(outputImage.c_str(), width, height, nbChannels, output.data(), nbChannels*width);
-  if (!errcode)
+  if (!errcode) {
-  {
     std::cout<<"Error while exporting the resulting image."<<std::endl;
     exit(errcode);
   }