MagickCore  6.9.13-9
Convert, Edit, Or Compose Bitmap Images
resample.c
1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % %
4 % %
5 % %
6 % RRRR EEEEE SSSSS AAA M M PPPP L EEEEE %
7 % R R E SS A A MM MM P P L E %
8 % RRRR EEE SSS AAAAA M M M PPPP L EEE %
9 % R R E SS A A M M P L E %
10 % R R EEEEE SSSSS A A M M P LLLLL EEEEE %
11 % %
12 % %
13 % MagickCore Pixel Resampling Methods %
14 % %
15 % Software Design %
16 % Cristy %
17 % Anthony Thyssen %
18 % August 2007 %
19 % %
20 % %
21 % Copyright 1999 ImageMagick Studio LLC, a non-profit organization %
22 % dedicated to making software imaging solutions freely available. %
23 % %
24 % You may not use this file except in compliance with the License. You may %
25 % obtain a copy of the License at %
26 % %
27 % https://imagemagick.org/script/license.php %
28 % %
29 % Unless required by applicable law or agreed to in writing, software %
30 % distributed under the License is distributed on an "AS IS" BASIS, %
31 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
32 % See the License for the specific language governing permissions and %
33 % limitations under the License. %
34 % %
35 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
36 %
37 %
38 */
39 ␌
40 /*
41  Include declarations.
42 */
43 #include "magick/studio.h"
44 #include "magick/artifact.h"
45 #include "magick/color-private.h"
46 #include "magick/cache.h"
47 #include "magick/draw.h"
48 #include "magick/exception-private.h"
49 #include "magick/gem.h"
50 #include "magick/image.h"
51 #include "magick/image-private.h"
52 #include "magick/log.h"
53 #include "magick/magick.h"
54 #include "magick/memory_.h"
55 #include "magick/pixel.h"
56 #include "magick/pixel-private.h"
57 #include "magick/quantum.h"
58 #include "magick/random_.h"
59 #include "magick/resample.h"
60 #include "magick/resize.h"
61 #include "magick/resize-private.h"
62 #include "magick/resource_.h"
63 #include "magick/transform.h"
64 #include "magick/signature-private.h"
65 #include "magick/token.h"
66 #include "magick/utility.h"
67 #include "magick/option.h"
68 /*
69  EWA Resampling Options
70 */
71 
72 /* select ONE resampling method */
73 #define EWA 1 /* Normal EWA handling - raw or clamped */
74  /* if 0 then use "High Quality EWA" */
75 #define EWA_CLAMP 1 /* EWA Clamping from Nicolas Robidoux */
76 
77 #define FILTER_LUT 1 /* Use a LUT rather then direct filter calls */
78 
79 /* output debugging information */
80 #define DEBUG_ELLIPSE 0 /* output ellipse info for debug */
81 #define DEBUG_HIT_MISS 0 /* output hit/miss pixels (as gnuplot commands) */
82 #define DEBUG_NO_PIXEL_HIT 0 /* Make pixels that fail to hit anything - RED */
83 
84 #if ! FILTER_DIRECT
85 #define WLUT_WIDTH 1024 /* size of the filter cache */
86 #endif
87 
88 /*
89  Typedef declarations.
90 */
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 elliptical area being resampled around center point */
128  double
129  A, B, C,
130  Vlimit, Ulimit, Uwidth, slope;
131 
132 #if FILTER_LUT
133  /* LUT of weights for filtered average in elliptical area */
134  double
135  filter_lut[WLUT_WIDTH];
136 #else
137  /* Use a Direct call to the filter functions */
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 efficient manner possible for the
302 % VirtualPixelMethod setting.
303 %
304 % The format of the ResamplePixelColor method is:
305 %
306 % MagickBooleanType ResamplePixelColor(ResampleFilter *resample_filter,
307 % const double u0,const double v0,MagickPixelPacket *pixel)
308 %
309 % A description of each parameter follows:
310 %
311 % o resample_filter: the resample filter.
312 %
313 % o u0,v0: A double representing the center of the area to resample,
314 % The distortion transformed x,y coordinate.
315 %
316 % o pixel: the resampled pixel is returned here.
317 %
318 */
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  (MagickRealType) resample_filter->average_pixel.opacity);
516  resample_filter->average_pixel.red *= weight;
517  resample_filter->average_pixel.green *= weight;
518  resample_filter->average_pixel.blue *= weight;
519  divisor_c = weight;
520 
521  /* background color */
522  weight = QuantumScale*((MagickRealType) QuantumRange-
523  (MagickRealType)
524  resample_filter->image->background_color.opacity);
525  resample_filter->average_pixel.red += weight*(MagickRealType)
526  resample_filter->image->background_color.red;
527  resample_filter->average_pixel.green += weight*(MagickRealType)
528  resample_filter->image->background_color.green;
529  resample_filter->average_pixel.blue += weight*(MagickRealType)
530  resample_filter->image->background_color.blue;
531  resample_filter->average_pixel.opacity += (MagickRealType)
532  resample_filter->image->background_color.opacity;
533  divisor_c += weight;
534 
535  /* alpha blend */
536  resample_filter->average_pixel.red /= divisor_c;
537  resample_filter->average_pixel.green /= divisor_c;
538  resample_filter->average_pixel.blue /= divisor_c;
539  resample_filter->average_pixel.opacity /= 2; /* 50% blend */
540 
541  }
542  }
543  *pixel=resample_filter->average_pixel;
544  break;
545  }
546  return(status);
547  }
548 
549  /*
550  Initialize weighted average data collection
551  */
552  hit = 0;
553  divisor_c = 0.0;
554  divisor_m = 0.0;
555  pixel->red = pixel->green = pixel->blue = 0.0;
556  if (pixel->matte != MagickFalse) pixel->opacity = 0.0;
557  if (pixel->colorspace == CMYKColorspace) pixel->index = 0.0;
558 
559  /*
560  Determine the parallelogram bounding box fitted to the ellipse
561  centered at u0,v0. This area is bounding by the lines...
562  */
563  v1 = (ssize_t)ceil(v0 - resample_filter->Vlimit); /* range of scan lines */
564  v2 = (ssize_t)floor(v0 + resample_filter->Vlimit);
565 
566  /* scan line start and width across the parallelogram */
567  u1 = u0 + (v1-v0)*resample_filter->slope - resample_filter->Uwidth;
568  uw = (ssize_t)(2.0*resample_filter->Uwidth)+1;
569 
570 #if DEBUG_ELLIPSE
571  (void) FormatLocaleFile(stderr, "v1=%ld; v2=%ld\n", (long)v1, (long)v2);
572  (void) FormatLocaleFile(stderr, "u1=%ld; uw=%ld\n", (long)u1, (long)uw);
573 #else
574 # define DEBUG_HIT_MISS 0 /* only valid if DEBUG_ELLIPSE is enabled */
575 #endif
576 
577  /*
578  Do weighted resampling of all pixels, within the scaled ellipse,
579  bound by a Parallelogram fitted to the ellipse.
580  */
581  DDQ = 2*resample_filter->A;
582  for( v=v1; v<=v2; v++ ) {
583 #if DEBUG_HIT_MISS
584  long uu = ceil(u1); /* actual pixel location (for debug only) */
585  (void) FormatLocaleFile(stderr, "# scan line from pixel %ld, %ld\n", (long)uu, (long)v);
586 #endif
587  u = (ssize_t)ceil(u1); /* first pixel in scanline */
588  u1 += resample_filter->slope; /* start of next scan line */
589 
590 
591  /* location of this first pixel, relative to u0,v0 */
592  U = (double)u-u0;
593  V = (double)v-v0;
594 
595  /* Q = ellipse quotient ( if Q<F then pixel is inside ellipse) */
596  Q = (resample_filter->A*U + resample_filter->B*V)*U + resample_filter->C*V*V;
597  DQ = resample_filter->A*(2.0*U+1) + resample_filter->B*V;
598 
599  /* get the scanline of pixels for this v */
600  pixels=GetCacheViewVirtualPixels(resample_filter->view,u,v,(size_t) uw,
601  1,resample_filter->exception);
602  if (pixels == (const PixelPacket *) NULL)
603  return(MagickFalse);
604  indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
605 
606  /* count up the weighted pixel colors */
607  for( u=0; u<uw; u++ ) {
608  weight = 0.0;
609 #if FILTER_LUT
610  /* Note that the ellipse has been pre-scaled so F = WLUT_WIDTH */
611  if ((Q >= 0.0) && ((int) Q < WLUT_WIDTH)) {
612  weight = resample_filter->filter_lut[(int)Q];
613 #else
614  /* Note that the ellipse has been pre-scaled so F = support^2 */
615  if ((Q >= 0.0) && (Q < (double)resample_filter->F)) {
616  weight = GetResizeFilterWeight(resample_filter->filter_def,
617  sqrt(Q)); /* a SquareRoot! Arrggghhhhh... */
618 #endif
619 
620  if (pixel->matte != MagickFalse)
621  pixel->opacity += weight*(MagickRealType) pixels->opacity;
622  divisor_m += weight;
623 
624  if (pixel->matte != MagickFalse)
625  weight *= QuantumScale*((MagickRealType)(QuantumRange-pixels->opacity));
626  pixel->red += weight*(MagickRealType) pixels->red;
627  pixel->green += weight*(MagickRealType) pixels->green;
628  pixel->blue += weight*(MagickRealType) pixels->blue;
629  if (pixel->colorspace == CMYKColorspace)
630  pixel->index += weight*(MagickRealType) (*indexes);
631  divisor_c += weight;
632 
633  hit++;
634 #if DEBUG_HIT_MISS
635  /* mark the pixel according to hit/miss of the ellipse */
636  (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
637  (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
638  (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
639  (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
640  } else {
641  (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
642  (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
643  (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
644  (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
645  }
646  uu++;
647 #else
648  }
649 #endif
650  pixels++;
651  indexes++;
652  Q += DQ;
653  DQ += DDQ;
654  }
655  }
656 #if DEBUG_ELLIPSE
657  (void) FormatLocaleFile(stderr, "Hit=%ld; Total=%ld;\n", (long)hit, (long)uw*(v2-v1) );
658 #endif
659 
660  /*
661  Result sanity check -- this should NOT happen
662  */
663  if ( hit == 0 || divisor_m <= MagickEpsilon || divisor_c <= MagickEpsilon ) {
664  /* not enough pixels, or bad weighting in resampling,
665  resort to direct interpolation */
666 #if DEBUG_NO_PIXEL_HIT
667  pixel->opacity = pixel->red = pixel->green = pixel->blue = 0;
668  pixel->red = QuantumRange; /* show pixels for which EWA fails */
669 #else
670  status=InterpolateMagickPixelPacket(resample_filter->image,
671  resample_filter->view,resample_filter->interpolate,u0,v0,pixel,
672  resample_filter->exception);
673 #endif
674  return status;
675  }
676 
677  /*
678  Finalize results of resampling
679  */
680  divisor_m = 1.0/divisor_m;
681  if (pixel->matte != MagickFalse)
682  pixel->opacity = (MagickRealType) ClampToQuantum(divisor_m*pixel->opacity);
683  divisor_c = 1.0/divisor_c;
684  pixel->red = (MagickRealType) ClampToQuantum(divisor_c*pixel->red);
685  pixel->green = (MagickRealType) ClampToQuantum(divisor_c*pixel->green);
686  pixel->blue = (MagickRealType) ClampToQuantum(divisor_c*pixel->blue);
687  if (pixel->colorspace == CMYKColorspace)
688  pixel->index = (MagickRealType) ClampToQuantum(divisor_c*pixel->index);
689  return(MagickTrue);
690 }
691 ␌
692 #if EWA && EWA_CLAMP
693 /*
694 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
695 % %
696 % %
697 % %
698 - C l a m p U p A x e s %
699 % %
700 % %
701 % %
702 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
703 %
704 % ClampUpAxes() function converts the input vectors into a major and
705 % minor axis unit vectors, and their magnitude. This allows us to
706 % ensure that the ellipse generated is never smaller than the unit
707 % circle and thus never too small for use in EWA resampling.
708 %
709 % This purely mathematical 'magic' was provided by Professor Nicolas
710 % Robidoux and his Masters student Chantal Racette.
711 %
712 % Reference: "We Recommend Singular Value Decomposition", David Austin
713 % http://www.ams.org/samplings/feature-column/fcarc-svd
714 %
715 % By generating major and minor axis vectors, we can actually use the
716 % ellipse in its "canonical form", by remapping the dx,dy of the
717 % sampled point into distances along the major and minor axis unit
718 % vectors.
719 %
720 % Reference: http://en.wikipedia.org/wiki/Ellipse#Canonical_form
721 */
722 static inline void ClampUpAxes(const double dux,
723  const double dvx,
724  const double duy,
725  const double dvy,
726  double *major_mag,
727  double *minor_mag,
728  double *major_unit_x,
729  double *major_unit_y,
730  double *minor_unit_x,
731  double *minor_unit_y)
732 {
733  /*
734  * ClampUpAxes takes an input 2x2 matrix
735  *
736  * [ a b ] = [ dux duy ]
737  * [ c d ] = [ dvx dvy ]
738  *
739  * and computes from it the major and minor axis vectors [major_x,
740  * major_y] and [minor_x,minor_y] of the smallest ellipse containing
741  * both the unit disk and the ellipse which is the image of the unit
742  * disk by the linear transformation
743  *
744  * [ dux duy ] [S] = [s]
745  * [ dvx dvy ] [T] = [t]
746  *
747  * (The vector [S,T] is the difference between a position in output
748  * space and [X,Y]; the vector [s,t] is the difference between a
749  * position in input space and [x,y].)
750  */
751  /*
752  * Output:
753  *
754  * major_mag is the half-length of the major axis of the "new"
755  * ellipse.
756  *
757  * minor_mag is the half-length of the minor axis of the "new"
758  * ellipse.
759  *
760  * major_unit_x is the x-coordinate of the major axis direction vector
761  * of both the "old" and "new" ellipses.
762  *
763  * major_unit_y is the y-coordinate of the major axis direction vector.
764  *
765  * minor_unit_x is the x-coordinate of the minor axis direction vector.
766  *
767  * minor_unit_y is the y-coordinate of the minor axis direction vector.
768  *
769  * Unit vectors are useful for computing projections, in particular,
770  * to compute the distance between a point in output space and the
771  * center of a unit disk in output space, using the position of the
772  * corresponding point [s,t] in input space. Following the clamping,
773  * the square of this distance is
774  *
775  * ( ( s * major_unit_x + t * major_unit_y ) / major_mag )^2
776  * +
777  * ( ( s * minor_unit_x + t * minor_unit_y ) / minor_mag )^2
778  *
779  * If such distances will be computed for many [s,t]'s, it makes
780  * sense to actually compute the reciprocal of major_mag and
781  * minor_mag and multiply them by the above unit lengths.
782  *
783  * Now, if you want to modify the input pair of tangent vectors so
784  * that it defines the modified ellipse, all you have to do is set
785  *
786  * newdux = major_mag * major_unit_x
787  * newdvx = major_mag * major_unit_y
788  * newduy = minor_mag * minor_unit_x = minor_mag * -major_unit_y
789  * newdvy = minor_mag * minor_unit_y = minor_mag * major_unit_x
790  *
791  * and use these tangent vectors as if they were the original ones.
792  * Usually, this is a drastic change in the tangent vectors even if
793  * the singular values are not clamped; for example, the minor axis
794  * vector always points in a direction which is 90 degrees
795  * counterclockwise from the direction of the major axis vector.
796  */
797  /*
798  * Discussion:
799  *
800  * GOAL: Fix things so that the pullback, in input space, of a disk
801  * of radius r in output space is an ellipse which contains, at
802  * least, a disc of radius r. (Make this hold for any r>0.)
803  *
804  * ESSENCE OF THE METHOD: Compute the product of the first two
805  * factors of an SVD of the linear transformation defining the
806  * ellipse and make sure that both its columns have norm at least 1.
807  * Because rotations and reflexions map disks to themselves, it is
808  * not necessary to compute the third (rightmost) factor of the SVD.
809  *
810  * DETAILS: Find the singular values and (unit) left singular
811  * vectors of Jinv, clampling up the singular values to 1, and
812  * multiply the unit left singular vectors by the new singular
813  * values in order to get the minor and major ellipse axis vectors.
814  *
815  * Image resampling context:
816  *
817  * The Jacobian matrix of the transformation at the output point
818  * under consideration is defined as follows:
819  *
820  * Consider the transformation (x,y) -> (X,Y) from input locations
821  * to output locations. (Anthony Thyssen, elsewhere in resample.c,
822  * uses the notation (u,v) -> (x,y).)
823  *
824  * The Jacobian matrix of the transformation at (x,y) is equal to
825  *
826  * J = [ A, B ] = [ dX/dx, dX/dy ]
827  * [ C, D ] [ dY/dx, dY/dy ]
828  *
829  * that is, the vector [A,C] is the tangent vector corresponding to
830  * input changes in the horizontal direction, and the vector [B,D]
831  * is the tangent vector corresponding to input changes in the
832  * vertical direction.
833  *
834  * In the context of resampling, it is natural to use the inverse
835  * Jacobian matrix Jinv because resampling is generally performed by
836  * pulling pixel locations in the output image back to locations in
837  * the input image. Jinv is
838  *
839  * Jinv = [ a, b ] = [ dx/dX, dx/dY ]
840  * [ c, d ] [ dy/dX, dy/dY ]
841  *
842  * Note: Jinv can be computed from J with the following matrix
843  * formula:
844  *
845  * Jinv = 1/(A*D-B*C) [ D, -B ]
846  * [ -C, A ]
847  *
848  * What we do is modify Jinv so that it generates an ellipse which
849  * is as close as possible to the original but which contains the
850  * unit disk. This can be accomplished as follows:
851  *
852  * Let
853  *
854  * Jinv = U Sigma V^T
855  *
856  * be an SVD decomposition of Jinv. (The SVD is not unique, but the
857  * final ellipse does not depend on the particular SVD.)
858  *
859  * We could clamp up the entries of the diagonal matrix Sigma so
860  * that they are at least 1, and then set
861  *
862  * Jinv = U newSigma V^T.
863  *
864  * However, we do not need to compute V for the following reason:
865  * V^T is an orthogonal matrix (that is, it represents a combination
866  * of rotations and reflexions) so that it maps the unit circle to
867  * itself. For this reason, the exact value of V does not affect the
868  * final ellipse, and we can choose V to be the identity
869  * matrix. This gives
870  *
871  * Jinv = U newSigma.
872  *
873  * In the end, we return the two diagonal entries of newSigma
874  * together with the two columns of U.
875  */
876  /*
877  * ClampUpAxes was written by Nicolas Robidoux and Chantal Racette
878  * of Laurentian University with insightful suggestions from Anthony
879  * Thyssen and funding from the National Science and Engineering
880  * Research Council of Canada. It is distinguished from its
881  * predecessors by its efficient handling of degenerate cases.
882  *
883  * The idea of clamping up the EWA ellipse's major and minor axes so
884  * that the result contains the reconstruction kernel filter support
885  * is taken from Andreas Gustafsson's Masters thesis "Interactive
886  * Image Warping", Helsinki University of Technology, Faculty of
887  * Information Technology, 59 pages, 1993 (see Section 3.6).
888  *
889  * The use of the SVD to clamp up the singular values of the
890  * Jacobian matrix of the pullback transformation for EWA resampling
891  * is taken from the astrophysicist Craig DeForest. It is
892  * implemented in his PDL::Transform code (PDL = Perl Data
893  * Language).
894  */
895  const double a = dux;
896  const double b = duy;
897  const double c = dvx;
898  const double d = dvy;
899  /*
900  * n is the matrix Jinv * transpose(Jinv). Eigenvalues of n are the
901  * squares of the singular values of Jinv.
902  */
903  const double aa = a*a;
904  const double bb = b*b;
905  const double cc = c*c;
906  const double dd = d*d;
907  /*
908  * Eigenvectors of n are left singular vectors of Jinv.
909  */
910  const double n11 = aa+bb;
911  const double n12 = a*c+b*d;
912  const double n21 = n12;
913  const double n22 = cc+dd;
914  const double det = a*d-b*c;
915  const double twice_det = det+det;
916  const double frobenius_squared = n11+n22;
917  const double discriminant =
918  (frobenius_squared+twice_det)*(frobenius_squared-twice_det);
919  /*
920  * In exact arithmetic, discriminant can't be negative. In floating
921  * point, it can, because of the bad conditioning of SVD
922  * decompositions done through the associated normal matrix.
923  */
924  const double sqrt_discriminant =
925  sqrt(discriminant > 0.0 ? discriminant : 0.0);
926  /*
927  * s1 is the largest singular value of the inverse Jacobian
928  * matrix. In other words, its reciprocal is the smallest singular
929  * value of the Jacobian matrix itself.
930  * If s1 = 0, both singular values are 0, and any orthogonal pair of
931  * left and right factors produces a singular decomposition of Jinv.
932  */
933  /*
934  * Initially, we only compute the squares of the singular values.
935  */
936  const double s1s1 = 0.5*(frobenius_squared+sqrt_discriminant);
937  /*
938  * s2 the smallest singular value of the inverse Jacobian
939  * matrix. Its reciprocal is the largest singular value of the
940  * Jacobian matrix itself.
941  */
942  const double s2s2 = 0.5*(frobenius_squared-sqrt_discriminant);
943  const double s1s1minusn11 = s1s1-n11;
944  const double s1s1minusn22 = s1s1-n22;
945  /*
946  * u1, the first column of the U factor of a singular decomposition
947  * of Jinv, is a (non-normalized) left singular vector corresponding
948  * to s1. It has entries u11 and u21. We compute u1 from the fact
949  * that it is an eigenvector of n corresponding to the eigenvalue
950  * s1^2.
951  */
952  const double s1s1minusn11_squared = s1s1minusn11*s1s1minusn11;
953  const double s1s1minusn22_squared = s1s1minusn22*s1s1minusn22;
954  /*
955  * The following selects the largest row of n-s1^2 I as the one
956  * which is used to find the eigenvector. If both s1^2-n11 and
957  * s1^2-n22 are zero, n-s1^2 I is the zero matrix. In that case,
958  * any vector is an eigenvector; in addition, norm below is equal to
959  * zero, and, in exact arithmetic, this is the only case in which
960  * norm = 0. So, setting u1 to the simple but arbitrary vector [1,0]
961  * if norm = 0 safely takes care of all cases.
962  */
963  const double temp_u11 =
964  ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? n12 : s1s1minusn22 );
965  const double temp_u21 =
966  ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? s1s1minusn11 : n21 );
967  const double norm = sqrt(temp_u11*temp_u11+temp_u21*temp_u21);
968  /*
969  * Finalize the entries of first left singular vector (associated
970  * with the largest singular value).
971  */
972  const double u11 = ( (norm>0.0) ? temp_u11/norm : 1.0 );
973  const double u21 = ( (norm>0.0) ? temp_u21/norm : 0.0 );
974  /*
975  * Clamp the singular values up to 1.
976  */
977  *major_mag = ( (s1s1<=1.0) ? 1.0 : sqrt(s1s1) );
978  *minor_mag = ( (s2s2<=1.0) ? 1.0 : sqrt(s2s2) );
979  /*
980  * Return the unit major and minor axis direction vectors.
981  */
982  *major_unit_x = u11;
983  *major_unit_y = u21;
984  *minor_unit_x = -u21;
985  *minor_unit_y = u11;
986 }
987 ␌
988 #endif
989 /*
990 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
991 % %
992 % %
993 % %
994 % S c a l e R e s a m p l e F i l t e r %
995 % %
996 % %
997 % %
998 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
999 %
1000 % ScaleResampleFilter() does all the calculations needed to resample an image
1001 % at a specific scale, defined by two scaling vectors. This not using
1002 % a orthogonal scaling, but two distorted scaling vectors, to allow the
1003 % generation of a angled ellipse.
1004 %
1005 % As only two derivative scaling vectors are used the center of the ellipse
1006 % must be the center of the lookup. That is any curvature that the
1007 % distortion may produce is discounted.
1008 %
1009 % The input vectors are produced by either finding the derivatives of the
1010 % distortion function, or the partial derivatives from a distortion mapping.
1011 % They do not need to be the orthogonal dx,dy scaling vectors, but can be
1012 % calculated from other derivatives. For example you could use dr,da/r
1013 % polar coordinate vector scaling vectors
1014 %
1015 % If u,v = DistortEquation(x,y) OR u = Fu(x,y); v = Fv(x,y)
1016 % Then the scaling vectors are determined from the derivatives...
1017 % du/dx, dv/dx and du/dy, dv/dy
1018 % If the resulting scaling vectors is orthogonally aligned then...
1019 % dv/dx = 0 and du/dy = 0
1020 % Producing an orthogonally aligned ellipse in source space for the area to
1021 % be resampled.
1022 %
1023 % Note that scaling vectors are different to argument order. Argument order
1024 % is the general order the derivatives are extracted from the distortion
1025 % equations, and not the scaling vectors. As such the middle two values
1026 % may be swapped from what you expect. Caution is advised.
1027 %
1028 % WARNING: It is assumed that any SetResampleFilter() method call will
1029 % always be performed before the ScaleResampleFilter() method, so that the
1030 % size of the ellipse will match the support for the resampling filter being
1031 % used.
1032 %
1033 % The format of the ScaleResampleFilter method is:
1034 %
1035 % void ScaleResampleFilter(const ResampleFilter *resample_filter,
1036 % const double dux,const double duy,const double dvx,const double dvy)
1037 %
1038 % A description of each parameter follows:
1039 %
1040 % o resample_filter: the resampling information defining the
1041 % image being resampled
1042 %
1043 % o dux,duy,dvx,dvy:
1044 % The derivatives or scaling vectors defining the EWA ellipse.
1045 % NOTE: watch the order, which is based on the order derivatives
1046 % are usually determined from distortion equations (see above).
1047 % The middle two values may need to be swapped if you are thinking
1048 % in terms of scaling vectors.
1049 %
1050 */
1051 MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
1052  const double dux,const double duy,const double dvx,const double dvy)
1053 {
1054  double A,B,C,F;
1055 
1056  assert(resample_filter != (ResampleFilter *) NULL);
1057  assert(resample_filter->signature == MagickCoreSignature);
1058 
1059  resample_filter->limit_reached = MagickFalse;
1060 
1061  /* A 'point' filter forces use of interpolation instead of area sampling */
1062  if ( resample_filter->filter == PointFilter )
1063  return; /* EWA turned off - nothing to do */
1064 
1065 #if DEBUG_ELLIPSE
1066  (void) FormatLocaleFile(stderr, "# -----\n" );
1067  (void) FormatLocaleFile(stderr, "dux=%lf; dvx=%lf; duy=%lf; dvy=%lf;\n",
1068  dux, dvx, duy, dvy);
1069 #endif
1070 
1071  /* Find Ellipse Coefficients such that
1072  A*u^2 + B*u*v + C*v^2 = F
1073  With u,v relative to point around which we are resampling.
1074  And the given scaling dx,dy vectors in u,v space
1075  du/dx,dv/dx and du/dy,dv/dy
1076  */
1077 #if EWA
1078  /* Direct conversion of derivatives into elliptical coefficients
1079  However when magnifying images, the scaling vectors will be small
1080  resulting in a ellipse that is too small to sample properly.
1081  As such we need to clamp the major/minor axis to a minimum of 1.0
1082  to prevent it getting too small.
1083  */
1084 #if EWA_CLAMP
1085  { double major_mag,
1086  minor_mag,
1087  major_x,
1088  major_y,
1089  minor_x,
1090  minor_y;
1091 
1092  ClampUpAxes(dux,dvx,duy,dvy, &major_mag, &minor_mag,
1093  &major_x, &major_y, &minor_x, &minor_y);
1094  major_x *= major_mag; major_y *= major_mag;
1095  minor_x *= minor_mag; minor_y *= minor_mag;
1096 #if DEBUG_ELLIPSE
1097  (void) FormatLocaleFile(stderr, "major_x=%lf; major_y=%lf; minor_x=%lf; minor_y=%lf;\n",
1098  major_x, major_y, minor_x, minor_y);
1099 #endif
1100  A = major_y*major_y+minor_y*minor_y;
1101  B = -2.0*(major_x*major_y+minor_x*minor_y);
1102  C = major_x*major_x+minor_x*minor_x;
1103  F = major_mag*minor_mag;
1104  F *= F; /* square it */
1105  }
1106 #else /* raw unclamped EWA */
1107  A = dvx*dvx+dvy*dvy;
1108  B = -2.0*(dux*dvx+duy*dvy);
1109  C = dux*dux+duy*duy;
1110  F = dux*dvy-duy*dvx;
1111  F *= F; /* square it */
1112 #endif /* EWA_CLAMP */
1113 
1114 #else /* HQ_EWA */
1115  /*
1116  This Paul Heckbert's "Higher Quality EWA" formula, from page 60 in his
1117  thesis, which adds a unit circle to the elliptical area so as to do both
1118  Reconstruction and Prefiltering of the pixels in the resampling. It also
1119  means it is always likely to have at least 4 pixels within the area of the
1120  ellipse, for weighted averaging. No scaling will result with F == 4.0 and
1121  a circle of radius 2.0, and F smaller than this means magnification is
1122  being used.
1123 
1124  NOTE: This method produces a very blurry result at near unity scale while
1125  producing perfect results for strong minification and magnifications.
1126 
1127  However filter support is fixed to 2.0 (no good for Windowed Sinc filters)
1128  */
1129  A = dvx*dvx+dvy*dvy+1;
1130  B = -2.0*(dux*dvx+duy*dvy);
1131  C = dux*dux+duy*duy+1;
1132  F = A*C - B*B/4;
1133 #endif
1134 
1135 #if DEBUG_ELLIPSE
1136  (void) FormatLocaleFile(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
1137 
1138  /* Figure out the various information directly about the ellipse.
1139  This information currently not needed at this time, but may be
1140  needed later for better limit determination.
1141 
1142  It is also good to have as a record for future debugging
1143  */
1144  { double alpha, beta, gamma, Major, Minor;
1145  double Eccentricity, Ellipse_Area, Ellipse_Angle;
1146 
1147  alpha = A+C;
1148  beta = A-C;
1149  gamma = sqrt(beta*beta + B*B );
1150 
1151  if ( alpha - gamma <= MagickEpsilon )
1152  Major= MagickMaximumValue;
1153  else
1154  Major= sqrt(2*F/(alpha - gamma));
1155  Minor = sqrt(2*F/(alpha + gamma));
1156 
1157  (void) FormatLocaleFile(stderr, "# Major=%lf; Minor=%lf\n", Major, Minor );
1158 
1159  /* other information about ellipse include... */
1160  Eccentricity = Major/Minor;
1161  Ellipse_Area = MagickPI*Major*Minor;
1162  Ellipse_Angle = atan2(B, A-C);
1163 
1164  (void) FormatLocaleFile(stderr, "# Angle=%lf Area=%lf\n",
1165  (double) RadiansToDegrees(Ellipse_Angle), Ellipse_Area);
1166  }
1167 #endif
1168 
1169  /* If one or both of the scaling vectors is impossibly large
1170  (producing a very large raw F value), we may as well not bother
1171  doing any form of resampling since resampled area is very large.
1172  In this case some alternative means of pixel sampling, such as
1173  the average of the whole image is needed to get a reasonable
1174  result. Calculate only as needed.
1175  */
1176  if ( (4*A*C - B*B) > MagickMaximumValue ) {
1177  resample_filter->limit_reached = MagickTrue;
1178  return;
1179  }
1180 
1181  /* Scale ellipse to match the filters support
1182  (that is, multiply F by the square of the support)
1183  Simpler to just multiply it by the support twice!
1184  */
1185  F *= resample_filter->support;
1186  F *= resample_filter->support;
1187 
1188  /* Orthogonal bounds of the ellipse */
1189  resample_filter->Ulimit = sqrt(C*F/(A*C-0.25*B*B));
1190  resample_filter->Vlimit = sqrt(A*F/(A*C-0.25*B*B));
1191 
1192  /* Horizontally aligned parallelogram fitted to Ellipse */
1193  resample_filter->Uwidth = sqrt(F/A); /* Half of the parallelogram width */
1194  resample_filter->slope = -B/(2.0*A); /* Reciprocal slope of the parallelogram */
1195 
1196 #if DEBUG_ELLIPSE
1197  (void) FormatLocaleFile(stderr, "Ulimit=%lf; Vlimit=%lf; UWidth=%lf; Slope=%lf;\n",
1198  resample_filter->Ulimit, resample_filter->Vlimit,
1199  resample_filter->Uwidth, resample_filter->slope );
1200 #endif
1201 
1202  /* Check the absolute area of the parallelogram involved.
1203  * This limit needs more work, as it is too slow for larger images
1204  * with tiled views of the horizon.
1205  */
1206  if ( (resample_filter->Uwidth * resample_filter->Vlimit)
1207  > (4.0*resample_filter->image_area)) {
1208  resample_filter->limit_reached = MagickTrue;
1209  return;
1210  }
1211 
1212  /* Scale ellipse formula to directly index the Filter Lookup Table */
1213  { double scale;
1214 #if FILTER_LUT
1215  /* scale so that F = WLUT_WIDTH; -- hardcoded */
1216  scale=(double) WLUT_WIDTH*PerceptibleReciprocal(F);
1217 #else
1218  /* scale so that F = resample_filter->F (support^2) */
1219  scale=resample_filter->F*PerceptibleReciprocal(F);
1220 #endif
1221  resample_filter->A = A*scale;
1222  resample_filter->B = B*scale;
1223  resample_filter->C = C*scale;
1224  }
1225 }
1226 ␌
1227 /*
1228 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1229 % %
1230 % %
1231 % %
1232 % S e t R e s a m p l e F i l t e r %
1233 % %
1234 % %
1235 % %
1236 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1237 %
1238 % SetResampleFilter() set the resampling filter lookup table based on a
1239 % specific filter. Note that the filter is used as a radial filter not as a
1240 % two pass orthogonally aligned resampling filter.
1241 %
1242 % The format of the SetResampleFilter method is:
1243 %
1244 % void SetResampleFilter(ResampleFilter *resample_filter,
1245 % const FilterTypes filter,const double blur)
1246 %
1247 % A description of each parameter follows:
1248 %
1249 % o resample_filter: resampling information structure
1250 %
1251 % o filter: the resize filter for elliptical weighting LUT
1252 %
1253 % o blur: filter blur factor (radial scaling) for elliptical weighting LUT
1254 %
1255 */
1256 MagickExport void SetResampleFilter(ResampleFilter *resample_filter,
1257  const FilterTypes filter,const double blur)
1258 {
1259  ResizeFilter
1260  *resize_filter;
1261 
1262  assert(resample_filter != (ResampleFilter *) NULL);
1263  assert(resample_filter->signature == MagickCoreSignature);
1264 
1265  resample_filter->do_interpolate = MagickFalse;
1266  resample_filter->filter = filter;
1267 
1268  /* Default cylindrical filter is a Cubic Keys filter */
1269  if ( filter == UndefinedFilter )
1270  resample_filter->filter = RobidouxFilter;
1271 
1272  if ( resample_filter->filter == PointFilter ) {
1273  resample_filter->do_interpolate = MagickTrue;
1274  return; /* EWA turned off - nothing more to do */
1275  }
1276 
1277  resize_filter = AcquireResizeFilter(resample_filter->image,
1278  resample_filter->filter,blur,MagickTrue,resample_filter->exception);
1279  if (resize_filter == (ResizeFilter *) NULL) {
1280  (void) ThrowMagickException(resample_filter->exception,GetMagickModule(),
1281  ModuleError, "UnableToSetFilteringValue",
1282  "Fall back to Interpolated 'Point' filter");
1283  resample_filter->filter = PointFilter;
1284  resample_filter->do_interpolate = MagickTrue;
1285  return; /* EWA turned off - nothing more to do */
1286  }
1287 
1288  /* Get the practical working support for the filter,
1289  * after any API call blur factors have been accounted for.
1290  */
1291 #if EWA
1292  resample_filter->support = GetResizeFilterSupport(resize_filter);
1293 #else
1294  resample_filter->support = 2.0; /* fixed support size for HQ-EWA */
1295 #endif
1296 
1297 #if FILTER_LUT
1298  /* Fill the LUT with the weights from the selected filter function */
1299  { int
1300  Q;
1301  double
1302  r_scale;
1303 
1304  /* Scale radius so the filter LUT covers the full support range */
1305  r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
1306  for(Q=0; Q<WLUT_WIDTH; Q++)
1307  resample_filter->filter_lut[Q] = (double)
1308  GetResizeFilterWeight(resize_filter,sqrt((double)Q)*r_scale);
1309 
1310  /* finished with the resize filter */
1311  resize_filter = DestroyResizeFilter(resize_filter);
1312  }
1313 #else
1314  /* save the filter and the scaled ellipse bounds needed for filter */
1315  resample_filter->filter_def = resize_filter;
1316  resample_filter->F = resample_filter->support*resample_filter->support;
1317 #endif
1318 
1319  /*
1320  Adjust the scaling of the default unit circle
1321  This assumes that any real scaling changes will always
1322  take place AFTER the filter method has been initialized.
1323  */
1324  ScaleResampleFilter(resample_filter, 1.0, 0.0, 0.0, 1.0);
1325 
1326 #if 0
1327  /*
1328  This is old code kept as a reference only. Basically it generates
1329  a Gaussian bell curve, with sigma = 0.5 if the support is 2.0
1330 
1331  Create Normal Gaussian 2D Filter Weighted Lookup Table.
1332  A normal EWA guassual lookup would use exp(Q*ALPHA)
1333  where Q = distance squared from 0.0 (center) to 1.0 (edge)
1334  and ALPHA = -4.0*ln(2.0) ==> -2.77258872223978123767
1335  The table is of length 1024, and equates to support radius of 2.0
1336  thus needs to be scaled by ALPHA*4/1024 and any blur factor squared
1337 
1338  The it comes from reference code provided by Fred Weinhaus.
1339  */
1340  r_scale = -2.77258872223978123767/(WLUT_WIDTH*blur*blur);
1341  for(Q=0; Q<WLUT_WIDTH; Q++)
1342  resample_filter->filter_lut[Q] = exp((double)Q*r_scale);
1343  resample_filter->support = WLUT_WIDTH;
1344 #endif
1345 
1346 #if FILTER_LUT
1347 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1348  #pragma omp single
1349 #endif
1350  {
1351  if (IsMagickTrue(GetImageArtifact(resample_filter->image,
1352  "resample:verbose")) )
1353  {
1354  int
1355  Q;
1356  double
1357  r_scale;
1358 
1359  /* Debug output of the filter weighting LUT
1360  Gnuplot the LUT data, the x scale index has been adjusted
1361  plot [0:2][-.2:1] "lut.dat" with lines
1362  The filter values should be normalized for comparison
1363  */
1364  printf("#\n");
1365  printf("# Resampling Filter LUT (%d values) for '%s' filter\n",
1366  WLUT_WIDTH, CommandOptionToMnemonic(MagickFilterOptions,
1367  resample_filter->filter) );
1368  printf("#\n");
1369  printf("# Note: values in table are using a squared radius lookup.\n");
1370  printf("# As such its distribution is not uniform.\n");
1371  printf("#\n");
1372  printf("# The X value is the support distance for the Y weight\n");
1373  printf("# so you can use gnuplot to plot this cylindrical filter\n");
1374  printf("# plot [0:2][-.2:1] \"lut.dat\" with lines\n");
1375  printf("#\n");
1376 
1377  /* Scale radius so the filter LUT covers the full support range */
1378  r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
1379  for(Q=0; Q<WLUT_WIDTH; Q++)
1380  printf("%8.*g %.*g\n",
1381  GetMagickPrecision(),sqrt((double)Q)*r_scale,
1382  GetMagickPrecision(),resample_filter->filter_lut[Q] );
1383  printf("\n\n"); /* generate a 'break' in gnuplot if multiple outputs */
1384  }
1385  /* Output the above once only for each image, and each setting
1386  (void) DeleteImageArtifact(resample_filter->image,"resample:verbose");
1387  */
1388  }
1389 #endif /* FILTER_LUT */
1390  return;
1391 }
1392 ␌
1393 /*
1394 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1395 % %
1396 % %
1397 % %
1398 % S e t R e s a m p l e F i l t e r I n t e r p o l a t e M e t h o d %
1399 % %
1400 % %
1401 % %
1402 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1403 %
1404 % SetResampleFilterInterpolateMethod() sets the resample filter interpolation
1405 % method.
1406 %
1407 % The format of the SetResampleFilterInterpolateMethod method is:
1408 %
1409 % MagickBooleanType SetResampleFilterInterpolateMethod(
1410 % ResampleFilter *resample_filter,const InterpolateMethod method)
1411 %
1412 % A description of each parameter follows:
1413 %
1414 % o resample_filter: the resample filter.
1415 %
1416 % o method: the interpolation method.
1417 %
1418 */
1419 MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
1420  ResampleFilter *resample_filter,const InterpolatePixelMethod method)
1421 {
1422  assert(resample_filter != (ResampleFilter *) NULL);
1423  assert(resample_filter->signature == MagickCoreSignature);
1424  assert(resample_filter->image != (Image *) NULL);
1425  if (IsEventLogging() != MagickFalse)
1426  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1427  resample_filter->image->filename);
1428  resample_filter->interpolate=method;
1429  return(MagickTrue);
1430 }
1431 ␌
1432 /*
1433 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1434 % %
1435 % %
1436 % %
1437 % S e t R e s a m p l e F i l t e r V i r t u a l P i x e l M e t h o d %
1438 % %
1439 % %
1440 % %
1441 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1442 %
1443 % SetResampleFilterVirtualPixelMethod() changes the virtual pixel method
1444 % associated with the specified resample filter.
1445 %
1446 % The format of the SetResampleFilterVirtualPixelMethod method is:
1447 %
1448 % MagickBooleanType SetResampleFilterVirtualPixelMethod(
1449 % ResampleFilter *resample_filter,const VirtualPixelMethod method)
1450 %
1451 % A description of each parameter follows:
1452 %
1453 % o resample_filter: the resample filter.
1454 %
1455 % o method: the virtual pixel method.
1456 %
1457 */
1458 MagickExport MagickBooleanType SetResampleFilterVirtualPixelMethod(
1459  ResampleFilter *resample_filter,const VirtualPixelMethod method)
1460 {
1461  assert(resample_filter != (ResampleFilter *) NULL);
1462  assert(resample_filter->signature == MagickCoreSignature);
1463  assert(resample_filter->image != (Image *) NULL);
1464  if (IsEventLogging() != MagickFalse)
1465  (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1466  resample_filter->image->filename);
1467  resample_filter->virtual_pixel=method;
1468  if (method != UndefinedVirtualPixelMethod)
1469  (void) SetCacheViewVirtualPixelMethod(resample_filter->view,method);
1470  return(MagickTrue);
1471 }
Definition: image.h:134