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