 /*
  * Khoros: $Id: rastfill.c,v 1.1 1991/05/10 15:41:33 khoros Exp $
  */

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

 /*
  * $Log: rastfill.c,v $
 * Revision 1.1  1991/05/10  15:41:33  khoros
 * Initial revision
 *
  */ 

/*
 *----------------------------------------------------------------------
 *
 * Copyright 1990, 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 1990 by UNM */

/*

	RASTFILL.C - Library of routines used to perform bilinear
                     interpolated raster fills.

	Written: Scott Wilson

        Modified: 11-Mar-90 SRW Fixed problem in assigning segments from
                  hit list. Fixed block-skip so that partial blocks are
                  filled too. Fixed bad pointer declaration in bilin_fill_f.
                  Fixed unsorted hit glitch.
                  22-Feb-91 SRW Fixed incorrect bounds check.
*/
#include "remote_gis/lvrast.h"	
#include <math.h>

bilin_fill_c(left,right,top,bot,ptr,nr,nc)
struct line *left,*right,*top,*bot;
unsigned char *ptr;
int nr,nc;
  {
    float minx,maxx,miny,maxy,intrp(),x1,x2,x3,x4,y1,y2,y3,y4;
    int i,j,k,xlo,ylo,xhi,yhi,n;
    struct line seglist[8];

    /* First find the max and min affected row and column numbers */
    x1 = left->wx1; x2 = right->wx1; x3 = right->wx2; x4 = left->wx2;
    y1 = left->wy1; y2 = right->wy1; y3 = right->wy2; y4 = left->wy2;
    minmax(x1,x2,x3,x4,&minx,&maxx);
    minmax(y1,y2,y3,y4,&miny,&maxy);
    xlo = floor((double)minx); xhi = ceil((double)maxx);
    ylo = floor((double)miny); yhi = ceil((double)maxy);

    /* Now step thru each affected row, filling in the points on that row. */
    if (ylo == yhi) yhi += 1;
    if (ylo < 0) ylo = 0; if (yhi > nr-1) yhi = nr-1;
    for (j=ylo; j<=yhi; j++)
      {
        n = fill_seglist(seglist,left,right,top,bot,j);
        for (k=0; k<n; k++)
          {
            xlo = seglist[k].x1;
            xhi = seglist[k].x2;
            if (xlo == xhi) xhi += 1;
            if (xlo < 0) xlo = 0;
            if (xhi > nc-1) xhi = nc-1;
            for (i=xlo; i<=xhi; i++)
                *(ptr+j*nc+i) = (unsigned char)(intrp(i,seglist+k));
          }
      }
  }

int
bilin_fill_s(left,right,top,bot,ptr,nr,nc)
struct line *left,*right,*top,*bot;
short *ptr;
int nr,nc;
  {
    float minx,maxx,miny,maxy,intrp(),x1,x2,x3,x4,y1,y2,y3,y4;
    int i,j,k,xlo,ylo,xhi,yhi,n;
    struct line seglist[8];

    /* First find the max and min affected row and column numbers */
    x1 = left->wx1; x2 = right->wx1; x3 = right->wx2; x4 = left->wx2;
    y1 = left->wy1; y2 = right->wy1; y3 = right->wy2; y4 = left->wy2;
    minmax(x1,x2,x3,x4,&minx,&maxx);
    minmax(y1,y2,y3,y4,&miny,&maxy);
    xlo = floor((double)minx); xhi = ceil((double)maxx);
    ylo = floor((double)miny); yhi = ceil((double)maxy);

    /* Now step thru each affected row, filling in the points on that row. */
    if (ylo == yhi) yhi += 1;
    if (ylo < 0) ylo = 0;
    if (yhi > nr-1) yhi = nr-1;
    for (j=ylo; j<=yhi; j++)
      {
        n = fill_seglist(seglist,left,right,top,bot,j);
        for (k=0; k<n; k++)
          {
            xlo = seglist[k].x1;
            xhi = seglist[k].x2;
            if (xlo == xhi) xhi += 1;
            if (xlo < 0) xlo = 0;
            if (xhi > nc-1) xhi = nc-1;
            for (i=xlo; i<=xhi; i++)
              *(ptr+j*nc+i) = intrp(i,seglist+k);
          }
      }
  }

bilin_fill_i(left,right,top,bot,ptr,nr,nc)
struct line *left,*right,*top,*bot;
int *ptr;
int nr,nc;
  {
    float minx,maxx,miny,maxy,intrp(),x1,x2,x3,x4,y1,y2,y3,y4;
    int i,j,k,xlo,ylo,xhi,yhi,n;
    struct line seglist[8];

    /* First find the max and min affected row and column numbers */
    x1 = left->wx1; x2 = right->wx1; x3 = right->wx2; x4 = left->wx2;
    y1 = left->wy1; y2 = right->wy1; y3 = right->wy2; y4 = left->wy2;
    minmax(x1,x2,x3,x4,&minx,&maxx);
    minmax(y1,y2,y3,y4,&miny,&maxy);
    xlo = floor((double)minx); xhi = ceil((double)maxx);
    ylo = floor((double)miny); yhi = ceil((double)maxy);

    /* Now step thru each affected row, filling in the points on that row. */
    if (ylo == yhi) yhi += 1;
    if (ylo < 0) ylo = 0;
    if (yhi > nr-1) yhi = nr-1;
    for (j=ylo; j<=yhi; j++)
      {
        n = fill_seglist(seglist,left,right,top,bot,j);
        for (k=0; k<n; k++)
          {
            xlo = seglist[k].x1;
            xhi = seglist[k].x2;
            if (xlo == xhi) xhi += 1;
            if (xlo < 0) xlo = 0;
            if (xhi > nc-1) xhi = nc-1;
            for (i=xlo; i<=xhi; i++)
                *(ptr+j*nc+i) = intrp(i,seglist+k);
          }
      }
  }

bilin_fill_f(left,right,top,bot,ptr,nr,nc)
struct line *left,*right,*top,*bot;
float *ptr;
int nr,nc;
  {
    float minx,maxx,miny,maxy,intrp(),x1,x2,x3,x4,y1,y2,y3,y4;
    int i,j,k,xlo,ylo,xhi,yhi,n;
    struct line seglist[8];

    /* First find the max and min affected row and column numbers */
    x1 = left->wx1; x2 = right->wx1; x3 = right->wx2; x4 = left->wx2;
    y1 = left->wy1; y2 = right->wy1; y3 = right->wy2; y4 = left->wy2;
    minmax(x1,x2,x3,x4,&minx,&maxx);
    minmax(y1,y2,y3,y4,&miny,&maxy);
    xlo = floor((double)minx); xhi = ceil((double)maxx);
    ylo = floor((double)miny); yhi = ceil((double)maxy);

    /* Now step thru each affected row, filling in the points on that row. */
    if (ylo == yhi) yhi += 1;
    if (ylo < 0) ylo = 0;
    if (yhi > nr-1) yhi = nr-1;
    for (j=ylo; j<=yhi; j++)
      {
        n = fill_seglist(seglist,left,right,top,bot,j);
        for (k=0; k<n; k++)
          {
            xlo = seglist[k].x1;
            xhi = seglist[k].x2;
            if (xlo == xhi) xhi += 1;
            if (xlo < 0) xlo = 0;
            if (xhi > nc-1) xhi = nc-1;
            for (i=xlo; i<=xhi; i++)
                *(ptr+j*nc+i) = intrp(i,seglist+k);
          }
      }
  }

fill_seglist(seglist,left,right,top,bot,y)
struct line seglist[],*left,*right,*top,*bot;
int y; /* Which scan line we're on */
  {
    /* Scan the line segments to see which ones intersect the raster line.
       Then build horizontal line segments to describe the regions. */
    struct intersection hits[8];
    int k,m,n;

    m = 0;

    n = hit_list(hits,left,right,top,bot,y);
    for (k=0; k<n-1; k++)
      {
            /* Build a segment from the next two intersections */
            seglist[m].x1 = floor((double)(hits[k].x));
            seglist[m].x2 = ceil((double)(hits[k+1].x));
            seglist[m].i1 = hits[k].inten;
            seglist[m].i2 = hits[k+1].inten;
            m++;
      }
    return(m);
  }

int
hit_list(hits,left,right,top,bot,y)
struct intersection hits[];
struct line *left,*right,*top,*bot;
int y;
  {
    /* Check each line segment  to see if it intersects the raster line
       at y. If so, place the location of the intersection and the
       interpolated value in the corresponding intersection structure. */
    float t,wx1,wy1,wx2,wy2;
    int m;
    extern int compar();

    m = 0;

    /* Check for horizontal line segment */
    wx1 = left->wx1; wx2 = left->wx2; wy1 = left->wy1; wy2 = left->wy2;
    if (wy1 == wy2 && wy1 == y)
      {
        hits[m].seg = left;
        hits[m].y = y;
        hits[m].x = left->wx1;
        hits[m].inten = left->i1;
        m++;
        hits[m].seg = left;
        hits[m].y = y;
        hits[m].x = left->wx2;
        hits[m].inten = left->i2;
        m++;
      }
    else if (between(wy1,y,wy2))
      /* Non-horizontal line segment */
      {
        t = (y-wy1)/(wy2-wy1);
        hits[m].seg = left;
        hits[m].y = y;
        hits[m].x = wx1 + t*(wx2-wx1);
        hits[m].inten = left->i1+t*(left->i2-left->i1);
        m++;
      }

    /* Check for horizontal line segment */
    wx1 = right->wx1; wx2 = right->wx2; wy1 = right->wy1; wy2 = right->wy2;
    if (wy1 == wy2 && wy1 == y)
      {
        hits[m].seg = right;
        hits[m].y = y;
        hits[m].x = right->wx1;
        hits[m].inten = right->i1;
        m++;
        hits[m].seg = right;
        hits[m].y = y;
        hits[m].x = right->wx2;
        hits[m].inten = right->i2;
        m++;
      }
    else if (between(wy1,y,wy2))
      /* Non-horizontal line segment */
      {
        t = (y-wy1)/(wy2-wy1);
        hits[m].seg = right;
        hits[m].y = y;
        hits[m].x = wx1 + t*(wx2-wx1);
        hits[m].inten = right->i1+t*(right->i2-right->i1);
        m++;
      }

    /* Check for horizontal line segment */
    wx1 = top->wx1; wx2 = top->wx2; wy1 = top->wy1; wy2 = top->wy2;
    if (wy1 == wy2 && wy1 == y)
      {
        hits[m].seg = top;
        hits[m].y = y;
        hits[m].x = top->wx1;
        hits[m].inten = top->i1;
        m++;
        hits[m].seg = top;
        hits[m].y = y;
        hits[m].x = top->wx2;
        hits[m].inten = top->i2;
        m++;
      }
    else if (between(wy1,y,wy2))
      /* Non-horizontal line segment */
      {
        t = (y-wy1)/(wy2-wy1);
        hits[m].seg = top;
        hits[m].y = y;
        hits[m].x = wx1 + t*(wx2-wx1);
        hits[m].inten = top->i1+t*(top->i2-top->i1);
        m++;
      }

    /* Check for horizontal line segment */
    wx1 = bot->wx1; wx2 = bot->wx2; wy1 = bot->wy1; wy2 = bot->wy2;
    if (wy1 == wy2 && wy1 == y)
      {
        hits[m].seg = bot;
        hits[m].y = y;
        hits[m].x = bot->wx1;
        hits[m].inten = bot->i1;
        m++;
        hits[m].seg = bot;
        hits[m].y = y;
        hits[m].x = bot->wx2;
        hits[m].inten = bot->i2;
        m++;
      }
    else if (between(wy1,y,wy2))
      /* Non-horizontal line segment */
      {
        t = (y-wy1)/(wy2-wy1);
        hits[m].seg = bot;
        hits[m].y = y;
        hits[m].x = wx1 + t*(wx2-wx1);
        hits[m].inten = bot->i1+t*(bot->i2-bot->i1);
        m++;
      }

    qsort(hits,m,sizeof(struct intersection),compar);

    return(m);
  }

int compar(a,b)
struct intersection *a,*b;
  {
    if (a->x > b->x) return(1);
    else if (a->x == b->x) return(0);
    else return(-1);
  }
