/*
 *  psychlops_FFTW_bridge.h
 *  Psychlops Standard Library (Universal)
 *
 *  Last Modified 2010/03/05 by Kenchi HOSOKAWA
 *  (C) 2010 Kenchi HOSOKAWA, Kazushi MARUYA, Takao SATO
 */


#ifndef HEADER_PSYCHLOPS_FFTW
#define HEADER_PSYCHLOPS_FFTW

#include "../../../psychlops_core.h"
#include "fftw3.h"


namespace Psychlops {

	double log2NormalDistibution(double log_x, double octave_mu, double octave_sigma);
	double cumulativeLog2NormalDistibution(double log_x, double octave_mu, double octave_sigma);


	class FFT2 {
	protected:
		int width_, height_;
		int left, top, x_zero, y_zero;
		fftw_complex *img_spc;
		fftw_complex *frq_spc;
		void construct_default();
	public:

		// Initialize
		FFT2();
		FFT2(int width, int height);
		FFT2(const Image &source);
		FFT2(const Matrix &source);
		FFT2(const Matrix &reali, const Matrix &imagi);
		~FFT2();
		void release();
		void set(int wid, int hei);
		void set(const Image &source);
		void set(const Matrix &source);
		void set(const Matrix &reali, const Matrix &imagi);

		// Accesser to elements
		int freqX(int x) const;
		int freqY(int y) const;

		double pix(int x, int y, double l);
		double pix(int x, int y, const Color &c);
		double getPix(int x, int y);

		double getDC();
		double setDC(double l);

		// Core FFT Execution
		void fft();
		void ifft();
		void normalizeFFT();
		//Matrix normalizeMichelson(const Matrix &matrix);
		//Matrix normalizeRMS(const Matrix &matrix);

		// Core Filter
		void filtering(const Matrix &filter);
		void filtering(double cutoff_lambda1, double cutoff_lambda2, double half_power_width = 0.0625);

		// Visualizer (read only)
		void getImage(Image &absolute);
		void getSpectrum(Image &absolute, double gamma);
		void getImage(Matrix &absolute);
		void getSpectrum(Matrix &absolute, double gamma = 1.0);
		void getImage(Matrix &reali, Matrix &imagi);
		void getSpectrum(Matrix &reali, Matrix &imagi);

		// Utilities
		void makeNoise();
	protected:
		static void getSpectrumExec(FFT2 &tmp, Image &result, double gamma);
		static void getRawSpectrumExec(FFT2 &tmp, Image &reali, Image &imagi, double gamma);
	public:
		static void getSpectrum(const Image &source, Image &result, double gamma);
		static void getSpectrum(const Matrix &source, Image &result, double gamma);
		static void getSpectrum(const Image &source, Image &reali, Image &imagi, double gamma);
		static void getSpectrum(const Matrix &source, Image &reali, Image &imagi, double gamma);
	protected:
		static void filterImageExec(FFT2 &tmp, Image &result, const Matrix &carnel);
	public:
		static void filterImage(const Image &source, Image &result, const Matrix &carnel);
		static void filterImage(const Matrix &source, Image &result, const Matrix &carnel);
		static void filterImage(const Image &source, Image &result, double cutoff_lambda1, double cutoff_lambda2, double half_power_width = 0.0625);
		static void filterImage(const Matrix &source, Image &result, double cutoff_lambda1, double cutoff_lambda2, double half_power_width = 0.0625);
		static Matrix makeFilter(int height, int width, double cutoff_lambda1, double cutoff_lambda2, double half_power_width = 0.0625);
		static Matrix makeFilterSector(int height, int width, double cutoff_lambda1, double cutoff_lambda2, double theta1, double theta2, double half_power_width = 0.0625, double half_power_width_ori = PI/90.0);
		static Matrix makeFilter(int height, int width, double cutoff_lambda1, double cutoff_lambda2, double cutoff_lambda3, double cutoff_lambda4, double theta1, double theta2, double half_power_width = 0.0625);

		static int originFor(int size);
	};


	class FFTW3D {
	protected:
		int width_, height_, depth_, image_size_;
		int left, top, front, x_zero, y_zero, z_zero;
		fftw_complex *img_spc;
		fftw_complex *frq_spc;


	public:
		FFTW3D();
		FFTW3D(int wid, int hei, int dep);
		FFTW3D(const Image *source, int framsize);
		FFTW3D(const Matrix *source, int framsize);
		~FFTW3D();
		void set(int wid, int hei, int dep);
		void set(const Image *source, int framsize);
		void set(const Matrix *source, int framsize);

		int freqX(int x) const;
		int freqY(int y) const;
		int freqT(int t) const;

		double pix(int x, int y, int z, double l);
		double pix(int x, int y, int z, const Color &c);
		double getPix(int x, int y, int z);

		double getDC();
		double setDC(double l);

		void fft();
		void ifft();
		//FFT2 filter(const Matrix &filter);

		void normalizeFFT();
		//Matrix normalizeMichelson(const Matrix &matrix);
		//Matrix normalizeRMS(const Matrix &matrix);
		//void getImage(Image &target);
		//void getSpectrum(Image &target, double gamma);
		//void getImage(Matrix &real, Matrix &imag);
		//void getSpectrum(Matrix &real, Matrix &imag);
		void drawCloudImage(Drawable &target = *Drawable::prime);
		void drawCloudSpectrum(Point c, double gamma, Drawable &target = *Drawable::prime);
/*
		void filtering(const Matrix &filter);
		void filtering(double cutoff_lambda1, double cutoff_lambda2, double slope);

	protected:
		static void getSpectrumExec(FFT2 &tmp, Image &result, double gamma);
	public:
		static void getSpectrum(const Image &source, Image &result, double gamma);
		static void getSpectrum(const Matrix &source, Image &result, double gamma);
	protected:
		static void filterImageExec(FFT2 &tmp, Image &result, const Matrix &carnel);
	public:
		static void filterImage(const Image &source, Image &result, const Matrix &carnel);
		static void filterImage(const Matrix &source, Image &result, const Matrix &carnel);
		static void filterImage(const Image &source, Image &result, double cutoff_lambda1, double cutoff_lambda2, double slope);
		static void filterImage(const Matrix &source, Image &result, double cutoff_lambda1, double cutoff_lambda2, double slope);
		static Matrix makeFilter(int height, int width, double cutoff_lambda1, double cutoff_lambda2, double slope);
*/
	};

}	/*	<- namespace Psycholops 	*/


#endif
