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