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 /*
  64   Typedef declarations.
  65 */
  66 #define WLUT_WIDTH 1024
  67 struct _ResampleFilter
  68 {
  69   CacheView
  70     *view;
  71
  72   Image
  73     *image;
  74
  75   ExceptionInfo
  76     *exception;
  77
  78   MagickBooleanType
  79     debug;
  80
  81   /* Information about image being resampled */
  82   ssize_t
  83     image_area;
  84
  85   InterpolatePixelMethod
  86     interpolate;
  87
  88   VirtualPixelMethod
  89     virtual_pixel;
  90
  91   FilterTypes
  92     filter;
  93
  94   /* processing settings needed */
  95   MagickBooleanType
  96     limit_reached,
  97     do_interpolate,
  98     average_defined;
  99
 100   MagickPixelPacket
 101     average_pixel;
 102
 103   /* current ellipitical area being resampled around center point */
 104   double
 105     A, B, C,
 106     Vlimit, Ulimit, Uwidth, slope;
 107
 108   /* LUT of weights for filtered average in elliptical area */
 109   double
 110     filter_lut[WLUT_WIDTH],
 111     support;
 112
 113   size_t
 114     signature;
 115 };
 116 \f
 117 /*
 118 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 119 %                                                                             %
 120 %                                                                             %
 121 %                                                                             %
 122 %   A c q u i r e R e s a m p l e I n f o                                     %
 123 %                                                                             %
 124 %                                                                             %
 125 %                                                                             %
 126 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 127 %
 128 %  AcquireResampleFilter() initializes the information resample needs do to a
 129 %  scaled lookup of a color from an image, using area sampling.
 130 %
 131 %  The algorithm is based on a Elliptical Weighted Average, where the pixels
 132 %  found in a large elliptical area is averaged together according to a
 133 %  weighting (filter) function.  For more details see "Fundamentals of Texture
 134 %  Mapping and Image Warping" a master's thesis by Paul.S.Heckbert, June 17,
 135 %  1989.  Available for free from, http://www.cs.cmu.edu/~ph/
 136 %
 137 %  As EWA resampling (or any sort of resampling) can require a lot of
 138 %  calculations to produce a distorted scaling of the source image for each
 139 %  output pixel, the ResampleFilter structure generated holds that information
 140 %  between individual image resampling.
 141 %
 142 %  This function will make the appropriate AcquireCacheView() calls
 143 %  to view the image, calling functions do not need to open a cache view.
 144 %
 145 %  Usage Example...
 146 %      resample_filter=AcquireResampleFilter(image,exception);
 147 %      for (y=0; y < (ssize_t) image->rows; y++) {
 148 %        for (x=0; x < (ssize_t) image->columns; x++) {
 149 %          X= ....;   Y= ....;
 150 %          ScaleResampleFilter(resample_filter, ... scaling vectors ...);
 151 %          (void) ResamplePixelColor(resample_filter,X,Y,&pixel);
 152 %          ... assign resampled pixel value ...
 153 %        }
 154 %      }
 155 %      DestroyResampleFilter(resample_filter);
 156 %
 157 %  The format of the AcquireResampleFilter method is:
 158 %
 159 %     ResampleFilter *AcquireResampleFilter(const Image *image,
 160 %       ExceptionInfo *exception)
 161 %
 162 %  A description of each parameter follows:
 163 %
 164 %    o image: the image.
 165 %
 166 %    o exception: return any errors or warnings in this structure.
 167 %
 168 */
 169 MagickExport ResampleFilter *AcquireResampleFilter(const Image *image,
 170   ExceptionInfo *exception)
 171 {
 172   register ResampleFilter
 173     *resample_filter;
 174
 175   assert(image != (Image *) NULL);
 176   assert(image->signature == MagickSignature);
 177   if (image->debug != MagickFalse)
 178     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 179   assert(exception != (ExceptionInfo *) NULL);
 180   assert(exception->signature == MagickSignature);
 181
 182   resample_filter=(ResampleFilter *) AcquireMagickMemory(
 183     sizeof(*resample_filter));
 184   if (resample_filter == (ResampleFilter *) NULL)
 185     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 186   (void) ResetMagickMemory(resample_filter,0,sizeof(*resample_filter));
 187
 188   resample_filter->image=ReferenceImage((Image *) image);
 189   resample_filter->view=AcquireCacheView(resample_filter->image);
 190   resample_filter->exception=exception;
 191
 192   resample_filter->debug=IsEventLogging();
 193   resample_filter->signature=MagickSignature;
 194
 195   resample_filter->image_area=(ssize_t) (resample_filter->image->columns*
 196     resample_filter->image->rows);
 197   resample_filter->average_defined = MagickFalse;
 198
 199   /* initialise the resampling filter settings */
 200   SetResampleFilter(resample_filter, resample_filter->image->filter,
 201     resample_filter->image->blur);
 202   resample_filter->interpolate = resample_filter->image->interpolate;
 203   resample_filter->virtual_pixel=GetImageVirtualPixelMethod(image);
 204
 205   /* init scale to a default of a unit circle */
 206   ScaleResampleFilter(resample_filter, 1.0, 0.0, 0.0, 1.0);
 207
 208   return(resample_filter);
 209 }
 210 \f
 211 /*
 212 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 213 %                                                                             %
 214 %                                                                             %
 215 %                                                                             %
 216 %   D e s t r o y R e s a m p l e I n f o                                     %
 217 %                                                                             %
 218 %                                                                             %
 219 %                                                                             %
 220 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 221 %
 222 %  DestroyResampleFilter() finalizes and cleans up the resampling
 223 %  resample_filter as returned by AcquireResampleFilter(), freeing any memory
 224 %  or other information as needed.
 225 %
 226 %  The format of the DestroyResampleFilter method is:
 227 %
 228 %      ResampleFilter *DestroyResampleFilter(ResampleFilter *resample_filter)
 229 %
 230 %  A description of each parameter follows:
 231 %
 232 %    o resample_filter: resampling information structure
 233 %
 234 */
 235 MagickExport ResampleFilter *DestroyResampleFilter(
 236   ResampleFilter *resample_filter)
 237 {
 238   assert(resample_filter != (ResampleFilter *) NULL);
 239   assert(resample_filter->signature == MagickSignature);
 240   assert(resample_filter->image != (Image *) NULL);
 241   if (resample_filter->debug != MagickFalse)
 242     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
 243       resample_filter->image->filename);
 244   resample_filter->view=DestroyCacheView(resample_filter->view);
 245   resample_filter->image=DestroyImage(resample_filter->image);
 246   resample_filter->signature=(~MagickSignature);
 247   resample_filter=(ResampleFilter *) RelinquishMagickMemory(resample_filter);
 248   return(resample_filter);
 249 }
 250 \f
 251 /*
 252 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 253 %                                                                             %
 254 %                                                                             %
 255 %                                                                             %
 256 %   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                         %
 257 %                                                                             %
 258 %                                                                             %
 259 %                                                                             %
 260 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 261 %
 262 %  InterpolateResampleFilter() applies bi-linear or tri-linear interpolation
 263 %  between a floating point coordinate and the pixels surrounding that
 264 %  coordinate.  No pixel area resampling, or scaling of the result is
 265 %  performed.
 266 %
 267 %  The format of the InterpolateResampleFilter method is:
 268 %
 269 %      MagickBooleanType InterpolateResampleFilter(
 270 %        ResampleInfo *resample_filter,const InterpolatePixelMethod method,
 271 %        const double x,const double y,MagickPixelPacket *pixel)
 272 %
 273 %  A description of each parameter follows:
 274 %
 275 %    o resample_filter: the resample filter.
 276 %
 277 %    o method: the pixel clor interpolation method.
 278 %
 279 %    o x,y: A double representing the current (x,y) position of the pixel.
 280 %
 281 %    o pixel: return the interpolated pixel here.
 282 %
 283 */
 284
 285 static inline double MagickMax(const double x,const double y)
 286 {
 287   if (x > y)
 288     return(x);
 289   return(y);
 290 }
 291
 292 static void BicubicInterpolate(const MagickPixelPacket *pixels,const double dx,
 293   MagickPixelPacket *pixel)
 294 {
 295   MagickRealType
 296     dx2,
 297     p,
 298     q,
 299     r,
 300     s;
 301
 302   dx2=dx*dx;
 303   p=(pixels[3].red-pixels[2].red)-(pixels[0].red-pixels[1].red);
 304   q=(pixels[0].red-pixels[1].red)-p;
 305   r=pixels[2].red-pixels[0].red;
 306   s=pixels[1].red;
 307   pixel->red=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
 308   p=(pixels[3].green-pixels[2].green)-(pixels[0].green-pixels[1].green);
 309   q=(pixels[0].green-pixels[1].green)-p;
 310   r=pixels[2].green-pixels[0].green;
 311   s=pixels[1].green;
 312   pixel->green=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
 313   p=(pixels[3].blue-pixels[2].blue)-(pixels[0].blue-pixels[1].blue);
 314   q=(pixels[0].blue-pixels[1].blue)-p;
 315   r=pixels[2].blue-pixels[0].blue;
 316   s=pixels[1].blue;
 317   pixel->blue=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
 318   p=(pixels[3].opacity-pixels[2].opacity)-(pixels[0].opacity-pixels[1].opacity);
 319   q=(pixels[0].opacity-pixels[1].opacity)-p;
 320   r=pixels[2].opacity-pixels[0].opacity;
 321   s=pixels[1].opacity;
 322   pixel->opacity=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
 323   if (pixel->colorspace == CMYKColorspace)
 324     {
 325       p=(pixels[3].index-pixels[2].index)-(pixels[0].index-pixels[1].index);
 326       q=(pixels[0].index-pixels[1].index)-p;
 327       r=pixels[2].index-pixels[0].index;
 328       s=pixels[1].index;
 329       pixel->index=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
 330     }
 331 }
 332
 333 static inline MagickRealType CubicWeightingFunction(const MagickRealType x)
 334 {
 335   MagickRealType
 336     alpha,
 337     gamma;
 338
 339   alpha=MagickMax(x+2.0,0.0);
 340   gamma=1.0*alpha*alpha*alpha;
 341   alpha=MagickMax(x+1.0,0.0);
 342   gamma-=4.0*alpha*alpha*alpha;
 343   alpha=MagickMax(x+0.0,0.0);
 344   gamma+=6.0*alpha*alpha*alpha;
 345   alpha=MagickMax(x-1.0,0.0);
 346   gamma-=4.0*alpha*alpha*alpha;
 347   return(gamma/6.0);
 348 }
 349
 350 static inline double MeshInterpolate(const PointInfo *delta,const double p,
 351   const double x,const double y)
 352 {
 353   return(delta->x*x+delta->y*y+(1.0-delta->x-delta->y)*p);
 354 }
 355
 356 static inline ssize_t NearestNeighbor(MagickRealType x)
 357 {
 358   if (x >= 0.0)
 359     return((ssize_t) (x+0.5));
 360   return((ssize_t) (x-0.5));
 361 }
 362
 363 static MagickBooleanType InterpolateResampleFilter(
 364   ResampleFilter *resample_filter,const InterpolatePixelMethod method,
 365   const double x,const double y,MagickPixelPacket *pixel)
 366 {
 367   MagickBooleanType
 368     status;
 369
 370   register const IndexPacket
 371     *indexes;
 372
 373   register const PixelPacket
 374     *p;
 375
 376   register ssize_t
 377     i;
 378
 379   assert(resample_filter != (ResampleFilter *) NULL);
 380   assert(resample_filter->signature == MagickSignature);
 381   status=MagickTrue;
 382   switch (method)
 383   {
 384     case AverageInterpolatePixel:
 385     {
 386       MagickPixelPacket
 387         pixels[16];
 388
 389       MagickRealType
 390         alpha[16],
 391         gamma;
 392
 393       p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x)-1,
 394         (ssize_t) floor(y)-1,4,4,resample_filter->exception);
 395       if (p == (const PixelPacket *) NULL)
 396         {
 397           status=MagickFalse;
 398           break;
 399         }
 400       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 401       for (i=0; i < 16L; i++)
 402       {
 403         GetMagickPixelPacket(resample_filter->image,pixels+i);
 404         SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
 405         alpha[i]=1.0;
 406         if (resample_filter->image->matte != MagickFalse)
 407           {
 408             alpha[i]=QuantumScale*((MagickRealType) GetAlphaPixelComponent(p));
 409             pixels[i].red*=alpha[i];
 410             pixels[i].green*=alpha[i];
 411             pixels[i].blue*=alpha[i];
 412             if (resample_filter->image->colorspace == CMYKColorspace)
 413               pixels[i].index*=alpha[i];
 414           }
 415         gamma=alpha[i];
 416         gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 417         pixel->red+=gamma*0.0625*pixels[i].red;
 418         pixel->green+=gamma*0.0625*pixels[i].green;
 419         pixel->blue+=gamma*0.0625*pixels[i].blue;
 420         pixel->opacity+=0.0625*pixels[i].opacity;
 421         if (resample_filter->image->colorspace == CMYKColorspace)
 422           pixel->index+=gamma*0.0625*pixels[i].index;
 423         p++;
 424       }
 425       break;
 426     }
 427     case BicubicInterpolatePixel:
 428     {
 429       MagickPixelPacket
 430         pixels[16],
 431         u[4];
 432
 433       MagickRealType
 434         alpha[16];
 435
 436       PointInfo
 437         delta;
 438
 439       p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x)-1,
 440         (ssize_t) floor(y)-1,4,4,resample_filter->exception);
 441       if (p == (const PixelPacket *) NULL)
 442         {
 443           status=MagickFalse;
 444           break;
 445         }
 446       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 447       for (i=0; i < 16L; i++)
 448       {
 449         GetMagickPixelPacket(resample_filter->image,pixels+i);
 450         SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
 451         alpha[i]=1.0;
 452         if (resample_filter->image->matte != MagickFalse)
 453           {
 454             alpha[i]=QuantumScale*((MagickRealType) GetAlphaPixelComponent(p));
 455             pixels[i].red*=alpha[i];
 456             pixels[i].green*=alpha[i];
 457             pixels[i].blue*=alpha[i];
 458             if (resample_filter->image->colorspace == CMYKColorspace)
 459               pixels[i].index*=alpha[i];
 460           }
 461         p++;
 462       }
 463       delta.x=x-floor(x);
 464       for (i=0; i < 4L; i++)
 465         BicubicInterpolate(pixels+4*i,delta.x,u+i);
 466       delta.y=y-floor(y);
 467       BicubicInterpolate(u,delta.y,pixel);
 468       break;
 469     }
 470     case BilinearInterpolatePixel:
 471     default:
 472     {
 473       MagickPixelPacket
 474         pixels[4];
 475
 476       MagickRealType
 477         alpha[4],
 478         gamma;
 479
 480       PointInfo
 481         delta,
 482         epsilon;
 483
 484       p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x),
 485         (ssize_t) floor(y),2,2,resample_filter->exception);
 486       if (p == (const PixelPacket *) NULL)
 487         {
 488           status=MagickFalse;
 489           break;
 490         }
 491       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 492       for (i=0; i < 4L; i++)
 493       {
 494         pixels[i].red=(MagickRealType) p[i].red;
 495         pixels[i].green=(MagickRealType) p[i].green;
 496         pixels[i].blue=(MagickRealType) p[i].blue;
 497         pixels[i].opacity=(MagickRealType) p[i].opacity;
 498         alpha[i]=1.0;
 499       }
 500       if (resample_filter->image->matte != MagickFalse)
 501         for (i=0; i < 4L; i++)
 502         {
 503           alpha[i]=QuantumScale*((MagickRealType) QuantumRange-p[i].opacity);
 504           pixels[i].red*=alpha[i];
 505           pixels[i].green*=alpha[i];
 506           pixels[i].blue*=alpha[i];
 507         }
 508       if (indexes != (IndexPacket *) NULL)
 509         for (i=0; i < 4L; i++)
 510         {
 511           pixels[i].index=(MagickRealType) indexes[i];
 512           if (resample_filter->image->colorspace == CMYKColorspace)
 513             pixels[i].index*=alpha[i];
 514         }
 515       delta.x=x-floor(x);
 516       delta.y=y-floor(y);
 517       epsilon.x=1.0-delta.x;
 518       epsilon.y=1.0-delta.y;
 519       gamma=((epsilon.y*(epsilon.x*alpha[0]+delta.x*alpha[1])+delta.y*
 520         (epsilon.x*alpha[2]+delta.x*alpha[3])));
 521       gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 522       pixel->red=gamma*(epsilon.y*(epsilon.x*pixels[0].red+delta.x*
 523         pixels[1].red)+delta.y*(epsilon.x*pixels[2].red+delta.x*pixels[3].red));
 524       pixel->green=gamma*(epsilon.y*(epsilon.x*pixels[0].green+delta.x*
 525         pixels[1].green)+delta.y*(epsilon.x*pixels[2].green+delta.x*
 526         pixels[3].green));
 527       pixel->blue=gamma*(epsilon.y*(epsilon.x*pixels[0].blue+delta.x*
 528         pixels[1].blue)+delta.y*(epsilon.x*pixels[2].blue+delta.x*
 529         pixels[3].blue));
 530       pixel->opacity=(epsilon.y*(epsilon.x*pixels[0].opacity+delta.x*
 531         pixels[1].opacity)+delta.y*(epsilon.x*pixels[2].opacity+delta.x*
 532         pixels[3].opacity));
 533       if (resample_filter->image->colorspace == CMYKColorspace)
 534         pixel->index=gamma*(epsilon.y*(epsilon.x*pixels[0].index+delta.x*
 535           pixels[1].index)+delta.y*(epsilon.x*pixels[2].index+delta.x*
 536           pixels[3].index));
 537       break;
 538     }
 539     case FilterInterpolatePixel:
 540     {
 541       CacheView
 542         *filter_view;
 543
 544       Image
 545         *excerpt_image,
 546         *filter_image;
 547
 548       MagickPixelPacket
 549         pixels[1];
 550
 551       RectangleInfo
 552         geometry;
 553
 554       geometry.width=4L;
 555       geometry.height=4L;
 556       geometry.x=(ssize_t) floor(x)-1L;
 557       geometry.y=(ssize_t) floor(y)-1L;
 558       excerpt_image=ExcerptImage(resample_filter->image,&geometry,
 559         resample_filter->exception);
 560       if (excerpt_image == (Image *) NULL)
 561         {
 562           status=MagickFalse;
 563           break;
 564         }
 565       filter_image=ResizeImage(excerpt_image,1,1,resample_filter->image->filter,
 566         resample_filter->image->blur,resample_filter->exception);
 567       excerpt_image=DestroyImage(excerpt_image);
 568       if (filter_image == (Image *) NULL)
 569         break;
 570       filter_view=AcquireCacheView(filter_image);
 571       p=GetCacheViewVirtualPixels(filter_view,0,0,1,1,
 572         resample_filter->exception);
 573       if (p != (const PixelPacket *) NULL)
 574         {
 575           indexes=GetVirtualIndexQueue(filter_image);
 576           GetMagickPixelPacket(resample_filter->image,pixels);
 577           SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
 578         }
 579       filter_view=DestroyCacheView(filter_view);
 580       filter_image=DestroyImage(filter_image);
 581       break;
 582     }
 583     case IntegerInterpolatePixel:
 584     {
 585       MagickPixelPacket
 586         pixels[1];
 587
 588       p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x),
 589         (ssize_t) floor(y),1,1,resample_filter->exception);
 590       if (p == (const PixelPacket *) NULL)
 591         {
 592           status=MagickFalse;
 593           break;
 594         }
 595       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 596       GetMagickPixelPacket(resample_filter->image,pixels);
 597       SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
 598       break;
 599     }
 600     case MeshInterpolatePixel:
 601     {
 602       MagickPixelPacket
 603         pixels[4];
 604
 605       MagickRealType
 606         alpha[4],
 607         gamma;
 608
 609       PointInfo
 610         delta,
 611         luminance;
 612
 613       p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x),
 614         (ssize_t) floor(y),2,2,resample_filter->exception);
 615       if (p == (const PixelPacket *) NULL)
 616         {
 617           status=MagickFalse;
 618           break;
 619         }
 620       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 621       for (i=0; i < 4L; i++)
 622       {
 623         GetMagickPixelPacket(resample_filter->image,pixels+i);
 624         SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
 625         alpha[i]=1.0;
 626         if (resample_filter->image->matte != MagickFalse)
 627           {
 628             alpha[i]=QuantumScale*((MagickRealType) GetAlphaPixelComponent(p));
 629             pixels[i].red*=alpha[i];
 630             pixels[i].green*=alpha[i];
 631             pixels[i].blue*=alpha[i];
 632             if (resample_filter->image->colorspace == CMYKColorspace)
 633               pixels[i].index*=alpha[i];
 634           }
 635         p++;
 636       }
 637       delta.x=x-floor(x);
 638       delta.y=y-floor(y);
 639       luminance.x=MagickPixelLuminance(pixels+0)-MagickPixelLuminance(pixels+3);
 640       luminance.y=MagickPixelLuminance(pixels+1)-MagickPixelLuminance(pixels+2);
 641       if (fabs(luminance.x) < fabs(luminance.y))
 642         {
 643           /*
 644             Diagonal 0-3 NW-SE.
 645           */
 646           if (delta.x <= delta.y)
 647             {
 648               /*
 649                 Bottom-left triangle  (pixel:2, diagonal: 0-3).
 650               */
 651               delta.y=1.0-delta.y;
 652               gamma=MeshInterpolate(&delta,alpha[2],alpha[3],alpha[0]);
 653               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 654               pixel->red=gamma*MeshInterpolate(&delta,pixels[2].red,
 655                 pixels[3].red,pixels[0].red);
 656               pixel->green=gamma*MeshInterpolate(&delta,pixels[2].green,
 657                 pixels[3].green,pixels[0].green);
 658               pixel->blue=gamma*MeshInterpolate(&delta,pixels[2].blue,
 659                 pixels[3].blue,pixels[0].blue);
 660               pixel->opacity=gamma*MeshInterpolate(&delta,pixels[2].opacity,
 661                 pixels[3].opacity,pixels[0].opacity);
 662               if (resample_filter->image->colorspace == CMYKColorspace)
 663                 pixel->index=gamma*MeshInterpolate(&delta,pixels[2].index,
 664                   pixels[3].index,pixels[0].index);
 665             }
 666           else
 667             {
 668               /*
 669                 Top-right triangle (pixel:1, diagonal: 0-3).
 670               */
 671               delta.x=1.0-delta.x;
 672               gamma=MeshInterpolate(&delta,alpha[1],alpha[0],alpha[3]);
 673               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 674               pixel->red=gamma*MeshInterpolate(&delta,pixels[1].red,
 675                 pixels[0].red,pixels[3].red);
 676               pixel->green=gamma*MeshInterpolate(&delta,pixels[1].green,
 677                 pixels[0].green,pixels[3].green);
 678               pixel->blue=gamma*MeshInterpolate(&delta,pixels[1].blue,
 679                 pixels[0].blue,pixels[3].blue);
 680               pixel->opacity=gamma*MeshInterpolate(&delta,pixels[1].opacity,
 681                 pixels[0].opacity,pixels[3].opacity);
 682               if (resample_filter->image->colorspace == CMYKColorspace)
 683                 pixel->index=gamma*MeshInterpolate(&delta,pixels[1].index,
 684                   pixels[0].index,pixels[3].index);
 685             }
 686         }
 687       else
 688         {
 689           /*
 690             Diagonal 1-2 NE-SW.
 691           */
 692           if (delta.x <= (1.0-delta.y))
 693             {
 694               /*
 695                 Top-left triangle (pixel 0, diagonal: 1-2).
 696               */
 697               gamma=MeshInterpolate(&delta,alpha[0],alpha[1],alpha[2]);
 698               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 699               pixel->red=gamma*MeshInterpolate(&delta,pixels[0].red,
 700                 pixels[1].red,pixels[2].red);
 701               pixel->green=gamma*MeshInterpolate(&delta,pixels[0].green,
 702                 pixels[1].green,pixels[2].green);
 703               pixel->blue=gamma*MeshInterpolate(&delta,pixels[0].blue,
 704                 pixels[1].blue,pixels[2].blue);
 705               pixel->opacity=gamma*MeshInterpolate(&delta,pixels[0].opacity,
 706                 pixels[1].opacity,pixels[2].opacity);
 707               if (resample_filter->image->colorspace == CMYKColorspace)
 708                 pixel->index=gamma*MeshInterpolate(&delta,pixels[0].index,
 709                   pixels[1].index,pixels[2].index);
 710             }
 711           else
 712             {
 713               /*
 714                 Bottom-right triangle (pixel: 3, diagonal: 1-2).
 715               */
 716               delta.x=1.0-delta.x;
 717               delta.y=1.0-delta.y;
 718               gamma=MeshInterpolate(&delta,alpha[3],alpha[2],alpha[1]);
 719               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 720               pixel->red=gamma*MeshInterpolate(&delta,pixels[3].red,
 721                 pixels[2].red,pixels[1].red);
 722               pixel->green=gamma*MeshInterpolate(&delta,pixels[3].green,
 723                 pixels[2].green,pixels[1].green);
 724               pixel->blue=gamma*MeshInterpolate(&delta,pixels[3].blue,
 725                 pixels[2].blue,pixels[1].blue);
 726               pixel->opacity=gamma*MeshInterpolate(&delta,pixels[3].opacity,
 727                 pixels[2].opacity,pixels[1].opacity);
 728               if (resample_filter->image->colorspace == CMYKColorspace)
 729                 pixel->index=gamma*MeshInterpolate(&delta,pixels[3].index,
 730                   pixels[2].index,pixels[1].index);
 731             }
 732         }
 733       break;
 734     }
 735     case NearestNeighborInterpolatePixel:
 736     {
 737       MagickPixelPacket
 738         pixels[1];
 739
 740       p=GetCacheViewVirtualPixels(resample_filter->view,NearestNeighbor(x),
 741         NearestNeighbor(y),1,1,resample_filter->exception);
 742       if (p == (const PixelPacket *) NULL)
 743         {
 744           status=MagickFalse;
 745           break;
 746         }
 747       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 748       GetMagickPixelPacket(resample_filter->image,pixels);
 749       SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
 750       break;
 751     }
 752     case SplineInterpolatePixel:
 753     {
 754       MagickPixelPacket
 755         pixels[16];
 756
 757       MagickRealType
 758         alpha[16],
 759         dx,
 760         dy,
 761         gamma;
 762
 763       PointInfo
 764         delta;
 765
 766       ssize_t
 767         j,
 768         n;
 769
 770       p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x)-1,
 771         (ssize_t) floor(y)-1,4,4,resample_filter->exception);
 772       if (p == (const PixelPacket *) NULL)
 773         {
 774           status=MagickFalse;
 775           break;
 776         }
 777       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 778       n=0;
 779       delta.x=x-floor(x);
 780       delta.y=y-floor(y);
 781       for (i=(-1); i < 3L; i++)
 782       {
 783         dy=CubicWeightingFunction((MagickRealType) i-delta.y);
 784         for (j=(-1); j < 3L; j++)
 785         {
 786           GetMagickPixelPacket(resample_filter->image,pixels+n);
 787           SetMagickPixelPacket(resample_filter->image,p,indexes+n,pixels+n);
 788           alpha[n]=1.0;
 789           if (resample_filter->image->matte != MagickFalse)
 790             {
 791               alpha[n]=QuantumScale*((MagickRealType)
 792                 GetAlphaPixelComponent(p));
 793               pixels[n].red*=alpha[n];
 794               pixels[n].green*=alpha[n];
 795               pixels[n].blue*=alpha[n];
 796               if (resample_filter->image->colorspace == CMYKColorspace)
 797                 pixels[n].index*=alpha[n];
 798             }
 799           dx=CubicWeightingFunction(delta.x-(MagickRealType) j);
 800           gamma=alpha[n];
 801           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 802           pixel->red+=gamma*dx*dy*pixels[n].red;
 803           pixel->green+=gamma*dx*dy*pixels[n].green;
 804           pixel->blue+=gamma*dx*dy*pixels[n].blue;
 805           if (resample_filter->image->matte != MagickFalse)
 806             pixel->opacity+=dx*dy*pixels[n].opacity;
 807           if (resample_filter->image->colorspace == CMYKColorspace)
 808             pixel->index+=gamma*dx*dy*pixels[n].index;
 809           n++;
 810           p++;
 811         }
 812       }
 813       break;
 814     }
 815   }
 816   return(status);
 817 }
 818 \f
 819 /*
 820 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 821 %                                                                             %
 822 %                                                                             %
 823 %                                                                             %
 824 %   R e s a m p l e P i x e l C o l o r                                       %
 825 %                                                                             %
 826 %                                                                             %
 827 %                                                                             %
 828 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 829 %
 830 %  ResamplePixelColor() samples the pixel values surrounding the location
 831 %  given using an elliptical weighted average, at the scale previously
 832 %  calculated, and in the most efficent manner possible for the
 833 %  VirtualPixelMethod setting.
 834 %
 835 %  The format of the ResamplePixelColor method is:
 836 %
 837 %     MagickBooleanType ResamplePixelColor(ResampleFilter *resample_filter,
 838 %       const double u0,const double v0,MagickPixelPacket *pixel)
 839 %
 840 %  A description of each parameter follows:
 841 %
 842 %    o resample_filter: the resample filter.
 843 %
 844 %    o u0,v0: A double representing the center of the area to resample,
 845 %        The distortion transformed transformed x,y coordinate.
 846 %
 847 %    o pixel: the resampled pixel is returned here.
 848 %
 849 */
 850 MagickExport MagickBooleanType ResamplePixelColor(
 851   ResampleFilter *resample_filter,const double u0,const double v0,
 852   MagickPixelPacket *pixel)
 853 {
 854   MagickBooleanType
 855     status;
 856
 857   ssize_t u,v, uw,v1,v2, hit;
 858   double u1;
 859   double U,V,Q,DQ,DDQ;
 860   double divisor_c,divisor_m;
 861   register double weight;
 862   register const PixelPacket *pixels;
 863   register const IndexPacket *indexes;
 864   assert(resample_filter != (ResampleFilter *) NULL);
 865   assert(resample_filter->signature == MagickSignature);
 866
 867   status=MagickTrue;
 868   GetMagickPixelPacket(resample_filter->image,pixel);
 869   if ( resample_filter->do_interpolate ) {
 870     status=InterpolateResampleFilter(resample_filter,
 871       resample_filter->interpolate,u0,v0,pixel);
 872     return(status);
 873   }
 874
 875   /*
 876     Does resample area Miss the image?
 877     And is that area a simple solid color - then return that color
 878   */
 879   hit = 0;
 880   switch ( resample_filter->virtual_pixel ) {
 881     case BackgroundVirtualPixelMethod:
 882     case ConstantVirtualPixelMethod:
 883     case TransparentVirtualPixelMethod:
 884     case BlackVirtualPixelMethod:
 885     case GrayVirtualPixelMethod:
 886     case WhiteVirtualPixelMethod:
 887     case MaskVirtualPixelMethod:
 888       if ( resample_filter->limit_reached
 889            || u0 + resample_filter->Ulimit < 0.0
 890            || u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 891            || v0 + resample_filter->Vlimit < 0.0
 892            || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows
 893            )
 894         hit++;
 895       break;
 896
 897     case UndefinedVirtualPixelMethod:
 898     case EdgeVirtualPixelMethod:
 899       if (    ( u0 + resample_filter->Ulimit < 0.0 && v0 + resample_filter->Vlimit < 0.0 )
 900            || ( u0 + resample_filter->Ulimit < 0.0
 901                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows )
 902            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 903                 && v0 + resample_filter->Vlimit < 0.0 )
 904            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 905                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows )
 906            )
 907         hit++;
 908       break;
 909     case HorizontalTileVirtualPixelMethod:
 910       if (    v0 + resample_filter->Vlimit < 0.0
 911            || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows
 912            )
 913         hit++;  /* outside the horizontally tiled images. */
 914       break;
 915     case VerticalTileVirtualPixelMethod:
 916       if (    u0 + resample_filter->Ulimit < 0.0
 917            || u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 918            )
 919         hit++;  /* outside the vertically tiled images. */
 920       break;
 921     case DitherVirtualPixelMethod:
 922       if (    ( u0 + resample_filter->Ulimit < -32.0 && v0 + resample_filter->Vlimit < -32.0 )
 923            || ( u0 + resample_filter->Ulimit < -32.0
 924                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+32.0 )
 925            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+32.0
 926                 && v0 + resample_filter->Vlimit < -32.0 )
 927            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+32.0
 928                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+32.0 )
 929            )
 930         hit++;
 931       break;
 932     case TileVirtualPixelMethod:
 933     case MirrorVirtualPixelMethod:
 934     case RandomVirtualPixelMethod:
 935     case HorizontalTileEdgeVirtualPixelMethod:
 936     case VerticalTileEdgeVirtualPixelMethod:
 937     case CheckerTileVirtualPixelMethod:
 938       /* resampling of area is always needed - no VP limits */
 939       break;
 940   }
 941   if ( hit ) {
 942     /* whole area is a solid color -- just return that color */
 943     status=InterpolateResampleFilter(resample_filter,IntegerInterpolatePixel,
 944       u0,v0,pixel);
 945     return(status);
 946   }
 947
 948   /*
 949     Scaling limits reached, return an 'averaged' result.
 950   */
 951   if ( resample_filter->limit_reached ) {
 952     switch ( resample_filter->virtual_pixel ) {
 953       /*  This is always handled by the above, so no need.
 954         case BackgroundVirtualPixelMethod:
 955         case ConstantVirtualPixelMethod:
 956         case TransparentVirtualPixelMethod:
 957         case GrayVirtualPixelMethod,
 958         case WhiteVirtualPixelMethod
 959         case MaskVirtualPixelMethod:
 960       */
 961       case UndefinedVirtualPixelMethod:
 962       case EdgeVirtualPixelMethod:
 963       case DitherVirtualPixelMethod:
 964       case HorizontalTileEdgeVirtualPixelMethod:
 965       case VerticalTileEdgeVirtualPixelMethod:
 966         /* We need an average edge pixel, from the correct edge!
 967            How should I calculate an average edge color?
 968            Just returning an averaged neighbourhood,
 969            works well in general, but falls down for TileEdge methods.
 970            This needs to be done properly!!!!!!
 971         */
 972         status=InterpolateResampleFilter(resample_filter,
 973           AverageInterpolatePixel,u0,v0,pixel);
 974         break;
 975       case HorizontalTileVirtualPixelMethod:
 976       case VerticalTileVirtualPixelMethod:
 977         /* just return the background pixel - Is there more direct way? */
 978         status=InterpolateResampleFilter(resample_filter,
 979            IntegerInterpolatePixel,(double)-1,(double)-1,pixel);
 980         break;
 981       case TileVirtualPixelMethod:
 982       case MirrorVirtualPixelMethod:
 983       case RandomVirtualPixelMethod:
 984       case CheckerTileVirtualPixelMethod:
 985       default:
 986         /* generate a average color of the WHOLE image */
 987         if ( resample_filter->average_defined == MagickFalse ) {
 988           Image
 989             *average_image;
 990
 991           CacheView
 992             *average_view;
 993
 994           GetMagickPixelPacket(resample_filter->image,
 995                 (MagickPixelPacket *)&(resample_filter->average_pixel));
 996           resample_filter->average_defined = MagickTrue;
 997
 998           /* Try to get an averaged pixel color of whole image */
 999           average_image=ResizeImage(resample_filter->image,1,1,BoxFilter,1.0,
1000            resample_filter->exception);
1001           if (average_image == (Image *) NULL)
1002             {
1003               *pixel=resample_filter->average_pixel; /* FAILED */
1004               break;
1005             }
1006           average_view=AcquireCacheView(average_image);
1007           pixels=(PixelPacket *)GetCacheViewVirtualPixels(average_view,0,0,1,1,
1008             resample_filter->exception);
1009           if (pixels == (const PixelPacket *) NULL) {
1010             average_view=DestroyCacheView(average_view);
1011             average_image=DestroyImage(average_image);
1012             *pixel=resample_filter->average_pixel; /* FAILED */
1013             break;
1014           }
1015           indexes=(IndexPacket *) GetCacheViewAuthenticIndexQueue(average_view);
1016           SetMagickPixelPacket(resample_filter->image,pixels,indexes,
1017             &(resample_filter->average_pixel));
1018           average_view=DestroyCacheView(average_view);
1019           average_image=DestroyImage(average_image);
1020 #if 0
1021           /* CheckerTile should average the image with background color */
1022           //if ( resample_filter->virtual_pixel == CheckerTileVirtualPixelMethod ) {
1023 #if 0
1024             resample_filter->average_pixel.red =
1025                       ( resample_filter->average_pixel.red +
1026                           resample_filter->image->background_color.red ) /2;
1027             resample_filter->average_pixel.green =
1028                       ( resample_filter->average_pixel.green +
1029                           resample_filter->image->background_color.green ) /2;
1030             resample_filter->average_pixel.blue =
1031                       ( resample_filter->average_pixel.blue +
1032                           resample_filter->image->background_color.blue ) /2;
1033             resample_filter->average_pixel.matte =
1034                       ( resample_filter->average_pixel.matte +
1035                           resample_filter->image->background_color.matte ) /2;
1036             resample_filter->average_pixel.black =
1037                       ( resample_filter->average_pixel.black +
1038                           resample_filter->image->background_color.black ) /2;
1039 #else
1040             resample_filter->average_pixel =
1041                           resample_filter->image->background_color;
1042 #endif
1043           }
1044 #endif
1045         }
1046         *pixel=resample_filter->average_pixel;
1047         break;
1048     }
1049     return(status);
1050   }
1051
1052   /*
1053     Initialize weighted average data collection
1054   */
1055   hit = 0;
1056   divisor_c = 0.0;
1057   divisor_m = 0.0;
1058   pixel->red = pixel->green = pixel->blue = 0.0;
1059   if (resample_filter->image->matte != MagickFalse) pixel->opacity = 0.0;
1060   if (resample_filter->image->colorspace == CMYKColorspace) pixel->index = 0.0;
1061
1062   /*
1063     Determine the parellelogram bounding box fitted to the ellipse
1064     centered at u0,v0.  This area is bounding by the lines...
1065         v = +/- sqrt(A)
1066         u = -By/2A  +/- sqrt(F/A)
1067     Which has been pre-calculated above.
1068   */
1069   v1 = (ssize_t)(v0 - resample_filter->Vlimit);               /* range of scan lines */
1070   v2 = (ssize_t)(v0 + resample_filter->Vlimit + 1);
1071
1072   u1 = u0 + (v1-v0)*resample_filter->slope - resample_filter->Uwidth; /* start of scanline for v=v1 */
1073   uw = (ssize_t)(2*resample_filter->Uwidth)+1;       /* width of parallelogram */
1074
1075   /*
1076     Do weighted resampling of all pixels,  within the scaled ellipse,
1077     bound by a Parellelogram fitted to the ellipse.
1078   */
1079   DDQ = 2*resample_filter->A;
1080   for( v=v1; v<=v2;  v++, u1+=resample_filter->slope ) {
1081     u = (ssize_t)u1;       /* first pixel in scanline  ( floor(u1) ) */
1082     U = (double)u-u0;   /* location of that pixel, relative to u0,v0 */
1083     V = (double)v-v0;
1084
1085     /* Q = ellipse quotent ( if Q<F then pixel is inside ellipse) */
1086     Q = U*(resample_filter->A*U + resample_filter->B*V) + resample_filter->C*V*V;
1087     DQ = resample_filter->A*(2.0*U+1) + resample_filter->B*V;
1088
1089     /* get the scanline of pixels for this v */
1090     pixels=GetCacheViewVirtualPixels(resample_filter->view,u,v,(size_t) uw,
1091       1,resample_filter->exception);
1092     if (pixels == (const PixelPacket *) NULL)
1093       return(MagickFalse);
1094     indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
1095
1096     /* count up the weighted pixel colors */
1097     for( u=0; u<uw; u++ ) {
1098       /* Note that the ellipse has been pre-scaled so F = WLUT_WIDTH */
1099       if ( Q < (double)WLUT_WIDTH ) {
1100         weight = resample_filter->filter_lut[(int)Q];
1101
1102         pixel->opacity  += weight*pixels->opacity;
1103         divisor_m += weight;
1104
1105         if (resample_filter->image->matte != MagickFalse)
1106           weight *= QuantumScale*((MagickRealType)(QuantumRange-pixels->opacity));
1107         pixel->red   += weight*pixels->red;
1108         pixel->green += weight*pixels->green;
1109         pixel->blue  += weight*pixels->blue;
1110         if (resample_filter->image->colorspace == CMYKColorspace)
1111           pixel->index += weight*(*indexes);
1112         divisor_c += weight;
1113
1114         hit++;
1115       }
1116       pixels++;
1117       indexes++;
1118       Q += DQ;
1119       DQ += DDQ;
1120     }
1121   }
1122
1123   /*
1124     Result sanity check -- this should NOT happen
1125   */
1126   if ( hit == 0 ) {  /* should be 4 */
1127     /* not enough pixels in resampling, resort to direct interpolation */
1128 #if 0
1129     pixel->opacity = pixel->red = pixel->green = pixel->blue = 0;
1130     pixel->red = QuantumRange; /* show pixels for which EWA fails */
1131 #else
1132     status=InterpolateResampleFilter(resample_filter,
1133       resample_filter->interpolate,u0,v0,pixel);
1134 #endif
1135     return status;
1136   }
1137
1138   /*
1139     Finialize results of resampling
1140   */
1141   divisor_m = 1.0/divisor_m;
1142   pixel->opacity = (MagickRealType) ClampToQuantum(divisor_m*pixel->opacity);
1143   divisor_c = 1.0/divisor_c;
1144   pixel->red   = (MagickRealType) ClampToQuantum(divisor_c*pixel->red);
1145   pixel->green = (MagickRealType) ClampToQuantum(divisor_c*pixel->green);
1146   pixel->blue  = (MagickRealType) ClampToQuantum(divisor_c*pixel->blue);
1147   if (resample_filter->image->colorspace == CMYKColorspace)
1148     pixel->index = (MagickRealType) ClampToQuantum(divisor_c*pixel->index);
1149   return(MagickTrue);
1150 }
1151 \f
1152 /*
1153 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1154 %                                                                             %
1155 %                                                                             %
1156 %                                                                             %
1157 %   S c a l e R e s a m p l e F i l t e r                                     %
1158 %                                                                             %
1159 %                                                                             %
1160 %                                                                             %
1161 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1162 %
1163 %  ScaleResampleFilter() does all the calculations needed to resample an image
1164 %  at a specific scale, defined by two scaling vectors.  This not using
1165 %  a orthogonal scaling, but two distorted scaling vectors, to allow the
1166 %  generation of a angled ellipse.
1167 %
1168 %  As only two deritive scaling vectors are used the center of the ellipse
1169 %  must be the center of the lookup.  That is any curvature that the
1170 %  distortion may produce is discounted.
1171 %
1172 %  The input vectors are produced by either finding the derivitives of the
1173 %  distortion function, or the partial derivitives from a distortion mapping.
1174 %  They do not need to be the orthogonal dx,dy scaling vectors, but can be
1175 %  calculated from other derivatives.  For example you could use  dr,da/r
1176 %  polar coordinate vector scaling vectors
1177 %
1178 %  If   u,v =  DistortEquation(x,y)
1179 %  Then the scaling vectors dx,dy (in u,v space) are the derivitives...
1180 %      du/dx, dv/dx     and    du/dy, dv/dy
1181 %  If the scaling is only othogonally aligned then...
1182 %      dv/dx = 0   and   du/dy  =  0
1183 %  Producing an othogonally alligned ellipse for the area to be resampled.
1184 %
1185 %  Note that scaling vectors are different to argument order.  Argument order
1186 %  is the general order the deritives are extracted from the distortion
1187 %  equations, EG: U(x,y), V(x,y).  Caution is advised if you are trying to
1188 %  define the ellipse directly from scaling vectors.
1189 %
1190 %  The format of the ScaleResampleFilter method is:
1191 %
1192 %     void ScaleResampleFilter(const ResampleFilter *resample_filter,
1193 %       const double dux,const double duy,const double dvx,const double dvy)
1194 %
1195 %  A description of each parameter follows:
1196 %
1197 %    o resample_filter: the resampling resample_filterrmation defining the
1198 %      image being resampled
1199 %
1200 %    o dux,duy,dvx,dvy:
1201 %         The partial derivitives or scaling vectors for resampling.
1202 %           dx = du/dx, dv/dx    and  dy = du/dy, dv/dy
1203 %
1204 %         The values are used to define the size and angle of the
1205 %         elliptical resampling area, centered on the lookup point.
1206 %
1207 */
1208 MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
1209   const double dux,const double duy,const double dvx,const double dvy)
1210 {
1211   double A,B,C,F;
1212
1213   assert(resample_filter != (ResampleFilter *) NULL);
1214   assert(resample_filter->signature == MagickSignature);
1215
1216   resample_filter->limit_reached = MagickFalse;
1217   resample_filter->do_interpolate = MagickFalse;
1218
1219   /* A 'point' filter forces use of interpolation instead of area sampling */
1220   if ( resample_filter->filter == PointFilter ) {
1221     resample_filter->do_interpolate = MagickTrue;
1222     return;
1223   }
1224
1225   /* Find Ellipse Coefficents such that
1226         A*u^2 + B*u*v + C*v^2 = F
1227      With u,v relative to point around which we are resampling.
1228      And the given scaling dx,dy vectors in u,v space
1229          du/dx,dv/dx   and  du/dy,dv/dy
1230   */
1231 #if 0
1232   /* Direct conversion of derivatives into elliptical coefficients
1233     No scaling will result in F == 1.0 and a unit circle.
1234     However if the ellipse becomes very small (magnification) or
1235     very long and thin, the ellipse may miss all source pixels!
1236   */
1237   A = dvx*dvx+dvy*dvy;
1238   B = -2.0*(dux*dvx+duy*dvy);
1239   C = dux*dux+duy*duy;
1240   F = dux*dvy+duy*dvx;
1241   F *= F; /* square it */
1242 #define F_UNITY 1.0
1243
1244 #else
1245   /*
1246     This Paul Heckbert's recomended "Higher Quality EWA" formula, from page
1247     60 in his thesis, which adds a unit circle to the elliptical area so as
1248     to do both Reconstruction and Prefiltering of the pixels in the
1249     resampling.  It also means it is always likely to have at least 4 pixels
1250     within the area of the ellipse, for weighted averaging.  No scaling will
1251     result with F == 4.0 and a circle of radius 2.0, and F smaller than this
1252     means magnification is being used.
1253   */
1254   A = dvx*dvx+dvy*dvy+1;
1255   B = -2.0*(dux*dvx+duy*dvy);
1256   C = dux*dux+duy*duy+1;
1257   F = A*C - B*B/4;
1258 #define F_UNITY 4.0
1259 #endif
1260
1261 /* DEBUGGING OUTPUT */
1262 #if 0
1263   fprintf(stderr, "dux=%lf; dvx=%lf;   duy=%lf; dvy%lf;\n",
1264        dux, dvx, duy, dvy);
1265   fprintf(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
1266 #endif
1267
1268 #if 0
1269   /* Figure out the Ellipses Major and Minor Axis, and other info.
1270      This information currently not needed at this time, but may be
1271      needed later for better limit determination.
1272
1273      It is also good to have as a record for future debugging
1274   */
1275   { double alpha, beta, gamma, Major, Minor;
1276     double Eccentricity, Ellipse_Area, Ellipse_angle;
1277     double max_horizontal_cross_section, max_vertical_cross_section;
1278     double parellelogram_slope;
1279
1280     alpha = A+C;
1281     beta  = A-C;
1282     gamma = sqrt(beta*beta + B*B );
1283
1284     if ( alpha - gamma <= MagickEpsilon )
1285       Major = MagickHuge;
1286     else
1287       Major = sqrt(2*F/(alpha - gamma));
1288     Minor = sqrt(2*F/(alpha + gamma));
1289
1290     fprintf(stderr, "\tMajor=%lf; Minor=%lf\n",
1291          Major, Minor );
1292
1293     /* other information about ellipse include... */
1294     Eccentricity = Major/Minor;
1295     Ellipse_Area = MagickPI*Major*Minor;
1296     Ellipse_angle =  atan2(B, A-C);
1297
1298     fprintf(stderr, "\tAngle=%lf Area=%lf\n",
1299          RadiansToDegrees(Ellipse_angle), Ellipse_Area );
1300
1301     /* Ellipse Orthogonal Bounds */
1302     /* max_horizontal_orthogonal = sqrt(4*A*F/(4*A*C-B*B)) */
1303     /* max_vertical_orthogonal   = sqrt(4*C*F/(4*A*C-B*B)) */
1304
1305     /* After optimization using the improved ellipse */
1306     /* Note how F cancels out divisor and the 4, leaving only A and C */
1307     max_horizontal_orthogonal = sqrt(A);
1308     max_vertical_orthogonal   = sqrt(C);
1309
1310     /* Parallelogram Bounds (axis intercepts) */
1311     max_horizontal_cross_section = sqrt(F/A);
1312     max_vertical_cross_section   = sqrt(F/C);
1313     parellelogram_slope = -B/(2*A);
1314   }
1315 #endif
1316
1317   /* Is default elliptical area, too small? Then the image is being magnified.
1318      Switch to using a orthogonal interpolation method,
1319      unless 'filter' interpolation was specifically requested.
1320   */
1321   if ( F <= F_UNITY && resample_filter->interpolate != FilterInterpolatePixel ) {
1322     resample_filter->do_interpolate = MagickTrue;
1323     return;
1324   }
1325
1326   /* If F is impossibly large, we may as well not bother doing any
1327      form of resampling, as you risk an near infinite resampled area.
1328      In this case some alturnative means of pixel sampling, such as
1329      the average of the whole image is needed to get a reasonable result.
1330   */
1331   if ( F > MagickHuge ) {
1332     resample_filter->limit_reached = MagickTrue;
1333     return;
1334   }
1335
1336   /* Othogonal bounds of the Improved Ellipse */
1337   resample_filter->Ulimit = sqrt(C)+1.0;   /* Horizontal Orthogonal Limit */
1338   resample_filter->Vlimit = sqrt(A)+1.0;   /* Vertical Orthogonal Limit */
1339
1340   /* Horizontally aligned Parallelogram fitted to ellipse */
1341   resample_filter->Uwidth = sqrt(F/A)+1.0;   /* Parallelogram Width */
1342   resample_filter->slope = -B/(2*A);        /* Slope of the parallelogram */
1343
1344
1345   /* Check the absolute area of the Parallogram involved...
1346    * This limit needs more work, as it gets too slow for
1347    * larger images involved with tiled views of the horizon. */
1348   if ( resample_filter->Uwidth * resample_filter->Vlimit
1349          > 5.0*resample_filter->image_area ) {
1350     resample_filter->limit_reached = MagickTrue;
1351     return;
1352   }
1353
1354   /* Scale ellipse formula to directly index the Filter Lookup Table */
1355   { register double scale;
1356     scale = (double)WLUT_WIDTH/(F*4);  /* 4 = cache_support^2 */
1357     resample_filter->A = A*scale;
1358     resample_filter->B = B*scale;
1359     resample_filter->C = C*scale;
1360     /* ..ple_filter->F = WLUT_WIDTH; -- hardcoded */
1361   }
1362 }
1363 \f
1364 /*
1365 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1366 %                                                                             %
1367 %                                                                             %
1368 %                                                                             %
1369 %   S e t R e s a m p l e F i l t e r                                         %
1370 %                                                                             %
1371 %                                                                             %
1372 %                                                                             %
1373 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1374 %
1375 %  SetResampleFilter() set the resampling filter lookup table based on a
1376 %  specific filter.  Note that the filter is used as a radial filter not as a
1377 %  two pass othogonally aligned resampling filter.
1378 %
1379 %  The default Filter, is Gaussian, which is the standard filter used by the
1380 %  original paper on the Elliptical Weighted Everage Algorithm. However other
1381 %  filters can also be used.
1382 %
1383 %  The format of the SetResampleFilter method is:
1384 %
1385 %    void SetResampleFilter(ResampleFilter *resample_filter,
1386 %      const FilterTypes filter,const double blur)
1387 %
1388 %  A description of each parameter follows:
1389 %
1390 %    o resample_filter: resampling resample_filterrmation structure
1391 %
1392 %    o filter: the resize filter for elliptical weighting LUT
1393 %
1394 %    o blur: filter blur factor (radial scaling) for elliptical weighting LUT
1395 %
1396 */
1397 MagickExport void SetResampleFilter(ResampleFilter *resample_filter,
1398   const FilterTypes filter,const double blur)
1399 {
1400   register int
1401      Q;
1402
1403   double
1404      r_scale;
1405
1406   ResizeFilter
1407      *resize_filter;
1408
1409   assert(resample_filter != (ResampleFilter *) NULL);
1410   assert(resample_filter->signature == MagickSignature);
1411
1412   resample_filter->filter = filter;
1413
1414   /* Scale radius so it equals 1.0, at edge of ellipse when a
1415      default blurring factor of 1.0 is used.
1416
1417      Note that these filters are being used as a radial filter, not as
1418      an othoginally alligned filter. How this effects results is still
1419      being worked out.
1420
1421      Future: Direct use of the resize filters in "resize.c" to set the lookup
1422      table, based on the filters working support window.
1423   */
1424   r_scale = sqrt(1.0/(double)WLUT_WIDTH);
1425   r_scale *= 2; /* for 2 pixel radius of Improved Elliptical Formula */
1426
1427   switch ( resample_filter->filter ) {
1428   case PointFilter:
1429     /* This equivelent to turning off the EWA algroithm.
1430        Only Interpolated lookup will be used.  */
1431     break;
1432   case UndefinedFilter:
1433     resample_filter->filter = GaussianFilter;
1434     /* FALL-THRU */
1435   default:
1436     /*
1437       Fill the LUT with a 1D resize filter function
1438       But make the Sinc/Bessel tapered window 2.0
1439       I also normalize the result so the filter is 1.0
1440     */
1441     resize_filter = AcquireResizeFilter(resample_filter->image,
1442          resample_filter->filter,blur,MagickTrue,resample_filter->exception);
1443     if (resize_filter != (ResizeFilter *) NULL) {
1444 #if 0
1445       /* At this time the filter support is completely ignored! */
1446       resample_filter->support = GetResizeFilterSupport(resize_filter);
1447       resample_filter->support /= blur; /* taken into account above */
1448       resample_filter->support *= resample_filter->support;
1449       resample_filter->support *= (double)WLUT_WIDTH/4;
1450       if ( resample_filter->support >= (double)WLUT_WIDTH )
1451         resample_filter->support = (double)WLUT_WIDTH; /* not used yet */
1452 #endif
1453       for(Q=0; Q<WLUT_WIDTH; Q++)
1454         resample_filter->filter_lut[Q] = (double)
1455              GetResizeFilterWeight(resize_filter,sqrt((double)Q)*r_scale);
1456       resize_filter = DestroyResizeFilter(resize_filter);
1457       break;
1458     }
1459     else {
1460       (void) ThrowMagickException(resample_filter->exception,GetMagickModule(),
1461            ModuleError, "UnableToSetFilteringValue",
1462            "Fall back to default EWA gaussian filter");
1463     }
1464 #if 0
1465   This is old code that is only
1466   /*
1467     Create Normal Gaussian 2D Filter Weighted Lookup Table.
1468     A normal EWA guassual lookup would use   exp(Q*ALPHA)
1469     where  Q = distance squared from 0.0 (center) to 1.0 (edge)
1470     and    ALPHA = -4.0*ln(2.0)  ==>  -2.77258872223978123767
1471     However the table is of length 1024, and equates to a radius of 2px
1472     thus needs to be scaled by  ALPHA*4/1024 and any blur factor squared
1473
1474     This is now known to be wrong -- very wrong!
1475   */
1476   /*r_scale = -2.77258872223978123767*4/WLUT_WIDTH/blur/blur;*/
1477   r_scale = -2.77258872223978123767/(WLUT_WIDTH*blur*blur);
1478   for(Q=0; Q<WLUT_WIDTH; Q++)
1479     resample_filter->filter_lut[Q] = exp((double)Q*r_scale);
1480   resample_filter->support = WLUT_WIDTH;
1481   break;
1482 #endif
1483   }
1484 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1485   /* if( GetOpenMPThreadId() == 0 ) */
1486 #endif
1487     if (GetImageArtifact(resample_filter->image,"resample:verbose")
1488           != (const char *) NULL)
1489       {
1490         /* Debug output of the filter weighting LUT
1491           Gnuplot the LUT with hoizontal adjusted to 'r' using...
1492             plot [0:2][-.2:1] "lut.dat" using (sqrt($0/1024)*2):1 with lines
1493           The filter values is normalized for comparision
1494         */
1495         printf("#\n");
1496         printf("# Resampling Filter LUT (%d values)\n", WLUT_WIDTH);
1497         printf("#\n");
1498         printf("# Note: values in table are using a squared radius lookup.\n");
1499         printf("# And the whole table represents the filters support.\n");
1500         printf("#\n");
1501         for(Q=0; Q<WLUT_WIDTH; Q++)
1502           printf("%8.*g %.*g\n",
1503                GetMagickPrecision(),sqrt((double)Q)*r_scale,
1504                GetMagickPrecision(),resample_filter->filter_lut[Q] );
1505       }
1506   return;
1507 }
1508 \f
1509 /*
1510 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1511 %                                                                             %
1512 %                                                                             %
1513 %                                                                             %
1514 %   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       %
1515 %                                                                             %
1516 %                                                                             %
1517 %                                                                             %
1518 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1519 %
1520 %  SetResampleFilterInterpolateMethod() changes the interpolation method
1521 %  associated with the specified resample filter.
1522 %
1523 %  The format of the SetResampleFilterInterpolateMethod method is:
1524 %
1525 %      MagickBooleanType SetResampleFilterInterpolateMethod(
1526 %        ResampleFilter *resample_filter,const InterpolateMethod method)
1527 %
1528 %  A description of each parameter follows:
1529 %
1530 %    o resample_filter: the resample filter.
1531 %
1532 %    o method: the interpolation method.
1533 %
1534 */
1535 MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
1536   ResampleFilter *resample_filter,const InterpolatePixelMethod method)
1537 {
1538   assert(resample_filter != (ResampleFilter *) NULL);
1539   assert(resample_filter->signature == MagickSignature);
1540   assert(resample_filter->image != (Image *) NULL);
1541
1542   if (resample_filter->debug != MagickFalse)
1543     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1544       resample_filter->image->filename);
1545
1546   resample_filter->interpolate=method;
1547
1548   return(MagickTrue);
1549 }
1550 \f
1551 /*
1552 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1553 %                                                                             %
1554 %                                                                             %
1555 %                                                                             %
1556 %   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     %
1557 %                                                                             %
1558 %                                                                             %
1559 %                                                                             %
1560 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1561 %
1562 %  SetResampleFilterVirtualPixelMethod() changes the virtual pixel method
1563 %  associated with the specified resample filter.
1564 %
1565 %  The format of the SetResampleFilterVirtualPixelMethod method is:
1566 %
1567 %      MagickBooleanType SetResampleFilterVirtualPixelMethod(
1568 %        ResampleFilter *resample_filter,const VirtualPixelMethod method)
1569 %
1570 %  A description of each parameter follows:
1571 %
1572 %    o resample_filter: the resample filter.
1573 %
1574 %    o method: the virtual pixel method.
1575 %
1576 */
1577 MagickExport MagickBooleanType SetResampleFilterVirtualPixelMethod(
1578   ResampleFilter *resample_filter,const VirtualPixelMethod method)
1579 {
1580   assert(resample_filter != (ResampleFilter *) NULL);
1581   assert(resample_filter->signature == MagickSignature);
1582   assert(resample_filter->image != (Image *) NULL);
1583   if (resample_filter->debug != MagickFalse)
1584     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1585       resample_filter->image->filename);
1586   resample_filter->virtual_pixel=method;
1587   if (method != UndefinedVirtualPixelMethod)
1588     (void) SetCacheViewVirtualPixelMethod(resample_filter->view,method);
1589   return(MagickTrue);
1590 }