Commit 3445f26d authored by Giovanni La Mura's avatar Giovanni La Mura
Browse files

Complete the implementation of write_matrix_as_ppm()

parent 9dc63b17
Loading
Loading
Loading
Loading
+18 −8
Original line number Diff line number Diff line
@@ -78,6 +78,7 @@ int write_matrix_as_ppm(
  const dcomplex *A, np_int m, np_int n, const std::string& file_name, const std::string& mode,
  int row_bin, int col_bin
) {
  int result = 0;
  double (*cabs)(const dcomplex& z) = [] (const dcomplex& z) -> double { return dcabs(z); };
  double (*function)(const dcomplex& z) = cabs;
  if (mode.compare("REAL") == 0) {
@@ -97,11 +98,13 @@ int write_matrix_as_ppm(
  
  ofstream ofs(file_name, ios::binary);
  // Header PPM: P5 = binary grayscale, width, height, max value
  ofs << "P5\n" << ppm_height << " " << ppm_width << "\n255\n";
  ofs << "P5\n" << ppm_width << " " << ppm_height << "\n255\n";

  double max_val = 0.0;
  double max_val = 0.0, min_val = 1.0e+09;
  double *logs = new double[ppm_size];

#pragma omp parallel for reduction(max: max_val) \
  reduction(min: min_val)
  for (np_int pi = 0; pi < ppm_size; ++pi) {
    int this_bin_size = bin_size;
    np_int vi = pi / ppm_width;
@@ -121,24 +124,31 @@ int write_matrix_as_ppm(
      this_bin_size *= width_spare;
    }
    double value = 0.0;
    for (np_int ai = first_i; ai < last_i; ai++) {
      for (np_int aj = first_j; aj < last_j; aj++) {
    np_int asize = (last_i - first_i) * (last_j - first_j);
    for (np_int aij = 0; aij < asize; ++aij) {
      np_int ai = first_i + aij / (last_j - first_j);
      np_int aj = first_j + aij % (last_j - first_j);
      value += function(A[n * ai + aj]);
    }
    }
    value /= this_bin_size;
    logs[ppm_width * vi + vj] = (value > 0.0) ? log10(value) : -1.0;
    if (value > max_val) max_val = value;
    if (value < min_val && min_val > 0.0) min_val = value;
  }
  if (min_val <= 0.0) result = 1; // Matrix contains non-positive elements: log-scale corruption!
  double lmax_val = (max_val > 0.0) ? std::log10(max_val) : 0.0;
  double lmin_val = (min_val > 0.0) ? std::log10(min_val) : 0.0;

  // Write normalized pixel values.
  for (int i = 0; i < mn; ++i) {
  for (int i = 0; i < ppm_size; ++i) {
    double pix_val = (logs[i] == 1000.0) ? 0 : 255.0 * (logs[i] - lmin_val) / (lmax_val - lmin_val);
    // Avoid 0-division errors for empty matrices
    unsigned char pixel = (max_val > 0.0) ? 
      static_cast<unsigned char>((1.0 + logs[i]) * 255.0) : 0;
      static_cast<unsigned char>(pix_val) : 0;
    ofs.put(pixel);
  }

  ofs.close();
  delete[] logs;
  return result;
}