#ifndef _IMAGE_HPP_
#define _IMAGE_HPP_

#include <cassert>
#include <fstream>
#include <iostream>
#include <sstream>
#include <string>
#include <vector>

using pixel_t = unsigned char;
using image_t = std::vector<pixel_t>;

// write image
std::string writePgm(const std::string &filename, int width, int height, 
        const image_t & data) 
{
    std::ofstream os(filename);
    if (not os)
        return "writePgm: failed to open output file";
    // write magic number
    os << "P2" << std::endl;
    // write size
    os << width << ' ' << height << std::endl;
    // write max value
    os << "255" << std::endl;
    // write data
    for (unsigned char pixel : data) 
        os << (unsigned)pixel << '\n';
    return "";
}

// read image
std::string readPgm(const std::string &filename, int & width, int & height, 
        image_t & data) 
{
    std::ifstream is(filename);
    if (not is)
        return "readPgm: failed to open input file";
    std::string tmp;
    // read magic number
    std::getline(is, tmp);
    if (tmp != "P2")
        return "readPnm: not a ASCII PGM file";
    // read comments
    std::getline(is, tmp);
    while (tmp[0] == '#')
        std::getline(is, tmp);
    // read size
    std::stringstream ss(tmp);
    ss >> width >> height;
    // read max value
    std::getline(is, tmp);
    // read data
    data = image_t(width*height);
    for (unsigned k=0; k<data.size(); k++)
    {
        int pixel;
        is >> pixel;
        data[k] = pixel;
    }
    return "";
}

// Returns a blurred subimage of data1 
// (using the specified radius for the convolution kernel).
// The size of data1 is width*height.
// The blurred region is, in data1, (x0, y0, x1, y1).
// Therefore the size of the subimage is (x1-x0)*(y1-y0)
image_t blur(const image_t & data1, int width, int height, int radius, int x0, 
        int y0, int x1, int y1)
{
    int width2 = x1-x0;
    int height2 = y1-y0;
    assert(data1.size() >= unsigned(width*height));
    assert(width2 > 0);
    assert(height2 > 0);
    image_t data2(width2*height2);
    for (int x=x0; x<x1; x++)
    {
        for (int y=y0; y<y1; y++)
        {
            unsigned s = 0;
            unsigned n = 0;
            for (int lx=-radius; lx<radius; lx++)
            {
                for (int ly=-radius; ly<radius; ly++)
                {
                    int xx = x+lx;
                    int yy = y+ly;
                    if (xx >= 0 and xx < width and yy >= 0 and yy < height)
                    {
                        s += data1[(y+ly)*width + (x+lx)];
                        n++;
                    }
                }
            }
            int pixel = n>0 ? s/n : data1[y*width + x];
            if (x==x0 or y==y0 or x==x1-1 or y==y1-1) pixel = 0; // draw border
            data2[(y-y0)*width2 + (x-x0)] = pixel;
        }
    }
    return data2;
}

image_t computeLaplacian(const image_t & data1, int width, int height, 
        double scaling)
{
    assert(data1.size() >= unsigned(width*height));
    // return the pixel value of (x,y) in [0,1]
    auto ind = [&data1,width](int x, int y) 
    {return double(data1[y*width + x])/255.0;};
    // compute laplacian image
    image_t data2(width*height);
    for (int x=1; x<width-1; x++)
    {
        for (int y=1; y<height-1; y++)
        {
            // compute laplacian value in [-4,4]
            double v = -4.0*ind(x,y) 
                + ind(x,y-1) + ind(x,y+1) + ind(x-1,y) + ind(x+1,y);
            // scale the value in [0,255]
            // and write this value in output image
            int i255 = 255.0 * (scaling*v+4.0)/8.0;
            data2[y*width + x] = std::min(255, std::max(0, i255));
        }
    }
    return data2;
}

image_t computeLaplacianOmp(const image_t & data1, int width, int height, 
        double scaling)
{
    assert(data1.size() >= unsigned(width*height));
    // return the pixel value of (x,y) in [0,1]
    auto ind = [&data1,width](int x, int y) 
    {return double(data1[y*width + x])/255.0;};
    // compute laplacian image
    image_t data2(width*height);
#pragma omp parallel for
    for (int x=1; x<width-1; x++)
    {
        for (int y=1; y<height-1; y++)
        {
            // compute laplacian value in [-4,4]
            double v = -4.0*ind(x,y) 
                + ind(x,y-1) + ind(x,y+1) + ind(x-1,y) + ind(x+1,y);
            // scale the value in [0,255]
            // and write this value in output image
            int i255 = 255.0 * (scaling*v+4.0)/8.0;
            data2[y*width + x] = std::min(255, std::max(0, i255));
        }
    }
    return data2;
}

#endif