 /*
  * Khoros: $Id: lvwarp.c,v 1.1 1991/05/10 15:43:20 khoros Exp $
  */

#if !defined(lint) && !defined(SABER)
static char rcsid[] = "Khoros: $Id: lvwarp.c,v 1.1 1991/05/10 15:43:20 khoros Exp $";
#endif

 /*
  * $Log: lvwarp.c,v $
 * Revision 1.1  1991/05/10  15:43:20  khoros
 * Initial revision
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * 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 too 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: lvwarp.c
 >>>>
 >>>>      Program Name: vwarp
 >>>>
 >>>> Date Last Updated: Thu Mar 28 00:21:15 1991 
 >>>>
 >>>>          Routines: lvwarp - the library call for vwarp
 >>>>
 >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>   <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<*/


#include "vinclude.h"


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


/****************************************************************
*
* Routine Name: lvwarp - library call for vwarp
*
* Purpose:
*    
*    Perform polynomial geometric transformation (warp)
*    
*    
* Input:
*    
*    image          pointer to xvimage structure to be processed
*    
*    coeffs1        pointer to xvimage structure containing X  polyno-
*                   mial coefficients image
*    
*    coeffs2        pointer to xvimage structure containing Y  polyno-
*                   mial coefficients image
*    
*    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
*    
*    12-Feb-91 Scott Wilson -  Fixed  incorrect  intensity  parameters
*    supplied to build the top and left line structures. MUCH BETTER!
*    
*    
****************************************************************/


/* -library_def */
int
lvwarp(image,coeffs1,coeffs2,xca,yca)
struct xvimage *image,*coeffs1,*coeffs2;
float xca,yca;
/* -library_def_end */

/* -library_code */
{
    struct line right,left,top,bot;
    int i,j,m1,m2;
    int nc,nr,*clistx1,*clisty1,*clistx2,*clisty2;
    float *clist1,*clist2;
    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();
    unsigned char *cptr,*cp;
    short *sptr,*sp;
    int   *iptr,*ip;
    float *fptr,*fp;
    int nb,bands,l;

    /* Check for valid coefficient types */
    switch(coeffs1->data_storage_type)
      {
        case VFF_TYP_1_BYTE:
        case VFF_TYP_2_BYTE:
        case VFF_TYP_4_BYTE:
        case VFF_TYP_FLOAT :
          break;
        default:
          fprintf(stderr,"vwarp: Invalid data storage type for X coefficients.\n");
          return(0);
      }
    switch(coeffs2->data_storage_type)
      {
        case VFF_TYP_1_BYTE:
        case VFF_TYP_2_BYTE:
        case VFF_TYP_4_BYTE:
        case VFF_TYP_FLOAT :
          break;
        default:
          fprintf(stderr,"vwarp: Invalid data storage type for Y coefficients.\n");
          return(0);
      }

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

    m1 = make_clist(&clist1,&clistx1,&clisty1,coeffs1);
    m2 = make_clist(&clist2,&clistx2,&clisty2,coeffs2);

    switch(image->data_storage_type)
      {
        case VFF_TYP_BIT:
          fprintf(stderr,"vwarp: 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,"vwarp: 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 */
                  eval_polys(xa,ya,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx1,&wy1,coeffs1,coeffs2);
                  eval_polys(xb,yb,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx2,&wy2,coeffs1,coeffs2);
                  eval_polys(xc,yc,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx3,&wy3,coeffs1,coeffs2);
                  eval_polys(xd,yd,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx4,&wy4,coeffs1,coeffs2);
                  /* 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,"vwarp: 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 */
                  eval_polys(xa,ya,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx1,&wy1,coeffs1,coeffs2);
                  eval_polys(xb,yb,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx2,&wy2,coeffs1,coeffs2);
                  eval_polys(xc,yc,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx3,&wy3,coeffs1,coeffs2);
                  eval_polys(xd,yd,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx4,&wy4,coeffs1,coeffs2);
                  /* 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,"vwarp: 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 */
                  eval_polys(xa,ya,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx1,&wy1,coeffs1,coeffs2);
                  eval_polys(xb,yb,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx2,&wy2,coeffs1,coeffs2);
                  eval_polys(xc,yc,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx3,&wy3,coeffs1,coeffs2);
                  eval_polys(xd,yd,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx4,&wy4,coeffs1,coeffs2);
                  /* 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,"vwarp: 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 */
                  eval_polys(xa,ya,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx1,&wy1,coeffs1,coeffs2);
                  eval_polys(xb,yb,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx2,&wy2,coeffs1,coeffs2);
                  eval_polys(xc,yc,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx3,&wy3,coeffs1,coeffs2);
                  eval_polys(xd,yd,clist1,clist2,
                             clistx1,clistx2,clisty1,clisty2,
                             m1,m2,&wx4,&wy4,coeffs1,coeffs2);
                  /* 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,"vwarp: Data type not supported!\n");
            free(clist1);  free(clist2);
            free(clistx1); free(clistx2);
            free(clisty1); free(clisty2);
            return(0);
          }
      }

    free(clist1);  free(clist2);
    free(clistx1); free(clistx2);
    free(clisty1); free(clisty2);

    return(1);

  }

eval_polys(x,y,clist1,clist2,clistx1,clistx2,clisty1,clisty2,
           m1,m2,xptr,yptr,coeffs1,coeffs2)
float x,y;
float *clist1,*clist2;
int *clistx1,*clistx2,*clisty1,*clisty2;
int m1,m2;
float *xptr,*yptr;
struct xvimage *coeffs1,*coeffs2;
  {
    /* Efficiently evalutate the polynomial for the specified input values */
    int i,ncc,ncr;
    float sum,sig,xpwrs[50],ypwrs[50]; /* Max of 50th order! */
 
    ncc = coeffs1->row_size > coeffs2->row_size ? 
          coeffs1->row_size : coeffs2->row_size;
    ncr = coeffs1->col_size > coeffs2->col_size ?
          coeffs1->col_size : coeffs2->col_size;

    /* Form powers of x */
    sig = 1.0;
    xpwrs[0] = 1.0;
    for (i=1; i<ncc; i++)
      {
        sig *= x;
        xpwrs[i] = sig;
      }

    /* Form powers of y */
    sig = 1.0;
    ypwrs[0] = 1.0;
    for (i=1; i<ncr; i++)
      {
        sig *= y;
        ypwrs[i] = sig;
      }

    /* Evalute the polynomial for the new x */
    sum = 0.0;
    for (i=0; i<m1; i++) sum += clist1[i]*xpwrs[clistx1[i]]*ypwrs[clisty1[i]];
    *xptr = sum;

    /* Evalute the polynomial for the new y */
    sum = 0.0;
    for (i=0; i<m2; i++) sum += clist2[i]*xpwrs[clistx2[i]]*ypwrs[clisty2[i]];
    *yptr = sum;
  }

int
make_clist(clist,clistx,clisty,coeffs)
float **clist;
int **clistx,**clisty;
struct xvimage *coeffs;
  {
    /* Make a list of non-zero matrix elements and thier locations */
    int i,j,k,m,n,ncr,ncc;
    float f;
    char *malloc();

    ncr = coeffs->col_size;
    ncc = coeffs->row_size;

    k = ncr*ncc;

    *clist = (float *)malloc(k*sizeof(float));
    *clistx = (int *)malloc(k*sizeof(int));
    *clisty = (int *)malloc(k*sizeof(int));

    if (*clist == NULL || *clistx == NULL || *clisty == NULL)
    {
        fprintf(stderr,"vwarp: unable to allocate memory.\n");
        return(0);
    }
                
    m = 0;  /* Offset into location arrays */
    n = 0;  /* Offset into imagedata array */
    for (i=0; i< ncr; i++)
      {
        for (j=0; j<ncc; j++)
          {
            switch (coeffs->data_storage_type)
              {
                case VFF_TYP_1_BYTE:
                  f = ((unsigned char *)(coeffs->imagedata))[n++];
                  break;
                case VFF_TYP_2_BYTE:
                  f = ((short *)(coeffs->imagedata))[n++];
                  break;
                case VFF_TYP_4_BYTE:
                  f = ((int *)(coeffs->imagedata))[n++];
                  break;
                case VFF_TYP_FLOAT:
                  f = ((float *)(coeffs->imagedata))[n++];
                  break;
                default:
                  fprintf(stderr,"make_clist: Invalid data type\n");
                  return(0);
                  break;
              }
            if (f != 0.0)
              {
                (*clist)[m] = f;
                (*clistx)[m] = j;
                (*clisty)[m] = i;
                m++;
              }
          }
      }
    return(m);
  }

/* -library_code_end */
