1// tResample.cpp 
2// 
3// Resample an image using various filers like nearest-neighbour, box, bilinear, and various bicubics. 
4// 
5// Copyright (c) 2020, 2024 Tristan Grimmer. 
6// Permission to use, copy, modify, and/or distribute this software for any purpose with or without fee is hereby 
7// granted, provided that the above copyright notice and this permission notice appear in all copies. 
8// 
9// THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE INCLUDING ALL 
10// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT, 
11// INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER IN 
12// AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR 
13// PERFORMANCE OF THIS SOFTWARE. 
14 
15#include "Image/tResample.h" 
16#include "System/tPrint.h" 
17using namespace tMath; 
18 
19 
20namespace tImage 
21{ 
22 enum class FilterDirection 
23 { 
24 Horizontal, 
25 Vertical 
26 }; 
27 
28 struct FilterParams 
29 { 
30 FilterParams() : RatioH(0.0f), RatioV(0.0f) { } 
31 union { float RatioH; float CubicCoeffB; float LanczosA; }; 
32 union { float RatioV; float CubicCoeffC; }; 
33 }; 
34 
35 tPixel4b KernelFilterNearest (const tPixel4b* src, int srcW, int srcH, float x, float y, FilterDirection, tResampleEdgeMode, const FilterParams&); 
36 tPixel4b KernelFilterBox (const tPixel4b* src, int srcW, int srcH, float x, float y, FilterDirection, tResampleEdgeMode, const FilterParams&); 
37 tPixel4b KernelFilterBilinear (const tPixel4b* src, int srcW, int srcH, float x, float y, FilterDirection, tResampleEdgeMode, const FilterParams&); 
38 tPixel4b KernelFilterBicubic (const tPixel4b* src, int srcW, int srcH, float x, float y, FilterDirection, tResampleEdgeMode, const FilterParams&); 
39 tPixel4b KernelFilterLanczos (const tPixel4b* src, int srcW, int srcH, float x, float y, FilterDirection, tResampleEdgeMode, const FilterParams&); 
40 
41 int GetSrcIndex(int idx, int count, tResampleEdgeMode); 
42 float ComputeCubicWeight(float x, float b, float c); 
43 float ComputeLanczosWeight(float x, float a); 
44} 
45 
46 
47const char* tImage::tResampleFilterNames[int(tResampleFilter::NumFilters)+1] = 
48{ 
49 "Nearest Neighbour", 
50 "Box", 
51 "Bilinear", 
52 "Bicubic Standard", 
53 "Bicubic CatmullRom", 
54 "Bicubic Mitchell", 
55 "Bicubic Cardinal", 
56 "Bicubic BSpline", 
57 "Lanczos Narrow", 
58 "Lanczos Normal", 
59 "Lanczos Wide", 
60 "None" 
61}; 
62 
63 
64const char* tImage::tResampleFilterNamesSimple[int(tResampleFilter::NumFilters)+1] = 
65{ 
66 "nearest", 
67 "box", 
68 "bilinear", 
69 "bicubic", 
70 "bicubic_catmullrom", 
71 "bicubic_mitchell", 
72 "bicubic_cardinal", 
73 "bicubic_bspline", 
74 "lanczos_narrow", 
75 "lanczos", 
76 "lanczos_wide", 
77 "none" 
78}; 
79 
80 
81const char* tImage::tResampleEdgeModeNames[int(tResampleEdgeMode::NumEdgeModes)+1] = 
82{ 
83 "Clamp", 
84 "Wrap", 
85 "None" 
86}; 
87 
88 
89const char* tImage::tResampleEdgeModeNamesSimple[int(tResampleEdgeMode::NumEdgeModes)+1] = 
90{ 
91 "clamp", 
92 "wrap", 
93 "none" 
94}; 
95 
96 
97inline int tImage::GetSrcIndex(int idx, int count, tResampleEdgeMode edgeMode) 
98{ 
99 tAssert(count > 0); 
100 switch (edgeMode) 
101 { 
102 case tResampleEdgeMode::Clamp: 
103 return tClamp(val: idx, min: 0, max: count-1); 
104 
105 case tResampleEdgeMode::Wrap: 
106 return tMod(n: idx, d: count); 
107 } 
108 
109 return -1; 
110} 
111 
112 
113bool tImage::Resample 
114( 
115 tPixel4b* src, int srcW, int srcH, 
116 tPixel4b* dst, int dstW, int dstH, 
117 tResampleFilter resampleFilter, 
118 tResampleEdgeMode edgeMode 
119) 
120{ 
121 if (!src || !dst || srcW<=0 || srcH<=0 || dstW<=0 || dstH<=0) 
122 return false; 
123 
124 if ((srcW == dstW) && (srcH == dstH)) 
125 { 
126 for (int p = 0; p < srcW*srcH; p++) 
127 dst[p] = src[p]; 
128 
129 return true; 
130 } 
131 
132 float ratioH = (dstW > 1) ? (float(srcW) - 1.0f) / float(dstW - 1) : 1.0f; 
133 float ratioV = (dstH > 1) ? (float(srcH) - 1.0f) / float(dstH - 1) : 1.0f; 
134 
135 // Decide what filer kernel to use. Different kernels may set different values in FilterParams. 
136 FilterParams params; 
137 typedef tPixel4b (*KernelFilterFn)(const tPixel4b* src, int srcW, int srcH, float x, float y, FilterDirection, tResampleEdgeMode, const FilterParams&); 
138 KernelFilterFn kernel; 
139 switch (resampleFilter) 
140 { 
141 case tResampleFilter::Nearest: 
142 kernel = KernelFilterNearest; 
143 break; 
144 
145 case tResampleFilter::Box: 
146 params.RatioH = ratioH; 
147 params.RatioV = ratioV; 
148 kernel = KernelFilterBox; 
149 break; 
150 
151 case tResampleFilter::Bilinear: 
152 kernel = KernelFilterBilinear; 
153 break; 
154 
155 case tResampleFilter::Bicubic_Standard: // Cardinal. B=0 C=3/4 
156 params.CubicCoeffB = 0.0f; 
157 params.CubicCoeffC = 3.0f/4.0f; 
158 kernel = KernelFilterBicubic; 
159 break; 
160 
161 case tResampleFilter::Bicubic_CatmullRom: // Cardinal. B=0 C=1/2 
162 params.CubicCoeffB = 0.0f; 
163 params.CubicCoeffC = 1.0f/2.0f; 
164 kernel = KernelFilterBicubic; 
165 break; 
166 
167 case tResampleFilter::Bicubic_Mitchell: // Balanced. B=1/3 C=1/3 
168 params.CubicCoeffB = 1.0f/3.0f; 
169 params.CubicCoeffC = 1.0f/3.0f; 
170 kernel = KernelFilterBicubic; 
171 break; 
172 
173 case tResampleFilter::Bicubic_Cardinal: // Pure Cardinal. B=0 C=1 
174 params.CubicCoeffB = 0.0f; 
175 params.CubicCoeffC = 1.0f; 
176 kernel = KernelFilterBicubic; 
177 break; 
178 
179 case tResampleFilter::Bicubic_BSpline: // Pure BSpline. Blurry. B=1 C=0 
180 params.CubicCoeffB = 1.0f; 
181 params.CubicCoeffC = 0.0f; 
182 kernel = KernelFilterBicubic; 
183 break; 
184 
185 case tResampleFilter::Lanczos_Narrow: // Lanczos. Ringy/Sharp. A=2 
186 params.LanczosA = 2.0f; 
187 kernel = KernelFilterLanczos; 
188 break; 
189 
190 case tResampleFilter::Lanczos_Normal: // Lanczos. Ringy/Sharp. A=3 
191 params.LanczosA = 3.0f; 
192 kernel = KernelFilterLanczos; 
193 break; 
194 
195 case tResampleFilter::Lanczos_Wide: // Lanczos. Ringy/Sharp. A=4 
196 params.LanczosA = 4.0f; 
197 kernel = KernelFilterLanczos; 
198 break; 
199 
200 case tResampleFilter::Invalid: 
201 default: 
202 return false; 
203 }  
204 
205 // By convention do horizontal first. Outer loop is for each src row. 
206 // hri stands for hozontal-resized-image. 
207 tPixel4b* hri = new tPixel4b[dstW*srcH]; 
208 for (int r = 0; r < srcH; r++) 
209 { 
210 // Fill in each dst pixel for the src row, 
211 float y = float(r); 
212 for (int c = 0; c < dstW; c++) 
213 { 
214 tPixel4b& dstPixel = hri[dstW*r + c]; 
215 float x = float(c) * ratioH; 
216 dstPixel = kernel(src, srcW, srcH, x, y, FilterDirection::Horizontal, edgeMode, params); 
217 } 
218 } 
219 
220 // Vertical resampling. Source is the horizontally resized image. 
221 for (int c = 0; c < dstW; c++) 
222 { 
223 float x = float(c); 
224 for (int r = 0; r < dstH; r++) 
225 { 
226 tPixel4b& dstPixel = dst[dstW*r + c]; 
227 float y = float(r) * ratioV; 
228 dstPixel = kernel(hri, dstW, srcH, x, y, FilterDirection::Vertical, edgeMode, params); 
229 } 
230 } 
231 
232 delete[] hri; 
233 return true; 
234} 
235 
236 
237tPixel4b tImage::KernelFilterNearest 
238( 
239 const tPixel4b* src, int srcW, int srcH, float x, float y, 
240 FilterDirection dir, tResampleEdgeMode edgeMode, const FilterParams& params 
241) 
242{ 
243 int ix = tClamp(val: int(x + 0.5f), min: 0, max: srcW-1); 
244 int iy = tClamp(val: int(y + 0.5f), min: 0, max: srcH-1); 
245 return src[srcW*iy + ix]; 
246} 
247 
248 
249tPixel4b tImage::KernelFilterBox 
250( 
251 const tPixel4b* src, int srcW, int srcH, float x, float y, 
252 FilterDirection dir, tResampleEdgeMode edgeMode, const FilterParams& params 
253) 
254{ 
255 float ratio = (dir == FilterDirection::Horizontal) ? params.RatioH : params.RatioV; 
256 int pixelDist = int(ratio + 1.0f); 
257 float maxDist = ratio; 
258 float weightTotal = 0.0f; 
259 tVector4 sampleTotal = tVector4::zero; 
260 
261 for (int ks = 1 - pixelDist; ks <= pixelDist; ks++) 
262 { 
263 int ix = (dir == FilterDirection::Horizontal) ? int(x) + ks : int(x); 
264 int iy = (dir == FilterDirection::Horizontal) ? int(y) : int(y) + ks; 
265 
266 float dist = tAbs( val: (dir == FilterDirection::Horizontal) ? x - float(ix) : y - float(iy) ); 
267 float weight = 0.0f; 
268 
269 int srcX = GetSrcIndex(idx: ix ,count: srcW, edgeMode); 
270 int srcY = GetSrcIndex(idx: iy ,count: srcH, edgeMode); 
271 tPixel4b srcPixel = src[srcW*srcY + srcX]; 
272 
273 if (ratio >= 1.0f) 
274 { 
275 weight = 1.0f - tMin(a: maxDist, b: dist)/maxDist; 
276 } 
277 else 
278 { 
279 if (dist >= (0.5f - ratio)) 
280 weight = 1.0f - dist; 
281 else 
282 return srcPixel; // Box is inside src pixel. Done. 
283 } 
284 
285 sampleTotal.x += srcPixel.R * weight; 
286 sampleTotal.y += srcPixel.G * weight; 
287 sampleTotal.z += srcPixel.B * weight; 
288 sampleTotal.w += srcPixel.A * weight; 
289 weightTotal += weight; 
290 } 
291 
292 // Renormalize sampleTotal back to [0, 256). 
293 sampleTotal /= weightTotal; 
294 return tPixel4b 
295 ( 
296 tClamp(val: int(tRound(v: sampleTotal.x)), min: 0, max: 255), 
297 tClamp(val: int(tRound(v: sampleTotal.y)), min: 0, max: 255), 
298 tClamp(val: int(tRound(v: sampleTotal.z)), min: 0, max: 255), 
299 tClamp(val: int(tRound(v: sampleTotal.w)), min: 0, max: 255) 
300 ); 
301} 
302 
303 
304tPixel4b tImage::KernelFilterBilinear 
305( 
306 const tPixel4b* src, int srcW, int srcH, float x, float y, 
307 FilterDirection dir, tResampleEdgeMode edgeMode, const FilterParams& params 
308) 
309{ 
310 int ix = int(x); 
311 int iy = int(y); 
312 
313 int srcXa = GetSrcIndex(idx: ix , count: srcW, edgeMode); 
314 int srcYa = GetSrcIndex(idx: iy , count: srcH, edgeMode); 
315 int srcXb = GetSrcIndex(idx: ix+1, count: srcW, edgeMode); 
316 int srcYb = GetSrcIndex(idx: iy+1, count: srcH, edgeMode); 
317 
318 tPixel4b a = src[srcW*srcYa + srcXa]; 
319 tPixel4b b = (dir == FilterDirection::Horizontal) ? 
320 src[srcW*srcYa + srcXb] : 
321 src[srcW*srcYb + srcXa]; 
322 
323 float weight = (dir == FilterDirection::Horizontal) ? float(x)-ix : float(y)-iy; 
324 
325 tVector4 av, bv, rv; 
326 a.GetDenorm(dest&: av); 
327 b.GetDenorm(dest&: bv); 
328 rv = av*(1.0f-weight) + bv*weight; 
329 tiClamp(val&: tiRound(v&: rv.x), min: 0.0f, max: 255.0f); 
330 tiClamp(val&: tiRound(v&: rv.y), min: 0.0f, max: 255.0f); 
331 tiClamp(val&: tiRound(v&: rv.z), min: 0.0f, max: 255.0f); 
332 tiClamp(val&: tiRound(v&: rv.w), min: 0.0f, max: 255.0f); 
333 
334 return tPixel4b(int(rv.x), int(rv.y), int(rv.z), int(rv.w)); 
335} 
336 
337 
338// This function is the cubic filter workhorse. It implements the weight function k(x) found 
339// at https://entropymine.com/imageworsener/bicubic/ 
340// If that site ever goes down, the original paper is from Mitchell and Netravali (1988). 
341float tImage::ComputeCubicWeight(float x, float b, float c) 
342{ 
343 float xa = tAbs(val: x); 
344 
345 // Case 3. Early exit the 'otherwise' case. 
346 if (xa >= 2.0f) 
347 return 0.0f; 
348 
349 // Common terms in the other two cases. 
350 float c6 = 6.0f*c; 
351 float xc = tCube(v: xa); 
352 float b12 = 12.0f*b; 
353 float xs = tSquare(v: xa); 
354 float r = 0.0f; 
355 
356 if (xa < 1.0f) 
357 { 
358 // Case 1. 
359 float b9 = 9.0f*b; 
360 float b2 = 2.0f*b; 
361 r = (12.0f-b9-c6)*xc + (-18.0f+b12+c6)*xs + (6.0f-b2); 
362 } 
363 else 
364 { 
365 // Case 2. 
366 float b6 = 6.0f*b; 
367 float c30 = 30.0f*c; 
368 float c48 = 48.0f*c; 
369 float b8 = 8.0f*b; 
370 float c24 = 24.0f*c; 
371 r = (-b-c6)*xc + (b6+c30)*xs + (-b12-c48)*xa + (b8+c24); 
372 } 
373 
374 return tClampMin(val: r/6.0f, min: 0.0f); 
375} 
376 
377 
378tPixel4b tImage::KernelFilterBicubic 
379( 
380 const tPixel4b* src, int srcW, int srcH, float x, float y, 
381 FilterDirection dir, tResampleEdgeMode edgeMode, const FilterParams& params 
382) 
383{ 
384 float weightTotal = 0.0f; 
385 tVector4 sampleTotal = tVector4::zero; 
386 
387 for (int ks = -2; ks < 2; ks++) 
388 { 
389 int ix = (dir == FilterDirection::Horizontal) ? int(x) + ks : int(x); 
390 int iy = (dir == FilterDirection::Horizontal) ? int(y) : int(y) + ks; 
391 
392 float diff = (dir == FilterDirection::Horizontal) ? x - float(ix) : y - float(iy); 
393 float weight = ComputeCubicWeight(x: diff, b: params.CubicCoeffB, c: params.CubicCoeffC); 
394 
395 int srcX = GetSrcIndex(idx: ix, count: srcW, edgeMode); 
396 int srcY = GetSrcIndex(idx: iy, count: srcH, edgeMode); 
397 tPixel4b srcPixel = src[srcW*srcY + srcX]; 
398 
399 sampleTotal.x += srcPixel.R * weight; 
400 sampleTotal.y += srcPixel.G * weight; 
401 sampleTotal.z += srcPixel.B * weight; 
402 sampleTotal.w += srcPixel.A * weight; 
403 weightTotal += weight; 
404 } 
405 
406 // Renormalize sampleTotal back to [0, 256). 
407 sampleTotal /= weightTotal; 
408 return tPixel4b 
409 ( 
410 tClamp(val: int(tRound(v: sampleTotal.x)), min: 0, max: 255), 
411 tClamp(val: int(tRound(v: sampleTotal.y)), min: 0, max: 255), 
412 tClamp(val: int(tRound(v: sampleTotal.z)), min: 0, max: 255), 
413 tClamp(val: int(tRound(v: sampleTotal.w)), min: 0, max: 255) 
414 ); 
415} 
416 
417 
418inline float tImage::ComputeLanczosWeight(float x, float a) 
419{ 
420 // See https://en.wikipedia.org/wiki/Lanczos_resampling 
421 return ((x >= -a) && (x <= a)) ? tSinc(x) * tSinc(x: x/a) : 0.0f; 
422} 
423 
424 
425tPixel4b tImage::KernelFilterLanczos 
426( 
427 const tPixel4b* src, int srcW, int srcH, float x, float y, 
428 FilterDirection dir, tResampleEdgeMode edgeMode, const FilterParams& params 
429) 
430{ 
431 int pixelDist = int(params.LanczosA); 
432 float weightTotal = 0.0f; 
433 tVector4 sampleTotal = tVector4::zero; 
434 
435 for (int ks = -pixelDist; ks < pixelDist; ks++) 
436 { 
437 int ix = (dir == FilterDirection::Horizontal) ? int(x) + ks : int(x); 
438 int iy = (dir == FilterDirection::Horizontal) ? int(y) : int(y) + ks; 
439 
440 float diff = (dir == FilterDirection::Horizontal) ? x - float(ix) : y - float(iy); 
441 float weight = ComputeLanczosWeight(x: tAbs(val: diff), a: params.LanczosA); 
442 
443 int srcX = GetSrcIndex(idx: ix, count: srcW, edgeMode); 
444 int srcY = GetSrcIndex(idx: iy, count: srcH, edgeMode); 
445 tPixel4b srcPixel = src[srcW*srcY + srcX]; 
446 
447 sampleTotal.x += srcPixel.R * weight; 
448 sampleTotal.y += srcPixel.G * weight; 
449 sampleTotal.z += srcPixel.B * weight; 
450 sampleTotal.w += srcPixel.A * weight; 
451 weightTotal += weight; 
452 } 
453 
454 // Renormalize totalSamples back to [0, 256). 
455 sampleTotal /= weightTotal; 
456 return tPixel4b 
457 ( 
458 tClamp(val: int(tRound(v: sampleTotal.x)), min: 0, max: 255), 
459 tClamp(val: int(tRound(v: sampleTotal.y)), min: 0, max: 255), 
460 tClamp(val: int(tRound(v: sampleTotal.z)), min: 0, max: 255), 
461 tClamp(val: int(tRound(v: sampleTotal.w)), min: 0, max: 255) 
462 ); 
463} 
464