1 // This code is in the public domain -- castanyo@yahoo.es
2
3 /** @file Filter.cpp
4 * @brief Image filters.
5 *
6 * Jonathan Blow articles:
7 * http://number-none.com/product/Mipmapping, Part 1/index.html
8 * http://number-none.com/product/Mipmapping, Part 2/index.html
9 *
10 * References from Thacher Ulrich:
11 * See _Graphics Gems III_ "General Filtered Image Rescaling", Dale A. Schumacher
12 * http://tog.acm.org/GraphicsGems/gemsiii/filter.c
13 *
14 * References from Paul Heckbert:
15 * A.V. Oppenheim, R.W. Schafer, Digital Signal Processing, Prentice-Hall, 1975
16 *
17 * R.W. Hamming, Digital Filters, Prentice-Hall, Englewood Cliffs, NJ, 1983
18 *
19 * W.K. Pratt, Digital Image Processing, John Wiley and Sons, 1978
20 *
21 * H.S. Hou, H.C. Andrews, "Cubic Splines for Image Interpolation and
22 * Digital Filtering", IEEE Trans. Acoustics, Speech, and Signal Proc.,
23 * vol. ASSP-26, no. 6, Dec. 1978, pp. 508-517
24 *
25 * Paul Heckbert's zoom library.
26 * http://www.xmission.com/~legalize/zoom.html
27 *
28 * Reconstruction Filters in Computer Graphics
29 * http://www.mentallandscape.com/Papers_siggraph88.pdf
30 *
31 * More references:
32 * http://www.worldserver.com/turk/computergraphics/ResamplingFilters.pdf
33 * http://www.dspguide.com/ch16.htm
34 */
35
36 #include "Filter.h"
37
38 #include <nvmath/Vector.h> // Vector4
39 #include <nvcore/Containers.h> // swap
40
41 using namespace nv;
42
43 namespace
44 {
45 // Sinc function.
sincf(const float x)46 inline static float sincf(const float x)
47 {
48 if (fabs(x) < NV_EPSILON) {
49 //return 1.0;
50 return 1.0f + x*x*(-1.0f/6.0f + x*x*1.0f/120.0f);
51 }
52 else {
53 return sin(x) / x;
54 }
55 }
56
57 // Bessel function of the first kind from Jon Blow's article.
58 // http://mathworld.wolfram.com/BesselFunctionoftheFirstKind.html
59 // http://en.wikipedia.org/wiki/Bessel_function
bessel0(float x)60 inline static float bessel0(float x)
61 {
62 const float EPSILON_RATIO = 1e-6f;
63 float xh, sum, pow, ds;
64 int k;
65
66 xh = 0.5f * x;
67 sum = 1.0f;
68 pow = 1.0f;
69 k = 0;
70 ds = 1.0;
71 while (ds > sum * EPSILON_RATIO) {
72 ++k;
73 pow = pow * (xh / k);
74 ds = pow * pow;
75 sum = sum + ds;
76 }
77
78 return sum;
79 }
80
81 /*// Alternative bessel function from Paul Heckbert.
82 static float _bessel0(float x)
83 {
84 const float EPSILON_RATIO = 1E-6;
85 float sum = 1.0f;
86 float y = x * x / 4.0f;
87 float t = y;
88 for(int i = 2; t > EPSILON_RATIO; i++) {
89 sum += t;
90 t *= y / float(i * i);
91 }
92 return sum;
93 }*/
94
95 } // namespace
96
97
Filter(float width)98 Filter::Filter(float width) : m_width(width)
99 {
100 }
101
~Filter()102 /*virtual*/ Filter::~Filter()
103 {
104 }
105
sampleDelta(float x,float scale) const106 float Filter::sampleDelta(float x, float scale) const
107 {
108 return evaluate((x + 0.5f)* scale);
109 }
110
sampleBox(float x,float scale,int samples) const111 float Filter::sampleBox(float x, float scale, int samples) const
112 {
113 float sum = 0;
114 float isamples = 1.0f / float(samples);
115
116 for(int s = 0; s < samples; s++)
117 {
118 float p = (x + (float(s) + 0.5f) * isamples) * scale;
119 float value = evaluate(p);
120 sum += value;
121 }
122
123 return sum * isamples;
124 }
125
sampleTriangle(float x,float scale,int samples) const126 float Filter::sampleTriangle(float x, float scale, int samples) const
127 {
128 float sum = 0;
129 float isamples = 1.0f / float(samples);
130
131 for(int s = 0; s < samples; s++)
132 {
133 float offset = (2 * float(s) + 1.0f) * isamples;
134 float p = (x + offset - 0.5f) * scale;
135 float value = evaluate(p);
136
137 float weight = offset;
138 if (weight > 1.0f) weight = 2.0f - weight;
139
140 sum += value * weight;
141 }
142
143 return 2 * sum * isamples;
144 }
145
146
147
148
149
BoxFilter()150 BoxFilter::BoxFilter() : Filter(0.5f) {}
BoxFilter(float width)151 BoxFilter::BoxFilter(float width) : Filter(width) {}
152
evaluate(float x) const153 float BoxFilter::evaluate(float x) const
154 {
155 if (fabs(x) <= m_width) return 1.0f;
156 else return 0.0f;
157 }
158
159
TriangleFilter()160 TriangleFilter::TriangleFilter() : Filter(1.0f) {}
TriangleFilter(float width)161 TriangleFilter::TriangleFilter(float width) : Filter(width) {}
162
evaluate(float x) const163 float TriangleFilter::evaluate(float x) const
164 {
165 x = fabs(x);
166 if( x < m_width ) return m_width - x;
167 return 0.0f;
168 }
169
170
QuadraticFilter()171 QuadraticFilter::QuadraticFilter() : Filter(1.5f) {}
172
evaluate(float x) const173 float QuadraticFilter::evaluate(float x) const
174 {
175 x = fabs(x);
176 if( x < 0.5f ) return 0.75f - x * x;
177 if( x < 1.5f ) {
178 float t = x - 1.5f;
179 return 0.5f * t * t;
180 }
181 return 0.0f;
182 }
183
184
CubicFilter()185 CubicFilter::CubicFilter() : Filter(1.0f) {}
186
evaluate(float x) const187 float CubicFilter::evaluate(float x) const
188 {
189 // f(t) = 2|t|^3 - 3|t|^2 + 1, -1 <= t <= 1
190 x = fabs(x);
191 if( x < 1.0f ) return((2.0f * x - 3.0f) * x * x + 1.0f);
192 return 0.0f;
193 }
194
195
BSplineFilter()196 BSplineFilter::BSplineFilter() : Filter(2.0f) {}
197
evaluate(float x) const198 float BSplineFilter::evaluate(float x) const
199 {
200 x = fabs(x);
201 if( x < 1.0f ) return (4.0f + x * x * (-6.0f + x * 3.0f)) / 6.0f;
202 if( x < 2.0f ) {
203 float t = 2.0f - x;
204 return t * t * t / 6.0f;
205 }
206 return 0.0f;
207 }
208
209
MitchellFilter()210 MitchellFilter::MitchellFilter() : Filter(2.0f) { setParameters(1.0f/3.0f, 1.0f/3.0f); }
211
evaluate(float x) const212 float MitchellFilter::evaluate(float x) const
213 {
214 x = fabs(x);
215 if( x < 1.0f ) return p0 + x * x * (p2 + x * p3);
216 if( x < 2.0f ) return q0 + x * (q1 + x * (q2 + x * q3));
217 return 0.0f;
218 }
219
setParameters(float b,float c)220 void MitchellFilter::setParameters(float b, float c)
221 {
222 p0 = (6.0f - 2.0f * b) / 6.0f;
223 p2 = (-18.0f + 12.0f * b + 6.0f * c) / 6.0f;
224 p3 = (12.0f - 9.0f * b - 6.0f * c) / 6.0f;
225 q0 = (8.0f * b + 24.0f * c) / 6.0f;
226 q1 = (-12.0f * b - 48.0f * c) / 6.0f;
227 q2 = (6.0f * b + 30.0f * c) / 6.0f;
228 q3 = (-b - 6.0f * c) / 6.0f;
229 }
230
231
LanczosFilter()232 LanczosFilter::LanczosFilter() : Filter(3.0f) {}
233
evaluate(float x) const234 float LanczosFilter::evaluate(float x) const
235 {
236 x = fabs(x);
237 if( x < 3.0f ) return sincf(PI * x) * sincf(PI * x / 3.0f);
238 return 0.0f;
239 }
240
241
SincFilter(float w)242 SincFilter::SincFilter(float w) : Filter(w) {}
243
evaluate(float x) const244 float SincFilter::evaluate(float x) const
245 {
246 return sincf(PI * x);
247 }
248
249
KaiserFilter(float w)250 KaiserFilter::KaiserFilter(float w) : Filter(w) { setParameters(4.0f, 1.0f); }
251
evaluate(float x) const252 float KaiserFilter::evaluate(float x) const
253 {
254 const float sinc_value = sincf(PI * x * stretch);
255 const float t = x / m_width;
256 if ((1 - t * t) >= 0) return sinc_value * bessel0(alpha * sqrtf(1 - t * t)) / bessel0(alpha);
257 else return 0;
258 }
259
setParameters(float alpha,float stretch)260 void KaiserFilter::setParameters(float alpha, float stretch)
261 {
262 this->alpha = alpha;
263 this->stretch = stretch;
264 }
265
266
267
268 /// Ctor.
Kernel1(const Filter & f,int iscale,int samples)269 Kernel1::Kernel1(const Filter & f, int iscale, int samples/*= 32*/)
270 {
271 nvDebugCheck(iscale > 1);
272 nvDebugCheck(samples > 0);
273
274 const float scale = 1.0f / iscale;
275
276 m_width = f.width() * iscale;
277 m_windowSize = (int)ceilf(2 * m_width);
278 m_data = new float[m_windowSize];
279
280 const float offset = float(m_windowSize) / 2;
281
282 float total = 0.0f;
283 for (int i = 0; i < m_windowSize; i++)
284 {
285 const float sample = f.sampleBox(i - offset, scale, samples);
286 m_data[i] = sample;
287 total += sample;
288 }
289
290 const float inv = 1.0f / total;
291 for (int i = 0; i < m_windowSize; i++)
292 {
293 m_data[i] *= inv;
294 }
295 }
296
297 /// Dtor.
~Kernel1()298 Kernel1::~Kernel1()
299 {
300 delete m_data;
301 }
302
303 /// Print the kernel for debugging purposes.
debugPrint()304 void Kernel1::debugPrint()
305 {
306 for (int i = 0; i < m_windowSize; i++) {
307 nvDebug("%d: %f\n", i, m_data[i]);
308 }
309 }
310
311
312
313 /// Ctor.
Kernel2(uint ws)314 Kernel2::Kernel2(uint ws) : m_windowSize(ws)
315 {
316 m_data = new float[m_windowSize * m_windowSize];
317 }
318
319 /// Copy ctor.
Kernel2(const Kernel2 & k)320 Kernel2::Kernel2(const Kernel2 & k) : m_windowSize(k.m_windowSize)
321 {
322 m_data = new float[m_windowSize * m_windowSize];
323 for (uint i = 0; i < m_windowSize * m_windowSize; i++) {
324 m_data[i] = k.m_data[i];
325 }
326 }
327
328
329 /// Dtor.
~Kernel2()330 Kernel2::~Kernel2()
331 {
332 delete m_data;
333 }
334
335 /// Normalize the filter.
normalize()336 void Kernel2::normalize()
337 {
338 float total = 0.0f;
339 for(uint i = 0; i < m_windowSize*m_windowSize; i++) {
340 total += fabs(m_data[i]);
341 }
342
343 float inv = 1.0f / total;
344 for(uint i = 0; i < m_windowSize*m_windowSize; i++) {
345 m_data[i] *= inv;
346 }
347 }
348
349 /// Transpose the kernel.
transpose()350 void Kernel2::transpose()
351 {
352 for(uint i = 0; i < m_windowSize; i++) {
353 for(uint j = i+1; j < m_windowSize; j++) {
354 swap(m_data[i*m_windowSize + j], m_data[j*m_windowSize + i]);
355 }
356 }
357 }
358
359 /// Init laplacian filter, usually used for sharpening.
initLaplacian()360 void Kernel2::initLaplacian()
361 {
362 nvDebugCheck(m_windowSize == 3);
363 // m_data[0] = -1; m_data[1] = -1; m_data[2] = -1;
364 // m_data[3] = -1; m_data[4] = +8; m_data[5] = -1;
365 // m_data[6] = -1; m_data[7] = -1; m_data[8] = -1;
366
367 m_data[0] = +0; m_data[1] = -1; m_data[2] = +0;
368 m_data[3] = -1; m_data[4] = +4; m_data[5] = -1;
369 m_data[6] = +0; m_data[7] = -1; m_data[8] = +0;
370
371 // m_data[0] = +1; m_data[1] = -2; m_data[2] = +1;
372 // m_data[3] = -2; m_data[4] = +4; m_data[5] = -2;
373 // m_data[6] = +1; m_data[7] = -2; m_data[8] = +1;
374 }
375
376
377 /// Init simple edge detection filter.
initEdgeDetection()378 void Kernel2::initEdgeDetection()
379 {
380 nvCheck(m_windowSize == 3);
381 m_data[0] = 0; m_data[1] = 0; m_data[2] = 0;
382 m_data[3] =-1; m_data[4] = 0; m_data[5] = 1;
383 m_data[6] = 0; m_data[7] = 0; m_data[8] = 0;
384 }
385
386 /// Init sobel filter.
initSobel()387 void Kernel2::initSobel()
388 {
389 if (m_windowSize == 3)
390 {
391 m_data[0] = -1; m_data[1] = 0; m_data[2] = 1;
392 m_data[3] = -2; m_data[4] = 0; m_data[5] = 2;
393 m_data[6] = -1; m_data[7] = 0; m_data[8] = 1;
394 }
395 else if (m_windowSize == 5)
396 {
397 float elements[] = {
398 -1, -2, 0, 2, 1,
399 -2, -3, 0, 3, 2,
400 -3, -4, 0, 4, 3,
401 -2, -3, 0, 3, 2,
402 -1, -2, 0, 2, 1
403 };
404
405 for (int i = 0; i < 5*5; i++) {
406 m_data[i] = elements[i];
407 }
408 }
409 else if (m_windowSize == 7)
410 {
411 float elements[] = {
412 -1, -2, -3, 0, 3, 2, 1,
413 -2, -3, -4, 0, 4, 3, 2,
414 -3, -4, -5, 0, 5, 4, 3,
415 -4, -5, -6, 0, 6, 5, 4,
416 -3, -4, -5, 0, 5, 4, 3,
417 -2, -3, -4, 0, 4, 3, 2,
418 -1, -2, -3, 0, 3, 2, 1
419 };
420
421 for (int i = 0; i < 7*7; i++) {
422 m_data[i] = elements[i];
423 }
424 }
425 else if (m_windowSize == 9)
426 {
427 float elements[] = {
428 -1, -2, -3, -4, 0, 4, 3, 2, 1,
429 -2, -3, -4, -5, 0, 5, 4, 3, 2,
430 -3, -4, -5, -6, 0, 6, 5, 4, 3,
431 -4, -5, -6, -7, 0, 7, 6, 5, 4,
432 -5, -6, -7, -8, 0, 8, 7, 6, 5,
433 -4, -5, -6, -7, 0, 7, 6, 5, 4,
434 -3, -4, -5, -6, 0, 6, 5, 4, 3,
435 -2, -3, -4, -5, 0, 5, 4, 3, 2,
436 -1, -2, -3, -4, 0, 4, 3, 2, 1
437 };
438
439 for (int i = 0; i < 9*9; i++) {
440 m_data[i] = elements[i];
441 }
442 }
443 }
444
445 /// Init prewitt filter.
initPrewitt()446 void Kernel2::initPrewitt()
447 {
448 if (m_windowSize == 3)
449 {
450 m_data[0] = -1; m_data[1] = 0; m_data[2] = -1;
451 m_data[3] = -1; m_data[4] = 0; m_data[5] = -1;
452 m_data[6] = -1; m_data[7] = 0; m_data[8] = -1;
453 }
454 else if (m_windowSize == 5)
455 {
456 // @@ Is this correct?
457 float elements[] = {
458 -2, -1, 0, 1, 2,
459 -2, -1, 0, 1, 2,
460 -2, -1, 0, 1, 2,
461 -2, -1, 0, 1, 2,
462 -2, -1, 0, 1, 2
463 };
464
465 for (int i = 0; i < 5*5; i++) {
466 m_data[i] = elements[i];
467 }
468 }
469 }
470
471 /// Init blended sobel filter.
initBlendedSobel(const Vector4 & scale)472 void Kernel2::initBlendedSobel(const Vector4 & scale)
473 {
474 nvCheck(m_windowSize == 9);
475
476 {
477 const float elements[] = {
478 -1, -2, -3, -4, 0, 4, 3, 2, 1,
479 -2, -3, -4, -5, 0, 5, 4, 3, 2,
480 -3, -4, -5, -6, 0, 6, 5, 4, 3,
481 -4, -5, -6, -7, 0, 7, 6, 5, 4,
482 -5, -6, -7, -8, 0, 8, 7, 6, 5,
483 -4, -5, -6, -7, 0, 7, 6, 5, 4,
484 -3, -4, -5, -6, 0, 6, 5, 4, 3,
485 -2, -3, -4, -5, 0, 5, 4, 3, 2,
486 -1, -2, -3, -4, 0, 4, 3, 2, 1
487 };
488
489 for (int i = 0; i < 9*9; i++) {
490 m_data[i] = elements[i] * scale.w();
491 }
492 }
493 {
494 const float elements[] = {
495 -1, -2, -3, 0, 3, 2, 1,
496 -2, -3, -4, 0, 4, 3, 2,
497 -3, -4, -5, 0, 5, 4, 3,
498 -4, -5, -6, 0, 6, 5, 4,
499 -3, -4, -5, 0, 5, 4, 3,
500 -2, -3, -4, 0, 4, 3, 2,
501 -1, -2, -3, 0, 3, 2, 1,
502 };
503
504 for (int i = 0; i < 7; i++) {
505 for (int e = 0; e < 7; e++) {
506 m_data[(i + 1) * 9 + e + 1] += elements[i * 7 + e] * scale.z();
507 }
508 }
509 }
510 {
511 const float elements[] = {
512 -1, -2, 0, 2, 1,
513 -2, -3, 0, 3, 2,
514 -3, -4, 0, 4, 3,
515 -2, -3, 0, 3, 2,
516 -1, -2, 0, 2, 1
517 };
518
519 for (int i = 0; i < 5; i++) {
520 for (int e = 0; e < 5; e++) {
521 m_data[(i + 2) * 9 + e + 2] += elements[i * 5 + e] * scale.y();
522 }
523 }
524 }
525 {
526 const float elements[] = {
527 -1, 0, 1,
528 -2, 0, 2,
529 -1, 0, 1,
530 };
531
532 for (int i = 0; i < 3; i++) {
533 for (int e = 0; e < 3; e++) {
534 m_data[(i + 3) * 9 + e + 3] += elements[i * 3 + e] * scale.x();
535 }
536 }
537 }
538 }
539
540
PolyphaseKernel(const Filter & f,uint srcLength,uint dstLength,int samples)541 PolyphaseKernel::PolyphaseKernel(const Filter & f, uint srcLength, uint dstLength, int samples/*= 32*/)
542 {
543 nvDebugCheck(samples > 0);
544
545 float scale = float(dstLength) / float(srcLength);
546 const float iscale = 1.0f / scale;
547
548 if (scale > 1) {
549 // Upsampling.
550 samples = 1;
551 scale = 1;
552 }
553
554 m_length = dstLength;
555 m_width = f.width() * iscale;
556 m_windowSize = (int)ceilf(m_width * 2) + 1;
557
558 m_data = new float[m_windowSize * m_length];
559 memset(m_data, 0, sizeof(float) * m_windowSize * m_length);
560
561 for (uint i = 0; i < m_length; i++)
562 {
563 const float center = (0.5f + i) * iscale;
564
565 const int left = (int)floorf(center - m_width);
566 const int right = (int)ceilf(center + m_width);
567 nvDebugCheck(right - left <= m_windowSize);
568
569 float total = 0.0f;
570 for (int j = 0; j < m_windowSize; j++)
571 {
572 const float sample = f.sampleBox(left + j - center, scale, samples);
573
574 m_data[i * m_windowSize + j] = sample;
575 total += sample;
576 }
577
578 // normalize weights.
579 for (int j = 0; j < m_windowSize; j++)
580 {
581 m_data[i * m_windowSize + j] /= total;
582 }
583 }
584 }
585
~PolyphaseKernel()586 PolyphaseKernel::~PolyphaseKernel()
587 {
588 delete [] m_data;
589 }
590
591
592 /// Print the kernel for debugging purposes.
debugPrint() const593 void PolyphaseKernel::debugPrint() const
594 {
595 for (uint i = 0; i < m_length; i++)
596 {
597 nvDebug("%d: ", i);
598 for (int j = 0; j < m_windowSize; j++)
599 {
600 nvDebug(" %6.4f", m_data[i * m_windowSize + j]);
601 }
602 nvDebug("\n");
603 }
604 }
605
606