/****************************************************************/
/*								*/
/* coor.c:	coordinate transformations and data		*/
/*		manipulation					*/
/*								*/
/*		Design: Walter Benzing 1994			*/
/*								*/
/*		NO RIGHTS RESERVED				*/
/****************************************************************/

/* standard includes */

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>

/* my includes */

#include "xdim.h"

/*
 * allocate a new buffer in act data
 */

int CO_NewData(int NewWidth, int NewHeight, int DataType, XDimInfo *XDim)
   {
    size_t Size;
    ViewportInfo *view;
    void *buffer;

    if(XDim->port >=0)
      {
       view = &((XDim->Views)[XDim->port]);
       view->changed |= DATA_AND_CHANGE;
       switch (DataType)
          {
           case MY_UCHAR:
           Size = NewWidth * NewHeight;
           break;
           default:
	   Size = NewWidth * NewHeight * sizeof(double);
	   break;
	  }
	if(Size != 0)
	  {
	   buffer = malloc(Size);
	   if(!buffer)
	     {
	      sprintf(XDim->error, "Couldn't allocate %d bytes of memory",
		      (int)Size);
	      return(ERROR);
	     }
	  }
	else
	  {
	   NewWidth = 0;
	   NewHeight = 0;
	  }
       if(view->charData)
         {
          free(view->charData);
          view->charData = NULL;
         }
       if(view->actData)
         {
          free(view->actData);
          view->actData = NULL;
         }
       switch (DataType)
          {
           case MY_UCHAR:
           view->charData = (unsigned char *) buffer;;
           break;
           default:
	   view->actData = (double *) buffer;;
	   break;
	  }
       view->Width = NewWidth;
       view->Height = NewHeight;
       view->dataType = DataType;
      }
    else
      {
       strcpy(XDim->error, "No viewport available");
       return(ERROR);
      }
    return(OK);
   }

/*
 * find Max, Min
 */

void CO_MinMax(double *max, double *min, int Size, double *Data)
   {
    int i;
    double Max, Min, *dp, val;

    dp = Data;
    Max = *dp;
    Min = *dp++;
    for(i=1; i<Size; i++)
       {
	val = *dp++;
	if(val<Min)
	  Min = val;
	else if(val>Max)
	  Max = val;
       }
    *max = Max;
    *min = Min;
   }

void CO_MinMaxChar(double *max, double *min, int Size, unsigned char *CData)
   {
    int i;
    unsigned char Max, Min, *cp, val;

    cp = CData;
    Max = *cp;
    Min = *cp++;
    for(i=1; i<Size; i++)
       {
	val = *cp++;
	if(val<Min)
	  Min = val;
	else if(val>Max)
	  Max = val;
       }
    *max = (double)Max;
    *min = (double)Min;
   }

/*
 * calculate the base for the axis
 */

double CO_Base(double Max)
   {
    double base;
    
    Max = fabs(Max);
    if(Max > 0.0)
      {
       base = log10(fabs(Max));
       base = floor(base/3.0+0.2)*3.0;
       base = exp(log(10.0)*base);
      }
    else
      base = 1.0;
    return(base);
   }

/*
 * Type cast to double
 */

int CO_Cast2Double(ViewportInfo *view, XDimInfo *XDim)
   {
    int i, Size;
    unsigned char *cp, *pp;
    double *dp;

    switch (view->dataType)
     {
      case MY_UCHAR:
      if(!view->charData)
	{
	 sprintf(XDim->error, "No %s Data available for cast to double!",
	    TypeNames[MY_UCHAR]);
	 return(ERROR);
	}

      /* CAUTION old data is removed with CO_NewData */
      
      pp = view->charData;
      view->charData = NULL;
      if(CO_NewData(view->Width, view->Height, MY_DOUBLE, XDim))
	 return(ERROR);
      cp = pp;
      dp = view->actData;
      Size = view->Width * view->Height;
      for(i=0; i<Size; i++)
	 {
	  *dp++ = (double)*cp++;
	 }
      free(pp);
      view->charData = NULL;
      break;
      default:
      break;
     }
    view->dataType = MY_DOUBLE;
    return(OK);
   }

/*
 * Type cast to unsigned char
 */

int CO_Cast2UChar(ViewportInfo *view, XDimInfo *XDim)
   {
    int i, Size;
    unsigned char *cp;
    double *dp, *pp;

    switch (view->dataType)
     {
      case MY_UCHAR:
      break;
      default:
      if(!view->actData)
	{
	 sprintf(XDim->error, "No %s Data available for cast to uchar!",
	    TypeNames[MY_DOUBLE]);
	 return(ERROR);
	}

      /* CAUTION old data is removed with CO_NewData */
      
      pp = view->actData;
      view->actData = NULL;
      if(CO_NewData(view->Width, view->Height, MY_UCHAR, XDim))
	 return(ERROR);
      cp = view->charData;
      dp = view->actData;
      Size = view->Width * view->Height;
      for(i=0; i<Size; i++)
	 {
	  *cp++ = (unsigned char)*dp++;
	 }
      free(pp);
      view->actData = NULL;
      break;
     }
    view->dataType = MY_UCHAR;
    return(OK);
   }

/*
 * format from double to integer, values
 * outside XDim->Max, XDim->Min are cut off
 */
   
int CO_Data2Int(ViewportInfo *view, XDimInfo *XDim)
   {
    int i, Size;
    myCoor *ip;
    double *dp, fact, Max, Min, val;
    unsigned char *cp;

    if(view->actData || view->charData)
      {
       Size = view->Width * view->Height;
       if(view->Data)
	 free(view->Data);
       view->Data = malloc(Size*sizeof(int));
       if(!view->Data)
	 {
	  sprintf(XDim->error, "Couldn't allocate %d bytes of memory",
	     Size*(int)sizeof(int));
	  return(ERROR);
	 }
       ip = view->Data;
       Min = view->Min;
       Max = view->Max;
       if(Max > Min)
	 fact = (double)(CO_ONE-1)/(Max - Min);
       else
	 fact = 1.0;
       view->FactZ = 1.0/fact;
       switch (view->dataType)
        {
         case MY_UCHAR:
	 cp = view->charData;
	 view->dMin = 0;
	 for(i=0; i<Size; i++)
	    {
	     val = (double)*cp++;
	     val = Mminmax(val,Min,Max);
	     *ip++ = (myCoor)(fact * (val - Min));
	    }
         break;
         default:
	 dp = view->actData;
	 if(Min > 0.0 || Max < 0.0)
	   {
	    view->dMin = 0;
	    for(i=0; i<Size; i++)
	       {
		val = *dp++;
		val = Mminmax(val,Min,Max);
		*ip++ = (myCoor)(fact * (val - Min));
	       }
	   }
	 else
	   {
	    view->dMin = fact*Min;
	    for(i=0; i<Size; i++)
	       {
		val = *dp++;
		val = Mminmax(val,Min,Max);
		*ip++ = (myCoor)(fact * val);
	       }
	   }
	 break;
	}
      }
    view->changed &= ~DATA_CHANGED;
    return(OK);
   }
/*
 * Return the Value depending on the Data Type
 */

double CO_GetVal(int x, int y, ViewportInfo *view)
   {
    double val;
    
    switch (view->dataType)
       {
        case MY_UCHAR:
        val = (double) *(view->charData+ view->Width*y+ x);
        break;
        default:
        val = *(view->actData+ view->Width*y+ x);
	break;
       }
    return(val);
   }

/*
 * do an interpolation
 */

double CO_Inter(double xs, double ys, ViewportInfo *view)
   {
    int del, dx, dy;
    double min,max,*dp, z1, z2, zn, fact;

    del = view->Width;
    min = view->Min;
    max = view->Max;
    dx = (int) xs;
    dy = (int) ys;
    if(view->dataType)
      zn = CO_GetVal(dx, dy, view);
    else
      {
       dp = view->actData + dy*del+ dx;
       zn = *dp;
      }
    zn = Mminmax(zn,min,max);
    fact = xs-floor(xs);
    if(dx<del-1)
      {
       if(view->dataType)
	 z1 = CO_GetVal(dx+1, dy, view);
       else
	 z1 = *(dp+1);
       z1 = Mminmax(z1,min,max);
       z1 = zn + fact*(z1-zn);
      }
    else
      z1 = zn;
    if(dy < view->Height-1)
      {
       if(view->dataType)
         {
          dy++;
	  zn = CO_GetVal(dx, dy, view);
	 }
       else
         {
	  dp += del;
	  zn = *dp;
         }
       zn = Mminmax(zn,min,max);
      }
    if(dx<del-1)
      {
       if(view->dataType)
	 z2 = CO_GetVal(dx+1, dy, view);
       else
         z2 = *(dp+1);
       z2 = Mminmax(z2,min,max);
       z2 = zn + fact*(z2-zn);
      }    
    else
      z2 = zn;
    fact = ys -floor(ys);
    zn = z1 + fact*(z2-z1);
    return(zn);
   }

/*
 * Start spline interpolation
 */

int CO_StartSpline(int len, double dx, double *a, double *z, XDimInfo *XDim)
   {
    int ret, i, len1, size;
    double fact, dx2, dx23, dx6,*b,*c,*d;
    
    size = len*sizeof(double);
    b = (double *)malloc(size);
    c = (double *)malloc(size);
    d = (double *)malloc(size);
    if(len < 3)
      {
       sprintf(XDim->error, "To do spline interpolation n must be > 2 !!!");
       free(a);
       ret = ERROR;
      }
    else if(b== NULL || c == NULL || d == NULL)
      {
       sprintf(XDim->error, "Couldn't allocate %d bytes of memory", 3*size);
       free(a);
       ret = ERROR;
      }
    else
      {
       dx2 = dx*dx;
       len1 = len-1;
       a[0] = 0.0;
       b[0] = 1.0;
       c[0] = 0.0;
       a[len1] = 0.0;
       b[len1] = 1.0;
       c[len1] = 0.0;
       d[0] = (z[1]-z[0])/dx2;
       d[len1] = (z[len1-1]-z[len1])/dx2;
       dx6 = dx2/6.0;
       dx23 = 2.0*dx2/3.0;
       for(i=1; i<len1; i++)
          {
           a[i] = dx6;
           b[i] = dx23;
           c[i] = dx6;
           d[i] = z[i+1]-2*z[i]+z[i-1];
          }
       for(i=1; i<len1; i++)
          {
           fact = -a[i]/b[i-1];
           b[i] += c[i-1]*fact;
           d[i] += d[i-1]*fact;
          }
       a[len1] = d[len1]/b[len1];
       a[0] = d[0]/b[0];
       for(i= len1-1; i>0; i--)
          {
           a[i] = (d[i]-c[i]*a[i+1])/b[i];
          } 
       ret = OK;
      }
    free(b);
    free(c);
    free(d);
    return(ret);
   }

/*
 * do spline interpolation
 */

double CO_Spline(double x, double dx, int len, double *xvec, double *yvec)
   {
    double y,xi,xi1, dx6;
    int frac;
    
    frac = ((int)x);
    xi1 = x- (double) frac;
    frac++;
    xi = x-(double) frac;
    dx6 = 1.0/dx/6.0;
    if(x < 0.0)
      {
       y = yvec[0];
      }
    else if(frac>=len)
      {
       y = yvec[len-1];
      }
    else
      {
       y = xvec[frac]*xi1*(xi1*xi1*dx6-dx/6.0)+yvec[frac]*xi1/dx-
           xvec[frac-1]*xi*(xi*xi*dx6-dx/6.0)-yvec[frac-1]*xi/dx;
      }
    return(y);
   }

/*
 * Cut a line through a grid
 */

int CO_NextStep(double x1, double y1, double *x, double *y, double *z,
	ViewportInfo *view)
   {
    int Next;
    double xs, ys, xn, yn, xn1, yn1;

    xn = *x;
    yn = *y;
    xn1 = floor(xn);
    yn1 = floor(yn);
    if(x1 >= xn1 && x1 < xn1+1.0 && y1 >= yn1 && y1 < yn1+1.0)
      Next = False;
    else
      Next = True;  
    xs = x1;
    ys = y1;
    if(xn1>=x1)
      {
       if(xn1>=xn)
         xn1 -= 1.0;
       if(xn1 >= x1)
         xs = xn1;
      }
    else
      {
       xn1 += 1.0;
       if(xn1<=x1)
          xs = xn1;
      }
    if(xn != x1)
      ys = y1 + (yn - y1)*(xs-x1)/(xn-x1);
    if(ys < yn1 || ys > yn1+1.0)
      {
       if(yn1>=y1)
         {
          if(yn1>=yn)
            yn1 -= 1.0;
          if(yn1>=y1)
	    ys = yn1;
	 }
       else
	 {
	  yn1 += 1.0;
	  if(yn1<=y1)
             ys = yn1;
	 }
       if(yn != y1)
	 xs = x1 + (xn - x1)*(ys-y1)/(yn-y1);
      }
    /* Now do Interpolation */
    
    *z = CO_Inter(xs, ys, view);
    *y = ys;
    *x = xs;
    return(Next);
   }

/*
 * Calculate Coordinates from screen x,y
 */

int CO_Point(int x, int y, double *xr, double *yr, double *zr,
	ViewportInfo *view, XDimInfo *XDim)
   {
    int next, found;
    double fact, Sub, cx, cy, cz;
    myCoor x0, y0, x1, y1, coMaxX, sinXY, cosXY, sinZ, cosZ, moveX, moveY;
    myCoor z1, MaxX, MaxY, coMaxY, zoomZ, zoomX, zoomY, height, width;
    myCoor x2, y2, x3, y3, dumx, dumy, Fact;
    double pi2, dAngle, cx1, cy1, cz1, cx2, cy2, Max, Min;
    int xa, ya, xas, yas, xb, yb, xbs, ybs, xc, yc, xcs, ycs, xd, yd, xds, yds;
    
    dAngle = 2.0*PI;
    pi2 = PI/2.0;
    sinXY = view->sinXY;
    cosXY = view->cosXY;
    sinZ = view->sinZ;
    cosZ = view->cosZ;
    moveX = view->moveX;
    moveY = view->moveY;
    zoomX = view->zoomX;
    zoomY = view->zoomY;
    zoomZ = view->zoomZ;
    MaxX = view->Width-1;
    MaxY = view->Height-1;
    width = view->pwidth;
    height = view->pheight;
    y = Mupdown(y,height);
    if(MaxY > MaxX)
      {
       if(MaxY >0)
	 Fact = MaxY;
       else
	 Fact = 1;
       coMaxY = CO_ONE >> 1;
       coMaxX = (MaxX<<(CO_ONESHIFT-1))/Fact;
      }
    else
      {
       if(MaxX >0)
	 Fact = MaxX;
       else
	 Fact = 1;
       coMaxX = CO_ONE >> 1;
       coMaxY = (MaxY<<(CO_ONESHIFT-1))/Fact;
      }


    if((dAngle -= pi2) <= view->rotXY)
      {
       x0 = -coMaxX;
       y0 = -coMaxY;
       xas = 0;
       yas = 0;
       x1 = coMaxX;
       y1 = -coMaxY;
       xbs = MaxX;
       ybs = 0;
       x2 = coMaxX;
       y2 = coMaxY;
       xcs = MaxX;
       ycs = MaxY;
       x3 = -coMaxX;
       y3 = coMaxY;
       xds = 0;
       yds = MaxY;
      }
    else if((dAngle -= pi2) <= view->rotXY)
      {
       x0 = -coMaxX;
       y0 = coMaxY;
       xas = 0;
       yas = MaxY;
       x1 = -coMaxX;
       y1 = -coMaxY;
       xbs = 0;
       ybs = 0;
       x2 = coMaxX;
       y2 = -coMaxY;
       xcs = MaxX;
       ycs = 0;
       x3 = coMaxX;
       y3 = coMaxY;
       xds = MaxX;
       yds = MaxY;
      }
    else if((dAngle -= pi2) <= view->rotXY)
      {
       x0 = coMaxX;
       y0 = coMaxY;
       xas = MaxX;
       yas = MaxY;
       x1 = -coMaxX;
       y1 = coMaxY;
       xbs = 0;
       ybs = MaxY;
       x2 = -coMaxX;
       y2 = -coMaxY;
       xcs = 0;
       ycs = 0;
       x3 = coMaxX;
       y3 = -coMaxY;
       xds = MaxX;
       yds = 0;
      }
    else
      {
       x0 = coMaxX;
       y0 = -coMaxY;
       xas = MaxX;
       yas = 0;
       x1 = coMaxX;
       y1 = coMaxY;
       xbs = MaxX;
       ybs = MaxY;
       x2 = -coMaxX;
       y2 = coMaxY;
       xcs = 0;
       ycs = MaxY;
       x3 = -coMaxX;
       y3 = -coMaxY;
       xds = 0;
       yds = 0;
      }
    CO_TransformXY(xa,ya,dumx,dumy,x0,y0,0,sinXY,cosXY,sinZ,cosZ,
	zoomX,zoomY,zoomZ,moveX,moveY,0);
    CO_TransformXY(xb,yb,dumx,dumy,x1,y1,0,sinXY,cosXY,sinZ,cosZ,
	zoomX,zoomY,zoomZ,moveX,moveY,0);
    CO_TransformXY(xc,yc,dumx,dumy,x2,y2,0,sinXY,cosXY,sinZ,cosZ,
	zoomX,zoomY,zoomZ,moveX,moveY,0);
    CO_TransformXY(xd,yd,dumx,dumy,x3,y3,0,sinXY,cosXY,sinZ,cosZ,
	zoomX,zoomY,zoomZ,moveX,moveY,0);
    if(x<xd || x> xb)
       return(False);;
    if(x < xa)
      {
       if(xa != xd)
         {
          cx1 = xds + ((double)((x-xd)*(xas-xds)))/(xa-xd);
          cy1 = yds + ((double)((x-xd)*(yas-yds)))/(xa-xd);
         }
       else
         {
          cx1 = xds;
          cy1 = yds;
         }
      }
    else
      {
       if(xb != xa)
         {
          cx1 = xas + ((double)((x-xa)*(xbs-xas)))/(xb-xa);
          cy1 = yas + ((double)((x-xa)*(ybs-yas)))/(xb-xa);
         }
       else
         {
          cx1 = xas;
          cy1 = yas;
         }
      }
    if(x < xc)
      {
       if(xc != xd)
         {
          cx2 = xds + ((double)((x-xd)*(xcs-xds)))/(xc-xd);
          cy2 = yds + ((double)((x-xd)*(ycs-yds)))/(xc-xd);
         }
       else
         {
          cx2 = xds;
          cy2 = yds;
         }
      }
    else
      {
       if(xb != xc)
         {
          cx2 = xcs + ((double)((x-xc)*(xbs-xcs)))/(xb-xc);
          cy2 = ycs + ((double)((x-xc)*(ybs-ycs)))/(xb-xc);
         }
       else
         {
          cx2 = xcs;
          cy2 = ycs;
         }
      }
    cz1 = CO_Inter(cx1, cy1, view);
    cx = cx1;
    cy = cy1;
    cz = cz1;
    Min = view->Min;
    Max = view->Max;
    if(Max > Min)
      fact = (CO_ONE-1)/(Max - Min);
    else
      fact = 1.0;
    x1 = (myCoor)((cx/Fact)*CO_ONE)-coMaxX;
    y1 = (myCoor)((cy/Fact)*CO_ONE)-coMaxY;
    if(Min <0.0 && Max > 0.0)
      Sub = 0.0;
    else
      Sub = Min;
    z1 = (myCoor)(fact * (cz - Sub));
    CO_TransformXY(xb,yb,dumx,dumy,x1,y1,z1,sinXY,cosXY,sinZ,cosZ,
	zoomX,zoomY,zoomZ,moveX,moveY,0);
    found = False;
    do
      {
       cx1 = cx;
       cy1 = cy;
       cz1 = cz;
       ya = yb;
       next=CO_NextStep(cx2,cy2,&cx,&cy,&cz, view);
       x1 = (myCoor)((cx/Fact)*CO_ONE)-coMaxX;
       y1 = (myCoor)((cy/Fact)*CO_ONE)-coMaxY;
       z1 = (myCoor)(fact * (cz - Sub));
       CO_TransformXY(xb,yb,dumx,dumy,x1,y1,z1,sinXY,cosXY,sinZ,cosZ,
	zoomX,zoomY,zoomZ,moveX,moveY,0);
       if((y <= yb && y >= ya) || (y <= ya && y >= yb))
         {
          found = True;
          if(ya!=yb)
            {
	     *xr = cx1 + ((cx-cx1)*(y-ya))/(yb-ya);
	     *yr = cy1 + ((cy-cy1)*(y-ya))/(yb-ya);
	     *zr = cz1 + ((cz-cz1)*(y-ya))/(yb-ya);
	    }
	  else
	    {
	     *xr = cx1;
	     *yr = cy1;
	     *zr = cz1;
	    }            
          break;
         }
      }
    while(next);
    return(found);
   }

/*
 * Calculate Axis Base values
 */

void CO_CalculateBases(ViewportInfo *view2)
   {
    switch (view2->dataType)
      {
       case MY_UCHAR:
       CO_MinMaxChar(&(view2->Max), &(view2->Min),
	  view2->Width * view2->Height,view2->charData);
       break;
       default:
       CO_MinMax(&(view2->Max), &(view2->Min),
	  view2->Width * view2->Height,view2->actData);
       break;
      }
    view2->baseX = CO_Base(view2->Width*view2->FactX
	+view2->OffX);
    view2->baseY = CO_Base(view2->Height*view2->FactY
	+view2->OffY);
    if(fabs(view2->Max) > fabs(view2->Min))
      view2->baseZ = CO_Base(view2->Max);
    else
      view2->baseZ = CO_Base(view2->Min);
   }
