1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
|
#ifdef FUN
/* this code is not compiled as such, but rather included in other code */
/*TABLE HYPER PRECISION
adds TYH bits of precision to table indexes without needing a
table that big, and without overusing 'double' variables;
does not add any visible difference, but has a small computational cost
*/
#if (SINC_TABLE_UNIT*KERNEL_WIDTH) < (1<<11)
#define TYH 4
#define SINC_TABLE(I) sinc_table[ (I) >> TYH]
#elif (SINC_TABLE_UNIT*KERNEL_WIDTH) < (1<<13)
#define TYH 2
#define SINC_TABLE(I) sinc_table[ (I) >> TYH]
#else
#define TYH 0
#define SINC_TABLE(I) sinc_table[I]
#endif
#if defined(KERNEL_sinc_fast) && ( SINC_TABLE_UNIT < (PIXEL_MAX-PIXEL_MIN) )
#warning the kernel is based on the sinc(); the function sinc is replaced by a table; but the table is not precise enough for your type of pixels
/* in this case you may want to interpolate in the table
#define MASK ((1<<TYH)-1)
#define SINC_TABLE(I) (sinc_table[ (I) >> TYH] * ((I)&MASK) + \
sinc_table[ 1+( (I) >> TYH)]* ((1<<TYH)- ((I)&MASK)) )
but it is unclear if this is really any faster than using the sinc()
function with doubles
*/
#endif
static void FUN(KERNEL) /*that is, resample_array_inv_conv_ ## KERNEL */
(const double *F,
const PIXEL_TYPE *src, int s_len, int s_stride,
PIXEL_TYPE *dst, int d_len, int d_stride)
{
int i,p=0,j,pj, firstj,lastj;
double x,px,c,s,v,dx;
#ifdef KERNEL_sinc_fast
int index, indexbound, increment;
#endif
px=F[1];
for(i=0;i<d_len;i++) {
v=0; c=0;
x=F[i];
/* FIXED! READ OUT OF ARRAY!!!
//nx=F[i+1]; */
dx=ABS(x-px);
if(dx<1) dx=1;
firstj=floor(x-dx*KERNEL_WIDTH);
lastj=ceil(x+dx*KERNEL_WIDTH);
j=firstj;
/* note that in the past I had reversed the usage of the filter;
this does not really make any difference as long as the filter is
symmetric... but I want to be more precise */
#ifdef KERNEL_sinc_fast
index=(( (((double)firstj-x) / dx) )* (double)(SINC_TABLE_UNIT << TYH));
increment= ( (SINC_TABLE_UNIT << TYH) / dx );
indexbound=( (KERNEL_WIDTH*SINC_TABLE_UNIT) << TYH);
/*roughly
index ~ -indexbound;
but if it is approximated to that, then the warping loses
sub-pixel precision, and there are visible artifacts
*/
/* this condition is always true
//if (index >= -indexbound) */
if(increment<1) {
printf("%s:%s: precision underflow! there is a jump of %f pixels in the mesh; this is too big w.r.t. the size %d of the kernel\n",
__FILE__,__FUNCTION__,dx, ((SINC_TABLE_UNIT << TYH)));
increment=1;
}
{ j++; index += increment; }
#endif
#ifndef KERNEL_sinc_fast
/* this is the convolution with a generic kernel */
for ( ; j<=lastj ; j++ ) {
s=KERNEL( ((double)j-x) / dx );
pj=j;
if ( pj<0 ) pj=0;
else if( pj>=s_len) pj=s_len-1;
v+=s* (double)src[pj* s_stride];
c+=s;
}
#else
/* this is how it is computed very fast for a kernel based on sinc, such as lanczos */
for ( ; index<=0 ; ) {
s=SINC_TABLE(-index );
pj=j;
if ( pj<0 ) pj=0;
else if( pj>=s_len) pj=s_len-1;
v+=s* (double)src[pj*s_stride];
c+=s;
j++;
index+=increment;
}
for ( ; index <= indexbound ; ) {
s=SINC_TABLE(index);
pj=j;
if ( pj<0 ) pj=0;
else if( pj>=s_len) pj=s_len-1;
v+=s* (double)src[pj*s_stride];
c+=s;
j++;
index+=increment;
}
#endif
if(c<0.0001 && c > -0.0001)
fprintf(stderr,"%s:%s:%d: too small c=%f\n",__FILE__,__FUNCTION__,__LINE__,c);
else
v/=c;
dst[p]=CLAMP(v,PIXEL_MIN,PIXEL_MAX);
p+= d_stride;
px=x;
}
}
#endif
|
