/******************************************************************************
 Warp 0.01 -- image filter plug-in for The Gimp program                        
 Copyright (C) 1996 John P. Beale                                              
 Mostly copied from LIC:Line Integral Convolution, Copyright (C) 1996 Tom Bech 
===============================================================================
 John's E-mail: beale@best.com
 Tom's  E-mail: tomb@ii.uib.no                                                 
 You can contact the original The Gimp authors at gimp@xcf.berkeley.edu        
===============================================================================
 This program is free software; you can redistribute it and/or modify it under 
 the terms of the GNU General Public License as published by the Free Software 
 Foundation; either version 2 of the License, or (at your option) any later    
 version.                                                                      
===============================================================================
 This program is distributed in the hope that it will be useful, but WITHOUT   
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 
 FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
===============================================================================
 You should have received a copy of the GNU General Public License along with  
 this program; if not, write to the Free Software Foundation, Inc., 675 Mass   
 Ave, Cambridge, MA 02139, USA.                                                
===============================================================================
 In other words, you can't sue me for whatever happens while using this ;)     
*******************************************************************************
 Changes (post 0.01):                                                          
*******************************************************************************
 This plug-in implements image warping                                         
******************************************************************************/

/*
@(GIMP)         = <plug-in warp "Effects/Warp">
@(GIMP_DEP)     = <warp.c>
@(GIMP_OBJ)     = <warp.o>
@(GIMP_LIB)     = <m c gimp>
@(GIMP_AUTHOR)  = <John Beale>
@(GIMP_EMAIL)   = <beale@best.com>
@(GIMP_DESC)    = <Implements vector image displacement (warping)>
@(GIMP_VERSION) = <0.01>
@(GIMP_URL)     = <http://www.best.com/~beale/>
*/

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "gimp.h"

/************/
/* Typedefs */
/************/

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#define D double

#define CHECKBOUNDS(x,y) (x>=0 && y>=0 && x<width && y<height)

#define EPSILON 1.0e-5

#define NUMX    40              /* Pseudo-random vector grid size */
#define NUMY    40

#define stepx   0.5
#define stepy   0.5

/***********************/
/* Some useful structs */
/***********************/

typedef struct
{
  double r,g,b;
} RGBPixel;

/*****************************/
/* Global variables and such */
/*****************************/

RGBPixel White = {1.0,1.0,1.0};
RGBPixel Black = {0.0,0.0,0.0};
RGBPixel BackGround;

double G[NUMX][NUMY][2];

double L=10.0,dx=2.0,dy=2.0,MAXV=2.5,ISteps=20.0;

Image input,output,effect;
unsigned char *dinput=NULL,*doutput=NULL,*temp_data=NULL, *scalarfield;
long progress_counter=0;            /* GIMP progress counter */

int width,height,channels,modulo;
int MainDialogID=-1;

long maxcounter;

long Effect_Toggles[3] = {0,0,1};
long Operator_Toggles[2] = {0,1};
long Convolve_Toggles[2] = {0,1};
long EffectImage=0;


long ScaleFac=10;     /* vector displacement scaling factor */
long VectorRot=-90;   /* pointwise vector rotation amount: -90,90 = twist */
long IntSteps=20;     /* number of iterations to do */
long DitherScale=2;
long MaxV=5;

/************************/
/* Convenience routines */
/************************/

/* RGBAdd(a,b): given two pointers to RGB values,
   form a = a + b
 */

void RGBAdd(RGBPixel *a,RGBPixel *b)
{
  a->r=a->r+b->r;
  a->g=a->g+b->g;
  a->b=a->b+b->b;
}

/* RGBMul(a,b): a = a*b where b is
   a scalar value
 */

void RGBMul(RGBPixel *a,double b)
{
  a->r=a->r*b;
  a->g=a->g*b;
  a->b=a->b*b;
}

/* force R,G,B to be in interval [0..1] */
void RGBClamp(RGBPixel *a)
{
  if (a->r>1.0) a->r=1.0;
  if (a->g>1.0) a->g=1.0;
  if (a->b>1.0) a->b=1.0;
  if (a->r<0.0) a->r=0.0;
  if (a->g<0.0) a->g=0.0;
  if (a->b<0.0) a->b=0.0;
}

void SetColor(RGBPixel *a,double r,double g,double b)
{
  a->r=r; a->g=g; a->b=b;
}

long int XYToIndex(int x,int y)
{
  return((long int)x*(long int)channels+(long int)y*(long int)modulo);
}

int CheckBounds(int x,int y)
{
  if (x<0 || y<0 || x>width-1 || y>height-1) return(1);
  else return(0);
}

unsigned char PeekMap(unsigned char *DispMap,int x,int y)
{
  long int index;

  index=(long int)x+(long int)width*(long int)y;

  return(DispMap[index]);
}

void PokeMap(unsigned char *DispMap,int x,int y,unsigned char value)
{
  long int index;

  index=(long int)x+(long int)width*(long int)y;

  DispMap[index]=value;
}

RGBPixel Peek(unsigned char *data,int x,int y)
{
  long int index=XYToIndex(x,y);
  RGBPixel color;

  color.r=((double)data[index])/255.0;
  color.g=((double)data[index+1])/255.0;
  color.b=((double)data[index+2])/255.0;
  return(color);
}

void Poke(unsigned char *data,int x,int y,RGBPixel *color)
{ 
  long int index=XYToIndex(x,y);

  data[index]=(unsigned char)(255.0*color->r);
  data[index+1]=(unsigned char)(255.0*color->g);
  data[index+2]=(unsigned char)(255.0*color->b);
}

/*************/
/* Main part */
/*************/

/***************************************************/
/* Compute the derivative in the x and y direction */
/* We use these convolution kernels:               */
/*     |1 0 -1|     |  1   2   1|                  */
/* DX: |2 0 -2| DY: |  0   0   0|                  */
/*     |1 0 -1|     | -1  -2  -1|                  */
/* (It's a varation of the Sobel kernels, really)  */
/***************************************************/

double gradx(unsigned char *image,int x,int y)
{
  double val=0;

  if (CHECKBOUNDS(x-1,y-1)) val=val+(D)PeekMap(image,x-1,y-1);
  if (CHECKBOUNDS(x+1,y-1)) val=val-(D)PeekMap(image,x+1,y-1);

  if (CHECKBOUNDS(x-1,y)) val=val+2*(D)PeekMap(image,x-1,y);
  if (CHECKBOUNDS(x+1,y)) val=val-2*(D)PeekMap(image,x+1,y);

  if (CHECKBOUNDS(x-1,y+1)) val=val+(D)PeekMap(image,x-1,y+1);
  if (CHECKBOUNDS(x+1,y+1)) val=val-(D)PeekMap(image,x+1,y+1);

  return(val);
}

double grady(unsigned char *image,int x,int y)
{
  double val=0;

  if (CHECKBOUNDS(x-1,y-1)) val=val+(D)PeekMap(image,x-1,y-1);
  if (CHECKBOUNDS(x,y-1)) val=val+2*(D)PeekMap(image,x,y-1);
  if (CHECKBOUNDS(x+1,y-1)) val=val+(D)PeekMap(image,x+1,y-1);

  if (CHECKBOUNDS(x-1,y+1)) val=val-(D)PeekMap(image,x-1,y+1);
  if (CHECKBOUNDS(x,y+1)) val=val-2*(D)PeekMap(image,x,y+1);
  if (CHECKBOUNDS(x+1,y+1)) val=val-(D)PeekMap(image,x+1,y+1);

  return(val);
}

/************************************/
/* A nice 2nd order cubic spline :) */
/************************************/

double cubic(double t)
{
  double at=fabs(t);
  if (at<1.0) return 2.0*at*at*at-3.0*at*at+1.0;
  else return 0.0;
}

double omega(double u,double v,int i,int j)
{
  while (i<0) i+=NUMX;
  while (j<0) j+=NUMY;
  i%=NUMX;
  j%=NUMY;
  return cubic(u)*cubic(v)*(G[i][j][0]*u+G[i][j][1]*v);
}

/*************************************************************/
/* The noise function (2D variant of Perlins noise function) */
/*************************************************************/

double noise(double x,double y)
{
  int i,sti=(int)floor(x/dx);
  int j,stj=(int)floor(y/dy);
  double sum=0.0;

  /* Calculate the double sum */
  /* ======================== */

  for (i=sti; i<=sti+1; i++)
	 {
		for (j=stj; j<=stj+1; j++)
		  sum+=omega((x-(double)i*dx)/dx,(y-(double)j*dy)/dy,i,j);
	 }

  return(sum);
}

/*************************************************/
/* Generates pseudo-random vectors with length 1 */
/*************************************************/

void generatevectors(void)
{
  double alpha;
  int i,j;

  for (i=0; i<NUMX; i++)
	 {
		for (j=0; j<NUMY; j++)
		  {
			 alpha = (double)(random()%1000)*(M_PI/500.0);
			 G[i][j][0] = cos(alpha);
			 G[i][j][1] = sin(alpha);
		  }
	 }
}

/* A simple triangle filter */
/* ======================== */

double filter(double u)
{
  double f=1.0-fabs(u)/L;
  if (f<0.0) f=0.0;
  return(f);
}


/* bilinear(): do a bilinear interpolation between the four values 
   in the two-by-two array of RGB pixels passed in "p". 
   x and y are real-valued interpolation quantities which
   should be between 0 and 1, but are forced (by wrapping)
   to that range if they aren't.
   */

static RGBPixel bilinear(double x, double y, RGBPixel *p)
{
  double   m0, m1;
  double   ix, iy;
  RGBPixel v;

  x = fmod(x, 1.0);
  y = fmod(y, 1.0);

  if (x < 0) x += 1.0;

  if (y < 0) y += 1.0;

  ix = 1.0 - x;
  iy = 1.0 - y;

  /* Red */
  /* === */

  m0 = ix * p[0].r + x * p[1].r;
  m1 = ix * p[2].r + x * p[3].r;

  v.r = iy * m0 + y * m1;

  /* Green */
  /* ===== */

  m0 = ix * p[0].g + x * p[1].g;
  m1 = ix * p[2].g + x * p[3].g;

  v.g = iy * m0 + y * m1;

  /* Blue */
  /* ==== */

  m0 = ix * p[0].b + x * p[1].b;
  m1 = ix * p[2].b + x * p[3].b;

  v.b = iy * m0 + y * m1;

  return(v);
} /* bilinear */


/* Get a bilinear-interpolated RGB pixel from the input image data global "dinput"
   based on the real-valued coordinates (u,v). If the point (u,v) is outside
   the image window, wrap it around until it is inside.
   */

void GetPixel(RGBPixel *p,double u,double v)
{
  register int x1, y1, x2, y2;
  static RGBPixel pp[4];
 
  x1 = (int)u;
  y1 = (int)v;

  if (x1 < 0) x1 = width - (-x1 % width);
  else        x1 = x1 % width;
  
  if (y1 < 0) y1 = height - (-y1 % height);
  else        y1 = y1 % height;

  x2 = (x1 + 1) % width;
  y2 = (y1 + 1) % height;
  /*
  pp[0] = Peek(dinput, x1, y1);
  pp[1] = Peek(dinput, x2, y1);
  pp[2] = Peek(dinput, x1, y2);
  pp[3] = Peek(dinput, x2, y2);
  */
  pp[0] = Peek(temp_data, x1, y1);
  pp[1] = Peek(temp_data, x2, y1);
  pp[2] = Peek(temp_data, x1, y2);
  pp[3] = Peek(temp_data, x2, y2);

  *p=bilinear(u,v,pp);
}

double Maximum(double a,double b,double c)
{
  double max=a;
  if (b>max) max=b;
  if (c>max) max=c;
  return(max);
}

double Minimum(double a,double b,double c)
{ 
  double min=a;
  if (b<min) min=b;
  if (c<min) min=c;
  return(min);
}

void RGB_To_Hue(RGBPixel *col,double *hue)
{
  double max,min,delta;

  max=Maximum(col->r,col->g,col->b);
  min=Minimum(col->r,col->g,col->b);

  if (max==min) *hue=-1.0;
  else
	 {
		delta=max-min;
		if (col->r==max)
		  {
			 *hue=(col->g-col->b)/delta;
		  }
		else if (col->g==max)
		  {
			 *hue=2.0+(col->b-col->r)/delta;
		  }
		else if (col->b==max)
		  {
			 *hue=4.0+(col->r-col->g)/delta;
		  }
		*hue=*hue*60.0;
		if (*hue<0.0) *hue=*hue+360.0;
	 }
}

void RGB_To_Saturation(RGBPixel *col,double *sat)
{
  double max,min,l;

  max=Maximum(col->r,col->g,col->b);
  min=Minimum(col->r,col->g,col->b);

  if (max==min) *sat=0.0;
  else
	 {
		l=(max+min)/2.0;
		if (l<=0.5) *sat=(max-min)/(max+min);
      else *sat=(max-min)/(2.0-max-min);
	 }
}

void RGB_To_Brightness(RGBPixel *col,double *bri)
{
  double max,min;

  max=Maximum(col->r,col->g,col->b);
  min=Minimum(col->r,col->g,col->b);

  *bri=(max+min)/2.0;
}

void RGBToHue(Image image,unsigned char **map)
{
  unsigned char *data,dval,*themap;
  int w,h;
  RGBPixel color;
  double val;
  long int maxc,cnt,index;

  data=gimp_image_data(image);

  w=gimp_image_width(image);
  h=gimp_image_height(image);
  maxc=(long int)w*(long int)h;

  themap=(unsigned char *)malloc((size_t)maxc*sizeof(unsigned char));

  for (cnt=0;cnt<maxc-1;cnt++)
	 {
		index=3*cnt;
		color.r=((double)data[index])/255.0;
		color.g=((double)data[index+1])/255.0;
		color.b=((double)data[index+2])/255.0;
		
		RGB_To_Hue(&color,&val);
		if (val==-1.0) val=0.0;
		
		dval=(unsigned char)(255.0*(val/360.0));
		themap[cnt]=dval;
	 }
  *map=themap;
}

void RGBToSaturation(Image image,unsigned char **map)
{
  unsigned char *data,dval,*themap;
  int w,h;
  RGBPixel color;
  double val;
  long int maxc,cnt,index;

  data=gimp_image_data(image);

  w=gimp_image_width(image);
  h=gimp_image_height(image);
  maxc=(long int)w*(long int)h;

  themap=(unsigned char *)malloc((size_t)maxc*sizeof(unsigned char));

  for (cnt=0;cnt<maxc-1;cnt++)
	 {
		index=3*cnt;
		color.r=((double)data[index])/255.0;
		color.g=((double)data[index+1])/255.0;
		color.b=((double)data[index+2])/255.0;
		
		RGB_To_Saturation(&color,&val);
		
		dval=(unsigned char)(255.0*val);
		themap[cnt]=dval;
	 }
  *map=themap;
}

void RGBToBrightness(Image image,unsigned char **map)
{
  unsigned char *data,dval,*themap;
  int w,h;
  RGBPixel color;
  double val;
  long int maxc,cnt,index;

  data=gimp_image_data(image);

  w=gimp_image_width(image);
  h=gimp_image_height(image);
  maxc=(long int)w*(long int)h;

  themap=(unsigned char *)malloc((size_t)maxc*sizeof(unsigned char));

  for (cnt=0;cnt<maxc-1;cnt++)
	 {
		index=3*cnt;
		color.r=((double)data[index])/255.0;
		color.g=((double)data[index+1])/255.0;
		color.b=((double)data[index+2])/255.0;
		
		RGB_To_Brightness(&color,&val);
		
		dval=(unsigned char)(255.0*val);
		themap[cnt]=dval;
	 }
  *map=themap;
}

/*************/
/* Main loop */
/*************/

/* ------------------------------------------------------------------
   find gradient of "scalarfield" at each point, generate a new image
   which is composed of the old image displaced by the given amount.
   Puts the warped output image in the global byte array "doutput".
  */

void ComputeWarp()
{
  int xcount,ycount;
  RGBPixel color;
  double vx,vy,tmp;
  double nx,ny,delta_x,delta_y;
  double delta;
  double r, theta;               /* displacement in polar coordinates */   

  for (ycount=0;ycount<height;ycount++) /* loop over rows of input image */
       {
	 
	 for (xcount=0;xcount<width;xcount++) /* loop over pixels in row */
	     {
	     	 /* Get direction vector at (x,y) and normalize it */
	       	 /* ============================================== */
	       	 
	         delta = (2*((D)rand()/RAND_MAX)-1)*DitherScale;
  	       	 vx=gradx(scalarfield,xcount,ycount) + delta;
  	       	 vy=grady(scalarfield,xcount,ycount) + delta;

 		 tmp=ScaleFac/100.0;
 	         vx*=tmp; 
		 vy*=tmp;
		 r = sqrt(vx*vx + vy*vy);
		 theta = atan2(vy,vx);
		 theta += (D)VectorRot * M_PI / 180;
		 delta_x = r * cos(theta);
		 delta_y = r * sin(theta);
		 
		 nx = xcount + delta_x;
		 ny = ycount + delta_y;

		 /* now, we have displacement vector (vx,vy) at 
                    image coordinates (xcount,ycount). Get the pixel
                    from that displacement. 
		  */

		 GetPixel(&color,nx,ny);
	 			 
	       	 Poke(doutput,xcount,ycount,&color);

	       	 progress_counter++;
	      if ((progress_counter % width)==0) gimp_do_progress(progress_counter,
							 maxcounter*IntSteps);
		  }
	 }
}


/* uses global temporary image array temp_data[] */
void DoWarp(void)
{
 int i;
 long p,pmax;

  if (Convolve_Toggles[0]==1) generatevectors();

  dx=dy=(double)VectorRot;
  MAXV=((double)MaxV)/10.0;
  ISteps=(double)IntSteps;

  width=gimp_image_width(input);
  height=gimp_image_height(input);
  maxcounter=(long int)width*(long int)height;
  channels=gimp_image_channels(input);
  modulo=channels*width;

  effect=gimp_get_input_image(EffectImage);

  if (Effect_Toggles[0]==1) RGBToHue(effect,&scalarfield);
  if (Effect_Toggles[1]==1) RGBToSaturation(effect,&scalarfield);
  else                      RGBToBrightness(effect,&scalarfield);

  gimp_free_image(effect);

  if (scalarfield==NULL)
    {
      gimp_message("Warp: Couldn't allocate temporary buffer - out of memory!\n");
      return;
    }

  output=gimp_get_output_image(0);

  pmax = width*height*3;
  /*  printf("pmax=%ld VectorRot=%ld\n",pmax,VectorRot); */

  temp_data = (unsigned char *) malloc(pmax*sizeof(unsigned char));
  
  if (!output || !temp_data)
	 {
		 gimp_message("Warp: Unable to create output image!\n");
		 return;
	 }

  gimp_init_progress("Warp");

  dinput=gimp_image_data(input);
  doutput=gimp_image_data(output);
 

 
  for (p=0; p<pmax;p++) {
    temp_data[p]=dinput[p];        
  }

  progress_counter = 0;
  for (i=0; i< IntSteps;i++) {
    ComputeWarp();
    gimp_display_image(output);
    gimp_update_image(output);
    for (p=0;p<pmax;p++)
     temp_data[p]=doutput[p];
  }

  free(temp_data);
  free(scalarfield);
}

/**************************/
/* Below is only UI stuff */
/**************************/

static void main_ok_callback(int item_id, void *client_data, void *call_data)
{
  gimp_close_dialog(MainDialogID, 1);
}

static void main_cancel_callback(int item_id, void *client_data, void *call_data)
{
  gimp_close_dialog(MainDialogID, 0); 
}

static void radio_callback (int item_ID, void *client_data, void *call_data)
{
  *((long*) client_data) = *((long*) call_data);
}

static void image_menu_callback (int item_ID, void *client_data, void *call_data)
{
  *((long*) client_data) = *((long*) call_data);
}

static void scale_callback (int item_ID,void *client_data,void *call_data)
{
  *((long*) client_data) = (long)(*((long*) call_data));
}

void CreateMainDialog(void)
{
  int MainGroupID,FrameID,RowID,ColID,HueID,SaturationID;
  int BrightnessID,MenuID;
  int ScaleFacID,VectorRotID,IntStepsID,DitherScaleID;

  MainDialogID = gimp_new_dialog("Warp");
  MainGroupID = gimp_new_column_group(MainDialogID, DEFAULT, NORMAL, "Warp Plug-In");

  RowID = gimp_new_row_group(MainDialogID, MainGroupID, NORMAL, "");

  FrameID = gimp_new_frame (MainDialogID, RowID, "Control channel");
  ColID = gimp_new_row_group(MainDialogID, FrameID, RADIO, "");

  HueID = gimp_new_radio_button (MainDialogID, ColID, "Hue");
  SaturationID = gimp_new_radio_button (MainDialogID, ColID, "Saturation");
  BrightnessID = gimp_new_radio_button (MainDialogID, ColID, "Brightness"); 


  RowID = gimp_new_row_group(MainDialogID, MainGroupID, NORMAL, "");

  FrameID = gimp_new_frame (MainDialogID, RowID, "Control Parameters");
  ColID = gimp_new_row_group(MainDialogID, FrameID, NORMAL, "");

  MenuID = gimp_new_image_menu (MainDialogID, ColID, IMAGE_CONSTRAIN_RGB, 
				"Control image:");

  ScaleFacID = gimp_new_scale (MainDialogID, ColID, 1, 100, 10, 0);
  gimp_new_label(MainDialogID, ScaleFacID, "Scale Factor:");

  VectorRotID = gimp_new_scale (MainDialogID, ColID, -90, 90, -90, 0);
  gimp_new_label(MainDialogID, VectorRotID, "Vector Rotation (degrees):");

  IntStepsID = gimp_new_scale (MainDialogID, ColID, 1, 90, 10, 0);
  gimp_new_label(MainDialogID, IntStepsID, "Number of cycles:");

  DitherScaleID = gimp_new_scale (MainDialogID, ColID, 0, 100, 2, 0);
  gimp_new_label(MainDialogID, DitherScaleID, "Dither value:");


  /* Default settings */
  /* ================ */

  gimp_change_item (MainDialogID, HueID, sizeof(long), &Effect_Toggles[0]);
  gimp_change_item (MainDialogID, SaturationID, sizeof(long), &Effect_Toggles[1]);
  gimp_change_item (MainDialogID, BrightnessID, sizeof(long), &Effect_Toggles[2]);

  /* Callbacks */
  /* ========= */
  
  gimp_add_callback (MainDialogID, HueID, radio_callback, &Effect_Toggles[0]);
  gimp_add_callback (MainDialogID, SaturationID, radio_callback, &Effect_Toggles[1]);
  gimp_add_callback (MainDialogID, BrightnessID, radio_callback, &Effect_Toggles[2]);
  gimp_add_callback (MainDialogID, MenuID, image_menu_callback, &EffectImage);

  gimp_add_callback (MainDialogID, ScaleFacID, scale_callback, &ScaleFac);
  gimp_add_callback (MainDialogID, VectorRotID, scale_callback, &VectorRot);
  gimp_add_callback (MainDialogID, IntStepsID, scale_callback, &IntSteps);
  gimp_add_callback (MainDialogID, DitherScaleID, scale_callback, &DitherScale);
	
  gimp_add_callback(MainDialogID, gimp_ok_item_id(MainDialogID), main_ok_callback, NULL);
  gimp_add_callback(MainDialogID, gimp_cancel_item_id(MainDialogID), main_cancel_callback, NULL);
}

int main(int argc,char **argv)
{
  /* Standard stuff straight out of Marc's tutorial :) */
  /* ================================================= */

  if (!gimp_init(argc,argv)) return(0);

  input=gimp_get_input_image(0);
  if (!input) return(0);

  /* We deal only with RGB images */
  /* ============================ */

  if (gimp_image_type(input)!=RGB_IMAGE) gimp_message("Warp: Uses RGB type images only!\n");
  else
	 {
		/* Create a nice dialog and show it */
		/* ================================ */

		CreateMainDialog();
 		if (gimp_show_dialog(MainDialogID)) DoWarp();
	 }

  /* Standard The End stuff.. */
  /* ======================== */

  gimp_free_image(input);
  gimp_free_image(output);

  gimp_quit();
  return(0);
}