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