//---------------------------------------------------------------------------
#define iloscWezlow 1000
#include <vcl.h>
#include <math.h>
#include "Matrix.h"
#include "nr.h"
#include "Polaczone.h"

#pragma hdrstop

//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
//BITMAPINFO bmi;
//unsigned char *bits;
Graphics::TBitmap *bmp;
Graphics::TBitmap* Dib;
matrix* Macierz;
int DispMode=2;
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner)
        : TForm(Owner)
{

}
//---------------------------------------------------------------------------

void __fastcall TForm1::WczytajClick(TObject *Sender)
{
if(OpenDialog1->Execute())
  {
   bmp = Image1->Picture->Bitmap;
   bmp->LoadFromFile(OpenDialog1->FileName);
 //  bmp->PixelFormat=pf8bit;    //!!!!!!!!!!!!!!!!
 //  bmp->Palette = GrayPalette();
   StatusBar1->SimpleText = OpenDialog1->FileName;
   Memo1->Lines->Append(OpenDialog1->FileName) ;// ();
   };
}
//---------------------------------------------------------------------------

void __fastcall TForm1::ZamknijClick(TObject *Sender)
{
 if(SaveDialog1->Execute()) Memo1->Lines->SaveToFile(SaveDialog1->FileName);

 Close();
}

//---------------------------------------------------------------------------
void __fastcall TForm1::PoliczFFTClick(TObject *Sender)
{
 int isig=1;

// matrix* Macierz;

 matrix* Macierz_Wizualizacja;

 Macierz=MatrixFromDib(bmp);

 Fft2(Macierz->x(), Macierz->n()/2, Macierz->m(), isig);

 Macierz_Wizualizacja=Macierz;

 SortMatrix(Macierz_Wizualizacja);

 try {
        FillDibWithMatrix(Image2->Picture->Bitmap,DispMode,100,Macierz_Wizualizacja); }
 catch(...){
        ShowMessage ("Nie pracuje");};

}

//-------------------------------------------------------------------------

matrix* __fastcall TForm1::PowerSpectrum(matrix* Matrix1)
{
 double *y;
 double real, imag;
 int indeks;

 int nCol = Matrix1->n()/2;
 int nRow = Matrix1->m();
 int N = nCol*nRow;

 y=new double[N];
 memset(y, 0,N*sizeof(double));

 indeks=0;
 double *tmp=Matrix1->x();
 for(int j=0;j<nRow;j++){
        for(int i=0;i<nCol;i++){
                real=*tmp++;    //real=*tmp  *tmp=*tmp+1
                imag=*tmp++;
                //Memo1->Lines->Add(FloatToStr(real));
                y[indeks]=(real)*(real)+(imag)*(imag);
                indeks++;
                }
        }

//for(int i=0;i<N;i++ )
//Memo1->Lines->Add(FloatToStr(y[i])); //4????Nie jestem pewna,czy jest ok

 matrix* spectrum=new matrix(nRow,nCol,y);
 delete[] y;
 return spectrum;
}

//---------------------------------------------------------------------------

double __fastcall TForm1::Maksimum (matrix* Matrix2)
{
 double maks = Matrix2->max();
  //Memo2->Lines->Add(FloatToStr(maks));
 return maks;
}


//---------------------------------------------------------------------------

double __fastcall TForm1::Poziomica (double s, double data[], int nx, int ny)
{
 double W = 0.0;

   for(int i=0;i<nx*ny;i++ )
          W+=((data[i])> s ? 1.0 :0.0);
//    Memo2->Lines->Add(FloatToStr(W));
 return W;
 }

//---------------------------------------------------------------------------

double __fastcall TForm1::ZnajdzBeta1 (matrix* Matrix, int mx, int my)
{
 double u[iloscWezlow];
 double v[iloscWezlow];
 double max = Maksimum(Matrix);
 double us=0, uks=0, vs=0, uvs=0;

 try {
      for(int k=0; k<iloscWezlow; k++){
        u[k] = log(max)-0.35-2.0*k/(1.0*iloscWezlow);
        v[k] = log(Poziomica(exp(u[k]), Matrix->x(), mx,my));
        };;}
 catch(...){
        ShowMessage ("PALANT!!!");};

 for(int k=0; k<iloscWezlow; k++){
        us+=u[k]/(1.0*iloscWezlow);
        uks+=u[k]*u[k]/(1.0*iloscWezlow);
        vs+=v[k]/(1.0*iloscWezlow);
        uvs+=u[k]*v[k]/(1.0*iloscWezlow);
        };
 double aa= (uvs-us*vs)/(uks-us*us);
 return aa;
 }
//---------------------------------------------------------------------------

double __fastcall TForm1::ZnajdzBeta2 (matrix* Matrix, int mx, int my)
{
 double u[iloscWezlow];
 double v[iloscWezlow];
 double max = Maksimum(Matrix);
 double us=0, uks=0, vs=0, uvs=0;

 try {
      for(int k=1; k<iloscWezlow-1; k++){
        u[k] = log(k*max/(1.0*iloscWezlow));
        v[k] = log( Poziomica ( exp(u[k]), Matrix->x(), mx, my ) );
        };
      }
 catch(...){
        ShowMessage ("PALANT!!!");};

 for(int k=1; k<iloscWezlow-1; k++){
        us+=u[k]/(1.0*iloscWezlow-2.0);
        uks+=u[k]*u[k]/(1.0*iloscWezlow-2.0);
        vs+=v[k]/(1.0*iloscWezlow-2.0);
        uvs+=u[k]*v[k]/(1.0*iloscWezlow-2.0);
        };
 double aa= (uvs-us*vs)/(uks-us*us);
 return aa;
 }


//---------------------------------------------------------------------------

matrix* __fastcall TForm1::MatrixFromDib(Graphics::TBitmap *bmp)
{
   bmp->PixelFormat=pf8bit;    //!!!!!!!!!!!!!!!!
   bmp->Palette = GrayPalette();
   unsigned int *data,*tmp;
   int W=bmp->Width;
   int H=bmp->Height;
   int newH,newW;
   newH=newW=1;
	while(H>newH) newH*=2;
 	while(W>newW) newW*=2;

   data = new unsigned int [newW*newH*2];
   memset(data, 0,newW*newH*2*sizeof(unsigned int));
	tmp=data;
   unsigned char* ptr;

   for (int j = 0; j<H ; j++){
        ptr = (unsigned char *)bmp->ScanLine[j];
   	for (int i = 0; i < W ; i++){
               	 *tmp++ = ptr[i];
                  tmp++ ;
                //Memo1->Lines->Add(FloatToStr(ptr[i])); //1OK
        }
        for(int k=W;k<newW;k++){
             	tmp++;
                tmp++;
        }

   }
   //for(int i=0;i<newH*newW*2;i++)
      // Memo1->Lines->Add(FloatToStr(data[i]));   //2OK
   matrix* m=new matrix(newH,newW*2,data);
   delete[] data;
   return m;
}

//---------------------------------------------------------------------------


void __fastcall TForm1::Fft2(double data[], int nx, int ny, int isign){
unsigned long nn[3];
nn[0]=0;                    ////????????????
nn[1]=(unsigned long) nx;
nn[2]=(unsigned long) ny;
fourn (data-1, nn, 2, isign);
if (isign==1){
        double c= double(1.0)/double(nn[1]*nn[2]);             //NORMALIZACJA!!!
        for( unsigned long i=0;i<2*nn[1]*nn[2];i++ ) data[i]*=c;
      };

//for (unsigned long i=0;i<2*nn[1]*nn[2];i++)
//Memo1->Lines->Add(FloatToStr(data[i]));   //3 OK KONTROLNE do data[]  Znak?
};


//---------------------------------------------------------------------------

void __fastcall TForm1::FillDibWithMatrix(Graphics::TBitmap* Dib,int DispMode,double gain,matrix* M){

if(!M  || !Dib ) return;

   int i,j;
   int W=M->n()/2;
   int H=M->m();
   double r,im,max,tmp;
   Dib->HandleType = bmDIB;
   Dib->PixelFormat = pf8bit;
   Dib->Palette = GrayPalette();
   Dib->Height = H;
   Dib->Width = W;
  // Image2->Picture->Bitmap = Dib;         ///?????to dodalam
   double *v = M->x();
   unsigned char *ptr;

   switch(DispMode){
      	case /*DIB_DISP_REAL*/ 0:
                         for (j = 0; j < H; j++){
                         ptr=(unsigned char *)Dib->ScanLine[j];
                         for (i = 0; i < W; i++){
                         ptr[i] = (unsigned char) floor( gain*fabs( *v++)) ;
                                v++;
      	                }
                  }
                break;
      case /*DIB_DISP_IMAG*/ 1:
                 v++;
 		for (j = 0; j < H; j++){
                         ptr=(unsigned char *)Dib->ScanLine[j];
                         for (i = 0; i < W; i++){
                         ptr[i] = (unsigned char) floor( gain*fabs( *v++ )) ;
                                v++;
                        }
                }
                break;

      case /*DIB_DISP_MODULE*/ 2:
		for (j = 0; j < H; j++){
                ptr=(unsigned char *)Dib->ScanLine[j];
                for (i = 0; i < W; i++){
                r = *v++;
                im = *v++;
                tmp = floor( gain*sqrt(r*r+im*im));
                                if(tmp>255)
				  ptr[i] = 255;
                                else
                                  ptr[i] = tmp;
              	      }
		}
            break;
      case /*DIB_DISP_POWER*/ 3:
       	    for (j = 0; j < H ; j++){
            ptr=(unsigned char *)Dib->ScanLine[j];
             	    	for (i = 0; i < W; i++){
                        r = *v++;
   	          	im = *v++;
                        tmp = floor( gain*(r*r+im*im) );
                                if(tmp>255)
				  ptr[i] = 255;
                                else
                                  ptr[i] = (unsigned char)tmp;

      	      }
	    }
            break;
            }
}

//-------------------------------------------------------------------------

HPALETTE __fastcall TForm1::GrayPalette()
 {
 HPALETTE hpal;
 WORD nColors=256;
 LOGPALETTE* logPal = (LOGPALETTE*) new unsigned char[sizeof(LOGPALETTE)+(nColors-1)*sizeof(PALETTEENTRY)];
 logPal->palVersion = 0x300;
 logPal->palNumEntries = nColors;
 for(int i=0; i<nColors;i++){
 logPal->palPalEntry[i].peRed   = i;
 logPal->palPalEntry[i].peGreen = i;
 logPal->palPalEntry[i].peBlue  = i;
 logPal->palPalEntry[i].peFlags = 0; }
 hpal = CreatePalette(logPal);
 delete[] logPal;
 return hpal;
 }

//------------------------------------------------------------------------

 void __fastcall TForm1::SortMatrix(matrix* Matrix) {
    int j;    int H=Matrix->m();
   int W=Matrix->n();
   double *p1,*p2,*tmp;

   // Najpierw przestawiamy kolumny
   tmp = new double[W/2];
   for (j = 0; j < H ; j++){
       p1 = Matrix->x() + j*W;
       p2 = p1+W/2 ;
        memcpy(tmp,p1,sizeof(double)*W/2);
        memcpy(p1,p2,sizeof(double)*W/2);
        memcpy(p2,tmp,sizeof(double)*W/2);
  }
   delete[] tmp;

  //Przestawiamy wiersze
   tmp = new double[W*H/2];
   p1= Matrix->x();
   p2= p1+W*H/2;
   memcpy(tmp,p1,sizeof(double)*W*H/2);
   memcpy(p1,p2,sizeof(double)*W*H/2);
   memcpy(p2,tmp,sizeof(double)*W*H/2);
   delete[] tmp;

}

//-------------------------------------------------------------------------
void __fastcall TForm1::Metoda1Click(TObject *Sender)  //Spektralny wymiar fraktalny metod¹ 1
{
 matrix* MacierzB1;

 matrix* PowSpec;

 MacierzB1=Macierz;

 PowSpec=PowerSpectrum(MacierzB1);

 double a1=ZnajdzBeta1 (PowSpec, PowSpec->n(), PowSpec->m());

 Edit1->Text=-2.0/a1;

 Memo1->Lines->Add(FloatToStr(-2.0/a1));
}
//---------------------------------------------------------------------------

void __fastcall TForm1::Metoda2Click(TObject *Sender)  //Spektralny wymiar fraktalny metod¹ 2
{
 matrix* MacierzB2;

 matrix* PowSpec;

 MacierzB2=Macierz;

 PowSpec=PowerSpectrum(MacierzB2);

 double a2=ZnajdzBeta2 (PowSpec, PowSpec->n(), PowSpec->m());

 Edit2->Text=-2.0/a2;

  Memo1->Lines->Add(FloatToStr(-2.0/a2));
}
//---------------------------------------------------------------------------