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