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