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