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