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   ssize_t
  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   size_t
 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 < (ssize_t) image->rows; y++) {
 147 %        for (x=0; x < (ssize_t) 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=(ssize_t) (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 ssize_t NearestNeighbor(MagickRealType x)
 356 {
 357   if (x >= 0.0)
 358     return((ssize_t) (x+0.5));
 359   return((ssize_t) (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 ssize_t
 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,(ssize_t) floor(x)-1,
 393         (ssize_t) 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) GetAlphaPixelComponent(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,(ssize_t) floor(x)-1,
 439         (ssize_t) 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) GetAlphaPixelComponent(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,(ssize_t) floor(x),
 484         (ssize_t) 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       CacheView
 541         *filter_view;
 542
 543       Image
 544         *excerpt_image,
 545         *filter_image;
 546
 547       MagickPixelPacket
 548         pixels[1];
 549
 550       RectangleInfo
 551         geometry;
 552
 553       geometry.width=4L;
 554       geometry.height=4L;
 555       geometry.x=(ssize_t) floor(x)-1L;
 556       geometry.y=(ssize_t) 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,(ssize_t) floor(x),
 588         (ssize_t) 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,(ssize_t) floor(x),
 613         (ssize_t) 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) GetAlphaPixelComponent(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       MagickPixelPacket
 754         pixels[16];
 755
 756       MagickRealType
 757         alpha[16],
 758         dx,
 759         dy,
 760         gamma;
 761
 762       PointInfo
 763         delta;
 764
 765       ssize_t
 766         j,
 767         n;
 768
 769       p=GetCacheViewVirtualPixels(resample_filter->view,(ssize_t) floor(x)-1,
 770         (ssize_t) 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)
 791                 GetAlphaPixelComponent(p));
 792               pixels[n].red*=alpha[n];
 793               pixels[n].green*=alpha[n];
 794               pixels[n].blue*=alpha[n];
 795               if (resample_filter->image->colorspace == CMYKColorspace)
 796                 pixels[n].index*=alpha[n];
 797             }
 798           dx=CubicWeightingFunction(delta.x-(MagickRealType) j);
 799           gamma=alpha[n];
 800           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
 801           pixel->red+=gamma*dx*dy*pixels[n].red;
 802           pixel->green+=gamma*dx*dy*pixels[n].green;
 803           pixel->blue+=gamma*dx*dy*pixels[n].blue;
 804           if (resample_filter->image->matte != MagickFalse)
 805             pixel->opacity+=dx*dy*pixels[n].opacity;
 806           if (resample_filter->image->colorspace == CMYKColorspace)
 807             pixel->index+=gamma*dx*dy*pixels[n].index;
 808           n++;
 809           p++;
 810         }
 811       }
 812       break;
 813     }
 814   }
 815   return(status);
 816 }
 817 \f
 818 /*
 819 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 820 %                                                                             %
 821 %                                                                             %
 822 %                                                                             %
 823 %   R e s a m p l e P i x e l C o l o r                                       %
 824 %                                                                             %
 825 %                                                                             %
 826 %                                                                             %
 827 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 828 %
 829 %  ResamplePixelColor() samples the pixel values surrounding the location
 830 %  given using an elliptical weighted average, at the scale previously
 831 %  calculated, and in the most efficent manner possible for the
 832 %  VirtualPixelMethod setting.
 833 %
 834 %  The format of the ResamplePixelColor method is:
 835 %
 836 %     MagickBooleanType ResamplePixelColor(ResampleFilter *resample_filter,
 837 %       const double u0,const double v0,MagickPixelPacket *pixel)
 838 %
 839 %  A description of each parameter follows:
 840 %
 841 %    o resample_filter: the resample filter.
 842 %
 843 %    o u0,v0: A double representing the center of the area to resample,
 844 %        The distortion transformed transformed x,y coordinate.
 845 %
 846 %    o pixel: the resampled pixel is returned here.
 847 %
 848 */
 849 MagickExport MagickBooleanType ResamplePixelColor(
 850   ResampleFilter *resample_filter,const double u0,const double v0,
 851   MagickPixelPacket *pixel)
 852 {
 853   MagickBooleanType
 854     status;
 855
 856   ssize_t u,v, uw,v1,v2, hit;
 857   double u1;
 858   double U,V,Q,DQ,DDQ;
 859   double divisor_c,divisor_m;
 860   register double weight;
 861   register const PixelPacket *pixels;
 862   register const IndexPacket *indexes;
 863   assert(resample_filter != (ResampleFilter *) NULL);
 864   assert(resample_filter->signature == MagickSignature);
 865
 866   status=MagickTrue;
 867   GetMagickPixelPacket(resample_filter->image,pixel);
 868   if ( resample_filter->do_interpolate ) {
 869     status=InterpolateResampleFilter(resample_filter,
 870       resample_filter->interpolate,u0,v0,pixel);
 871     return(status);
 872   }
 873
 874   /*
 875     Does resample area Miss the image?
 876     And is that area a simple solid color - then return that color
 877   */
 878   hit = 0;
 879   switch ( resample_filter->virtual_pixel ) {
 880     case BackgroundVirtualPixelMethod:
 881     case ConstantVirtualPixelMethod:
 882     case TransparentVirtualPixelMethod:
 883     case BlackVirtualPixelMethod:
 884     case GrayVirtualPixelMethod:
 885     case WhiteVirtualPixelMethod:
 886     case MaskVirtualPixelMethod:
 887       if ( resample_filter->limit_reached
 888            || u0 + resample_filter->sqrtC < 0.0
 889            || u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
 890            || v0 + resample_filter->sqrtA < 0.0
 891            || v0 - resample_filter->sqrtA > (double) resample_filter->image->rows
 892            )
 893         hit++;
 894       break;
 895
 896     case UndefinedVirtualPixelMethod:
 897     case EdgeVirtualPixelMethod:
 898       if (    ( u0 + resample_filter->sqrtC < 0.0 && v0 + resample_filter->sqrtA < 0.0 )
 899            || ( u0 + resample_filter->sqrtC < 0.0
 900                 && v0 - resample_filter->sqrtA > (double) resample_filter->image->rows )
 901            || ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
 902                 && v0 + resample_filter->sqrtA < 0.0 )
 903            || ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
 904                 && v0 - resample_filter->sqrtA > (double) resample_filter->image->rows )
 905            )
 906         hit++;
 907       break;
 908     case HorizontalTileVirtualPixelMethod:
 909       if (    v0 + resample_filter->sqrtA < 0.0
 910            || v0 - resample_filter->sqrtA > (double) resample_filter->image->rows
 911            )
 912         hit++;  /* outside the horizontally tiled images. */
 913       break;
 914     case VerticalTileVirtualPixelMethod:
 915       if (    u0 + resample_filter->sqrtC < 0.0
 916            || u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
 917            )
 918         hit++;  /* outside the vertically tiled images. */
 919       break;
 920     case DitherVirtualPixelMethod:
 921       if (    ( u0 + resample_filter->sqrtC < -32.0 && v0 + resample_filter->sqrtA < -32.0 )
 922            || ( u0 + resample_filter->sqrtC < -32.0
 923                 && v0 - resample_filter->sqrtA > (double) resample_filter->image->rows+32.0 )
 924            || ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns+32.0
 925                 && v0 + resample_filter->sqrtA < -32.0 )
 926            || ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns+32.0
 927                 && v0 - resample_filter->sqrtA > (double) resample_filter->image->rows+32.0 )
 928            )
 929         hit++;
 930       break;
 931     case TileVirtualPixelMethod:
 932     case MirrorVirtualPixelMethod:
 933     case RandomVirtualPixelMethod:
 934     case HorizontalTileEdgeVirtualPixelMethod:
 935     case VerticalTileEdgeVirtualPixelMethod:
 936     case CheckerTileVirtualPixelMethod:
 937       /* resampling of area is always needed - no VP limits */
 938       break;
 939   }
 940   if ( hit ) {
 941     /* whole area is a solid color -- just return that color */
 942     status=InterpolateResampleFilter(resample_filter,IntegerInterpolatePixel,
 943       u0,v0,pixel);
 944     return(status);
 945   }
 946
 947   /*
 948     Scaling limits reached, return an 'averaged' result.
 949   */
 950   if ( resample_filter->limit_reached ) {
 951     switch ( resample_filter->virtual_pixel ) {
 952       /*  This is always handled by the above, so no need.
 953         case BackgroundVirtualPixelMethod:
 954         case ConstantVirtualPixelMethod:
 955         case TransparentVirtualPixelMethod:
 956         case GrayVirtualPixelMethod,
 957         case WhiteVirtualPixelMethod
 958         case MaskVirtualPixelMethod:
 959       */
 960       case UndefinedVirtualPixelMethod:
 961       case EdgeVirtualPixelMethod:
 962       case DitherVirtualPixelMethod:
 963       case HorizontalTileEdgeVirtualPixelMethod:
 964       case VerticalTileEdgeVirtualPixelMethod:
 965         /* We need an average edge pixel, for the right edge!
 966            How should I calculate an average edge color?
 967            Just returning an averaged neighbourhood,
 968            works well in general, but falls down for TileEdge methods.
 969            This needs to be done properly!!!!!!
 970         */
 971         status=InterpolateResampleFilter(resample_filter,
 972           AverageInterpolatePixel,u0,v0,pixel);
 973         break;
 974       case HorizontalTileVirtualPixelMethod:
 975       case VerticalTileVirtualPixelMethod:
 976         /* just return the background pixel - Is there more direct way? */
 977         status=InterpolateResampleFilter(resample_filter,
 978            IntegerInterpolatePixel,(double)-1,(double)-1,pixel);
 979         break;
 980       case TileVirtualPixelMethod:
 981       case MirrorVirtualPixelMethod:
 982       case RandomVirtualPixelMethod:
 983       case CheckerTileVirtualPixelMethod:
 984       default:
 985         /* generate a average color of the WHOLE image */
 986         if ( resample_filter->average_defined == MagickFalse ) {
 987           Image
 988             *average_image;
 989
 990           CacheView
 991             *average_view;
 992
 993           GetMagickPixelPacket(resample_filter->image,
 994                 (MagickPixelPacket *)&(resample_filter->average_pixel));
 995           resample_filter->average_defined = MagickTrue;
 996
 997           /* Try to get an averaged pixel color of whole image */
 998           average_image=ResizeImage(resample_filter->image,1,1,BoxFilter,1.0,
 999            resample_filter->exception);
1000           if (average_image == (Image *) NULL)
1001             {
1002               *pixel=resample_filter->average_pixel; /* FAILED */
1003               break;
1004             }
1005           average_view=AcquireCacheView(average_image);
1006           pixels=(PixelPacket *)GetCacheViewVirtualPixels(average_view,0,0,1,1,
1007             resample_filter->exception);
1008           if (pixels == (const PixelPacket *) NULL) {
1009             average_view=DestroyCacheView(average_view);
1010             average_image=DestroyImage(average_image);
1011             *pixel=resample_filter->average_pixel; /* FAILED */
1012             break;
1013           }
1014           indexes=(IndexPacket *) GetCacheViewAuthenticIndexQueue(average_view);
1015           SetMagickPixelPacket(resample_filter->image,pixels,indexes,
1016             &(resample_filter->average_pixel));
1017           average_view=DestroyCacheView(average_view);
1018           average_image=DestroyImage(average_image);
1019 #if 0
1020           /* CheckerTile should average the image with background color */
1021           //if ( resample_filter->virtual_pixel == CheckerTileVirtualPixelMethod ) {
1022 #if 0
1023             resample_filter->average_pixel.red =
1024                       ( resample_filter->average_pixel.red +
1025                           resample_filter->image->background_color.red ) /2;
1026             resample_filter->average_pixel.green =
1027                       ( resample_filter->average_pixel.green +
1028                           resample_filter->image->background_color.green ) /2;
1029             resample_filter->average_pixel.blue =
1030                       ( resample_filter->average_pixel.blue +
1031                           resample_filter->image->background_color.blue ) /2;
1032             resample_filter->average_pixel.matte =
1033                       ( resample_filter->average_pixel.matte +
1034                           resample_filter->image->background_color.matte ) /2;
1035             resample_filter->average_pixel.black =
1036                       ( resample_filter->average_pixel.black +
1037                           resample_filter->image->background_color.black ) /2;
1038 #else
1039             resample_filter->average_pixel =
1040                           resample_filter->image->background_color;
1041 #endif
1042           }
1043 #endif
1044         }
1045         *pixel=resample_filter->average_pixel;
1046         break;
1047     }
1048     return(status);
1049   }
1050
1051   /*
1052     Initialize weighted average data collection
1053   */
1054   hit = 0;
1055   divisor_c = 0.0;
1056   divisor_m = 0.0;
1057   pixel->red = pixel->green = pixel->blue = 0.0;
1058   if (resample_filter->image->matte != MagickFalse) pixel->opacity = 0.0;
1059   if (resample_filter->image->colorspace == CMYKColorspace) pixel->index = 0.0;
1060
1061   /*
1062     Determine the parellelogram bounding box fitted to the ellipse
1063     centered at u0,v0.  This area is bounding by the lines...
1064         v = +/- sqrt(A)
1065         u = -By/2A  +/- sqrt(F/A)
1066     Which has been pre-calculated above.
1067   */
1068   v1 = (ssize_t)(v0 - resample_filter->sqrtA);               /* range of scan lines */
1069   v2 = (ssize_t)(v0 + resample_filter->sqrtA + 1);
1070
1071   u1 = u0 + (v1-v0)*resample_filter->slope - resample_filter->sqrtU; /* start of scanline for v=v1 */
1072   uw = (ssize_t)(2*resample_filter->sqrtU)+1;       /* width of parallelogram */
1073
1074   /*
1075     Do weighted resampling of all pixels,  within the scaled ellipse,
1076     bound by a Parellelogram fitted to the ellipse.
1077   */
1078   DDQ = 2*resample_filter->A;
1079   for( v=v1; v<=v2;  v++, u1+=resample_filter->slope ) {
1080     u = (ssize_t)u1;       /* first pixel in scanline  ( floor(u1) ) */
1081     U = (double)u-u0;   /* location of that pixel, relative to u0,v0 */
1082     V = (double)v-v0;
1083
1084     /* Q = ellipse quotent ( if Q<F then pixel is inside ellipse) */
1085     Q = U*(resample_filter->A*U + resample_filter->B*V) + resample_filter->C*V*V;
1086     DQ = resample_filter->A*(2.0*U+1) + resample_filter->B*V;
1087
1088     /* get the scanline of pixels for this v */
1089     pixels=GetCacheViewVirtualPixels(resample_filter->view,u,v,(size_t) uw,
1090       1,resample_filter->exception);
1091     if (pixels == (const PixelPacket *) NULL)
1092       return(MagickFalse);
1093     indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
1094
1095     /* count up the weighted pixel colors */
1096     for( u=0; u<uw; u++ ) {
1097       /* Note that the ellipse has been pre-scaled so F = WLUT_WIDTH */
1098       if ( Q < (double)WLUT_WIDTH ) {
1099         weight = resample_filter->filter_lut[(int)Q];
1100
1101         pixel->opacity  += weight*pixels->opacity;
1102         divisor_m += weight;
1103
1104         if (resample_filter->image->matte != MagickFalse)
1105           weight *= QuantumScale*((MagickRealType)(QuantumRange-pixels->opacity));
1106         pixel->red   += weight*pixels->red;
1107         pixel->green += weight*pixels->green;
1108         pixel->blue  += weight*pixels->blue;
1109         if (resample_filter->image->colorspace == CMYKColorspace)
1110           pixel->index += weight*(*indexes);
1111         divisor_c += weight;
1112
1113         hit++;
1114       }
1115       pixels++;
1116       indexes++;
1117       Q += DQ;
1118       DQ += DDQ;
1119     }
1120   }
1121
1122   /*
1123     Result sanity check -- this should NOT happen
1124   */
1125   if ( hit < 4 || divisor_c < 1.0 ) {
1126     /* not enough pixels in resampling, resort to direct interpolation */
1127     status=InterpolateResampleFilter(resample_filter,
1128       resample_filter->interpolate,u0,v0,pixel);
1129     return status;
1130   }
1131
1132   /*
1133     Finialize results of resampling
1134   */
1135   divisor_m = 1.0/divisor_m;
1136   pixel->opacity = (MagickRealType) ClampToQuantum(divisor_m*pixel->opacity);
1137   divisor_c = 1.0/divisor_c;
1138   pixel->red   = (MagickRealType) ClampToQuantum(divisor_c*pixel->red);
1139   pixel->green = (MagickRealType) ClampToQuantum(divisor_c*pixel->green);
1140   pixel->blue  = (MagickRealType) ClampToQuantum(divisor_c*pixel->blue);
1141   if (resample_filter->image->colorspace == CMYKColorspace)
1142     pixel->index = (MagickRealType) ClampToQuantum(divisor_c*pixel->index);
1143   return(MagickTrue);
1144 }
1145 \f
1146 /*
1147 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1148 %                                                                             %
1149 %                                                                             %
1150 %                                                                             %
1151 %   S c a l e R e s a m p l e F i l t e r                                     %
1152 %                                                                             %
1153 %                                                                             %
1154 %                                                                             %
1155 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1156 %
1157 %  ScaleResampleFilter() does all the calculations needed to resample an image
1158 %  at a specific scale, defined by two scaling vectors.  This not using
1159 %  a orthogonal scaling, but two distorted scaling vectors, to allow the
1160 %  generation of a angled ellipse.
1161 %
1162 %  As only two deritive scaling vectors are used the center of the ellipse
1163 %  must be the center of the lookup.  That is any curvature that the
1164 %  distortion may produce is discounted.
1165 %
1166 %  The input vectors are produced by either finding the derivitives of the
1167 %  distortion function, or the partial derivitives from a distortion mapping.
1168 %  They do not need to be the orthogonal dx,dy scaling vectors, but can be
1169 %  calculated from other derivatives.  For example you could use  dr,da/r
1170 %  polar coordinate vector scaling vectors
1171 %
1172 %  If   u,v =  DistortEquation(x,y)
1173 %  Then the scaling vectors dx,dy (in u,v space) are the derivitives...
1174 %      du/dx, dv/dx     and    du/dy, dv/dy
1175 %  If the scaling is only othogonally aligned then...
1176 %      dv/dx = 0   and   du/dy  =  0
1177 %  Producing an othogonally alligned ellipse for the area to be resampled.
1178 %
1179 %  Note that scaling vectors are different to argument order.  Argument order
1180 %  is the general order the deritives are extracted from the distortion
1181 %  equations, EG: U(x,y), V(x,y).  Caution is advised if you are trying to
1182 %  define the ellipse directly from scaling vectors.
1183 %
1184 %  The format of the ScaleResampleFilter method is:
1185 %
1186 %     void ScaleResampleFilter(const ResampleFilter *resample_filter,
1187 %       const double dux,const double duy,const double dvx,const double dvy)
1188 %
1189 %  A description of each parameter follows:
1190 %
1191 %    o resample_filter: the resampling resample_filterrmation defining the
1192 %      image being resampled
1193 %
1194 %    o dux,duy,dvx,dvy:
1195 %         The partial derivitives or scaling vectors for resampling.
1196 %           dx = du/dx, dv/dx    and  dy = du/dy, dv/dy
1197 %
1198 %         The values are used to define the size and angle of the
1199 %         elliptical resampling area, centered on the lookup point.
1200 %
1201 */
1202 MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
1203   const double dux,const double duy,const double dvx,const double dvy)
1204 {
1205   double A,B,C,F, area;
1206
1207   assert(resample_filter != (ResampleFilter *) NULL);
1208   assert(resample_filter->signature == MagickSignature);
1209
1210   resample_filter->limit_reached = MagickFalse;
1211   resample_filter->do_interpolate = MagickFalse;
1212
1213   /* A 'point' filter forces use of interpolation instead of area sampling */
1214   if ( resample_filter->filter == PointFilter ) {
1215     resample_filter->do_interpolate = MagickTrue;
1216     return;
1217   }
1218
1219   /* Find Ellipse Coefficents such that
1220         A*u^2 + B*u*v + C*v^2 = F
1221      With u,v relative to point around which we are resampling.
1222      And the given scaling dx,dy vectors in u,v space
1223          du/dx,dv/dx   and  du/dy,dv/dy
1224   */
1225 #if 0
1226   /* Direct conversions of derivatives to elliptical coefficients
1227      No scaling will result in F == 1.0 and a unit circle.
1228   */
1229   A = dvx*dvx+dvy*dvy;
1230   B = (dux*dvx+duy*dvy)*-2.0;
1231   C = dux*dux+duy*duy;
1232   F = dux*dvy+duy*dvx;
1233   F *= F;
1234 #define F_UNITY 1.0
1235 #else
1236   /* This Paul Heckbert's recomended "Higher Quality EWA" formula, from page
1237      60 in his thesis, which adds a unit circle to the elliptical area so are
1238      to do both Reconstruction and Prefiltering of the pixels in the
1239      resampling.  It also means it is likely to have at least 4 pixels within
1240      the area of the ellipse, for weighted averaging.
1241      No scaling will result if F == 4.0 and a circle of radius 2.0
1242   */
1243   A = dvx*dvx+dvy*dvy+1;
1244   B = (dux*dvx+duy*dvy)*-2.0;
1245   C = dux*dux+duy*duy+1;
1246   F = A*C - B*B/4;
1247 #define F_UNITY 4.0
1248 #endif
1249
1250 /* DEBUGGING OUTPUT */
1251 #if 0
1252   fprintf(stderr, "dux=%lf; dvx=%lf;   duy=%lf; dvy%lf;\n",
1253        dux, dvx, duy, dvy);
1254   fprintf(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
1255 #endif
1256
1257 #if 0
1258   /* Figure out the Ellipses Major and Minor Axis, and other info.
1259      This information currently not needed at this time, but may be
1260      needed later for better limit determination.
1261   */
1262   { double alpha, beta, gamma, Major, Minor;
1263     double Eccentricity, Ellipse_Area, Ellipse_angle;
1264     double max_horizontal_cross_section, max_vertical_cross_section;
1265     alpha = A+C;
1266     beta  = A-C;
1267     gamma = sqrt(beta*beta + B*B );
1268
1269     if ( alpha - gamma <= MagickEpsilon )
1270       Major = MagickHuge;
1271     else
1272       Major = sqrt(2*F/(alpha - gamma));
1273     Minor = sqrt(2*F/(alpha + gamma));
1274
1275     fprintf(stderr, "\tMajor=%lf; Minor=%lf\n",
1276          Major, Minor );
1277
1278     /* other information about ellipse include... */
1279     Eccentricity = Major/Minor;
1280     Ellipse_Area = MagickPI*Major*Minor;
1281     Ellipse_angle =  atan2(B, A-C);
1282
1283     fprintf(stderr, "\tAngle=%lf Area=%lf\n",
1284          RadiansToDegrees(Ellipse_angle), Ellipse_Area );
1285
1286     /* Ellipse limits */
1287
1288     /* orthogonal rectangle - improved ellipse */
1289     max_horizontal_orthogonal = sqrt(A); /* = sqrt(4*A*F/(4*A*C-B*B)) */
1290     max_vertical_orthogonal   = sqrt(C); /* = sqrt(4*C*F/(4*A*C-B*B)) */
1291
1292     /* parallelogram bounds -- what we are using */
1293     max_horizontal_cross_section = sqrt(F/A);
1294     max_vertical_cross_section   = sqrt(F/C);
1295   }
1296 #endif
1297
1298   /* Is default elliptical area, too small? Image being magnified?
1299      Switch to doing pure 'point' interpolation of the pixel.
1300      That is turn off  EWA Resampling.
1301   */
1302   if ( F <= F_UNITY ) {
1303     resample_filter->do_interpolate = MagickTrue;
1304     return;
1305   }
1306
1307
1308   /* If F is impossibly large, we may as well not bother doing any
1309    * form of resampling, as you risk an infinite resampled area.
1310   */
1311   if ( F > MagickHuge ) {
1312     resample_filter->limit_reached = MagickTrue;
1313     return;
1314   }
1315
1316   /* Othogonal bounds of the ellipse */
1317   resample_filter->sqrtA = sqrt(A)+1.0;     /* Vertical Orthogonal Limit */
1318   resample_filter->sqrtC = sqrt(C)+1.0;     /* Horizontal Orthogonal Limit */
1319
1320   /* Horizontally aligned Parallelogram fitted to ellipse */
1321   resample_filter->sqrtU = sqrt(F/A)+1.0;   /* Parallelogram Width */
1322   resample_filter->slope = -B/(2*A);        /* Slope of the parallelogram */
1323
1324   /* The size of the area of the parallelogram we will be sampling */
1325   area = 4 * resample_filter->sqrtA * resample_filter->sqrtU;
1326
1327   /* Absolute limit on the area to be resampled
1328    * This limit needs more work, as it gets too slow for
1329    * larger images involved with tiled views of the horizon. */
1330   if ( area > 20.0*resample_filter->image_area ) {
1331     resample_filter->limit_reached = MagickTrue;
1332     return;
1333   }
1334
1335   /* Scale ellipse formula to directly fit the Filter Lookup Table */
1336   { register double scale;
1337     scale = (double)WLUT_WIDTH/F;
1338     resample_filter->A = A*scale;
1339     resample_filter->B = B*scale;
1340     resample_filter->C = C*scale;
1341     /* ..ple_filter->F = WLUT_WIDTH; -- hardcoded */
1342   }
1343 }
1344 \f
1345 /*
1346 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1347 %                                                                             %
1348 %                                                                             %
1349 %                                                                             %
1350 %   S e t R e s a m p l e F i l t e r                                         %
1351 %                                                                             %
1352 %                                                                             %
1353 %                                                                             %
1354 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1355 %
1356 %  SetResampleFilter() set the resampling filter lookup table based on a
1357 %  specific filter.  Note that the filter is used as a radial filter not as a
1358 %  two pass othogonally aligned resampling filter.
1359 %
1360 %  The default Filter, is Gaussian, which is the standard filter used by the
1361 %  original paper on the Elliptical Weighted Everage Algorithm. However other
1362 %  filters can also be used.
1363 %
1364 %  The format of the SetResampleFilter method is:
1365 %
1366 %    void SetResampleFilter(ResampleFilter *resample_filter,
1367 %      const FilterTypes filter,const double blur)
1368 %
1369 %  A description of each parameter follows:
1370 %
1371 %    o resample_filter: resampling resample_filterrmation structure
1372 %
1373 %    o filter: the resize filter for elliptical weighting LUT
1374 %
1375 %    o blur: filter blur factor (radial scaling) for elliptical weighting LUT
1376 %
1377 */
1378 MagickExport void SetResampleFilter(ResampleFilter *resample_filter,
1379   const FilterTypes filter,const double blur)
1380 {
1381   register int
1382      Q;
1383
1384   double
1385      r_scale;
1386
1387   ResizeFilter
1388      *resize_filter;
1389
1390   assert(resample_filter != (ResampleFilter *) NULL);
1391   assert(resample_filter->signature == MagickSignature);
1392
1393   resample_filter->filter = filter;
1394
1395   /* Scale radius so it equals 1.0, at edge of ellipse when a
1396      default blurring factor of 1.0 is used.
1397
1398      Note that these filters are being used as a radial filter, not as
1399      an othoginally alligned filter. How this effects results is still
1400      being worked out.
1401
1402      Future: Direct use of teh resize filters in "resize.c" to set the lookup
1403      table, based on the filters working support window.
1404   */
1405   r_scale = sqrt(1.0/(double)WLUT_WIDTH)/blur;
1406   r_scale *= 2; /* for 2 pixel radius of Improved Elliptical Formula */
1407
1408   switch ( filter ) {
1409   case PointFilter:
1410     /* This equivelent to turning off the EWA algroithm.
1411        Only Interpolated lookup will be used.  */
1412     break;
1413   default:
1414     /*
1415       Fill the LUT with a 1D resize filter function
1416       But make the Sinc/Bessel tapered window 2.0
1417       I also normalize the result so the filter is 1.0
1418     */
1419     resize_filter = AcquireResizeFilter(resample_filter->image,filter,
1420          (MagickRealType)1.0,MagickTrue,resample_filter->exception);
1421     if (resize_filter != (ResizeFilter *) NULL) {
1422       resample_filter->support = GetResizeFilterSupport(resize_filter);
1423       resample_filter->support /= blur; /* taken into account above */
1424       resample_filter->support *= resample_filter->support;
1425       resample_filter->support *= (double)WLUT_WIDTH/4;
1426       if ( resample_filter->support >= (double)WLUT_WIDTH )
1427            resample_filter->support = (double)WLUT_WIDTH;  /* hack */
1428       for(Q=0; Q<WLUT_WIDTH; Q++)
1429         if ( (double) Q < resample_filter->support )
1430           resample_filter->filter_lut[Q] = (double)
1431                GetResizeFilterWeight(resize_filter,sqrt((double)Q)*r_scale);
1432         else
1433           resample_filter->filter_lut[Q] = 0.0;
1434       resize_filter = DestroyResizeFilter(resize_filter);
1435       break;
1436     }
1437     else {
1438       (void) ThrowMagickException(resample_filter->exception,GetMagickModule(),
1439            ModuleError, "UnableToSetFilteringValue",
1440            "Fall back to default EWA gaussian filter");
1441     }
1442     /* FALLTHRU - on exception */
1443   /*case GaussianFilter:*/
1444   case UndefinedFilter:
1445     /*
1446       Create Normal Gaussian 2D Filter Weighted Lookup Table.
1447       A normal EWA guassual lookup would use   exp(Q*ALPHA)
1448       where  Q = distantce squared from 0.0 (center) to 1.0 (edge)
1449       and    ALPHA = -4.0*ln(2.0)  ==>  -2.77258872223978123767
1450       However the table is of length 1024, and equates to a radius of 2px
1451       thus needs to be scaled by  ALPHA*4/1024 and any blur factor squared
1452     */
1453     /*r_scale = -2.77258872223978123767*4/WLUT_WIDTH/blur/blur;*/
1454     r_scale = -2.77258872223978123767/WLUT_WIDTH/blur/blur;
1455     for(Q=0; Q<WLUT_WIDTH; Q++)
1456       resample_filter->filter_lut[Q] = exp((double)Q*r_scale);
1457     resample_filter->support = WLUT_WIDTH;
1458     break;
1459   }
1460   if (GetImageArtifact(resample_filter->image,"resample:verbose")
1461         != (const char *) NULL)
1462     /* Debug output of the filter weighting LUT
1463       Gnuplot the LUT with hoizontal adjusted to 'r' using...
1464         plot [0:2][-.2:1] "lut.dat" using (sqrt($0/1024)*2):1 with lines
1465       THe filter values is normalized for comparision
1466     */
1467     for(Q=0; Q<WLUT_WIDTH; Q++)
1468       printf("%lf\n", resample_filter->filter_lut[Q]
1469                         /resample_filter->filter_lut[0] );
1470   return;
1471 }
1472 \f
1473 /*
1474 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1475 %                                                                             %
1476 %                                                                             %
1477 %                                                                             %
1478 %   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       %
1479 %                                                                             %
1480 %                                                                             %
1481 %                                                                             %
1482 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1483 %
1484 %  SetResampleFilterInterpolateMethod() changes the interpolation method
1485 %  associated with the specified resample filter.
1486 %
1487 %  The format of the SetResampleFilterInterpolateMethod method is:
1488 %
1489 %      MagickBooleanType SetResampleFilterInterpolateMethod(
1490 %        ResampleFilter *resample_filter,const InterpolateMethod method)
1491 %
1492 %  A description of each parameter follows:
1493 %
1494 %    o resample_filter: the resample filter.
1495 %
1496 %    o method: the interpolation method.
1497 %
1498 */
1499 MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
1500   ResampleFilter *resample_filter,const InterpolatePixelMethod method)
1501 {
1502   assert(resample_filter != (ResampleFilter *) NULL);
1503   assert(resample_filter->signature == MagickSignature);
1504   assert(resample_filter->image != (Image *) NULL);
1505   if (resample_filter->debug != MagickFalse)
1506     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1507       resample_filter->image->filename);
1508   resample_filter->interpolate=method;
1509   return(MagickTrue);
1510 }
1511 \f
1512 /*
1513 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1514 %                                                                             %
1515 %                                                                             %
1516 %                                                                             %
1517 %   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     %
1518 %                                                                             %
1519 %                                                                             %
1520 %                                                                             %
1521 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1522 %
1523 %  SetResampleFilterVirtualPixelMethod() changes the virtual pixel method
1524 %  associated with the specified resample filter.
1525 %
1526 %  The format of the SetResampleFilterVirtualPixelMethod method is:
1527 %
1528 %      MagickBooleanType SetResampleFilterVirtualPixelMethod(
1529 %        ResampleFilter *resample_filter,const VirtualPixelMethod method)
1530 %
1531 %  A description of each parameter follows:
1532 %
1533 %    o resample_filter: the resample filter.
1534 %
1535 %    o method: the virtual pixel method.
1536 %
1537 */
1538 MagickExport MagickBooleanType SetResampleFilterVirtualPixelMethod(
1539   ResampleFilter *resample_filter,const VirtualPixelMethod method)
1540 {
1541   assert(resample_filter != (ResampleFilter *) NULL);
1542   assert(resample_filter->signature == MagickSignature);
1543   assert(resample_filter->image != (Image *) NULL);
1544   if (resample_filter->debug != MagickFalse)
1545     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1546       resample_filter->image->filename);
1547   resample_filter->virtual_pixel=method;
1548   if (method != UndefinedVirtualPixelMethod)
1549     (void) SetCacheViewVirtualPixelMethod(resample_filter->view,method);
1550   return(MagickTrue);
1551 }