MagickCore 6.9.13
Loading...
Searching...
No Matches
resample.c
1/*
2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3% %
4% %
5% %
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 %
11% %
12% %
13% MagickCore Pixel Resampling Methods %
14% %
15% Software Design %
16% Cristy %
17% Anthony Thyssen %
18% August 2007 %
19% %
20% %
21% Copyright 1999 ImageMagick Studio LLC, a non-profit organization %
22% dedicated to making software imaging solutions freely available. %
23% %
24% You may not use this file except in compliance with the License. You may %
25% obtain a copy of the License at %
26% %
27% https://imagemagick.org/script/license.php %
28% %
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. %
34% %
35%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36%
37%
38*/
39␌
40/*
41 Include declarations.
42*/
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/resource_.h"
63#include "magick/transform.h"
64#include "magick/signature-private.h"
65#include "magick/token.h"
66#include "magick/utility.h"
67#include "magick/option.h"
68/*
69 EWA Resampling Options
70*/
71
72/* select ONE resampling method */
73#define EWA 1 /* Normal EWA handling - raw or clamped */
74 /* if 0 then use "High Quality EWA" */
75#define EWA_CLAMP 1 /* EWA Clamping from Nicolas Robidoux */
76
77#define FILTER_LUT 1 /* Use a LUT rather then direct filter calls */
78
79/* output debugging information */
80#define DEBUG_ELLIPSE 0 /* output ellipse info for debug */
81#define DEBUG_HIT_MISS 0 /* output hit/miss pixels (as gnuplot commands) */
82#define DEBUG_NO_PIXEL_HIT 0 /* Make pixels that fail to hit anything - RED */
83
84#if ! FILTER_DIRECT
85#define WLUT_WIDTH 1024 /* size of the filter cache */
86#endif
87
88/*
89 Typedef declarations.
90*/
91struct _ResampleFilter
92{
93 CacheView
94 *view;
95
96 Image
97 *image;
98
99 ExceptionInfo
100 *exception;
101
102 MagickBooleanType
103 debug;
104
105 /* Information about image being resampled */
106 ssize_t
107 image_area;
108
109 InterpolatePixelMethod
110 interpolate;
111
112 VirtualPixelMethod
113 virtual_pixel;
114
115 FilterTypes
116 filter;
117
118 /* processing settings needed */
119 MagickBooleanType
120 limit_reached,
121 do_interpolate,
122 average_defined;
123
124 MagickPixelPacket
125 average_pixel;
126
127 /* current elliptical area being resampled around center point */
128 double
129 A, B, C,
130 Vlimit, Ulimit, Uwidth, slope;
131
132#if FILTER_LUT
133 /* LUT of weights for filtered average in elliptical area */
134 double
135 filter_lut[WLUT_WIDTH];
136#else
137 /* Use a Direct call to the filter functions */
138 ResizeFilter
139 *filter_def;
140
141 double
142 F;
143#endif
144
145 /* the practical working support of the filter */
146 double
147 support;
148
149 size_t
150 signature;
151};
152␌
153/*
154%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
155% %
156% %
157% %
158% A c q u i r e R e s a m p l e I n f o %
159% %
160% %
161% %
162%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
163%
164% AcquireResampleFilter() initializes the information resample needs do to a
165% scaled lookup of a color from an image, using area sampling.
166%
167% The algorithm is based on a Elliptical Weighted Average, where the pixels
168% found in a large elliptical area is averaged together according to a
169% weighting (filter) function. For more details see "Fundamentals of Texture
170% Mapping and Image Warping" a master's thesis by Paul.S.Heckbert, June 17,
171% 1989. Available for free from, http://www.cs.cmu.edu/~ph/
172%
173% As EWA resampling (or any sort of resampling) can require a lot of
174% calculations to produce a distorted scaling of the source image for each
175% output pixel, the ResampleFilter structure generated holds that information
176% between individual image resampling.
177%
178% This function will make the appropriate AcquireVirtualCacheView() calls
179% to view the image, calling functions do not need to open a cache view.
180%
181% Usage Example...
182% resample_filter=AcquireResampleFilter(image,exception);
183% SetResampleFilter(resample_filter, GaussianFilter, 1.0);
184% for (y=0; y < (ssize_t) image->rows; y++) {
185% for (x=0; x < (ssize_t) image->columns; x++) {
186% u= ....; v= ....;
187% ScaleResampleFilter(resample_filter, ... scaling vectors ...);
188% (void) ResamplePixelColor(resample_filter,u,v,&pixel);
189% ... assign resampled pixel value ...
190% }
191% }
192% DestroyResampleFilter(resample_filter);
193%
194% The format of the AcquireResampleFilter method is:
195%
196% ResampleFilter *AcquireResampleFilter(const Image *image,
197% ExceptionInfo *exception)
198%
199% A description of each parameter follows:
200%
201% o image: the image.
202%
203% o exception: return any errors or warnings in this structure.
204%
205*/
206MagickExport ResampleFilter *AcquireResampleFilter(const Image *image,
207 ExceptionInfo *exception)
208{
209 ResampleFilter
210 *resample_filter;
211
212 assert(image != (Image *) NULL);
213 assert(image->signature == MagickCoreSignature);
214 assert(exception != (ExceptionInfo *) NULL);
215 assert(exception->signature == MagickCoreSignature);
216 if (IsEventLogging() != MagickFalse)
217 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
218
219 resample_filter=(ResampleFilter *) AcquireQuantumMemory(1,
220 sizeof(*resample_filter));
221 if (resample_filter == (ResampleFilter *) NULL)
222 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
223 (void) memset(resample_filter,0,sizeof(*resample_filter));
224
225 resample_filter->exception=exception;
226 resample_filter->image=ReferenceImage((Image *) image);
227 resample_filter->view=AcquireVirtualCacheView(resample_filter->image,exception);
228
229 resample_filter->debug=IsEventLogging();
230 resample_filter->signature=MagickCoreSignature;
231
232 resample_filter->image_area=(ssize_t) (image->columns*image->rows);
233 resample_filter->average_defined = MagickFalse;
234
235 /* initialise the resampling filter settings */
236 SetResampleFilter(resample_filter, image->filter, image->blur);
237 (void) SetResampleFilterInterpolateMethod(resample_filter,
238 image->interpolate);
239 (void) SetResampleFilterVirtualPixelMethod(resample_filter,
240 GetImageVirtualPixelMethod(image));
241
242 return(resample_filter);
243}
244␌
245/*
246%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
247% %
248% %
249% %
250% D e s t r o y R e s a m p l e I n f o %
251% %
252% %
253% %
254%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
255%
256% DestroyResampleFilter() finalizes and cleans up the resampling
257% resample_filter as returned by AcquireResampleFilter(), freeing any memory
258% or other information as needed.
259%
260% The format of the DestroyResampleFilter method is:
261%
262% ResampleFilter *DestroyResampleFilter(ResampleFilter *resample_filter)
263%
264% A description of each parameter follows:
265%
266% o resample_filter: resampling information structure
267%
268*/
269MagickExport ResampleFilter *DestroyResampleFilter(
270 ResampleFilter *resample_filter)
271{
272 assert(resample_filter != (ResampleFilter *) NULL);
273 assert(resample_filter->signature == MagickCoreSignature);
274 assert(resample_filter->image != (Image *) NULL);
275 if (IsEventLogging() != MagickFalse)
276 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
277 resample_filter->image->filename);
278 resample_filter->view=DestroyCacheView(resample_filter->view);
279 resample_filter->image=DestroyImage(resample_filter->image);
280#if ! FILTER_LUT
281 resample_filter->filter_def=DestroyResizeFilter(resample_filter->filter_def);
282#endif
283 resample_filter->signature=(~MagickCoreSignature);
284 resample_filter=(ResampleFilter *) RelinquishMagickMemory(resample_filter);
285 return(resample_filter);
286}
287␌
288/*
289%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
290% %
291% %
292% %
293% R e s a m p l e P i x e l C o l o r %
294% %
295% %
296% %
297%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
298%
299% ResamplePixelColor() samples the pixel values surrounding the location
300% given using an elliptical weighted average, at the scale previously
301% calculated, and in the most efficient manner possible for the
302% VirtualPixelMethod setting.
303%
304% The format of the ResamplePixelColor method is:
305%
306% MagickBooleanType ResamplePixelColor(ResampleFilter *resample_filter,
307% const double u0,const double v0,MagickPixelPacket *pixel)
308%
309% A description of each parameter follows:
310%
311% o resample_filter: the resample filter.
312%
313% o u0,v0: A double representing the center of the area to resample,
314% The distortion transformed x,y coordinate.
315%
316% o pixel: the resampled pixel is returned here.
317%
318*/
319MagickExport MagickBooleanType ResamplePixelColor(
320 ResampleFilter *resample_filter,const double u0,const double v0,
321 MagickPixelPacket *pixel)
322{
323 MagickBooleanType
324 status;
325
326 ssize_t u,v, v1, v2, uw, hit;
327 double u1;
328 double U,V,Q,DQ,DDQ;
329 double divisor_c,divisor_m;
330 double weight;
331 const PixelPacket *pixels;
332 const IndexPacket *indexes;
333 assert(resample_filter != (ResampleFilter *) NULL);
334 assert(resample_filter->signature == MagickCoreSignature);
335
336 status=MagickTrue;
337 /* GetMagickPixelPacket(resample_filter->image,pixel); */
338 if ( resample_filter->do_interpolate ) {
339 status=InterpolateMagickPixelPacket(resample_filter->image,
340 resample_filter->view,resample_filter->interpolate,u0,v0,pixel,
341 resample_filter->exception);
342 return(status);
343 }
344
345#if DEBUG_ELLIPSE
346 (void) FormatLocaleFile(stderr, "u0=%lf; v0=%lf;\n", u0, v0);
347#endif
348
349 /*
350 Does resample area Miss the image Proper?
351 If and that area a simple solid color - then simply return that color!
352 This saves a lot of calculation when resampling outside the bounds of
353 the source image.
354
355 However it probably should be expanded to image bounds plus the filters
356 scaled support size.
357 */
358 hit = 0;
359 switch ( resample_filter->virtual_pixel ) {
360 case BackgroundVirtualPixelMethod:
361 case ConstantVirtualPixelMethod:
362 case TransparentVirtualPixelMethod:
363 case BlackVirtualPixelMethod:
364 case GrayVirtualPixelMethod:
365 case WhiteVirtualPixelMethod:
366 case MaskVirtualPixelMethod:
367 if ( resample_filter->limit_reached
368 || u0 + resample_filter->Ulimit < 0.0
369 || u0 - resample_filter->Ulimit > (double) resample_filter->image->columns-1.0
370 || v0 + resample_filter->Vlimit < 0.0
371 || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows-1.0
372 )
373 hit++;
374 break;
375
376 case UndefinedVirtualPixelMethod:
377 case EdgeVirtualPixelMethod:
378 if ( ( u0 + resample_filter->Ulimit < 0.0 && v0 + resample_filter->Vlimit < 0.0 )
379 || ( u0 + resample_filter->Ulimit < 0.0
380 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows-1.0 )
381 || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns-1.0
382 && v0 + resample_filter->Vlimit < 0.0 )
383 || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns-1.0
384 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows-1.0 )
385 )
386 hit++;
387 break;
388 case HorizontalTileVirtualPixelMethod:
389 if ( v0 + resample_filter->Vlimit < 0.0
390 || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows-1.0
391 )
392 hit++; /* outside the horizontally tiled images. */
393 break;
394 case VerticalTileVirtualPixelMethod:
395 if ( u0 + resample_filter->Ulimit < 0.0
396 || u0 - resample_filter->Ulimit > (double) resample_filter->image->columns-1.0
397 )
398 hit++; /* outside the vertically tiled images. */
399 break;
400 case DitherVirtualPixelMethod:
401 if ( ( u0 + resample_filter->Ulimit < -32.0 && v0 + resample_filter->Vlimit < -32.0 )
402 || ( u0 + resample_filter->Ulimit < -32.0
403 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+31.0 )
404 || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+31.0
405 && v0 + resample_filter->Vlimit < -32.0 )
406 || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+31.0
407 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+31.0 )
408 )
409 hit++;
410 break;
411 case TileVirtualPixelMethod:
412 case MirrorVirtualPixelMethod:
413 case RandomVirtualPixelMethod:
414 case HorizontalTileEdgeVirtualPixelMethod:
415 case VerticalTileEdgeVirtualPixelMethod:
416 case CheckerTileVirtualPixelMethod:
417 /* resampling of area is always needed - no VP limits */
418 break;
419 }
420 if ( hit ) {
421 /* The area being resampled is simply a solid color
422 * just return a single lookup color.
423 *
424 * Should this return the users requested interpolated color?
425 */
426 status=InterpolateMagickPixelPacket(resample_filter->image,
427 resample_filter->view,IntegerInterpolatePixel,u0,v0,pixel,
428 resample_filter->exception);
429 return(status);
430 }
431
432 /*
433 When Scaling limits reached, return an 'averaged' result.
434 */
435 if ( resample_filter->limit_reached ) {
436 switch ( resample_filter->virtual_pixel ) {
437 /* This is always handled by the above, so no need.
438 case BackgroundVirtualPixelMethod:
439 case ConstantVirtualPixelMethod:
440 case TransparentVirtualPixelMethod:
441 case GrayVirtualPixelMethod,
442 case WhiteVirtualPixelMethod
443 case MaskVirtualPixelMethod:
444 */
445 case UndefinedVirtualPixelMethod:
446 case EdgeVirtualPixelMethod:
447 case DitherVirtualPixelMethod:
448 case HorizontalTileEdgeVirtualPixelMethod:
449 case VerticalTileEdgeVirtualPixelMethod:
450 /* We need an average edge pixel, from the correct edge!
451 How should I calculate an average edge color?
452 Just returning an averaged neighbourhood,
453 works well in general, but falls down for TileEdge methods.
454 This needs to be done properly!!!!!!
455 */
456 status=InterpolateMagickPixelPacket(resample_filter->image,
457 resample_filter->view,AverageInterpolatePixel,u0,v0,pixel,
458 resample_filter->exception);
459 break;
460 case HorizontalTileVirtualPixelMethod:
461 case VerticalTileVirtualPixelMethod:
462 /* just return the background pixel - Is there a better way? */
463 status=InterpolateMagickPixelPacket(resample_filter->image,
464 resample_filter->view,IntegerInterpolatePixel,-1.0,-1.0,pixel,
465 resample_filter->exception);
466 break;
467 case TileVirtualPixelMethod:
468 case MirrorVirtualPixelMethod:
469 case RandomVirtualPixelMethod:
470 case CheckerTileVirtualPixelMethod:
471 default:
472 /* generate a average color of the WHOLE image */
473 if ( resample_filter->average_defined == MagickFalse ) {
474 Image
475 *average_image;
476
477 CacheView
478 *average_view;
479
480 GetMagickPixelPacket(resample_filter->image,(MagickPixelPacket *)
481 &resample_filter->average_pixel);
482 resample_filter->average_defined=MagickTrue;
483
484 /* Try to get an averaged pixel color of whole image */
485 average_image=ResizeImage(resample_filter->image,1,1,BoxFilter,1.0,
486 resample_filter->exception);
487 if (average_image == (Image *) NULL)
488 {
489 *pixel=resample_filter->average_pixel; /* FAILED */
490 break;
491 }
492 average_view=AcquireVirtualCacheView(average_image,
493 &average_image->exception);
494 pixels=(PixelPacket *)GetCacheViewVirtualPixels(average_view,0,0,1,1,
495 resample_filter->exception);
496 if (pixels == (const PixelPacket *) NULL) {
497 average_view=DestroyCacheView(average_view);
498 average_image=DestroyImage(average_image);
499 *pixel=resample_filter->average_pixel; /* FAILED */
500 break;
501 }
502 indexes=(IndexPacket *) GetCacheViewAuthenticIndexQueue(average_view);
503 SetMagickPixelPacket(resample_filter->image,pixels,indexes,
504 &(resample_filter->average_pixel));
505 average_view=DestroyCacheView(average_view);
506 average_image=DestroyImage(average_image);
507
508 if ( resample_filter->virtual_pixel == CheckerTileVirtualPixelMethod )
509 {
510 /* CheckerTile is a alpha blend of the image's average pixel
511 color and the current background color */
512
513 /* image's average pixel color */
514 weight = QuantumScale*((MagickRealType) QuantumRange-
515 (MagickRealType) resample_filter->average_pixel.opacity);
516 resample_filter->average_pixel.red *= weight;
517 resample_filter->average_pixel.green *= weight;
518 resample_filter->average_pixel.blue *= weight;
519 divisor_c = weight;
520
521 /* background color */
522 weight = QuantumScale*((MagickRealType) QuantumRange-
523 (MagickRealType)
524 resample_filter->image->background_color.opacity);
525 resample_filter->average_pixel.red += weight*(MagickRealType)
526 resample_filter->image->background_color.red;
527 resample_filter->average_pixel.green += weight*(MagickRealType)
528 resample_filter->image->background_color.green;
529 resample_filter->average_pixel.blue += weight*(MagickRealType)
530 resample_filter->image->background_color.blue;
531 resample_filter->average_pixel.opacity += (MagickRealType)
532 resample_filter->image->background_color.opacity;
533 divisor_c += weight;
534
535 /* alpha blend */
536 resample_filter->average_pixel.red /= divisor_c;
537 resample_filter->average_pixel.green /= divisor_c;
538 resample_filter->average_pixel.blue /= divisor_c;
539 resample_filter->average_pixel.opacity /= 2; /* 50% blend */
540
541 }
542 }
543 *pixel=resample_filter->average_pixel;
544 break;
545 }
546 return(status);
547 }
548
549 /*
550 Initialize weighted average data collection
551 */
552 hit = 0;
553 divisor_c = 0.0;
554 divisor_m = 0.0;
555 pixel->red = pixel->green = pixel->blue = 0.0;
556 if (pixel->matte != MagickFalse) pixel->opacity = 0.0;
557 if (pixel->colorspace == CMYKColorspace) pixel->index = 0.0;
558
559 /*
560 Determine the parallelogram bounding box fitted to the ellipse
561 centered at u0,v0. This area is bounding by the lines...
562 */
563 v1 = (ssize_t)ceil(v0 - resample_filter->Vlimit); /* range of scan lines */
564 v2 = (ssize_t)floor(v0 + resample_filter->Vlimit);
565
566 /* scan line start and width across the parallelogram */
567 u1 = u0 + (v1-v0)*resample_filter->slope - resample_filter->Uwidth;
568 uw = (ssize_t)(2.0*resample_filter->Uwidth)+1;
569
570#if DEBUG_ELLIPSE
571 (void) FormatLocaleFile(stderr, "v1=%ld; v2=%ld\n", (long)v1, (long)v2);
572 (void) FormatLocaleFile(stderr, "u1=%ld; uw=%ld\n", (long)u1, (long)uw);
573#else
574# define DEBUG_HIT_MISS 0 /* only valid if DEBUG_ELLIPSE is enabled */
575#endif
576
577 /*
578 Do weighted resampling of all pixels, within the scaled ellipse,
579 bound by a Parallelogram fitted to the ellipse.
580 */
581 DDQ = 2*resample_filter->A;
582 for( v=v1; v<=v2; v++ ) {
583#if DEBUG_HIT_MISS
584 long uu = ceil(u1); /* actual pixel location (for debug only) */
585 (void) FormatLocaleFile(stderr, "# scan line from pixel %ld, %ld\n", (long)uu, (long)v);
586#endif
587 u = (ssize_t)ceil(u1); /* first pixel in scanline */
588 u1 += resample_filter->slope; /* start of next scan line */
589
590
591 /* location of this first pixel, relative to u0,v0 */
592 U = (double)u-u0;
593 V = (double)v-v0;
594
595 /* Q = ellipse quotient ( if Q<F then pixel is inside ellipse) */
596 Q = (resample_filter->A*U + resample_filter->B*V)*U + resample_filter->C*V*V;
597 DQ = resample_filter->A*(2.0*U+1) + resample_filter->B*V;
598
599 /* get the scanline of pixels for this v */
600 pixels=GetCacheViewVirtualPixels(resample_filter->view,u,v,(size_t) uw,
601 1,resample_filter->exception);
602 if (pixels == (const PixelPacket *) NULL)
603 return(MagickFalse);
604 indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
605
606 /* count up the weighted pixel colors */
607 for( u=0; u<uw; u++ ) {
608 weight = 0.0;
609#if FILTER_LUT
610 /* Note that the ellipse has been pre-scaled so F = WLUT_WIDTH */
611 if (((int) Q >= 0) && ((int) Q < WLUT_WIDTH)) {
612 weight = resample_filter->filter_lut[(int) Q];
613#else
614 /* Note that the ellipse has been pre-scaled so F = support^2 */
615 if ((Q >= 0.0) && (Q < resample_filter->F)) {
616 weight = GetResizeFilterWeight(resample_filter->filter_def,
617 sqrt(Q)); /* a SquareRoot! Arrggghhhhh... */
618#endif
619
620 if (pixel->matte != MagickFalse)
621 pixel->opacity += weight*(MagickRealType) pixels->opacity;
622 divisor_m += weight;
623
624 if (pixel->matte != MagickFalse)
625 weight *= QuantumScale*((MagickRealType)(QuantumRange-pixels->opacity));
626 pixel->red += weight*(MagickRealType) pixels->red;
627 pixel->green += weight*(MagickRealType) pixels->green;
628 pixel->blue += weight*(MagickRealType) pixels->blue;
629 if (pixel->colorspace == CMYKColorspace)
630 pixel->index += weight*(MagickRealType) (*indexes);
631 divisor_c += weight;
632
633 hit++;
634#if DEBUG_HIT_MISS
635 /* mark the pixel according to hit/miss of the ellipse */
636 (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
637 (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
638 (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
639 (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
640 } else {
641 (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
642 (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
643 (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
644 (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
645 }
646 uu++;
647#else
648 }
649#endif
650 pixels++;
651 indexes++;
652 Q += DQ;
653 DQ += DDQ;
654 }
655 }
656#if DEBUG_ELLIPSE
657 (void) FormatLocaleFile(stderr, "Hit=%ld; Total=%ld;\n", (long)hit, (long)uw*(v2-v1) );
658#endif
659
660 /*
661 Result sanity check -- this should NOT happen
662 */
663 if ( hit == 0 || divisor_m <= MagickEpsilon || divisor_c <= MagickEpsilon ) {
664 /* not enough pixels, or bad weighting in resampling,
665 resort to direct interpolation */
666#if DEBUG_NO_PIXEL_HIT
667 pixel->opacity = pixel->red = pixel->green = pixel->blue = 0;
668 pixel->red = QuantumRange; /* show pixels for which EWA fails */
669#else
670 status=InterpolateMagickPixelPacket(resample_filter->image,
671 resample_filter->view,resample_filter->interpolate,u0,v0,pixel,
672 resample_filter->exception);
673#endif
674 return status;
675 }
676
677 /*
678 Finalize results of resampling
679 */
680 divisor_m = 1.0/divisor_m;
681 if (pixel->matte != MagickFalse)
682 pixel->opacity = (MagickRealType) ClampToQuantum(divisor_m*pixel->opacity);
683 divisor_c = 1.0/divisor_c;
684 pixel->red = (MagickRealType) ClampToQuantum(divisor_c*pixel->red);
685 pixel->green = (MagickRealType) ClampToQuantum(divisor_c*pixel->green);
686 pixel->blue = (MagickRealType) ClampToQuantum(divisor_c*pixel->blue);
687 if (pixel->colorspace == CMYKColorspace)
688 pixel->index = (MagickRealType) ClampToQuantum(divisor_c*pixel->index);
689 return(MagickTrue);
690}
691␌
692#if EWA && EWA_CLAMP
693/*
694%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
695% %
696% %
697% %
698- C l a m p U p A x e s %
699% %
700% %
701% %
702%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
703%
704% ClampUpAxes() function converts the input vectors into a major and
705% minor axis unit vectors, and their magnitude. This allows us to
706% ensure that the ellipse generated is never smaller than the unit
707% circle and thus never too small for use in EWA resampling.
708%
709% This purely mathematical 'magic' was provided by Professor Nicolas
710% Robidoux and his Masters student Chantal Racette.
711%
712% Reference: "We Recommend Singular Value Decomposition", David Austin
713% http://www.ams.org/samplings/feature-column/fcarc-svd
714%
715% By generating major and minor axis vectors, we can actually use the
716% ellipse in its "canonical form", by remapping the dx,dy of the
717% sampled point into distances along the major and minor axis unit
718% vectors.
719%
720% Reference: http://en.wikipedia.org/wiki/Ellipse#Canonical_form
721*/
722static inline void ClampUpAxes(const double dux,
723 const double dvx,
724 const double duy,
725 const double dvy,
726 double *major_mag,
727 double *minor_mag,
728 double *major_unit_x,
729 double *major_unit_y,
730 double *minor_unit_x,
731 double *minor_unit_y)
732{
733 /*
734 * ClampUpAxes takes an input 2x2 matrix
735 *
736 * [ a b ] = [ dux duy ]
737 * [ c d ] = [ dvx dvy ]
738 *
739 * and computes from it the major and minor axis vectors [major_x,
740 * major_y] and [minor_x,minor_y] of the smallest ellipse containing
741 * both the unit disk and the ellipse which is the image of the unit
742 * disk by the linear transformation
743 *
744 * [ dux duy ] [S] = [s]
745 * [ dvx dvy ] [T] = [t]
746 *
747 * (The vector [S,T] is the difference between a position in output
748 * space and [X,Y]; the vector [s,t] is the difference between a
749 * position in input space and [x,y].)
750 */
751 /*
752 * Output:
753 *
754 * major_mag is the half-length of the major axis of the "new"
755 * ellipse.
756 *
757 * minor_mag is the half-length of the minor axis of the "new"
758 * ellipse.
759 *
760 * major_unit_x is the x-coordinate of the major axis direction vector
761 * of both the "old" and "new" ellipses.
762 *
763 * major_unit_y is the y-coordinate of the major axis direction vector.
764 *
765 * minor_unit_x is the x-coordinate of the minor axis direction vector.
766 *
767 * minor_unit_y is the y-coordinate of the minor axis direction vector.
768 *
769 * Unit vectors are useful for computing projections, in particular,
770 * to compute the distance between a point in output space and the
771 * center of a unit disk in output space, using the position of the
772 * corresponding point [s,t] in input space. Following the clamping,
773 * the square of this distance is
774 *
775 * ( ( s * major_unit_x + t * major_unit_y ) / major_mag )^2
776 * +
777 * ( ( s * minor_unit_x + t * minor_unit_y ) / minor_mag )^2
778 *
779 * If such distances will be computed for many [s,t]'s, it makes
780 * sense to actually compute the reciprocal of major_mag and
781 * minor_mag and multiply them by the above unit lengths.
782 *
783 * Now, if you want to modify the input pair of tangent vectors so
784 * that it defines the modified ellipse, all you have to do is set
785 *
786 * newdux = major_mag * major_unit_x
787 * newdvx = major_mag * major_unit_y
788 * newduy = minor_mag * minor_unit_x = minor_mag * -major_unit_y
789 * newdvy = minor_mag * minor_unit_y = minor_mag * major_unit_x
790 *
791 * and use these tangent vectors as if they were the original ones.
792 * Usually, this is a drastic change in the tangent vectors even if
793 * the singular values are not clamped; for example, the minor axis
794 * vector always points in a direction which is 90 degrees
795 * counterclockwise from the direction of the major axis vector.
796 */
797 /*
798 * Discussion:
799 *
800 * GOAL: Fix things so that the pullback, in input space, of a disk
801 * of radius r in output space is an ellipse which contains, at
802 * least, a disc of radius r. (Make this hold for any r>0.)
803 *
804 * ESSENCE OF THE METHOD: Compute the product of the first two
805 * factors of an SVD of the linear transformation defining the
806 * ellipse and make sure that both its columns have norm at least 1.
807 * Because rotations and reflexions map disks to themselves, it is
808 * not necessary to compute the third (rightmost) factor of the SVD.
809 *
810 * DETAILS: Find the singular values and (unit) left singular
811 * vectors of Jinv, clampling up the singular values to 1, and
812 * multiply the unit left singular vectors by the new singular
813 * values in order to get the minor and major ellipse axis vectors.
814 *
815 * Image resampling context:
816 *
817 * The Jacobian matrix of the transformation at the output point
818 * under consideration is defined as follows:
819 *
820 * Consider the transformation (x,y) -> (X,Y) from input locations
821 * to output locations. (Anthony Thyssen, elsewhere in resample.c,
822 * uses the notation (u,v) -> (x,y).)
823 *
824 * The Jacobian matrix of the transformation at (x,y) is equal to
825 *
826 * J = [ A, B ] = [ dX/dx, dX/dy ]
827 * [ C, D ] [ dY/dx, dY/dy ]
828 *
829 * that is, the vector [A,C] is the tangent vector corresponding to
830 * input changes in the horizontal direction, and the vector [B,D]
831 * is the tangent vector corresponding to input changes in the
832 * vertical direction.
833 *
834 * In the context of resampling, it is natural to use the inverse
835 * Jacobian matrix Jinv because resampling is generally performed by
836 * pulling pixel locations in the output image back to locations in
837 * the input image. Jinv is
838 *
839 * Jinv = [ a, b ] = [ dx/dX, dx/dY ]
840 * [ c, d ] [ dy/dX, dy/dY ]
841 *
842 * Note: Jinv can be computed from J with the following matrix
843 * formula:
844 *
845 * Jinv = 1/(A*D-B*C) [ D, -B ]
846 * [ -C, A ]
847 *
848 * What we do is modify Jinv so that it generates an ellipse which
849 * is as close as possible to the original but which contains the
850 * unit disk. This can be accomplished as follows:
851 *
852 * Let
853 *
854 * Jinv = U Sigma V^T
855 *
856 * be an SVD decomposition of Jinv. (The SVD is not unique, but the
857 * final ellipse does not depend on the particular SVD.)
858 *
859 * We could clamp up the entries of the diagonal matrix Sigma so
860 * that they are at least 1, and then set
861 *
862 * Jinv = U newSigma V^T.
863 *
864 * However, we do not need to compute V for the following reason:
865 * V^T is an orthogonal matrix (that is, it represents a combination
866 * of rotations and reflexions) so that it maps the unit circle to
867 * itself. For this reason, the exact value of V does not affect the
868 * final ellipse, and we can choose V to be the identity
869 * matrix. This gives
870 *
871 * Jinv = U newSigma.
872 *
873 * In the end, we return the two diagonal entries of newSigma
874 * together with the two columns of U.
875 */
876 /*
877 * ClampUpAxes was written by Nicolas Robidoux and Chantal Racette
878 * of Laurentian University with insightful suggestions from Anthony
879 * Thyssen and funding from the National Science and Engineering
880 * Research Council of Canada. It is distinguished from its
881 * predecessors by its efficient handling of degenerate cases.
882 *
883 * The idea of clamping up the EWA ellipse's major and minor axes so
884 * that the result contains the reconstruction kernel filter support
885 * is taken from Andreas Gustafsson's Masters thesis "Interactive
886 * Image Warping", Helsinki University of Technology, Faculty of
887 * Information Technology, 59 pages, 1993 (see Section 3.6).
888 *
889 * The use of the SVD to clamp up the singular values of the
890 * Jacobian matrix of the pullback transformation for EWA resampling
891 * is taken from the astrophysicist Craig DeForest. It is
892 * implemented in his PDL::Transform code (PDL = Perl Data
893 * Language).
894 */
895 const double a = dux;
896 const double b = duy;
897 const double c = dvx;
898 const double d = dvy;
899 /*
900 * n is the matrix Jinv * transpose(Jinv). Eigenvalues of n are the
901 * squares of the singular values of Jinv.
902 */
903 const double aa = a*a;
904 const double bb = b*b;
905 const double cc = c*c;
906 const double dd = d*d;
907 /*
908 * Eigenvectors of n are left singular vectors of Jinv.
909 */
910 const double n11 = aa+bb;
911 const double n12 = a*c+b*d;
912 const double n21 = n12;
913 const double n22 = cc+dd;
914 const double det = a*d-b*c;
915 const double twice_det = det+det;
916 const double frobenius_squared = n11+n22;
917 const double discriminant =
918 (frobenius_squared+twice_det)*(frobenius_squared-twice_det);
919 /*
920 * In exact arithmetic, discriminant can't be negative. In floating
921 * point, it can, because of the bad conditioning of SVD
922 * decompositions done through the associated normal matrix.
923 */
924 const double sqrt_discriminant =
925 sqrt(discriminant > 0.0 ? discriminant : 0.0);
926 /*
927 * s1 is the largest singular value of the inverse Jacobian
928 * matrix. In other words, its reciprocal is the smallest singular
929 * value of the Jacobian matrix itself.
930 * If s1 = 0, both singular values are 0, and any orthogonal pair of
931 * left and right factors produces a singular decomposition of Jinv.
932 */
933 /*
934 * Initially, we only compute the squares of the singular values.
935 */
936 const double s1s1 = 0.5*(frobenius_squared+sqrt_discriminant);
937 /*
938 * s2 the smallest singular value of the inverse Jacobian
939 * matrix. Its reciprocal is the largest singular value of the
940 * Jacobian matrix itself.
941 */
942 const double s2s2 = 0.5*(frobenius_squared-sqrt_discriminant);
943 const double s1s1minusn11 = s1s1-n11;
944 const double s1s1minusn22 = s1s1-n22;
945 /*
946 * u1, the first column of the U factor of a singular decomposition
947 * of Jinv, is a (non-normalized) left singular vector corresponding
948 * to s1. It has entries u11 and u21. We compute u1 from the fact
949 * that it is an eigenvector of n corresponding to the eigenvalue
950 * s1^2.
951 */
952 const double s1s1minusn11_squared = s1s1minusn11*s1s1minusn11;
953 const double s1s1minusn22_squared = s1s1minusn22*s1s1minusn22;
954 /*
955 * The following selects the largest row of n-s1^2 I as the one
956 * which is used to find the eigenvector. If both s1^2-n11 and
957 * s1^2-n22 are zero, n-s1^2 I is the zero matrix. In that case,
958 * any vector is an eigenvector; in addition, norm below is equal to
959 * zero, and, in exact arithmetic, this is the only case in which
960 * norm = 0. So, setting u1 to the simple but arbitrary vector [1,0]
961 * if norm = 0 safely takes care of all cases.
962 */
963 const double temp_u11 =
964 ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? n12 : s1s1minusn22 );
965 const double temp_u21 =
966 ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? s1s1minusn11 : n21 );
967 const double norm = sqrt(temp_u11*temp_u11+temp_u21*temp_u21);
968 /*
969 * Finalize the entries of first left singular vector (associated
970 * with the largest singular value).
971 */
972 const double u11 = ( (norm>0.0) ? temp_u11/norm : 1.0 );
973 const double u21 = ( (norm>0.0) ? temp_u21/norm : 0.0 );
974 /*
975 * Clamp the singular values up to 1.
976 */
977 *major_mag = ( (s1s1<=1.0) ? 1.0 : sqrt(s1s1) );
978 *minor_mag = ( (s2s2<=1.0) ? 1.0 : sqrt(s2s2) );
979 /*
980 * Return the unit major and minor axis direction vectors.
981 */
982 *major_unit_x = u11;
983 *major_unit_y = u21;
984 *minor_unit_x = -u21;
985 *minor_unit_y = u11;
986}
987␌
988#endif
989/*
990%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
991% %
992% %
993% %
994% S c a l e R e s a m p l e F i l t e r %
995% %
996% %
997% %
998%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
999%
1000% ScaleResampleFilter() does all the calculations needed to resample an image
1001% at a specific scale, defined by two scaling vectors. This not using
1002% a orthogonal scaling, but two distorted scaling vectors, to allow the
1003% generation of a angled ellipse.
1004%
1005% As only two derivative scaling vectors are used the center of the ellipse
1006% must be the center of the lookup. That is any curvature that the
1007% distortion may produce is discounted.
1008%
1009% The input vectors are produced by either finding the derivatives of the
1010% distortion function, or the partial derivatives from a distortion mapping.
1011% They do not need to be the orthogonal dx,dy scaling vectors, but can be
1012% calculated from other derivatives. For example you could use dr,da/r
1013% polar coordinate vector scaling vectors
1014%
1015% If u,v = DistortEquation(x,y) OR u = Fu(x,y); v = Fv(x,y)
1016% Then the scaling vectors are determined from the derivatives...
1017% du/dx, dv/dx and du/dy, dv/dy
1018% If the resulting scaling vectors is orthogonally aligned then...
1019% dv/dx = 0 and du/dy = 0
1020% Producing an orthogonally aligned ellipse in source space for the area to
1021% be resampled.
1022%
1023% Note that scaling vectors are different to argument order. Argument order
1024% is the general order the derivatives are extracted from the distortion
1025% equations, and not the scaling vectors. As such the middle two values
1026% may be swapped from what you expect. Caution is advised.
1027%
1028% WARNING: It is assumed that any SetResampleFilter() method call will
1029% always be performed before the ScaleResampleFilter() method, so that the
1030% size of the ellipse will match the support for the resampling filter being
1031% used.
1032%
1033% The format of the ScaleResampleFilter method is:
1034%
1035% void ScaleResampleFilter(const ResampleFilter *resample_filter,
1036% const double dux,const double duy,const double dvx,const double dvy)
1037%
1038% A description of each parameter follows:
1039%
1040% o resample_filter: the resampling information defining the
1041% image being resampled
1042%
1043% o dux,duy,dvx,dvy:
1044% The derivatives or scaling vectors defining the EWA ellipse.
1045% NOTE: watch the order, which is based on the order derivatives
1046% are usually determined from distortion equations (see above).
1047% The middle two values may need to be swapped if you are thinking
1048% in terms of scaling vectors.
1049%
1050*/
1051MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
1052 const double dux,const double duy,const double dvx,const double dvy)
1053{
1054 double A,B,C,F;
1055
1056 assert(resample_filter != (ResampleFilter *) NULL);
1057 assert(resample_filter->signature == MagickCoreSignature);
1058
1059 resample_filter->limit_reached = MagickFalse;
1060
1061 /* A 'point' filter forces use of interpolation instead of area sampling */
1062 if ( resample_filter->filter == PointFilter )
1063 return; /* EWA turned off - nothing to do */
1064
1065#if DEBUG_ELLIPSE
1066 (void) FormatLocaleFile(stderr, "# -----\n" );
1067 (void) FormatLocaleFile(stderr, "dux=%lf; dvx=%lf; duy=%lf; dvy=%lf;\n",
1068 dux, dvx, duy, dvy);
1069#endif
1070
1071 /* Find Ellipse Coefficients such that
1072 A*u^2 + B*u*v + C*v^2 = F
1073 With u,v relative to point around which we are resampling.
1074 And the given scaling dx,dy vectors in u,v space
1075 du/dx,dv/dx and du/dy,dv/dy
1076 */
1077#if EWA
1078 /* Direct conversion of derivatives into elliptical coefficients
1079 However when magnifying images, the scaling vectors will be small
1080 resulting in a ellipse that is too small to sample properly.
1081 As such we need to clamp the major/minor axis to a minimum of 1.0
1082 to prevent it getting too small.
1083 */
1084#if EWA_CLAMP
1085 { double major_mag,
1086 minor_mag,
1087 major_x,
1088 major_y,
1089 minor_x,
1090 minor_y;
1091
1092 ClampUpAxes(dux,dvx,duy,dvy, &major_mag, &minor_mag,
1093 &major_x, &major_y, &minor_x, &minor_y);
1094 major_x *= major_mag; major_y *= major_mag;
1095 minor_x *= minor_mag; minor_y *= minor_mag;
1096#if DEBUG_ELLIPSE
1097 (void) FormatLocaleFile(stderr, "major_x=%lf; major_y=%lf; minor_x=%lf; minor_y=%lf;\n",
1098 major_x, major_y, minor_x, minor_y);
1099#endif
1100 A = major_y*major_y+minor_y*minor_y;
1101 B = -2.0*(major_x*major_y+minor_x*minor_y);
1102 C = major_x*major_x+minor_x*minor_x;
1103 F = major_mag*minor_mag;
1104 F *= F; /* square it */
1105 }
1106#else /* raw unclamped EWA */
1107 A = dvx*dvx+dvy*dvy;
1108 B = -2.0*(dux*dvx+duy*dvy);
1109 C = dux*dux+duy*duy;
1110 F = dux*dvy-duy*dvx;
1111 F *= F; /* square it */
1112#endif /* EWA_CLAMP */
1113
1114#else /* HQ_EWA */
1115 /*
1116 This Paul Heckbert's "Higher Quality EWA" formula, from page 60 in his
1117 thesis, which adds a unit circle to the elliptical area so as to do both
1118 Reconstruction and Prefiltering of the pixels in the resampling. It also
1119 means it is always likely to have at least 4 pixels within the area of the
1120 ellipse, for weighted averaging. No scaling will result with F == 4.0 and
1121 a circle of radius 2.0, and F smaller than this means magnification is
1122 being used.
1123
1124 NOTE: This method produces a very blurry result at near unity scale while
1125 producing perfect results for strong minification and magnifications.
1126
1127 However filter support is fixed to 2.0 (no good for Windowed Sinc filters)
1128 */
1129 A = dvx*dvx+dvy*dvy+1;
1130 B = -2.0*(dux*dvx+duy*dvy);
1131 C = dux*dux+duy*duy+1;
1132 F = A*C - B*B/4;
1133#endif
1134
1135#if DEBUG_ELLIPSE
1136 (void) FormatLocaleFile(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
1137
1138 /* Figure out the various information directly about the ellipse.
1139 This information currently not needed at this time, but may be
1140 needed later for better limit determination.
1141
1142 It is also good to have as a record for future debugging
1143 */
1144 { double alpha, beta, gamma, Major, Minor;
1145 double Eccentricity, Ellipse_Area, Ellipse_Angle;
1146
1147 alpha = A+C;
1148 beta = A-C;
1149 gamma = sqrt(beta*beta + B*B );
1150
1151 if ( alpha - gamma <= MagickEpsilon )
1152 Major= MagickMaximumValue;
1153 else
1154 Major= sqrt(2*F/(alpha - gamma));
1155 Minor = sqrt(2*F/(alpha + gamma));
1156
1157 (void) FormatLocaleFile(stderr, "# Major=%lf; Minor=%lf\n", Major, Minor );
1158
1159 /* other information about ellipse include... */
1160 Eccentricity = Major/Minor;
1161 Ellipse_Area = MagickPI*Major*Minor;
1162 Ellipse_Angle = atan2(B, A-C);
1163
1164 (void) FormatLocaleFile(stderr, "# Angle=%lf Area=%lf\n",
1165 (double) RadiansToDegrees(Ellipse_Angle), Ellipse_Area);
1166 }
1167#endif
1168
1169 /* If one or both of the scaling vectors is impossibly large
1170 (producing a very large raw F value), we may as well not bother
1171 doing any form of resampling since resampled area is very large.
1172 In this case some alternative means of pixel sampling, such as
1173 the average of the whole image is needed to get a reasonable
1174 result. Calculate only as needed.
1175 */
1176 if ( (4*A*C - B*B) > MagickMaximumValue ) {
1177 resample_filter->limit_reached = MagickTrue;
1178 return;
1179 }
1180
1181 /* Scale ellipse to match the filters support
1182 (that is, multiply F by the square of the support)
1183 Simpler to just multiply it by the support twice!
1184 */
1185 F *= resample_filter->support;
1186 F *= resample_filter->support;
1187
1188 /* Orthogonal bounds of the ellipse */
1189 resample_filter->Ulimit = sqrt(C*F/(A*C-0.25*B*B));
1190 resample_filter->Vlimit = sqrt(A*F/(A*C-0.25*B*B));
1191
1192 /* Horizontally aligned parallelogram fitted to Ellipse */
1193 resample_filter->Uwidth = sqrt(F/A); /* Half of the parallelogram width */
1194 resample_filter->slope = -B/(2.0*A); /* Reciprocal slope of the parallelogram */
1195
1196#if DEBUG_ELLIPSE
1197 (void) FormatLocaleFile(stderr, "Ulimit=%lf; Vlimit=%lf; UWidth=%lf; Slope=%lf;\n",
1198 resample_filter->Ulimit, resample_filter->Vlimit,
1199 resample_filter->Uwidth, resample_filter->slope );
1200#endif
1201
1202 /* Check the absolute area of the parallelogram involved.
1203 * This limit needs more work, as it is too slow for larger images
1204 * with tiled views of the horizon.
1205 */
1206 if ( (resample_filter->Uwidth * resample_filter->Vlimit)
1207 > (4.0*resample_filter->image_area)) {
1208 resample_filter->limit_reached = MagickTrue;
1209 return;
1210 }
1211
1212 /* Scale ellipse formula to directly index the Filter Lookup Table */
1213 { double scale;
1214#if FILTER_LUT
1215 /* scale so that F = WLUT_WIDTH; -- hardcoded */
1216 scale=(double) WLUT_WIDTH*PerceptibleReciprocal(F);
1217#else
1218 /* scale so that F = resample_filter->F (support^2) */
1219 scale=resample_filter->F*PerceptibleReciprocal(F);
1220#endif
1221 resample_filter->A = A*scale;
1222 resample_filter->B = B*scale;
1223 resample_filter->C = C*scale;
1224 }
1225}
1226␌
1227/*
1228%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1229% %
1230% %
1231% %
1232% S e t R e s a m p l e F i l t e r %
1233% %
1234% %
1235% %
1236%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1237%
1238% SetResampleFilter() set the resampling filter lookup table based on a
1239% specific filter. Note that the filter is used as a radial filter not as a
1240% two pass orthogonally aligned resampling filter.
1241%
1242% The format of the SetResampleFilter method is:
1243%
1244% void SetResampleFilter(ResampleFilter *resample_filter,
1245% const FilterTypes filter,const double blur)
1246%
1247% A description of each parameter follows:
1248%
1249% o resample_filter: resampling information structure
1250%
1251% o filter: the resize filter for elliptical weighting LUT
1252%
1253% o blur: filter blur factor (radial scaling) for elliptical weighting LUT
1254%
1255*/
1256MagickExport void SetResampleFilter(ResampleFilter *resample_filter,
1257 const FilterTypes filter,const double blur)
1258{
1259 ResizeFilter
1260 *resize_filter;
1261
1262 assert(resample_filter != (ResampleFilter *) NULL);
1263 assert(resample_filter->signature == MagickCoreSignature);
1264
1265 resample_filter->do_interpolate = MagickFalse;
1266 resample_filter->filter = filter;
1267
1268 /* Default cylindrical filter is a Cubic Keys filter */
1269 if ( filter == UndefinedFilter )
1270 resample_filter->filter = RobidouxFilter;
1271
1272 if ( resample_filter->filter == PointFilter ) {
1273 resample_filter->do_interpolate = MagickTrue;
1274 return; /* EWA turned off - nothing more to do */
1275 }
1276
1277 resize_filter = AcquireResizeFilter(resample_filter->image,
1278 resample_filter->filter,blur,MagickTrue,resample_filter->exception);
1279 if (resize_filter == (ResizeFilter *) NULL) {
1280 (void) ThrowMagickException(resample_filter->exception,GetMagickModule(),
1281 ModuleError, "UnableToSetFilteringValue",
1282 "Fall back to Interpolated 'Point' filter");
1283 resample_filter->filter = PointFilter;
1284 resample_filter->do_interpolate = MagickTrue;
1285 return; /* EWA turned off - nothing more to do */
1286 }
1287
1288 /* Get the practical working support for the filter,
1289 * after any API call blur factors have been accounted for.
1290 */
1291#if EWA
1292 resample_filter->support = GetResizeFilterSupport(resize_filter);
1293#else
1294 resample_filter->support = 2.0; /* fixed support size for HQ-EWA */
1295#endif
1296
1297#if FILTER_LUT
1298 /* Fill the LUT with the weights from the selected filter function */
1299 { int
1300 Q;
1301 double
1302 r_scale;
1303
1304 /* Scale radius so the filter LUT covers the full support range */
1305 r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
1306 for(Q=0; Q<WLUT_WIDTH; Q++)
1307 resample_filter->filter_lut[Q] = (double)
1308 GetResizeFilterWeight(resize_filter,sqrt((double)Q)*r_scale);
1309
1310 /* finished with the resize filter */
1311 resize_filter = DestroyResizeFilter(resize_filter);
1312 }
1313#else
1314 /* save the filter and the scaled ellipse bounds needed for filter */
1315 resample_filter->filter_def = resize_filter;
1316 resample_filter->F = resample_filter->support*resample_filter->support;
1317#endif
1318
1319 /*
1320 Adjust the scaling of the default unit circle
1321 This assumes that any real scaling changes will always
1322 take place AFTER the filter method has been initialized.
1323 */
1324 ScaleResampleFilter(resample_filter, 1.0, 0.0, 0.0, 1.0);
1325
1326#if 0
1327 /*
1328 This is old code kept as a reference only. Basically it generates
1329 a Gaussian bell curve, with sigma = 0.5 if the support is 2.0
1330
1331 Create Normal Gaussian 2D Filter Weighted Lookup Table.
1332 A normal EWA guassual lookup would use exp(Q*ALPHA)
1333 where Q = distance squared from 0.0 (center) to 1.0 (edge)
1334 and ALPHA = -4.0*ln(2.0) ==> -2.77258872223978123767
1335 The table is of length 1024, and equates to support radius of 2.0
1336 thus needs to be scaled by ALPHA*4/1024 and any blur factor squared
1337
1338 The it comes from reference code provided by Fred Weinhaus.
1339 */
1340 r_scale = -2.77258872223978123767/(WLUT_WIDTH*blur*blur);
1341 for(Q=0; Q<WLUT_WIDTH; Q++)
1342 resample_filter->filter_lut[Q] = exp((double)Q*r_scale);
1343 resample_filter->support = WLUT_WIDTH;
1344#endif
1345
1346#if FILTER_LUT
1347#if defined(MAGICKCORE_OPENMP_SUPPORT)
1348 #pragma omp single
1349#endif
1350 {
1351 if (IsMagickTrue(GetImageArtifact(resample_filter->image,
1352 "resample:verbose")) )
1353 {
1354 int
1355 Q;
1356 double
1357 r_scale;
1358
1359 /* Debug output of the filter weighting LUT
1360 Gnuplot the LUT data, the x scale index has been adjusted
1361 plot [0:2][-.2:1] "lut.dat" with lines
1362 The filter values should be normalized for comparison
1363 */
1364 printf("#\n");
1365 printf("# Resampling Filter LUT (%d values) for '%s' filter\n",
1366 WLUT_WIDTH, CommandOptionToMnemonic(MagickFilterOptions,
1367 resample_filter->filter) );
1368 printf("#\n");
1369 printf("# Note: values in table are using a squared radius lookup.\n");
1370 printf("# As such its distribution is not uniform.\n");
1371 printf("#\n");
1372 printf("# The X value is the support distance for the Y weight\n");
1373 printf("# so you can use gnuplot to plot this cylindrical filter\n");
1374 printf("# plot [0:2][-.2:1] \"lut.dat\" with lines\n");
1375 printf("#\n");
1376
1377 /* Scale radius so the filter LUT covers the full support range */
1378 r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
1379 for(Q=0; Q<WLUT_WIDTH; Q++)
1380 printf("%8.*g %.*g\n",
1381 GetMagickPrecision(),sqrt((double)Q)*r_scale,
1382 GetMagickPrecision(),resample_filter->filter_lut[Q] );
1383 printf("\n\n"); /* generate a 'break' in gnuplot if multiple outputs */
1384 }
1385 /* Output the above once only for each image, and each setting
1386 (void) DeleteImageArtifact(resample_filter->image,"resample:verbose");
1387 */
1388 }
1389#endif /* FILTER_LUT */
1390 return;
1391}
1392␌
1393/*
1394%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1395% %
1396% %
1397% %
1398% 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 %
1399% %
1400% %
1401% %
1402%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1403%
1404% SetResampleFilterInterpolateMethod() sets the resample filter interpolation
1405% method.
1406%
1407% The format of the SetResampleFilterInterpolateMethod method is:
1408%
1409% MagickBooleanType SetResampleFilterInterpolateMethod(
1410% ResampleFilter *resample_filter,const InterpolateMethod method)
1411%
1412% A description of each parameter follows:
1413%
1414% o resample_filter: the resample filter.
1415%
1416% o method: the interpolation method.
1417%
1418*/
1419MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
1420 ResampleFilter *resample_filter,const InterpolatePixelMethod method)
1421{
1422 assert(resample_filter != (ResampleFilter *) NULL);
1423 assert(resample_filter->signature == MagickCoreSignature);
1424 assert(resample_filter->image != (Image *) NULL);
1425 if (IsEventLogging() != MagickFalse)
1426 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1427 resample_filter->image->filename);
1428 resample_filter->interpolate=method;
1429 return(MagickTrue);
1430}
1431␌
1432/*
1433%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1434% %
1435% %
1436% %
1437% 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 %
1438% %
1439% %
1440% %
1441%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1442%
1443% SetResampleFilterVirtualPixelMethod() changes the virtual pixel method
1444% associated with the specified resample filter.
1445%
1446% The format of the SetResampleFilterVirtualPixelMethod method is:
1447%
1448% MagickBooleanType SetResampleFilterVirtualPixelMethod(
1449% ResampleFilter *resample_filter,const VirtualPixelMethod method)
1450%
1451% A description of each parameter follows:
1452%
1453% o resample_filter: the resample filter.
1454%
1455% o method: the virtual pixel method.
1456%
1457*/
1458MagickExport MagickBooleanType SetResampleFilterVirtualPixelMethod(
1459 ResampleFilter *resample_filter,const VirtualPixelMethod method)
1460{
1461 assert(resample_filter != (ResampleFilter *) NULL);
1462 assert(resample_filter->signature == MagickCoreSignature);
1463 assert(resample_filter->image != (Image *) NULL);
1464 if (IsEventLogging() != MagickFalse)
1465 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1466 resample_filter->image->filename);
1467 resample_filter->virtual_pixel=method;
1468 if (method != UndefinedVirtualPixelMethod)
1469 (void) SetCacheViewVirtualPixelMethod(resample_filter->view,method);
1470 return(MagickTrue);
1471}
_CacheView
Definition cache-view.c:66
_ExceptionInfo
Definition exception.h:103
_Image
Definition image.h:134
_MagickPixelPacket
Definition pixel.h:105
_PixelPacket
Definition pixel.h:132
_ResampleFilter
Definition resample.c:92
_ResizeFilter
Definition resize.c:88