/**************************************************************************
  Copyright (C) 1992 Guy Moreillon
  All rights reserved.

  This software may be freely copied, modified, and redistributed
  provided that this copyright notice is preserved on all copies.

  You may not distribute this software, in whole or in part, as part of
  any commercial product without the express consent of the authors.

  There is no warranty or other guarantee of fitness of this software
  for any purpose.  It is provided solely "as is".
**************************************************************************/
/**************************************************************************

				RAD
		       (Interactive Radiosity)
		    	    Version 1.0
			       1992		    	    
		    	    
		  
		    	  Guy Moreillon
			  Angelo Mangili
		    	  
**************************************************************************/

/**************************************************************************
  Fichier	: soft_zoom.c
  Description	: Routines de zoom d'image
**************************************************************************/

#include <math.h>
#include "soft_zoom.h"


#define EPSILON		0.0001

#define SHIFT  12
#define ONE    (1<<SHIFT)


typedef struct FILTER {
    int   n, w_tot;
    short *data;
    short *w;
} FILTER;


typedef struct RGBA {
   char A;
   char B;
   char G;
   char R;
} RGBA;



static int   filter_shape = BOX;
static float blurfactor   = 1.0;
static int   firsted      = 0;


/* --- prototypes --- */

FILTER *make_filter(short *in_buffer,
                    int in_size_x, 
		    int out_size_x,
		    int *maxn);

get_row(unsigned long *image, 
        short         *row,
	int           size_x,
	int           y,
        int           z);

put_row(unsigned long *image, 
        short         *row,
	int size_x,
	int y,
	int z);


float filterrad();
float filterinteg();
float mitchell();

/* --- */


void soft_zoom(unsigned long  *in_image, 
               unsigned long  *out_image,
               int in_size_x,  int in_size_y, int in_size_z,
               int out_size_x, int out_size_y)
{
    int     y, z, current_y;
    FILTER  *filter, *x_filter, *y_filter;
    short   **filter_rows;
    short   *row;
    int     *accrow;
    int     nrows;
    int     fy;
    int     i, max;
    short   *in_buffer, *out_buffer;

    in_buffer  = (short *)malloc(in_size_x  *sizeof(short));
    out_buffer = (short *)malloc(out_size_x *sizeof(short));
    accrow     = (int *)malloc(out_size_x * sizeof(int));

    x_filter = make_filter(in_buffer, in_size_x, out_size_x, &nrows);
    y_filter = make_filter(0,         in_size_y, out_size_y, &nrows);
    filter_rows = (short **) malloc(nrows * sizeof(short *));
    for(i = 0; i < nrows; i++) {
       row = (short *)malloc(out_size_x *sizeof(short));
       filter_rows[i] = row;
    }

    for (z = 0; z < in_size_z; z++) {
       current_y = 0;
       for(y = 0; y < out_size_y; y++) {
          filter = y_filter +y;
          max = ((int)filter->data) /sizeof(short) +(filter->n-1);

          while (current_y <= max) {
             get_row(in_image, in_buffer, in_size_x, current_y++, z);
             row = filter_rows[0];
             for(i=0; i<(nrows-1); i++) filter_rows[i] = filter_rows[i+1];
             filter_rows[nrows-1] = row;
             applyx_filter(filter_rows[nrows-1], x_filter, out_size_x);
	  }
          if (filter->n == 1) {
             put_row(out_image, filter_rows[nrows-1], out_size_x, y, z);
          } else {
             bzero(accrow, out_size_x * sizeof(int));
             for(fy = 0; fy < filter->n; fy++) 
                add_row(accrow, filter_rows[fy +(nrows-1) -(filter->n-1)],
	               filter->w[fy], out_size_x);
             div_row(accrow, out_buffer, filter->w_tot, out_size_x);
             put_row(out_image, out_buffer, out_size_x, y, z);
	  }
       }
    }
}


add_row(iptr, sptr, w, n)
   int   *iptr;
   short *sptr;
   int   w, n;
{
   while(n--) {
      if(n>=7) {
	  *iptr++ += (w * *sptr++);
	  *iptr++ += (w * *sptr++);
	  *iptr++ += (w * *sptr++);
	  *iptr++ += (w * *sptr++);
	  *iptr++ += (w * *sptr++);
	  *iptr++ += (w * *sptr++);
	  *iptr++ += (w * *sptr++);
	  *iptr++ += (w * *sptr++);
	  n -= 7;
      } else {
	  *iptr++ += (w * *sptr++);
      }
   }
}


div_row(iptr, sptr, tot, n)
   int   *iptr;
   short *sptr;
   int   tot, n;
{
    while(n--) {
       if(n>=7) {
	   *sptr++ = (*iptr++)/tot;
	   *sptr++ = (*iptr++)/tot;
	   *sptr++ = (*iptr++)/tot;
	   *sptr++ = (*iptr++)/tot;
	   *sptr++ = (*iptr++)/tot;
	   *sptr++ = (*iptr++)/tot;
	   *sptr++ = (*iptr++)/tot;
	   *sptr++ = (*iptr++)/tot;
	   n -= 7;
       } else {
	   *sptr++ = (*iptr++)/tot;
       }
    }
}


#define ITEM(X, Y, TYPE) \
   ((RGBA *)image +((Y *size_x) +X))->TYPE


get_row(unsigned long *image, 
        short         *row,
	int           size_x, 
	int           y, 
        int           z)
{
    int  i;

    switch (z) {
       case 0:
	  for (i = 0; i < size_x; i++) row[i] = (short)ITEM(i, y, R); 
	  break;
       case 1:
	  for (i = 0; i < size_x; i++) row[i] = (short)ITEM(i, y, G); 
	  break;
       case 2:
	  for (i = 0; i < size_x; i++) row[i] = (short)ITEM(i, y, B);
	  break;
       case 3:
	  for (i = 0; i < size_x; i++) row[i] = (short)ITEM(i, y, A);
	  break;
       default:
	  printf("tmp_zoom: invalid z size\n");
	  exit(0);
    }
}


put_row(unsigned long *image, 
        short         *row,
	int size_x,
	int y,
	int z)
{
   int i;

    switch (z) {
       case 0:
	  for (i = 0; i < size_x; i++) ITEM(i, y, R) = (char)row[i];
	  break;
       case 1:
	  for (i = 0; i < size_x; i++) ITEM(i, y, G) = (char)row[i];
	  break;
       case 2:
	  for (i = 0; i < size_x; i++) ITEM(i, y, B) = (char)row[i];
	  break;
       case 3:
	  for (i = 0; i < size_x; i++) ITEM(i, y, A) = (char)row[i];
	  break;
       default:
	  printf("tmp_zoom: invalid z size\n");
	  exit(0);
    }
}

#undef ITEM(X, Y, TYPE)



FILTER *make_filter(short *in_buffer,
                    int in_size_x, 
		    int out_size_x,
		    int *maxn)
{
    FILTER *f, *filter;
    int    x, n;
    float  bmin, bmax, bcent, brad;
    float  fmin, fmax, acent, arad;
    int    amin, amax;
    float  cover;

    f = filter = (FILTER *) malloc(out_size_x *sizeof(FILTER));
    *maxn = 0;
    for(x = 0; x < out_size_x; x++) {
       if (out_size_x < in_size_x) {
	  brad  = filterrad()     /out_size_x;
	  bcent = ((float)x +0.5) /out_size_x;
	  amin  = (int)floor((bcent -brad) *in_size_x +EPSILON);
	  amax  = (int)floor((bcent +brad) *in_size_x -EPSILON);
	  if(amin <  0)         amin = 0;
	  if(amax >= in_size_x) amax = in_size_x -1;
	  f->n     = 1 +amax -amin;
	  f->data  = in_buffer +amin;
	  f->w     = (short *)malloc(f->n *sizeof(short));
	  f->w_tot = 0;
	  for (n = 0; n < f->n; n++) {
	     bmin      = out_size_x *((((float)amin +n)   /in_size_x) -bcent);
	     bmax      = out_size_x *((((float)amin +n +1)/in_size_x) -bcent);
	     cover     = filterinteg(bmax, 1) -filterinteg(bmin, 0);  
	     f->w[n]   = (ONE *cover) +0.5;
	     f->w_tot += f->w[n];
	  }
       } else {
	  arad = filterrad() /in_size_x;
	  bmin = ((float)x    ) /out_size_x;
	  bmax = ((float)x+1.0) /out_size_x;
	  amin = (int)floor((bmin-arad) *in_size_x +0.5 +EPSILON);
	  amax = (int)floor((bmax+arad) *in_size_x -0.5 -EPSILON);
	  if(amin <  0)         amin = 0;
	  if(amax >= in_size_x) amax = in_size_x-1;
	  f->n     = 1 + amax - amin;
	  f->data  = in_buffer + amin;
	  f->w     = (short *)malloc(f->n * sizeof(short));
	  f->w_tot = 0;
	  for (n=0; n<f->n; n++) {
	     acent     = (amin +n +0.5) /in_size_x;
	     fmin      = in_size_x *(bmin-acent);
	     fmax      = in_size_x *(bmax-acent);
	     cover     = filterinteg(fmax, 1) -filterinteg(fmin, 0);  
	     f->w[n]   = (ONE *cover) +0.5;
	     f->w_tot += f->w[n];
	  }
       }
       if(f->n>*maxn) *maxn = f->n;
       f++;
    }
    return (filter);
}


applyx_filter(out_buffer, xfilt, out_size_x)
   register short *out_buffer;
   register FILTER *xfilt;
   register int out_size_x;
{
    register short *w;
    register short *dptr;
    register int n, val;

    while(out_size_x--) {
   if((n=xfilt->n) == 1) {
       *out_buffer++ = *xfilt->data;
   } else {
       w = xfilt->w;
       dptr = xfilt->data;
       val = 0;
       n = xfilt->n;
       while(n--) 
            val += *w++ * *dptr++;
       *out_buffer++ = val / xfilt->w_tot;
   }
   xfilt++;
    }
}


float filterrad()
{
    switch(filter_shape) {
       case BOX:
          return (0.5 *blurfactor);
          break;
       case TRIANGLE:
          return (1.0 *blurfactor);
          break;
       case QUADRATIC:
          return (1.0 *blurfactor);
          break;
      case MITCHELL:
          return (2.0*blurfactor);
          break;
      default:
          printf("soft_zoom: bad choice in filterrad()\n");
          exit(0);
    }
}


float quadinteg(x)
   float x;
{
   if(x<-1.0)
       return 0.0;
   if(x<-0.5)
       return 2.0*((1.0/3.0)*(x*x*x+1.0) + x*x + x);
   else
       return -(2.0/3.0)*x*x*x + x + 0.5;
}


float filterinteg(x,side)
   float x;
   int side;
{
    float val;

    x = x/blurfactor;
    switch(filter_shape) {
   case BOX:
       if (x<-0.5)     return 0.0;
       else if (x>0.5) return 1.0;
       else            return x+0.5;
       break;
   case TRIANGLE:
       if(x<-1.0) return 0.0;
       else if(x>1.0) return 1.0;
       else if(x<0.0) {
         val = x+1.0;
         return 0.5*val*val;
       } else {
         val = 1.0-x;
         return 1.0-0.5*val*val;
       }
       break;
   case QUADRATIC:
       if(x<0.0) return quadinteg(x);
       else return 1.0-quadinteg(-x);
       break;
   case MITCHELL:
       if(side == 0) return 0.0;
       return mitchell(x);
    }
}


float p0, p2, p3, q0, q1, q2, q3;


/*
 * see Mitchell&Netravali, "Reconstruction Filters in Computer Graphics",
 * SIGGRAPH 88.  Mitchell code provided by Paul Heckbert.
 */
float mitchell(x)   /* Mitchell & Netravali's two-param cubic */
   float x;
{
    if(!firsted) {
       mitchellinit(1.0/3.0,1.0/3.0);
       firsted = 1;
    }
    if (x<-2.) return 0.0;
    if (x<-1.) return 2.0*(q0-x*(q1-x*(q2-x*q3)));
    if (x<0.) return 2.0*(p0+x*x*(p2-x*p3));
    if (x<1.) return 2.0*(p0+x*x*(p2+x*p3));
    if (x<2.) return 2.0*(q0+x*(q1+x*(q2+x*q3)));
    return 0.0;
}

mitchellinit(b,c)
   float b, c;
{

    p0 = (  6. -  2.*b        ) / 6.;
    p2 = (-18. + 12.*b +  6.*c) / 6.;
    p3 = ( 12. -  9.*b -  6.*c) / 6.;
    q0 = (        8.*b + 24.*c) / 6.;
    q1 = (     - 12.*b - 48.*c) / 6.;
    q2 = (        6.*b + 30.*c) / 6.;
    q3 = (     -     b -  6.*c) / 6.;
}


