add of moments extraction from rawls files

main/CMakeLists.txt

@@ -7,6 +7,7 @@ add_executable(rawls_merge_median_incr rawls_merge_median_incr.cpp)
 add_executable(rawls_merge_MON_incr rawls_merge_MON_incr.cpp)
 add_executable(rawls_merge_MON_pct_incr rawls_merge_MON_pct_incr.cpp)
 add_executable(rawls_update rawls_update.cpp)
+add_executable(extract_stats_images extract_stats_images.cpp)
 
 target_link_libraries(rawls_convert LINK_PUBLIC rawls)
 target_link_libraries(rawls_merge_mean LINK_PUBLIC rawls)
@@ -16,6 +17,7 @@ target_link_libraries(rawls_merge_mean_incr_v2 LINK_PUBLIC rawls)
 target_link_libraries(rawls_merge_median_incr LINK_PUBLIC rawls)
 target_link_libraries(rawls_merge_MON_incr LINK_PUBLIC rawls)
 target_link_libraries(rawls_merge_MON_pct_incr LINK_PUBLIC rawls)
+target_link_libraries(extract_stats_images LINK_PUBLIC rawls)
 target_link_libraries(rawls_update LINK_PUBLIC rawls rawls_v1)
 
 set_property(TARGET rawls_merge_mean PROPERTY CXX_STANDARD 17)
@@ -26,4 +28,5 @@ set_property(TARGET rawls_merge_mean_incr_v2 PROPERTY CXX_STANDARD 17)
 set_property(TARGET rawls_merge_median_incr PROPERTY CXX_STANDARD 17)
 set_property(TARGET rawls_merge_MON_incr PROPERTY CXX_STANDARD 17)
 set_property(TARGET rawls_merge_MON_pct_incr PROPERTY CXX_STANDARD 17)
+set_property(TARGET extract_stats_images PROPERTY CXX_STANDARD 17)
 set_property(TARGET rawls_update PROPERTY CXX_STANDARD 17)

main/extract_stats_images.cpp

@@ -0,0 +1,267 @@
+#include <stdio.h>
+#include <string.h>
+#include <sstream>
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <tuple>
+#include <cmath>
+#include <numeric>  
+
+#include <algorithm>
+#include <filesystem>
+
+#include "rawls.h"
+
+struct Point {
+    unsigned x;
+    unsigned y;
+};
+
+struct Tile {
+    Point p1;
+    Point p2;
+};
+
+
+void writeProgress(float progress, bool moveUp = false){
+    int barWidth = 200;
+
+    if (moveUp){
+        // move up line
+        std::cout << "\e[A";
+        std::cout.flush();
+    }
+
+    std::cout << "[";
+    int pos = barWidth * progress;
+    for (int i = 0; i < barWidth; ++i) {
+        if (i < pos) std::cout << "=";
+        else if (i == pos) std::cout << ">";
+        else std::cout << " ";
+    }
+    std::cout << "] " << int(progress * 100.0) << " %\r";
+    std::cout.flush();
+}
+
+float getEstimator(std::string estimator, std::vector<float> values) {
+
+    if (estimator == "median") {
+        
+        std::sort(values.begin(), values.end());
+
+        unsigned size = values.size();
+
+        if (size % 2 == 0)
+        {
+            return (values[size / 2 - 1] + values[size / 2]) / 2;
+        }
+        else 
+        {
+            return values[size / 2];
+        }
+    } else if (estimator == "mean") {
+
+        return std::accumulate(values.begin(), values.end(), 0.0) / values.size(); 
+
+    } else if (estimator == "var") {
+        // Calculate the mean
+        const float mean = std::accumulate(values.begin(), values.end(), 0.0) / values.size();
+
+        // Now calculate the variance
+        auto variance_func = [&mean](float accumulator, const float& val) {
+            return accumulator + pow(val - mean, 2);
+        };
+
+        return std::accumulate(values.begin(), values.end(), 0.0, variance_func) / values.size();
+
+    } else if (estimator == "std") {
+
+        return sqrt(getEstimator("var", values));
+
+    } else if (estimator == "skewness") {
+
+        unsigned size = values.size();
+
+        float mean = getEstimator("mean", values);
+        float std = getEstimator("std", values);
+
+        // Now calculate the sum of pow 3
+        auto order3_func = [&mean, &size](float accumulator, const float& val) {
+            return accumulator + pow(val - mean, 3);
+        };
+
+        float order3 = std::accumulate(values.begin(), values.end(), 0.0, order3_func);
+
+        return order3 / ((size -1 ) * pow(std, 3));
+
+    } else if (estimator == "kurtosis") {
+        
+        unsigned size = values.size();
+
+        float mean = getEstimator("mean", values);
+
+        // Now calculate the sum of pow 4
+        auto order4_func = [&mean, &size](float accumulator, const float& val) {
+            return accumulator + pow(val - mean, 4);
+        };
+
+        float order4 = std::accumulate(values.begin(), values.end(), 0.0, order4_func);
+
+        // Now calculate the sum of pow 2
+        auto order2_func = [&mean, &size](float accumulator, const float& val) {
+            return accumulator + pow(val - mean, 2);
+        };
+
+        float order2 = std::accumulate(values.begin(), values.end(), 0.0, order2_func);
+
+        return size * (order4 / pow(order2, 2));
+    }
+
+    // by default
+    return 0.;
+}
+
+int main(int argc, char *argv[]){
+
+    std::string folderName;
+    std::string estimator;
+    unsigned blockHeight;
+    unsigned blockWidth;
+    std::string outfileName;
+
+    for (int i = 1; i < argc; ++i) {
+        if (!strcmp(argv[i], "--folder") || !strcmp(argv[i], "-folder")) {
+            folderName = argv[++i];
+        } else if (!strcmp(argv[i], "--estimator") || !strcmp(argv[i], "-estimator")) {
+            estimator = argv[++i];
+        } else if (!strcmp(argv[i], "--bwidth") || !strcmp(argv[i], "-bwidth")) {
+            blockHeight = atoi(argv[++i]);
+        } else if (!strcmp(argv[i], "--bheight") || !strcmp(argv[i], "-bheight")) {
+            blockWidth = atoi(argv[++i]);
+        } else if (!strcmp(argv[i], "--outfile") || !strcmp(argv[i], "-outfile")) {
+            outfileName = argv[++i];
+        }
+    }
+
+    std::vector<std::string> imagesPath;
+
+    for (const auto & entry : std::filesystem::directory_iterator(folderName)){
+        std::string imageName = entry.path().string();
+        if (rawls::HasExtension(imageName, ".rawls") || rawls::HasExtension(imageName, ".rawls_20")){
+            imagesPath.push_back(imageName);
+        }
+    }
+
+    std::sort(imagesPath.begin(), imagesPath.end());
+
+    std::tuple<unsigned, unsigned, unsigned> data = rawls::getDimensionsRAWLS(imagesPath.at(0));
+
+    unsigned outputWidth = std::get<0>(data);
+    unsigned outputHeight = std::get<1>(data);
+    unsigned nbChanels = std::get<2>(data);
+
+    // new buffer size as new output buffer image (default 3 channels)
+    float* outputBuffer = new float[outputHeight * outputWidth * nbChanels];
+
+    // get all tiles to apply
+    unsigned nWidth = ceil(outputWidth / (float)blockWidth);
+    unsigned nHeight = ceil(outputHeight / (float)blockHeight);
+
+    std::vector<Tile> tiles;
+
+    for (unsigned i = 0; i < nWidth; i++) {
+        for (unsigned j = 0; j < nHeight; j++) {
+
+            unsigned x1 = i * blockWidth;
+            unsigned y1 = j * blockHeight;
+
+            unsigned x2 = i * blockWidth + blockWidth;
+            unsigned y2 = j * blockHeight + blockHeight;
+
+            x2 = x2 > outputWidth ? outputWidth: x2;
+            y2 = y2 > outputHeight ? outputHeight: y2;
+            
+            Point p1 = {x1, y1};
+            Point p2 = {x2, y2};
+
+            Tile tile = {p1, p2};
+            tiles.push_back(tile);
+        }
+    }
+
+    unsigned nsamples = imagesPath.size();
+    unsigned nloop = tiles.size() * nsamples;
+    unsigned nloopCounter = 0;
+
+    for (unsigned t_index = 0; t_index < tiles.size(); t_index++){
+
+        Tile tile = tiles.at(t_index);
+
+        //std::cout << "Tile: (" << tile.p1.x << ", " << tile.p1.y << ")" << " => " << "(" << tile.p2.x << ", " << tile.p2.y << ")" << std::endl;
+
+        unsigned nvalues = (tile.p2.x - tile.p1.x) * (tile.p2.y - tile.p1.y) * 3;
+
+        std::vector<std::vector<float>> rgbValues(nvalues);
+
+        for (unsigned i = 0; i < nsamples; i++) {
+            
+            try {
+                
+                float* RGBpixels = rawls::getPixelsRAWLS(imagesPath.at(i));
+
+                unsigned index = 0;
+                for (int y = tile.p1.y; y < tile.p2.y; ++y) {
+                    for (int x = tile.p1.x; x < tile.p2.x; ++x) {
+            
+                        rgbValues.at(index).push_back(RGBpixels[3 * (y * outputWidth + x) + 0]);
+                        rgbValues.at(index + 1).push_back(RGBpixels[3 * (y * outputWidth + x) + 1]);
+                        rgbValues.at(index + 2).push_back(RGBpixels[3 * (y * outputWidth + x) + 2]);
+
+                        index += 3;
+                    
+                    }
+                }
+
+                delete RGBpixels;
+
+            } catch(std::exception& e){
+                std::cout << "Error occurs when reading file" << std::endl;
+            }
+
+            // display progress
+            nloopCounter += 1;
+            writeProgress(nloopCounter / (float)nloop);
+        }
+
+        // extract stat and add predicted value into output buffer
+        unsigned index = 0;
+
+        for (int y = tile.p1.y; y < tile.p2.y; ++y) {
+            for (int x = tile.p1.x; x < tile.p2.x; ++x) {
+                
+                outputBuffer[3 * (y * outputWidth + x) + 0] = getEstimator(estimator, rgbValues.at(index + 0));
+                outputBuffer[3 * (y * outputWidth + x) + 1] = getEstimator(estimator, rgbValues.at(index + 1));
+                outputBuffer[3 * (y * outputWidth + x) + 2] = getEstimator(estimator, rgbValues.at(index + 2));
+
+                index += 3;
+            }
+        }
+    }
+
+    // Save here new rawls image
+    std::string comments = rawls::getCommentsRAWLS(imagesPath.at(0));
+
+    bool success = rawls::saveAsRAWLS(outputWidth, outputHeight, nbChanels, comments, outputBuffer, outfileName);
+
+    if (success) {
+        std::cout << "New image saved into " << outfileName << std::endl;
+    }
+    else
+    {
+        std::cout << "Error while saving current image " << outfileName << std::endl;
+    }
+
+    delete outputBuffer;
+
+}