Implement write_matrix_as_ppm() function in utils module

   this program in the COPYING file. If not, see: <https://www.gnu.org/licenses/>.

 */

/*! \file utils.h

/**

 * \file utils.h

 *

 * \brief Definition of auxiliary code utilities.

 */

#ifndef INCLUDE_UTILS_H_

#define INCLUDE_UTILS_H_

/*! \brief Obtain an estimate of the RAM overhead factor for optimization.

/**

 * \brief Obtain an estimate of the RAM overhead factor for optimization.

 *

 * Code speed-up optimization usually comes at the cost of increased memory

 * requirements. While this should not generally be a serious issue, it can

 */

double get_ram_overhead();

/*! \brief Write a double complex matrix to a text file.

/**

 * \brief Write a double complex matrix to a text file.

 *

 * \param af: `VirtualAsciiFile *` Pointer to an existing VirtualAsciiFile.

 * \param mat: `dcomplex **` Pointer to the matrix.

 * \param rows: `int` Number of rows in the matrix.

 * \param columns: `int` Number of columns in the matrix.

 * \param format: `const string&` Format of the line (default is \" %5d %5d (%17.8lE,%17.8lE)\n\")

 * \param first_index: `int` Index of the first element (default is 1, i.e. base 1 FORTRAN array notation)

 * \param af: `VirtualAsciiFile *` Pointer to an existing VirtualAsciiFile. [IN]

 * \param mat: `dcomplex **` Pointer to the matrix. [IN]

 * \param rows: `int` Number of rows in the matrix. [IN]

 * \param columns: `int` Number of columns in the matrix. [IN]

 * \param format: `const string&` Format of the line (default is \" %5d %5d (%17.8lE,%17.8lE)\n\"). [IN]

 * \param first_index: `int` Index of the first element (default is 1, i.e. base 1 FORTRAN array notation). [IN]

 * \return result: `int` An exit code (0 if successful).

 */

int write_dcomplex_matrix(

			  VirtualAsciiFile *af, dcomplex **mat, int rows,

			  int columns, const std::string& format=" %5d %5d (%17.8lE,%17.8lE)\n",

			  int first_index=1

  VirtualAsciiFile *af, dcomplex **mat, int rows, int columns,

  const std::string& format=" %5d %5d (%17.8lE,%17.8lE)\n", int first_index=1

);

/**

 * \brief Draw a PPM representation of a matrix.

 *

 * \param A: `const dcomplex *` Vector representation of the matrix. [IN]

 * \param m: `int` Number of rows in the matrix. [IN]

 * \param n: `int` Number of columns in the matrix. [IN]

 * \param file_name: `const string&` Reference to a string with the name of the output file. [IN]

 * \param mode: `const string&` Reference to a string for the value to be represented ("MAG" to

 * draw magnitudes, "RE" to draw real parts, "IM" to draw imaginary parts. Optional, default is

 * "MAG"). [IN]

 * \param row_bin: `int` Number of rows to bin together (optional, default is 1). [IN]

 * \param col_bin: `int` Number of columns to bin together (optional, default is 1). [IN]

 * \return result: `int` An exit code (0 if successful).

 */

int write_matrix_as_ppm(

  const dcomplex *A, int m, int n, const std::string& file_name, const std::string& mode="MAG",

  int row_bin=1, int col_bin=1

);

#endif

 * \brief Implementation of auxiliary code utilities.

 */

#include <fstream>

#include <hdf5.h>

#ifndef INCLUDE_TYPES_H_

#include "../include/types.h"

  }

  return result;

}

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

) {

  double (*cabs)(const dcomplex& z) = [] (const dcomplex& z) -> double { return dcabs(z); };

  double (*function)(const dcomplex& z) = cabs;

  if (mode.compare("REAL") == 0) {

    function = ℜ

  } else if (mode.compare("IMAG") == 0) {

    function = &imag;

  }

  np_int mn = m * n;

  np_int ppm_height = (row_bin > 1) ? m / row_bin : m;

  np_int ppm_width = (col_bin > 1) ? n / col_bin : n;

  const int height_spare = m % row_bin;

  const int width_spare = n % col_bin;

  if (height_spare > 0) ++ppm_height;

  if (width_spare > 0) ++ppm_width;

  np_int ppm_size = ppm_width * ppm_height;

  int bin_size = row_bin * col_bin;

  ofstream ofs(file_name, ios::binary);

  // Header PPM: P5 = binary grayscale, width, height, max value

  ofs << "P5\n" << ppm_height << " " << ppm_width << "\n255\n";

  double max_val = 0.0;

  double *logs = new double[ppm_size];

  for (np_int pi = 0; pi < ppm_size; ++pi) {

    int this_bin_size = bin_size;

    np_int vi = pi / ppm_width;

    np_int vj = pi % ppm_width;

    np_int first_i = vi * row_bin;

    np_int last_i = first_i + row_bin;

    if (last_i > m) {

      last_i = m;

      this_bin_size /= row_bin;

      this_bin_size *= height_spare;

    }

    np_int first_j = vj * col_bin;

    np_int last_j = first_j + col_bin;

    if (last_j > n) {

      last_j = n;

      this_bin_size /= col_bin;

      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++) {

	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;

  }

  // Write normalized pixel values.

  for (int i = 0; i < mn; ++i) {

    // 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;

    ofs.put(pixel);

  }

  ofs.close();

  delete[] logs;

}