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