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