1 /++
2 	Image resizing support for [arsd.color.MemoryImage]. Handles up and down scaling.
3 	See [imageResize] for the main function, all others are lower level if you need
4 	more control.
5 
6 
7 	Note that this focuses more on quality than speed. You can tweak the `filterScale`
8 	argument to speed things up at the expense of quality though (lower number = faster).
9 
10 	I've found:
11 
12 	---
13 	auto size = calculateSizeKeepingAspectRatio(i.width, i.height, maxWidth, maxHeight);
14 	if(size.width != i.width || size.height != i.height) {
15 		i = imageResize(i, size.width, size.height, null, 1.0, 0.6);
16 	}
17 	---
18 
19 	Gives decent results balancing quality and speed. Compiling with ldc or gdc can also
20 	speed up your program.
21 
22 
23 
24 	Authors:
25 		Originally written in C by Rich Geldreich, ported to D by ketmar.
26 	License:
27 		Public Domain / Unlicense - http://unlicense.org/
28 +/
29 module arsd.imageresize;
30 
31 import arsd.color;
32 
33 // ////////////////////////////////////////////////////////////////////////// //
34 // Separable filtering image rescaler v2.21, Rich Geldreich - richgel99@gmail.com
35 //
36 // This is free and unencumbered software released into the public domain.
37 //
38 // Anyone is free to copy, modify, publish, use, compile, sell, or
39 // distribute this software, either in source code form or as a compiled
40 // binary, for any purpose, commercial or non-commercial, and by any
41 // means.
42 //
43 // In jurisdictions that recognize copyright laws, the author or authors
44 // of this software dedicate any and all copyright interest in the
45 // software to the public domain. We make this dedication for the benefit
46 // of the public at large and to the detriment of our heirs and
47 // successors. We intend this dedication to be an overt act of
48 // relinquishment in perpetuity of all present and future rights to this
49 // software under copyright law.
50 //
51 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
52 // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
53 // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
54 // IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR
55 // OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE,
56 // ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR
57 // OTHER DEALINGS IN THE SOFTWARE.
58 //
59 // For more information, please refer to <http://unlicense.org/>
60 //
61 // Feb. 1996: Creation, losely based on a heavily bugfixed version of Schumacher's resampler in Graphics Gems 3.
62 // Oct. 2000: Ported to C++, tweaks.
63 // May 2001: Continous to discrete mapping, box filter tweaks.
64 // March 9, 2002: Kaiser filter grabbed from Jonathan Blow's GD magazine mipmap sample code.
65 // Sept. 8, 2002: Comments cleaned up a bit.
66 // Dec. 31, 2008: v2.2: Bit more cleanup, released as public domain.
67 // June 4, 2012: v2.21: Switched to unlicense.org, integrated GCC fixes supplied by Peter Nagy <petern@crytek.com>, Anteru at anteru.net, and clay@coge.net,
68 // added Codeblocks project (for testing with MinGW and GCC), VS2008 static code analysis pass.
69 // float or double
70 private:
71 
72 //version = iresample_debug;
73 
74 
75 // ////////////////////////////////////////////////////////////////////////// //
76 public enum ImageResizeDefaultFilter = "lanczos4"; /// Default filter for image resampler.
77 public enum ImageResizeMaxDimension = 65536; /// Maximum image width/height for image resampler.
78 
79 
80 // ////////////////////////////////////////////////////////////////////////// //
81 /// Number of known image resizer filters.
82 public @property int imageResizeFilterCount () { pragma(inline, true); return NumFilters; }
83 
84 /// Get filter name. Will return `null` for invalid index.
85 public string imageResizeFilterName (long idx) { pragma(inline, true); return (idx >= 0 && idx < NumFilters ? gFilters.ptr[cast(uint)idx].name : null); }
86 
87 /// Find filter index by name. Will use default filter for invalid names.
88 public int imageResizeFindFilter (const(char)[] name, const(char)[] defaultFilter=ImageResizeDefaultFilter) {
89   int res = resamplerFindFilterInternal(name);
90   if (res >= 0) return res;
91   res = resamplerFindFilterInternal(defaultFilter);
92   if (res >= 0) return res;
93   res = resamplerFindFilterInternal("lanczos4");
94   assert(res >= 0);
95   return res;
96 }
97 
98 /++
99 	Calculates a new size that fits inside the maximums while keeping the original aspect ratio.
100 
101 	History:
102 		Added March 18, 2021 (dub v9.4)
103 +/
104 public Size calculateSizeKeepingAspectRatio(int currentWidth, int currentHeight, int maxWidth, int maxHeight) {
105 	if(currentWidth <= maxWidth && currentHeight <= maxHeight)
106 		return Size(currentWidth, currentHeight);
107 
108 	float shrinkage = 1.0;
109 
110 	if(currentWidth > maxWidth) {
111 		shrinkage = cast(float) maxWidth / currentWidth;
112 	}
113 	if(currentHeight > maxHeight) {
114 		auto shrinkage2 = cast(float) maxHeight / currentHeight;
115 		if(shrinkage2 < shrinkage)
116 			shrinkage = shrinkage2;
117 	}
118 
119 	return Size(cast(int) (currentWidth * shrinkage), cast(int) (currentHeight * shrinkage));
120 }
121 
122 // ////////////////////////////////////////////////////////////////////////// //
123 /// Resize image.
124 public TrueColorImage imageResize(int Components=4) (MemoryImage msrcimg, int dstwdt, int dsthgt, const(char)[] filter=null, float gamma=1.0f, float filterScale=1.0f) {
125   static assert(Components == 1 || Components == 3 || Components == 4, "invalid number of components in color");
126   return imageResize!Components(msrcimg, dstwdt, dsthgt, imageResizeFindFilter(filter), gamma, filterScale);
127 }
128 
129 /// ditto
130 public TrueColorImage imageResize(int Components=4) (MemoryImage msrcimg, int dstwdt, int dsthgt, int filter, float gamma=1.0f, float filterScale=1.0f) {
131   static assert(Components == 1 || Components == 3 || Components == 4, "invalid number of components in color");
132   if (msrcimg is null || msrcimg.width < 1 || msrcimg.height < 1 || msrcimg.width > ImageResizeMaxDimension || msrcimg.height > ImageResizeMaxDimension) {
133     throw new Exception("invalid source image");
134   }
135   if (dstwdt < 1 || dsthgt < 1 || dstwdt > ImageResizeMaxDimension || dsthgt > ImageResizeMaxDimension) throw new Exception("invalid destination image size");
136   auto resimg = new TrueColorImage(dstwdt, dsthgt);
137   scope(failure) .destroy(resimg);
138   if (auto tc = cast(TrueColorImage)msrcimg) {
139     imageResize!Components(
140       delegate (Color[] destrow, int y) { destrow[] = tc.imageData.colors[y*tc.width..(y+1)*tc.width]; },
141       delegate (int y, const(Color)[] row) { resimg.imageData.colors[y*resimg.width..(y+1)*resimg.width] = row[]; },
142       msrcimg.width, msrcimg.height, dstwdt, dsthgt, filter, gamma, filterScale
143     );
144   } else {
145     imageResize!Components(
146       delegate (Color[] destrow, int y) { foreach (immutable x, ref c; destrow) c = msrcimg.getPixel(cast(int)x, y); },
147       delegate (int y, const(Color)[] row) { resimg.imageData.colors[y*resimg.width..(y+1)*resimg.width] = row[]; },
148       msrcimg.width, msrcimg.height, dstwdt, dsthgt, filter, gamma, filterScale
149     );
150   }
151   return resimg;
152 }
153 
154 
155 private {
156   enum Linear2srgbTableSize = 4096;
157   enum InvLinear2srgbTableSize = cast(float)(1.0f/Linear2srgbTableSize);
158   float[256] srgb2linear = void;
159   ubyte[Linear2srgbTableSize] linear2srgb = void;
160   float lastGamma = float.nan;
161 }
162 
163 /// Resize image.
164 /// Partial gamma correction looks better on mips; set to 1.0 to disable gamma correction.
165 /// Filter scale: values < 1.0 cause aliasing, but create sharper looking mips (0.75f, for example).
166 public void imageResize(int Components=4) (
167   scope void delegate (Color[] destrow, int y) srcGetRow,
168   scope void delegate (int y, const(Color)[] row) dstPutRow,
169   int srcwdt, int srchgt, int dstwdt, int dsthgt,
170   int filter=-1, float gamma=1.0f, float filterScale=1.0f
171 ) {
172   static assert(Components == 1 || Components == 3 || Components == 4, "invalid number of components in color");
173   assert(srcGetRow !is null);
174   assert(dstPutRow !is null);
175 
176   if (srcwdt < 1 || srchgt < 1 || dstwdt < 1 || dsthgt < 1 ||
177       srcwdt > ImageResizeMaxDimension || srchgt > ImageResizeMaxDimension ||
178       dstwdt > ImageResizeMaxDimension || dsthgt > ImageResizeMaxDimension) throw new Exception("invalid image size");
179 
180   if (filter < 0 || filter >= NumFilters) {
181     filter = resamplerFindFilterInternal(ImageResizeDefaultFilter);
182     if (filter < 0) {
183       filter = resamplerFindFilterInternal("lanczos4");
184     }
185   }
186   assert(filter >= 0 && filter < NumFilters);
187 
188 
189   if (lastGamma != gamma) {
190     version(iresample_debug) { import core.stdc.stdio; stderr.fprintf("creating translation tables for gamma %f (previous gamma is %f)\n", gamma, lastGamma); }
191     foreach (immutable i, ref v; srgb2linear[]) {
192       import std.math : pow;
193       v = cast(float)pow(cast(int)i*1.0f/255.0f, gamma);
194     }
195     immutable float invSourceGamma = 1.0f/gamma;
196     foreach (immutable i, ref v; linear2srgb[]) {
197       import std.math : pow;
198       int k = cast(int)(255.0f*pow(cast(int)i*InvLinear2srgbTableSize, invSourceGamma)+0.5f);
199       if (k < 0) k = 0; else if (k > 255) k = 255;
200       v = cast(ubyte)k;
201     }
202     lastGamma = gamma;
203   }
204   version(iresample_debug) { import core.stdc.stdio; stderr.fprintf("filter is %d\n", filter); }
205 
206   ImageResampleWorker[Components] resamplers;
207   float[][Components] samples;
208   Color[] srcrow, dstrow;
209   scope(exit) {
210     foreach (ref rsm; resamplers[]) .destroy(rsm);
211     foreach (ref smr; samples[]) .destroy(smr);
212   }
213 
214   // now create a ImageResampleWorker instance for each component to process
215   // the first instance will create new contributor tables, which are shared by the resamplers
216   // used for the other components (a memory and slight cache efficiency optimization).
217   resamplers[0] = new ImageResampleWorker(srcwdt, srchgt, dstwdt, dsthgt, ImageResampleWorker.BoundaryClamp, 0.0f, 1.0f, filter, null, null, filterScale, filterScale);
218   samples[0].length = srcwdt;
219   srcrow.length = srcwdt;
220   dstrow.length = dstwdt;
221   foreach (immutable i; 1..Components) {
222     resamplers[i] = new ImageResampleWorker(srcwdt, srchgt, dstwdt, dsthgt, ImageResampleWorker.BoundaryClamp, 0.0f, 1.0f, filter, resamplers[0].getClistX(), resamplers[0].getClistY(), filterScale, filterScale);
223     samples[i].length = srcwdt;
224   }
225 
226   int dsty = 0;
227   foreach (immutable int srcy; 0..srchgt) {
228     // get row components
229     srcGetRow(srcrow, srcy);
230     {
231       auto scp = srcrow.ptr;
232       foreach (immutable x; 0..srcwdt) {
233         auto sc = *scp++;
234         samples.ptr[0].ptr[x] = srgb2linear.ptr[sc.r]; // first component
235         static if (Components > 1) samples.ptr[1].ptr[x] = srgb2linear.ptr[sc.g]; // second component
236         static if (Components > 2) samples.ptr[2].ptr[x] = srgb2linear.ptr[sc.b]; // thirs component
237         static if (Components == 4) samples.ptr[3].ptr[x] = sc.a*(1.0f/255.0f); // fourth component is alpha, and it is already linear
238       }
239     }
240 
241     foreach (immutable c; 0..Components) if (!resamplers.ptr[c].putLine(samples.ptr[c].ptr)) assert(0, "out of memory");
242 
243     for (;;) {
244       int compIdx = 0;
245       for (; compIdx < Components; ++compIdx) {
246         const(float)* outsmp = resamplers.ptr[compIdx].getLine();
247         if (outsmp is null) break;
248         auto dsc = dstrow.ptr;
249         // alpha?
250         static if (Components == 4) {
251           if (compIdx == 3) {
252             foreach (immutable x; 0..dstwdt) {
253               dsc.a = Color.clampToByte(cast(int)(255.0f*(*outsmp++)+0.5f));
254               ++dsc;
255             }
256             continue;
257           }
258         }
259         // color
260         auto dsb = (cast(ubyte*)dsc)+compIdx;
261         foreach (immutable x; 0..dstwdt) {
262           int j = cast(int)(Linear2srgbTableSize*(*outsmp++)+0.5f);
263           if (j < 0) j = 0; else if (j >= Linear2srgbTableSize) j = Linear2srgbTableSize-1;
264           *dsb = linear2srgb.ptr[j];
265           dsb += 4;
266         }
267       }
268       if (compIdx < Components) break;
269       // fill destination line
270       assert(dsty < dsthgt);
271       static if (Components != 4) {
272         auto dsc = dstrow.ptr;
273         foreach (immutable x; 0..dstwdt) {
274           static if (Components == 1) dsc.g = dsc.b = dsc.r;
275           dsc.a = 255;
276           ++dsc;
277         }
278       }
279       //version(iresample_debug) { import core.stdc.stdio; stderr.fprintf("writing dest row %d with %u components\n", dsty, Components); }
280       dstPutRow(dsty, dstrow);
281       ++dsty;
282     }
283   }
284 }
285 
286 
287 // ////////////////////////////////////////////////////////////////////////// //
288 public final class ImageResampleWorker {
289 nothrow @trusted @nogc:
290 public:
291   alias ResampleReal = float;
292   alias Sample = ResampleReal;
293 
294   static struct Contrib {
295     ResampleReal weight;
296     ushort pixel;
297   }
298 
299   static struct ContribList {
300     ushort n;
301     Contrib* p;
302   }
303 
304   alias BoundaryOp = int;
305   enum /*Boundary_Op*/ {
306     BoundaryWrap = 0,
307     BoundaryReflect = 1,
308     BoundaryClamp = 2,
309   }
310 
311   alias Status = int;
312   enum /*Status*/ {
313     StatusOkay = 0,
314     StatusOutOfMemory = 1,
315     StatusBadFilterName = 2,
316     StatusScanBufferFull = 3,
317   }
318 
319 private:
320   alias FilterFunc = ResampleReal function (ResampleReal) nothrow @trusted @nogc;
321 
322   int mIntermediateX;
323 
324   int mResampleSrcX;
325   int mResampleSrcY;
326   int mResampleDstX;
327   int mResampleDstY;
328 
329   BoundaryOp mBoundaryOp;
330 
331   Sample* mPdstBuf;
332   Sample* mPtmpBuf;
333 
334   ContribList* mPclistX;
335   ContribList* mPclistY;
336 
337   bool mClistXForced;
338   bool mClistYForced;
339 
340   bool mDelayXResample;
341 
342   int* mPsrcYCount;
343   ubyte* mPsrcYFlag;
344 
345   // The maximum number of scanlines that can be buffered at one time.
346   enum MaxScanBufSize = ImageResizeMaxDimension;
347 
348   static struct ScanBuf {
349     int[MaxScanBufSize] scanBufY;
350     Sample*[MaxScanBufSize] scanBufL;
351   }
352 
353   ScanBuf* mPscanBuf;
354 
355   int mCurSrcY;
356   int mCurDstY;
357 
358   Status mStatus;
359 
360   // The make_clist() method generates, for all destination samples,
361   // the list of all source samples with non-zero weighted contributions.
362   ContribList* makeClist(
363     int srcX, int dstX, BoundaryOp boundaryOp,
364     FilterFunc Pfilter,
365     ResampleReal filterSupport,
366     ResampleReal filterScale,
367     ResampleReal srcOfs)
368   {
369     import core.stdc.stdlib : calloc, free;
370     import std.math : floor, ceil;
371 
372     static struct ContribBounds {
373       // The center of the range in DISCRETE coordinates (pixel center = 0.0f).
374       ResampleReal center;
375       int left, right;
376     }
377 
378     ContribList* Pcontrib, PcontribRes;
379     Contrib* Pcpool;
380     Contrib* PcpoolNext;
381     ContribBounds* PcontribBounds;
382 
383     if ((Pcontrib = cast(ContribList*)calloc(dstX, ContribList.sizeof)) is null) return null;
384     scope(exit) if (Pcontrib !is null) free(Pcontrib);
385 
386     PcontribBounds = cast(ContribBounds*)calloc(dstX, ContribBounds.sizeof);
387     if (PcontribBounds is null) return null;
388     scope(exit) free(PcontribBounds);
389 
390     enum ResampleReal NUDGE = 0.5f;
391     immutable ResampleReal ooFilterScale = 1.0f/filterScale;
392     immutable ResampleReal xscale = dstX/cast(ResampleReal)srcX;
393 
394     if (xscale < 1.0f) {
395       int total = 0;
396       // Handle case when there are fewer destination samples than source samples (downsampling/minification).
397       // stretched half width of filter
398       immutable ResampleReal halfWidth = (filterSupport/xscale)*filterScale;
399       // Find the range of source sample(s) that will contribute to each destination sample.
400       foreach (immutable i; 0..dstX) {
401         // Convert from discrete to continuous coordinates, scale, then convert back to discrete.
402         ResampleReal center = (cast(ResampleReal)i+NUDGE)/xscale;
403         center -= NUDGE;
404         center += srcOfs;
405         immutable int left = castToInt(cast(ResampleReal)floor(center-halfWidth));
406         immutable int right = castToInt(cast(ResampleReal)ceil(center+halfWidth));
407         PcontribBounds[i].center = center;
408         PcontribBounds[i].left = left;
409         PcontribBounds[i].right = right;
410         total += (right-left+1);
411       }
412 
413       // Allocate memory for contributors.
414       if (total == 0 || ((Pcpool = cast(Contrib*)calloc(total, Contrib.sizeof)) is null)) return null;
415       //scope(failure) free(Pcpool);
416       //immutable int total = n;
417 
418       PcpoolNext = Pcpool;
419 
420       // Create the list of source samples which contribute to each destination sample.
421       foreach (immutable i; 0..dstX) {
422         int maxK = -1;
423         ResampleReal maxW = -1e+20f;
424 
425         ResampleReal center = PcontribBounds[i].center;
426         immutable int left = PcontribBounds[i].left;
427         immutable int right = PcontribBounds[i].right;
428 
429         Pcontrib[i].n = 0;
430         Pcontrib[i].p = PcpoolNext;
431         PcpoolNext += (right-left+1);
432         assert(PcpoolNext-Pcpool <= total);
433 
434         ResampleReal totalWeight0 = 0;
435         foreach (immutable j; left..right+1) totalWeight0 += Pfilter((center-cast(ResampleReal)j)*xscale*ooFilterScale);
436         immutable ResampleReal norm = cast(ResampleReal)(1.0f/totalWeight0);
437 
438         ResampleReal totalWeight1 = 0;
439         foreach (immutable j; left..right+1) {
440           immutable ResampleReal weight = Pfilter((center-cast(ResampleReal)j)*xscale*ooFilterScale)*norm;
441           if (weight == 0.0f) continue;
442           immutable int n = reflect(j, srcX, boundaryOp);
443           // Increment the number of source samples which contribute to the current destination sample.
444           immutable int k = Pcontrib[i].n++;
445           Pcontrib[i].p[k].pixel = cast(ushort)(n); // store src sample number
446           Pcontrib[i].p[k].weight = weight; // store src sample weight
447           totalWeight1 += weight; // total weight of all contributors
448           if (weight > maxW) {
449             maxW = weight;
450             maxK = k;
451           }
452         }
453         //assert(Pcontrib[i].n);
454         //assert(max_k != -1);
455         if (maxK == -1 || Pcontrib[i].n == 0) return null;
456         if (totalWeight1 != 1.0f) Pcontrib[i].p[maxK].weight += 1.0f-totalWeight1;
457       }
458     } else {
459       int total = 0;
460       // Handle case when there are more destination samples than source samples (upsampling).
461       immutable ResampleReal halfWidth = filterSupport*filterScale;
462       // Find the source sample(s) that contribute to each destination sample.
463       foreach (immutable i; 0..dstX) {
464         // Convert from discrete to continuous coordinates, scale, then convert back to discrete.
465         ResampleReal center = (cast(ResampleReal)i+NUDGE)/xscale;
466         center -= NUDGE;
467         center += srcOfs;
468         immutable int left = castToInt(cast(ResampleReal)floor(center-halfWidth));
469         immutable int right = castToInt(cast(ResampleReal)ceil(center+halfWidth));
470         PcontribBounds[i].center = center;
471         PcontribBounds[i].left = left;
472         PcontribBounds[i].right = right;
473         total += (right-left+1);
474       }
475 
476       // Allocate memory for contributors.
477       if (total == 0 || ((Pcpool = cast(Contrib*)calloc(total, Contrib.sizeof)) is null)) return null;
478       //scope(failure) free(Pcpool);
479 
480       PcpoolNext = Pcpool;
481 
482       // Create the list of source samples which contribute to each destination sample.
483       foreach (immutable i; 0..dstX) {
484         int maxK = -1;
485         ResampleReal maxW = -1e+20f;
486 
487         ResampleReal center = PcontribBounds[i].center;
488         immutable int left = PcontribBounds[i].left;
489         immutable int right = PcontribBounds[i].right;
490 
491         Pcontrib[i].n = 0;
492         Pcontrib[i].p = PcpoolNext;
493         PcpoolNext += (right-left+1);
494         assert(PcpoolNext-Pcpool <= total);
495 
496         ResampleReal totalWeight0 = 0;
497         foreach (immutable j; left..right+1) totalWeight0 += Pfilter((center-cast(ResampleReal)j)*ooFilterScale);
498         immutable ResampleReal norm = cast(ResampleReal)(1.0f/totalWeight0);
499 
500         ResampleReal totalWeight1 = 0;
501         foreach (immutable j; left..right+1) {
502           immutable ResampleReal weight = Pfilter((center-cast(ResampleReal)j)*ooFilterScale)*norm;
503           if (weight == 0.0f) continue;
504           immutable int n = reflect(j, srcX, boundaryOp);
505           // Increment the number of source samples which contribute to the current destination sample.
506           immutable int k = Pcontrib[i].n++;
507           Pcontrib[i].p[k].pixel = cast(ushort)(n); // store src sample number
508           Pcontrib[i].p[k].weight = weight; // store src sample weight
509           totalWeight1 += weight; // total weight of all contributors
510           if (weight > maxW) {
511             maxW = weight;
512             maxK = k;
513           }
514         }
515         //assert(Pcontrib[i].n);
516         //assert(max_k != -1);
517         if (maxK == -1 || Pcontrib[i].n == 0) return null;
518         if (totalWeight1 != 1.0f) Pcontrib[i].p[maxK].weight += 1.0f-totalWeight1;
519       }
520     }
521     // don't free return value
522     PcontribRes = Pcontrib;
523     Pcontrib = null;
524     return PcontribRes;
525   }
526 
527   static int countOps (const(ContribList)* Pclist, int k) {
528     int t = 0;
529     foreach (immutable i; 0..k) t += Pclist[i].n;
530     return t;
531   }
532 
533   private ResampleReal mLo;
534   private ResampleReal mHi;
535 
536   ResampleReal clampSample (ResampleReal f) const {
537     pragma(inline, true);
538     if (f < mLo) f = mLo; else if (f > mHi) f = mHi;
539     return f;
540   }
541 
542 public:
543   // src_x/src_y - Input dimensions
544   // dst_x/dst_y - Output dimensions
545   // boundary_op - How to sample pixels near the image boundaries
546   // sample_low/sample_high - Clamp output samples to specified range, or disable clamping if sample_low >= sample_high
547   // Pclist_x/Pclist_y - Optional pointers to contributor lists from another instance of a ImageResampleWorker
548   // src_x_ofs/src_y_ofs - Offset input image by specified amount (fractional values okay)
549   this(
550     int srcX, int srcY,
551     int dstX, int dstY,
552     BoundaryOp boundaryOp=BoundaryClamp,
553     ResampleReal sampleLow=0.0f, ResampleReal sampleHigh=0.0f,
554     int PfilterIndex=-1,
555     ContribList* PclistX=null,
556     ContribList* PclistY=null,
557     ResampleReal filterXScale=1.0f,
558     ResampleReal filterYScale=1.0f,
559     ResampleReal srcXOfs=0.0f,
560     ResampleReal srcYOfs=0.0f)
561   {
562     import core.stdc.stdlib : calloc, malloc;
563 
564     int i, j;
565     ResampleReal support;
566     FilterFunc func;
567 
568     assert(srcX > 0);
569     assert(srcY > 0);
570     assert(dstX > 0);
571     assert(dstY > 0);
572 
573     mLo = sampleLow;
574     mHi = sampleHigh;
575 
576     mDelayXResample = false;
577     mIntermediateX = 0;
578     mPdstBuf = null;
579     mPtmpBuf = null;
580     mClistXForced = false;
581     mPclistX = null;
582     mClistYForced = false;
583     mPclistY = null;
584     mPsrcYCount = null;
585     mPsrcYFlag = null;
586     mPscanBuf = null;
587     mStatus = StatusOkay;
588 
589     mResampleSrcX = srcX;
590     mResampleSrcY = srcY;
591     mResampleDstX = dstX;
592     mResampleDstY = dstY;
593 
594     mBoundaryOp = boundaryOp;
595 
596     if ((mPdstBuf = cast(Sample*)malloc(mResampleDstX*Sample.sizeof)) is null) {
597       mStatus = StatusOutOfMemory;
598       return;
599     }
600 
601     if (PfilterIndex < 0 || PfilterIndex >= NumFilters) {
602       PfilterIndex = resamplerFindFilterInternal(ImageResizeDefaultFilter);
603       if (PfilterIndex < 0 || PfilterIndex >= NumFilters) {
604         mStatus = StatusBadFilterName;
605         return;
606       }
607     }
608 
609     func = gFilters[PfilterIndex].func;
610     support = gFilters[PfilterIndex].support;
611 
612     // Create contributor lists, unless the user supplied custom lists.
613     if (PclistX is null) {
614       mPclistX = makeClist(mResampleSrcX, mResampleDstX, mBoundaryOp, func, support, filterXScale, srcXOfs);
615       if (mPclistX is null) {
616         mStatus = StatusOutOfMemory;
617         return;
618       }
619     } else {
620       mPclistX = PclistX;
621       mClistXForced = true;
622     }
623 
624     if (PclistY is null) {
625       mPclistY = makeClist(mResampleSrcY, mResampleDstY, mBoundaryOp, func, support, filterYScale, srcYOfs);
626       if (mPclistY is null) {
627         mStatus = StatusOutOfMemory;
628         return;
629       }
630     } else {
631       mPclistY = PclistY;
632       mClistYForced = true;
633     }
634 
635     if ((mPsrcYCount = cast(int*)calloc(mResampleSrcY, int.sizeof)) is null) {
636       mStatus = StatusOutOfMemory;
637       return;
638     }
639 
640     if ((mPsrcYFlag = cast(ubyte*)calloc(mResampleSrcY, ubyte.sizeof)) is null) {
641       mStatus = StatusOutOfMemory;
642       return;
643     }
644 
645     // Count how many times each source line contributes to a destination line.
646     for (i = 0; i < mResampleDstY; ++i) {
647       for (j = 0; j < mPclistY[i].n; ++j) {
648         ++mPsrcYCount[resamplerRangeCheck(mPclistY[i].p[j].pixel, mResampleSrcY)];
649       }
650     }
651 
652     if ((mPscanBuf = cast(ScanBuf*)malloc(ScanBuf.sizeof)) is null) {
653       mStatus = StatusOutOfMemory;
654       return;
655     }
656 
657     for (i = 0; i < MaxScanBufSize; ++i) {
658       mPscanBuf.scanBufY.ptr[i] = -1;
659       mPscanBuf.scanBufL.ptr[i] = null;
660     }
661 
662     mCurSrcY = mCurDstY = 0;
663     {
664       // Determine which axis to resample first by comparing the number of multiplies required
665       // for each possibility.
666       int xOps = countOps(mPclistX, mResampleDstX);
667       int yOps = countOps(mPclistY, mResampleDstY);
668 
669       // Hack 10/2000: Weight Y axis ops a little more than X axis ops.
670       // (Y axis ops use more cache resources.)
671       int xyOps = xOps*mResampleSrcY+(4*yOps*mResampleDstX)/3;
672       int yxOps = (4*yOps*mResampleSrcX)/3+xOps*mResampleDstY;
673 
674       // Now check which resample order is better. In case of a tie, choose the order
675       // which buffers the least amount of data.
676       if (xyOps > yxOps || (xyOps == yxOps && mResampleSrcX < mResampleDstX)) {
677         mDelayXResample = true;
678         mIntermediateX = mResampleSrcX;
679       } else {
680         mDelayXResample = false;
681         mIntermediateX = mResampleDstX;
682       }
683     }
684 
685     if (mDelayXResample) {
686       if ((mPtmpBuf = cast(Sample*)malloc(mIntermediateX*Sample.sizeof)) is null) {
687         mStatus = StatusOutOfMemory;
688         return;
689       }
690     }
691   }
692 
693   ~this () {
694      import core.stdc.stdlib : free;
695 
696      if (mPdstBuf !is null) {
697        free(mPdstBuf);
698        mPdstBuf = null;
699      }
700 
701      if (mPtmpBuf !is null) {
702        free(mPtmpBuf);
703        mPtmpBuf = null;
704      }
705 
706      // Don't deallocate a contibutor list if the user passed us one of their own.
707      if (mPclistX !is null && !mClistXForced) {
708        free(mPclistX.p);
709        free(mPclistX);
710        mPclistX = null;
711      }
712      if (mPclistY !is null && !mClistYForced) {
713        free(mPclistY.p);
714        free(mPclistY);
715        mPclistY = null;
716      }
717 
718      if (mPsrcYCount !is null) {
719        free(mPsrcYCount);
720        mPsrcYCount = null;
721      }
722 
723      if (mPsrcYFlag !is null) {
724        free(mPsrcYFlag);
725        mPsrcYFlag = null;
726      }
727 
728      if (mPscanBuf !is null) {
729        foreach (immutable i; 0..MaxScanBufSize) if (mPscanBuf.scanBufL.ptr[i] !is null) free(mPscanBuf.scanBufL.ptr[i]);
730        free(mPscanBuf);
731        mPscanBuf = null;
732      }
733   }
734 
735   // Reinits resampler so it can handle another frame.
736   void restart () {
737     import core.stdc.stdlib : free;
738     if (StatusOkay != mStatus) return;
739     mCurSrcY = mCurDstY = 0;
740     foreach (immutable i; 0..mResampleSrcY) {
741       mPsrcYCount[i] = 0;
742       mPsrcYFlag[i] = false;
743     }
744     foreach (immutable i; 0..mResampleDstY) {
745       foreach (immutable j; 0..mPclistY[i].n) {
746         ++mPsrcYCount[resamplerRangeCheck(mPclistY[i].p[j].pixel, mResampleSrcY)];
747       }
748     }
749     foreach (immutable i; 0..MaxScanBufSize) {
750       mPscanBuf.scanBufY.ptr[i] = -1;
751       free(mPscanBuf.scanBufL.ptr[i]);
752       mPscanBuf.scanBufL.ptr[i] = null;
753     }
754   }
755 
756   // false on out of memory.
757   bool putLine (const(Sample)* Psrc) {
758     int i;
759 
760     if (mCurSrcY >= mResampleSrcY) return false;
761 
762     // Does this source line contribute to any destination line? if not, exit now.
763     if (!mPsrcYCount[resamplerRangeCheck(mCurSrcY, mResampleSrcY)]) {
764       ++mCurSrcY;
765       return true;
766     }
767 
768     // Find an empty slot in the scanline buffer. (FIXME: Perf. is terrible here with extreme scaling ratios.)
769     for (i = 0; i < MaxScanBufSize; ++i) if (mPscanBuf.scanBufY.ptr[i] == -1) break;
770 
771     // If the buffer is full, exit with an error.
772     if (i == MaxScanBufSize) {
773       mStatus = StatusScanBufferFull;
774       return false;
775     }
776 
777     mPsrcYFlag[resamplerRangeCheck(mCurSrcY, mResampleSrcY)] = true;
778     mPscanBuf.scanBufY.ptr[i] = mCurSrcY;
779 
780     // Does this slot have any memory allocated to it?
781     if (!mPscanBuf.scanBufL.ptr[i]) {
782       import core.stdc.stdlib : malloc;
783       if ((mPscanBuf.scanBufL.ptr[i] = cast(Sample*)malloc(mIntermediateX*Sample.sizeof)) is null) {
784         mStatus = StatusOutOfMemory;
785         return false;
786       }
787     }
788 
789     // Resampling on the X axis first?
790     if (mDelayXResample) {
791       import core.stdc.string : memcpy;
792       assert(mIntermediateX == mResampleSrcX);
793       // Y-X resampling order
794       memcpy(mPscanBuf.scanBufL.ptr[i], Psrc, mIntermediateX*Sample.sizeof);
795     } else {
796       assert(mIntermediateX == mResampleDstX);
797       // X-Y resampling order
798       resampleX(mPscanBuf.scanBufL.ptr[i], Psrc);
799     }
800 
801     ++mCurSrcY;
802 
803     return true;
804   }
805 
806   // null if no scanlines are currently available (give the resampler more scanlines!)
807   const(Sample)* getLine () {
808     // if all the destination lines have been generated, then always return null
809     if (mCurDstY == mResampleDstY) return null;
810     // check to see if all the required contributors are present, if not, return null
811     foreach (immutable i; 0..mPclistY[mCurDstY].n) {
812       if (!mPsrcYFlag[resamplerRangeCheck(mPclistY[mCurDstY].p[i].pixel, mResampleSrcY)]) return null;
813     }
814     resampleY(mPdstBuf);
815     ++mCurDstY;
816     return mPdstBuf;
817   }
818 
819   @property Status status () const { pragma(inline, true); return mStatus; }
820 
821   // returned contributor lists can be shared with another ImageResampleWorker
822   void getClists (ContribList** ptrClistX, ContribList** ptrClistY) {
823     if (ptrClistX !is null) *ptrClistX = mPclistX;
824     if (ptrClistY !is null) *ptrClistY = mPclistY;
825   }
826 
827   @property ContribList* getClistX () { pragma(inline, true); return mPclistX; }
828   @property ContribList* getClistY () { pragma(inline, true); return mPclistY; }
829 
830   // filter accessors
831   static @property auto filters () {
832     static struct FilterRange {
833     pure nothrow @trusted @nogc:
834       int idx;
835       @property bool empty () const { pragma(inline, true); return (idx >= NumFilters); }
836       @property string front () const { pragma(inline, true); return (idx < NumFilters ? gFilters[idx].name : null); }
837       void popFront () { if (idx < NumFilters) ++idx; }
838       int length () const { return cast(int)NumFilters; }
839       alias opDollar = length;
840     }
841     return FilterRange();
842   }
843 
844 private:
845   /* Ensure that the contributing source sample is
846   * within bounds. If not, reflect, clamp, or wrap.
847   */
848   int reflect (in int j, in int srcX, in BoundaryOp boundaryOp) {
849     int n;
850     if (j < 0) {
851       if (boundaryOp == BoundaryReflect) {
852         n = -j;
853         if (n >= srcX) n = srcX-1;
854       } else if (boundaryOp == BoundaryWrap) {
855         n = posmod(j, srcX);
856       } else {
857         n = 0;
858       }
859     } else if (j >= srcX) {
860       if (boundaryOp == BoundaryReflect) {
861         n = (srcX-j)+(srcX-1);
862         if (n < 0) n = 0;
863       } else if (boundaryOp == BoundaryWrap) {
864         n = posmod(j, srcX);
865       } else {
866         n = srcX-1;
867       }
868     } else {
869       n = j;
870     }
871     return n;
872   }
873 
874   void resampleX (Sample* Pdst, const(Sample)* Psrc) {
875     assert(Pdst);
876     assert(Psrc);
877 
878     Sample total;
879     ContribList *Pclist = mPclistX;
880     Contrib *p;
881 
882     for (int i = mResampleDstX; i > 0; --i, ++Pclist) {
883       int j = void;
884       for (j = Pclist.n, p = Pclist.p, total = 0; j > 0; --j, ++p) total += Psrc[p.pixel]*p.weight;
885       *Pdst++ = total;
886     }
887   }
888 
889   void scaleYMov (Sample* Ptmp, const(Sample)* Psrc, ResampleReal weight, int dstX) {
890     // Not += because temp buf wasn't cleared.
891     for (int i = dstX; i > 0; --i) *Ptmp++ = *Psrc++*weight;
892   }
893 
894   void scaleYAdd (Sample* Ptmp, const(Sample)* Psrc, ResampleReal weight, int dstX) {
895     for (int i = dstX; i > 0; --i) (*Ptmp++) += *Psrc++*weight;
896   }
897 
898   void clamp (Sample* Pdst, int n) {
899     while (n > 0) {
900       *Pdst = clampSample(*Pdst);
901       ++Pdst;
902       --n;
903     }
904   }
905 
906   void resampleY (Sample* Pdst) {
907     Sample* Psrc;
908     ContribList* Pclist = &mPclistY[mCurDstY];
909 
910     Sample* Ptmp = mDelayXResample ? mPtmpBuf : Pdst;
911     assert(Ptmp);
912 
913     // process each contributor
914     foreach (immutable i; 0..Pclist.n) {
915       // locate the contributor's location in the scan buffer -- the contributor must always be found!
916       int j = void;
917       for (j = 0; j < MaxScanBufSize; ++j) if (mPscanBuf.scanBufY.ptr[j] == Pclist.p[i].pixel) break;
918       assert(j < MaxScanBufSize);
919       Psrc = mPscanBuf.scanBufL.ptr[j];
920       if (!i) {
921         scaleYMov(Ptmp, Psrc, Pclist.p[i].weight, mIntermediateX);
922       } else {
923         scaleYAdd(Ptmp, Psrc, Pclist.p[i].weight, mIntermediateX);
924       }
925 
926       /* If this source line doesn't contribute to any
927        * more destination lines then mark the scanline buffer slot
928        * which holds this source line as free.
929        * (The max. number of slots used depends on the Y
930        * axis sampling factor and the scaled filter width.)
931        */
932 
933       if (--mPsrcYCount[resamplerRangeCheck(Pclist.p[i].pixel, mResampleSrcY)] == 0) {
934         mPsrcYFlag[resamplerRangeCheck(Pclist.p[i].pixel, mResampleSrcY)] = false;
935         mPscanBuf.scanBufY.ptr[j] = -1;
936       }
937     }
938 
939     // now generate the destination line
940     if (mDelayXResample) {
941       // X was resampling delayed until after Y resampling
942       assert(Pdst != Ptmp);
943       resampleX(Pdst, Ptmp);
944     } else {
945       assert(Pdst == Ptmp);
946     }
947 
948     if (mLo < mHi) clamp(Pdst, mResampleDstX);
949   }
950 }
951 
952 
953 // ////////////////////////////////////////////////////////////////////////// //
954 private nothrow @trusted @nogc:
955 int resamplerRangeCheck (int v, int h) {
956   version(assert) {
957     //import std.conv : to;
958     //assert(v >= 0 && v < h, "invalid v ("~to!string(v)~"), should be in [0.."~to!string(h)~")");
959     assert(v >= 0 && v < h); // alas, @nogc
960     return v;
961   } else {
962     pragma(inline, true);
963     return v;
964   }
965 }
966 
967 enum M_PI = 3.14159265358979323846;
968 
969 // Float to int cast with truncation.
970 int castToInt (ImageResampleWorker.ResampleReal i) { pragma(inline, true); return cast(int)i; }
971 
972 // (x mod y) with special handling for negative x values.
973 int posmod (int x, int y) {
974   pragma(inline, true);
975   if (x >= 0) {
976     return (x%y);
977   } else {
978     int m = (-x)%y;
979     if (m != 0) m = y-m;
980     return m;
981   }
982 }
983 
984 // To add your own filter, insert the new function below and update the filter table.
985 // There is no need to make the filter function particularly fast, because it's
986 // only called during initializing to create the X and Y axis contributor tables.
987 
988 /* pulse/Fourier window */
989 enum BoxFilterSupport = 0.5f;
990 ImageResampleWorker.ResampleReal boxFilter (ImageResampleWorker.ResampleReal t) {
991   // make_clist() calls the filter function with t inverted (pos = left, neg = right)
992   if (t >= -0.5f && t < 0.5f) return 1.0f; else return 0.0f;
993 }
994 
995 /* box (*) box, bilinear/triangle */
996 enum TentFilterSupport = 1.0f;
997 ImageResampleWorker.ResampleReal tentFilter (ImageResampleWorker.ResampleReal t) {
998   if (t < 0.0f) t = -t;
999   if (t < 1.0f) return 1.0f-t; else return 0.0f;
1000 }
1001 
1002 /* box (*) box (*) box */
1003 enum BellSupport = 1.5f;
1004 ImageResampleWorker.ResampleReal bellFilter (ImageResampleWorker.ResampleReal t) {
1005   if (t < 0.0f) t = -t;
1006   if (t < 0.5f) return (0.75f-(t*t));
1007   if (t < 1.5f) { t = (t-1.5f); return (0.5f*(t*t)); }
1008   return (0.0f);
1009 }
1010 
1011 /* box (*) box (*) box (*) box */
1012 enum BSplineSupport = 2.0f;
1013 ImageResampleWorker.ResampleReal BSplineFilter (ImageResampleWorker.ResampleReal t) {
1014   if (t < 0.0f) t = -t;
1015   if (t < 1.0f) { immutable ImageResampleWorker.ResampleReal tt = t*t; return ((0.5f*tt*t)-tt+(2.0f/3.0f)); }
1016   if (t < 2.0f) { t = 2.0f-t; return ((1.0f/6.0f)*(t*t*t)); }
1017   return 0.0f;
1018 }
1019 
1020 // Dodgson, N., "Quadratic Interpolation for Image Resampling"
1021 enum QuadraticSupport = 1.5f;
1022 ImageResampleWorker.ResampleReal quadratic (ImageResampleWorker.ResampleReal t, in ImageResampleWorker.ResampleReal R) {
1023   pragma(inline, true);
1024   if (t < 0.0f) t = -t;
1025   if (t < QuadraticSupport) {
1026     immutable ImageResampleWorker.ResampleReal tt = t*t;
1027     if (t <= 0.5f) return (-2.0f*R)*tt+0.5f*(R+1.0f);
1028     return (R*tt)+(-2.0f*R-0.5f)*t+(3.0f/4.0f)*(R+1.0f);
1029   }
1030   return 0.0f;
1031 }
1032 
1033 ImageResampleWorker.ResampleReal quadraticInterpFilter (ImageResampleWorker.ResampleReal t) {
1034   return quadratic(t, 1.0f);
1035 }
1036 
1037 ImageResampleWorker.ResampleReal quadraticApproxFilter (ImageResampleWorker.ResampleReal t) {
1038   return quadratic(t, 0.5f);
1039 }
1040 
1041 ImageResampleWorker.ResampleReal quadraticMixFilter (ImageResampleWorker.ResampleReal t) {
1042   return quadratic(t, 0.8f);
1043 }
1044 
1045 // Mitchell, D. and A. Netravali, "Reconstruction Filters in Computer Graphics."
1046 // Computer Graphics, Vol. 22, No. 4, pp. 221-228.
1047 // (B, C)
1048 // (1/3, 1/3)  - Defaults recommended by Mitchell and Netravali
1049 // (1, 0)    - Equivalent to the Cubic B-Spline
1050 // (0, 0.5)   - Equivalent to the Catmull-Rom Spline
1051 // (0, C)   - The family of Cardinal Cubic Splines
1052 // (B, 0)   - Duff's tensioned B-Splines.
1053 ImageResampleWorker.ResampleReal mitchell (ImageResampleWorker.ResampleReal t, in ImageResampleWorker.ResampleReal B, in ImageResampleWorker.ResampleReal C) {
1054   ImageResampleWorker.ResampleReal tt = t*t;
1055   if (t < 0.0f) t = -t;
1056   if (t < 1.0f) {
1057     t = (((12.0f-9.0f*B-6.0f*C)*(t*tt))+
1058          ((-18.0f+12.0f*B+6.0f*C)*tt)+
1059          (6.0f-2.0f*B));
1060     return (t/6.0f);
1061   }
1062   if (t < 2.0f) {
1063     t = (((-1.0f*B-6.0f*C)*(t*tt))+
1064          ((6.0f*B+30.0f*C)*tt)+
1065          ((-12.0f*B-48.0f*C)*t)+
1066          (8.0f*B+24.0f*C));
1067     return (t/6.0f);
1068   }
1069   return 0.0f;
1070 }
1071 
1072 enum MitchellSupport = 2.0f;
1073 ImageResampleWorker.ResampleReal mitchellFilter (ImageResampleWorker.ResampleReal t) {
1074   return mitchell(t, 1.0f/3.0f, 1.0f/3.0f);
1075 }
1076 
1077 enum CatmullRomSupport = 2.0f;
1078 ImageResampleWorker.ResampleReal catmullRomFilter (ImageResampleWorker.ResampleReal t) {
1079   return mitchell(t, 0.0f, 0.5f);
1080 }
1081 
1082 double sinc (double x) {
1083   pragma(inline, true);
1084   import std.math : sin;
1085   x *= M_PI;
1086   if (x < 0.01f && x > -0.01f) return 1.0f+x*x*(-1.0f/6.0f+x*x*1.0f/120.0f);
1087   return sin(x)/x;
1088 }
1089 
1090 ImageResampleWorker.ResampleReal clean (double t) {
1091   pragma(inline, true);
1092   import std.math : abs;
1093   enum EPSILON = cast(ImageResampleWorker.ResampleReal)0.0000125f;
1094   if (abs(t) < EPSILON) return 0.0f;
1095   return cast(ImageResampleWorker.ResampleReal)t;
1096 }
1097 
1098 //static double blackman_window(double x)
1099 //{
1100 //  return 0.42f+0.50f*cos(M_PI*x)+0.08f*cos(2.0f*M_PI*x);
1101 //}
1102 
1103 double blackmanExactWindow (double x) {
1104   pragma(inline, true);
1105   import std.math : cos;
1106   return 0.42659071f+0.49656062f*cos(M_PI*x)+0.07684867f*cos(2.0f*M_PI*x);
1107 }
1108 
1109 enum BlackmanSupport = 3.0f;
1110 ImageResampleWorker.ResampleReal blackmanFilter (ImageResampleWorker.ResampleReal t) {
1111   if (t < 0.0f) t = -t;
1112   if (t < 3.0f) {
1113     //return clean(sinc(t)*blackman_window(t/3.0f));
1114     return clean(sinc(t)*blackmanExactWindow(t/3.0f));
1115   }
1116   return (0.0f);
1117 }
1118 
1119 // with blackman window
1120 enum GaussianSupport = 1.25f;
1121 ImageResampleWorker.ResampleReal gaussianFilter (ImageResampleWorker.ResampleReal t) {
1122   import std.math : exp, sqrt;
1123   if (t < 0) t = -t;
1124   if (t < GaussianSupport) return clean(exp(-2.0f*t*t)*sqrt(2.0f/M_PI)*blackmanExactWindow(t/GaussianSupport));
1125   return 0.0f;
1126 }
1127 
1128 // Windowed sinc -- see "Jimm Blinn's Corner: Dirty Pixels" pg. 26.
1129 enum Lanczos3Support = 3.0f;
1130 ImageResampleWorker.ResampleReal lanczos3Filter (ImageResampleWorker.ResampleReal t) {
1131   if (t < 0.0f) t = -t;
1132   if (t < 3.0f) return clean(sinc(t)*sinc(t/3.0f));
1133   return (0.0f);
1134 }
1135 
1136 enum Lanczos4Support = 4.0f;
1137 ImageResampleWorker.ResampleReal lanczos4Filter (ImageResampleWorker.ResampleReal t) {
1138   if (t < 0.0f) t = -t;
1139   if (t < 4.0f) return clean(sinc(t)*sinc(t/4.0f));
1140   return (0.0f);
1141 }
1142 
1143 enum Lanczos6Support = 6.0f;
1144 ImageResampleWorker.ResampleReal lanczos6Filter (ImageResampleWorker.ResampleReal t) {
1145   if (t < 0.0f) t = -t;
1146   if (t < 6.0f) return clean(sinc(t)*sinc(t/6.0f));
1147   return (0.0f);
1148 }
1149 
1150 enum Lanczos12Support = 12.0f;
1151 ImageResampleWorker.ResampleReal lanczos12Filter (ImageResampleWorker.ResampleReal t) {
1152   if (t < 0.0f) t = -t;
1153   if (t < 12.0f) return clean(sinc(t)*sinc(t/12.0f));
1154   return (0.0f);
1155 }
1156 
1157 double bessel0 (double x) {
1158   enum EpsilonRatio = cast(double)1E-16;
1159   double xh = 0.5*x;
1160   double sum = 1.0;
1161   double pow = 1.0;
1162   int k = 0;
1163   double ds = 1.0;
1164   // FIXME: Shouldn't this stop after X iterations for max. safety?
1165   while (ds > sum*EpsilonRatio) {
1166     ++k;
1167     pow = pow*(xh/k);
1168     ds = pow*pow;
1169     sum = sum+ds;
1170   }
1171   return sum;
1172 }
1173 
1174 enum KaiserAlpha = cast(ImageResampleWorker.ResampleReal)4.0;
1175 double kaiser (double alpha, double halfWidth, double x) {
1176   pragma(inline, true);
1177   import std.math : sqrt;
1178   immutable double ratio = (x/halfWidth);
1179   return bessel0(alpha*sqrt(1-ratio*ratio))/bessel0(alpha);
1180 }
1181 
1182 enum KaiserSupport = 3;
1183 static ImageResampleWorker.ResampleReal kaiserFilter (ImageResampleWorker.ResampleReal t) {
1184   if (t < 0.0f) t = -t;
1185   if (t < KaiserSupport) {
1186     import std.math : exp, log;
1187     // db atten
1188     immutable ImageResampleWorker.ResampleReal att = 40.0f;
1189     immutable ImageResampleWorker.ResampleReal alpha = cast(ImageResampleWorker.ResampleReal)(exp(log(cast(double)0.58417*(att-20.96))*0.4)+0.07886*(att-20.96));
1190     //const ImageResampleWorker.Resample_Real alpha = KAISER_ALPHA;
1191     return cast(ImageResampleWorker.ResampleReal)clean(sinc(t)*kaiser(alpha, KaiserSupport, t));
1192   }
1193   return 0.0f;
1194 }
1195 
1196 // filters[] is a list of all the available filter functions.
1197 struct FilterInfo {
1198   string name;
1199   ImageResampleWorker.FilterFunc func;
1200   ImageResampleWorker.ResampleReal support;
1201 }
1202 
1203 static immutable FilterInfo[16] gFilters = [
1204    FilterInfo("box",              &boxFilter,             BoxFilterSupport),
1205    FilterInfo("tent",             &tentFilter,            TentFilterSupport),
1206    FilterInfo("bell",             &bellFilter,            BellSupport),
1207    FilterInfo("bspline",          &BSplineFilter,         BSplineSupport),
1208    FilterInfo("mitchell",         &mitchellFilter,        MitchellSupport),
1209    FilterInfo("lanczos3",         &lanczos3Filter,        Lanczos3Support),
1210    FilterInfo("blackman",         &blackmanFilter,        BlackmanSupport),
1211    FilterInfo("lanczos4",         &lanczos4Filter,        Lanczos4Support),
1212    FilterInfo("lanczos6",         &lanczos6Filter,        Lanczos6Support),
1213    FilterInfo("lanczos12",        &lanczos12Filter,       Lanczos12Support),
1214    FilterInfo("kaiser",           &kaiserFilter,          KaiserSupport),
1215    FilterInfo("gaussian",         &gaussianFilter,        GaussianSupport),
1216    FilterInfo("catmullrom",       &catmullRomFilter,      CatmullRomSupport),
1217    FilterInfo("quadratic_interp", &quadraticInterpFilter, QuadraticSupport),
1218    FilterInfo("quadratic_approx", &quadraticApproxFilter, QuadraticSupport),
1219    FilterInfo("quadratic_mix",    &quadraticMixFilter,    QuadraticSupport),
1220 ];
1221 
1222 enum NumFilters = cast(int)gFilters.length;
1223 
1224 
1225 bool rsmStringEqu (const(char)[] s0, const(char)[] s1) {
1226   for (;;) {
1227     if (s0.length && (s0.ptr[0] <= ' ' || s0.ptr[0] == '_')) { s0 = s0[1..$]; continue; }
1228     if (s1.length && (s1.ptr[0] <= ' ' || s1.ptr[0] == '_')) { s1 = s1[1..$]; continue; }
1229     if (s0.length == 0) {
1230       while (s1.length && (s1.ptr[0] <= ' ' || s1.ptr[0] == '_')) s1 = s1[1..$];
1231       return (s1.length == 0);
1232     }
1233     if (s1.length == 0) {
1234       while (s0.length && (s0.ptr[0] <= ' ' || s0.ptr[0] == '_')) s0 = s0[1..$];
1235       return (s0.length == 0);
1236     }
1237     assert(s0.length && s1.length);
1238     char c0 = s0.ptr[0];
1239     char c1 = s1.ptr[0];
1240     if (c0 >= 'A' && c0 <= 'Z') c0 += 32; // poor man's tolower
1241     if (c1 >= 'A' && c1 <= 'Z') c1 += 32; // poor man's tolower
1242     if (c0 != c1) return false;
1243     s0 = s0[1..$];
1244     s1 = s1[1..$];
1245   }
1246 }
1247 
1248 
1249 int resamplerFindFilterInternal (const(char)[] name) {
1250   if (name.length) {
1251     foreach (immutable idx, const ref fi; gFilters[]) if (rsmStringEqu(name, fi.name)) return cast(int)idx;
1252   }
1253   return -1;
1254 }