 /*
  * Khoros: $Id: lvexpwarp.c,v 1.2 1991/12/18 08:56:14 dkhoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lvexpwarp.c,v 1.2 1991/12/18 08:56:14 dkhoros Exp $";
#endif

 /*
  * $Log: lvexpwarp.c,v $
 * Revision 1.2  1991/12/18  08:56:14  dkhoros
 * HellPatch3
 *
  */

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1991, University of New Mexico.  All rights reserved.
 * Permission to copy and modify this software and its documen-
 * tation only for internal use in your organization is hereby
 * granted, provided that this notice is retained thereon and
 * on all copies.  UNM makes no representations as to the sui-
 * tability and operability of this software for any purpose.
 * It is provided "as is" without express or implied warranty.
 * 
 * UNM DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
 * INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FIT-
 * NESS.  IN NO EVENT SHALL UNM BE LIABLE FOR ANY SPECIAL,
 * INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY OTHER DAMAGES WHAT-
 * SOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
 * IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS
 * ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR PER-
 * FORMANCE OF THIS SOFTWARE.
 * 
 * No other rights, including, for example, the right to redis-
 * tribute this software and its documentation or the right to
 * prepare derivative works, are granted unless specifically
 * provided in a separate license agreement.
 *---------------------------------------------------------------------
 */

#include "unmcopyright.h"        /* Copyright 1991 by UNM */

/*>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>  <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 >>>>
 >>>>         File Name: lvexpwarp.c
 >>>>
 >>>>      Program Name: vexpwarp
 >>>>
 >>>> Date Last Updated: Mon Nov 11 11:28:44 1991 
 >>>>
 >>>>          Routines: lvexpwarp - the library call for vexpwarp
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


/* -library_includes */
#include "remote_gis/lvrast.h"
/* -library_includes_end */


/****************************************************************
*
* Routine Name: lvexpwarp - library call for vexpwarp
*
* Purpose:
*    
*    Perform  geometric  transformation  (warp)  using   symbolic
*    expression
*    
*    

* Input:
*    
*    image          pointer to xvimage structure to be processed
*    
*    xmap           string containing the x map function
*    
*    ymap           string containing the y map function
*    
*    xca            X coordinate of center of action
*    
*    yca            Y coordinate of center of action
*    
*    

* Output:
*    
*    image          holds the result of the warp operation.  The  out-
*                   put data type is the same as the input data type.
*    
*    Return Value:  1 on success, 0 on failure.
*    
*    

*
* Written By: Scott Wilson
*    
*    9-Feb-91 Scott Wilson - Use new xvexpr multiple stmt syntax
*    
*    11-Nov-91 Scott Wilson - Correct bad pointer  parameter  passing,
*    updated to 17-Oct-91 version from IBF toolbox
*    
*    

****************************************************************/


/* -library_def */
int
lvexpwarp(image,xmap,ymap,xca,yca)
struct xvimage *image;
char *xmap,*ymap;
float xca,yca;
/* -library_def_end */

/* -library_code */
{
    struct line right,left,top,bot;
    int i,j;
    int nc,nr;
    int x1,x2,x3,x4,y1,y2,y3,y4;
    float wx1,wx2,wx3,wx4,wy1,wy2,wy3,wy4;
    float xa,xb,xc,xd,ya,yb,yc,yd;
    char *malloc(),s[2048];
    unsigned char *cptr,*cp;
    short *sptr,*sp;
    int   *iptr,*ip;
    float *fptr,*fp;
    int nb,bands,l;

    nc = image->row_size;
    nr = image->col_size;
    nb = image->num_data_bands;

    xve_init();

    switch(image->data_storage_type)
      {
        case VFF_TYP_BIT:
          fprintf(stderr,"lvexpwarp: Cannot warp a BIT image!\n");
          return(0);
          break;
        case VFF_TYP_1_BYTE:
          cptr = (unsigned char *)malloc(nb*nr*nc*sizeof(unsigned char));
          if (cptr == NULL)
            { 
              fprintf(stderr,"lvexpwarp: Could not allocate enough memory!\n");
              return(0);
            }
          bzero((char *)cptr,nb*nr*nc*sizeof(unsigned char));
          cp = (unsigned char *)image->imagedata;
          for (i=0; i<nr-1; i++)
            {
              for (j=0; j<nc-1; j++)
                {
                  /* Locate the input points to be transformed */
                  x1 = j;    y1 = i;    /* Upper left */
                  x2 = x1+1; y2 = y1;   /* Upper right */
                  x3 = x1+1; y3 = y1+1; /* Lower right */
                  x4 = x1;   y4 = y1+1; /* Lower left */
                  /* Take care of center-of-action */
                  xa = x1 - xca; xb = x2 - xca; xc = x3 - xca; xd = x4 - xca;
                  ya = y1 - yca; yb = y2 - yca; yc = y3 - yca; yd = y4 - yca;
                  /* Compute location in output plane */
                  sprintf(s,"x=%e;y=%e;%s",xa,ya,xmap);
                  xve_eval_float(1,s,&wx1,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xa,ya,ymap);
                  xve_eval_float(1,s,&wy1,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xb,yb,xmap);
                  xve_eval_float(1,s,&wx2,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xb,yb,ymap);
                  xve_eval_float(1,s,&wy2,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xc,yc,xmap);
                  xve_eval_float(1,s,&wx3,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xc,yc,ymap);
                  xve_eval_float(1,s,&wy3,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xd,yd,xmap);
                  xve_eval_float(1,s,&wx4,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xd,yd,ymap);
                  xve_eval_float(1,s,&wy4,NULL);
                  /* Take care of center-of-action */
                  wx1 += xca; wx2 += xca; wx3 += xca; wx4 += xca;
                  wy1 += yca; wy2 += yca; wy3 += yca; wy4 += yca;
                  for (bands=0; bands<nb; bands++)
                    {
                      l = bands*nr*nc;
                      /* Build line segments */
                      make_line(&left ,x1,y1,x4,y4,wx1,wy1,wx4,wy4,
                                (float)(*(cp+l+nc*y1+x1)),
                                (float)(*(cp+l+nc*y4+x4)));
                      make_line(&right,x2,y2,x3,y3,wx2,wy2,wx3,wy3,
                                (float)(*(cp+l+nc*y2+x2)),
                                (float)(*(cp+l+nc*y3+x3)));
                      make_line(&top  ,x1,y1,x2,y2,wx1,wy1,wx2,wy2,
                                (float)(*(cp+l+nc*y1+x1)),
                                (float)(*(cp+l+nc*y2+x2)));
                      make_line(&bot  ,x4,y4,x3,y3,wx4,wy4,wx3,wy3,
                                (float)(*(cp+l+nc*y4+x4)),
                                (float)(*(cp+l+nc*y3+x3)));
                      /* Fill the output plane polygon */
                      bilin_fill_c(&left,&right,&top,&bot,cptr+l,nr,nc);
                    }
                }
            }
          /* Fix up the header information */
          free(image->imagedata);
          image->imagedata = (char *)cptr;
          break;
        case VFF_TYP_2_BYTE:
          sptr = (short *)malloc(nb*nr*nc*sizeof(short));
          if (sptr == NULL)
            { 
              fprintf(stderr,"lvexpwarp: Could not allocate enough memory!\n");
              return(0);
            }
          bzero((char *)sptr,nb*nr*nc*sizeof(short));
          sp = (short *)image->imagedata;
          for (i=0; i<nr-1; i++)
            {
              for (j=0; j<nc-1; j++)
                {
                  /* Locate the input points to be transformed */
                  x1 = j;    y1 = i;    /* Upper left */
                  x2 = x1+1; y2 = y1;   /* Upper right */
                  x3 = x1+1; y3 = y1+1; /* Lower right */
                  x4 = x1;   y4 = y1+1; /* Lower left */
                  /* Take care of center-of-action */
                  xa = x1 - xca; xb = x2 - xca; xc = x3 - xca; xd = x4 - xca;
                  ya = y1 - yca; yb = y2 - yca; yc = y3 - yca; yd = y4 - yca;
                  /* Compute location in output plane */
                  sprintf(s,"x=%e;y=%e;%s",xa,ya,xmap);
                  xve_eval_float(1,s,&wx1,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xa,ya,ymap);
                  xve_eval_float(1,s,&wy1,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xb,yb,xmap);
                  xve_eval_float(1,s,&wx2,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xb,yb,ymap);
                  xve_eval_float(1,s,&wy2,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xc,yc,xmap);
                  xve_eval_float(1,s,&wx3,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xc,yc,ymap);
                  xve_eval_float(1,s,&wy3,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xd,yd,xmap);
                  xve_eval_float(1,s,&wx4,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xd,yd,ymap);
                  xve_eval_float(1,s,&wy4,NULL);
                  /* Take care of center-of-action */
                  wx1 += xca; wx2 += xca; wx3 += xca; wx4 += xca;
                  wy1 += yca; wy2 += yca; wy3 += yca; wy4 += yca;
                  for (bands=0; bands<nb; bands++)
                    {
                      l = bands*nr*nc;
                      /* Build line segments */
                      make_line(&left ,x1,y1,x4,y4,wx1,wy1,wx4,wy4,
                                (float)(*(sp+l+nc*y1+x1)),
                                (float)(*(sp+l+nc*y4+x4)));
                      make_line(&right,x2,y2,x3,y3,wx2,wy2,wx3,wy3,
                                (float)(*(sp+l+nc*y2+x2)),
                                (float)(*(sp+l+nc*y3+x3)));
                      make_line(&top  ,x1,y1,x2,y2,wx1,wy1,wx2,wy2,
                                (float)(*(sp+l+nc*y1+x1)),
                                (float)(*(sp+l+nc*y2+x2)));
                      make_line(&bot  ,x4,y4,x3,y3,wx4,wy4,wx3,wy3,
                                (float)(*(sp+l+nc*y4+x4)),
                                (float)(*(sp+l+nc*y3+x3)));
                      /* Fill the output plane polygon */
                      bilin_fill_s(&left,&right,&top,&bot,sptr+l,nr,nc);
                    }
                }
            }
          /* Fix up the header information */
          free(image->imagedata);
          image->imagedata = (char *)sptr;
          break;
        case VFF_TYP_4_BYTE:
          iptr = (int *)malloc(nb*nr*nc*sizeof(int));
          if (iptr == NULL)
            { 
              fprintf(stderr,"lvexpwarp: Could not allocate enough memory!\n");
              return(0);
            }
          bzero((char *)iptr,nb*nr*nc*sizeof(int));
          ip = (int *)image->imagedata;
          for (i=0; i<nr-1; i++)
            {
              for (j=0; j<nc-1; j++)
                {
                  /* Locate the input points to be transformed */
                  x1 = j;    y1 = i;    /* Upper left */
                  x2 = x1+1; y2 = y1;   /* Upper right */
                  x3 = x1+1; y3 = y1+1; /* Lower right */
                  x4 = x1;   y4 = y1+1; /* Lower left */
                  /* Take care of center-of-action */
                  xa = x1 - xca; xb = x2 - xca; xc = x3 - xca; xd = x4 - xca;
                  ya = y1 - yca; yb = y2 - yca; yc = y3 - yca; yd = y4 - yca;
                  /* Compute location in output plane */
                  sprintf(s,"x=%e;y=%e;%s",xa,ya,xmap);
                  xve_eval_float(1,s,&wx1,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xa,ya,ymap);
                  xve_eval_float(1,s,&wy1,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xb,yb,xmap);
                  xve_eval_float(1,s,&wx2,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xb,yb,ymap);
                  xve_eval_float(1,s,&wy2,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xc,yc,xmap);
                  xve_eval_float(1,s,&wx3,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xc,yc,ymap);
                  xve_eval_float(1,s,&wy3,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xd,yd,xmap);
                  xve_eval_float(1,s,&wx4,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xd,yd,ymap);
                  xve_eval_float(1,s,&wy4,NULL);
                  /* Take care of center-of-action */
                  wx1 += xca; wx2 += xca; wx3 += xca; wx4 += xca;
                  wy1 += yca; wy2 += yca; wy3 += yca; wy4 += yca;
                  for (bands=0; bands<nb; bands++)
                    {
                      l = bands*nr*nc;
                      /* Build line segments */
                      make_line(&left ,x1,y1,x4,y4,wx1,wy1,wx4,wy4,
                                (float)(*(ip+l+nc*y1+x1)),
                                (float)(*(ip+l+nc*y4+x4)));
                      make_line(&right,x2,y2,x3,y3,wx2,wy2,wx3,wy3,
                                (float)(*(ip+l+nc*y2+x2)),
                                (float)(*(ip+l+nc*y3+x3)));
                      make_line(&top  ,x1,y1,x2,y2,wx1,wy1,wx2,wy2,
                                (float)(*(ip+l+nc*y1+x1)),
                                (float)(*(ip+l+nc*y2+x2)));
                      make_line(&bot  ,x4,y4,x3,y3,wx4,wy4,wx3,wy3,
                                (float)(*(ip+l+nc*y4+x4)),
                                (float)(*(ip+l+nc*y3+x3)));
                      /* Fill the output plane polygon */
                      bilin_fill_i(&left,&right,&top,&bot,iptr+l,nr,nc);
                    }
                }
            }
          /* Fix up the header information */
          free(image->imagedata);
          image->imagedata = (char *)iptr;
          break;
        case VFF_TYP_FLOAT:
          fptr = (float *)malloc(nb*nr*nc*sizeof(float));
          if (fptr == NULL)
            { 
              fprintf(stderr,"lvexpwarp: Could not allocate enough memory!\n");
              return(0);
            }
          bzero((char *)fptr,nb*nr*nc*sizeof(float));
          fp = (float *)image->imagedata;
          for (i=0; i<nr-1; i++)
            {
              for (j=0; j<nc-1; j++)
                {
                  /* Locate the input points to be transformed */
                  x1 = j;    y1 = i;    /* Upper left */
                  x2 = x1+1; y2 = y1;   /* Upper right */
                  x3 = x1+1; y3 = y1+1; /* Lower right */
                  x4 = x1;   y4 = y1+1; /* Lower left */
                  /* Take care of center-of-action */
                  xa = x1 - xca; xb = x2 - xca; xc = x3 - xca; xd = x4 - xca;
                  ya = y1 - yca; yb = y2 - yca; yc = y3 - yca; yd = y4 - yca;
                  /* Compute location in output plane */
                  sprintf(s,"x=%e;y=%e;%s",xa,ya,xmap);
                  xve_eval_float(1,s,&wx1,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xa,ya,ymap);
                  xve_eval_float(1,s,&wy1,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xb,yb,xmap);
                  xve_eval_float(1,s,&wx2,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xb,yb,ymap);
                  xve_eval_float(1,s,&wy2,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xc,yc,xmap);
                  xve_eval_float(1,s,&wx3,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xc,yc,ymap);
                  xve_eval_float(1,s,&wy3,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xd,yd,xmap);
                  xve_eval_float(1,s,&wx4,NULL);
                  sprintf(s,"x=%e;y=%e;%s",xd,yd,ymap);
                  xve_eval_float(1,s,&wy4,NULL);
                  /* Take care of center-of-action */
                  wx1 += xca; wx2 += xca; wx3 += xca; wx4 += xca;
                  wy1 += yca; wy2 += yca; wy3 += yca; wy4 += yca;
                  for(bands=0; bands<nb; bands++)
                    {
                      l = bands*nr*nc;
                      /* Build line segments */
                      make_line(&left ,x1,y1,x4,y4,wx1,wy1,wx4,wy4,
                                *(fp+l+nc*y1+x1),
                                *(fp+l+nc*y4+x4));
                      make_line(&right,x2,y2,x3,y3,wx2,wy2,wx3,wy3,
                                *(fp+l+nc*y2+x2),
                                *(fp+l+nc*y3+x3));
                      make_line(&top  ,x1,y1,x2,y2,wx1,wy1,wx2,wy2,
                                *(fp+l+nc*y1+x1),
                                *(fp+l+nc*y2+x2));
                      make_line(&bot  ,x4,y4,x3,y3,wx4,wy4,wx3,wy3,
                                *(fp+l+nc*y4+x4),
                                *(fp+l+nc*y3+x3));
                      /* Fill the output plane polygon */
                      bilin_fill_f(&left,&right,&top,&bot,fptr+l,nr,nc);
                    }
                }
            }
          /* Fix up the header information */
          free(image->imagedata);
          image->imagedata = (char *)fptr;
          break;
        case VFF_TYP_COMPLEX:
          {
            fprintf(stderr,"lvexpwarp: Data type not supported!\n");
            return(0);
          }
        case VFF_TYP_DOUBLE:
          {
            fprintf(stderr,"lvexpwarp: Data type not supported!\n");
            return(0);
          }
        case VFF_TYP_DCOMPLEX:
          {
            fprintf(stderr,"lvexpwarp: Data type not supported!\n");
            return(0);
          }
      }
    return(1);
  }
/* -library_code_end */
