extract_stats_images_all_reduced.cpp 17 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530
  1. #include <stdio.h>
  2. #include <string.h>
  3. #include <sstream>
  4. #include <iostream>
  5. #include <fstream>
  6. #include <vector>
  7. #include <tuple>
  8. #include <cmath>
  9. #include <numeric>
  10. #include <map>
  11. #include <algorithm>
  12. #include <filesystem>
  13. #include <unistd.h>
  14. #include <bits/stdc++.h>
  15. #include <iostream>
  16. #include <sys/stat.h>
  17. #include <sys/types.h>
  18. #include "rawls.h"
  19. struct Point {
  20. unsigned x;
  21. unsigned y;
  22. };
  23. struct Tile {
  24. Point p1;
  25. Point p2;
  26. };
  27. std::vector<std::string> split(const std::string& s, char delimiter)
  28. {
  29. std::vector<std::string> tokens;
  30. std::string token;
  31. std::istringstream tokenStream(s);
  32. while (std::getline(tokenStream, token, delimiter))
  33. {
  34. tokens.push_back(token);
  35. }
  36. return tokens;
  37. }
  38. void writeProgress(float progress, bool moveUp = false){
  39. int barWidth = 150;
  40. if (moveUp){
  41. // move up line
  42. std::cout << "\e[A";
  43. std::cout.flush();
  44. }
  45. std::cout << "[";
  46. int pos = barWidth * progress;
  47. for (int i = 0; i < barWidth; ++i) {
  48. if (i < pos) std::cout << "=";
  49. else if (i == pos) std::cout << ">";
  50. else std::cout << " ";
  51. }
  52. std::cout << "] " << int(progress * 100.0) << " %\r";
  53. std::cout.flush();
  54. }
  55. float getMedian(std::vector<float> &values) {
  56. std::sort(values.begin(), values.end());
  57. unsigned size = values.size();
  58. if (size % 2 == 0)
  59. {
  60. return (values[size / 2 - 1] + values[size / 2]) / 2;
  61. }
  62. else
  63. {
  64. return values[size / 2];
  65. }
  66. }
  67. float getVariance(float &mean, std::vector<float> &values){
  68. // Now calculate the variance
  69. auto variance_func = [&mean](float accumulator, const float& val) {
  70. return accumulator + pow(val - mean, 2);
  71. };
  72. return std::accumulate(values.begin(), values.end(), 0.0, variance_func) / values.size();
  73. }
  74. float getSkewness(float &mean, float &std, std::vector<float> &values) {
  75. unsigned size = values.size();
  76. // Now calculate the sum of pow 3
  77. auto order3_func = [&mean, &std](float accumulator, const float& val) {
  78. return accumulator + pow((val - mean) / std, 3);
  79. };
  80. float order3 = std::accumulate(values.begin(), values.end(), 0.0, order3_func);
  81. return order3 / size;
  82. }
  83. float getKurtosis(float &mean, float &std, std::vector<float> &values){
  84. unsigned size = values.size();
  85. // Now calculate the sum of pow 4
  86. auto order4_func = [&mean, &std](float accumulator, const float& val) {
  87. return accumulator + pow((val - mean) / std, 4);
  88. };
  89. float order4 = std::accumulate(values.begin(), values.end(), 0.0, order4_func);
  90. return order4 / size;
  91. }
  92. float getMode(std::vector<float> &values) {
  93. std::vector<float> pvalues;
  94. for (unsigned i = 0; i < values.size(); i++){
  95. pvalues.push_back(roundf(values.at(i) * 100) / 100.0);
  96. }
  97. typedef std::map<float,unsigned int> CounterMap;
  98. CounterMap counts;
  99. for (int i = 0; i < pvalues.size(); ++i)
  100. {
  101. CounterMap::iterator it(counts.find(pvalues[i]));
  102. if (it != counts.end()){
  103. it->second++;
  104. } else {
  105. counts[pvalues[i]] = 1;
  106. }
  107. }
  108. // Create a map iterator and point to beginning of map
  109. std::map<float, unsigned int>::iterator it = counts.begin();
  110. unsigned noccurences = 0;
  111. float modeValue = 0.;
  112. // Iterate over the map using Iterator till end.
  113. while (it != counts.end())
  114. {
  115. // Accessing KEY from element pointed by it.
  116. float potentialMode = it->first;
  117. // Accessing VALUE from element pointed by it.
  118. unsigned count = it->second;
  119. if (count > noccurences) {
  120. noccurences = count;
  121. modeValue = potentialMode;
  122. }
  123. // Increment the Iterator to point to next entry
  124. it++;
  125. }
  126. return modeValue;
  127. }
  128. // float getEstimator(std::string estimator, std::vector<float> values) {
  129. // // another version of scripts in order to quick compute data
  130. // if (estimator == "median") {
  131. // std::sort(values.begin(), values.end());
  132. // unsigned size = values.size();
  133. // if (size % 2 == 0)
  134. // {
  135. // return (values[size / 2 - 1] + values[size / 2]) / 2;
  136. // }
  137. // else
  138. // {
  139. // return values[size / 2];
  140. // }
  141. // } else if (estimator == "mean") {
  142. // return std::accumulate(values.begin(), values.end(), 0.0) / values.size();
  143. // } else if (estimator == "var") {
  144. // // Calculate the mean
  145. // const float mean = std::accumulate(values.begin(), values.end(), 0.0) / values.size();
  146. // // Now calculate the variance
  147. // auto variance_func = [&mean](float accumulator, const float& val) {
  148. // return accumulator + pow(val - mean, 2);
  149. // };
  150. // return std::accumulate(values.begin(), values.end(), 0.0, variance_func) / values.size();
  151. // } else if (estimator == "std") {
  152. // return sqrt(getEstimator("var", values));
  153. // } else if (estimator == "skewness") {
  154. // unsigned size = values.size();
  155. // float mean = getEstimator("mean", values);
  156. // float std = getEstimator("std", values);
  157. // // Now calculate the sum of pow 3
  158. // auto order3_func = [&mean, &std](float accumulator, const float& val) {
  159. // return accumulator + pow((val - mean) / std, 3);
  160. // };
  161. // float order3 = std::accumulate(values.begin(), values.end(), 0.0, order3_func);
  162. // return order3 / size;
  163. // } else if (estimator == "kurtosis") {
  164. // unsigned size = values.size();
  165. // float mean = getEstimator("mean", values);
  166. // float std = getEstimator("std", values);
  167. // // Now calculate the sum of pow 4
  168. // auto order4_func = [&mean, &std](float accumulator, const float& val) {
  169. // return accumulator + pow((val - mean) / std, 4);
  170. // };
  171. // float order4 = std::accumulate(values.begin(), values.end(), 0.0, order4_func);
  172. // return order4 / size;
  173. // } else if (estimator == "mode") {
  174. // std::vector<float> pvalues;
  175. // for (unsigned i = 0; i < values.size(); i++){
  176. // pvalues.push_back(roundf(values.at(i) * 100) / 100.0);
  177. // }
  178. // typedef std::map<float,unsigned int> CounterMap;
  179. // CounterMap counts;
  180. // for (int i = 0; i < pvalues.size(); ++i)
  181. // {
  182. // CounterMap::iterator it(counts.find(pvalues[i]));
  183. // if (it != counts.end()){
  184. // it->second++;
  185. // } else {
  186. // counts[pvalues[i]] = 1;
  187. // }
  188. // }
  189. // // Create a map iterator and point to beginning of map
  190. // std::map<float, unsigned int>::iterator it = counts.begin();
  191. // unsigned noccurences = 0;
  192. // float modeValue = 0.;
  193. // // Iterate over the map using Iterator till end.
  194. // while (it != counts.end())
  195. // {
  196. // // Accessing KEY from element pointed by it.
  197. // float potentialMode = it->first;
  198. // // Accessing VALUE from element pointed by it.
  199. // unsigned count = it->second;
  200. // if (count > noccurences) {
  201. // noccurences = count;
  202. // modeValue = potentialMode;
  203. // }
  204. // // Increment the Iterator to point to next entry
  205. // it++;
  206. // }
  207. // return modeValue;
  208. // }
  209. // // by default
  210. // return 0.;
  211. // }
  212. int main(int argc, char *argv[]){
  213. std::string folderName;
  214. std::vector<std::string> estimators = {"mean", "median", "var", "std", "skewness", "kurtosis", "mode"};
  215. unsigned blockHeight;
  216. unsigned blockWidth;
  217. unsigned nfiles = 10000;
  218. std::string outputFolder;
  219. for (int i = 1; i < argc; ++i) {
  220. if (!strcmp(argv[i], "--folder") || !strcmp(argv[i], "-folder")) {
  221. folderName = argv[++i];
  222. } else if (!strcmp(argv[i], "--bwidth") || !strcmp(argv[i], "-bwidth")) {
  223. blockHeight = atoi(argv[++i]);
  224. } else if (!strcmp(argv[i], "--bheight") || !strcmp(argv[i], "-bheight")) {
  225. blockWidth = atoi(argv[++i]);
  226. } else if (!strcmp(argv[i], "--output") || !strcmp(argv[i], "-output")) {
  227. outputFolder = argv[++i];
  228. } else if (!strcmp(argv[i], "--nfiles") || !strcmp(argv[i], "-nfiles")) {
  229. nfiles = atoi(argv[++i]);
  230. }
  231. }
  232. // create outputs directory
  233. mkdir(outputFolder.c_str(), 0755);
  234. auto elements = split(folderName, '/');
  235. std::string sceneName = elements.at(elements.size() - 1);
  236. for (int i = 0; i < estimators.size(); i++) {
  237. mkdir((outputFolder + "/" + estimators[i]).c_str(), 0755);
  238. mkdir((outputFolder + "/" + estimators[i] + "/" + sceneName).c_str(), 0755);
  239. }
  240. // get all files path
  241. std::vector<std::string> imagesPath;
  242. for (const auto & entry : std::filesystem::directory_iterator(folderName)){
  243. std::string imageName = entry.path().string();
  244. if (rawls::HasExtension(imageName, ".rawls") || rawls::HasExtension(imageName, ".rawls_20")){
  245. imagesPath.push_back(imageName);
  246. }
  247. }
  248. if (imagesPath.size() != nfiles) {
  249. return 0;
  250. }
  251. std::sort(imagesPath.begin(), imagesPath.end());
  252. std::tuple<unsigned, unsigned, unsigned> data = rawls::getDimensionsRAWLS(imagesPath.at(0));
  253. unsigned outputWidth = std::get<0>(data);
  254. unsigned outputHeight = std::get<1>(data);
  255. unsigned nbChanels = std::get<2>(data);
  256. std::vector<float*> outputBuffers;
  257. std::vector<std::string> outputFiles;
  258. std::vector<std::string> selectedEstimators;
  259. // new buffer size as new output buffer image (default 3 channels)
  260. for (int i = 0; i < estimators.size(); i++) {
  261. std::string outputFile = outputFolder + "/" + estimators[i] + "/" + sceneName + "/" + sceneName + ".rawls";
  262. std::ifstream ifile;
  263. ifile.open(outputFile);
  264. if(!ifile) {
  265. // create new buffer entry
  266. selectedEstimators.push_back(estimators[i]);
  267. outputFiles.push_back(outputFile);
  268. outputBuffers.push_back(new float[outputHeight * outputWidth * nbChanels]);
  269. } else {
  270. ifile.close();
  271. }
  272. }
  273. // get all tiles to apply
  274. unsigned nWidth = ceil(outputWidth / (float)blockWidth);
  275. unsigned nHeight = ceil(outputHeight / (float)blockHeight);
  276. std::vector<Tile> tiles;
  277. for (unsigned i = 0; i < nWidth; i++) {
  278. for (unsigned j = 0; j < nHeight; j++) {
  279. unsigned x1 = i * blockWidth;
  280. unsigned y1 = j * blockHeight;
  281. unsigned x2 = i * blockWidth + blockWidth;
  282. unsigned y2 = j * blockHeight + blockHeight;
  283. x2 = x2 > outputWidth ? outputWidth: x2;
  284. y2 = y2 > outputHeight ? outputHeight: y2;
  285. Point p1 = {x1, y1};
  286. Point p2 = {x2, y2};
  287. Tile tile = {p1, p2};
  288. tiles.push_back(tile);
  289. }
  290. }
  291. unsigned nsamples = imagesPath.size();
  292. unsigned nloop = tiles.size() * nsamples;
  293. unsigned nloopCounter = 0;
  294. for (unsigned t_index = 0; t_index < tiles.size(); t_index++){
  295. Tile tile = tiles.at(t_index);
  296. //std::cout << "Tile: (" << tile.p1.x << ", " << tile.p1.y << ")" << " => " << "(" << tile.p2.x << ", " << tile.p2.y << ")" << std::endl;
  297. unsigned nvalues = (tile.p2.x - tile.p1.x) * (tile.p2.y - tile.p1.y) * 3;
  298. std::vector<std::vector<float>> rgbValues(nvalues);
  299. for (unsigned i = 0; i < nsamples; i++) {
  300. try {
  301. float* RGBpixels = rawls::getPixelsRAWLS(imagesPath.at(i));
  302. std::cout << "Read image n°" << i << " / " << nsamples << " for tile n°" << t_index << " / " << tiles.size() << std::endl;
  303. unsigned index = 0;
  304. for (int y = tile.p1.y; y < tile.p2.y; ++y) {
  305. for (int x = tile.p1.x; x < tile.p2.x; ++x) {
  306. rgbValues.at(index).push_back(RGBpixels[3 * (y * outputWidth + x) + 0]);
  307. rgbValues.at(index + 1).push_back(RGBpixels[3 * (y * outputWidth + x) + 1]);
  308. rgbValues.at(index + 2).push_back(RGBpixels[3 * (y * outputWidth + x) + 2]);
  309. index += 3;
  310. }
  311. }
  312. delete RGBpixels;
  313. } catch(std::exception& e){
  314. std::cout << "Error occurs when reading file" << std::endl;
  315. }
  316. // display progress
  317. nloopCounter += 1;
  318. writeProgress(nloopCounter / (float)nloop);
  319. }
  320. // for (int i = 0; i < outputFiles.size(); i++) {
  321. // extract stat and add predicted value into output buffer
  322. unsigned index = 0;
  323. for (int y = tile.p1.y; y < tile.p2.y; ++y) {
  324. for (int x = tile.p1.x; x < tile.p2.x; ++x) {
  325. // Here we will compute each estimator in specific order
  326. // => {"median", "var", "std", "skewness", "kurtosis", "mode"}
  327. auto rvalues = rgbValues.at(index + 0);
  328. auto gvalues = rgbValues.at(index + 1);
  329. auto bvalues = rgbValues.at(index + 2);
  330. // Index [0] : MEAN
  331. float rmean = std::accumulate(rvalues.begin(), rvalues.end(), 0.0) / rvalues.size();
  332. float gmean = std::accumulate(gvalues.begin(), gvalues.end(), 0.0) / gvalues.size();
  333. float bmean = std::accumulate(bvalues.begin(), bvalues.end(), 0.0) / bvalues.size();
  334. outputBuffers.at(0)[3 * (y * outputWidth + x) + 0] = rmean;
  335. outputBuffers.at(0)[3 * (y * outputWidth + x) + 1] = gmean;
  336. outputBuffers.at(0)[3 * (y * outputWidth + x) + 2] = bmean;
  337. // Index [1] : MEDIAN
  338. float rmedian = getMedian(rvalues);
  339. float gmedian = getMedian(gvalues);
  340. float bmedian = getMedian(bvalues);
  341. outputBuffers.at(1)[3 * (y * outputWidth + x) + 0] = rmedian;
  342. outputBuffers.at(1)[3 * (y * outputWidth + x) + 1] = gmedian;
  343. outputBuffers.at(1)[3 * (y * outputWidth + x) + 2] = bmedian;
  344. // Index [2] : VARIANCE
  345. float rvariance = getVariance(rmean, rvalues);
  346. float gvariance = getVariance(gmean, gvalues);
  347. float bvariance = getVariance(bmean, bvalues);
  348. outputBuffers.at(2)[3 * (y * outputWidth + x) + 0] = rvariance;
  349. outputBuffers.at(2)[3 * (y * outputWidth + x) + 1] = gvariance;
  350. outputBuffers.at(2)[3 * (y * outputWidth + x) + 2] = bvariance;
  351. // Index [3] : STD
  352. float rstd = sqrt(rvariance);
  353. float gstd = sqrt(gvariance);
  354. float bstd = sqrt(bvariance);
  355. outputBuffers.at(3)[3 * (y * outputWidth + x) + 0] = rstd;
  356. outputBuffers.at(3)[3 * (y * outputWidth + x) + 1] = gstd;
  357. outputBuffers.at(3)[3 * (y * outputWidth + x) + 2] = bstd;
  358. // Index [4] : SKEWNESS
  359. float rskew = getSkewness(rmean, rstd, rvalues);
  360. float gskew = getSkewness(gmean, gstd, gvalues);
  361. float bskew = getSkewness(bmean, bstd, bvalues);
  362. outputBuffers.at(4)[3 * (y * outputWidth + x) + 0] = rskew;
  363. outputBuffers.at(4)[3 * (y * outputWidth + x) + 1] = gskew;
  364. outputBuffers.at(4)[3 * (y * outputWidth + x) + 2] = bskew;
  365. // Index [5] : KURTOSIS
  366. float rkurtosis = getKurtosis(rmean, rstd, rvalues);
  367. float gkurtosis = getKurtosis(gmean, gstd, gvalues);
  368. float bkurtosis = getKurtosis(bmean, bstd, bvalues);
  369. outputBuffers.at(5)[3 * (y * outputWidth + x) + 0] = rkurtosis;
  370. outputBuffers.at(5)[3 * (y * outputWidth + x) + 1] = gkurtosis;
  371. outputBuffers.at(5)[3 * (y * outputWidth + x) + 2] = bkurtosis;
  372. // Index [6] : MODE (TODO: check computation time and if very necessary)
  373. float rmode = getMode(rvalues);
  374. float gmode = getMode(gvalues);
  375. float bmode = getMode(bvalues);
  376. outputBuffers.at(6)[3 * (y * outputWidth + x) + 0] = rmode;
  377. outputBuffers.at(6)[3 * (y * outputWidth + x) + 1] = gmode;
  378. outputBuffers.at(6)[3 * (y * outputWidth + x) + 2] = bmode;
  379. index += 3;
  380. }
  381. }
  382. // }
  383. }
  384. // Save here new rawls image
  385. std::string comments = rawls::getCommentsRAWLS(imagesPath.at(0));
  386. for (int i = 0; i < outputFiles.size(); i++) {
  387. // construct specific outfile name
  388. bool success = rawls::saveAsRAWLS(outputWidth, outputHeight, nbChanels, comments, outputBuffers[i], outputFiles[i]);
  389. if (success) {
  390. std::cout << "New image saved into " << outputFiles[i] << std::endl;
  391. }
  392. else
  393. {
  394. std::cout << "Error while saving current image " << outputFiles[i] << std::endl;
  395. }
  396. delete outputBuffers[i];
  397. }
  398. }