//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop

#include "imagepro.h"
#include <values.h>
#include <math.h>
#include <time.h>

//Copies TImage object data to unsigned char *
void imagepro::CopyImageToPtr(TImage *Image, unsigned char *Ptr) {
  /*unsigned char *ptr;
  for(int y = 0; y < H; y++){
    ptr = (unsigned char*)Image->Picture->Bitmap->ScanLine[y];
    memcpy(Ptr,ptr,W);
    Ptr += W;
  } */
    unsigned char *PtrRed = Ptr;
    unsigned char *PtrGreen = Ptr + W*H;
    unsigned char *PtrBlue = Ptr + 2*W*H;
    RGBTRIPLE *ptrc;
    for(int y=0, x; y<H; y++){
       ptrc = (RGBTRIPLE *)Image->Picture->Bitmap->ScanLine[y];
       for(x=0; x<W; x++){
           *PtrRed = (unsigned char)(*ptrc).rgbtRed;
           *PtrBlue = (unsigned char)(*ptrc).rgbtBlue;
           *PtrGreen = (unsigned char)(*ptrc).rgbtGreen;
           PtrRed++;
           PtrBlue++;
           PtrGreen++;
           ptrc++;
       }
    }
}

//Copies data of unsigned char * to TImage object
void imagepro::CopyPtrToImage(unsigned char *Ptr, TImage *Image) {
  /*unsigned char *ptr;
  for(int y = 0; y < H; y++){
    ptr = (unsigned char*)Image->Picture->Bitmap->ScanLine[y];
    memcpy(ptr,Ptr,W);
    Ptr += W;
  } */

    RGBTRIPLE *ptrc;
    RGBTRIPLE *x1 = new RGBTRIPLE;
    for(int y = 0; y < H; y++){
        ptrc = (RGBTRIPLE *)Image->Picture->Bitmap->ScanLine[y];
        for(int x = 0; x < W; x++){
            //if(GrayToColor != NULL)
            //    *ptrc = (RGBTRIPLE)GrayToColor[*Ptr];
            //else{
                x1->rgbtBlue = *Ptr;
                x1->rgbtGreen = *Ptr;
                x1->rgbtRed = *Ptr;
                *ptrc = *x1;
            //}
            Ptr++;
            ptrc++;
        }
    }
    delete x1;
}

//---------------------------------------------------------------------------
imagepro::imagepro() {
  nOfPts = 0;
  nOfIterations = 0;
  fcounter = 0;
}

imagepro::~imagepro() {

}
//---------------------------------------------------------------------------
//Loads a BMP image from disk
void imagepro::LoadBMPImage(TImage *currImage, AnsiString fname) {
  Graphics::TBitmap *timg;
  timg = new Graphics::TBitmap;
  timg->LoadFromFile(fname);
  timg->PixelFormat = pf24bit;
  currImage->Width = timg->Width;
  currImage->Height = timg->Height;
  currImage->Picture->Bitmap = timg;
  H = timg->Height;
  W = timg->Width;
  img = new unsigned char[W*H*3];
  wrk = new unsigned char[W*H*3];

  CopyImageToPtr(currImage, img);
  CopyImageToPtr(currImage, wrk);

  delete timg;
}
//---------------------------------------------------------------------------
//
void imagepro::SobelOperator(TImage *dImage) {
  int sum, minv, maxv;

  int *M1;
  M1 = new int[W*H];
  int *M2;
  M2 = new int[W*H];
  S = new int[W*H];

  for(int i = 0; i < (W*H); i++) {
    M1[i] = 0;
    M2[i] = 0;
  }

  for (int i = 1; i < W-1; i++)
    for (int j = 1; j < H-1; j++) {
      sum = 0;
      sum += wrk[(j-1)*W+(i-1)];
      sum += wrk[(j-1)*W+i]*2;
      sum += wrk[(j-1)*W+(i+1)];
      sum += -wrk[(j+1)*W+(i-1)];
      sum += -wrk[(j+1)*W+i]*2;
      sum += -wrk[(j+1)*W+(i+1)];
      M1[j*W+i] = (int)((float)sum/4);
    }

  for (int i = 1; i < W-1; i++)
    for (int j = 1; j < H-1; j++) {
      sum = 0;
      sum += -wrk[(j-1)*W+(i-1)];
      sum += wrk[(j-1)*W+(i+1)];
      sum += -wrk[j*W+(i-1)]*2;
      sum += wrk[(j)*W+(i+1)]*2;
      sum += -wrk[(j+1)*W+(i-1)];
      sum += wrk[(j+1)*W+(i+1)];
      M2[j*W+i] = (int)((float)sum/4);
    }

  minv = 255;
  maxv = 0;
  for (int i = 0; i < W; i++)
    for (int j = 0; j < H; j++) {
      S[j*W+i] = (int)sqrt(M1[j*W+i]*M1[j*W+i] + M2[j*W+i]*M2[j*W+i]);
      if (S[j*W+i] < minv) minv = S[j*W+i];
      if (S[j*W+i] > maxv) maxv = S[j*W+i];
    }

  dS = new unsigned char[W*H];
  for (int i = 0; i < W; i++)
    for (int j = 0; j < H; j++) {
        dS[j*W+i] = (unsigned char)((float)(S[j*W+i]-minv)/(maxv-minv)*255);
        if (dS[j*W+i] > 5)
          dS[j*W+i] = 255;
        else
          dS[j*W+i] = 0;
    }

  wrk = dS;
  CopyPtrToImage(wrk, dImage);
  dImage->Refresh();
  delete M1, M2;
}

//---------------------------------------------------------------------------
//blurring the image
void imagepro::BlurImage(TImage *dImage)
{
  B = new unsigned char[W*H];
  int sum;
  int f_width = 5;
  float F[5] = {1, 3, 9, 3, 1};


  //filter each row
  for (int j = 0; j < H; j++) {
    for (int i = 0; i < W; i++) {
      if (!((i < f_width/2)||(i > W-f_width/2))) {
        sum = 0;
        for (int k = 0; k < f_width; k++)
          sum += wrk[j*W+(i-f_width/2+k)]*F[k];
        B[j*W+i] = (int)((float)sum/17);
      }
    }
  }
  //filter each column
  for (int i = 0; i < W; i++) {
    for (int j = 0; j < H; j++) {
      if (!((j < f_width/2)||(j > H-f_width/2))) {
        sum = 0;
        for (int k = 0; k < f_width; k++)
          sum += B[(j-f_width/2+k)*W+i]*F[k];
        B[j*W+i] = (int)((float)sum/17);
      }
    }
  }

  wrk = B;
  CopyPtrToImage(wrk, dImage);
  dImage->Refresh();
}
//---------------------------------------------------------------------------
int imagepro::ReturnIntensity(int x, int y)
{
  if ((x >= 0)&&(x < W)&&(y >= 0)&&(y < H))
    return dS[y*W + x];
  else
    return 0;
}
//---------------------------------------------------------------------------
int imagepro::ReturnImageIntensity(int x, int y)
{
  return img[y*W + x];
}
//---------------------------------------------------------------------------
void imagepro::FirstPointToList(int x, int y, TListBox *lbox)
{
  nOfPts = 1;
  xL = new int[nOfPts];
  yL = new int[nOfPts];
  xL[0] = x;
  yL[0] = y;
  lbox->Items->Add("x = "+(AnsiString)x+"   y = "+(AnsiString)y);
}
//---------------------------------------------------------------------------
void imagepro::AddPointToList(int x, int y, TListBox *lbox)
{
  nOfPts++;
  int *xT;
  int *yT;
  xT = new int[nOfPts];
  yT = new int[nOfPts];
  for (int i = 0; i < nOfPts-1; i++) {
    xT[i] = xL[i];
    yT[i] = yL[i];
  }
  delete xL;
  delete yL;

  xT[nOfPts-1] = x;
  yT[nOfPts-1] = y;
  xL = new int[nOfPts];
  yL = new int[nOfPts];
  for (int i = 0; i < nOfPts; i++) {
    xL[i] = xT[i];
    yL[i] = yT[i];
  }
  delete xT;
  delete yT;
  lbox->Items->Add("x = "+(AnsiString)x+"   y = "+(AnsiString)y);
}
//---------------------------------------------------------------------------
bool imagepro::AddPointInBetweenList(int pos)
{
  int newXpos;
  int newYpos;
  if (pos == nOfPts-1) {
    newXpos = (xL[pos] + xL[0])/2;
    newYpos = (yL[pos] + yL[0])/2;
  } else {
    newXpos = (xL[pos] + xL[pos+1])/2;
    newYpos = (yL[pos] + yL[pos+1])/2;
  }
  bool noadd = true;

  if ((newXpos >= 0)&&(newXpos < W)&&(newYpos >= 0)&&(newYpos < H)) {
    int *xT;
    xT = new int[nOfPts+1];
    int *yT;
    yT = new int[nOfPts+1];

    for (int i = 0; i <= pos; i++) {
      xT[i] = xL[i];
      yT[i] = yL[i];
    }
    xT[pos+1] = newXpos;
    yT[pos+1] = newYpos;
    for (int i = pos+1; i < nOfPts; i++) {
      xT[i+1] = xL[i];
      yT[i+1] = yL[i];
    }
    delete xL, yL;
    nOfPts++;
    xL = new int[nOfPts];
    yL = new int[nOfPts];
    for (int i = 0; i < nOfPts; i++) {
      xL[i] = xT[i];
      yL[i] = yT[i];
    }
    delete xT, yT;
    noadd = false;
  }
  //returns true if point is NOT added in the list
  return noadd;
}
//---------------------------------------------------------------------------
void imagepro::DrawX(int x, int y, int s, TImage *dImage)
{
  TCanvas *can;
  can = dImage->Picture->Bitmap->Canvas;
  can->Brush->Color = clRed;
  can->Pen->Color = clRed;
  can->Pen->Width = 2;
  can->MoveTo(x-s,y-s);
  can->LineTo(x+s,y+s);
  can->MoveTo(x+s,y-s);
  can->LineTo(x-s,y+s);
}
//---------------------------------------------------------------------------
void imagepro::DrawLine(int x1, int x2, int y1, int y2, TImage *dImage)
{
  TCanvas *can;
  can = dImage->Picture->Bitmap->Canvas;
  can->Brush->Color = clRed;
  can->Pen->Color = clRed;
  can->MoveTo(x1,y1);
  can->LineTo(x2,y2);
}
//---------------------------------------------------------------------------
void imagepro::RemoveLastPixelFromList(TListBox *lbox, TImage *dImage)
{
  nOfPts--;
  lbox->Items->Delete(nOfPts);
  CopyPtrToImage(wrk, dImage);
  dImage->Refresh();
  for (int i = 0; i < nOfPts; i++)
    DrawX(xL[i],yL[i],3,dImage);
}
//---------------------------------------------------------------------------
void imagepro::CopyList(int *source, int *target, int n)
{
  for (int i = 0; i < n; i++)
    target[i] = source[i];
}
//---------------------------------------------------------------------------
void imagepro::SaveBitmapImage(TImage *Image, AnsiString fname)
{
  if (fcounter < 100) fname = fname + "0";
  if (fcounter < 10) fname = fname + "0";
  fname = fname + (AnsiString)fcounter + ".bmp";
  Image->Picture->Bitmap->SaveToFile(fname);
  fcounter++;
  Application->ProcessMessages();
}
//---------------------------------------------------------------------------
void imagepro::Algorithm(TImage *dImage, double alpha, double beta, double gamma, int nbhood, int distance)
{
  bool noNewPoints = false;
  int AlgCount = 0;

  while ((!noNewPoints)&&(AlgCount < 5)) {
    //compute snake optimisation
    SnakeIterations(dImage, alpha, beta, gamma, nbhood);
    /*//refresh image display
      CopyPtrToImage(wrk, dImage);
      dImage->Refresh();
      for (int i = 0; i < nOfPts; i++)
        DrawX(xL[i],yL[i],3,dImage);
    //save currrent display to file
      SaveBitmapImage(dImage);*/
    Application->ProcessMessages();
    //calculate new points if needed
    noNewPoints = CalculateNewListPoints(distance);
    //noNewPoints = true;
    AlgCount++;
  }

}
//---------------------------------------------------------------------------
void imagepro::SnakeIterations(TImage *dispImage, double alpha, double beta, double gamma, int nbhood)
{
  int StopSnake = 0.8*nOfPts; //number of the points that should be static
  int nOfStaticPoints;  //number of points that are static
  bool StopSnakeFlag = false; //stop the snake algorithm
  bool StopSnakeFlagReal = false; //stop the snake algorithm
  AnsiString FName = "Frame";

  //temporary (iteration) point list
  int *xLt;
  xLt = new int[nOfPts];
  CopyList(xL, xLt, nOfPts);
  int *yLt;
  yLt = new int[nOfPts];
  CopyList(yL, yLt, nOfPts);

  //temporary (comparison) point list
  int *xLtt;
  xLtt = new int[nOfPts];
  int *yLtt;
  yLtt = new int[nOfPts];

  double ave; //the average distance between points
  double energySum;  //the sum of energies

  //calculate the average distance between points in the list
  ave = aveDCount(xL, yL, nOfPts);
  //list of Enegry values
  double *E;
  E = new double[nOfPts];
  for (int i = 0; i < nOfPts; i++)
    E[i] = (alpha*energyCont(ave,i,xL,yL,0,0,nOfPts) + beta*energyCurv(i,xL,yL,0,0,nOfPts) + gamma*energyIm(i,xL,yL,0,0));

  //MAIN SNAKE ITERATIONS
  while ((!StopSnakeFlagReal)&&(nOfIterations < 50)) {

    if (StopSnakeFlag)
      StopSnakeFlagReal = true;

    //make the temporary (comparison) list
    CopyList(xL, xLtt, nOfPts);
    CopyList(yL, yLtt, nOfPts);

    //loop through all the points in the list
    for (int i = 0; i < nOfPts; i++) {
      //loop through neighbourhood
      xLt[i] = xL[i];
      yLt[i] = yL[i];
      for (int y = -nbhood; y <= nbhood; y++)
        for (int x = -nbhood; x <= nbhood; x++)
          if (!((x == 0)&&(y == 0))) {
            energySum = (alpha*energyCont(ave,i,xL,yL,x,y,nOfPts) + beta*energyCurv(i,xL,yL,x,y,nOfPts) + gamma*energyIm(i,xL,yL,x,y));
            //check if enegry is lower
            if (energySum < E[i]) {
              E[i] = energySum;
              xLt[i] = xL[i] + x;
              yLt[i] = yL[i] + y;
            }
          }
      //update the list
      xL[i] = xLt[i];
      yL[i] = yLt[i];
      //calculate the average distance between points in the list
      ave = aveDCount(xL, yL, nOfPts);
      //calculate the new energz of neighbour points (all of points, really)
      for (int i = 0; i < nOfPts; i++)
        E[i] = (alpha*energyCont(ave,i,xL,yL,0,0,nOfPts) + beta*energyCurv(i,xL,yL,0,0,nOfPts) + gamma*energyIm(i,xL,yL,0,0));
    }

    //check how many points did change position
    nOfStaticPoints = 0;
    for (int i = 0; i < nOfPts; i++) {
      if ((xL[i] == xLtt[i])&&(yL[i] == yLtt[i]))
        nOfStaticPoints++;
    }
    //continue snake yes/no
    if (nOfStaticPoints >= StopSnake)
      StopSnakeFlag = true;

    //increase number of iterations
    nOfIterations++;

    //displaying and saving the image
      dispImage->Picture->Bitmap->PixelFormat = pf24bit;
      CopyPtrToImage(wrk, dispImage);
      dispImage->Refresh();
      //for (int i = 0; i < nOfPts; i++)
      //  DrawX(xL[i],yL[i],3,dispImage);
      for (int i = 0; i < nOfPts-1; i++)
        DrawLine(xL[i],xL[i+1],yL[i],yL[i+1],dispImage);
      DrawLine(xL[nOfPts-1],xL[0],yL[nOfPts-1],yL[0],dispImage);
      Application->ProcessMessages();
      SaveBitmapImage(dispImage, FName);

  }
  //end - MAIN SNAKE ITERATIONS

  delete xLt, yLt, xLtt, yLtt;
}
//---------------------------------------------------------------------------
//adds new points to the list of points if necessary
bool imagepro::CalculateNewListPoints(int min_distance)
{
  double edist;
  bool *results;
  bool result = false;
  results = new bool[nOfPts];
  for (int i = 0; i < nOfPts; i++)
    results[i] = false;
  int cnt = 0;

  //loop through all list points except the last one
  for (int i = 0; i < nOfPts-1; i++) {
    edist = sqrt((xL[i]-xL[i+1])*(xL[i]-xL[i+1]) + (yL[i]-yL[i+1])*(yL[i]-yL[i+1]));
    if (edist > 2*min_distance) {
      results[cnt] = AddPointInBetweenList(i);
      i++;
    }
    cnt++;
  }
  //the last list point
  edist = sqrt((xL[nOfPts-1]-xL[0])*(xL[nOfPts-1]-xL[0]) + (yL[nOfPts-1]-yL[0])*(yL[nOfPts-1]-yL[0]));
  if (edist > 2*min_distance)
    results[nOfPts-1] = AddPointInBetweenList(nOfPts-1);

  //returns true if there are no new points added in the list
  for (int i = 0; i < nOfPts; i++)
    if (results[i] == true)
      result = true;

  delete results;
  //returns true if there are no new points added in the list
  return result;
}

//---------------------------------------------------------------------------
//calculates the average distance between points
double imagepro::aveDCount(int* pX, int* pY, int l) {
  double sum = 0;
  for (int i = 0; i < (l-1); i++)
    sum += sqrt(pow(*(pX+i)-*(pX+i+1),2) + pow(*(pY+i)-*(pY+i+1),2));
  return (sum += sqrt(pow(*(pX+l)-*(pX),2) + pow(*(pY+l)-*(pY),2)))/l;
}

//---------------------------------------------------------------------------
//calculates the Continuity Energy
double imagepro::energyCont(double ave, int i, int* pX, int* pY, int xs, int ys, int l) {
 if(i>0)
   return pow(ave-sqrt(pow(*(pX+i)+xs-*(pX+i-1),2)+pow(*(pY+i)+ys-*(pY+i-1),2)),2);
 else
   return pow(ave-sqrt(pow(*(pX)+xs-*(pX+l),2)+pow(*(pY)+ys-*(pY+l),2)),2);
}
//---------------------------------------------------------------------------
//calculates the Curvature Energy
double imagepro::energyCurv(int i, int* pX, int* pY, int xs, int ys, int l) {
  if(i>0 && i<(l-1))
    return pow(*(pX+i-1)-2*(*(pX+i)+xs)+*(pX+i+1),2)+pow(*(pY+i-1)-2*(*(pY+i)+ys)+*(pY+i+1),2);
  else
    if(i == (l-1))
      return pow(*(pX+l-2)-2*(*(pX+l-1)+xs)+*(pX),2)+pow(*(pY+l-2)-2*(*(pY+l-1)+ys)+*(pY),2);
    else
      return pow(*(pX+l-1)-2*(*(pX)+xs)+*(pX+1),2)+pow(*(pY+l-1)-2*(*(pY)+ys)+*(pY+1),2);
}
//---------------------------------------------------------------------------
//calculates the Image (edge) Energy
double imagepro::energyIm(int i, int* pX, int* pY, int xs, int ys) {
  double temp;/*
  temp = 4*(ReturnImageIntensity(*(pX+i),*(pY+i)))-(ReturnImageIntensity(*(pX+i)+1,*(pY+i)));
  temp += -(ReturnImageIntensity(*(pX+i),*(pY+i))+1)-(ReturnImageIntensity(*(pX+i)-1,*(pY+i)));
  temp += -(ReturnImageIntensity(*(pX+i),*(pY+i)-1));
  temp *= -temp;*/
  temp = -pow(ReturnIntensity(*(pX+i)+xs,*(pY+i)+ys), 2);
  return temp;
}
//---------------------------------------------------------------------------
void imagepro::DisplaySnakeResult(TImage *dImage, unsigned char *ptr, int mode)
{
  CopyPtrToImage(ptr, dImage);
  dImage->Refresh();

  if (mode == 1) {
    for (int i = 0; i < nOfPts; i++)
      DrawX(xL[i],yL[i],3,dImage);
  }
  if (mode == 2) {
    for (int i = 0; i < nOfPts-1; i++)
      DrawLine(xL[i],xL[i+1],yL[i],yL[i+1],dImage);
    DrawLine(xL[nOfPts-1],xL[0],yL[nOfPts-1],yL[0],dImage);
  }
  Application->ProcessMessages();
  SaveBitmapImage(dImage, "Snap");
}
//---------------------------------------------------------------------------

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

//---------------------------------------------------------------------------
#pragma package(smart_init)