2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
6 % RRRR EEEEE SSSSS AAA M M PPPP L EEEEE %
7 % R R E SS A A MM MM P P L E %
8 % RRRR EEE SSS AAAAA M M M PPPP L EEE %
9 % R R E SS A A M M P L E %
10 % R R EEEEE SSSSS A A M M P LLLLL EEEEE %
13 % MagickCore Pixel Resampling Methods %
21 % Copyright 1999-2010 ImageMagick Studio LLC, a non-profit organization %
22 % dedicated to making software imaging solutions freely available. %
24 % You may not use this file except in compliance with the License. You may %
25 % obtain a copy of the License at %
27 % http://www.imagemagick.org/script/license.php %
29 % Unless required by applicable law or agreed to in writing, software %
30 % distributed under the License is distributed on an "AS IS" BASIS, %
31 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
32 % See the License for the specific language governing permissions and %
33 % limitations under the License. %
35 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
43 #include "magick/studio.h"
44 #include "magick/artifact.h"
45 #include "magick/color-private.h"
46 #include "magick/cache.h"
47 #include "magick/draw.h"
48 #include "magick/exception-private.h"
49 #include "magick/gem.h"
50 #include "magick/image.h"
51 #include "magick/image-private.h"
52 #include "magick/log.h"
53 #include "magick/memory_.h"
54 #include "magick/pixel-private.h"
55 #include "magick/quantum.h"
56 #include "magick/random_.h"
57 #include "magick/resample.h"
58 #include "magick/resize.h"
59 #include "magick/resize-private.h"
60 #include "magick/transform.h"
61 #include "magick/signature-private.h"
65 #define WLUT_WIDTH 1024
66 struct _ResampleFilter
80 /* Information about image being resampled */
84 InterpolatePixelMethod
93 /* processing settings needed */
102 /* current ellipitical area being resampled around center point */
105 sqrtA, sqrtC, sqrtU, slope;
107 /* LUT of weights for filtered average in elliptical area */
109 filter_lut[WLUT_WIDTH],
117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
121 % A c q u i r e R e s a m p l e I n f o %
125 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
127 % AcquireResampleFilter() initializes the information resample needs do to a
128 % scaled lookup of a color from an image, using area sampling.
130 % The algorithm is based on a Elliptical Weighted Average, where the pixels
131 % found in a large elliptical area is averaged together according to a
132 % weighting (filter) function. For more details see "Fundamentals of Texture
133 % Mapping and Image Warping" a master's thesis by Paul.S.Heckbert, June 17,
134 % 1989. Available for free from, http://www.cs.cmu.edu/~ph/
136 % As EWA resampling (or any sort of resampling) can require a lot of
137 % calculations to produce a distorted scaling of the source image for each
138 % output pixel, the ResampleFilter structure generated holds that information
139 % between individual image resampling.
141 % This function will make the appropriate AcquireCacheView() calls
142 % to view the image, calling functions do not need to open a cache view.
145 % resample_filter=AcquireResampleFilter(image,exception);
146 % for (y=0; y < (long) image->rows; y++) {
147 % for (x=0; x < (long) image->columns; x++) {
149 % ScaleResampleFilter(resample_filter, ... scaling vectors ...);
150 % (void) ResamplePixelColor(resample_filter,X,Y,&pixel);
151 % ... assign resampled pixel value ...
154 % DestroyResampleFilter(resample_filter);
156 % The format of the AcquireResampleFilter method is:
158 % ResampleFilter *AcquireResampleFilter(const Image *image,
159 % ExceptionInfo *exception)
161 % A description of each parameter follows:
163 % o image: the image.
165 % o exception: return any errors or warnings in this structure.
168 MagickExport ResampleFilter *AcquireResampleFilter(const Image *image,
169 ExceptionInfo *exception)
171 register ResampleFilter
174 assert(image != (Image *) NULL);
175 assert(image->signature == MagickSignature);
176 if (image->debug != MagickFalse)
177 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
178 assert(exception != (ExceptionInfo *) NULL);
179 assert(exception->signature == MagickSignature);
181 resample_filter=(ResampleFilter *) AcquireMagickMemory(
182 sizeof(*resample_filter));
183 if (resample_filter == (ResampleFilter *) NULL)
184 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
185 (void) ResetMagickMemory(resample_filter,0,sizeof(*resample_filter));
187 resample_filter->image=ReferenceImage((Image *) image);
188 resample_filter->view=AcquireCacheView(resample_filter->image);
189 resample_filter->exception=exception;
191 resample_filter->debug=IsEventLogging();
192 resample_filter->signature=MagickSignature;
194 resample_filter->image_area = (long) resample_filter->image->columns *
195 resample_filter->image->rows;
196 resample_filter->average_defined = MagickFalse;
198 /* initialise the resampling filter settings */
199 SetResampleFilter(resample_filter, resample_filter->image->filter,
200 resample_filter->image->blur);
201 resample_filter->interpolate = resample_filter->image->interpolate;
202 resample_filter->virtual_pixel=GetImageVirtualPixelMethod(image);
204 /* init scale to a default of a unit circle */
205 ScaleResampleFilter(resample_filter, 1.0, 0.0, 0.0, 1.0);
207 return(resample_filter);
211 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
215 % D e s t r o y R e s a m p l e I n f o %
219 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
221 % DestroyResampleFilter() finalizes and cleans up the resampling
222 % resample_filter as returned by AcquireResampleFilter(), freeing any memory
223 % or other information as needed.
225 % The format of the DestroyResampleFilter method is:
227 % ResampleFilter *DestroyResampleFilter(ResampleFilter *resample_filter)
229 % A description of each parameter follows:
231 % o resample_filter: resampling information structure
234 MagickExport ResampleFilter *DestroyResampleFilter(
235 ResampleFilter *resample_filter)
237 assert(resample_filter != (ResampleFilter *) NULL);
238 assert(resample_filter->signature == MagickSignature);
239 assert(resample_filter->image != (Image *) NULL);
240 if (resample_filter->debug != MagickFalse)
241 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
242 resample_filter->image->filename);
243 resample_filter->view=DestroyCacheView(resample_filter->view);
244 resample_filter->image=DestroyImage(resample_filter->image);
245 resample_filter->signature=(~MagickSignature);
246 resample_filter=(ResampleFilter *) RelinquishMagickMemory(resample_filter);
247 return(resample_filter);
251 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
255 % I n t e r p o l a t e R e s a m p l e F i l t e r %
259 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
261 % InterpolateResampleFilter() applies bi-linear or tri-linear interpolation
262 % between a floating point coordinate and the pixels surrounding that
263 % coordinate. No pixel area resampling, or scaling of the result is
266 % The format of the InterpolateResampleFilter method is:
268 % MagickBooleanType InterpolateResampleFilter(
269 % ResampleInfo *resample_filter,const InterpolatePixelMethod method,
270 % const double x,const double y,MagickPixelPacket *pixel)
272 % A description of each parameter follows:
274 % o resample_filter: the resample filter.
276 % o method: the pixel clor interpolation method.
278 % o x,y: A double representing the current (x,y) position of the pixel.
280 % o pixel: return the interpolated pixel here.
284 static inline double MagickMax(const double x,const double y)
291 static void BicubicInterpolate(const MagickPixelPacket *pixels,const double dx,
292 MagickPixelPacket *pixel)
302 p=(pixels[3].red-pixels[2].red)-(pixels[0].red-pixels[1].red);
303 q=(pixels[0].red-pixels[1].red)-p;
304 r=pixels[2].red-pixels[0].red;
306 pixel->red=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
307 p=(pixels[3].green-pixels[2].green)-(pixels[0].green-pixels[1].green);
308 q=(pixels[0].green-pixels[1].green)-p;
309 r=pixels[2].green-pixels[0].green;
311 pixel->green=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
312 p=(pixels[3].blue-pixels[2].blue)-(pixels[0].blue-pixels[1].blue);
313 q=(pixels[0].blue-pixels[1].blue)-p;
314 r=pixels[2].blue-pixels[0].blue;
316 pixel->blue=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
317 p=(pixels[3].opacity-pixels[2].opacity)-(pixels[0].opacity-pixels[1].opacity);
318 q=(pixels[0].opacity-pixels[1].opacity)-p;
319 r=pixels[2].opacity-pixels[0].opacity;
321 pixel->opacity=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
322 if (pixel->colorspace == CMYKColorspace)
324 p=(pixels[3].index-pixels[2].index)-(pixels[0].index-pixels[1].index);
325 q=(pixels[0].index-pixels[1].index)-p;
326 r=pixels[2].index-pixels[0].index;
328 pixel->index=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
332 static inline MagickRealType CubicWeightingFunction(const MagickRealType x)
338 alpha=MagickMax(x+2.0,0.0);
339 gamma=1.0*alpha*alpha*alpha;
340 alpha=MagickMax(x+1.0,0.0);
341 gamma-=4.0*alpha*alpha*alpha;
342 alpha=MagickMax(x+0.0,0.0);
343 gamma+=6.0*alpha*alpha*alpha;
344 alpha=MagickMax(x-1.0,0.0);
345 gamma-=4.0*alpha*alpha*alpha;
349 static inline double MeshInterpolate(const PointInfo *delta,const double p,
350 const double x,const double y)
352 return(delta->x*x+delta->y*y+(1.0-delta->x-delta->y)*p);
355 static inline long NearestNeighbor(MagickRealType x)
358 return((long) (x+0.5));
359 return((long) (x-0.5));
362 static MagickBooleanType InterpolateResampleFilter(
363 ResampleFilter *resample_filter,const InterpolatePixelMethod method,
364 const double x,const double y,MagickPixelPacket *pixel)
369 register const IndexPacket
372 register const PixelPacket
378 assert(resample_filter != (ResampleFilter *) NULL);
379 assert(resample_filter->signature == MagickSignature);
383 case AverageInterpolatePixel:
392 p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x)-1,(long)
393 floor(y)-1,4,4,resample_filter->exception);
394 if (p == (const PixelPacket *) NULL)
399 indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
400 for (i=0; i < 16L; i++)
402 GetMagickPixelPacket(resample_filter->image,pixels+i);
403 SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
405 if (resample_filter->image->matte != MagickFalse)
407 alpha[i]=QuantumScale*((MagickRealType) QuantumRange-GetOpacitySample(p));
408 pixels[i].red*=alpha[i];
409 pixels[i].green*=alpha[i];
410 pixels[i].blue*=alpha[i];
411 if (resample_filter->image->colorspace == CMYKColorspace)
412 pixels[i].index*=alpha[i];
415 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
416 pixel->red+=gamma*0.0625*pixels[i].red;
417 pixel->green+=gamma*0.0625*pixels[i].green;
418 pixel->blue+=gamma*0.0625*pixels[i].blue;
419 pixel->opacity+=0.0625*pixels[i].opacity;
420 if (resample_filter->image->colorspace == CMYKColorspace)
421 pixel->index+=gamma*0.0625*pixels[i].index;
426 case BicubicInterpolatePixel:
438 p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x)-1,(long)
439 floor(y)-1,4,4,resample_filter->exception);
440 if (p == (const PixelPacket *) NULL)
445 indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
446 for (i=0; i < 16L; i++)
448 GetMagickPixelPacket(resample_filter->image,pixels+i);
449 SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
451 if (resample_filter->image->matte != MagickFalse)
453 alpha[i]=QuantumScale*((MagickRealType) QuantumRange-GetOpacitySample(p));
454 pixels[i].red*=alpha[i];
455 pixels[i].green*=alpha[i];
456 pixels[i].blue*=alpha[i];
457 if (resample_filter->image->colorspace == CMYKColorspace)
458 pixels[i].index*=alpha[i];
463 for (i=0; i < 4L; i++)
464 BicubicInterpolate(pixels+4*i,delta.x,u+i);
466 BicubicInterpolate(u,delta.y,pixel);
469 case BilinearInterpolatePixel:
483 p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x),(long)
484 floor(y),2,2,resample_filter->exception);
485 if (p == (const PixelPacket *) NULL)
490 indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
491 for (i=0; i < 4L; i++)
493 pixels[i].red=(MagickRealType) p[i].red;
494 pixels[i].green=(MagickRealType) p[i].green;
495 pixels[i].blue=(MagickRealType) p[i].blue;
496 pixels[i].opacity=(MagickRealType) p[i].opacity;
499 if (resample_filter->image->matte != MagickFalse)
500 for (i=0; i < 4L; i++)
502 alpha[i]=QuantumScale*((MagickRealType) QuantumRange-p[i].opacity);
503 pixels[i].red*=alpha[i];
504 pixels[i].green*=alpha[i];
505 pixels[i].blue*=alpha[i];
507 if (indexes != (IndexPacket *) NULL)
508 for (i=0; i < 4L; i++)
510 pixels[i].index=(MagickRealType) indexes[i];
511 if (resample_filter->image->colorspace == CMYKColorspace)
512 pixels[i].index*=alpha[i];
516 epsilon.x=1.0-delta.x;
517 epsilon.y=1.0-delta.y;
518 gamma=((epsilon.y*(epsilon.x*alpha[0]+delta.x*alpha[1])+delta.y*
519 (epsilon.x*alpha[2]+delta.x*alpha[3])));
520 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
521 pixel->red=gamma*(epsilon.y*(epsilon.x*pixels[0].red+delta.x*
522 pixels[1].red)+delta.y*(epsilon.x*pixels[2].red+delta.x*pixels[3].red));
523 pixel->green=gamma*(epsilon.y*(epsilon.x*pixels[0].green+delta.x*
524 pixels[1].green)+delta.y*(epsilon.x*pixels[2].green+delta.x*
526 pixel->blue=gamma*(epsilon.y*(epsilon.x*pixels[0].blue+delta.x*
527 pixels[1].blue)+delta.y*(epsilon.x*pixels[2].blue+delta.x*
529 pixel->opacity=(epsilon.y*(epsilon.x*pixels[0].opacity+delta.x*
530 pixels[1].opacity)+delta.y*(epsilon.x*pixels[2].opacity+delta.x*
532 if (resample_filter->image->colorspace == CMYKColorspace)
533 pixel->index=gamma*(epsilon.y*(epsilon.x*pixels[0].index+delta.x*
534 pixels[1].index)+delta.y*(epsilon.x*pixels[2].index+delta.x*
538 case FilterInterpolatePixel:
555 geometry.x=(long) floor(x)-1L;
556 geometry.y=(long) floor(y)-1L;
557 excerpt_image=ExcerptImage(resample_filter->image,&geometry,
558 resample_filter->exception);
559 if (excerpt_image == (Image *) NULL)
564 filter_image=ResizeImage(excerpt_image,1,1,resample_filter->image->filter,
565 resample_filter->image->blur,resample_filter->exception);
566 excerpt_image=DestroyImage(excerpt_image);
567 if (filter_image == (Image *) NULL)
569 filter_view=AcquireCacheView(filter_image);
570 p=GetCacheViewVirtualPixels(filter_view,0,0,1,1,
571 resample_filter->exception);
572 if (p != (const PixelPacket *) NULL)
574 indexes=GetVirtualIndexQueue(filter_image);
575 GetMagickPixelPacket(resample_filter->image,pixels);
576 SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
578 filter_view=DestroyCacheView(filter_view);
579 filter_image=DestroyImage(filter_image);
582 case IntegerInterpolatePixel:
587 p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x),(long)
588 floor(y),1,1,resample_filter->exception);
589 if (p == (const PixelPacket *) NULL)
594 indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
595 GetMagickPixelPacket(resample_filter->image,pixels);
596 SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
599 case MeshInterpolatePixel:
612 p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x),(long)
613 floor(y),2,2,resample_filter->exception);
614 if (p == (const PixelPacket *) NULL)
619 indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
620 for (i=0; i < 4L; i++)
622 GetMagickPixelPacket(resample_filter->image,pixels+i);
623 SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
625 if (resample_filter->image->matte != MagickFalse)
627 alpha[i]=QuantumScale*((MagickRealType) QuantumRange-GetOpacitySample(p));
628 pixels[i].red*=alpha[i];
629 pixels[i].green*=alpha[i];
630 pixels[i].blue*=alpha[i];
631 if (resample_filter->image->colorspace == CMYKColorspace)
632 pixels[i].index*=alpha[i];
638 luminance.x=MagickPixelLuminance(pixels+0)-MagickPixelLuminance(pixels+3);
639 luminance.y=MagickPixelLuminance(pixels+1)-MagickPixelLuminance(pixels+2);
640 if (fabs(luminance.x) < fabs(luminance.y))
645 if (delta.x <= delta.y)
648 Bottom-left triangle (pixel:2, diagonal: 0-3).
651 gamma=MeshInterpolate(&delta,alpha[2],alpha[3],alpha[0]);
652 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
653 pixel->red=gamma*MeshInterpolate(&delta,pixels[2].red,
654 pixels[3].red,pixels[0].red);
655 pixel->green=gamma*MeshInterpolate(&delta,pixels[2].green,
656 pixels[3].green,pixels[0].green);
657 pixel->blue=gamma*MeshInterpolate(&delta,pixels[2].blue,
658 pixels[3].blue,pixels[0].blue);
659 pixel->opacity=gamma*MeshInterpolate(&delta,pixels[2].opacity,
660 pixels[3].opacity,pixels[0].opacity);
661 if (resample_filter->image->colorspace == CMYKColorspace)
662 pixel->index=gamma*MeshInterpolate(&delta,pixels[2].index,
663 pixels[3].index,pixels[0].index);
668 Top-right triangle (pixel:1, diagonal: 0-3).
671 gamma=MeshInterpolate(&delta,alpha[1],alpha[0],alpha[3]);
672 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
673 pixel->red=gamma*MeshInterpolate(&delta,pixels[1].red,
674 pixels[0].red,pixels[3].red);
675 pixel->green=gamma*MeshInterpolate(&delta,pixels[1].green,
676 pixels[0].green,pixels[3].green);
677 pixel->blue=gamma*MeshInterpolate(&delta,pixels[1].blue,
678 pixels[0].blue,pixels[3].blue);
679 pixel->opacity=gamma*MeshInterpolate(&delta,pixels[1].opacity,
680 pixels[0].opacity,pixels[3].opacity);
681 if (resample_filter->image->colorspace == CMYKColorspace)
682 pixel->index=gamma*MeshInterpolate(&delta,pixels[1].index,
683 pixels[0].index,pixels[3].index);
691 if (delta.x <= (1.0-delta.y))
694 Top-left triangle (pixel 0, diagonal: 1-2).
696 gamma=MeshInterpolate(&delta,alpha[0],alpha[1],alpha[2]);
697 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
698 pixel->red=gamma*MeshInterpolate(&delta,pixels[0].red,
699 pixels[1].red,pixels[2].red);
700 pixel->green=gamma*MeshInterpolate(&delta,pixels[0].green,
701 pixels[1].green,pixels[2].green);
702 pixel->blue=gamma*MeshInterpolate(&delta,pixels[0].blue,
703 pixels[1].blue,pixels[2].blue);
704 pixel->opacity=gamma*MeshInterpolate(&delta,pixels[0].opacity,
705 pixels[1].opacity,pixels[2].opacity);
706 if (resample_filter->image->colorspace == CMYKColorspace)
707 pixel->index=gamma*MeshInterpolate(&delta,pixels[0].index,
708 pixels[1].index,pixels[2].index);
713 Bottom-right triangle (pixel: 3, diagonal: 1-2).
717 gamma=MeshInterpolate(&delta,alpha[3],alpha[2],alpha[1]);
718 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
719 pixel->red=gamma*MeshInterpolate(&delta,pixels[3].red,
720 pixels[2].red,pixels[1].red);
721 pixel->green=gamma*MeshInterpolate(&delta,pixels[3].green,
722 pixels[2].green,pixels[1].green);
723 pixel->blue=gamma*MeshInterpolate(&delta,pixels[3].blue,
724 pixels[2].blue,pixels[1].blue);
725 pixel->opacity=gamma*MeshInterpolate(&delta,pixels[3].opacity,
726 pixels[2].opacity,pixels[1].opacity);
727 if (resample_filter->image->colorspace == CMYKColorspace)
728 pixel->index=gamma*MeshInterpolate(&delta,pixels[3].index,
729 pixels[2].index,pixels[1].index);
734 case NearestNeighborInterpolatePixel:
739 p=GetCacheViewVirtualPixels(resample_filter->view,NearestNeighbor(x),
740 NearestNeighbor(y),1,1,resample_filter->exception);
741 if (p == (const PixelPacket *) NULL)
746 indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
747 GetMagickPixelPacket(resample_filter->image,pixels);
748 SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
751 case SplineInterpolatePixel:
769 p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x)-1,(long)
770 floor(y)-1,4,4,resample_filter->exception);
771 if (p == (const PixelPacket *) NULL)
776 indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
780 for (i=(-1); i < 3L; i++)
782 dy=CubicWeightingFunction((MagickRealType) i-delta.y);
783 for (j=(-1); j < 3L; j++)
785 GetMagickPixelPacket(resample_filter->image,pixels+n);
786 SetMagickPixelPacket(resample_filter->image,p,indexes+n,pixels+n);
788 if (resample_filter->image->matte != MagickFalse)
790 alpha[n]=QuantumScale*((MagickRealType) QuantumRange-GetOpacitySample(p));
791 pixels[n].red*=alpha[n];
792 pixels[n].green*=alpha[n];
793 pixels[n].blue*=alpha[n];
794 if (resample_filter->image->colorspace == CMYKColorspace)
795 pixels[n].index*=alpha[n];
797 dx=CubicWeightingFunction(delta.x-(MagickRealType) j);
799 gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
800 pixel->red+=gamma*dx*dy*pixels[n].red;
801 pixel->green+=gamma*dx*dy*pixels[n].green;
802 pixel->blue+=gamma*dx*dy*pixels[n].blue;
803 if (resample_filter->image->matte != MagickFalse)
804 pixel->opacity+=dx*dy*pixels[n].opacity;
805 if (resample_filter->image->colorspace == CMYKColorspace)
806 pixel->index+=gamma*dx*dy*pixels[n].index;
818 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
822 % R e s a m p l e P i x e l C o l o r %
826 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
828 % ResamplePixelColor() samples the pixel values surrounding the location
829 % given using an elliptical weighted average, at the scale previously
830 % calculated, and in the most efficent manner possible for the
831 % VirtualPixelMethod setting.
833 % The format of the ResamplePixelColor method is:
835 % MagickBooleanType ResamplePixelColor(ResampleFilter *resample_filter,
836 % const double u0,const double v0,MagickPixelPacket *pixel)
838 % A description of each parameter follows:
840 % o resample_filter: the resample filter.
842 % o u0,v0: A double representing the center of the area to resample,
843 % The distortion transformed transformed x,y coordinate.
845 % o pixel: the resampled pixel is returned here.
848 MagickExport MagickBooleanType ResamplePixelColor(
849 ResampleFilter *resample_filter,const double u0,const double v0,
850 MagickPixelPacket *pixel)
855 long u,v, uw,v1,v2, hit;
858 double divisor_c,divisor_m;
859 register double weight;
860 register const PixelPacket *pixels;
861 register const IndexPacket *indexes;
862 assert(resample_filter != (ResampleFilter *) NULL);
863 assert(resample_filter->signature == MagickSignature);
866 GetMagickPixelPacket(resample_filter->image,pixel);
867 if ( resample_filter->do_interpolate ) {
868 status=InterpolateResampleFilter(resample_filter,
869 resample_filter->interpolate,u0,v0,pixel);
874 Does resample area Miss the image?
875 And is that area a simple solid color - then return that color
878 switch ( resample_filter->virtual_pixel ) {
879 case BackgroundVirtualPixelMethod:
880 case ConstantVirtualPixelMethod:
881 case TransparentVirtualPixelMethod:
882 case BlackVirtualPixelMethod:
883 case GrayVirtualPixelMethod:
884 case WhiteVirtualPixelMethod:
885 case MaskVirtualPixelMethod:
886 if ( resample_filter->limit_reached
887 || u0 + resample_filter->sqrtC < 0.0
888 || u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
889 || v0 + resample_filter->sqrtA < 0.0
890 || v0 - resample_filter->sqrtA > (double) resample_filter->image->rows
895 case UndefinedVirtualPixelMethod:
896 case EdgeVirtualPixelMethod:
897 if ( ( u0 + resample_filter->sqrtC < 0.0 && v0 + resample_filter->sqrtA < 0.0 )
898 || ( u0 + resample_filter->sqrtC < 0.0
899 && v0 - resample_filter->sqrtA > (double) resample_filter->image->rows )
900 || ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
901 && v0 + resample_filter->sqrtA < 0.0 )
902 || ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
903 && v0 - resample_filter->sqrtA > (double) resample_filter->image->rows )
907 case HorizontalTileVirtualPixelMethod:
908 if ( v0 + resample_filter->sqrtA < 0.0
909 || v0 - resample_filter->sqrtA > (double) resample_filter->image->rows
911 hit++; /* outside the horizontally tiled images. */
913 case VerticalTileVirtualPixelMethod:
914 if ( u0 + resample_filter->sqrtC < 0.0
915 || u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
917 hit++; /* outside the vertically tiled images. */
919 case DitherVirtualPixelMethod:
920 if ( ( u0 + resample_filter->sqrtC < -32.0 && v0 + resample_filter->sqrtA < -32.0 )
921 || ( u0 + resample_filter->sqrtC < -32.0
922 && v0 - resample_filter->sqrtA > (double) resample_filter->image->rows+32.0 )
923 || ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns+32.0
924 && v0 + resample_filter->sqrtA < -32.0 )
925 || ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns+32.0
926 && v0 - resample_filter->sqrtA > (double) resample_filter->image->rows+32.0 )
930 case TileVirtualPixelMethod:
931 case MirrorVirtualPixelMethod:
932 case RandomVirtualPixelMethod:
933 case HorizontalTileEdgeVirtualPixelMethod:
934 case VerticalTileEdgeVirtualPixelMethod:
935 case CheckerTileVirtualPixelMethod:
936 /* resampling of area is always needed - no VP limits */
940 /* whole area is a solid color -- just return that color */
941 status=InterpolateResampleFilter(resample_filter,IntegerInterpolatePixel,
947 Scaling limits reached, return an 'averaged' result.
949 if ( resample_filter->limit_reached ) {
950 switch ( resample_filter->virtual_pixel ) {
951 /* This is always handled by the above, so no need.
952 case BackgroundVirtualPixelMethod:
953 case ConstantVirtualPixelMethod:
954 case TransparentVirtualPixelMethod:
955 case GrayVirtualPixelMethod,
956 case WhiteVirtualPixelMethod
957 case MaskVirtualPixelMethod:
959 case UndefinedVirtualPixelMethod:
960 case EdgeVirtualPixelMethod:
961 case DitherVirtualPixelMethod:
962 case HorizontalTileEdgeVirtualPixelMethod:
963 case VerticalTileEdgeVirtualPixelMethod:
964 /* We need an average edge pixel, for the right edge!
965 How should I calculate an average edge color?
966 Just returning an averaged neighbourhood,
967 works well in general, but falls down for TileEdge methods.
968 This needs to be done properly!!!!!!
970 status=InterpolateResampleFilter(resample_filter,
971 AverageInterpolatePixel,u0,v0,pixel);
973 case HorizontalTileVirtualPixelMethod:
974 case VerticalTileVirtualPixelMethod:
975 /* just return the background pixel - Is there more direct way? */
976 status=InterpolateResampleFilter(resample_filter,
977 IntegerInterpolatePixel,(double)-1,(double)-1,pixel);
979 case TileVirtualPixelMethod:
980 case MirrorVirtualPixelMethod:
981 case RandomVirtualPixelMethod:
982 case CheckerTileVirtualPixelMethod:
984 /* generate a average color of the WHOLE image */
985 if ( resample_filter->average_defined == MagickFalse ) {
992 GetMagickPixelPacket(resample_filter->image,
993 (MagickPixelPacket *)&(resample_filter->average_pixel));
994 resample_filter->average_defined = MagickTrue;
996 /* Try to get an averaged pixel color of whole image */
997 average_image=ResizeImage(resample_filter->image,1,1,BoxFilter,1.0,
998 resample_filter->exception);
999 if (average_image == (Image *) NULL)
1001 *pixel=resample_filter->average_pixel; /* FAILED */
1004 average_view=AcquireCacheView(average_image);
1005 pixels=(PixelPacket *)GetCacheViewVirtualPixels(average_view,0,0,1,1,
1006 resample_filter->exception);
1007 if (pixels == (const PixelPacket *) NULL) {
1008 average_view=DestroyCacheView(average_view);
1009 average_image=DestroyImage(average_image);
1010 *pixel=resample_filter->average_pixel; /* FAILED */
1013 indexes=(IndexPacket *) GetCacheViewAuthenticIndexQueue(average_view);
1014 SetMagickPixelPacket(resample_filter->image,pixels,indexes,
1015 &(resample_filter->average_pixel));
1016 average_view=DestroyCacheView(average_view);
1017 average_image=DestroyImage(average_image);
1019 /* CheckerTile should average the image with background color */
1020 //if ( resample_filter->virtual_pixel == CheckerTileVirtualPixelMethod ) {
1022 resample_filter->average_pixel.red =
1023 ( resample_filter->average_pixel.red +
1024 resample_filter->image->background_color.red ) /2;
1025 resample_filter->average_pixel.green =
1026 ( resample_filter->average_pixel.green +
1027 resample_filter->image->background_color.green ) /2;
1028 resample_filter->average_pixel.blue =
1029 ( resample_filter->average_pixel.blue +
1030 resample_filter->image->background_color.blue ) /2;
1031 resample_filter->average_pixel.matte =
1032 ( resample_filter->average_pixel.matte +
1033 resample_filter->image->background_color.matte ) /2;
1034 resample_filter->average_pixel.black =
1035 ( resample_filter->average_pixel.black +
1036 resample_filter->image->background_color.black ) /2;
1038 resample_filter->average_pixel =
1039 resample_filter->image->background_color;
1044 *pixel=resample_filter->average_pixel;
1051 Initialize weighted average data collection
1056 pixel->red = pixel->green = pixel->blue = 0.0;
1057 if (resample_filter->image->matte != MagickFalse) pixel->opacity = 0.0;
1058 if (resample_filter->image->colorspace == CMYKColorspace) pixel->index = 0.0;
1061 Determine the parellelogram bounding box fitted to the ellipse
1062 centered at u0,v0. This area is bounding by the lines...
1064 u = -By/2A +/- sqrt(F/A)
1065 Which has been pre-calculated above.
1067 v1 = (long)(v0 - resample_filter->sqrtA); /* range of scan lines */
1068 v2 = (long)(v0 + resample_filter->sqrtA + 1);
1070 u1 = u0 + (v1-v0)*resample_filter->slope - resample_filter->sqrtU; /* start of scanline for v=v1 */
1071 uw = (long)(2*resample_filter->sqrtU)+1; /* width of parallelogram */
1074 Do weighted resampling of all pixels, within the scaled ellipse,
1075 bound by a Parellelogram fitted to the ellipse.
1077 DDQ = 2*resample_filter->A;
1078 for( v=v1; v<=v2; v++, u1+=resample_filter->slope ) {
1079 u = (long)u1; /* first pixel in scanline ( floor(u1) ) */
1080 U = (double)u-u0; /* location of that pixel, relative to u0,v0 */
1083 /* Q = ellipse quotent ( if Q<F then pixel is inside ellipse) */
1084 Q = U*(resample_filter->A*U + resample_filter->B*V) + resample_filter->C*V*V;
1085 DQ = resample_filter->A*(2.0*U+1) + resample_filter->B*V;
1087 /* get the scanline of pixels for this v */
1088 pixels=GetCacheViewVirtualPixels(resample_filter->view,u,v,(unsigned long) uw,
1089 1,resample_filter->exception);
1090 if (pixels == (const PixelPacket *) NULL)
1091 return(MagickFalse);
1092 indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
1094 /* count up the weighted pixel colors */
1095 for( u=0; u<uw; u++ ) {
1096 /* Note that the ellipse has been pre-scaled so F = WLUT_WIDTH */
1097 if ( Q < (double)WLUT_WIDTH ) {
1098 weight = resample_filter->filter_lut[(int)Q];
1100 pixel->opacity += weight*pixels->opacity;
1101 divisor_m += weight;
1103 if (resample_filter->image->matte != MagickFalse)
1104 weight *= QuantumScale*((MagickRealType)(QuantumRange-pixels->opacity));
1105 pixel->red += weight*pixels->red;
1106 pixel->green += weight*pixels->green;
1107 pixel->blue += weight*pixels->blue;
1108 if (resample_filter->image->colorspace == CMYKColorspace)
1109 pixel->index += weight*(*indexes);
1110 divisor_c += weight;
1122 Result sanity check -- this should NOT happen
1124 if ( hit < 4 || divisor_c < 1.0 ) {
1125 /* not enough pixels in resampling, resort to direct interpolation */
1126 status=InterpolateResampleFilter(resample_filter,
1127 resample_filter->interpolate,u0,v0,pixel);
1132 Finialize results of resampling
1134 divisor_m = 1.0/divisor_m;
1135 pixel->opacity = (MagickRealType) RoundToQuantum(divisor_m*pixel->opacity);
1136 divisor_c = 1.0/divisor_c;
1137 pixel->red = (MagickRealType) RoundToQuantum(divisor_c*pixel->red);
1138 pixel->green = (MagickRealType) RoundToQuantum(divisor_c*pixel->green);
1139 pixel->blue = (MagickRealType) RoundToQuantum(divisor_c*pixel->blue);
1140 if (resample_filter->image->colorspace == CMYKColorspace)
1141 pixel->index = (MagickRealType) RoundToQuantum(divisor_c*pixel->index);
1146 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1150 % S c a l e R e s a m p l e F i l t e r %
1154 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1156 % ScaleResampleFilter() does all the calculations needed to resample an image
1157 % at a specific scale, defined by two scaling vectors. This not using
1158 % a orthogonal scaling, but two distorted scaling vectors, to allow the
1159 % generation of a angled ellipse.
1161 % As only two deritive scaling vectors are used the center of the ellipse
1162 % must be the center of the lookup. That is any curvature that the
1163 % distortion may produce is discounted.
1165 % The input vectors are produced by either finding the derivitives of the
1166 % distortion function, or the partial derivitives from a distortion mapping.
1167 % They do not need to be the orthogonal dx,dy scaling vectors, but can be
1168 % calculated from other derivatives. For example you could use dr,da/r
1169 % polar coordinate vector scaling vectors
1171 % If u,v = DistortEquation(x,y)
1172 % Then the scaling vectors dx,dy (in u,v space) are the derivitives...
1173 % du/dx, dv/dx and du/dy, dv/dy
1174 % If the scaling is only othogonally aligned then...
1175 % dv/dx = 0 and du/dy = 0
1176 % Producing an othogonally alligned ellipse for the area to be resampled.
1178 % Note that scaling vectors are different to argument order. Argument order
1179 % is the general order the deritives are extracted from the distortion
1180 % equations, EG: U(x,y), V(x,y). Caution is advised if you are trying to
1181 % define the ellipse directly from scaling vectors.
1183 % The format of the ScaleResampleFilter method is:
1185 % void ScaleResampleFilter(const ResampleFilter *resample_filter,
1186 % const double dux,const double duy,const double dvx,const double dvy)
1188 % A description of each parameter follows:
1190 % o resample_filter: the resampling resample_filterrmation defining the
1191 % image being resampled
1193 % o dux,duy,dvx,dvy:
1194 % The partial derivitives or scaling vectors for resampling.
1195 % dx = du/dx, dv/dx and dy = du/dy, dv/dy
1197 % The values are used to define the size and angle of the
1198 % elliptical resampling area, centered on the lookup point.
1201 MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
1202 const double dux,const double duy,const double dvx,const double dvy)
1204 double A,B,C,F, area;
1206 assert(resample_filter != (ResampleFilter *) NULL);
1207 assert(resample_filter->signature == MagickSignature);
1209 resample_filter->limit_reached = MagickFalse;
1210 resample_filter->do_interpolate = MagickFalse;
1212 /* A 'point' filter forces use of interpolation instead of area sampling */
1213 if ( resample_filter->filter == PointFilter ) {
1214 resample_filter->do_interpolate = MagickTrue;
1218 /* Find Ellipse Coefficents such that
1219 A*u^2 + B*u*v + C*v^2 = F
1220 With u,v relative to point around which we are resampling.
1221 And the given scaling dx,dy vectors in u,v space
1222 du/dx,dv/dx and du/dy,dv/dy
1225 /* Direct conversions of derivatives to elliptical coefficients
1226 No scaling will result in F == 1.0 and a unit circle.
1228 A = dvx*dvx+dvy*dvy;
1229 B = (dux*dvx+duy*dvy)*-2.0;
1230 C = dux*dux+duy*duy;
1231 F = dux*dvy+duy*dvx;
1235 /* This Paul Heckbert's recomended "Higher Quality EWA" formula, from page
1236 60 in his thesis, which adds a unit circle to the elliptical area so are
1237 to do both Reconstruction and Prefiltering of the pixels in the
1238 resampling. It also means it is likely to have at least 4 pixels within
1239 the area of the ellipse, for weighted averaging.
1240 No scaling will result if F == 4.0 and a circle of radius 2.0
1242 A = dvx*dvx+dvy*dvy+1;
1243 B = (dux*dvx+duy*dvy)*-2.0;
1244 C = dux*dux+duy*duy+1;
1249 /* DEBUGGING OUTPUT */
1251 fprintf(stderr, "dux=%lf; dvx=%lf; duy=%lf; dvy%lf;\n",
1252 dux, dvx, duy, dvy);
1253 fprintf(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
1257 /* Figure out the Ellipses Major and Minor Axis, and other info.
1258 This information currently not needed at this time, but may be
1259 needed later for better limit determination.
1261 { double alpha, beta, gamma, Major, Minor;
1262 double Eccentricity, Ellipse_Area, Ellipse_angle;
1263 double max_horizontal_cross_section, max_vertical_cross_section;
1266 gamma = sqrt(beta*beta + B*B );
1268 if ( alpha - gamma <= MagickEpsilon )
1271 Major = sqrt(2*F/(alpha - gamma));
1272 Minor = sqrt(2*F/(alpha + gamma));
1274 fprintf(stderr, "\tMajor=%lf; Minor=%lf\n",
1277 /* other information about ellipse include... */
1278 Eccentricity = Major/Minor;
1279 Ellipse_Area = MagickPI*Major*Minor;
1280 Ellipse_angle = atan2(B, A-C);
1282 fprintf(stderr, "\tAngle=%lf Area=%lf\n",
1283 RadiansToDegrees(Ellipse_angle), Ellipse_Area );
1285 /* Ellipse limits */
1287 /* orthogonal rectangle - improved ellipse */
1288 max_horizontal_orthogonal = sqrt(A); /* = sqrt(4*A*F/(4*A*C-B*B)) */
1289 max_vertical_orthogonal = sqrt(C); /* = sqrt(4*C*F/(4*A*C-B*B)) */
1291 /* parallelogram bounds -- what we are using */
1292 max_horizontal_cross_section = sqrt(F/A);
1293 max_vertical_cross_section = sqrt(F/C);
1297 /* Is default elliptical area, too small? Image being magnified?
1298 Switch to doing pure 'point' interpolation of the pixel.
1299 That is turn off EWA Resampling.
1301 if ( F <= F_UNITY ) {
1302 resample_filter->do_interpolate = MagickTrue;
1307 /* If F is impossibly large, we may as well not bother doing any
1308 * form of resampling, as you risk an infinite resampled area.
1310 if ( F > MagickHuge ) {
1311 resample_filter->limit_reached = MagickTrue;
1315 /* Othogonal bounds of the ellipse */
1316 resample_filter->sqrtA = sqrt(A)+1.0; /* Vertical Orthogonal Limit */
1317 resample_filter->sqrtC = sqrt(C)+1.0; /* Horizontal Orthogonal Limit */
1319 /* Horizontally aligned Parallelogram fitted to ellipse */
1320 resample_filter->sqrtU = sqrt(F/A)+1.0; /* Parallelogram Width */
1321 resample_filter->slope = -B/(2*A); /* Slope of the parallelogram */
1323 /* The size of the area of the parallelogram we will be sampling */
1324 area = 4 * resample_filter->sqrtA * resample_filter->sqrtU;
1326 /* Absolute limit on the area to be resampled
1327 * This limit needs more work, as it gets too slow for
1328 * larger images involved with tiled views of the horizon. */
1329 if ( area > 20.0*resample_filter->image_area ) {
1330 resample_filter->limit_reached = MagickTrue;
1334 /* Scale ellipse formula to directly fit the Filter Lookup Table */
1335 { register double scale;
1336 scale = (double)WLUT_WIDTH/F;
1337 resample_filter->A = A*scale;
1338 resample_filter->B = B*scale;
1339 resample_filter->C = C*scale;
1340 /* ..ple_filter->F = WLUT_WIDTH; -- hardcoded */
1345 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1349 % S e t R e s a m p l e F i l t e r %
1353 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1355 % SetResampleFilter() set the resampling filter lookup table based on a
1356 % specific filter. Note that the filter is used as a radial filter not as a
1357 % two pass othogonally aligned resampling filter.
1359 % The default Filter, is Gaussian, which is the standard filter used by the
1360 % original paper on the Elliptical Weighted Everage Algorithm. However other
1361 % filters can also be used.
1363 % The format of the SetResampleFilter method is:
1365 % void SetResampleFilter(ResampleFilter *resample_filter,
1366 % const FilterTypes filter,const double blur)
1368 % A description of each parameter follows:
1370 % o resample_filter: resampling resample_filterrmation structure
1372 % o filter: the resize filter for elliptical weighting LUT
1374 % o blur: filter blur factor (radial scaling) for elliptical weighting LUT
1377 MagickExport void SetResampleFilter(ResampleFilter *resample_filter,
1378 const FilterTypes filter,const double blur)
1389 assert(resample_filter != (ResampleFilter *) NULL);
1390 assert(resample_filter->signature == MagickSignature);
1392 resample_filter->filter = filter;
1394 /* Scale radius so it equals 1.0, at edge of ellipse when a
1395 default blurring factor of 1.0 is used.
1397 Note that these filters are being used as a radial filter, not as
1398 an othoginally alligned filter. How this effects results is still
1401 Future: Direct use of teh resize filters in "resize.c" to set the lookup
1402 table, based on the filters working support window.
1404 r_scale = sqrt(1.0/(double)WLUT_WIDTH)/blur;
1405 r_scale *= 2; /* for 2 pixel radius of Improved Elliptical Formula */
1409 /* This equivelent to turning off the EWA algroithm.
1410 Only Interpolated lookup will be used. */
1414 Fill the LUT with a 1D resize filter function
1415 But make the Sinc/Bessel tapered window 2.0
1416 I also normalize the result so the filter is 1.0
1418 resize_filter = AcquireResizeFilter(resample_filter->image,filter,
1419 (MagickRealType)1.0,MagickTrue,resample_filter->exception);
1420 if (resize_filter != (ResizeFilter *) NULL) {
1421 resample_filter->support = GetResizeFilterSupport(resize_filter);
1422 resample_filter->support /= blur; /* taken into account above */
1423 resample_filter->support *= resample_filter->support;
1424 resample_filter->support *= (double)WLUT_WIDTH/4;
1425 if ( resample_filter->support >= (double)WLUT_WIDTH )
1426 resample_filter->support = (double)WLUT_WIDTH; /* hack */
1427 for(Q=0; Q<WLUT_WIDTH; Q++)
1428 if ( (double) Q < resample_filter->support )
1429 resample_filter->filter_lut[Q] = (double)
1430 GetResizeFilterWeight(resize_filter,sqrt((double)Q)*r_scale);
1432 resample_filter->filter_lut[Q] = 0.0;
1433 resize_filter = DestroyResizeFilter(resize_filter);
1437 (void) ThrowMagickException(resample_filter->exception,GetMagickModule(),
1438 ModuleError, "UnableToSetFilteringValue",
1439 "Fall back to default EWA gaussian filter");
1441 /* FALLTHRU - on exception */
1442 /*case GaussianFilter:*/
1443 case UndefinedFilter:
1445 Create Normal Gaussian 2D Filter Weighted Lookup Table.
1446 A normal EWA guassual lookup would use exp(Q*ALPHA)
1447 where Q = distantce squared from 0.0 (center) to 1.0 (edge)
1448 and ALPHA = -4.0*ln(2.0) ==> -2.77258872223978123767
1449 However the table is of length 1024, and equates to a radius of 2px
1450 thus needs to be scaled by ALPHA*4/1024 and any blur factor squared
1452 /*r_scale = -2.77258872223978123767*4/WLUT_WIDTH/blur/blur;*/
1453 r_scale = -2.77258872223978123767/WLUT_WIDTH/blur/blur;
1454 for(Q=0; Q<WLUT_WIDTH; Q++)
1455 resample_filter->filter_lut[Q] = exp((double)Q*r_scale);
1456 resample_filter->support = WLUT_WIDTH;
1459 if (GetImageArtifact(resample_filter->image,"resample:verbose")
1460 != (const char *) NULL)
1461 /* Debug output of the filter weighting LUT
1462 Gnuplot the LUT with hoizontal adjusted to 'r' using...
1463 plot [0:2][-.2:1] "lut.dat" using (sqrt($0/1024)*2):1 with lines
1464 THe filter values is normalized for comparision
1466 for(Q=0; Q<WLUT_WIDTH; Q++)
1467 printf("%lf\n", resample_filter->filter_lut[Q]
1468 /resample_filter->filter_lut[0] );
1473 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1477 % S e t R e s a m p l e F i l t e r I n t e r p o l a t e M e t h o d %
1481 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1483 % SetResampleFilterInterpolateMethod() changes the interpolation method
1484 % associated with the specified resample filter.
1486 % The format of the SetResampleFilterInterpolateMethod method is:
1488 % MagickBooleanType SetResampleFilterInterpolateMethod(
1489 % ResampleFilter *resample_filter,const InterpolateMethod method)
1491 % A description of each parameter follows:
1493 % o resample_filter: the resample filter.
1495 % o method: the interpolation method.
1498 MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
1499 ResampleFilter *resample_filter,const InterpolatePixelMethod method)
1501 assert(resample_filter != (ResampleFilter *) NULL);
1502 assert(resample_filter->signature == MagickSignature);
1503 assert(resample_filter->image != (Image *) NULL);
1504 if (resample_filter->debug != MagickFalse)
1505 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1506 resample_filter->image->filename);
1507 resample_filter->interpolate=method;
1512 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1516 % S e t R e s a m p l e F i l t e r V i r t u a l P i x e l M e t h o d %
1520 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1522 % SetResampleFilterVirtualPixelMethod() changes the virtual pixel method
1523 % associated with the specified resample filter.
1525 % The format of the SetResampleFilterVirtualPixelMethod method is:
1527 % MagickBooleanType SetResampleFilterVirtualPixelMethod(
1528 % ResampleFilter *resample_filter,const VirtualPixelMethod method)
1530 % A description of each parameter follows:
1532 % o resample_filter: the resample filter.
1534 % o method: the virtual pixel method.
1537 MagickExport MagickBooleanType SetResampleFilterVirtualPixelMethod(
1538 ResampleFilter *resample_filter,const VirtualPixelMethod method)
1540 assert(resample_filter != (ResampleFilter *) NULL);
1541 assert(resample_filter->signature == MagickSignature);
1542 assert(resample_filter->image != (Image *) NULL);
1543 if (resample_filter->debug != MagickFalse)
1544 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1545 resample_filter->image->filename);
1546 resample_filter->virtual_pixel=method;
1547 (void) SetCacheViewVirtualPixelMethod(resample_filter->view,method);