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