#define cimg_display 1
#include <iostream>
#include <sstream>
#include <string>
#include <math.h>
#include <gsl/gsl_linalg.h>
#include <time.h>
#include "CImg.h"

using namespace cimg_library;
using namespace std;

unsigned char w( unsigned char z ) {
  if (z < 128){
    return z;
  }else {
    return 255-z;
  }
}


int main() {
  //[TODO ]IO reading a file list
  // ugly placeholder

  CImg<unsigned char> imageList[12];

  for (int i = 0; i<12 ; i++){
    string filename;
    stringstream index;
    if (i<9){
      index << i+1;
      filename = "../pics/candy/Img0" + index.str() + ".jpg";
    }else {
      index << i+1;
      filename = "../pics/candy/Img" + index.str() + ".jpg";
    }
    imageList[i] = CImg<unsigned char>(filename.c_str());
  }

  int sample_index_x[25]; // N samples
  int sample_index_y[25];
  srand( time( NULL ) );
  
  //hardcoded shutter speed
  double shutter_speed[12]={13.0, 6.0, 3.0,  1.6, 
                            0.77, 0.5, 0.25, 0.125, 
                            0.0625, 0.03125, 0.015625,
                            0.0078125
                             };
//  double speed = 13.0;

  const int P = 12;
  const int N = 25;
  const int n = 256; 

  // get width and height for sampling
  // cout << "spec: " << imageList[0].spectrum() << endl;
  // cout << "height: " << imageList[0].width() << endl;
  for( int i = 0 ; i < N ; i++ ){
    sample_index_x[i] = rand() % imageList[0].width();
    sample_index_y[i] = rand() % imageList[0].height();
    printf( "x : %d\ty : %d\n" , sample_index_x[i] , sample_index_y[i] );
  }
  // shutter speeds
  for( int i = 0 ; i < P ; i++ ){
    shutter_speed[i] = log( shutter_speed[i] );
    //speed /= 2;
  }

  //  N   P
  //(25 * 12) + 254 + 1 = 555
  // 25 + 256 = 281

  double _A_R [555*281]={0.0};
  double _A_G [555*281]={0.0};
  double _A_B [555*281]={0.0};
  double _b_R [555]={0.0};
  double _b_G [555]={0.0};
  double _b_B [555]={0.0};
//  double _x [281];

  gsl_matrix      *V_matrix;
  gsl_vector      *S_vector, *work;

  V_matrix = gsl_matrix_alloc( 256 + N , 256 + N );
  S_vector = gsl_vector_alloc( 256 + N );
  work = gsl_vector_calloc( 256 + N );

  //solving linear algebra
  int k = 0;
  for (int i =0; i< 25;i++){			// for all sample
    for (int j = 0; j< 12;j++){ 		// for all pics
      unsigned char R_val, G_val, B_val;     
      R_val = imageList[j](sample_index_x[i], sample_index_y[i], 0);
      G_val = imageList[j](sample_index_x[i], sample_index_y[i], 1);
      B_val = imageList[j](sample_index_x[i], sample_index_y[i], 2);

      unsigned char wij;
      // R
      wij = w(R_val);
      _A_R[k*281 + R_val] = wij; 
      _A_R[k*281 + n + i] = -wij; 
      _b_R[k] = wij*shutter_speed[j]; 

      // G
      wij = w(G_val);
      _A_G[k*281 + G_val] = wij; 
      _A_G[k*281 + n + i] = -wij; 
      _b_G[k] = wij*shutter_speed[j]; 

      // B
      wij = w(B_val);
      _A_B[k*281 + B_val] = wij; 
      _A_B[k*281 + n + i] = -wij; 
      _b_B[k] = wij*shutter_speed[j]; 

      // increase k
      k++;
    }
  }

  _A_R[k*281+128] = 1.0;
  _A_G[k*281+128] = 1.0;
  _A_B[k*281+128] = 1.0;
  k++;
  
  // smooth equations
  for (unsigned char i =0; i< 254; i++){
    //R
    _A_R[k*281+i] = 1*w(i+1);
    _A_R[k*281+i+1] = -2*w(i+1);
    _A_R[k*281+i+2] = 1*w(i+1);

    //G
    _A_G[k*281+i] = 1*w(i+1);
    _A_G[k*281+i+1] = -2*w(i+1);
    _A_G[k*281+i+2] = 1*w(i+1);

    //B
    _A_B[k*281+i] = 1*w(i+1);
    _A_B[k*281+i+1] = -2*w(i+1);
    _A_B[k*281+i+2] = 1*w(i+1);

    k++;
  }

  //solve equation and output file
  //-----------------------------

  FILE	*R_out, *G_out, *B_out;
  R_out = fopen("R_curve.txt","w");
  G_out = fopen("G_curve.txt","w");
  B_out = fopen("B_curve.txt","w");


  gsl_matrix_view A_view;
  gsl_vector_view b_view;

  //R
  A_view = gsl_matrix_view_array( _A_R , N * P + 1 + 254 , 256 + N );
  b_view = gsl_vector_view_array( _b_R , N * P + 1 + 254 );
  gsl_linalg_SV_decomp( &A_view.matrix , V_matrix , S_vector , work );
  gsl_vector_set_all( work , 0.0 );
  gsl_linalg_SV_solve( &A_view.matrix , V_matrix , S_vector , &b_view.vector , work );
/*
  for (int i=0;i < 555*281; i++){
    cout << _A_R[i] << " ";
    if (i% 281 ==0)
      cout << endl;
  }
*/
//  for (int i =0; i<n;i++){
//    cout << gsl_vector_ptr(work,0)[i] << " " ;
//  }
  for( int i = 0 ; i < n ; i++ ){
    fprintf( R_out , "%lf\n" , gsl_vector_get( work , i ));
  }  
  gsl_vector_set_all( work , 0.0 );
  gsl_vector_set_all( S_vector , 0.0 );
  gsl_matrix_set_all( V_matrix , 0.0 );

  //G
  A_view = gsl_matrix_view_array( _A_G , N * P + 1 + 254 , 256 + N );
  b_view = gsl_vector_view_array( _b_G , N * P + 1 + 254 );
  gsl_linalg_SV_decomp( &A_view.matrix , V_matrix , S_vector , work );
  gsl_vector_set_all( work , 0.0 );
  gsl_linalg_SV_solve( &A_view.matrix , V_matrix , S_vector , &b_view.vector , work );
  for( int i = 0 ; i < n ; i++ ){
    fprintf( G_out , "%lf\n" ,  gsl_vector_get( work , i ));
  }  
  gsl_vector_set_all( work , 0.0 );
  gsl_vector_set_all( S_vector , 0.0 );
  gsl_matrix_set_all( V_matrix , 0.0 );

  //B
  A_view = gsl_matrix_view_array( _A_B , N * P + 1 + 254 , 256 + N );
  b_view = gsl_vector_view_array( _b_B , N * P + 1 + 254 );
  gsl_linalg_SV_decomp( &A_view.matrix , V_matrix , S_vector , work );
  gsl_vector_set_all( work , 0.0 );
  gsl_linalg_SV_solve( &A_view.matrix , V_matrix , S_vector , &b_view.vector , work );
  for( int i = 0 ; i < n ; i++ ){
    fprintf( B_out , "%lf\n" , gsl_vector_get( work , i ));
  }  



  fclose(R_out);
  fclose(G_out);
  fclose(B_out);

  return 0;
}