main.cpp 5.53 KB
Newer Older
1
#include <algorithm>
Lars Kuehne's avatar
Lars Kuehne committed
2 3
#include <cassert>
#include <chrono>
plgruener's avatar
plgruener committed
4
#include <cmath>
Lars Kuehne's avatar
Lars Kuehne committed
5
#include <fstream>
6
#include <iostream>
Lars Kuehne's avatar
Lars Kuehne committed
7
#include <random>
plgruener's avatar
plgruener committed
8
#include <sstream>
Christoph Saffer's avatar
Christoph Saffer committed
9 10
#include <vector>

11
#include <boost/program_options.hpp>
Christoph Saffer's avatar
Christoph Saffer committed
12

plgruener's avatar
plgruener committed
13
#include "colormaps.hpp"
14
#include "compute.hpp"
plgruener's avatar
plgruener committed
15

plgruener's avatar
plgruener committed
16 17


18
int main(int argc, char* argv[]) {
plgruener's avatar
plgruener committed
19
  // get arguments from CLI
20 21
  // these can be input by user
  int num_iterations;
plgruener's avatar
plgruener committed
22
  float threshold;
23 24 25 26 27 28 29 30 31 32 33 34 35
  float alphamin;
  float alphamax;
  int alpha_num_intervals;
  float betamin;
  float betamax;
  int beta_num_intervals;

  namespace po = boost::program_options;
  try {
    po::options_description desc("Options");
    desc.add_options()
      ("help", "Help message")
      ("iterations,n", po::value<int>(&num_iterations)->default_value(100), "Number of iterations")
plgruener's avatar
plgruener committed
36
      ("threshold,s", po::value<float>(&threshold)->default_value(1), " Threshold above that computation is stopped")
37 38 39 40 41 42 43 44 45 46 47 48 49 50 51
      ("amin,a", po::value<float>(&alphamin)->default_value(0), " α lower bound")
      ("amax,A", po::value<float>(&alphamax)->default_value(1), " α upper bound")
      ("alphas,w", po::value<int>(&alpha_num_intervals)->default_value(100), "α resolution/width of image")
      ("betas,h", po::value<int>(&beta_num_intervals)->default_value(100), "β resolution/height of image")
      ("bmin,b", po::value<float>(&betamin)->default_value(0), " β lower bound")
      ("bmax,B", po::value<float>(&betamax)->default_value(1), " β upper bound")
      ;

    po::variables_map vm;
    po::store(po::parse_command_line(argc, argv, desc), vm);

    if (vm.count("help")) {
      std::cout << desc << std::endl;
      return 0;
    }
52

53 54 55 56 57 58 59 60 61
    po::notify(vm);

    // check if our integers are >0, else throw invalid-argument-error
    if ( (num_iterations < 1) || (alpha_num_intervals < 1) || (beta_num_intervals <1) ) {
      throw po::validation_error(po::validation_error::invalid_option_value);
    }

  } catch (po::error &e) {
    std::cerr << e.what() << std::endl;
plgruener's avatar
plgruener committed
62 63
    return -1;
  }
64 65

  // these are computed
66
  int alpha_num_params = alpha_num_intervals + 1;
Lars Kuehne's avatar
Lars Kuehne committed
67
  float alpha_interval_size = (alphamax - alphamin) / (alpha_num_intervals);
68
  int beta_num_params = beta_num_intervals + 1;
Lars Kuehne's avatar
Lars Kuehne committed
69
  float beta_interval_size = (betamax - betamin) / (beta_num_intervals);
plgruener's avatar
plgruener committed
70

71
  // fill Parametervectors:
Christoph Saffer's avatar
Christoph Saffer committed
72
  aligned_vector<float> alphas(alpha_num_params);
Lars Kuehne's avatar
Lars Kuehne committed
73 74 75 76
  float* alphasp = alphas.data();
#pragma omp simd aligned(alphasp : 64)
  for (int i = 0; i < alpha_num_params; i++) {
    alphasp[i] = alphamin + i * alpha_interval_size;
77
  }
Christoph Saffer's avatar
Christoph Saffer committed
78
  aligned_vector<float> betas(beta_num_params);
Lars Kuehne's avatar
Lars Kuehne committed
79 80 81 82
  float* betasp = betas.data();
#pragma omp simd aligned(betasp : 64)
  for (int i = 0; i < beta_num_params; i++) {
    betasp[i] = betamin + i * beta_interval_size;
83 84
  }

Lars Kuehne's avatar
Lars Kuehne committed
85 86 87
  aligned_vector<float> x_start = {0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9};
  // aligned_vector<float> y_start = {0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9};
  aligned_vector<float> y_start = {0, 0, 0, 0, 0, 0, 0, 0, 0};
88

89
  aligned_vector<float> result(alpha_num_params * beta_num_params);
90
  aligned_vector<int> colors(alpha_num_params * beta_num_params);
91
  aligned_vector<int> colors_rgb(3 * alpha_num_params * beta_num_params);
Lars Kuehne's avatar
Lars Kuehne committed
92

plgruener's avatar
plgruener committed
93 94
  auto time_start = std::chrono::system_clock::now();

95
#pragma omp parallel for schedule(dynamic)
Lars Kuehne's avatar
Lars Kuehne committed
96 97
  for (int a = 0; a < alpha_num_params; a++) {
    for (int b = beta_num_params - 1; b >= 0; b--) {
98 99
      result[(beta_num_params - b - 1) * alpha_num_params + a] =
          compute(alphas[a], betas[b], x_start, y_start, num_iterations, threshold);
100
    }
101
  }
plgruener's avatar
plgruener committed
102
  auto time_end = std::chrono::system_clock::now();
Lars Kuehne's avatar
Lars Kuehne committed
103 104
  float const elapsed_seconds =
      std::chrono::duration<float>(time_end - time_start).count();
plgruener's avatar
plgruener committed
105
  std::cout << "TIME: " << elapsed_seconds << std::endl;
plgruener's avatar
plgruener committed
106

107 108 109 110 111 112 113 114 115 116
  // Output result into .csv
  std::string file_result = "result.csv";
  std::ofstream ostrm_csv(file_result);
  ostrm_csv << "alpha beta value\n";
  for (int a = 0; a < alpha_num_params; a++) {
    for (int b = beta_num_params - 1; b >= 0; b--) {
      ostrm_csv << alphas[a] << ' ' << betas[b] << ' ' << result[(beta_num_params - b - 1) * alpha_num_params + a] << std::endl;
    }
  }

117
  // transform floats into grayscale-colors:
118
  for (int i = 0; i < alpha_num_params*beta_num_params; ++i) {
plgruener's avatar
plgruener committed
119
    if (result[i] > threshold) {
120
      colors[i] = 255;
121 122 123
      colors_rgb[3*i] = 255;
      colors_rgb[3*i+1] = 255;
      colors_rgb[3*i+2] = 255;
124
    } else {
plgruener's avatar
plgruener committed
125
      colors[i] = (int) 254 * result[i]/threshold;
plgruener's avatar
plgruener committed
126
      // RGB color gradient: viridis from matplotlib
plgruener's avatar
plgruener committed
127
      int idx = floor(255*result[i]/threshold);
plgruener's avatar
plgruener committed
128 129 130 131 132
      
      int r = floor(255*viridis[idx][0]);
      int g = floor(255*viridis[idx][1]);
      int b = floor(255*viridis[idx][2]);

133 134 135
      colors_rgb[3*i] = r;  // red
      colors_rgb[3*i+1] = g;// green
      colors_rgb[3*i+2] = b;  // blue
136 137 138 139 140 141
    }
  }

  // Output pixel vector into PGM File
  std::string filename = "picture.pgm";
  std::ofstream ostrm(filename);
142
  ostrm << "P5" << '\n';
143 144
  ostrm << alpha_num_params << ' ' << beta_num_params << '\n';
  ostrm << 255 << '\n';  // max. gray value
145
  // 1 byte per pixel
146 147 148 149 150 151 152 153 154 155 156 157 158 159
  for (int i=0;i<colors.size();++i){
  ostrm.write(reinterpret_cast<char*>(&colors[i]),1);//colors.size());
  }
  // Output pixel vector into PPM File
  std::string filename_rgb = "picture_rgb.ppm";
  std::ofstream ostrm_rgb(filename_rgb);
  ostrm_rgb << "P6" << '\n';
  ostrm_rgb << alpha_num_params << ' ' << beta_num_params << '\n';
  ostrm_rgb << 255 << '\n';  // max. gray value
  // 3 byte per pixel
  for (int i=0;i<colors_rgb.size();++i){
    // das geht theoretisch auch mit dem ganzen Array, aber praktisch nicht?!
  ostrm_rgb.write(reinterpret_cast<char*>(&colors_rgb[i]),1);//colors_rgb.size());
  }
plgruener's avatar
plgruener committed
160
}