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-2011 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/magick.h"
54 #include "magick/memory_.h"
55 #include "magick/pixel.h"
56 #include "magick/pixel-private.h"
57 #include "magick/quantum.h"
58 #include "magick/random_.h"
59 #include "magick/resample.h"
60 #include "magick/resize.h"
61 #include "magick/resize-private.h"
62 #include "magick/transform.h"
63 #include "magick/signature-private.h"
64 #include "magick/utility.h"
66 EWA Resampling Options
69 /* select ONE resampling method */
70 #define EWA 1 /* Normal EWA handling - raw or clamped */
71 /* if 0 then use "High Quality EWA" */
72 #define EWA_CLAMP 1 /* EWA Clamping from Nicolas Robidoux */
74 #define FILTER_LUT 1 /* Use a LUT rather then direct filter calls */
76 /* output debugging information */
77 #define DEBUG_ELLIPSE 0 /* output ellipse info for debug */
78 #define DEBUG_HIT_MISS 0 /* output hit/miss pixels (as gnuplot commands) */
79 #define DEBUG_NO_PIXEL_HIT 0 /* Make pixels that fail to hit anything - RED */
82 #define WLUT_WIDTH 1024 /* size of the filter cache */
88 struct _ResampleFilter
102 /* Information about image being resampled */
106 InterpolatePixelMethod
115 /* processing settings needed */
124 /* current ellipitical area being resampled around center point */
127 Vlimit, Ulimit, Uwidth, slope;
130 /* LUT of weights for filtered average in elliptical area */
132 filter_lut[WLUT_WIDTH];
134 /* Use a Direct call to the filter functions */
142 /* the practical working support of the filter */
151 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
155 % A c q u i r e R e s a m p l e I n f o %
159 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
161 % AcquireResampleFilter() initializes the information resample needs do to a
162 % scaled lookup of a color from an image, using area sampling.
164 % The algorithm is based on a Elliptical Weighted Average, where the pixels
165 % found in a large elliptical area is averaged together according to a
166 % weighting (filter) function. For more details see "Fundamentals of Texture
167 % Mapping and Image Warping" a master's thesis by Paul.S.Heckbert, June 17,
168 % 1989. Available for free from, http://www.cs.cmu.edu/~ph/
170 % As EWA resampling (or any sort of resampling) can require a lot of
171 % calculations to produce a distorted scaling of the source image for each
172 % output pixel, the ResampleFilter structure generated holds that information
173 % between individual image resampling.
175 % This function will make the appropriate AcquireCacheView() calls
176 % to view the image, calling functions do not need to open a cache view.
179 % resample_filter=AcquireResampleFilter(image,exception);
180 % SetResampleFilter(resample_filter, GaussianFilter, 1.0);
181 % for (y=0; y < (ssize_t) image->rows; y++) {
182 % for (x=0; x < (ssize_t) image->columns; x++) {
184 % ScaleResampleFilter(resample_filter, ... scaling vectors ...);
185 % (void) ResamplePixelColor(resample_filter,u,v,&pixel);
186 % ... assign resampled pixel value ...
189 % DestroyResampleFilter(resample_filter);
191 % The format of the AcquireResampleFilter method is:
193 % ResampleFilter *AcquireResampleFilter(const Image *image,
194 % ExceptionInfo *exception)
196 % A description of each parameter follows:
198 % o image: the image.
200 % o exception: return any errors or warnings in this structure.
203 MagickExport ResampleFilter *AcquireResampleFilter(const Image *image,
204 ExceptionInfo *exception)
206 register ResampleFilter
209 assert(image != (Image *) NULL);
210 assert(image->signature == MagickSignature);
211 if (image->debug != MagickFalse)
212 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
213 assert(exception != (ExceptionInfo *) NULL);
214 assert(exception->signature == MagickSignature);
216 resample_filter=(ResampleFilter *) AcquireMagickMemory(
217 sizeof(*resample_filter));
218 if (resample_filter == (ResampleFilter *) NULL)
219 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
220 (void) ResetMagickMemory(resample_filter,0,sizeof(*resample_filter));
222 resample_filter->exception=exception;
223 resample_filter->image=ReferenceImage((Image *) image);
224 resample_filter->view=AcquireCacheView(resample_filter->image);
226 resample_filter->debug=IsEventLogging();
227 resample_filter->signature=MagickSignature;
229 resample_filter->image_area=(ssize_t) (image->columns*image->rows);
230 resample_filter->average_defined = MagickFalse;
232 /* initialise the resampling filter settings */
233 SetResampleFilter(resample_filter, image->filter, image->blur);
234 (void) SetResampleFilterInterpolateMethod(resample_filter,
236 (void) SetResampleFilterVirtualPixelMethod(resample_filter,
237 GetImageVirtualPixelMethod(image));
239 return(resample_filter);
243 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
247 % D e s t r o y R e s a m p l e I n f o %
251 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
253 % DestroyResampleFilter() finalizes and cleans up the resampling
254 % resample_filter as returned by AcquireResampleFilter(), freeing any memory
255 % or other information as needed.
257 % The format of the DestroyResampleFilter method is:
259 % ResampleFilter *DestroyResampleFilter(ResampleFilter *resample_filter)
261 % A description of each parameter follows:
263 % o resample_filter: resampling information structure
266 MagickExport ResampleFilter *DestroyResampleFilter(
267 ResampleFilter *resample_filter)
269 assert(resample_filter != (ResampleFilter *) NULL);
270 assert(resample_filter->signature == MagickSignature);
271 assert(resample_filter->image != (Image *) NULL);
272 if (resample_filter->debug != MagickFalse)
273 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
274 resample_filter->image->filename);
275 resample_filter->view=DestroyCacheView(resample_filter->view);
276 resample_filter->image=DestroyImage(resample_filter->image);
278 resample_filter->filter_def=DestroyResizeFilter(resample_filter->filter_def);
280 resample_filter->signature=(~MagickSignature);
281 resample_filter=(ResampleFilter *) RelinquishMagickMemory(resample_filter);
282 return(resample_filter);
286 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
290 % R e s a m p l e P i x e l C o l o r %
294 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
296 % ResamplePixelColor() samples the pixel values surrounding the location
297 % given using an elliptical weighted average, at the scale previously
298 % calculated, and in the most efficent manner possible for the
299 % VirtualPixelMethod setting.
301 % The format of the ResamplePixelColor method is:
303 % MagickBooleanType ResamplePixelColor(ResampleFilter *resample_filter,
304 % const double u0,const double v0,MagickPixelPacket *pixel)
306 % A description of each parameter follows:
308 % o resample_filter: the resample filter.
310 % o u0,v0: A double representing the center of the area to resample,
311 % The distortion transformed transformed x,y coordinate.
313 % o pixel: the resampled pixel is returned here.
316 MagickExport MagickBooleanType ResamplePixelColor(
317 ResampleFilter *resample_filter,const double u0,const double v0,
318 MagickPixelPacket *pixel)
323 ssize_t u,v, v1, v2, uw, hit;
326 double divisor_c,divisor_m;
327 register double weight;
328 register const PixelPacket *pixels;
329 register const IndexPacket *indexes;
330 assert(resample_filter != (ResampleFilter *) NULL);
331 assert(resample_filter->signature == MagickSignature);
334 /* GetMagickPixelPacket(resample_filter->image,pixel); */
335 if ( resample_filter->do_interpolate ) {
336 status=InterpolateMagickPixelPacket(resample_filter->image,
337 resample_filter->view,resample_filter->interpolate,u0,v0,pixel,
338 resample_filter->exception);
343 (void) FormatLocaleFile(stderr, "u0=%lf; v0=%lf;\n", u0, v0);
347 Does resample area Miss the image?
348 And is that area a simple solid color - then return that color
351 switch ( resample_filter->virtual_pixel ) {
352 case BackgroundVirtualPixelMethod:
353 case ConstantVirtualPixelMethod:
354 case TransparentVirtualPixelMethod:
355 case BlackVirtualPixelMethod:
356 case GrayVirtualPixelMethod:
357 case WhiteVirtualPixelMethod:
358 case MaskVirtualPixelMethod:
359 if ( resample_filter->limit_reached
360 || u0 + resample_filter->Ulimit < 0.0
361 || u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
362 || v0 + resample_filter->Vlimit < 0.0
363 || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows
368 case UndefinedVirtualPixelMethod:
369 case EdgeVirtualPixelMethod:
370 if ( ( u0 + resample_filter->Ulimit < 0.0 && v0 + resample_filter->Vlimit < 0.0 )
371 || ( u0 + resample_filter->Ulimit < 0.0
372 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows )
373 || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
374 && v0 + resample_filter->Vlimit < 0.0 )
375 || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
376 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows )
380 case HorizontalTileVirtualPixelMethod:
381 if ( v0 + resample_filter->Vlimit < 0.0
382 || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows
384 hit++; /* outside the horizontally tiled images. */
386 case VerticalTileVirtualPixelMethod:
387 if ( u0 + resample_filter->Ulimit < 0.0
388 || u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
390 hit++; /* outside the vertically tiled images. */
392 case DitherVirtualPixelMethod:
393 if ( ( u0 + resample_filter->Ulimit < -32.0 && v0 + resample_filter->Vlimit < -32.0 )
394 || ( u0 + resample_filter->Ulimit < -32.0
395 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+32.0 )
396 || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+32.0
397 && v0 + resample_filter->Vlimit < -32.0 )
398 || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+32.0
399 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+32.0 )
403 case TileVirtualPixelMethod:
404 case MirrorVirtualPixelMethod:
405 case RandomVirtualPixelMethod:
406 case HorizontalTileEdgeVirtualPixelMethod:
407 case VerticalTileEdgeVirtualPixelMethod:
408 case CheckerTileVirtualPixelMethod:
409 /* resampling of area is always needed - no VP limits */
413 /* whole area is a solid color -- just return that color */
414 status=InterpolateMagickPixelPacket(resample_filter->image,
415 resample_filter->view,IntegerInterpolatePixel,u0,v0,pixel,
416 resample_filter->exception);
421 Scaling limits reached, return an 'averaged' result.
423 if ( resample_filter->limit_reached ) {
424 switch ( resample_filter->virtual_pixel ) {
425 /* This is always handled by the above, so no need.
426 case BackgroundVirtualPixelMethod:
427 case ConstantVirtualPixelMethod:
428 case TransparentVirtualPixelMethod:
429 case GrayVirtualPixelMethod,
430 case WhiteVirtualPixelMethod
431 case MaskVirtualPixelMethod:
433 case UndefinedVirtualPixelMethod:
434 case EdgeVirtualPixelMethod:
435 case DitherVirtualPixelMethod:
436 case HorizontalTileEdgeVirtualPixelMethod:
437 case VerticalTileEdgeVirtualPixelMethod:
438 /* We need an average edge pixel, from the correct edge!
439 How should I calculate an average edge color?
440 Just returning an averaged neighbourhood,
441 works well in general, but falls down for TileEdge methods.
442 This needs to be done properly!!!!!!
444 status=InterpolateMagickPixelPacket(resample_filter->image,
445 resample_filter->view,AverageInterpolatePixel,u0,v0,pixel,
446 resample_filter->exception);
448 case HorizontalTileVirtualPixelMethod:
449 case VerticalTileVirtualPixelMethod:
450 /* just return the background pixel - Is there more direct way? */
451 status=InterpolateMagickPixelPacket(resample_filter->image,
452 resample_filter->view,IntegerInterpolatePixel,-1.0,-1.0,pixel,
453 resample_filter->exception);
455 case TileVirtualPixelMethod:
456 case MirrorVirtualPixelMethod:
457 case RandomVirtualPixelMethod:
458 case CheckerTileVirtualPixelMethod:
460 /* generate a average color of the WHOLE image */
461 if ( resample_filter->average_defined == MagickFalse ) {
468 GetMagickPixelPacket(resample_filter->image,(MagickPixelPacket *)
469 &resample_filter->average_pixel);
470 resample_filter->average_defined=MagickTrue;
472 /* Try to get an averaged pixel color of whole image */
473 average_image=ResizeImage(resample_filter->image,1,1,BoxFilter,1.0,
474 resample_filter->exception);
475 if (average_image == (Image *) NULL)
477 *pixel=resample_filter->average_pixel; /* FAILED */
480 average_view=AcquireCacheView(average_image);
481 pixels=(PixelPacket *)GetCacheViewVirtualPixels(average_view,0,0,1,1,
482 resample_filter->exception);
483 if (pixels == (const PixelPacket *) NULL) {
484 average_view=DestroyCacheView(average_view);
485 average_image=DestroyImage(average_image);
486 *pixel=resample_filter->average_pixel; /* FAILED */
489 indexes=(IndexPacket *) GetCacheViewAuthenticIndexQueue(average_view);
490 SetMagickPixelPacket(resample_filter->image,pixels,indexes,
491 &(resample_filter->average_pixel));
492 average_view=DestroyCacheView(average_view);
493 average_image=DestroyImage(average_image);
495 if ( resample_filter->virtual_pixel == CheckerTileVirtualPixelMethod )
497 /* CheckerTile is avergae of image average half background */
498 /* FUTURE: replace with a 50% blend of both pixels */
500 weight = QuantumScale*((MagickRealType)(QuantumRange-
501 resample_filter->average_pixel.opacity));
502 resample_filter->average_pixel.red *= weight;
503 resample_filter->average_pixel.green *= weight;
504 resample_filter->average_pixel.blue *= weight;
507 weight = QuantumScale*((MagickRealType)(QuantumRange-
508 resample_filter->image->background_color.opacity));
509 resample_filter->average_pixel.red +=
510 weight*resample_filter->image->background_color.red;
511 resample_filter->average_pixel.green +=
512 weight*resample_filter->image->background_color.green;
513 resample_filter->average_pixel.blue +=
514 weight*resample_filter->image->background_color.blue;
515 resample_filter->average_pixel.opacity +=
516 resample_filter->image->background_color.opacity;
519 resample_filter->average_pixel.red /= divisor_c;
520 resample_filter->average_pixel.green /= divisor_c;
521 resample_filter->average_pixel.blue /= divisor_c;
522 resample_filter->average_pixel.opacity /= 2;
526 *pixel=resample_filter->average_pixel;
533 Initialize weighted average data collection
538 pixel->red = pixel->green = pixel->blue = 0.0;
539 if (pixel->matte != MagickFalse) pixel->opacity = 0.0;
540 if (pixel->colorspace == CMYKColorspace) pixel->index = 0.0;
543 Determine the parellelogram bounding box fitted to the ellipse
544 centered at u0,v0. This area is bounding by the lines...
546 v1 = (ssize_t)ceil(v0 - resample_filter->Vlimit); /* range of scan lines */
547 v2 = (ssize_t)floor(v0 + resample_filter->Vlimit);
549 /* scan line start and width accross the parallelogram */
550 u1 = u0 + (v1-v0)*resample_filter->slope - resample_filter->Uwidth;
551 uw = (ssize_t)(2.0*resample_filter->Uwidth)+1;
554 (void) FormatLocaleFile(stderr, "v1=%ld; v2=%ld\n", (long)v1, (long)v2);
555 (void) FormatLocaleFile(stderr, "u1=%ld; uw=%ld\n", (long)u1, (long)uw);
557 # define DEBUG_HIT_MISS 0 /* only valid if DEBUG_ELLIPSE is enabled */
561 Do weighted resampling of all pixels, within the scaled ellipse,
562 bound by a Parellelogram fitted to the ellipse.
564 DDQ = 2*resample_filter->A;
565 for( v=v1; v<=v2; v++ ) {
567 long uu = ceil(u1); /* actual pixel location (for debug only) */
568 (void) FormatLocaleFile(stderr, "# scan line from pixel %ld, %ld\n", (long)uu, (long)v);
570 u = (ssize_t)ceil(u1); /* first pixel in scanline */
571 u1 += resample_filter->slope; /* start of next scan line */
574 /* location of this first pixel, relative to u0,v0 */
578 /* Q = ellipse quotent ( if Q<F then pixel is inside ellipse) */
579 Q = (resample_filter->A*U + resample_filter->B*V)*U + resample_filter->C*V*V;
580 DQ = resample_filter->A*(2.0*U+1) + resample_filter->B*V;
582 /* get the scanline of pixels for this v */
583 pixels=GetCacheViewVirtualPixels(resample_filter->view,u,v,(size_t) uw,
584 1,resample_filter->exception);
585 if (pixels == (const PixelPacket *) NULL)
587 indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
589 /* count up the weighted pixel colors */
590 for( u=0; u<uw; u++ ) {
592 /* Note that the ellipse has been pre-scaled so F = WLUT_WIDTH */
593 if ( Q < (double)WLUT_WIDTH ) {
594 weight = resample_filter->filter_lut[(int)Q];
596 /* Note that the ellipse has been pre-scaled so F = support^2 */
597 if ( Q < (double)resample_filter->F ) {
598 weight = GetResizeFilterWeight(resample_filter->filter_def,
599 sqrt(Q)); /* a SquareRoot! Arrggghhhhh... */
602 pixel->opacity += weight*pixels->opacity;
605 if (pixel->matte != MagickFalse)
606 weight *= QuantumScale*((MagickRealType)(QuantumRange-pixels->opacity));
607 pixel->red += weight*pixels->red;
608 pixel->green += weight*pixels->green;
609 pixel->blue += weight*pixels->blue;
610 if (pixel->colorspace == CMYKColorspace)
611 pixel->index += weight*(*indexes);
616 /* mark the pixel according to hit/miss of the ellipse */
617 (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
618 (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
619 (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
620 (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
622 (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
623 (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
624 (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
625 (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
638 (void) FormatLocaleFile(stderr, "Hit=%ld; Total=%ld;\n", (long)hit, (long)uw*(v2-v1) );
642 Result sanity check -- this should NOT happen
645 /* not enough pixels in resampling, resort to direct interpolation */
646 #if DEBUG_NO_PIXEL_HIT
647 pixel->opacity = pixel->red = pixel->green = pixel->blue = 0;
648 pixel->red = QuantumRange; /* show pixels for which EWA fails */
650 status=InterpolateMagickPixelPacket(resample_filter->image,
651 resample_filter->view,resample_filter->interpolate,u0,v0,pixel,
652 resample_filter->exception);
658 Finialize results of resampling
660 divisor_m = 1.0/divisor_m;
661 pixel->opacity = (MagickRealType) ClampToQuantum(divisor_m*pixel->opacity);
662 divisor_c = 1.0/divisor_c;
663 pixel->red = (MagickRealType) ClampToQuantum(divisor_c*pixel->red);
664 pixel->green = (MagickRealType) ClampToQuantum(divisor_c*pixel->green);
665 pixel->blue = (MagickRealType) ClampToQuantum(divisor_c*pixel->blue);
666 if (pixel->colorspace == CMYKColorspace)
667 pixel->index = (MagickRealType) ClampToQuantum(divisor_c*pixel->index);
673 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
677 - C l a m p U p A x e s %
681 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
683 % ClampUpAxes() function converts the input vectors into a major and
684 % minor axis unit vectors, and their magnitude. This allows us to
685 % ensure that the ellipse generated is never smaller than the unit
686 % circle and thus never too small for use in EWA resampling.
688 % This purely mathematical 'magic' was provided by Professor Nicolas
689 % Robidoux and his Masters student Chantal Racette.
691 % Reference: "We Recommend Singular Value Decomposition", David Austin
692 % http://www.ams.org/samplings/feature-column/fcarc-svd
694 % By generating major and minor axis vectors, we can actually use the
695 % ellipse in its "canonical form", by remapping the dx,dy of the
696 % sampled point into distances along the major and minor axis unit
699 % Reference: http://en.wikipedia.org/wiki/Ellipse#Canonical_form
701 static inline void ClampUpAxes(const double dux,
707 double *major_unit_x,
708 double *major_unit_y,
709 double *minor_unit_x,
710 double *minor_unit_y)
713 * ClampUpAxes takes an input 2x2 matrix
715 * [ a b ] = [ dux duy ]
716 * [ c d ] = [ dvx dvy ]
718 * and computes from it the major and minor axis vectors [major_x,
719 * major_y] and [minor_x,minor_y] of the smallest ellipse containing
720 * both the unit disk and the ellipse which is the image of the unit
721 * disk by the linear transformation
723 * [ dux duy ] [S] = [s]
724 * [ dvx dvy ] [T] = [t]
726 * (The vector [S,T] is the difference between a position in output
727 * space and [X,Y]; the vector [s,t] is the difference between a
728 * position in input space and [x,y].)
733 * major_mag is the half-length of the major axis of the "new"
736 * minor_mag is the half-length of the minor axis of the "new"
739 * major_unit_x is the x-coordinate of the major axis direction vector
740 * of both the "old" and "new" ellipses.
742 * major_unit_y is the y-coordinate of the major axis direction vector.
744 * minor_unit_x is the x-coordinate of the minor axis direction vector.
746 * minor_unit_y is the y-coordinate of the minor axis direction vector.
748 * Unit vectors are useful for computing projections, in particular,
749 * to compute the distance between a point in output space and the
750 * center of a unit disk in output space, using the position of the
751 * corresponding point [s,t] in input space. Following the clamping,
752 * the square of this distance is
754 * ( ( s * major_unit_x + t * major_unit_y ) / major_mag )^2
756 * ( ( s * minor_unit_x + t * minor_unit_y ) / minor_mag )^2
758 * If such distances will be computed for many [s,t]'s, it makes
759 * sense to actually compute the reciprocal of major_mag and
760 * minor_mag and multiply them by the above unit lengths.
762 * Now, if you want to modify the input pair of tangent vectors so
763 * that it defines the modified ellipse, all you have to do is set
765 * newdux = major_mag * major_unit_x
766 * newdvx = major_mag * major_unit_y
767 * newduy = minor_mag * minor_unit_x = minor_mag * -major_unit_y
768 * newdvy = minor_mag * minor_unit_y = minor_mag * major_unit_x
770 * and use these tangent vectors as if they were the original ones.
771 * Usually, this is a drastic change in the tangent vectors even if
772 * the singular values are not clamped; for example, the minor axis
773 * vector always points in a direction which is 90 degrees
774 * counterclockwise from the direction of the major axis vector.
779 * GOAL: Fix things so that the pullback, in input space, of a disk
780 * of radius r in output space is an ellipse which contains, at
781 * least, a disc of radius r. (Make this hold for any r>0.)
783 * ESSENCE OF THE METHOD: Compute the product of the first two
784 * factors of an SVD of the linear transformation defining the
785 * ellipse and make sure that both its columns have norm at least 1.
786 * Because rotations and reflexions map disks to themselves, it is
787 * not necessary to compute the third (rightmost) factor of the SVD.
789 * DETAILS: Find the singular values and (unit) left singular
790 * vectors of Jinv, clampling up the singular values to 1, and
791 * multiply the unit left singular vectors by the new singular
792 * values in order to get the minor and major ellipse axis vectors.
794 * Image resampling context:
796 * The Jacobian matrix of the transformation at the output point
797 * under consideration is defined as follows:
799 * Consider the transformation (x,y) -> (X,Y) from input locations
800 * to output locations. (Anthony Thyssen, elsewhere in resample.c,
801 * uses the notation (u,v) -> (x,y).)
803 * The Jacobian matrix of the transformation at (x,y) is equal to
805 * J = [ A, B ] = [ dX/dx, dX/dy ]
806 * [ C, D ] [ dY/dx, dY/dy ]
808 * that is, the vector [A,C] is the tangent vector corresponding to
809 * input changes in the horizontal direction, and the vector [B,D]
810 * is the tangent vector corresponding to input changes in the
811 * vertical direction.
813 * In the context of resampling, it is natural to use the inverse
814 * Jacobian matrix Jinv because resampling is generally performed by
815 * pulling pixel locations in the output image back to locations in
816 * the input image. Jinv is
818 * Jinv = [ a, b ] = [ dx/dX, dx/dY ]
819 * [ c, d ] [ dy/dX, dy/dY ]
821 * Note: Jinv can be computed from J with the following matrix
824 * Jinv = 1/(A*D-B*C) [ D, -B ]
827 * What we do is modify Jinv so that it generates an ellipse which
828 * is as close as possible to the original but which contains the
829 * unit disk. This can be accomplished as follows:
835 * be an SVD decomposition of Jinv. (The SVD is not unique, but the
836 * final ellipse does not depend on the particular SVD.)
838 * We could clamp up the entries of the diagonal matrix Sigma so
839 * that they are at least 1, and then set
841 * Jinv = U newSigma V^T.
843 * However, we do not need to compute V for the following reason:
844 * V^T is an orthogonal matrix (that is, it represents a combination
845 * of rotations and reflexions) so that it maps the unit circle to
846 * itself. For this reason, the exact value of V does not affect the
847 * final ellipse, and we can choose V to be the identity
852 * In the end, we return the two diagonal entries of newSigma
853 * together with the two columns of U.
856 * ClampUpAxes was written by Nicolas Robidoux and Chantal Racette
857 * of Laurentian University with insightful suggestions from Anthony
858 * Thyssen and funding from the National Science and Engineering
859 * Research Council of Canada. It is distinguished from its
860 * predecessors by its efficient handling of degenerate cases.
862 * The idea of clamping up the EWA ellipse's major and minor axes so
863 * that the result contains the reconstruction kernel filter support
864 * is taken from Andreas Gustaffson's Masters thesis "Interactive
865 * Image Warping", Helsinki University of Technology, Faculty of
866 * Information Technology, 59 pages, 1993 (see Section 3.6).
868 * The use of the SVD to clamp up the singular values of the
869 * Jacobian matrix of the pullback transformation for EWA resampling
870 * is taken from the astrophysicist Craig DeForest. It is
871 * implemented in his PDL::Transform code (PDL = Perl Data
874 const double a = dux;
875 const double b = duy;
876 const double c = dvx;
877 const double d = dvy;
879 * n is the matrix Jinv * transpose(Jinv). Eigenvalues of n are the
880 * squares of the singular values of Jinv.
882 const double aa = a*a;
883 const double bb = b*b;
884 const double cc = c*c;
885 const double dd = d*d;
887 * Eigenvectors of n are left singular vectors of Jinv.
889 const double n11 = aa+bb;
890 const double n12 = a*c+b*d;
891 const double n21 = n12;
892 const double n22 = cc+dd;
893 const double det = a*d-b*c;
894 const double twice_det = det+det;
895 const double frobenius_squared = n11+n22;
896 const double discriminant =
897 (frobenius_squared+twice_det)*(frobenius_squared-twice_det);
898 const double sqrt_discriminant = sqrt(discriminant);
900 * s1 is the largest singular value of the inverse Jacobian
901 * matrix. In other words, its reciprocal is the smallest singular
902 * value of the Jacobian matrix itself.
903 * If s1 = 0, both singular values are 0, and any orthogonal pair of
904 * left and right factors produces a singular decomposition of Jinv.
907 * Initially, we only compute the squares of the singular values.
909 const double s1s1 = 0.5*(frobenius_squared+sqrt_discriminant);
911 * s2 the smallest singular value of the inverse Jacobian
912 * matrix. Its reciprocal is the largest singular value of the
913 * Jacobian matrix itself.
915 const double s2s2 = 0.5*(frobenius_squared-sqrt_discriminant);
916 const double s1s1minusn11 = s1s1-n11;
917 const double s1s1minusn22 = s1s1-n22;
919 * u1, the first column of the U factor of a singular decomposition
920 * of Jinv, is a (non-normalized) left singular vector corresponding
921 * to s1. It has entries u11 and u21. We compute u1 from the fact
922 * that it is an eigenvector of n corresponding to the eigenvalue
925 const double s1s1minusn11_squared = s1s1minusn11*s1s1minusn11;
926 const double s1s1minusn22_squared = s1s1minusn22*s1s1minusn22;
928 * The following selects the largest row of n-s1^2 I as the one
929 * which is used to find the eigenvector. If both s1^2-n11 and
930 * s1^2-n22 are zero, n-s1^2 I is the zero matrix. In that case,
931 * any vector is an eigenvector; in addition, norm below is equal to
932 * zero, and, in exact arithmetic, this is the only case in which
933 * norm = 0. So, setting u1 to the simple but arbitrary vector [1,0]
934 * if norm = 0 safely takes care of all cases.
936 const double temp_u11 =
937 ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? n12 : s1s1minusn22 );
938 const double temp_u21 =
939 ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? s1s1minusn11 : n21 );
940 const double norm = sqrt(temp_u11*temp_u11+temp_u21*temp_u21);
942 * Finalize the entries of first left singular vector (associated
943 * with the largest singular value).
945 const double u11 = ( (norm>0.0) ? temp_u11/norm : 1.0 );
946 const double u21 = ( (norm>0.0) ? temp_u21/norm : 0.0 );
948 * Clamp the singular values up to 1.
950 *major_mag = ( (s1s1<=1.0) ? 1.0 : sqrt(s1s1) );
951 *minor_mag = ( (s2s2<=1.0) ? 1.0 : sqrt(s2s2) );
953 * Return the unit major and minor axis direction vectors.
957 *minor_unit_x = -u21;
963 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
967 % S c a l e R e s a m p l e F i l t e r %
971 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
973 % ScaleResampleFilter() does all the calculations needed to resample an image
974 % at a specific scale, defined by two scaling vectors. This not using
975 % a orthogonal scaling, but two distorted scaling vectors, to allow the
976 % generation of a angled ellipse.
978 % As only two deritive scaling vectors are used the center of the ellipse
979 % must be the center of the lookup. That is any curvature that the
980 % distortion may produce is discounted.
982 % The input vectors are produced by either finding the derivitives of the
983 % distortion function, or the partial derivitives from a distortion mapping.
984 % They do not need to be the orthogonal dx,dy scaling vectors, but can be
985 % calculated from other derivatives. For example you could use dr,da/r
986 % polar coordinate vector scaling vectors
988 % If u,v = DistortEquation(x,y) OR u = Fu(x,y); v = Fv(x,y)
989 % Then the scaling vectors are determined from the deritives...
990 % du/dx, dv/dx and du/dy, dv/dy
991 % If the resulting scaling vectors is othogonally aligned then...
992 % dv/dx = 0 and du/dy = 0
993 % Producing an othogonally alligned ellipse in source space for the area to
996 % Note that scaling vectors are different to argument order. Argument order
997 % is the general order the deritives are extracted from the distortion
998 % equations, and not the scaling vectors. As such the middle two vaules
999 % may be swapped from what you expect. Caution is advised.
1001 % WARNING: It is assumed that any SetResampleFilter() method call will
1002 % always be performed before the ScaleResampleFilter() method, so that the
1003 % size of the ellipse will match the support for the resampling filter being
1006 % The format of the ScaleResampleFilter method is:
1008 % void ScaleResampleFilter(const ResampleFilter *resample_filter,
1009 % const double dux,const double duy,const double dvx,const double dvy)
1011 % A description of each parameter follows:
1013 % o resample_filter: the resampling resample_filterrmation defining the
1014 % image being resampled
1016 % o dux,duy,dvx,dvy:
1017 % The deritives or scaling vectors defining the EWA ellipse.
1018 % NOTE: watch the order, which is based on the order deritives
1019 % are usally determined from distortion equations (see above).
1020 % The middle two values may need to be swapped if you are thinking
1021 % in terms of scaling vectors.
1024 MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
1025 const double dux,const double duy,const double dvx,const double dvy)
1029 assert(resample_filter != (ResampleFilter *) NULL);
1030 assert(resample_filter->signature == MagickSignature);
1032 resample_filter->limit_reached = MagickFalse;
1034 /* A 'point' filter forces use of interpolation instead of area sampling */
1035 if ( resample_filter->filter == PointFilter )
1036 return; /* EWA turned off - nothing to do */
1039 (void) FormatLocaleFile(stderr, "# -----\n" );
1040 (void) FormatLocaleFile(stderr, "dux=%lf; dvx=%lf; duy=%lf; dvy=%lf;\n",
1041 dux, dvx, duy, dvy);
1044 /* Find Ellipse Coefficents such that
1045 A*u^2 + B*u*v + C*v^2 = F
1046 With u,v relative to point around which we are resampling.
1047 And the given scaling dx,dy vectors in u,v space
1048 du/dx,dv/dx and du/dy,dv/dy
1051 /* Direct conversion of derivatives into elliptical coefficients
1052 However when magnifying images, the scaling vectors will be small
1053 resulting in a ellipse that is too small to sample properly.
1054 As such we need to clamp the major/minor axis to a minumum of 1.0
1055 to prevent it getting too small.
1065 ClampUpAxes(dux,dvx,duy,dvy, &major_mag, &minor_mag,
1066 &major_x, &major_y, &minor_x, &minor_y);
1067 major_x *= major_mag; major_y *= major_mag;
1068 minor_x *= minor_mag; minor_y *= minor_mag;
1070 (void) FormatLocaleFile(stderr, "major_x=%lf; major_y=%lf; minor_x=%lf; minor_y=%lf;\n",
1071 major_x, major_y, minor_x, minor_y);
1073 A = major_y*major_y+minor_y*minor_y;
1074 B = -2.0*(major_x*major_y+minor_x*minor_y);
1075 C = major_x*major_x+minor_x*minor_x;
1076 F = major_mag*minor_mag;
1077 F *= F; /* square it */
1079 #else /* raw unclamped EWA */
1080 A = dvx*dvx+dvy*dvy;
1081 B = -2.0*(dux*dvx+duy*dvy);
1082 C = dux*dux+duy*duy;
1083 F = dux*dvy-duy*dvx;
1084 F *= F; /* square it */
1085 #endif /* EWA_CLAMP */
1089 This Paul Heckbert's "Higher Quality EWA" formula, from page 60 in his
1090 thesis, which adds a unit circle to the elliptical area so as to do both
1091 Reconstruction and Prefiltering of the pixels in the resampling. It also
1092 means it is always likely to have at least 4 pixels within the area of the
1093 ellipse, for weighted averaging. No scaling will result with F == 4.0 and
1094 a circle of radius 2.0, and F smaller than this means magnification is
1097 NOTE: This method produces a very blury result at near unity scale while
1098 producing perfect results for strong minitification and magnifications.
1100 However filter support is fixed to 2.0 (no good for Windowed Sinc filters)
1102 A = dvx*dvx+dvy*dvy+1;
1103 B = -2.0*(dux*dvx+duy*dvy);
1104 C = dux*dux+duy*duy+1;
1109 (void) FormatLocaleFile(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
1111 /* Figure out the various information directly about the ellipse.
1112 This information currently not needed at this time, but may be
1113 needed later for better limit determination.
1115 It is also good to have as a record for future debugging
1117 { double alpha, beta, gamma, Major, Minor;
1118 double Eccentricity, Ellipse_Area, Ellipse_Angle;
1122 gamma = sqrt(beta*beta + B*B );
1124 if ( alpha - gamma <= MagickEpsilon )
1127 Major = sqrt(2*F/(alpha - gamma));
1128 Minor = sqrt(2*F/(alpha + gamma));
1130 (void) FormatLocaleFile(stderr, "# Major=%lf; Minor=%lf\n", Major, Minor );
1132 /* other information about ellipse include... */
1133 Eccentricity = Major/Minor;
1134 Ellipse_Area = MagickPI*Major*Minor;
1135 Ellipse_Angle = atan2(B, A-C);
1137 (void) FormatLocaleFile(stderr, "# Angle=%lf Area=%lf\n",
1138 RadiansToDegrees(Ellipse_Angle), Ellipse_Area);
1142 /* If one or both of the scaling vectors is impossibly large
1143 (producing a very large raw F value), we may as well not bother
1144 doing any form of resampling since resampled area is very large.
1145 In this case some alternative means of pixel sampling, such as
1146 the average of the whole image is needed to get a reasonable
1147 result. Calculate only as needed.
1149 if ( (4*A*C - B*B) > MagickHuge ) {
1150 resample_filter->limit_reached = MagickTrue;
1154 /* Scale ellipse to match the filters support
1155 (that is, multiply F by the square of the support).
1157 F *= resample_filter->support;
1158 F *= resample_filter->support;
1160 /* Orthogonal bounds of the ellipse */
1161 resample_filter->Ulimit = sqrt(C*F/(A*C-0.25*B*B));
1162 resample_filter->Vlimit = sqrt(A*F/(A*C-0.25*B*B));
1164 /* Horizontally aligned parallelogram fitted to Ellipse */
1165 resample_filter->Uwidth = sqrt(F/A); /* Half of the parallelogram width */
1166 resample_filter->slope = -B/(2.0*A); /* Reciprocal slope of the parallelogram */
1169 (void) FormatLocaleFile(stderr, "Ulimit=%lf; Vlimit=%lf; UWidth=%lf; Slope=%lf;\n",
1170 resample_filter->Ulimit, resample_filter->Vlimit,
1171 resample_filter->Uwidth, resample_filter->slope );
1174 /* Check the absolute area of the parallelogram involved.
1175 * This limit needs more work, as it is too slow for larger images
1176 * with tiled views of the horizon.
1178 if ( (resample_filter->Uwidth * resample_filter->Vlimit)
1179 > (4.0*resample_filter->image_area)) {
1180 resample_filter->limit_reached = MagickTrue;
1184 /* Scale ellipse formula to directly index the Filter Lookup Table */
1185 { register double scale;
1187 /* scale so that F = WLUT_WIDTH; -- hardcoded */
1188 scale = (double)WLUT_WIDTH/F;
1190 /* scale so that F = resample_filter->F (support^2) */
1191 scale = resample_filter->F/F;
1193 resample_filter->A = A*scale;
1194 resample_filter->B = B*scale;
1195 resample_filter->C = C*scale;
1200 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1204 % S e t R e s a m p l e F i l t e r %
1208 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1210 % SetResampleFilter() set the resampling filter lookup table based on a
1211 % specific filter. Note that the filter is used as a radial filter not as a
1212 % two pass othogonally aligned resampling filter.
1214 % The default Filter, is Gaussian, which is the standard filter used by the
1215 % original paper on the Elliptical Weighted Everage Algorithm. However other
1216 % filters can also be used.
1218 % The format of the SetResampleFilter method is:
1220 % void SetResampleFilter(ResampleFilter *resample_filter,
1221 % const FilterTypes filter,const double blur)
1223 % A description of each parameter follows:
1225 % o resample_filter: resampling resample_filterrmation structure
1227 % o filter: the resize filter for elliptical weighting LUT
1229 % o blur: filter blur factor (radial scaling) for elliptical weighting LUT
1232 MagickExport void SetResampleFilter(ResampleFilter *resample_filter,
1233 const FilterTypes filter,const double blur)
1238 assert(resample_filter != (ResampleFilter *) NULL);
1239 assert(resample_filter->signature == MagickSignature);
1241 resample_filter->do_interpolate = MagickFalse;
1242 resample_filter->filter = filter;
1244 if ( filter == PointFilter )
1246 resample_filter->do_interpolate = MagickTrue;
1247 return; /* EWA turned off - nothing more to do */
1250 /* Set a default cylindrical filter of a 'low blur' Jinc windowed Jinc */
1251 if ( filter == UndefinedFilter )
1252 resample_filter->filter = RobidouxFilter;
1254 resize_filter = AcquireResizeFilter(resample_filter->image,
1255 resample_filter->filter,blur,MagickTrue,resample_filter->exception);
1256 if (resize_filter == (ResizeFilter *) NULL)
1258 (void) ThrowMagickException(resample_filter->exception,GetMagickModule(),
1259 ModuleError, "UnableToSetFilteringValue",
1260 "Fall back to default EWA gaussian filter");
1261 resample_filter->filter = PointFilter;
1264 /* Get the practical working support for the filter,
1265 * after any API call blur factors have been accoded for.
1268 resample_filter->support = GetResizeFilterSupport(resize_filter);
1270 resample_filter->support = 2.0; /* fixed support size for HQ-EWA */
1274 /* Fill the LUT with the weights from the selected filter function */
1279 /* Scale radius so the filter LUT covers the full support range */
1280 r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
1281 for(Q=0; Q<WLUT_WIDTH; Q++)
1282 resample_filter->filter_lut[Q] = (double)
1283 GetResizeFilterWeight(resize_filter,sqrt((double)Q)*r_scale);
1285 /* finished with the resize filter */
1286 resize_filter = DestroyResizeFilter(resize_filter);
1289 /* save the filter and the scaled ellipse bounds needed for filter */
1290 resample_filter->filter_def = resize_filter;
1291 resample_filter->F = resample_filter->support*resample_filter->support;
1295 Adjust the scaling of the default unit circle
1296 This assumes that any real scaling changes will always
1297 take place AFTER the filter method has been initialized.
1299 ScaleResampleFilter(resample_filter, 1.0, 0.0, 0.0, 1.0);
1302 /* This is old code kept as a reference only. It is very wrong,
1303 and I don't understand exactly what it was attempting to do.
1306 Create Normal Gaussian 2D Filter Weighted Lookup Table.
1307 A normal EWA guassual lookup would use exp(Q*ALPHA)
1308 where Q = distance squared from 0.0 (center) to 1.0 (edge)
1309 and ALPHA = -4.0*ln(2.0) ==> -2.77258872223978123767
1310 The table is of length 1024, and equates to support radius of 2.0
1311 thus needs to be scaled by ALPHA*4/1024 and any blur factor squared
1313 The above came from some reference code provided by Fred Weinhaus
1314 and seems to have been a guess that was appropriate for its use
1315 in a 3d perspective landscape mapping program.
1317 r_scale = -2.77258872223978123767/(WLUT_WIDTH*blur*blur);
1318 for(Q=0; Q<WLUT_WIDTH; Q++)
1319 resample_filter->filter_lut[Q] = exp((double)Q*r_scale);
1320 resample_filter->support = WLUT_WIDTH;
1325 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1333 /* Scale radius so the filter LUT covers the full support range */
1334 r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
1335 if (IsMagickTrue(GetImageArtifact(resample_filter->image,"resample:verbose")) )
1337 /* Debug output of the filter weighting LUT
1338 Gnuplot the LUT with hoizontal adjusted to 'r' using...
1339 plot [0:2][-.2:1] "lut.dat" using (sqrt($0/1024)*2):1 with lines
1340 The filter values is normalized for comparision
1343 printf("# Resampling Filter LUT (%d values)\n", WLUT_WIDTH);
1345 printf("# Note: values in table are using a squared radius lookup.\n");
1346 printf("# And the whole table represents the filters support.\n");
1347 printf("\n"); /* generates a 'break' in gnuplot if multiple outputs */
1348 for(Q=0; Q<WLUT_WIDTH; Q++)
1349 printf("%8.*g %.*g\n",
1350 GetMagickPrecision(),sqrt((double)Q)*r_scale,
1351 GetMagickPrecision(),resample_filter->filter_lut[Q] );
1353 /* output the above once only for each image, and each setting */
1354 (void) DeleteImageArtifact(resample_filter->image,"resample:verbose");
1356 #endif /* FILTER_LUT */
1361 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1365 % 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 %
1369 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1371 % SetResampleFilterInterpolateMethod() sets the resample filter interpolation
1374 % The format of the SetResampleFilterInterpolateMethod method is:
1376 % MagickBooleanType SetResampleFilterInterpolateMethod(
1377 % ResampleFilter *resample_filter,const InterpolateMethod method)
1379 % A description of each parameter follows:
1381 % o resample_filter: the resample filter.
1383 % o method: the interpolation method.
1386 MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
1387 ResampleFilter *resample_filter,const InterpolatePixelMethod method)
1389 assert(resample_filter != (ResampleFilter *) NULL);
1390 assert(resample_filter->signature == MagickSignature);
1391 assert(resample_filter->image != (Image *) NULL);
1392 if (resample_filter->debug != MagickFalse)
1393 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1394 resample_filter->image->filename);
1395 resample_filter->interpolate=method;
1400 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1404 % 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 %
1408 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1410 % SetResampleFilterVirtualPixelMethod() changes the virtual pixel method
1411 % associated with the specified resample filter.
1413 % The format of the SetResampleFilterVirtualPixelMethod method is:
1415 % MagickBooleanType SetResampleFilterVirtualPixelMethod(
1416 % ResampleFilter *resample_filter,const VirtualPixelMethod method)
1418 % A description of each parameter follows:
1420 % o resample_filter: the resample filter.
1422 % o method: the virtual pixel method.
1425 MagickExport MagickBooleanType SetResampleFilterVirtualPixelMethod(
1426 ResampleFilter *resample_filter,const VirtualPixelMethod method)
1428 assert(resample_filter != (ResampleFilter *) NULL);
1429 assert(resample_filter->signature == MagickSignature);
1430 assert(resample_filter->image != (Image *) NULL);
1431 if (resample_filter->debug != MagickFalse)
1432 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1433 resample_filter->image->filename);
1434 resample_filter->virtual_pixel=method;
1435 if (method != UndefinedVirtualPixelMethod)
1436 (void) SetCacheViewVirtualPixelMethod(resample_filter->view,method);