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-2011 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 magnitude.  This allows us to
1248 % 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 % 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 %
1262 % Reference: http://en.wikipedia.org/wiki/Ellipse#Canonical_form
1263 */
1264 static inline void ClampUpAxes(const double dux,
1265                                const double dvx,
1266                                const double duy,
1267                                const double dvy,
1268                                double *major_mag,
1269                                double *minor_mag,
1270                                double *major_unit_x,
1271                                double *major_unit_y,
1272                                double *minor_unit_x,
1273                                double *minor_unit_y)
1274 {
1275   /*
1276    * ClampUpAxes takes an input 2x2 matrix
1277    *
1278    * [ a b ] = [ dux duy ]
1279    * [ c d ] = [ dvx dvy ]
1280    *
1281    * and computes from it the major and minor axis vectors [major_x,
1282    * major_y] and [minor_x,minor_y] of the smallest ellipse containing
1283    * both the unit disk and the ellipse which is the image of the unit
1284    * disk by the linear transformation
1285    *
1286    * [ dux duy ] [S] = [s]
1287    * [ dvx dvy ] [T] = [t]
1288    *
1289    * (The vector [S,T] is the difference between a position in output
1290    * space and [X,Y]; the vector [s,t] is the difference between a
1291    * position in input space and [x,y].)
1292    */
1293   /*
1294    * Outputs:
1295    *
1296    * major_mag is the half-length of the major axis of the "new"
1297    * ellipse.
1298    *
1299    * minor_mag is the half-length of the minor axis of the "new"
1300    * ellipse.
1301    *
1302    * major_unit_x is the x-coordinate of the major axis direction vector
1303    * of both the "old" and "new" ellipses.
1304    *
1305    * major_unit_y is the y-coordinate of the major axis direction vector.
1306    *
1307    * minor_unit_x is the x-coordinate of the minor axis direction vector.
1308    *
1309    * minor_unit_y is the y-coordinate of the minor axis direction vector.
1310    *
1311    * Unit vectors are useful for computing projections, in particular,
1312    * to compute the distance between a point in output space and the
1313    * center (of a disk) from the position of the corresponding point
1314    * in input space.
1315    *
1316    * Now, if you want to modify the input pair of tangent vectors so
1317    * that it defines the modified ellipse, all you have to do is set
1318    *
1319    * newdux = major_mag * major_unit_x
1320    * newdvx = major_mag * major_unit_y
1321    * newduy = minor_mag * minor_unit_x = minor_mag * -major_unit_y
1322    * newdvy = minor_mag * minor_unit_y = minor_mag *  major_unit_x
1323    *
1324    * and use these tangent vectors as if they were the original ones.
1325    * Usually, this is a drastic change in the tangent vectors even if
1326    * the singular values are not clamped; for example, the minor axis
1327    * vector always points in a direction which is 90 degrees
1328    * counterclockwise from the direction of the major axis vector.
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    * ESSENCE OF THE METHOD: Compute the product of the first two
1338    * factors of an SVD 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 third (rightmost) factor of the SVD.
1342    *
1343    * DETAILS: Find the singular values and (unit) left singular
1344    * vectors of Jinv, clampling up the singular values to 1, and
1345    * multiply the unit left singular vectors by the new singular
1346    * values in order to get the minor and major ellipse axis vectors.
1347    *
1348    * Image resampling context:
1349    *
1350    * The Jacobian matrix of the transformation at the output point
1351    * under consideration is defined as follows:
1352    *
1353    * Consider the transformation (x,y) -> (X,Y) from input locations
1354    * to output locations. (Anthony Thyssen, elsewhere in resample.c,
1355    * uses the notation (u,v) -> (x,y).)
1356    *
1357    * The Jacobian matrix of the transformation at (x,y) is equal to
1358    *
1359    *   J = [ A, B ] = [ dX/dx, dX/dy ]
1360    *       [ C, D ]   [ dY/dx, dY/dy ]
1361    *
1362    * that is, the vector [A,C] is the tangent vector corresponding to
1363    * input changes in the horizontal direction, and the vector [B,D]
1364    * is the tangent vector corresponding to input changes in the
1365    * vertical direction.
1366    *
1367    * In the context of resampling, it is natural to use the inverse
1368    * Jacobian matrix Jinv because resampling is generally performed by
1369    * pulling pixel locations in the output image back to locations in
1370    * the input image. Jinv is
1371    *
1372    *   Jinb = [ a, b ] = [ dx/dX, dx/dY ]
1373    *          [ c, d ]   [ dy/dX, dy/dY ]
1374    *
1375    * Note: Jinv can be computed from J with the following matrix
1376    * formula:
1377    *
1378    *   Jinv = 1/(A*D-B*C) [  D, -B ]
1379    *                      [ -C,  A ]
1380    *
1381    * What we do is modify Jinv so that it generates an ellipse which
1382    * is as close as possible to the original but which contains the
1383    * unit disk. This can be accomplished 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, but the
1390    * final ellipse does not depend on the particular SVD.)
1391    *
1392    * We could clamp up the entries of the diagonal matrix Sigma so
1393    * that they are at least 1, and then set
1394    *
1395    *   Jinv = U newSigma V^T.
1396    *
1397    * However, we do not need to compute V for the following reason:
1398    * V^T is an orthogonal matrix (that is, it represents a combination
1399    * of rotations and reflexions) so that it maps the unit circle to
1400    * itself. For this reason, the exact value of V does not affect the
1401    * final ellipse, and we can choose V to be the identity
1402    * matrix. This gives
1403    *
1404    *   Jinv = U newSigma.
1405    *
1406    * In the end, we return the two diagonal entries of newSigma
1407    * together with the two columns of U.
1408    */
1409   /*
1410    * ClampUpAxes was written by Nicolas Robidoux and Chantal Racette
1411    * of Laurentian University with insightful suggestions from Anthony
1412    * Thyssen and funding from the National Science and Engineering
1413    * Research Council of Canada. It is distinguished from its
1414    * predecessors by its efficient handling of degenerate cases.
1415    *
1416    * The idea of clamping up the EWA ellipse's major and minor axes so
1417    * that the result contains the reconstruction kernel filter support
1418    * is taken from Andreas Gustaffson's Masters thesis "Interactive
1419    * Image Warping", Helsinki University of Technology, Faculty of
1420    * Information Technology, 59 pages, 1993 (see Section 3.6).
1421    *
1422    * The use of the SVD to clamp up the singular values of the
1423    * Jacobian matrix of the pullback transformation for EWA resampling
1424    * is taken from the astrophysicist Craig DeForest.  It is
1425    * implemented in his PDL::Transform code (PDL = Perl Data
1426    * Language).
1427    */
1428   const double a = dux;
1429   const double b = duy;
1430   const double c = dvx;
1431   const double d = dvy;
1432   /*
1433    * n is the matrix Jinv * transpose(Jinv). Eigenvalues of n are the
1434    * squares of the singular values of Jinv.
1435    */
1436   const double aa = a*a;
1437   const double bb = b*b;
1438   const double cc = c*c;
1439   const double dd = d*d;
1440   /*
1441    * Eigenvectors of n are left singular vectors of Jinv.
1442    */
1443   const double n11 = aa+bb;
1444   const double n12 = a*c+b*d;
1445   const double n21 = n12;
1446   const double n22 = cc+dd;
1447   const double det = a*d-b*c;
1448   const double twice_det = det+det;
1449   const double frobenius_squared = n11+n22;
1450   const double discriminant =
1451     (frobenius_squared+twice_det)*(frobenius_squared-twice_det);
1452   const double sqrt_discriminant = sqrt(discriminant);
1453   /*
1454    * s1 is the largest singular value of the inverse Jacobian
1455    * matrix. In other words, its reciprocal is the smallest singular
1456    * value of the Jacobian matrix itself.
1457    * If s1 = 0, both singular values are 0, and any orthogonal pair of
1458    * left and right factors produces a singular decomposition of Jinv.
1459    */
1460   /*
1461    * Initially, we only compute the squares of the singular values.
1462    */
1463   const double s1s1 = 0.5*(frobenius_squared+sqrt_discriminant);
1464   /*
1465    * s2 the smallest singular value of the inverse Jacobian
1466    * matrix. Its reciprocal is the largest singular value of the
1467    * Jacobian matrix itself.
1468    */
1469   const double s2s2 = 0.5*(frobenius_squared-sqrt_discriminant);
1470   const double s1s1minusn11 = s1s1-n11;
1471   const double s1s1minusn22 = s1s1-n22;
1472   /*
1473    * u1, the first column of the U factor of a singular decomposition
1474    * of Jinv, is a (non-normalized) left singular vector corresponding
1475    * to s1. It has entries u11 and u21. We compute u1 from the fact
1476    * that it is an eigenvector of n corresponding to the eigenvalue
1477    * s1^2.
1478    */
1479   const double s1s1minusn11_squared = s1s1minusn11*s1s1minusn11;
1480   const double s1s1minusn22_squared = s1s1minusn22*s1s1minusn22;
1481   /*
1482    * The following selects the largest row of n-s1^2 I as the one
1483    * which is used to find the eigenvector. If both s1^2-n11 and
1484    * s1^2-n22 are zero, n-s1^2 I is the zero matrix.  In that case,
1485    * any vector is an eigenvector; in addition, norm below is equal to
1486    * zero, and, in exact arithmetic, this is the only case in which
1487    * norm = 0. So, setting u1 to the simple but arbitrary vector [1,0]
1488    * if norm = 0 safely takes care of all cases.
1489    */
1490   const double temp_u11 =
1491     ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? n12 : s1s1minusn22 );
1492   const double temp_u21 =
1493     ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? s1s1minusn11 : n21 );
1494   const double norm = sqrt(temp_u11*temp_u11+temp_u21*temp_u21);
1495   /*
1496    * Finalize the entries of first left singular vector (associated
1497    * with the largest singular value).
1498    */
1499   const double u11 = ( (norm>0.0) ? temp_u11/norm : 1.0 );
1500   const double u21 = ( (norm>0.0) ? temp_u21/norm : 0.0 );
1501   /*
1502    * Clamp the singular values up to 1.
1503    */
1504   *major_mag = ( (s1s1<=1.0) ? 1.0 : sqrt(s1s1) );
1505   *minor_mag = ( (s2s2<=1.0) ? 1.0 : sqrt(s2s2) );
1506   /*
1507    * Return the unit major and minor axis direction vectors.
1508    */
1509   *major_unit_x = u11;
1510   *major_unit_y = u21;
1511   *minor_unit_x = -u21;
1512   *minor_unit_y = u11;
1513 }
1514 \f
1515 #endif
1516 /*
1517 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1518 %                                                                             %
1519 %                                                                             %
1520 %                                                                             %
1521 %   S c a l e R e s a m p l e F i l t e r                                     %
1522 %                                                                             %
1523 %                                                                             %
1524 %                                                                             %
1525 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1526 %
1527 %  ScaleResampleFilter() does all the calculations needed to resample an image
1528 %  at a specific scale, defined by two scaling vectors.  This not using
1529 %  a orthogonal scaling, but two distorted scaling vectors, to allow the
1530 %  generation of a angled ellipse.
1531 %
1532 %  As only two deritive scaling vectors are used the center of the ellipse
1533 %  must be the center of the lookup.  That is any curvature that the
1534 %  distortion may produce is discounted.
1535 %
1536 %  The input vectors are produced by either finding the derivitives of the
1537 %  distortion function, or the partial derivitives from a distortion mapping.
1538 %  They do not need to be the orthogonal dx,dy scaling vectors, but can be
1539 %  calculated from other derivatives.  For example you could use  dr,da/r
1540 %  polar coordinate vector scaling vectors
1541 %
1542 %  If   u,v =  DistortEquation(x,y)   OR   u = Fu(x,y); v = Fv(x,y)
1543 %  Then the scaling vectors are determined from the deritives...
1544 %      du/dx, dv/dx     and    du/dy, dv/dy
1545 %  If the resulting scaling vectors is othogonally aligned then...
1546 %      dv/dx = 0   and   du/dy  =  0
1547 %  Producing an othogonally alligned ellipse in source space for the area to
1548 %  be resampled.
1549 %
1550 %  Note that scaling vectors are different to argument order.  Argument order
1551 %  is the general order the deritives are extracted from the distortion
1552 %  equations, and not the scaling vectors. As such the middle two vaules
1553 %  may be swapped from what you expect.  Caution is advised.
1554 %
1555 %  WARNING: It is assumed that any SetResampleFilter() method call will
1556 %  always be performed before the ScaleResampleFilter() method, so that the
1557 %  size of the ellipse will match the support for the resampling filter being
1558 %  used.
1559 %
1560 %  The format of the ScaleResampleFilter method is:
1561 %
1562 %     void ScaleResampleFilter(const ResampleFilter *resample_filter,
1563 %       const double dux,const double duy,const double dvx,const double dvy)
1564 %
1565 %  A description of each parameter follows:
1566 %
1567 %    o resample_filter: the resampling resample_filterrmation defining the
1568 %      image being resampled
1569 %
1570 %    o dux,duy,dvx,dvy:
1571 %         The deritives or scaling vectors defining the EWA ellipse.
1572 %         NOTE: watch the order, which is based on the order deritives
1573 %         are usally determined from distortion equations (see above).
1574 %         The middle two values may need to be swapped if you are thinking
1575 %         in terms of scaling vectors.
1576 %
1577 */
1578 MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
1579   const double dux,const double duy,const double dvx,const double dvy)
1580 {
1581   double A,B,C,F;
1582
1583   assert(resample_filter != (ResampleFilter *) NULL);
1584   assert(resample_filter->signature == MagickSignature);
1585
1586   resample_filter->limit_reached = MagickFalse;
1587
1588   /* A 'point' filter forces use of interpolation instead of area sampling */
1589   if ( resample_filter->filter == PointFilter )
1590     return; /* EWA turned off - nothing to do */
1591
1592 #if DEBUG_ELLIPSE
1593   fprintf(stderr, "# -----\n" );
1594   fprintf(stderr, "dux=%lf; dvx=%lf;   duy=%lf; dvy=%lf;\n",
1595        dux, dvx, duy, dvy);
1596 #endif
1597
1598   /* Find Ellipse Coefficents such that
1599         A*u^2 + B*u*v + C*v^2 = F
1600      With u,v relative to point around which we are resampling.
1601      And the given scaling dx,dy vectors in u,v space
1602          du/dx,dv/dx   and  du/dy,dv/dy
1603   */
1604 #if EWA
1605   /* Direct conversion of derivatives into elliptical coefficients
1606      However when magnifying images, the scaling vectors will be small
1607      resulting in a ellipse that is too small to sample properly.
1608      As such we need to clamp the major/minor axis to a minumum of 1.0
1609      to prevent it getting too small.
1610   */
1611 #if EWA_CLAMP
1612   { double major_mag,
1613            minor_mag,
1614            major_x,
1615            major_y,
1616            minor_x,
1617            minor_y;
1618
1619   ClampUpAxes(dux,dvx,duy,dvy, &major_mag, &minor_mag,
1620                 &major_x, &major_y, &minor_x, &minor_y);
1621   major_x *= major_mag;  major_y *= major_mag;
1622   minor_x *= minor_mag;  minor_y *= minor_mag;
1623 #if DEBUG_ELLIPSE
1624   fprintf(stderr, "major_x=%lf; major_y=%lf;  minor_x=%lf; minor_y=%lf;\n",
1625         major_x, major_y, minor_x, minor_y);
1626 #endif
1627   A = major_y*major_y+minor_y*minor_y;
1628   B = -2.0*(major_x*major_y+minor_x*minor_y);
1629   C = major_x*major_x+minor_x*minor_x;
1630   F = major_mag*minor_mag;
1631   F *= F; /* square it */
1632   }
1633 #else /* raw unclamped EWA */
1634   A = dvx*dvx+dvy*dvy;
1635   B = -2.0*(dux*dvx+duy*dvy);
1636   C = dux*dux+duy*duy;
1637   F = dux*dvy-duy*dvx;
1638   F *= F; /* square it */
1639 #endif /* EWA_CLAMP */
1640
1641 #else /* HQ_EWA */
1642   /*
1643     This Paul Heckbert's "Higher Quality EWA" formula, from page 60 in his
1644     thesis, which adds a unit circle to the elliptical area so as to do both
1645     Reconstruction and Prefiltering of the pixels in the resampling.  It also
1646     means it is always likely to have at least 4 pixels within the area of the
1647     ellipse, for weighted averaging.  No scaling will result with F == 4.0 and
1648     a circle of radius 2.0, and F smaller than this means magnification is
1649     being used.
1650
1651     NOTE: This method produces a very blury result at near unity scale while
1652     producing perfect results for strong minitification and magnifications.
1653
1654     However filter support is fixed to 2.0 (no good for Windowed Sinc filters)
1655   */
1656   A = dvx*dvx+dvy*dvy+1;
1657   B = -2.0*(dux*dvx+duy*dvy);
1658   C = dux*dux+duy*duy+1;
1659   F = A*C - B*B/4;
1660 #endif
1661
1662 #if DEBUG_ELLIPSE
1663   fprintf(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
1664
1665   /* Figure out the various information directly about the ellipse.
1666      This information currently not needed at this time, but may be
1667      needed later for better limit determination.
1668
1669      It is also good to have as a record for future debugging
1670   */
1671   { double alpha, beta, gamma, Major, Minor;
1672     double Eccentricity, Ellipse_Area, Ellipse_Angle;
1673
1674     alpha = A+C;
1675     beta  = A-C;
1676     gamma = sqrt(beta*beta + B*B );
1677
1678     if ( alpha - gamma <= MagickEpsilon )
1679       Major = MagickHuge;
1680     else
1681       Major = sqrt(2*F/(alpha - gamma));
1682     Minor = sqrt(2*F/(alpha + gamma));
1683
1684     fprintf(stderr, "# Major=%lf; Minor=%lf\n", Major, Minor );
1685
1686     /* other information about ellipse include... */
1687     Eccentricity = Major/Minor;
1688     Ellipse_Area = MagickPI*Major*Minor;
1689     Ellipse_Angle = atan2(B, A-C);
1690
1691     fprintf(stderr, "# Angle=%lf   Area=%lf\n",
1692          RadiansToDegrees(Ellipse_Angle), Ellipse_Area);
1693   }
1694 #endif
1695
1696   /* If one or both of the scaling vectors is impossibly large
1697      (producing a very large raw F value), we may as well not bother
1698      doing any form of resampling since resampled area is very large.
1699      In this case some alternative means of pixel sampling, such as
1700      the average of the whole image is needed to get a reasonable
1701      result. Calculate only as needed.
1702   */
1703   if ( (4*A*C - B*B) > MagickHuge ) {
1704     resample_filter->limit_reached = MagickTrue;
1705     return;
1706   }
1707
1708   /* Scale ellipse to match the filters support
1709      (that is, multiply F by the square of the support).
1710   */
1711   F *= resample_filter->support;
1712   F *= resample_filter->support;
1713
1714   /* Orthogonal bounds of the ellipse */
1715   resample_filter->Ulimit = sqrt(4*C*F/(4*A*C-B*B));
1716   resample_filter->Vlimit = sqrt(4*A*F/(4*A*C-B*B));
1717
1718   /* Horizontally aligned parallelogram fitted to Ellipse */
1719   resample_filter->Uwidth = sqrt(F/A); /* Half of the parallelogram width */
1720   resample_filter->slope = -B/(2*A); /* Reciprocal slope of the parallelogram */
1721
1722 #if DEBUG_ELLIPSE
1723   fprintf(stderr, "Ulimit=%lf; Vlimit=%lf; UWidth=%lf; Slope=%lf;\n",
1724            resample_filter->Ulimit, resample_filter->Vlimit,
1725            resample_filter->Uwidth, resample_filter->slope );
1726 #endif
1727
1728   /* Check the absolute area of the parallelogram involved.
1729    * This limit needs more work, as it is too slow for larger images
1730    * with tiled views of the horizon.
1731   */
1732   if ( (resample_filter->Uwidth * resample_filter->Vlimit)
1733          > (4.0*resample_filter->image_area)) {
1734     resample_filter->limit_reached = MagickTrue;
1735     return;
1736   }
1737
1738   /* Scale ellipse formula to directly index the Filter Lookup Table */
1739   { register double scale;
1740 #if FILTER_LUT
1741     /* scale so that F = WLUT_WIDTH; -- hardcoded */
1742     scale = (double)WLUT_WIDTH/F;
1743 #else
1744     /* scale so that F = resample_filter->F (support^2) */
1745     scale = resample_filter->F/F;
1746 #endif
1747     resample_filter->A = A*scale;
1748     resample_filter->B = B*scale;
1749     resample_filter->C = C*scale;
1750   }
1751 }
1752 \f
1753 /*
1754 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1755 %                                                                             %
1756 %                                                                             %
1757 %                                                                             %
1758 %   S e t R e s a m p l e F i l t e r                                         %
1759 %                                                                             %
1760 %                                                                             %
1761 %                                                                             %
1762 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1763 %
1764 %  SetResampleFilter() set the resampling filter lookup table based on a
1765 %  specific filter.  Note that the filter is used as a radial filter not as a
1766 %  two pass othogonally aligned resampling filter.
1767 %
1768 %  The default Filter, is Gaussian, which is the standard filter used by the
1769 %  original paper on the Elliptical Weighted Everage Algorithm. However other
1770 %  filters can also be used.
1771 %
1772 %  The format of the SetResampleFilter method is:
1773 %
1774 %    void SetResampleFilter(ResampleFilter *resample_filter,
1775 %      const FilterTypes filter,const double blur)
1776 %
1777 %  A description of each parameter follows:
1778 %
1779 %    o resample_filter: resampling resample_filterrmation structure
1780 %
1781 %    o filter: the resize filter for elliptical weighting LUT
1782 %
1783 %    o blur: filter blur factor (radial scaling) for elliptical weighting LUT
1784 %
1785 */
1786 MagickExport void SetResampleFilter(ResampleFilter *resample_filter,
1787   const FilterTypes filter,const double blur)
1788 {
1789   ResizeFilter
1790      *resize_filter;
1791
1792   assert(resample_filter != (ResampleFilter *) NULL);
1793   assert(resample_filter->signature == MagickSignature);
1794
1795   resample_filter->do_interpolate = MagickFalse;
1796   resample_filter->filter = filter;
1797
1798   if ( filter == PointFilter )
1799     {
1800       resample_filter->do_interpolate = MagickTrue;
1801       return;  /* EWA turned off - nothing more to do */
1802     }
1803
1804   /* Set a default cylindrical filter of a 'low blur' Jinc windowed Jinc */
1805   if ( filter == UndefinedFilter )
1806     resample_filter->filter = RobidouxFilter;
1807
1808   resize_filter = AcquireResizeFilter(resample_filter->image,
1809        resample_filter->filter,blur,MagickTrue,resample_filter->exception);
1810   if (resize_filter == (ResizeFilter *) NULL)
1811     {
1812       (void) ThrowMagickException(resample_filter->exception,GetMagickModule(),
1813            ModuleError, "UnableToSetFilteringValue",
1814            "Fall back to default EWA gaussian filter");
1815       resample_filter->filter = PointFilter;
1816     }
1817
1818   /* Get the practical working support for the filter,
1819    * after any API call blur factors have been accoded for.
1820    */
1821 #if EWA
1822   resample_filter->support = GetResizeFilterSupport(resize_filter);
1823 #else
1824   resample_filter->support = 2.0;  /* fixed support size for HQ-EWA */
1825 #endif
1826
1827 #if FILTER_LUT
1828   /* Fill the LUT with the weights from the selected filter function */
1829   { register int
1830        Q;
1831     double
1832        r_scale;
1833     /* Scale radius so the filter LUT covers the full support range */
1834     r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
1835     for(Q=0; Q<WLUT_WIDTH; Q++)
1836       resample_filter->filter_lut[Q] = (double)
1837            GetResizeFilterWeight(resize_filter,sqrt((double)Q)*r_scale);
1838
1839     /* finished with the resize filter */
1840     resize_filter = DestroyResizeFilter(resize_filter);
1841   }
1842 #else
1843   /* save the filter and the scaled ellipse bounds needed for filter */
1844   resample_filter->filter_def = resize_filter;
1845   resample_filter->F = resample_filter->support*resample_filter->support;
1846 #endif
1847
1848   /*
1849     Adjust the scaling of the default unit circle
1850     This assumes that any real scaling changes will always
1851     take place AFTER the filter method has been initialized.
1852   */
1853   ScaleResampleFilter(resample_filter, 1.0, 0.0, 0.0, 1.0);
1854
1855 #if 0
1856   /* This is old code kept as a reference only.  It is very wrong,
1857      and I don't understand exactly what it was attempting to do.
1858   */
1859   /*
1860     Create Normal Gaussian 2D Filter Weighted Lookup Table.
1861     A normal EWA guassual lookup would use   exp(Q*ALPHA)
1862     where  Q = distance squared from 0.0 (center) to 1.0 (edge)
1863     and    ALPHA = -4.0*ln(2.0)  ==>  -2.77258872223978123767
1864     The table is of length 1024, and equates to support radius of 2.0
1865     thus needs to be scaled by  ALPHA*4/1024 and any blur factor squared
1866
1867     The above came from some reference code provided by Fred Weinhaus
1868     and seems to have been a guess that was appropriate for its use
1869     in a 3d perspective landscape mapping program.
1870   */
1871   r_scale = -2.77258872223978123767/(WLUT_WIDTH*blur*blur);
1872   for(Q=0; Q<WLUT_WIDTH; Q++)
1873     resample_filter->filter_lut[Q] = exp((double)Q*r_scale);
1874   resample_filter->support = WLUT_WIDTH;
1875   break;
1876 #endif
1877
1878 #if FILTER_LUT
1879 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1880   #pragma omp single
1881 #endif
1882   { register int
1883        Q;
1884     double
1885        r_scale;
1886
1887     /* Scale radius so the filter LUT covers the full support range */
1888     r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
1889     if (IsMagickTrue(GetImageArtifact(resample_filter->image,"resample:verbose")) )
1890       {
1891         /* Debug output of the filter weighting LUT
1892           Gnuplot the LUT with hoizontal adjusted to 'r' using...
1893             plot [0:2][-.2:1] "lut.dat" using (sqrt($0/1024)*2):1 with lines
1894           The filter values is normalized for comparision
1895         */
1896         printf("#\n");
1897         printf("# Resampling Filter LUT (%d values)\n", WLUT_WIDTH);
1898         printf("#\n");
1899         printf("# Note: values in table are using a squared radius lookup.\n");
1900         printf("# And the whole table represents the filters support.\n");
1901         printf("\n"); /* generates a 'break' in gnuplot if multiple outputs */
1902         for(Q=0; Q<WLUT_WIDTH; Q++)
1903           printf("%8.*g %.*g\n",
1904                GetMagickPrecision(),sqrt((double)Q)*r_scale,
1905                GetMagickPrecision(),resample_filter->filter_lut[Q] );
1906       }
1907     /* output the above once only for each image, and each setting */
1908     (void) DeleteImageArtifact(resample_filter->image,"resample:verbose");
1909   }
1910 #endif /* FILTER_LUT */
1911   return;
1912 }
1913 \f
1914 /*
1915 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1916 %                                                                             %
1917 %                                                                             %
1918 %                                                                             %
1919 %   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       %
1920 %                                                                             %
1921 %                                                                             %
1922 %                                                                             %
1923 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1924 %
1925 %  SetResampleFilterInterpolateMethod() changes the interpolation method
1926 %  associated with the specified resample filter.
1927 %
1928 %  The format of the SetResampleFilterInterpolateMethod method is:
1929 %
1930 %      MagickBooleanType SetResampleFilterInterpolateMethod(
1931 %        ResampleFilter *resample_filter,const InterpolateMethod method)
1932 %
1933 %  A description of each parameter follows:
1934 %
1935 %    o resample_filter: the resample filter.
1936 %
1937 %    o method: the interpolation method.
1938 %
1939 */
1940 MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
1941   ResampleFilter *resample_filter,const InterpolatePixelMethod method)
1942 {
1943   assert(resample_filter != (ResampleFilter *) NULL);
1944   assert(resample_filter->signature == MagickSignature);
1945   assert(resample_filter->image != (Image *) NULL);
1946
1947   if (resample_filter->debug != MagickFalse)
1948     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1949       resample_filter->image->filename);
1950
1951   resample_filter->interpolate=method;
1952
1953   return(MagickTrue);
1954 }
1955 \f
1956 /*
1957 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1958 %                                                                             %
1959 %                                                                             %
1960 %                                                                             %
1961 %   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     %
1962 %                                                                             %
1963 %                                                                             %
1964 %                                                                             %
1965 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1966 %
1967 %  SetResampleFilterVirtualPixelMethod() changes the virtual pixel method
1968 %  associated with the specified resample filter.
1969 %
1970 %  The format of the SetResampleFilterVirtualPixelMethod method is:
1971 %
1972 %      MagickBooleanType SetResampleFilterVirtualPixelMethod(
1973 %        ResampleFilter *resample_filter,const VirtualPixelMethod method)
1974 %
1975 %  A description of each parameter follows:
1976 %
1977 %    o resample_filter: the resample filter.
1978 %
1979 %    o method: the virtual pixel method.
1980 %
1981 */
1982 MagickExport MagickBooleanType SetResampleFilterVirtualPixelMethod(
1983   ResampleFilter *resample_filter,const VirtualPixelMethod method)
1984 {
1985   assert(resample_filter != (ResampleFilter *) NULL);
1986   assert(resample_filter->signature == MagickSignature);
1987   assert(resample_filter->image != (Image *) NULL);
1988   if (resample_filter->debug != MagickFalse)
1989     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1990       resample_filter->image->filename);
1991   resample_filter->virtual_pixel=method;
1992   if (method != UndefinedVirtualPixelMethod)
1993     (void) SetCacheViewVirtualPixelMethod(resample_filter->view,method);
1994   return(MagickTrue);
1995 }