/*
 * This software is copyrighted as noted below.  It may be freely copied,
 * modified, and redistributed, provided that the copyright notice is
 * preserved on all copies.
 *
 * There is no warranty or other guarantee of fitness for this software,
 * it is provided solely "as is".  Bug reports or fixes may be sent
 * to the author, who may or may not act on them as he desires.
 *
 * You may not include this software in a program or other software product
 * without supplying the source, or without informing the end-user that the
 * source is available for no extra charge.
 *
 * If you modify this software, you should include a notice giving the
 * name of the person performing the modification, the date of modification,
 * and the reason for such modification.
 */
/*
 * slowzoom.c - perform a better zoom ( but slower ).
 *
 * Author:      Raul Rivero
 *              Mathematics Dept.
 *              University of Oviedo
 * Date:        Mon 5 Jul 1993
 * Copyright (c) 1993, Raul Rivero
 *
 */
/*
 * This code is strongly based on the "Filtered Image Rescaling", by Dale
 * Schumacher [Graphics Gems III, 8-16]. The filter used is "the bell".
 *
 * The core is the same but some parts have been re-estructure. This
 * code doesn't try to be clear, it tries to be efficient.
 *
 * Thanks, Dale!.
 *
 * Dale Schumacher
 * Minnesota
 * ( dal@syntel.mn.org )
 *
 */

#include <lug.h>
#include <lugfnts.h>

#define bell_support                    (1.5)

typedef struct {
  int      pixel;
  double   weight;
}CONTRIB;

typedef struct {
  int      n;
  CONTRIB  *p;
}CLIST;


/*
 * Functions prototype.
 */
double
bell_filter(
#ifdef USE_PROTOTYPES
        double
#endif
);

CLIST *
calculate_filter(
#ifdef USE_PROTOTYPES
        register,
        register,
        double
#endif
);

byte *
slow_zoom(
#ifdef USE_PROTOTYPES
        register byte *,
        int,
        int,
        int,
        int,
        CLIST *,
        CLIST *
#endif
);


extern int LUGverbose;

slow_adjust_bitmap( inbitmap, outbitmap, newx, newy )
bitmap_hdr *inbitmap, *outbitmap;
int newx, newy;
{
  int i;
  CLIST *xcontrib, *ycontrib;

  if ( inbitmap->magic != LUGUSED )
    error( 19 );

  /* Correct values ? */
  if ( newx < 1 || newy < 1 )
    error( 22 );

  /* If the same proportions ... */
  if ( newx == inbitmap->xsize && newy == inbitmap->ysize ) {
    copy_bitmap( inbitmap, outbitmap );
    return 0;
  }

  /*
   * Ok, we need a resampling. But we accept only true
   * color or gray scaled images.
   */
  if ( inbitmap->depth < 24 && !isagrayscaled(inbitmap) )
    error( 23 );

  /* Fill the new header */
  outbitmap->magic = LUGUSED;
  outbitmap->xsize = newx;
  outbitmap->ysize = newy;
  outbitmap->depth = inbitmap->depth;
  outbitmap->colors = inbitmap->colors;

  /* If it's a mapped image, copy the cmap */
  if ( outbitmap->depth < 24 )
    outbitmap->cmap = (byte *) create_bw_pallete();

  /*
   * All preprocess do it. Let's go !.
   */

  /* Get the horizontal mask */
  VPRINTF( stderr, "Calculating masks: horizontal" ); VFLUSH( stderr );
  xcontrib = calculate_filter( inbitmap->xsize, newx, bell_support );
  /* Get the vertical mask */
  VPRINTF( stderr, ", vertical" ); VFLUSH( stderr );
  ycontrib = calculate_filter( inbitmap->ysize, newy, bell_support );
  VPRINTF( stderr, ". Done.\n");

  /* Zooms each component ( and applies the filter ) */
  VPRINTF( stderr, "Zoooming R" ); VFLUSH( stderr );
  outbitmap->r = slow_zoom( inbitmap->r, inbitmap->xsize, inbitmap->ysize,
                                         newx, newy, xcontrib, ycontrib );
  if ( outbitmap->depth > 8 ) {
    VPRINTF( stderr, "Zoooming G" ); VFLUSH( stderr );
    outbitmap->g = slow_zoom( inbitmap->g, inbitmap->xsize, inbitmap->ysize,
                                           newx, newy, xcontrib, ycontrib );

    VPRINTF( stderr, "Zoooming B" ); VFLUSH( stderr );
    outbitmap->b = slow_zoom( inbitmap->b, inbitmap->xsize, inbitmap->ysize,
                                           newx, newy, xcontrib, ycontrib );

  }

  /* Free the masks */
  for ( i = 0; i < newx; i++ )
    Free( xcontrib[i].p );
  for ( i = 0; i < newy; i++ )
    Free( ycontrib[i].p );
  Free(xcontrib);
  Free(ycontrib);
}

double
bell_filter(t)          /* box (*) box (*) box */
double t;
{
  if ( t < 0 )
    t = -t;
  if ( t < .5 )
    return( .75 - (t * t) );
  if ( t < 1.5 ) {
    t -= 1.5;
    return( .5 * (t * t) );
  }
  return( 0. );
}

CLIST *
calculate_filter( oldsize, newsize, filter_width )
register oldsize, newsize;
double filter_width;
{
  register int i, j, n, k;
  double scale, fscale;
  double weight;
  double left, center, right;
  double aux;
  CLIST *contrib;
  int oldsize2 = oldsize << 1;
  int contributor_size;

  /* Get the zoom scaled values */
  scale = (double) newsize / (double) oldsize;
  fscale = 1.0 / scale;

  /* Pre-calculate filter contributions for a row */
  contrib = (CLIST *) Malloc( newsize * sizeof(CLIST) );
  if( scale < 1.0 ) {
    /*
     * We are minification the current axis.
     */
    /* Scale the filter */
    filter_width *= fscale; 
    /* Get the memory size for each contributor */
    contributor_size  = (int) (filter_width * 2 + 2) * sizeof(CONTRIB);
    /* Pre-calculate ... */
    for( i = 0; i < newsize; i++ ) {
      contrib[i].p = (CONTRIB *) Malloc( contributor_size );
      /* Build the borders ... */
      center = (double) i * fscale;
      aux    = center - filter_width;
      left   = CEILING( aux );
      aux    = center + filter_width;
      right  = FLOOR( aux );
      /* ... and build the filter. */
      for( j = left; j <= right; j++ ) {
        /* Get the current weight */
        weight = center - (double) j;
        weight = bell_filter(weight * scale) * scale;
        if( j < 0 ) {
          n = -j;
        }else if( j >= oldsize ) {
                 n = oldsize2 - j - 1;
              }else {
                 n = j;
              }
        k = contrib[i].n++;
        /*** fprintf( stdout, "%3d-%3d  %3d-%3d\n", newsize, i,
                      (int) (filter_width * 2 + 1), k ); ***/
        contrib[i].p[k].pixel  = n;
        contrib[i].p[k].weight = weight;
      }
    }
  }else {
    /*
     * We are magnification the current axis.
     */
    /* Get the memory size for each contributor */
    contributor_size = (int) (filter_width * 2 + 2) * sizeof(CONTRIB);
    /* Pre-calculate ... */
    for( i = 0; i < newsize; i++ ) {
      contrib[i].p = (CONTRIB *) Malloc( contributor_size );
      /* Build the borders ... */
      center = (double) i * fscale;
      aux = center - filter_width;
      left = CEILING( aux );
      aux = center + filter_width;
      right = FLOOR( aux );
      /* ... and build the filter. */
      for( j = left; j <= right; j++ ) {
        weight = center - (double) j;
        weight = bell_filter( weight );
        if( j < 0 ) {
          n = -j;
        }else if( j >= oldsize ) {
                 n = oldsize2 - j - 1;
              }else {
                 n = j;
              }
        k = contrib[i].n++;
        contrib[i].p[k].pixel  = n;
        contrib[i].p[k].weight = weight;
      }
    }
  }

  return contrib;
}

byte *
slow_zoom( base, oldx, oldy, newx, newy, xcontrib, ycontrib )
register byte *base;
int oldx, oldy, newx, newy;
CLIST *xcontrib, *ycontrib;
{
  int i, j, k;
  register byte *tmpixel, *outpixel;
  byte *tmp;
  byte *out;
  double weight;

  /*
   * First, we resample horizontally so we need a
   * temporary buffer to store it.
   */
  tmpixel = tmp = (byte *) Malloc( newx * oldy );

  /*
   * Apply the [hozizontal] filter to zoom horizontally.
   */
  VPRINTF( stderr, ", horizontally" ); VFLUSH( stderr );
  for( i = 0; i < oldy; i++ ) {
    for( j = 0; j < newx; j++) {
      weight = 0.0;
      for( k = 0; k < xcontrib[j].n; k++) {
        weight += base[xcontrib[j].p[k].pixel] * xcontrib[j].p[k].weight;
      }
      *tmpixel++ = CORRECT( weight );
    }
    base += oldx;
  }

  /* Get memory fo our output zoomed image */
  out = (byte *) Malloc( newx * newy );

  /*
   * Apply the [vertical] filter to zoom vertically.
   */
  VPRINTF( stderr, ", vertically" ); VFLUSH( stderr );
  tmpixel = tmp;
  for( i = 0; i < newx; i++, tmpixel++ ) {
    outpixel = out + i;
    for( j = 0; j < newy; j++ ) {
      weight = 0.0;
      for( k = 0; k < ycontrib[j].n; k++ ) {
        weight += *(tmpixel + newx*ycontrib[j].p[k].pixel)
                  * ycontrib[j].p[k].weight;
      }
      *outpixel = CORRECT( weight );
      outpixel += newx;
    }
  }
  VPRINTF( stderr, ", done\n" );

  /* Free the auxiliar buffer */
  Free( tmp );

  return out;
}