#include "tmachine.h"

#include "tpixelgr.h"

#include "quickputP.h"

//#include "tspecialstyleid.h"

#include "tsystem.h"

#include "tcolorstyles.h"

#include "tpixelutils.h"

//#include "tstopwatch.h"

#ifndef TNZCORE_LIGHT

#include "tpalette.h"

#include "trastercm.h"

#include "tropcm.h"

#endif

using namespace TConsts;

#ifdef _WIN32

#include <emmintrin.h> // per SSE2</emmintrin.h>

#endif

#include <memory></memory>

//===========================================================================

/*

Versione con estensione dell'ultimo pixel e con default_value

per pixel "fuori" dall'immagine di ingresso.

Sistemi di coordinate a meno di una traslazione:

UV: coordinate dell'immagine di partenza

ST: coordinate di filtro (il raggio del filtro e' intero in ST)

FG: coordinate del filtro discretizzato

XY: coordinate dell'immagine di arrivo

Tra UV e ST c'e' una traslazione intera finche' non c'e' ingrandimento

e non c'e' blur. Il blur aggiunge uno scale, e altrettanto fa l'ingrandimento.

Tra ST e FG c'e' uno scale per la risoluzione del filtro.

Oggetti:

out : pixel di output (centro di ST e FG)

ref : pixel di riferimento dell'immagine di input

pix : pixel contribuente

Notazione per le coordinate:

obj_x  : coordinate intere di obj

obj_x_ : coordinate float  di obj

Notazione per le coppie di coordinate:

obj1_obj2_x  : coordinate intere di obj1 rispetto a obj2

obj1_obj2_x_ : coordinate float  di obj1 rispetto a obj2

Matrici affini:

aff_xy2uv  : matrice di trasformazione delle coordinate da XY a UV

aff0_xy2uv : stessa matrice con la parte di shift messa a 0

Una tantum:

aff_uv2xy  = aff

aff_xy2uv  = aff_inv (aff_uv2xy)

aff0_uv2xy = aff_place (00, 00, aff_uv2xy)

  vedi il codice, comunque ottimizzo una rotazione seguita da una scalatura

  anisotropa. Cerco i fattori di scala facendo radici di somme di quadrati.

  Mi regolo sui fattori di scala come se dovessi considerare la scalatura

  da sola. In questo modo tutto si riporta alla vecchia maniera se i fattori

  di scala sono uguali. Se sono diversi il risultato e' comunque esatto per

  rotazioni di multipli di 90 gradi, e non ha discontinuita'.

aff0_uv2st = aff_mult (aff0_xy2st, aff0_uv2xy)

aff0_st2fg = aff_scale (filter_resolution, Aff_I)

aff0_uv2fg = aff_mult (aff0_st2fg, aff0_uv2st)

pix_ref_uv[] = tutti quelli che servono (vedi sotto)

pix_ref_fg_  = AFF_M_V (aff0_uv2fg, pix_ref_uv)

pix_ref_fg[] = ROUND (pix_ref_fg_)

Ciclo su out_xy:

out_uv_ = AFF_M_V (aff_xy2uv, out_xy)

ref_uv = INT_LE (out_uv_)

ref_out_uv_ = ref_uv - out_uv_

ref_out_fg_ = AFF_M_V (aff0_uv2fg, ref_out_uv_)

ref_out_fg  = ROUND (ref_out_fg_)

Ciclo sui pix:

pix_out_fg = pix_ref_fg + ref_out_fg

weight = filter[pix_out_f] * filter[pix_out_g]

Per sapere quali sono i pix che servono:

-filter_fg_radius < pix_out_fg < filter_fg_radius

min_pix_out_uv_ < pix_out_uv_ < max_pix_out_uv_

min_pix_out_uv_ < pix_ref_uv_ + ref_out_uv_ < max_pix_out_uv_

min_pix_out_uv_ + out_ref_uv_ < pix_ref_uv_ < max_pix_out_uv_ + out_ref_uv_

min_pix_out_uv_ < pix_ref_uv_ < max_pix_out_uv_ + 1

Ciclo su tutti quelli che soddisfano questa condizione

0 <= out_ref_uv_ < 1

-1 < ref_out_uv_ <= 0

min_ref_out_fg_ <= ref_out_fg_ <= max_ref_out_fg_

min_ref_out_fg  <= ref_out_fg  <= max_ref_out_fg

-filter_fg_radius < pix_out_fg < filter_fg_radius

-filter_fg_radius < pix_ref_fg + ref_out_fg < filter_fg_radius

-filter_fg_radius - ref_out_fg < pix_ref_fg < filter_fg_radius - ref_out_fg

-filter_fg_radius - max_ref_out_fg < pix_ref_fg < filter_fg_radius - min_ref_out_fg

Scarto quelli che non soddisfano questa condizione

Come e' fatto il filtro:

  TOP                          filter_array[filter_array_size-1]

   |  filter[max_filter_fg]

   |  filter[max_pix_out_fg]

   |  filter[0]

   |  filter[min_pix_out_fg]

  BOT filter[min_filter_fg] == filter_array[0]

*/

//------------------------------------------------------------------------------

//---------------------------------------------------------------------------

#if !defined(TNZ_LITTLE_ENDIAN)

TNZ_LITTLE_ENDIAN undefined !!

#endif

	//2^36 * 1.5,  (52-_shiftamt=36) uses limited precision to floor

	const double _double2fixmagic = 68719476736.0 * 1.5;

//16.16 fixed point representation

const TINT32 _shiftamt = 16;

#if TNZ_LITTLE_ENDIAN

#define iexp_ 1

#define iman_ 0

#else

#define iexp_ 0

#define iman_ 1

#endif

inline TINT32 Double2Int(double val)

{

	val = val + _double2fixmagic;

	return ((TINT32 *)&val)[iman_] >> _shiftamt;

}

#define DOUBLE_TO_INT32(D) (d2iaux = D, d2iaux += _double2fixmagic, (((TINT32 *)&(d2iaux))[iman_] >> _shiftamt))

//#define USE_DOUBLE_TO_INT

//===========================================================================

inline double sinc0(double x, int a)

{

	return sin((M_PI / (a)) * (x)) / ((M_PI / (a)) * (x));

}

inline double sinc(double x, int a)

{

	return (x) == 0.0 ? 1.0 : sin((M_PI / (a)) * (x)) / ((M_PI / (a)) * (x));

}

inline UCHAR TO8BIT(float X)

{

	return (((X) < 0.0F) ? 0 : (((X) > 255.0F) ? 255 : tround(X)));

}

const UCHAR BORDER_GR8 = 255;

const UCHAR GREY_GR8 = 127;

#ifdef USE_INLINE_FUNS

//---------------------------------------------------------------------------

inline double aff0MV1(const TAffine &aff, double v1, double v2)

{

	return aff.a11 * v1 + aff.a12 * v2;

}

//---------------------------------------------------------------------------

inline double affMV1(const TAffine &aff, double v1, double v2)

{

	return aff.a11 * v1 + aff.a12 * v2 + aff.a13;

}

//---------------------------------------------------------------------------

inline double aff0MV2(const TAffine &aff, double v1, double v2)

{

	return aff.a21 * v1 + aff.a22 * v2;

}

//---------------------------------------------------------------------------

inline double affMV2(const TAffine &aff, double v1, double v2)

{

	return aff.a21 * v1 + aff.a22 * v2 + aff.a23;

}

#else // !USE_INLINE_FUNS

#ifndef USE_DOUBLE_TO_INT

#define ROUND(x) ((int)(((int)(-0.9F) == 0 && (x) < 0.0F) ? ((x)-0.5F) : ((x) + 0.5F)))

#define ROUNDP(x) ((int)((x) + 0.5F))

#define FLOOR(x) ((int)(x) > (x) ? (int)(x)-1 : (int)(x))

#define CEIL(x) ((int)(x) < (x) ? (int)(x) + 1 : (int)(x))

#else

#define ROUND(x) (DOUBLE_TO_INT32(((int)(-0.9F) == 0 && (x) < 0.0F) ? ((x)-0.5F) : ((x) + 0.5F)))

#define ROUNDP(x) (DOUBLE_TO_INT32((x) + 0.5F))

#define FLOOR(x) (DOUBLE_TO_INT32(x) > (x) ? DOUBLE_TO_INT32(x) - 1 : DOUBLE_TO_INT32(x))

#define CEIL(x) (DOUBLE_TO_INT32(x) < (x) ? DOUBLE_TO_INT32(x) + 1 : DOUBLE_TO_INT32(x))

#endif

#define INTLE(x) (FLOOR(x))

#define INTGT(x) (FLOOR(x) + 1)

#define INTLT(x) (CEIL(x) - 1)

#define INTGE(x) (CEIL(x))

#define NOT_LESS_THAN(MIN, X) \

	{                         \

		if ((X) < (MIN))      \

			(X) = (MIN);      \

	}

#define NOT_MORE_THAN(MAX, X) \

	{                         \

		if ((X) > (MAX))      \

			(X) = (MAX);      \

	}

#define tround ROUND

#define troundp ROUNDP

#define tfloor FLOOR

#define tceil CEIL

#define intLE INTLE

#define intGT INTGT

#define intLT INTLT

#define intGE INTGE

#define notLessThan NOT_LESS_THAN

#define notMoreThan NOT_MORE_THAN

#define AFF0_M_V_1(AFF, V1, V2) ((AFF).a11 * (V1) + (AFF).a12 * (V2))

#define AFF0_M_V_2(AFF, V1, V2) ((AFF).a21 * (V1) + (AFF).a22 * (V2))

#define AFF_M_V_1(AFF, V1, V2) ((AFF).a11 * (V1) + (AFF).a12 * (V2) + (AFF).a13)

#define AFF_M_V_2(AFF, V1, V2) ((AFF).a21 * (V1) + (AFF).a22 * (V2) + (AFF).a23)

#define aff0MV1 AFF0_M_V_1

#define aff0MV2 AFF0_M_V_2

#define affMV1 AFF_M_V_1

#define affMV2 AFF_M_V_2

#endif // USE_INLINE_FUNS

//---------------------------------------------------------------------------

struct FILTER {

	int first, last;

	float *w;

	float *w_base;

};

struct NOCALC {

	int first, last;

};

//---------------------------------------------------------------------------

inline int get_filter_radius(TRop::ResampleFilterType flt_type)

{

	switch (flt_type) {

	case TRop::Triangle:

		return 1;

	case TRop::Mitchell:

		return 2;

	case TRop::Cubic5:

		return 2;

	case TRop::Cubic75:

		return 2;

	case TRop::Cubic1:

		return 2;

	case TRop::Hann2:

		return 2;

	case TRop::Hann3:

		return 3;

	case TRop::Hamming2:

		return 2;

	case TRop::Hamming3:

		return 3;

	case TRop::Lanczos2:

		return 2;

	case TRop::Lanczos3:

		return 3;

	case TRop::Gauss:

		return 2;

	default:

		assert(!"bad filter type");

	}

	return 0;

}

//---------------------------------------------------------------------------

//!Equivalent to aff * TRectD(u0, v0, u1, v0).

inline void minmax(double u0, double v0,

				   double u1, double v1, const TAffine &aff,

				   double &x0, double &y0,

				   double &x1, double &y1)

{

	double xmin, ymin;

	double xmax, ymax;

	double x_a, y_a;

	double x_b, y_b;

	double x_c, y_c;

	double x_d, y_d;

	x_a = affMV1(aff, u0, v0);

	y_a = affMV2(aff, u0, v0);

	x_b = affMV1(aff, u1, v0);

	y_b = affMV2(aff, u1, v0);

	x_c = affMV1(aff, u1, v1);

	y_c = affMV2(aff, u1, v1);

	x_d = affMV1(aff, u0, v1);

	y_d = affMV2(aff, u0, v1);

	xmin = std::min(x_a, x_b);

	xmax = std::max(x_a, x_b);

	xmin = std::min(xmin, x_c);

	xmax = std::max(xmax, x_c);

	xmin = std::min(xmin, x_d);

	xmax = std::max(xmax, x_d);

	ymin = std::min(y_a, y_b);

	ymax = std::max(y_a, y_b);

	ymin = std::min(ymin, y_c);

	ymax = std::max(ymax, y_c);

	ymin = std::min(ymin, y_d);

	ymax = std::max(ymax, y_d);

	x0 = xmin;

	y0 = ymin;

	x1 = xmax;

	y1 = ymax;

}

/*---------------------------------------------------------------------------*/

/*

inline bool trivial_rot (TAffine inv, int *dudx, int *dudy,

                                    int *dvdx, int *dvdy)

{

*dudx = 0;

*dudy = 0;

*dvdx = 0;

*dvdy = 0;

if (! (inv.a12 == 0 && inv.a21 == 0 || inv.a11 == 0 && inv.a22 == 0))

  return false;

if (! (inv.a11 == 1 || inv.a11 == 0 || inv.a11 == -1))

  return false;

if (! (inv.a12 == 1 || inv.a12 == 0 || inv.a12 == -1))

  return false;

if (! (inv.a21 == 1 || inv.a21 == 0 || inv.a21 == -1))

  return false;

if (! (inv.a22 == 1 || inv.a22 == 0 || inv.a22 == -1))

  return false;

*dudx = (int)inv.a11;

*dudy = (int)inv.a12;

*dvdx = (int)inv.a21;

*dvdy = (int)inv.a22;

return true;

}

*/

//-----------------------------------------------------------------------------

//

// see Mitchell&Netravali, "Reconstruction Filters in Computer Graphics",

// SIGGRAPH 88.  Mitchell code provided by Paul Heckbert.

//

//

//-----------------------------------------------------------------------------

static double p0, p2, p3, q0, q1, q2, q3;

inline void mitchellinit(double b, double c)

{

	p0 = (6.0 - 2.0 * b) / 6.0;

	p2 = (-18.0 + 12.0 * b + 6.0 * c) / 6.0;

	p3 = (12.0 - 9.0 * b - 6.0 * c) / 6.0;

	q0 = (8.0 * b + 24.0 * c) / 6.0;

	q1 = (-12.0 * b - 48.0 * c) / 6.0;

	q2 = (6.0 * b + 30.0 * c) / 6.0;

	q3 = (-b - 6.0 * c) / 6.0;

}

const int fltradMitchell = 2;

static inline double flt_mitchell(double x) /*Mitchell & Netravali's two-param cubic*/

{

	static int mitfirsted;

	if (!mitfirsted) {

		mitchellinit(1.0 / 3.0, 1.0 / 3.0);

		mitfirsted = 1;

	}

	if (x < -2.0)

		return 0.0;

	if (x < -1.0)

		return (q0 - x * (q1 - x * (q2 - x * q3)));

	if (x < 0.0)

		return (p0 + x * x * (p2 - x * p3));

	if (x < 1.0)

		return (p0 + x * x * (p2 + x * p3));

	if (x < 2.0)

		return (q0 + x * (q1 + x * (q2 + x * q3)));

	return 0.0;

}

//-----------------------------------------------------------------------------

const int fltradTriangle = 1;

static inline double flt_triangle(double x)

{

	if (x < -1.0)

		return 0.0;

	if (x < 0.0)

		return 1.0 + x;

	if (x < 1.0)

		return 1.0 - x;

	return 0.0;

}

//-----------------------------------------------------------------------------

const int fltradCubic5 = 2;

static inline double flt_cubic_5(double x)

{

	if (x < 0.0)

		x = -x;

	if (x < 1.0)

		return 2.5 * x * x * x - 3.5 * x * x + 1;

	if (x < 2.0)

		return 0.5 * x * x * x - 2.5 * x * x + 4 * x - 2;

	return 0.0;

}

//-----------------------------------------------------------------------------

const int fltradCubic75 = 2;

static inline double flt_cubic_75(double x)

{

	if (x < 0.0)

		x = -x;

	if (x < 1.0)

		return 2.75 * x * x * x - 3.75 * x * x + 1;

	if (x < 2.0)

		return 0.75 * x * x * x - 3.75 * x * x + 6 * x - 3;

	return 0.0;

}

//-----------------------------------------------------------------------------

const int fltradCubic1 = 2;

static inline double flt_cubic_1(double x)

{

	if (x < 0.0)

		x = -x;

	if (x < 1.0)

		return 3 * x * x * x - 4 * x * x + 1;

	if (x < 2.0)

		return x * x * x - 5 * x * x + 8 * x - 4;

	return 0.0;

}

//-----------------------------------------------------------------------------

const int fltradHann2 = 2;

static inline double flt_hann2(double x)

{

	if (x <= -2.0)

		return 0.0;

	if (x < 2.0)

		return sinc(x, 1) * (0.5 + 0.5 * cos(M_PI_2 * x));

	return 0.0;

}

//-----------------------------------------------------------------------------

const int fltradHann3 = 3;

static inline double flt_hann3(double x)

{

	if (x <= -3.0)

		return 0.0;

	if (x < 3.0)

		return sinc(x, 1) * (0.5 + 0.5 * cos(M_PI_3 * x));

	return 0.0;

}

//-----------------------------------------------------------------------------

const int fltradHamming2 = 2;

static inline double flt_hamming2(double x)

{

	if (x <= -2.0)

		return 0.0;

	if (x < 2.0)

		return sinc(x, 1) * (0.54 + 0.46 * cos(M_PI_2 * x));

	return 0.0;

}

//-----------------------------------------------------------------------------

const int fltradHamming3 = 3;

static inline double flt_hamming3(double x)

{

	if (x <= -3.0)

		return 0.0;

	if (x < 3.0)

		return sinc(x, 1) * (0.54 + 0.46 * cos(M_PI_3 * x));

	return 0.0;

}

//-----------------------------------------------------------------------------

const int fltradLanczos2 = 2;

static inline double flt_lanczos2(double x)

{

	if (x <= -2.0)

		return 0.0;

	if (x < 2.0)

		return sinc(x, 1) * sinc(x, 2);

	return 0.0;

}

//-----------------------------------------------------------------------------

const int fltradLanczos3 = 3;

static inline double flt_lanczos3(double x)

{

	if (x <= -3.0)

		return 0.0;

	if (x < 3.0)

		return sinc(x, 1) * sinc(x, 3);

	return 0.0;

}

//-----------------------------------------------------------------------------

const int fltradGauss = 2;

static inline double flt_gauss(double x)

{

	if (x <= -2.0)

		return 0.0;

	if (x < 2.0)

		return exp(-M_PI * x * x);

	return 0.0; /* exp(-M_PI*2*2)~=3.5*10^-6 */

}

//-----------------------------------------------------------------------------

const int fltradW1 = 2;

static inline double flt_w_1(double x)

{

	if (x < 0.0)

		x = -x;

	if (x < 0.5)

		return 1 - 0.5 * x;

	if (x < 1.0)

		return 1.5 - 1.5 * x;

	if (x < 1.5)

		return 0.5 - 0.5 * x;

	if (x < 2.0)

		return 0.5 * x - 1.0;

	return 0.0;

}

//-----------------------------------------------------------------------------

static inline void get_flt_fun_rad(TRop::ResampleFilterType flt_type,

								   double (**flt_fun)(double), double &flt_rad)

{

	double (*fun)(double);

	double rad;

	switch (flt_type) {

	case TRop::Triangle: fun = flt_triangle; rad = fltradTriangle; break;

	case TRop::Mitchell: fun = flt_mitchell; rad = fltradMitchell; break;

	case TRop::Cubic5  : fun = flt_cubic_5 ; rad = fltradCubic5  ; break;

	case TRop::Cubic75 : fun = flt_cubic_75; rad = fltradCubic75 ; break;

	case TRop::Cubic1  : fun = flt_cubic_1 ; rad = fltradCubic1  ; break;

	case TRop::Hann2   : fun = flt_hann2   ; rad = fltradHann2   ; break;

	case TRop::Hann3   : fun = flt_hann3   ; rad = fltradHann3   ; break;

	case TRop::Hamming2: fun = flt_hamming2; rad = fltradHamming2; break;

	case TRop::Hamming3: fun = flt_hamming3; rad = fltradHamming3; break;

	case TRop::Lanczos2: fun = flt_lanczos2; rad = fltradLanczos2; break;

	case TRop::Lanczos3: fun = flt_lanczos3; rad = fltradLanczos3; break;

	case TRop::Gauss   : fun = flt_gauss   ; rad = fltradGauss   ; break;

	case 101           : fun = flt_w_1     ; rad = fltradW1      ; break;

	default            : fun = flt_triangle; rad = fltradTriangle; break;

	}

	if (flt_fun)

		*flt_fun = fun;

	flt_rad = rad;

}

//---------------------------------------------------------------------------

static FILTER *create_filter(TRop::ResampleFilterType flt_type, double blur,

							 double dx_du, double delta_x, int lx,

							 double &xrad, int &umin, int &umax, int &uwidth)

{

	double (*flt_fun)(double);

	FILTER *filter, *f;

	double du_dx;

	int x;

	double u_;

	int u, ulo, uhi, ulomin, uhimax, m, n, nmax;

	double flt_rad, rad_u, rad_x, nodedist_u, nodefreq_u, sum, norm, w;

#ifdef USE_DOUBLE_TO_INT

	double d2iaux;

#endif

	get_flt_fun_rad(flt_type, &flt_fun, flt_rad);

	du_dx = 1 / dx_du;

	if (dx_du > 1)

		nodedist_u = blur; /* magnification */

	else

		nodedist_u = du_dx * blur; /* minification */

	rad_u = flt_rad * nodedist_u;

	rad_x = rad_u * dx_du;

	nodefreq_u = 1 / nodedist_u;

	/*

mu = lu - 1;

*/

	filter = new FILTER[lx];

	nmax = 0;

	ulomin = c_maxint - 1;

	uhimax = c_minint + 1;

	for (x = 0; x < lx; x++) {

		f = filter + x;

		u_ = (x - delta_x) * du_dx;

		ulo = intGT(u_ - rad_u);

		uhi = intLT(u_ + rad_u);

		/*

  NOT_LESS_THAN( 0, ulo)

  NOT_MORE_THAN(mu, uhi)

*/

		m = uhi - ulo + 1;

		if (m > 0) {

			f->w_base = new float[m];

			f->w = f->w_base - ulo;

			for (sum = 0.0, u = ulo; u <= uhi; u++) {

				w = (*flt_fun)((u - u_) * nodefreq_u);

				sum += w;

				f->w[u] = (float)w;

			}

			for (; ulo <= uhi; ulo++)

				if (f->w[ulo])

					break;

			for (; uhi >= ulo; uhi--)

				if (f->w[uhi])

					break;

			if (ulo < ulomin)

				ulomin = ulo;

			if (uhi > uhimax)

				uhimax = uhi;

			n = uhi - ulo + 1;

			if (n > nmax)

				nmax = n;

			f->first = ulo;

			f->last = uhi;

			norm = 1 / sum;

			for (u = ulo; u <= uhi; u++)

				f->w[u] *= (float)norm;

		} else {

			f->w_base = 0;

			f->first = ulo;

			f->last = uhi;

		}

	}

	xrad = rad_x;

	umin = ulomin;

	umax = uhimax;

	uwidth = nmax;

	return filter;

}

//-----------------------------------------------------------------------------

static NOCALC *create_nocalc(TRop::ResampleFilterType flt_type, double blur,

							 double dx_du, double delta_x, int lx,

							 int umin, int umax, int &xwidth)

{

	/*

Il nocalc serve a stabilire che un insieme di pixel u (di ingresso)

non ha bisogno di essere calcolato, perche tutti i pixel x (di uscita)

su cui questo insieme si distribuisce non hanno bisogno di essere calcolati.

Il significato del nocalc a coordinata x e':

se arrivati a x si e' trovata una sequenza di width pixel x che

non e' necessario calcolare, allora non e' necessario calcolare tutti i

pixel u da nocalc->first a nocalc->last.

Per primo va calcolata la width. Deve essere garantito che tutti i pixel u

siano coperti dal vettore di nocalc.

Un pixel u viene usato da un intervallo di x largo quanto il filtro (in x),

cioe' un intervallo aperto (-radx_,radx_) intorno a x(u).

Aggiungendo una unita' x a questa larghezza si ha una larghezza tale che

se tutti i pixel x sono nocalc, un intervallo largo 1 in unita' x

di pixel u non necessita di essere calcolato.

Vogliamo che ulo_ <= first <= last < uhi_ con uhi_ = ulo_ + u(1).

Devono essere nocalc almeno gli x in (x(ulo_)-radx_, x(uhi_)+radx_).

Poniamo x = x(uhi_)+radx_-1.

uhi_ = u(x-radx_+1)

ulo_ = u(x-radx_)

x(ulo_)-radx_ = x-2*radx_   ma questo punto e' escluso, quindi l'intero GT e'

x - width + 1 = INT_LE (x-2*radx_+1)

1 - INT_LE (-2*radx_+1) = width

1 + INT_GE (2*radx_-1) = width

INT_GE (2*radx_) = width

Pero' per sicurezza facciamo

INT_GT (2*radx_) = width

*/

	NOCALC *nocalc;

	int width;

	double flt_rad;

	double rad_x;

	double du_dx;

	double ulo_, uhi_;

	int ulo, uhi;

	int x;

#ifdef USE_DOUBLE_TO_INT

	double d2iaux;

#endif

	du_dx = 1 / dx_du;

	get_flt_fun_rad(flt_type, 0, flt_rad);

	if (dx_du > 1) /* sto ingrandendo */

		rad_x = flt_rad * blur * dx_du;

	else

		rad_x = flt_rad * blur;

	rad_x += 0.5; /* ?!?!?!?!? */

	width = intGT(2 * rad_x + 1);

	nocalc = new NOCALC[lx + width - 1];

	for (x = 0; x < lx + width - 1; x++) {

		ulo_ = (x - rad_x - delta_x) * du_dx;