granicus.if.org Git - imagemagick/blob - magick/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-2010 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 "magick/studio.h"
  44 #include "magick/artifact.h"
  45 #include "magick/color-private.h"
  46 #include "magick/cache.h"
  47 #include "magick/draw.h"
  48 #include "magick/exception-private.h"
  49 #include "magick/gem.h"
  50 #include "magick/image.h"
  51 #include "magick/image-private.h"
  52 #include "magick/log.h"
  53 #include "magick/magick.h"
  54 #include "magick/memory_.h"
  55 #include "magick/pixel-private.h"
  56 #include "magick/quantum.h"
  57 #include "magick/random_.h"
  58 #include "magick/resample.h"
  59 #include "magick/resize.h"
  60 #include "magick/resize-private.h"
  61 #include "magick/transform.h"
  62 #include "magick/signature-private.h"
  63 #include "magick/utility.h"
  64 /*
  65   EWA Resampling Options
  66 */
  67
  68 /* select ONE resampling method */
  69 #define EWA 1                 /* Normal EWA handling - raw or clamped */
  70                               /* if 0 then use "High Quality EWA" */
  71 #define EWA_CLAMP 1           /* EWA Clamping from Nicolas Robidoux */
  72
  73 #define FILTER_LUT 1          /* Use a LUT rather then direct filter calls */
  74
  75 /* output debugging information */
  76 #define DEBUG_ELLIPSE 0       /* output ellipse info for debug */
  77 #define DEBUG_HIT_MISS 0      /* output hit/miss pixels (as gnuplot commands) */
  78 #define DEBUG_NO_PIXEL_HIT 0  /* Make pixels that fail to hit anything - RED */
  79
  80 #if ! FILTER_DIRECT
  81 #define WLUT_WIDTH 1024       /* size of the filter cache */
  82 #endif
  83
  84 /*
  85   Typedef declarations.
  86 */
  87 struct _ResampleFilter
  88 {
  89   CacheView
  90     *view;
  91
  92   Image
  93     *image;
  94
  95   ExceptionInfo
  96     *exception;
  97
  98   MagickBooleanType
  99     debug;
 100
 101   /* Information about image being resampled */
 102   ssize_t
 103     image_area;
 104
 105   InterpolatePixelMethod
 106     interpolate;
 107
 108   VirtualPixelMethod
 109     virtual_pixel;
 110
 111   FilterTypes
 112     filter;
 113
 114   /* processing settings needed */
 115   MagickBooleanType
 116     limit_reached,
 117     do_interpolate,
 118     average_defined;
 119
 120   MagickPixelPacket
 121     average_pixel;
 122
 123   /* current ellipitical area being resampled around center point */
 124   double
 125     A, B, C,
 126     Vlimit, Ulimit, Uwidth, slope;
 127
 128 #if FILTER_LUT
 129   /* LUT of weights for filtered average in elliptical area */
 130   double
 131     filter_lut[WLUT_WIDTH];
 132 #else
 133   /* Use a Direct call to the filter functions */
 134   ResizeFilter
 135     *filter_def;
 136
 137   double
 138     F;
 139 #endif
 140
 141   /* the practical working support of the filter */
 142   double
 143     support;
 144
 145   size_t
 146     signature;
 147 };
 148 \f
 149 /*
 150 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 151 %                                                                             %
 152 %                                                                             %
 153 %                                                                             %
 154 %   A c q u i r e R e s a m p l e I n f o                                     %
 155 %                                                                             %
 156 %                                                                             %
 157 %                                                                             %
 158 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 159 %
 160 %  AcquireResampleFilter() initializes the information resample needs do to a
 161 %  scaled lookup of a color from an image, using area sampling.
 162 %
 163 %  The algorithm is based on a Elliptical Weighted Average, where the pixels
 164 %  found in a large elliptical area is averaged together according to a
 165 %  weighting (filter) function.  For more details see "Fundamentals of Texture
 166 %  Mapping and Image Warping" a master's thesis by Paul.S.Heckbert, June 17,
 167 %  1989.  Available for free from, http://www.cs.cmu.edu/~ph/
 168 %
 169 %  As EWA resampling (or any sort of resampling) can require a lot of
 170 %  calculations to produce a distorted scaling of the source image for each
 171 %  output pixel, the ResampleFilter structure generated holds that information
 172 %  between individual image resampling.
 173 %
 174 %  This function will make the appropriate AcquireCacheView() calls
 175 %  to view the image, calling functions do not need to open a cache view.
 176 %
 177 %  Usage Example...
 178 %      resample_filter=AcquireResampleFilter(image,exception);
 179 %      SetResampleFilter(resample_filter, GaussianFilter, 1.0);
 180 %      for (y=0; y < (ssize_t) image->rows; y++) {
 181 %        for (x=0; x < (ssize_t) image->columns; x++) {
 182 %          u= ....;   v= ....;
 183 %          ScaleResampleFilter(resample_filter, ... scaling vectors ...);
 184 %          (void) ResamplePixelColor(resample_filter,u,v,&pixel);
 185 %          ... assign resampled pixel value ...
 186 %        }
 187 %      }
 188 %      DestroyResampleFilter(resample_filter);
 189 %
 190 %  The format of the AcquireResampleFilter method is:
 191 %
 192 %     ResampleFilter *AcquireResampleFilter(const Image *image,
 193 %       ExceptionInfo *exception)
 194 %
 195 %  A description of each parameter follows:
 196 %
 197 %    o image: the image.
 198 %
 199 %    o exception: return any errors or warnings in this structure.
 200 %
 201 */
 202 MagickExport ResampleFilter *AcquireResampleFilter(const Image *image,
 203   ExceptionInfo *exception)
 204 {
 205   register ResampleFilter
 206     *resample_filter;
 207
 208   assert(image != (Image *) NULL);
 209   assert(image->signature == MagickSignature);
 210   if (image->debug != MagickFalse)
 211     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 212   assert(exception != (ExceptionInfo *) NULL);
 213   assert(exception->signature == MagickSignature);
 214
 215   resample_filter=(ResampleFilter *) AcquireMagickMemory(
 216     sizeof(*resample_filter));
 217   if (resample_filter == (ResampleFilter *) NULL)
 218     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 219   (void) ResetMagickMemory(resample_filter,0,sizeof(*resample_filter));
 220
 221   resample_filter->image=ReferenceImage((Image *) image);
 222   resample_filter->view=AcquireCacheView(resample_filter->image);
 223   resample_filter->exception=exception;
 224
 225   resample_filter->debug=IsEventLogging();
 226   resample_filter->signature=MagickSignature;
 227
 228   resample_filter->image_area=(ssize_t) (image->columns*image->rows);
 229   resample_filter->average_defined = MagickFalse;
 230
 231   /* initialise the resampling filter settings */
 232   SetResampleFilter(resample_filter, image->filter, image->blur);
 233   (void) SetResampleFilterInterpolateMethod(resample_filter,
 234     image->interpolate);
 235   (void) SetResampleFilterVirtualPixelMethod(resample_filter,
 236     GetImageVirtualPixelMethod(image));
 237
 238   return(resample_filter);
 239 }
 240 \f
 241 /*
 242 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 243 %                                                                             %
 244 %                                                                             %
 245 %                                                                             %
 246 %   D e s t r o y R e s a m p l e I n f o                                     %
 247 %                                                                             %
 248 %                                                                             %
 249 %                                                                             %
 250 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 251 %
 252 %  DestroyResampleFilter() finalizes and cleans up the resampling
 253 %  resample_filter as returned by AcquireResampleFilter(), freeing any memory
 254 %  or other information as needed.
 255 %
 256 %  The format of the DestroyResampleFilter method is:
 257 %
 258 %      ResampleFilter *DestroyResampleFilter(ResampleFilter *resample_filter)
 259 %
 260 %  A description of each parameter follows:
 261 %
 262 %    o resample_filter: resampling information structure
 263 %
 264 */
 265 MagickExport ResampleFilter *DestroyResampleFilter(
 266   ResampleFilter *resample_filter)
 267 {
 268   assert(resample_filter != (ResampleFilter *) NULL);
 269   assert(resample_filter->signature == MagickSignature);
 270   assert(resample_filter->image != (Image *) NULL);
 271   if (resample_filter->debug != MagickFalse)
 272     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
 273       resample_filter->image->filename);
 274   resample_filter->view=DestroyCacheView(resample_filter->view);
 275   resample_filter->image=DestroyImage(resample_filter->image);
 276 #if ! FILTER_LUT
 277   resample_filter->filter_def=DestroyResizeFilter(resample_filter->filter_def);
 278 #endif
 279   resample_filter->signature=(~MagickSignature);
 280   resample_filter=(ResampleFilter *) RelinquishMagickMemory(resample_filter);
 281   return(resample_filter);
 282 }
 283 \f
 284 /*
 285 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 286 %                                                                             %
 287 %                                                                             %
 288 %                                                                             %
 289 %   I n t e r p o l a t e R e s a m p l e F i l t e r                         %
 290 %                                                                             %
 291 %                                                                             %
 292 %                                                                             %
 293 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 294 %
 295 %  InterpolateResampleFilter() applies bi-linear or tri-linear interpolation
 296 %  between a floating point coordinate and the pixels surrounding that
 297 %  coordinate.  No pixel area resampling, or scaling of the result is
 298 %  performed.
 299 %
 300 %  The format of the InterpolateResampleFilter method is:
 301 %
 302 %      MagickBooleanType InterpolateResampleFilter(
 303 %        ResampleInfo *resample_filter,const InterpolatePixelMethod method,
 304 %        const double x,const double y,MagickPixelPacket *pixel)
 305 %
 306 %  A description of each parameter follows:
 307 %
 308 %    o resample_filter: the resample filter.
 309 %
 310 %    o method: the pixel clor interpolation method.
 311 %
 312 %    o x,y: A double representing the current (x,y) position of the pixel.
 313 %
 314 %    o pixel: return the interpolated pixel here.
 315 %
 316 */
 317
 318 static inline double MagickMax(const double x,const double y)
 319 {
 320   if (x > y)
 321     return(x);
 322   return(y);
 323 }
 324
 325 static void BicubicInterpolate(const MagickPixelPacket *pixels,const double dx,
 326   MagickPixelPacket *pixel)
 327 {
 328   MagickRealType
 329     dx2,
 330     p,
 331     q,
 332     r,
 333     s;
 334
 335   dx2=dx*dx;
 336   p=(pixels[3].red-pixels[2].red)-(pixels[0].red-pixels[1].red);
 337   q=(pixels[0].red-pixels[1].red)-p;
 338   r=pixels[2].red-pixels[0].red;
 339   s=pixels[1].red;
 340   pixel->red=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
 341   p=(pixels[3].green-pixels[2].green)-(pixels[0].green-pixels[1].green);
 342   q=(pixels[0].green-pixels[1].green)-p;
 343   r=pixels[2].green-pixels[0].green;
 344   s=pixels[1].green;
 345   pixel->green=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
 346   p=(pixels[3].blue-pixels[2].blue)-(pixels[0].blue-pixels[1].blue);
 347   q=(pixels[0].blue-pixels[1].blue)-p;
 348   r=pixels[2].blue-pixels[0].blue;
 349   s=pixels[1].blue;
 350   pixel->blue=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
 351   p=(pixels[3].opacity-pixels[2].opacity)-(pixels[0].opacity-pixels[1].opacity);
 352   q=(pixels[0].opacity-pixels[1].opacity)-p;
 353   r=pixels[2].opacity-pixels[0].opacity;
 354   s=pixels[1].opacity;
 355   pixel->opacity=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
 356   if (pixel->colorspace == CMYKColorspace)
 357     {
 358       p=(pixels[3].index-pixels[2].index)-(pixels[0].index-pixels[1].index);
 359       q=(pixels[0].index-pixels[1].index)-p;
 360       r=pixels[2].index-pixels[0].index;
 361       s=pixels[1].index;
 362       pixel->index=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
 363     }
 364 }
 365
 366 static inline MagickRealType CubicWeightingFunction(const MagickRealType x)
 367 {
 368   MagickRealType
 369     alpha,
 370     gamma;
 371
 372   alpha=MagickMax(x+2.0,0.0);
 373   gamma=1.0*alpha*alpha*alpha;
 374   alpha=MagickMax(x+1.0,0.0);
 375   gamma-=4.0*alpha*alpha*alpha;
 376   alpha=MagickMax(x+0.0,0.0);
 377   gamma+=6.0*alpha*alpha*alpha;
 378   alpha=MagickMax(x-1.0,0.0);
 379   gamma-=4.0*alpha*alpha*alpha;
 380   return(gamma/6.0);
 381 }
 382
 383 static inline double MeshInterpolate(const PointInfo *delta,const double p,
 384   const double x,const double y)
 385 {
 386   return(delta->x*x+delta->y*y+(1.0-delta->x-delta->y)*p);
 387 }
 388
 389 static inline ssize_t NearestNeighbor(MagickRealType x)
 390 {
 391   if (x >= 0.0)
 392     return((ssize_t) (x+0.5));
 393   return((ssize_t) (x-0.5));
 394 }
 395
 396 static MagickBooleanType InterpolateResampleFilter(
 397   ResampleFilter *resample_filter,const InterpolatePixelMethod method,
 398   const double x,const double y,MagickPixelPacket *pixel)
 399 {
 400   MagickBooleanType
 401     status;
 402
 403   register const IndexPacket
 404     *indexes;
 405
 406   register const PixelPacket
 407     *p;
 408
 409   register ssize_t
 410     i;
 411
 412   assert(resample_filter != (ResampleFilter *) NULL);
 413   assert(resample_filter->signature == MagickSignature);
 414   status=MagickTrue;
 415
 416   switch (method)
 417   {
 418     case AverageInterpolatePixel:
 419     {
 420       MagickPixelPacket
 421         pixels[16];
 422
 423       MagickRealType
 424         alpha[16],
 425         gamma;
 426
 427       p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x)-1,
 428         (ssize_t) floor(y)-1,4,4,resample_filter->exception);
 429       if (p == (const PixelPacket *) NULL)
 430         {
 431           status=MagickFalse;
 432           break;
 433         }
 434       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 435       for (i=0; i < 16L; i++)
 436       {
 437         GetMagickPixelPacket(resample_filter->image,pixels+i);
 438         SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
 439         alpha[i]=1.0;
 440         if (resample_filter->image->matte != MagickFalse)
 441           {
 442             alpha[i]=QuantumScale*((MagickRealType) GetAlphaPixelComponent(p));
 443             pixels[i].red*=alpha[i];
 444             pixels[i].green*=alpha[i];
 445             pixels[i].blue*=alpha[i];
 446             if (resample_filter->image->colorspace == CMYKColorspace)
 447               pixels[i].index*=alpha[i];
 448           }
 449         gamma=alpha[i];
 450         gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 451         pixel->red+=gamma*0.0625*pixels[i].red;
 452         pixel->green+=gamma*0.0625*pixels[i].green;
 453         pixel->blue+=gamma*0.0625*pixels[i].blue;
 454         pixel->opacity+=0.0625*pixels[i].opacity;
 455         if (resample_filter->image->colorspace == CMYKColorspace)
 456           pixel->index+=gamma*0.0625*pixels[i].index;
 457         p++;
 458       }
 459       break;
 460     }
 461     case BicubicInterpolatePixel:
 462     {
 463       MagickPixelPacket
 464         pixels[16],
 465         u[4];
 466
 467       MagickRealType
 468         alpha[16];
 469
 470       PointInfo
 471         delta;
 472
 473       p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x)-1,
 474         (ssize_t) floor(y)-1,4,4,resample_filter->exception);
 475       if (p == (const PixelPacket *) NULL)
 476         {
 477           status=MagickFalse;
 478           break;
 479         }
 480       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 481       for (i=0; i < 16L; i++)
 482       {
 483         GetMagickPixelPacket(resample_filter->image,pixels+i);
 484         SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
 485         alpha[i]=1.0;
 486         if (resample_filter->image->matte != MagickFalse)
 487           {
 488             alpha[i]=QuantumScale*((MagickRealType) GetAlphaPixelComponent(p));
 489             pixels[i].red*=alpha[i];
 490             pixels[i].green*=alpha[i];
 491             pixels[i].blue*=alpha[i];
 492             if (resample_filter->image->colorspace == CMYKColorspace)
 493               pixels[i].index*=alpha[i];
 494           }
 495         p++;
 496       }
 497       delta.x=x-floor(x);
 498       for (i=0; i < 4L; i++)
 499         BicubicInterpolate(pixels+4*i,delta.x,u+i);
 500       delta.y=y-floor(y);
 501       BicubicInterpolate(u,delta.y,pixel);
 502       break;
 503     }
 504     case BilinearInterpolatePixel:
 505     default:
 506     {
 507       MagickPixelPacket
 508         pixels[4];
 509
 510       MagickRealType
 511         alpha[4],
 512         gamma;
 513
 514       PointInfo
 515         delta,
 516         epsilon;
 517
 518       p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x),
 519         (ssize_t) floor(y),2,2,resample_filter->exception);
 520       if (p == (const PixelPacket *) NULL)
 521         {
 522           status=MagickFalse;
 523           break;
 524         }
 525       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 526       for (i=0; i < 4L; i++)
 527       {
 528         pixels[i].red=(MagickRealType) p[i].red;
 529         pixels[i].green=(MagickRealType) p[i].green;
 530         pixels[i].blue=(MagickRealType) p[i].blue;
 531         pixels[i].opacity=(MagickRealType) p[i].opacity;
 532         alpha[i]=1.0;
 533       }
 534       if (resample_filter->image->matte != MagickFalse)
 535         for (i=0; i < 4L; i++)
 536         {
 537           alpha[i]=QuantumScale*((MagickRealType) QuantumRange-p[i].opacity);
 538           pixels[i].red*=alpha[i];
 539           pixels[i].green*=alpha[i];
 540           pixels[i].blue*=alpha[i];
 541         }
 542       if (indexes != (IndexPacket *) NULL)
 543         for (i=0; i < 4L; i++)
 544         {
 545           pixels[i].index=(MagickRealType) indexes[i];
 546           if (resample_filter->image->colorspace == CMYKColorspace)
 547             pixels[i].index*=alpha[i];
 548         }
 549       delta.x=x-floor(x);
 550       delta.y=y-floor(y);
 551       epsilon.x=1.0-delta.x;
 552       epsilon.y=1.0-delta.y;
 553       gamma=((epsilon.y*(epsilon.x*alpha[0]+delta.x*alpha[1])+delta.y*
 554         (epsilon.x*alpha[2]+delta.x*alpha[3])));
 555       gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 556       pixel->red=gamma*(epsilon.y*(epsilon.x*pixels[0].red+delta.x*
 557         pixels[1].red)+delta.y*(epsilon.x*pixels[2].red+delta.x*pixels[3].red));
 558       pixel->green=gamma*(epsilon.y*(epsilon.x*pixels[0].green+delta.x*
 559         pixels[1].green)+delta.y*(epsilon.x*pixels[2].green+delta.x*
 560         pixels[3].green));
 561       pixel->blue=gamma*(epsilon.y*(epsilon.x*pixels[0].blue+delta.x*
 562         pixels[1].blue)+delta.y*(epsilon.x*pixels[2].blue+delta.x*
 563         pixels[3].blue));
 564       pixel->opacity=(epsilon.y*(epsilon.x*pixels[0].opacity+delta.x*
 565         pixels[1].opacity)+delta.y*(epsilon.x*pixels[2].opacity+delta.x*
 566         pixels[3].opacity));
 567       if (resample_filter->image->colorspace == CMYKColorspace)
 568         pixel->index=gamma*(epsilon.y*(epsilon.x*pixels[0].index+delta.x*
 569           pixels[1].index)+delta.y*(epsilon.x*pixels[2].index+delta.x*
 570           pixels[3].index));
 571       break;
 572     }
 573     case FilterInterpolatePixel:
 574     {
 575       CacheView
 576         *filter_view;
 577
 578       Image
 579         *excerpt_image,
 580         *filter_image;
 581
 582       MagickPixelPacket
 583         pixels[1];
 584
 585       RectangleInfo
 586         geometry;
 587
 588       geometry.width=4L;
 589       geometry.height=4L;
 590       geometry.x=(ssize_t) floor(x)-1L;
 591       geometry.y=(ssize_t) floor(y)-1L;
 592       excerpt_image=ExcerptImage(resample_filter->image,&geometry,
 593         resample_filter->exception);
 594       if (excerpt_image == (Image *) NULL)
 595         {
 596           status=MagickFalse;
 597           break;
 598         }
 599       filter_image=ResizeImage(excerpt_image,1,1,resample_filter->image->filter,
 600         resample_filter->image->blur,resample_filter->exception);
 601       excerpt_image=DestroyImage(excerpt_image);
 602       if (filter_image == (Image *) NULL)
 603         break;
 604       filter_view=AcquireCacheView(filter_image);
 605       p=GetCacheViewVirtualPixels(filter_view,0,0,1,1,
 606         resample_filter->exception);
 607       if (p != (const PixelPacket *) NULL)
 608         {
 609           indexes=GetVirtualIndexQueue(filter_image);
 610           GetMagickPixelPacket(resample_filter->image,pixels);
 611           SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
 612         }
 613       filter_view=DestroyCacheView(filter_view);
 614       filter_image=DestroyImage(filter_image);
 615       break;
 616     }
 617     case IntegerInterpolatePixel:
 618     {
 619       MagickPixelPacket
 620         pixels[1];
 621
 622       p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x),
 623         (ssize_t) floor(y),1,1,resample_filter->exception);
 624       if (p == (const PixelPacket *) NULL)
 625         {
 626           status=MagickFalse;
 627           break;
 628         }
 629       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 630       GetMagickPixelPacket(resample_filter->image,pixels);
 631       SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
 632       break;
 633     }
 634     case MeshInterpolatePixel:
 635     {
 636       MagickPixelPacket
 637         pixels[4];
 638
 639       MagickRealType
 640         alpha[4],
 641         gamma;
 642
 643       PointInfo
 644         delta,
 645         luminance;
 646
 647       p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x),
 648         (ssize_t) floor(y),2,2,resample_filter->exception);
 649       if (p == (const PixelPacket *) NULL)
 650         {
 651           status=MagickFalse;
 652           break;
 653         }
 654       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 655       for (i=0; i < 4L; i++)
 656       {
 657         GetMagickPixelPacket(resample_filter->image,pixels+i);
 658         SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
 659         alpha[i]=1.0;
 660         if (resample_filter->image->matte != MagickFalse)
 661           {
 662             alpha[i]=QuantumScale*((MagickRealType) GetAlphaPixelComponent(p));
 663             pixels[i].red*=alpha[i];
 664             pixels[i].green*=alpha[i];
 665             pixels[i].blue*=alpha[i];
 666             if (resample_filter->image->colorspace == CMYKColorspace)
 667               pixels[i].index*=alpha[i];
 668           }
 669         p++;
 670       }
 671       delta.x=x-floor(x);
 672       delta.y=y-floor(y);
 673       luminance.x=MagickPixelLuminance(pixels+0)-MagickPixelLuminance(pixels+3);
 674       luminance.y=MagickPixelLuminance(pixels+1)-MagickPixelLuminance(pixels+2);
 675       if (fabs(luminance.x) < fabs(luminance.y))
 676         {
 677           /*
 678             Diagonal 0-3 NW-SE.
 679           */
 680           if (delta.x <= delta.y)
 681             {
 682               /*
 683                 Bottom-left triangle  (pixel:2, diagonal: 0-3).
 684               */
 685               delta.y=1.0-delta.y;
 686               gamma=MeshInterpolate(&delta,alpha[2],alpha[3],alpha[0]);
 687               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 688               pixel->red=gamma*MeshInterpolate(&delta,pixels[2].red,
 689                 pixels[3].red,pixels[0].red);
 690               pixel->green=gamma*MeshInterpolate(&delta,pixels[2].green,
 691                 pixels[3].green,pixels[0].green);
 692               pixel->blue=gamma*MeshInterpolate(&delta,pixels[2].blue,
 693                 pixels[3].blue,pixels[0].blue);
 694               pixel->opacity=gamma*MeshInterpolate(&delta,pixels[2].opacity,
 695                 pixels[3].opacity,pixels[0].opacity);
 696               if (resample_filter->image->colorspace == CMYKColorspace)
 697                 pixel->index=gamma*MeshInterpolate(&delta,pixels[2].index,
 698                   pixels[3].index,pixels[0].index);
 699             }
 700           else
 701             {
 702               /*
 703                 Top-right triangle (pixel:1, diagonal: 0-3).
 704               */
 705               delta.x=1.0-delta.x;
 706               gamma=MeshInterpolate(&delta,alpha[1],alpha[0],alpha[3]);
 707               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 708               pixel->red=gamma*MeshInterpolate(&delta,pixels[1].red,
 709                 pixels[0].red,pixels[3].red);
 710               pixel->green=gamma*MeshInterpolate(&delta,pixels[1].green,
 711                 pixels[0].green,pixels[3].green);
 712               pixel->blue=gamma*MeshInterpolate(&delta,pixels[1].blue,
 713                 pixels[0].blue,pixels[3].blue);
 714               pixel->opacity=gamma*MeshInterpolate(&delta,pixels[1].opacity,
 715                 pixels[0].opacity,pixels[3].opacity);
 716               if (resample_filter->image->colorspace == CMYKColorspace)
 717                 pixel->index=gamma*MeshInterpolate(&delta,pixels[1].index,
 718                   pixels[0].index,pixels[3].index);
 719             }
 720         }
 721       else
 722         {
 723           /*
 724             Diagonal 1-2 NE-SW.
 725           */
 726           if (delta.x <= (1.0-delta.y))
 727             {
 728               /*
 729                 Top-left triangle (pixel 0, diagonal: 1-2).
 730               */
 731               gamma=MeshInterpolate(&delta,alpha[0],alpha[1],alpha[2]);
 732               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 733               pixel->red=gamma*MeshInterpolate(&delta,pixels[0].red,
 734                 pixels[1].red,pixels[2].red);
 735               pixel->green=gamma*MeshInterpolate(&delta,pixels[0].green,
 736                 pixels[1].green,pixels[2].green);
 737               pixel->blue=gamma*MeshInterpolate(&delta,pixels[0].blue,
 738                 pixels[1].blue,pixels[2].blue);
 739               pixel->opacity=gamma*MeshInterpolate(&delta,pixels[0].opacity,
 740                 pixels[1].opacity,pixels[2].opacity);
 741               if (resample_filter->image->colorspace == CMYKColorspace)
 742                 pixel->index=gamma*MeshInterpolate(&delta,pixels[0].index,
 743                   pixels[1].index,pixels[2].index);
 744             }
 745           else
 746             {
 747               /*
 748                 Bottom-right triangle (pixel: 3, diagonal: 1-2).
 749               */
 750               delta.x=1.0-delta.x;
 751               delta.y=1.0-delta.y;
 752               gamma=MeshInterpolate(&delta,alpha[3],alpha[2],alpha[1]);
 753               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 754               pixel->red=gamma*MeshInterpolate(&delta,pixels[3].red,
 755                 pixels[2].red,pixels[1].red);
 756               pixel->green=gamma*MeshInterpolate(&delta,pixels[3].green,
 757                 pixels[2].green,pixels[1].green);
 758               pixel->blue=gamma*MeshInterpolate(&delta,pixels[3].blue,
 759                 pixels[2].blue,pixels[1].blue);
 760               pixel->opacity=gamma*MeshInterpolate(&delta,pixels[3].opacity,
 761                 pixels[2].opacity,pixels[1].opacity);
 762               if (resample_filter->image->colorspace == CMYKColorspace)
 763                 pixel->index=gamma*MeshInterpolate(&delta,pixels[3].index,
 764                   pixels[2].index,pixels[1].index);
 765             }
 766         }
 767       break;
 768     }
 769     case NearestNeighborInterpolatePixel:
 770     {
 771       MagickPixelPacket
 772         pixels[1];
 773
 774       p=GetCacheViewVirtualPixels(resample_filter->view,NearestNeighbor(x),
 775         NearestNeighbor(y),1,1,resample_filter->exception);
 776       if (p == (const PixelPacket *) NULL)
 777         {
 778           status=MagickFalse;
 779           break;
 780         }
 781       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 782       GetMagickPixelPacket(resample_filter->image,pixels);
 783       SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
 784       break;
 785     }
 786     case SplineInterpolatePixel:
 787     {
 788       MagickPixelPacket
 789         pixels[16];
 790
 791       MagickRealType
 792         alpha[16],
 793         dx,
 794         dy,
 795         gamma;
 796
 797       PointInfo
 798         delta;
 799
 800       ssize_t
 801         j,
 802         n;
 803
 804       p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x)-1,
 805         (ssize_t) floor(y)-1,4,4,resample_filter->exception);
 806       if (p == (const PixelPacket *) NULL)
 807         {
 808           status=MagickFalse;
 809           break;
 810         }
 811       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 812       n=0;
 813       delta.x=x-floor(x);
 814       delta.y=y-floor(y);
 815       for (i=(-1); i < 3L; i++)
 816       {
 817         dy=CubicWeightingFunction((MagickRealType) i-delta.y);
 818         for (j=(-1); j < 3L; j++)
 819         {
 820           GetMagickPixelPacket(resample_filter->image,pixels+n);
 821           SetMagickPixelPacket(resample_filter->image,p,indexes+n,pixels+n);
 822           alpha[n]=1.0;
 823           if (resample_filter->image->matte != MagickFalse)
 824             {
 825               alpha[n]=QuantumScale*((MagickRealType)
 826                 GetAlphaPixelComponent(p));
 827               pixels[n].red*=alpha[n];
 828               pixels[n].green*=alpha[n];
 829               pixels[n].blue*=alpha[n];
 830               if (resample_filter->image->colorspace == CMYKColorspace)
 831                 pixels[n].index*=alpha[n];
 832             }
 833           dx=CubicWeightingFunction(delta.x-(MagickRealType) j);
 834           gamma=alpha[n];
 835           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 836           pixel->red+=gamma*dx*dy*pixels[n].red;
 837           pixel->green+=gamma*dx*dy*pixels[n].green;
 838           pixel->blue+=gamma*dx*dy*pixels[n].blue;
 839           if (resample_filter->image->matte != MagickFalse)
 840             pixel->opacity+=dx*dy*pixels[n].opacity;
 841           if (resample_filter->image->colorspace == CMYKColorspace)
 842             pixel->index+=gamma*dx*dy*pixels[n].index;
 843           n++;
 844           p++;
 845         }
 846       }
 847       break;
 848     }
 849   }
 850   return(status);
 851 }
 852 \f
 853 /*
 854 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 855 %                                                                             %
 856 %                                                                             %
 857 %                                                                             %
 858 %   R e s a m p l e P i x e l C o l o r                                       %
 859 %                                                                             %
 860 %                                                                             %
 861 %                                                                             %
 862 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 863 %
 864 %  ResamplePixelColor() samples the pixel values surrounding the location
 865 %  given using an elliptical weighted average, at the scale previously
 866 %  calculated, and in the most efficent manner possible for the
 867 %  VirtualPixelMethod setting.
 868 %
 869 %  The format of the ResamplePixelColor method is:
 870 %
 871 %     MagickBooleanType ResamplePixelColor(ResampleFilter *resample_filter,
 872 %       const double u0,const double v0,MagickPixelPacket *pixel)
 873 %
 874 %  A description of each parameter follows:
 875 %
 876 %    o resample_filter: the resample filter.
 877 %
 878 %    o u0,v0: A double representing the center of the area to resample,
 879 %        The distortion transformed transformed x,y coordinate.
 880 %
 881 %    o pixel: the resampled pixel is returned here.
 882 %
 883 */
 884 MagickExport MagickBooleanType ResamplePixelColor(
 885   ResampleFilter *resample_filter,const double u0,const double v0,
 886   MagickPixelPacket *pixel)
 887 {
 888   MagickBooleanType
 889     status;
 890
 891   ssize_t u,v, v1, v2, uw, hit;
 892   double u1;
 893   double U,V,Q,DQ,DDQ;
 894   double divisor_c,divisor_m;
 895   register double weight;
 896   register const PixelPacket *pixels;
 897   register const IndexPacket *indexes;
 898   assert(resample_filter != (ResampleFilter *) NULL);
 899   assert(resample_filter->signature == MagickSignature);
 900
 901   status=MagickTrue;
 902   GetMagickPixelPacket(resample_filter->image,pixel);
 903   if ( resample_filter->do_interpolate ) {
 904     status=InterpolateResampleFilter(resample_filter,
 905       resample_filter->interpolate,u0,v0,pixel);
 906     return(status);
 907   }
 908
 909 #if DEBUG_ELLIPSE
 910   fprintf(stderr, "u0=%lf; v0=%lf;\n", u0, v0);
 911 #endif
 912
 913   /*
 914     Does resample area Miss the image?
 915     And is that area a simple solid color - then return that color
 916   */
 917   hit = 0;
 918   switch ( resample_filter->virtual_pixel ) {
 919     case BackgroundVirtualPixelMethod:
 920     case ConstantVirtualPixelMethod:
 921     case TransparentVirtualPixelMethod:
 922     case BlackVirtualPixelMethod:
 923     case GrayVirtualPixelMethod:
 924     case WhiteVirtualPixelMethod:
 925     case MaskVirtualPixelMethod:
 926       if ( resample_filter->limit_reached
 927            || u0 + resample_filter->Ulimit < 0.0
 928            || u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 929            || v0 + resample_filter->Vlimit < 0.0
 930            || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows
 931            )
 932         hit++;
 933       break;
 934
 935     case UndefinedVirtualPixelMethod:
 936     case EdgeVirtualPixelMethod:
 937       if (    ( u0 + resample_filter->Ulimit < 0.0 && v0 + resample_filter->Vlimit < 0.0 )
 938            || ( u0 + resample_filter->Ulimit < 0.0
 939                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows )
 940            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 941                 && v0 + resample_filter->Vlimit < 0.0 )
 942            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 943                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows )
 944            )
 945         hit++;
 946       break;
 947     case HorizontalTileVirtualPixelMethod:
 948       if (    v0 + resample_filter->Vlimit < 0.0
 949            || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows
 950            )
 951         hit++;  /* outside the horizontally tiled images. */
 952       break;
 953     case VerticalTileVirtualPixelMethod:
 954       if (    u0 + resample_filter->Ulimit < 0.0
 955            || u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 956            )
 957         hit++;  /* outside the vertically tiled images. */
 958       break;
 959     case DitherVirtualPixelMethod:
 960       if (    ( u0 + resample_filter->Ulimit < -32.0 && v0 + resample_filter->Vlimit < -32.0 )
 961            || ( u0 + resample_filter->Ulimit < -32.0
 962                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+32.0 )
 963            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+32.0
 964                 && v0 + resample_filter->Vlimit < -32.0 )
 965            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+32.0
 966                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+32.0 )
 967            )
 968         hit++;
 969       break;
 970     case TileVirtualPixelMethod:
 971     case MirrorVirtualPixelMethod:
 972     case RandomVirtualPixelMethod:
 973     case HorizontalTileEdgeVirtualPixelMethod:
 974     case VerticalTileEdgeVirtualPixelMethod:
 975     case CheckerTileVirtualPixelMethod:
 976       /* resampling of area is always needed - no VP limits */
 977       break;
 978   }
 979   if ( hit ) {
 980     /* whole area is a solid color -- just return that color */
 981     status=InterpolateResampleFilter(resample_filter,IntegerInterpolatePixel,
 982       u0,v0,pixel);
 983     return(status);
 984   }
 985
 986   /*
 987     Scaling limits reached, return an 'averaged' result.
 988   */
 989   if ( resample_filter->limit_reached ) {
 990     switch ( resample_filter->virtual_pixel ) {
 991       /*  This is always handled by the above, so no need.
 992         case BackgroundVirtualPixelMethod:
 993         case ConstantVirtualPixelMethod:
 994         case TransparentVirtualPixelMethod:
 995         case GrayVirtualPixelMethod,
 996         case WhiteVirtualPixelMethod
 997         case MaskVirtualPixelMethod:
 998       */
 999       case UndefinedVirtualPixelMethod:
1000       case EdgeVirtualPixelMethod:
1001       case DitherVirtualPixelMethod:
1002       case HorizontalTileEdgeVirtualPixelMethod:
1003       case VerticalTileEdgeVirtualPixelMethod:
1004         /* We need an average edge pixel, from the correct edge!
1005            How should I calculate an average edge color?
1006            Just returning an averaged neighbourhood,
1007            works well in general, but falls down for TileEdge methods.
1008            This needs to be done properly!!!!!!
1009         */
1010         status=InterpolateResampleFilter(resample_filter,
1011           AverageInterpolatePixel,u0,v0,pixel);
1012         break;
1013       case HorizontalTileVirtualPixelMethod:
1014       case VerticalTileVirtualPixelMethod:
1015         /* just return the background pixel - Is there more direct way? */
1016         status=InterpolateResampleFilter(resample_filter,
1017            IntegerInterpolatePixel,(double)-1,(double)-1,pixel);
1018         break;
1019       case TileVirtualPixelMethod:
1020       case MirrorVirtualPixelMethod:
1021       case RandomVirtualPixelMethod:
1022       case CheckerTileVirtualPixelMethod:
1023       default:
1024         /* generate a average color of the WHOLE image */
1025         if ( resample_filter->average_defined == MagickFalse ) {
1026           Image
1027             *average_image;
1028
1029           CacheView
1030             *average_view;
1031
1032           GetMagickPixelPacket(resample_filter->image,(MagickPixelPacket *)
1033             &resample_filter->average_pixel);
1034           resample_filter->average_defined=MagickTrue;
1035
1036           /* Try to get an averaged pixel color of whole image */
1037           average_image=ResizeImage(resample_filter->image,1,1,BoxFilter,1.0,
1038             resample_filter->exception);
1039           if (average_image == (Image *) NULL)
1040             {
1041               *pixel=resample_filter->average_pixel; /* FAILED */
1042               break;
1043             }
1044           average_view=AcquireCacheView(average_image);
1045           pixels=(PixelPacket *)GetCacheViewVirtualPixels(average_view,0,0,1,1,
1046             resample_filter->exception);
1047           if (pixels == (const PixelPacket *) NULL) {
1048             average_view=DestroyCacheView(average_view);
1049             average_image=DestroyImage(average_image);
1050             *pixel=resample_filter->average_pixel; /* FAILED */
1051             break;
1052           }
1053           indexes=(IndexPacket *) GetCacheViewAuthenticIndexQueue(average_view);
1054           SetMagickPixelPacket(resample_filter->image,pixels,indexes,
1055             &(resample_filter->average_pixel));
1056           average_view=DestroyCacheView(average_view);
1057           average_image=DestroyImage(average_image);
1058
1059           if ( resample_filter->virtual_pixel == CheckerTileVirtualPixelMethod )
1060             {
1061               /* CheckerTile is avergae of image average half background */
1062               /* FUTURE: replace with a 50% blend of both pixels */
1063
1064               weight = QuantumScale*((MagickRealType)(QuantumRange-
1065                           resample_filter->average_pixel.opacity));
1066               resample_filter->average_pixel.red *= weight;
1067               resample_filter->average_pixel.green *= weight;
1068               resample_filter->average_pixel.blue *= weight;
1069               divisor_c = weight;
1070
1071               weight = QuantumScale*((MagickRealType)(QuantumRange-
1072                           resample_filter->image->background_color.opacity));
1073               resample_filter->average_pixel.red +=
1074                       weight*resample_filter->image->background_color.red;
1075               resample_filter->average_pixel.green +=
1076                       weight*resample_filter->image->background_color.green;
1077               resample_filter->average_pixel.blue +=
1078                       weight*resample_filter->image->background_color.blue;
1079               resample_filter->average_pixel.opacity +=
1080                       resample_filter->image->background_color.opacity;
1081               divisor_c += weight;
1082
1083               resample_filter->average_pixel.red /= divisor_c;
1084               resample_filter->average_pixel.green /= divisor_c;
1085               resample_filter->average_pixel.blue /= divisor_c;
1086               resample_filter->average_pixel.opacity /= 2;
1087
1088             }
1089         }
1090         *pixel=resample_filter->average_pixel;
1091         break;
1092     }
1093     return(status);
1094   }
1095
1096   /*
1097     Initialize weighted average data collection
1098   */
1099   hit = 0;
1100   divisor_c = 0.0;
1101   divisor_m = 0.0;
1102   pixel->red = pixel->green = pixel->blue = 0.0;
1103   if (resample_filter->image->matte != MagickFalse) pixel->opacity = 0.0;
1104   if (resample_filter->image->colorspace == CMYKColorspace) pixel->index = 0.0;
1105
1106   /*
1107     Determine the parellelogram bounding box fitted to the ellipse
1108     centered at u0,v0.  This area is bounding by the lines...
1109   */
1110   v1 = (ssize_t)ceil(v0 - resample_filter->Vlimit);  /* range of scan lines */
1111   v2 = (ssize_t)floor(v0 + resample_filter->Vlimit);
1112
1113   /* scan line start and width accross the parallelogram */
1114   u1 = u0 + (v1-v0)*resample_filter->slope - resample_filter->Uwidth;
1115   uw = (ssize_t)(2.0*resample_filter->Uwidth)+1;
1116
1117 #if DEBUG_ELLIPSE
1118   fprintf(stderr, "v1=%ld; v2=%ld\n", (long)v1, (long)v2);
1119   fprintf(stderr, "u1=%ld; uw=%ld\n", (long)u1, (long)uw);
1120 #else
1121 # define DEBUG_HIT_MISS 0 /* only valid if DEBUG_ELLIPSE is enabled */
1122 #endif
1123
1124   /*
1125     Do weighted resampling of all pixels,  within the scaled ellipse,
1126     bound by a Parellelogram fitted to the ellipse.
1127   */
1128   DDQ = 2*resample_filter->A;
1129   for( v=v1; v<=v2;  v++ ) {
1130 #if DEBUG_HIT_MISS
1131     long uu = ceil(u1);   /* actual pixel location (for debug only) */
1132     fprintf(stderr, "# scan line from pixel %ld, %ld\n", (long)uu, (long)v);
1133 #endif
1134     u = (ssize_t)ceil(u1);        /* first pixel in scanline */
1135     u1 += resample_filter->slope; /* start of next scan line */
1136
1137
1138     /* location of this first pixel, relative to u0,v0 */
1139     U = (double)u-u0;
1140     V = (double)v-v0;
1141
1142     /* Q = ellipse quotent ( if Q<F then pixel is inside ellipse) */
1143     Q = (resample_filter->A*U + resample_filter->B*V)*U + resample_filter->C*V*V;
1144     DQ = resample_filter->A*(2.0*U+1) + resample_filter->B*V;
1145
1146     /* get the scanline of pixels for this v */
1147     pixels=GetCacheViewVirtualPixels(resample_filter->view,u,v,(size_t) uw,
1148       1,resample_filter->exception);
1149     if (pixels == (const PixelPacket *) NULL)
1150       return(MagickFalse);
1151     indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
1152
1153     /* count up the weighted pixel colors */
1154     for( u=0; u<uw; u++ ) {
1155 #if FILTER_LUT
1156       /* Note that the ellipse has been pre-scaled so F = WLUT_WIDTH */
1157       if ( Q < (double)WLUT_WIDTH ) {
1158         weight = resample_filter->filter_lut[(int)Q];
1159 #else
1160       /* Note that the ellipse has been pre-scaled so F = support^2 */
1161       if ( Q < (double)resample_filter->F ) {
1162         weight = GetResizeFilterWeight(resample_filter->filter_def,
1163              sqrt(Q));    /* a SquareRoot!  Arrggghhhhh... */
1164 #endif
1165
1166         pixel->opacity  += weight*pixels->opacity;
1167         divisor_m += weight;
1168
1169         if (resample_filter->image->matte != MagickFalse)
1170           weight *= QuantumScale*((MagickRealType)(QuantumRange-pixels->opacity));
1171         pixel->red   += weight*pixels->red;
1172         pixel->green += weight*pixels->green;
1173         pixel->blue  += weight*pixels->blue;
1174         if (resample_filter->image->colorspace == CMYKColorspace)
1175           pixel->index += weight*(*indexes);
1176         divisor_c += weight;
1177
1178         hit++;
1179 #if DEBUG_HIT_MISS
1180         /* mark the pixel according to hit/miss of the ellipse */
1181         fprintf(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
1182                      (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
1183         fprintf(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
1184                      (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
1185       } else {
1186         fprintf(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
1187                      (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
1188         fprintf(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
1189                      (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
1190       }
1191       uu++;
1192 #else
1193       }
1194 #endif
1195       pixels++;
1196       indexes++;
1197       Q += DQ;
1198       DQ += DDQ;
1199     }
1200   }
1201 #if DEBUG_ELLIPSE
1202   fprintf(stderr, "Hit=%ld;  Total=%ld;\n", (long)hit, (long)uw*(v2-v1) );
1203 #endif
1204
1205   /*
1206     Result sanity check -- this should NOT happen
1207   */
1208   if ( hit == 0 ) {
1209     /* not enough pixels in resampling, resort to direct interpolation */
1210 #if DEBUG_NO_PIXEL_HIT
1211     pixel->opacity = pixel->red = pixel->green = pixel->blue = 0;
1212     pixel->red = QuantumRange; /* show pixels for which EWA fails */
1213 #else
1214     status=InterpolateResampleFilter(resample_filter,
1215       resample_filter->interpolate,u0,v0,pixel);
1216 #endif
1217     return status;
1218   }
1219
1220   /*
1221     Finialize results of resampling
1222   */
1223   divisor_m = 1.0/divisor_m;
1224   pixel->opacity = (MagickRealType) ClampToQuantum(divisor_m*pixel->opacity);
1225   divisor_c = 1.0/divisor_c;
1226   pixel->red   = (MagickRealType) ClampToQuantum(divisor_c*pixel->red);
1227   pixel->green = (MagickRealType) ClampToQuantum(divisor_c*pixel->green);
1228   pixel->blue  = (MagickRealType) ClampToQuantum(divisor_c*pixel->blue);
1229   if (resample_filter->image->colorspace == CMYKColorspace)
1230     pixel->index = (MagickRealType) ClampToQuantum(divisor_c*pixel->index);
1231   return(MagickTrue);
1232 }
1233 \f
1234 #if EWA && EWA_CLAMP
1235 /*
1236 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1237 %                                                                             %
1238 %                                                                             %
1239 %                                                                             %
1240 -   C l a m p U p A x e s                                                     %
1241 %                                                                             %
1242 %                                                                             %
1243 %                                                                             %
1244 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1245 %
1246 % ClampUpAxes() function converts the input vectors into a major and
1247 % minor axis unit vectors, and their magnatude.  This form allows us
1248 % to ensure that the ellipse generated is never smaller than the unit
1249 % circle and thus never too small for use in EWA resampling.
1250 %
1251 % This purely mathematical 'magic' was provided by Professor Nicolas
1252 % Robidoux and his Masters student Chantal Racette.
1253 %
1254 % See Reference: "We Recommend Singular Value Decomposition", David Austin
1255 %   http://www.ams.org/samplings/feature-column/fcarc-svd
1256 %
1257 % By generating Major and Minor Axis vectors, we can actually use the
1258 % ellipse in its "canonical form", by remapping the dx,dy of the
1259 % sampled point into distances along the major and minor axis unit
1260 % vectors.
1261 %    http://en.wikipedia.org/wiki/Ellipse#Canonical_form
1262 */
1263 static inline void ClampUpAxes(const double dux,
1264                                const double dvx,
1265                                const double duy,
1266                                const double dvy,
1267                                double *major_mag,
1268                                double *minor_mag,
1269                                double *major_unit_x,
1270                                double *major_unit_y,
1271                                double *minor_unit_x,
1272                                double *minor_unit_y)
1273 {
1274   /*
1275    * ClampUpAxes takes an input 2x2 matrix
1276    *
1277    * [ a b ] = [ dux duy ]
1278    * [ c d ] = [ dvx dvy ]
1279    *
1280    * and computes from it the major and minor axis vectors [major_x,
1281    * major_y] and [minor_x,minor_y] of the smallest ellipse containing
1282    * both the unit disk and the ellipse which is the image of the unit
1283    * disk by the linear transformation
1284    *
1285    * [ dux duy ] [S] = [s]
1286    * [ dvx dvy ] [T] = [t]
1287    *
1288    * (The vector [S,T] is the difference between a position in output
1289    * space and [X,Y]; the vector [s,t] is the difference between a
1290    * position in input space and [x,y].)
1291    */
1292   /*
1293    * Outputs:
1294    *
1295    * major_mag is the half-length of the major axis of the "new"
1296    * ellipse (in input space).
1297    *
1298    * minor_mag is the half-length of the minor axis of the "new"
1299    * ellipse (in input space).
1300    *
1301    * major_unit_x is the x-coordinate of the major axis direction vector
1302    * of both the "old" and "new" ellipses.
1303    *
1304    * major_unit_y is the y-coordinate of the major axis direction vector.
1305    *
1306    * minor_unit_x is the x-coordinate of the minor axis direction vector.
1307    *
1308    * minor_unit_y is the y-coordinate of the minor axis direction vector.
1309    *
1310    * Unit vectors are useful for computing projections, in particular,
1311    * to compute the distance between a point in output space and the
1312    * center (of a disk) from the position of the corresponding point
1313    * in input space.
1314    *
1315    * Now, if you want to modify the input pair of tangent vectors so
1316    * that it defines the modified ellipse, all you have to do is set
1317    *
1318    * newdux = major_mag * major_unit_x
1319    * newdvx = major_mag * major_unit_y
1320    * newduy = minor_mag * minor_unit_x = minor_mag * -major_unit_y
1321    * newdvy = minor_mag * minor_unit_y = minor_mag *  major_unit_x
1322    *
1323    * and use these new tangent vectors "as if" they were the original
1324    * ones.  Most of the time this is a rather drastic change in the
1325    * tangent vectors (even if the singular values are large enough not
1326    * to be clampled). A technical explanation of why things still work
1327    * is found at the end of the discussion below.
1328    *
1329    */
1330   /*
1331    * Discussion:
1332    *
1333    * GOAL: Fix things so that the pullback, in input space, of a disk
1334    * of radius r in output space is an ellipse which contains, at
1335    * least, a disc of radius r. (Make this hold for any r>0.)
1336    *
1337    * SUMMARY OF THE METHOD: Compute the non-unitary factor of the left
1338    * polar decomposition of the linear transformation defining the
1339    * ellipse and make sure that both its columns have norm at least 1.
1340    * Because rotations and reflexions map disks to themselves, it is
1341    * not necessary to compute the other factor of the polar
1342    * decomposition.
1343    *
1344    * DETAILS: Find the singular values and (unit) left singular
1345    * vectors of Jinv, clampling up the singular values to 1, and
1346    * multiplying the unit left singular vectors by the new singular
1347    * values in order to get the minor and major ellipse axis vectors.
1348    *
1349    * Inputs:
1350    *
1351    * The Jacobian matrix of the transformation at the output point
1352    * under consideration is defined as follows:
1353    *
1354    * Consider the transformation (x,y) -> (X,Y) from input locations
1355    * to output locations. (Anthony Thyssen, elsewhere in resample.c,
1356    * uses the notation (u,v) -> (x,y) instead of (x,y) -> (X,Y).)
1357    *
1358    * The Jacobian matrix J is equal to
1359    *
1360    *   [ A, B ] = [ dX/dx, dX/dy ]
1361    *   [ C, D ] = [ dY/dx, dY/dy ]
1362    *
1363    * Consequently, the vector [A,C] is the tangent vector
1364    * corresponding to input changes in the horizontal direction, and
1365    * the vector [B,D] is the tangent vector corresponding to input
1366    * changes in the vertical direction.
1367    *
1368    * In the context of resampling, it is more natural to use the
1369    * inverse Jacobian matrix Jinv. Jinv is
1370    *
1371    *   [ a, b ] = [ dx/dX, dx/dY ]
1372    *   [ c, d ] = [ dy/dX, dy/dY ]
1373    *
1374    * Note: Jinv can be computed from J with the following matrix
1375    * formula:
1376    *
1377    *   Jinv = 1/(A*D-B*C) [  D, -B ]
1378    *                      [ -C,  A ]
1379    *
1380    * What we (implicitly) want to do is replace Jinv by a new Jinv
1381    * which generates an ellipse which is as close as possible to the
1382    * original but which contains the unit disk. This is accomplished
1383    * as follows:
1384    *
1385    * Let
1386    *
1387    *   Jinv = U Sigma V^T
1388    *
1389    * be an SVD decomposition of Jinv. (The SVD is not unique. The
1390    * final ellipse does not depend on the particular SVD. It only
1391    * depends on the hermitian factor of the left polar decomposition,
1392    * which is unique.) In principle, what we want is to clamp up the
1393    * entries of the diagonal matrix Sigma so that they are at least 1,
1394    * and then set
1395    *
1396    *   Jinv = U newSigma V^T.
1397    *
1398    * However, we do not need to compute V^T for the following reason:
1399    * V is an orthogonal matrix (that is, it represents a combination
1400    * of a rotation and a reflexion). Consequently, V maps the unit
1401    * circle to itself. For this reason, the exact value of V does not
1402    * affect the final ellipse, and we choose the identity matrix.
1403    * That is, we simply set
1404    *
1405    *   Jinv = U newSigma,
1406    *
1407    * omitting the V^T factor altogether. In the end, we return the two
1408    * diagonal entries of newSigma together with the two columns of U,
1409    * for a total of six returned quantities.
1410    */
1411   /*
1412    * ClampUpAxes was written by Nicolas Robidoux and Chantal Racette
1413    * of Laurentian University with funding from the National Science
1414    * and Engineering Research Council of Canada.
1415    *
1416    * The idea of using the SVD to clamp the singular values of the
1417    * linear part of the affine approximation of the pullback
1418    * transformation comes from the astrophysicist Craig DeForest, who
1419    * implemented it for use with (approximate) Gaussian filtering in
1420    * his PDL::Transform code (PDL = Perl Data Language).
1421    *
1422    * The only new math in the following is the selection of the
1423    * largest row of the eigen matrix system in order to stabilize the
1424    * computation in near rank-deficient cases, and the corresponding
1425    * efficient repair of degenerate cases using the norm of this
1426    * largest row. Omitting the "V^T" factor of the SVD is also a new
1427    * "trick." It corresponds to moving from the SVD to the left polar
1428    * decomposition.
1429    */
1430   const double a = dux;
1431   const double b = duy;
1432   const double c = dvx;
1433   const double d = dvy;
1434   /*
1435    * n is the matrix Jinv * transpose(Jinv). Eigenvalues of n are the
1436    * squares of the singular values of Jinv.
1437    */
1438   const double aa = a*a;
1439   const double bb = b*b;
1440   const double cc = c*c;
1441   const double dd = d*d;
1442   /*
1443    * Eigenvectors of n are left singular vectors of Jinv.
1444    */
1445   const double n11 = aa+bb;
1446   const double n12 = a*c+b*d;
1447   const double n21 = n12;
1448   const double n22 = cc+dd;
1449   const double det = a*d-b*c;
1450   const double twice_det = det+det;
1451   const double frobenius_squared = n11+n22;
1452   const double discriminant =
1453     (frobenius_squared+twice_det)*(frobenius_squared-twice_det);
1454   const double sqrt_discriminant = sqrt(discriminant);
1455   /*
1456    * s1 is the largest singular value of the inverse Jacobian
1457    * matrix. In other words, its reciprocal is the smallest singular
1458    * value of the Jacobian matrix itself.
1459    * If s1 = 0, both singular values are 0, and any orthogonal pair of
1460    * left and right factors produces a singular decomposition of Jinv.
1461    */
1462   /*
1463    * Initially, we only compute the squares of the singular values.
1464    */
1465   const double s1s1 = 0.5*(frobenius_squared+sqrt_discriminant);
1466   /*
1467    * s2 the smallest singular value of the inverse Jacobian
1468    * matrix. Its reciprocal is the largest singular value of the
1469    * Jacobian matrix itself.
1470    */
1471   const double s2s2 = 0.5*(frobenius_squared-sqrt_discriminant);
1472   const double s1s1minusn11 = s1s1-n11;
1473   const double s1s1minusn22 = s1s1-n22;
1474   /*
1475    * u1, the first column of the U factor of a singular decomposition
1476    * of Jinv, is a (non-normalized) left singular vector corresponding
1477    * to s1. It has entries u11 and u21. We compute u1 from the fact
1478    * that it is an eigenvector of n corresponding to the eigenvalue
1479    * s1^2.
1480    */
1481   const double s1s1minusn11_squared = s1s1minusn11*s1s1minusn11;
1482   const double s1s1minusn22_squared = s1s1minusn22*s1s1minusn22;
1483   /*
1484    * The following selects the largest row of n-s1^2 I as the one
1485    * which is used to find the eigenvector. If both s1^2-n11 and
1486    * s1^2-n22 are zero, n-s1^2 I is the zero matrix.  In that case,
1487    * any vector is an eigenvector; in addition, norm below is equal to
1488    * zero, and, in exact arithmetic, this is the only case in which
1489    * norm = 0. So, setting u1 to the simple but arbitrary vector [1,0]
1490    * if norm = 0 safely takes care of all cases.
1491    */
1492   const double temp_u11 =
1493     ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? n12 : s1s1minusn22 );
1494   const double temp_u21 =
1495     ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? s1s1minusn11 : n21 );
1496   const double norm = sqrt(temp_u11*temp_u11+temp_u21*temp_u21);
1497   /*
1498    * Finalize the entries of first left singular vector (associated
1499    * with the largest singular value).
1500    */
1501   const double u11 = ( (norm>0.0) ? temp_u11/norm : 1.0 );
1502   const double u21 = ( (norm>0.0) ? temp_u21/norm : 0.0 );
1503   /*
1504    * Clamp the singular values up to 1.
1505    */
1506   *major_mag = ( (s1s1<=1.0) ? 1.0 : sqrt(s1s1) );
1507   *minor_mag = ( (s2s2<=1.0) ? 1.0 : sqrt(s2s2) );
1508   /*
1509    * Return the unit major and minor axis direction vectors.
1510    */
1511   *major_unit_x = u11;
1512   *major_unit_y = u21;
1513   *minor_unit_x = -u21;
1514   *minor_unit_y = u11;
1515 }
1516 \f
1517 #endif
1518 /*
1519 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1520 %                                                                             %
1521 %                                                                             %
1522 %                                                                             %
1523 %   S c a l e R e s a m p l e F i l t e r                                     %
1524 %                                                                             %
1525 %                                                                             %
1526 %                                                                             %
1527 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1528 %
1529 %  ScaleResampleFilter() does all the calculations needed to resample an image
1530 %  at a specific scale, defined by two scaling vectors.  This not using
1531 %  a orthogonal scaling, but two distorted scaling vectors, to allow the
1532 %  generation of a angled ellipse.
1533 %
1534 %  As only two deritive scaling vectors are used the center of the ellipse
1535 %  must be the center of the lookup.  That is any curvature that the
1536 %  distortion may produce is discounted.
1537 %
1538 %  The input vectors are produced by either finding the derivitives of the
1539 %  distortion function, or the partial derivitives from a distortion mapping.
1540 %  They do not need to be the orthogonal dx,dy scaling vectors, but can be
1541 %  calculated from other derivatives.  For example you could use  dr,da/r
1542 %  polar coordinate vector scaling vectors
1543 %
1544 %  If   u,v =  DistortEquation(x,y)   OR   u = Fu(x,y); v = Fv(x,y)
1545 %  Then the scaling vectors are determined from the deritives...
1546 %      du/dx, dv/dx     and    du/dy, dv/dy
1547 %  If the resulting scaling vectors is othogonally aligned then...
1548 %      dv/dx = 0   and   du/dy  =  0
1549 %  Producing an othogonally alligned ellipse in source space for the area to
1550 %  be resampled.
1551 %
1552 %  Note that scaling vectors are different to argument order.  Argument order
1553 %  is the general order the deritives are extracted from the distortion
1554 %  equations, and not the scaling vectors. As such the middle two vaules
1555 %  may be swapped from what you expect.  Caution is advised.
1556 %
1557 %  WARNING: It is assumed that any SetResampleFilter() method call will
1558 %  always be performed before the ScaleResampleFilter() method, so that the
1559 %  size of the ellipse will match the support for the resampling filter being
1560 %  used.
1561 %
1562 %  The format of the ScaleResampleFilter method is:
1563 %
1564 %     void ScaleResampleFilter(const ResampleFilter *resample_filter,
1565 %       const double dux,const double duy,const double dvx,const double dvy)
1566 %
1567 %  A description of each parameter follows:
1568 %
1569 %    o resample_filter: the resampling resample_filterrmation defining the
1570 %      image being resampled
1571 %
1572 %    o dux,duy,dvx,dvy:
1573 %         The deritives or scaling vectors defining the EWA ellipse.
1574 %         NOTE: watch the order, which is based on the order deritives
1575 %         are usally determined from distortion equations (see above).
1576 %         The middle two values may need to be swapped if you are thinking
1577 %         in terms of scaling vectors.
1578 %
1579 */
1580 MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
1581   const double dux,const double duy,const double dvx,const double dvy)
1582 {
1583   double A,B,C,F;
1584
1585   assert(resample_filter != (ResampleFilter *) NULL);
1586   assert(resample_filter->signature == MagickSignature);
1587
1588   resample_filter->limit_reached = MagickFalse;
1589
1590   /* A 'point' filter forces use of interpolation instead of area sampling */
1591   if ( resample_filter->filter == PointFilter )
1592     return; /* EWA turned off - nothing to do */
1593
1594 #if DEBUG_ELLIPSE
1595   fprintf(stderr, "# -----\n" );
1596   fprintf(stderr, "dux=%lf; dvx=%lf;   duy=%lf; dvy=%lf;\n",
1597        dux, dvx, duy, dvy);
1598 #endif
1599
1600   /* Find Ellipse Coefficents such that
1601         A*u^2 + B*u*v + C*v^2 = F
1602      With u,v relative to point around which we are resampling.
1603      And the given scaling dx,dy vectors in u,v space
1604          du/dx,dv/dx   and  du/dy,dv/dy
1605   */
1606 #if EWA
1607   /* Direct conversion of derivatives into elliptical coefficients
1608      However when magnifying images, the scaling vectors will be small
1609      resulting in a ellipse that is too small to sample properly.
1610      As such we need to clamp the major/minor axis to a minumum of 1.0
1611      to prevent it getting too small.
1612   */
1613 #if EWA_CLAMP
1614   { double major_mag,
1615            minor_mag,
1616            major_x,
1617            major_y,
1618            minor_x,
1619            minor_y;
1620
1621   ClampUpAxes(dux,dvx,duy,dvy, &major_mag, &minor_mag,
1622                 &major_x, &major_y, &minor_x, &minor_y);
1623   major_x *= major_mag;  major_y *= major_mag;
1624   minor_x *= minor_mag;  minor_y *= minor_mag;
1625 #if DEBUG_ELLIPSE
1626   fprintf(stderr, "major_x=%lf; major_y=%lf;  minor_x=%lf; minor_y=%lf;\n",
1627         major_x, major_y, minor_x, minor_y);
1628 #endif
1629   A = major_y*major_y+minor_y*minor_y;
1630   B = -2.0*(major_x*major_y+minor_x*minor_y);
1631   C = major_x*major_x+minor_x*minor_x;
1632   F = major_mag*minor_mag;
1633   F *= F; /* square it */
1634   }
1635 #else /* raw unclamped EWA */
1636   A = dvx*dvx+dvy*dvy;
1637   B = -2.0*(dux*dvx+duy*dvy);
1638   C = dux*dux+duy*duy;
1639   F = dux*dvy-duy*dvx;
1640   F *= F; /* square it */
1641 #endif /* EWA_CLAMP */
1642
1643 #else /* HQ_EWA */
1644   /*
1645     This Paul Heckbert's "Higher Quality EWA" formula, from page 60 in his
1646     thesis, which adds a unit circle to the elliptical area so as to do both
1647     Reconstruction and Prefiltering of the pixels in the resampling.  It also
1648     means it is always likely to have at least 4 pixels within the area of the
1649     ellipse, for weighted averaging.  No scaling will result with F == 4.0 and
1650     a circle of radius 2.0, and F smaller than this means magnification is
1651     being used.
1652
1653     NOTE: This method produces a very blury result at near unity scale while
1654     producing perfect results for strong minitification and magnifications.
1655
1656     However filter support is fixed to 2.0 (no good for Windowed Sinc filters)
1657   */
1658   A = dvx*dvx+dvy*dvy+1;
1659   B = -2.0*(dux*dvx+duy*dvy);
1660   C = dux*dux+duy*duy+1;
1661   F = A*C - B*B/4;
1662 #endif
1663
1664 #if DEBUG_ELLIPSE
1665   fprintf(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
1666
1667   /* Figure out the various information directly about the ellipse.
1668      This information currently not needed at this time, but may be
1669      needed later for better limit determination.
1670
1671      It is also good to have as a record for future debugging
1672   */
1673   { double alpha, beta, gamma, Major, Minor;
1674     double Eccentricity, Ellipse_Area, Ellipse_Angle;
1675
1676     alpha = A+C;
1677     beta  = A-C;
1678     gamma = sqrt(beta*beta + B*B );
1679
1680     if ( alpha - gamma <= MagickEpsilon )
1681       Major = MagickHuge;
1682     else
1683       Major = sqrt(2*F/(alpha - gamma));
1684     Minor = sqrt(2*F/(alpha + gamma));
1685
1686     fprintf(stderr, "# Major=%lf; Minor=%lf\n", Major, Minor );
1687
1688     /* other information about ellipse include... */
1689     Eccentricity = Major/Minor;
1690     Ellipse_Area = MagickPI*Major*Minor;
1691     Ellipse_Angle = atan2(B, A-C);
1692
1693     fprintf(stderr, "# Angle=%lf   Area=%lf\n",
1694          RadiansToDegrees(Ellipse_Angle), Ellipse_Area);
1695   }
1696 #endif
1697
1698   /* If one or both of the scaling vectors is impossibly large
1699      (producing a very large raw F value), we may as well not bother
1700      doing any form of resampling since resampled area is very large.
1701      In this case some alternative means of pixel sampling, such as
1702      the average of the whole image is needed to get a reasonable
1703      result. Calculate only as needed.
1704   */
1705   if ( (4*A*C - B*B) > MagickHuge ) {
1706     resample_filter->limit_reached = MagickTrue;
1707     return;
1708   }
1709
1710   /* Scale ellipse to match the filters support
1711      (that is, multiply F by the square of the support).
1712   */
1713   F *= resample_filter->support;
1714   F *= resample_filter->support;
1715
1716   /* Orthogonal bounds of the ellipse */
1717   resample_filter->Ulimit = sqrt(4*C*F/(4*A*C-B*B));
1718   resample_filter->Vlimit = sqrt(4*A*F/(4*A*C-B*B));
1719
1720   /* Horizontally aligned parallelogram fitted to Ellipse */
1721   resample_filter->Uwidth = sqrt(F/A); /* Half of the parallelogram width */
1722   resample_filter->slope = -B/(2*A); /* Reciprocal slope of the parallelogram */
1723
1724 #if DEBUG_ELLIPSE
1725   fprintf(stderr, "Ulimit=%lf; Vlimit=%lf; UWidth=%lf; Slope=%lf;\n",
1726            resample_filter->Ulimit, resample_filter->Vlimit,
1727            resample_filter->Uwidth, resample_filter->slope );
1728 #endif
1729
1730   /* Check the absolute area of the parallelogram involved.
1731    * This limit needs more work, as it is too slow for larger images
1732    * with tiled views of the horizon.
1733   */
1734   if ( (resample_filter->Uwidth * resample_filter->Vlimit)
1735          > (4.0*resample_filter->image_area)) {
1736     resample_filter->limit_reached = MagickTrue;
1737     return;
1738   }
1739
1740   /* Scale ellipse formula to directly index the Filter Lookup Table */
1741   { register double scale;
1742 #if FILTER_LUT
1743     /* scale so that F = WLUT_WIDTH; -- hardcoded */
1744     scale = (double)WLUT_WIDTH/F;
1745 #else
1746     /* scale so that F = resample_filter->F (support^2) */
1747     scale = resample_filter->F/F;
1748 #endif
1749     resample_filter->A = A*scale;
1750     resample_filter->B = B*scale;
1751     resample_filter->C = C*scale;
1752   }
1753 }
1754 \f
1755 /*
1756 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1757 %                                                                             %
1758 %                                                                             %
1759 %                                                                             %
1760 %   S e t R e s a m p l e F i l t e r                                         %
1761 %                                                                             %
1762 %                                                                             %
1763 %                                                                             %
1764 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1765 %
1766 %  SetResampleFilter() set the resampling filter lookup table based on a
1767 %  specific filter.  Note that the filter is used as a radial filter not as a
1768 %  two pass othogonally aligned resampling filter.
1769 %
1770 %  The default Filter, is Gaussian, which is the standard filter used by the
1771 %  original paper on the Elliptical Weighted Everage Algorithm. However other
1772 %  filters can also be used.
1773 %
1774 %  The format of the SetResampleFilter method is:
1775 %
1776 %    void SetResampleFilter(ResampleFilter *resample_filter,
1777 %      const FilterTypes filter,const double blur)
1778 %
1779 %  A description of each parameter follows:
1780 %
1781 %    o resample_filter: resampling resample_filterrmation structure
1782 %
1783 %    o filter: the resize filter for elliptical weighting LUT
1784 %
1785 %    o blur: filter blur factor (radial scaling) for elliptical weighting LUT
1786 %
1787 */
1788 MagickExport void SetResampleFilter(ResampleFilter *resample_filter,
1789   const FilterTypes filter,const double blur)
1790 {
1791   ResizeFilter
1792      *resize_filter;
1793
1794   assert(resample_filter != (ResampleFilter *) NULL);
1795   assert(resample_filter->signature == MagickSignature);
1796
1797   resample_filter->do_interpolate = MagickFalse;
1798   resample_filter->filter = filter;
1799
1800   if ( filter == PointFilter )
1801     {
1802       resample_filter->do_interpolate = MagickTrue;
1803       return;  /* EWA turned off - nothing more to do */
1804     }
1805
1806   /* Set a default cylindrical filter of a 'low blur' Jinc windowed Jinc */
1807   if ( filter == UndefinedFilter )
1808     resample_filter->filter = RobidouxFilter;
1809
1810   resize_filter = AcquireResizeFilter(resample_filter->image,
1811        resample_filter->filter,blur,MagickTrue,resample_filter->exception);
1812   if (resize_filter == (ResizeFilter *) NULL)
1813     {
1814       (void) ThrowMagickException(resample_filter->exception,GetMagickModule(),
1815            ModuleError, "UnableToSetFilteringValue",
1816            "Fall back to default EWA gaussian filter");
1817       resample_filter->filter = PointFilter;
1818     }
1819
1820   /* Get the practical working support for the filter,
1821    * after any API call blur factors have been accoded for.
1822    */
1823 #if EWA
1824   resample_filter->support = GetResizeFilterSupport(resize_filter);
1825 #else
1826   resample_filter->support = 2.0;  /* fixed support size for HQ-EWA */
1827 #endif
1828
1829 #if FILTER_LUT
1830   /* Fill the LUT with the weights from the selected filter function */
1831   { register int
1832        Q;
1833     double
1834        r_scale;
1835     /* Scale radius so the filter LUT covers the full support range */
1836     r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
1837     for(Q=0; Q<WLUT_WIDTH; Q++)
1838       resample_filter->filter_lut[Q] = (double)
1839            GetResizeFilterWeight(resize_filter,sqrt((double)Q)*r_scale);
1840
1841     /* finished with the resize filter */
1842     resize_filter = DestroyResizeFilter(resize_filter);
1843   }
1844 #else
1845   /* save the filter and the scaled ellipse bounds needed for filter */
1846   resample_filter->filter_def = resize_filter;
1847   resample_filter->F = resample_filter->support*resample_filter->support;
1848 #endif
1849
1850   /*
1851     Adjust the scaling of the default unit circle
1852     This assumes that any real scaling changes will always
1853     take place AFTER the filter method has been initialized.
1854   */
1855   ScaleResampleFilter(resample_filter, 1.0, 0.0, 0.0, 1.0);
1856
1857 #if 0
1858   /* This is old code kept as a reference only.  It is very wrong,
1859      and I don't understand exactly what it was attempting to do.
1860   */
1861   /*
1862     Create Normal Gaussian 2D Filter Weighted Lookup Table.
1863     A normal EWA guassual lookup would use   exp(Q*ALPHA)
1864     where  Q = distance squared from 0.0 (center) to 1.0 (edge)
1865     and    ALPHA = -4.0*ln(2.0)  ==>  -2.77258872223978123767
1866     The table is of length 1024, and equates to support radius of 2.0
1867     thus needs to be scaled by  ALPHA*4/1024 and any blur factor squared
1868
1869     The above came from some reference code provided by Fred Weinhaus
1870     and seems to have been a guess that was appropriate for its use
1871     in a 3d perspective landscape mapping program.
1872   */
1873   r_scale = -2.77258872223978123767/(WLUT_WIDTH*blur*blur);
1874   for(Q=0; Q<WLUT_WIDTH; Q++)
1875     resample_filter->filter_lut[Q] = exp((double)Q*r_scale);
1876   resample_filter->support = WLUT_WIDTH;
1877   break;
1878 #endif
1879
1880 #if FILTER_LUT
1881 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1882   #pragma omp single
1883 #endif
1884   { register int
1885        Q;
1886     double
1887        r_scale;
1888
1889     /* Scale radius so the filter LUT covers the full support range */
1890     r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
1891     if (IsMagickTrue(GetImageArtifact(resample_filter->image,"resample:verbose")) )
1892       {
1893         /* Debug output of the filter weighting LUT
1894           Gnuplot the LUT with hoizontal adjusted to 'r' using...
1895             plot [0:2][-.2:1] "lut.dat" using (sqrt($0/1024)*2):1 with lines
1896           The filter values is normalized for comparision
1897         */
1898         printf("#\n");
1899         printf("# Resampling Filter LUT (%d values)\n", WLUT_WIDTH);
1900         printf("#\n");
1901         printf("# Note: values in table are using a squared radius lookup.\n");
1902         printf("# And the whole table represents the filters support.\n");
1903         printf("\n"); /* generates a 'break' in gnuplot if multiple outputs */
1904         for(Q=0; Q<WLUT_WIDTH; Q++)
1905           printf("%8.*g %.*g\n",
1906                GetMagickPrecision(),sqrt((double)Q)*r_scale,
1907                GetMagickPrecision(),resample_filter->filter_lut[Q] );
1908       }
1909     /* output the above once only for each image, and each setting */
1910     (void) DeleteImageArtifact(resample_filter->image,"resample:verbose");
1911   }
1912 #endif /* FILTER_LUT */
1913   return;
1914 }
1915 \f
1916 /*
1917 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1918 %                                                                             %
1919 %                                                                             %
1920 %                                                                             %
1921 %   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       %
1922 %                                                                             %
1923 %                                                                             %
1924 %                                                                             %
1925 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1926 %
1927 %  SetResampleFilterInterpolateMethod() changes the interpolation method
1928 %  associated with the specified resample filter.
1929 %
1930 %  The format of the SetResampleFilterInterpolateMethod method is:
1931 %
1932 %      MagickBooleanType SetResampleFilterInterpolateMethod(
1933 %        ResampleFilter *resample_filter,const InterpolateMethod method)
1934 %
1935 %  A description of each parameter follows:
1936 %
1937 %    o resample_filter: the resample filter.
1938 %
1939 %    o method: the interpolation method.
1940 %
1941 */
1942 MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
1943   ResampleFilter *resample_filter,const InterpolatePixelMethod method)
1944 {
1945   assert(resample_filter != (ResampleFilter *) NULL);
1946   assert(resample_filter->signature == MagickSignature);
1947   assert(resample_filter->image != (Image *) NULL);
1948
1949   if (resample_filter->debug != MagickFalse)
1950     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1951       resample_filter->image->filename);
1952
1953   resample_filter->interpolate=method;
1954
1955   return(MagickTrue);
1956 }
1957 \f
1958 /*
1959 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1960 %                                                                             %
1961 %                                                                             %
1962 %                                                                             %
1963 %   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     %
1964 %                                                                             %
1965 %                                                                             %
1966 %                                                                             %
1967 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1968 %
1969 %  SetResampleFilterVirtualPixelMethod() changes the virtual pixel method
1970 %  associated with the specified resample filter.
1971 %
1972 %  The format of the SetResampleFilterVirtualPixelMethod method is:
1973 %
1974 %      MagickBooleanType SetResampleFilterVirtualPixelMethod(
1975 %        ResampleFilter *resample_filter,const VirtualPixelMethod method)
1976 %
1977 %  A description of each parameter follows:
1978 %
1979 %    o resample_filter: the resample filter.
1980 %
1981 %    o method: the virtual pixel method.
1982 %
1983 */
1984 MagickExport MagickBooleanType SetResampleFilterVirtualPixelMethod(
1985   ResampleFilter *resample_filter,const VirtualPixelMethod method)
1986 {
1987   assert(resample_filter != (ResampleFilter *) NULL);
1988   assert(resample_filter->signature == MagickSignature);
1989   assert(resample_filter->image != (Image *) NULL);
1990   if (resample_filter->debug != MagickFalse)
1991     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1992       resample_filter->image->filename);
1993   resample_filter->virtual_pixel=method;
1994   if (method != UndefinedVirtualPixelMethod)
1995     (void) SetCacheViewVirtualPixelMethod(resample_filter->view,method);
1996   return(MagickTrue);
1997 }