#include "diffraction.h"
#include "winexception.h"
#include "config.h"
#include "fourier.h"
#include "ressources.h"
#include "resstring.h"
#include "parameters.h"
#include <math.h>
#include <string.h>
#include <stdio.h>
#include <stdlib.h>
#include <complex>

#define ABS(a) (arg = (a), arg > 0 ? arg : -arg)
#define ROUND(a) (arg = (a), (ceil(arg) - floor(arg) <= 0.5 ? ceil(arg) : floor(arg)) )
#define PIXEL_WHITE (-1)
#define PIXEL_BLACK (0)
#define TEST_WB(a) (arg = (a), (arg == PIXEL_WHITE || arg == PIXEL_BLACK) ? arg : 0)
#define SQR(a) (arg = (a), arg*arg)

Diffraction::Diffraction(HWND hwnd, CRaster *bmp, ColorMap *cmap):
    _hwnd(hwnd), _ndim(2), _isign(1), _bmp(bmp), _cmap(cmap),
    _four(NULL), _data(NULL), _dataAbs(NULL), _dataMin(0), _dataMax(0)
{
    _nn = new unsigned long[_ndim];
    _nn[0] = (unsigned long) (bmp->Height);
    _nn[1] = (unsigned long) (bmp->Width);

    if(_nn[0] != 0 && _nn[1] != 0)
      {
      _data = (double*) malloc(sizeof(double)*2*_nn[0]*_nn[1]);
      _dataAbs = (double*) malloc(sizeof(double)*_nn[0]*_nn[1]);
      }
}

Diffraction::Diffraction(const Diffraction &rhs):
    _hwnd(rhs._hwnd), _ndim(rhs._ndim), _isign(rhs._isign),
    _bmp(rhs._bmp), _cmap(rhs._cmap), _data(NULL), _dataAbs(NULL),
    _dataMin(rhs._dataMin), _dataMax(rhs._dataMax)
{
    _nn = new unsigned long[_ndim];
    _nn[0] = rhs._nn[0];
    _nn[1] = rhs._nn[1];
    
    if(_nn[0] != 0 && _nn[1] != 0)
      {
      _data = (double*) malloc(sizeof(double)*2*_nn[0]*_nn[1]);
      _dataAbs = (double*) malloc(sizeof(double)*_nn[0]*_nn[1]);
      
      memcpy(_data, rhs._data, sizeof(double)*2*_nn[0]*_nn[1]);
      memcpy(_dataAbs, rhs._dataAbs, sizeof(double)*_nn[0]*_nn[1]);
      }           
}

Diffraction::~Diffraction()
{
    if(_nn != NULL)      delete[] _nn;
    if(_data != NULL)    free(_data);
    if(_dataAbs != NULL) free(_dataAbs);
    if(_four != NULL)    delete[] _four;
}

Diffraction &Diffraction::operator=(const Diffraction &rhs)
{
if(&rhs != this)
  {
    _hwnd = rhs._hwnd;
    _isign = rhs._isign;
    _bmp = rhs._bmp;
    _cmap = rhs._cmap;
    _four = rhs._four;
    _dataMin = rhs._dataMin;
    _dataMax = rhs._dataMax;
    _data = NULL;
    _dataAbs = NULL;
    
    if(_ndim != rhs._ndim || _nn[0] != rhs._nn[0] || _nn[1] != rhs._nn[1])
      {
      _ndim = rhs._ndim;
      
      delete[] _nn;
      _nn = new unsigned long[_ndim];
      _nn[0] = rhs._nn[0];
      _nn[1] = rhs._nn[1];
      
      if(_nn[0] != 0 && _nn[1] != 0)
        {
        double *CopyData;
        CopyData = (double*) realloc(_data, sizeof(double)*2*_nn[0]*_nn[1]);
        
        if(!CopyData)
          throw WinException(ResString(ERROR_MEMORY_REALLOC));
        else
          _data = CopyData;
          
        CopyData = (double*) realloc(_dataAbs, sizeof(double)*_nn[0]*_nn[1]);
        
        if(!CopyData)
          throw WinException(ResString(ERROR_MEMORY_REALLOC));
        else
          _dataAbs = CopyData;  
        }
      }
  }

    return (*this);
}

// Calcul de la valeur absolue d'un vecteur de complexes
// de taille 2*N pour chacune de ses composantes.
// Retourne les valeurs absolues dans un vecteur de
// taille N.
void Diffraction::AbsoluteValue()
{
    unsigned long n, N(_nn[0]*_nn[1]);
    double x, y;
    
    for(n = 0; n < N; n++)
       {
       x = fabs(_data[2*n]);
       y = fabs(_data[2*n+1]);
  
       _dataAbs[n] = sqrt(x*x + y*y);
       }
}

// Cherche les minimums et maximum dans
// un vecteur de taille N.
void Diffraction::FindMinMax()
{
    unsigned long n, N(_nn[0]*_nn[1]);
    _dataMin = _dataMax = _dataAbs[0];
    
    for(n = 1; n < N; n++)
       {
       if(_dataAbs[n] < _dataMin)
         _dataMin = _dataAbs[n];
       else if(_dataAbs[n] > _dataMax)
         _dataMax = _dataAbs[n];
       }
}

// Lecture des données de l'image
// -1: pixel blanc
// 0: pixel noir
void Diffraction::CreateData()
{
    static double arg;
    unsigned long n, m;
    
    switch((_bmp->BPP))    
      {
      case 8: // Codage 8 bits
            for(n = 0; n < _nn[0]*_nn[1]; n++)
               {
               _data[2*n] = TEST_WB((double)(_bmp->Raster[n]))+1; // Partie réelle
               _data[2*n+1] = 0.0; // Partie imaginaire nulle
               }
      break;
      
      case 24: // Codage 24 bits 
            for(n = 0; n < _nn[0]*_nn[1]; n++)
               {
               _data[2*n] = TEST_WB((double)(_bmp->Raster[3*n]))+1; // Lit une donnée sur trois pour la partie réelle
               _data[2*n+1] = 0.0; // Partie imaginaire nulle
               }
      break;
      
      default:
            {
            char errmsg[1024] = "Le fichier contient ";
            char cBPP[8];
            sprintf(cBPP, "%d", _bmp->BPP);
            
            strcat(errmsg, cBPP);
            strcat(errmsg, " bits/pixel au lieu de 8 ou 24 bits/pixel.");
            
            throw WinException(errmsg);
            }
      }      
}

// Calcul de la solution numérique
// de l'image de diffraction.
// (Approximation de Fraunhofer)
void Diffraction::NumericalSolution(const Parameters *parameters)
{
    static const double lambda(parameters->lambda()); // [m] Longueur d'onde de l'onde plane incidente
    static const double k(2.0*(parameters->pi())/(parameters->lambda()));
    static const double z(parameters->screenDistance());  // [m] Distance à l'écran
    static const double dx(parameters->dxStep()), dy(parameters->dyStep());
    static const std::complex<double> I(0.0, 1.0); // I^2 = -1
    double X, Y, temp;
    double Xmax(_nn[0]*dx/2.0), Ymax(_nn[1]*dy/2.0);
    unsigned long n, count(0), N(_nn[0]*_nn[1]);
    std::complex<double> LinearPhase, TotalPhase;
    
    // Transformée de Fourier.
    // data (de taille 2*nn[0]*nn[1]) est arrangé ligne par ligne,
    // chaque réel étant suivit de la partie imaginaire
    // ! Requiert une longueur de data qui soit une
    // puissance de 2 !
    if(_four) _four->Fourn();
    else throw WinException(ResString(ERROR_FOURIER_INIT));
    
    // Calcul de la phase linéaire
    LinearPhase = exp(I*k*z) / (I*lambda*z);
    
    for(n = 0; n < N; n++)
       {
       // Coordonnées (X, Y) du nème noeud
       // de la grille spatiale
       X = ((double)(n % _nn[0]))*dx - Xmax;
       Y = ((double)count)*dy - Ymax;
       
       if((n % _nn[1]) == 0) count++;
       
       // Calcul de la phase totale
       // phase linéaire * phase quadratique
       TotalPhase = LinearPhase * exp(I*k/(2.0*z) * (X*X + Y*Y));
       
       // Solution numérique
       // Multiplication complexe de la transformée de
       // Fourier par la phase totale.
       temp = _data[2*n];
       _data[2*n]   = _data[2*n]*real(TotalPhase) - _data[2*n+1]*imag(TotalPhase);
       _data[2*n+1] = temp*imag(TotalPhase) + _data[2*n+1]*real(TotalPhase);
       }
}

// Shift dans les données pour un meilleur affichage
// graphique.
void Diffraction::ShiftData()
{
    unsigned long m, n;
    double *newData;
    newData = (double*) malloc(sizeof(double)*_nn[0]*_nn[1]);
    
    if(!newData)
      throw WinException(ResString(ERROR_MEMORY_ALLOC));
      
    // Décalage horizontal
    for(n = 0; n < _nn[0]; n++)
       for(m = 0; m < _nn[1]/2; m++)
          {
          newData[n*_nn[1]+m] = _dataAbs[n*_nn[1]+m+_nn[1]/2];
          newData[n*_nn[1]+m+_nn[1]/2] = _dataAbs[n*_nn[1]+m];
          }  
          
    // Copie
    memcpy(_dataAbs, newData, sizeof(double)*_nn[0]*_nn[1]);
    
    // Décalage vertical
    for(n = 0; n < _nn[0]/2; n++)
       for(m = 0;m < _nn[1]; m++)
          {
          newData[n*_nn[1]+m] = _dataAbs[(n+_nn[0]/2)*_nn[1]+m];
          newData[(n+_nn[0]/2)*_nn[1]+m] = _dataAbs[n*_nn[1]+m];
          }
          
    // Copie
    memcpy(_dataAbs, newData, sizeof(double)*_nn[0]*_nn[1]);
    
    free(newData);
}

// Remplacement des données dans 'CRaster _bmp'
void Diffraction::ReplaceData(const Parameters *parameters, bool doShift)
{
    static double arg;
    unsigned long n, m, temp, N(_nn[0]*_nn[1]);
    unsigned long cmapSize((_cmap->GetSize())/3);
    int Height(((_bmp->Height) > 0) ? (_bmp->Height) : -(_bmp->Height));
    bool HeightPos((_bmp->Height) > 0);
    int Width(_bmp->Width);
    int newBPP(24);
    int newBytesPerRow;
    double mcmap, pcmap;
    
    newBytesPerRow = (Width*newBPP)/8;
    newBytesPerRow += (4 - newBytesPerRow%4) % 4;
    
    // Recherche des valeurs minimales et maximales
    FindMinMax();

    // Shift pour un meilleur affichage graphique
    if(doShift) ShiftData();
    
    // Calcul de la correspondance linéaire
    // entre les valeurs des données de 'dataAbs'
    // et la carte de couleurs 'cmap'
    if(fabs(_dataMax - _dataMin) <= (parameters->tiny()))
      {
      mcmap = 0.0;
      pcmap = (_dataMin < 0) ? _dataMin : -_dataMin;
      }
    else
      {
      mcmap = (-(double)(cmapSize-1)) / (_dataMin - _dataMax);
      pcmap = ((double)cmapSize*_dataMin) / (_dataMin - _dataMax);
      }   
      
    // Mise à jour des variables de CRaster *bmp
    if((_bmp->BPP) != newBPP)
      {
      _bmp->BPP = newBPP;
      _bmp->BytesPerRow = newBytesPerRow;
      
      BITMAPINFOHEADER *bmih;
      bmih = &(_bmp->pbmi->bmiHeader);
      bmih->biBitCount = newBPP;
      bmih->biSizeImage = newBytesPerRow;
      
      // Mise à jour de Raster
      char *newRaster;
      newRaster = (char*) realloc(_bmp->Raster, sizeof(char)*newBytesPerRow*Height);
      
      if(!newRaster)
        throw WinException(ResString(ERROR_MEMORY_REALLOC));
      else
        (_bmp->Raster) = newRaster;
      }
      
      // Remplissage de l'image par les données
      for(n = 0; n < Height; n++)
         {
         char *pPixel;
         pPixel = (char*)(_bmp->Raster) + newBytesPerRow *
                (HeightPos ? n : Height-n-1);
   
         for(m = 0; m < Width; m++)
            {
            temp = 3*(unsigned long)floor(mcmap*_dataAbs[n*Width+m]+pcmap);
                
            *(pPixel+3*m)   = (char) ROUND((_cmap->map[temp+2])*255.0);
            *(pPixel+3*m+1) = (char) ROUND((_cmap->map[temp+1])*255.0);
            *(pPixel+3*m+2) = (char) ROUND((_cmap->map[temp])*255.0);
            }
         }
}

// Recharge la carte des couleurs
void Diffraction::ReloadMap(const Parameters *parameters)
{
    if(_dataAbs != NULL) ReplaceData(parameters, false);
}

// Recharge une nouvelle carte des couleurs
void Diffraction::ReloadMap(const Parameters *parameters, ColorMap *newMap)
{
    _cmap = newMap;
    if(_dataAbs != NULL) ReplaceData(parameters, false);
}

// Calcul de la figure de diffraction
// à partir d'un diaphragme fournit
// par un fichier image bitmap (8 ou 24 bits uniquement).
// (Approximation de Fraunhofer)
//
// !! Les dimensions requises de l'image
// sont celle d'un carré dont la longueur est une puissance de 2.
// Ceci est requis pas l'algorithme calculant
// la transformée de Fourier (FFT).
void Diffraction::DoDiffraction(const Parameters *parameters)
{
    if(_bmp != NULL)
      {
      // Lecture des données du diaphragme
      CreateData();
    
      // Initialisation de la transformée de Fourier
      _four = new Fourier(_data, _nn, _ndim, _isign);
    
      // Calcul de la solution numérique
      // de l'image de diffraction
      NumericalSolution(parameters);
    
      // Data contient à-priori des composantes complexes non-nulles
      // Calcul de la valeur absolue
      AbsoluteValue();
    
      // Remplacement des données dans 'CRaster _bmp'
      ReplaceData(parameters, true);
      }
}
