1// tResample.cpp 
2// 
3// Resample an image using various filers like nearest-neighbour, box, bilinear, and various bicubics. 
4// 
5// Copyright (c) 2020 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 tPixel KernelFilterNearest (const tPixel* src, int srcW, int srcH, float x, float y, FilterDirection, tResampleEdgeMode, const FilterParams&); 
36 tPixel KernelFilterBox (const tPixel* src, int srcW, int srcH, float x, float y, FilterDirection, tResampleEdgeMode, const FilterParams&); 
37 tPixel KernelFilterBilinear (const tPixel* src, int srcW, int srcH, float x, float y, FilterDirection, tResampleEdgeMode, const FilterParams&); 
38 tPixel KernelFilterBicubic (const tPixel* src, int srcW, int srcH, float x, float y, FilterDirection, tResampleEdgeMode, const FilterParams&); 
39 tPixel KernelFilterLanczos (const tPixel* 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)] = 
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}; 
61 
62 
63const char* tImage::tResampleEdgeModeNames[int(tResampleEdgeMode::NumEdgeModes)] = 
64{ 
65 "Clamp", 
66 "Wrap" 
67}; 
68 
69 
70inline int tImage::GetSrcIndex(int idx, int count, tResampleEdgeMode edgeMode) 
71{ 
72 tAssert(count > 0); 
73 switch (edgeMode) 
74 { 
75 case tResampleEdgeMode::Clamp: 
76 return tClamp(idx, 0, count-1); 
77 
78 case tResampleEdgeMode::Wrap: 
79 return tMod(idx, count); 
80 } 
81 
82 return -1; 
83} 
84 
85 
86bool tImage::Resample 
87( 
88 tPixel* src, int srcW, int srcH, 
89 tPixel* dst, int dstW, int dstH, 
90 tResampleFilter resampleFilter, 
91 tResampleEdgeMode edgeMode 
92) 
93{ 
94 if (!src || !dst || srcW<=0 || srcH<=0 || dstW<=0 || dstH<=0) 
95 return false; 
96 
97 if ((srcW == dstW) && (srcH == dstH)) 
98 { 
99 for (int p = 0; p < srcW*srcH; p++) 
100 dst[p] = src[p]; 
101 
102 return true; 
103 } 
104 
105 float ratioH = (dstW > 1) ? (float(srcW) - 1.0f) / float(dstW - 1) : 1.0f; 
106 float ratioV = (dstH > 1) ? (float(srcH) - 1.0f) / float(dstH - 1) : 1.0f; 
107 
108 // Decide what filer kernel to use. Different kernels may set different values in FilterParams. 
109 FilterParams params; 
110 typedef tPixel (*KernelFilterFn)(const tPixel* src, int srcW, int srcH, float x, float y, FilterDirection, tResampleEdgeMode, const FilterParams&); 
111 KernelFilterFn kernel; 
112 switch (resampleFilter) 
113 { 
114 case tResampleFilter::Nearest: 
115 kernel = KernelFilterNearest; 
116 break; 
117 
118 case tResampleFilter::Box: 
119 params.RatioH = ratioH; 
120 params.RatioV = ratioV; 
121 kernel = KernelFilterBox; 
122 break; 
123 
124 case tResampleFilter::Bilinear: 
125 kernel = KernelFilterBilinear; 
126 break; 
127 
128 case tResampleFilter::Bicubic_Standard: // Cardinal. B=0 C=3/4 
129 params.CubicCoeffB = 0.0f; 
130 params.CubicCoeffC = 3.0f/4.0f; 
131 kernel = KernelFilterBicubic; 
132 break; 
133 
134 case tResampleFilter::Bicubic_CatmullRom: // Cardinal. B=0 C=1/2 
135 params.CubicCoeffB = 0.0f; 
136 params.CubicCoeffC = 1.0f/2.0f; 
137 kernel = KernelFilterBicubic; 
138 break; 
139 
140 case tResampleFilter::Bicubic_Mitchell: // Balanced. B=1/3 C=1/3 
141 params.CubicCoeffB = 1.0f/3.0f; 
142 params.CubicCoeffC = 1.0f/3.0f; 
143 kernel = KernelFilterBicubic; 
144 break; 
145 
146 case tResampleFilter::Bicubic_Cardinal: // Pure Cardinal. B=0 C=1 
147 params.CubicCoeffB = 0.0f; 
148 params.CubicCoeffC = 1.0f; 
149 kernel = KernelFilterBicubic; 
150 break; 
151 
152 case tResampleFilter::Bicubic_BSpline: // Pure BSpline. Blurry. B=1 C=0 
153 params.CubicCoeffB = 1.0f; 
154 params.CubicCoeffC = 0.0f; 
155 kernel = KernelFilterBicubic; 
156 break; 
157 
158 case tResampleFilter::Lanczos_Narrow: // Lanczos. Ringy/Sharp. A=2 
159 params.LanczosA = 2.0f; 
160 kernel = KernelFilterLanczos; 
161 break; 
162 
163 case tResampleFilter::Lanczos_Normal: // Lanczos. Ringy/Sharp. A=3 
164 params.LanczosA = 3.0f; 
165 kernel = KernelFilterLanczos; 
166 break; 
167 
168 case tResampleFilter::Lanczos_Wide: // Lanczos. Ringy/Sharp. A=4 
169 params.LanczosA = 4.0f; 
170 kernel = KernelFilterLanczos; 
171 break; 
172 
173 case tResampleFilter::Invalid: 
174 default: 
175 return false; 
176 }  
177 
178 // By convention do horizontal first. Outer loop is for each src row. 
179 // hri stands for hozontal-resized-image. 
180 tPixel* hri = new tPixel[dstW*srcH]; 
181 for (int r = 0; r < srcH; r++) 
182 { 
183 // Fill in each dst pixel for the src row, 
184 float y = float(r); 
185 for (int c = 0; c < dstW; c++) 
186 { 
187 tPixel& dstPixel = hri[dstW*r + c]; 
188 float x = float(c) * ratioH; 
189 dstPixel = kernel(src, srcW, srcH, x, y, FilterDirection::Horizontal, edgeMode, params); 
190 } 
191 } 
192 
193 // Vertical resampling. Source is the horizontally resized image. 
194 for (int c = 0; c < dstW; c++) 
195 { 
196 float x = float(c); 
197 for (int r = 0; r < dstH; r++) 
198 { 
199 tPixel& dstPixel = dst[dstW*r + c]; 
200 float y = float(r) * ratioV; 
201 dstPixel = kernel(hri, dstW, srcH, x, y, FilterDirection::Vertical, edgeMode, params); 
202 } 
203 } 
204 
205 delete[] hri; 
206 return true; 
207} 
208 
209 
210tPixel tImage::KernelFilterNearest 
211( 
212 const tPixel* src, int srcW, int srcH, float x, float y, 
213 FilterDirection dir, tResampleEdgeMode edgeMode, const FilterParams& params 
214) 
215{ 
216 int ix = tClamp(int(x + 0.5f), 0, srcW-1); 
217 int iy = tClamp(int(y + 0.5f), 0, srcH-1); 
218 return src[srcW*iy + ix]; 
219} 
220 
221 
222tPixel tImage::KernelFilterBox 
223( 
224 const tPixel* src, int srcW, int srcH, float x, float y, 
225 FilterDirection dir, tResampleEdgeMode edgeMode, const FilterParams& params 
226) 
227{ 
228 float ratio = (dir == FilterDirection::Horizontal) ? params.RatioH : params.RatioV; 
229 int pixelDist = int(ratio + 1.0f); 
230 float maxDist = ratio; 
231 float weightTotal = 0.0f; 
232 tVector4 sampleTotal = tVector4::zero; 
233 
234 for (int ks = 1 - pixelDist; ks <= pixelDist; ks++) 
235 { 
236 int ix = (dir == FilterDirection::Horizontal) ? int(x) + ks : int(x); 
237 int iy = (dir == FilterDirection::Horizontal) ? int(y) : int(y) + ks; 
238 
239 float dist = tAbs( (dir == FilterDirection::Horizontal) ? x - float(ix) : y - float(iy) ); 
240 float weight = 0.0f; 
241 
242 int srcX = GetSrcIndex(ix ,srcW, edgeMode); 
243 int srcY = GetSrcIndex(iy ,srcH, edgeMode); 
244 tPixel srcPixel = src[srcW*srcY + srcX]; 
245 
246 if (ratio >= 1.0f) 
247 { 
248 weight = 1.0f - tMin(maxDist, dist)/maxDist; 
249 } 
250 else 
251 { 
252 if (dist >= (0.5f - ratio)) 
253 weight = 1.0f - dist; 
254 else 
255 return srcPixel; // Box is inside src pixel. Done. 
256 } 
257 
258 sampleTotal.x += srcPixel.R * weight; 
259 sampleTotal.y += srcPixel.G * weight; 
260 sampleTotal.z += srcPixel.B * weight; 
261 sampleTotal.w += srcPixel.A * weight; 
262 weightTotal += weight; 
263 } 
264 
265 // Renormalize sampleTotal back to [0, 256). 
266 sampleTotal /= weightTotal; 
267 return tPixel 
268 ( 
269 tClamp(int(tRound(sampleTotal.x)), 0, 255), 
270 tClamp(int(tRound(sampleTotal.y)), 0, 255), 
271 tClamp(int(tRound(sampleTotal.z)), 0, 255), 
272 tClamp(int(tRound(sampleTotal.w)), 0, 255) 
273 ); 
274} 
275 
276 
277tPixel tImage::KernelFilterBilinear 
278( 
279 const tPixel* src, int srcW, int srcH, float x, float y, 
280 FilterDirection dir, tResampleEdgeMode edgeMode, const FilterParams& params 
281) 
282{ 
283 int ix = int(x); 
284 int iy = int(y); 
285 
286 int srcXa = GetSrcIndex(ix , srcW, edgeMode); 
287 int srcYa = GetSrcIndex(iy , srcH, edgeMode); 
288 int srcXb = GetSrcIndex(ix+1, srcW, edgeMode); 
289 int srcYb = GetSrcIndex(iy+1, srcH, edgeMode); 
290 
291 tPixel a = src[srcW*srcYa + srcXa]; 
292 tPixel b = (dir == FilterDirection::Horizontal) ? 
293 src[srcW*srcYa + srcXb] : 
294 src[srcW*srcYb + srcXa]; 
295 
296 float weight = (dir == FilterDirection::Horizontal) ? float(x)-ix : float(y)-iy; 
297 
298 tVector4 av, bv, rv; 
299 a.GetDenorm(av); 
300 b.GetDenorm(bv); 
301 rv = av*(1.0f-weight) + bv*weight; 
302 tiClamp(tiRound(rv.x), 0.0f, 255.0f); 
303 tiClamp(tiRound(rv.y), 0.0f, 255.0f); 
304 tiClamp(tiRound(rv.z), 0.0f, 255.0f); 
305 tiClamp(tiRound(rv.w), 0.0f, 255.0f); 
306 
307 return tPixel(int(rv.x), int(rv.y), int(rv.z), int(rv.w)); 
308} 
309 
310 
311// This function is the cubic filter workhorse. It implements the weight function k(x) found 
312// at https://entropymine.com/imageworsener/bicubic/ 
313// If that site ever goes down, the original paper is from Mitchell and Netravali (1988). 
314float tImage::ComputeCubicWeight(float x, float b, float c) 
315{ 
316 float xa = tAbs(x); 
317 
318 // Case 3. Early exit the 'otherwise' case. 
319 if (xa >= 2.0f) 
320 return 0.0f; 
321 
322 // Common terms in the other two cases. 
323 float c6 = 6.0f*c; 
324 float xc = tCube(xa); 
325 float b12 = 12.0f*b; 
326 float xs = tSquare(xa); 
327 float r = 0.0f; 
328 
329 if (xa < 1.0f) 
330 { 
331 // Case 1. 
332 float b9 = 9.0f*b; 
333 float b2 = 2.0f*b; 
334 r = (12.0f-b9-c6)*xc + (-18.0f+b12+c6)*xs + (6.0f-b2); 
335 } 
336 else 
337 { 
338 // Case 2. 
339 float b6 = 6.0f*b; 
340 float c30 = 30.0f*c; 
341 float c48 = 48.0f*c; 
342 float b8 = 8.0f*b; 
343 float c24 = 24.0f*c; 
344 r = (-b-c6)*xc + (b6+c30)*xs + (-b12-c48)*xa + (b8+c24); 
345 } 
346 
347 return tClampMin(r/6.0f, 0.0f); 
348} 
349 
350 
351tPixel tImage::KernelFilterBicubic 
352( 
353 const tPixel* src, int srcW, int srcH, float x, float y, 
354 FilterDirection dir, tResampleEdgeMode edgeMode, const FilterParams& params 
355) 
356{ 
357 float weightTotal = 0.0f; 
358 tVector4 sampleTotal = tVector4::zero; 
359 
360 for (int ks = -2; ks < 2; ks++) 
361 { 
362 int ix = (dir == FilterDirection::Horizontal) ? int(x) + ks : int(x); 
363 int iy = (dir == FilterDirection::Horizontal) ? int(y) : int(y) + ks; 
364 
365 float diff = (dir == FilterDirection::Horizontal) ? x - float(ix) : y - float(iy); 
366 float weight = ComputeCubicWeight(diff, params.CubicCoeffB, params.CubicCoeffC); 
367 
368 int srcX = GetSrcIndex(ix, srcW, edgeMode); 
369 int srcY = GetSrcIndex(iy, srcH, edgeMode); 
370 tPixel srcPixel = src[srcW*srcY + srcX]; 
371 
372 sampleTotal.x += srcPixel.R * weight; 
373 sampleTotal.y += srcPixel.G * weight; 
374 sampleTotal.z += srcPixel.B * weight; 
375 sampleTotal.w += srcPixel.A * weight; 
376 weightTotal += weight; 
377 } 
378 
379 // Renormalize sampleTotal back to [0, 256). 
380 sampleTotal /= weightTotal; 
381 return tPixel 
382 ( 
383 tClamp(int(tRound(sampleTotal.x)), 0, 255), 
384 tClamp(int(tRound(sampleTotal.y)), 0, 255), 
385 tClamp(int(tRound(sampleTotal.z)), 0, 255), 
386 tClamp(int(tRound(sampleTotal.w)), 0, 255) 
387 ); 
388} 
389 
390 
391inline float tImage::ComputeLanczosWeight(float x, float a) 
392{ 
393 // See https://en.wikipedia.org/wiki/Lanczos_resampling 
394 return ((x >= -a) && (x <= a)) ? tSinc(x) * tSinc(x/a) : 0.0f; 
395} 
396 
397 
398tPixel tImage::KernelFilterLanczos 
399( 
400 const tPixel* src, int srcW, int srcH, float x, float y, 
401 FilterDirection dir, tResampleEdgeMode edgeMode, const FilterParams& params 
402) 
403{ 
404 int pixelDist = int(params.LanczosA); 
405 float weightTotal = 0.0f; 
406 tVector4 sampleTotal = tVector4::zero; 
407 
408 for (int ks = -pixelDist; ks < pixelDist; ks++) 
409 { 
410 int ix = (dir == FilterDirection::Horizontal) ? int(x) + ks : int(x); 
411 int iy = (dir == FilterDirection::Horizontal) ? int(y) : int(y) + ks; 
412 
413 float diff = (dir == FilterDirection::Horizontal) ? x - float(ix) : y - float(iy); 
414 float weight = ComputeLanczosWeight(tAbs(diff), params.LanczosA); 
415 
416 int srcX = GetSrcIndex(ix, srcW, edgeMode); 
417 int srcY = GetSrcIndex(iy, srcH, edgeMode); 
418 tPixel srcPixel = src[srcW*srcY + srcX]; 
419 
420 sampleTotal.x += srcPixel.R * weight; 
421 sampleTotal.y += srcPixel.G * weight; 
422 sampleTotal.z += srcPixel.B * weight; 
423 sampleTotal.w += srcPixel.A * weight; 
424 weightTotal += weight; 
425 } 
426 
427 // Renormalize totalSamples back to [0, 256). 
428 sampleTotal /= weightTotal; 
429 return tPixel 
430 ( 
431 tClamp(int(tRound(sampleTotal.x)), 0, 255), 
432 tClamp(int(tRound(sampleTotal.y)), 0, 255), 
433 tClamp(int(tRound(sampleTotal.z)), 0, 255), 
434 tClamp(int(tRound(sampleTotal.w)), 0, 255) 
435 ); 
436} 
437