granicus.if.org Git - imagemagick/blob - magick/resample.c

   1 /*
   2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   3 %                                                                             %
   4 %                                                                             %
   5 %                                                                             %
   6 %           RRRR    EEEEE   SSSSS   AAA   M   M  PPPP   L      EEEEE          %
   7 %           R   R   E       SS     A   A  MM MM  P   P  L      E              %
   8 %           RRRR    EEE      SSS   AAAAA  M M M  PPPP   L      EEE            %
   9 %           R R     E          SS  A   A  M   M  P      L      E              %
  10 %           R  R    EEEEE   SSSSS  A   A  M   M  P      LLLLL  EEEEE          %
  11 %                                                                             %
  12 %                                                                             %
  13 %                      MagickCore Pixel Resampling Methods                    %
  14 %                                                                             %
  15 %                              Software Design                                %
  16 %                                John Cristy                                  %
  17 %                              Anthony Thyssen                                %
  18 %                                August 2007                                  %
  19 %                                                                             %
  20 %                                                                             %
  21 %  Copyright 1999-2011 ImageMagick Studio LLC, a non-profit organization      %
  22 %  dedicated to making software imaging solutions freely available.           %
  23 %                                                                             %
  24 %  You may not use this file except in compliance with the License.  You may  %
  25 %  obtain a copy of the License at                                            %
  26 %                                                                             %
  27 %    http://www.imagemagick.org/script/license.php                            %
  28 %                                                                             %
  29 %  Unless required by applicable law or agreed to in writing, software        %
  30 %  distributed under the License is distributed on an "AS IS" BASIS,          %
  31 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
  32 %  See the License for the specific language governing permissions and        %
  33 %  limitations under the License.                                             %
  34 %                                                                             %
  35 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  36 %
  37 %
  38 */
  39 \f
  40 /*
  41   Include declarations.
  42 */
  43 #include "magick/studio.h"
  44 #include "magick/artifact.h"
  45 #include "magick/color-private.h"
  46 #include "magick/cache.h"
  47 #include "magick/draw.h"
  48 #include "magick/exception-private.h"
  49 #include "magick/gem.h"
  50 #include "magick/image.h"
  51 #include "magick/image-private.h"
  52 #include "magick/log.h"
  53 #include "magick/magick.h"
  54 #include "magick/memory_.h"
  55 #include "magick/pixel.h"
  56 #include "magick/pixel-private.h"
  57 #include "magick/quantum.h"
  58 #include "magick/random_.h"
  59 #include "magick/resample.h"
  60 #include "magick/resize.h"
  61 #include "magick/resize-private.h"
  62 #include "magick/transform.h"
  63 #include "magick/signature-private.h"
  64 #include "magick/utility.h"
  65 /*
  66   EWA Resampling Options
  67 */
  68
  69 /* select ONE resampling method */
  70 #define EWA 1                 /* Normal EWA handling - raw or clamped */
  71                               /* if 0 then use "High Quality EWA" */
  72 #define EWA_CLAMP 1           /* EWA Clamping from Nicolas Robidoux */
  73
  74 #define FILTER_LUT 1          /* Use a LUT rather then direct filter calls */
  75
  76 /* output debugging information */
  77 #define DEBUG_ELLIPSE 0       /* output ellipse info for debug */
  78 #define DEBUG_HIT_MISS 0      /* output hit/miss pixels (as gnuplot commands) */
  79 #define DEBUG_NO_PIXEL_HIT 0  /* Make pixels that fail to hit anything - RED */
  80
  81 #if ! FILTER_DIRECT
  82 #define WLUT_WIDTH 1024       /* size of the filter cache */
  83 #endif
  84
  85 /*
  86   Typedef declarations.
  87 */
  88 struct _ResampleFilter
  89 {
  90   CacheView
  91     *view;
  92
  93   Image
  94     *image;
  95
  96   ExceptionInfo
  97     *exception;
  98
  99   MagickBooleanType
 100     debug;
 101
 102   /* Information about image being resampled */
 103   ssize_t
 104     image_area;
 105
 106   InterpolatePixelMethod
 107     interpolate;
 108
 109   VirtualPixelMethod
 110     virtual_pixel;
 111
 112   FilterTypes
 113     filter;
 114
 115   /* processing settings needed */
 116   MagickBooleanType
 117     limit_reached,
 118     do_interpolate,
 119     average_defined;
 120
 121   MagickPixelPacket
 122     average_pixel;
 123
 124   /* current ellipitical area being resampled around center point */
 125   double
 126     A, B, C,
 127     Vlimit, Ulimit, Uwidth, slope;
 128
 129 #if FILTER_LUT
 130   /* LUT of weights for filtered average in elliptical area */
 131   double
 132     filter_lut[WLUT_WIDTH];
 133 #else
 134   /* Use a Direct call to the filter functions */
 135   ResizeFilter
 136     *filter_def;
 137
 138   double
 139     F;
 140 #endif
 141
 142   /* the practical working support of the filter */
 143   double
 144     support;
 145
 146   size_t
 147     signature;
 148 };
 149 \f
 150 /*
 151 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 152 %                                                                             %
 153 %                                                                             %
 154 %                                                                             %
 155 %   A c q u i r e R e s a m p l e I n f o                                     %
 156 %                                                                             %
 157 %                                                                             %
 158 %                                                                             %
 159 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 160 %
 161 %  AcquireResampleFilter() initializes the information resample needs do to a
 162 %  scaled lookup of a color from an image, using area sampling.
 163 %
 164 %  The algorithm is based on a Elliptical Weighted Average, where the pixels
 165 %  found in a large elliptical area is averaged together according to a
 166 %  weighting (filter) function.  For more details see "Fundamentals of Texture
 167 %  Mapping and Image Warping" a master's thesis by Paul.S.Heckbert, June 17,
 168 %  1989.  Available for free from, http://www.cs.cmu.edu/~ph/
 169 %
 170 %  As EWA resampling (or any sort of resampling) can require a lot of
 171 %  calculations to produce a distorted scaling of the source image for each
 172 %  output pixel, the ResampleFilter structure generated holds that information
 173 %  between individual image resampling.
 174 %
 175 %  This function will make the appropriate AcquireCacheView() calls
 176 %  to view the image, calling functions do not need to open a cache view.
 177 %
 178 %  Usage Example...
 179 %      resample_filter=AcquireResampleFilter(image,exception);
 180 %      SetResampleFilter(resample_filter, GaussianFilter, 1.0);
 181 %      for (y=0; y < (ssize_t) image->rows; y++) {
 182 %        for (x=0; x < (ssize_t) image->columns; x++) {
 183 %          u= ....;   v= ....;
 184 %          ScaleResampleFilter(resample_filter, ... scaling vectors ...);
 185 %          (void) ResamplePixelColor(resample_filter,u,v,&pixel);
 186 %          ... assign resampled pixel value ...
 187 %        }
 188 %      }
 189 %      DestroyResampleFilter(resample_filter);
 190 %
 191 %  The format of the AcquireResampleFilter method is:
 192 %
 193 %     ResampleFilter *AcquireResampleFilter(const Image *image,
 194 %       ExceptionInfo *exception)
 195 %
 196 %  A description of each parameter follows:
 197 %
 198 %    o image: the image.
 199 %
 200 %    o exception: return any errors or warnings in this structure.
 201 %
 202 */
 203 MagickExport ResampleFilter *AcquireResampleFilter(const Image *image,
 204   ExceptionInfo *exception)
 205 {
 206   register ResampleFilter
 207     *resample_filter;
 208
 209   assert(image != (Image *) NULL);
 210   assert(image->signature == MagickSignature);
 211   if (image->debug != MagickFalse)
 212     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
 213   assert(exception != (ExceptionInfo *) NULL);
 214   assert(exception->signature == MagickSignature);
 215
 216   resample_filter=(ResampleFilter *) AcquireMagickMemory(
 217     sizeof(*resample_filter));
 218   if (resample_filter == (ResampleFilter *) NULL)
 219     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
 220   (void) ResetMagickMemory(resample_filter,0,sizeof(*resample_filter));
 221
 222   resample_filter->exception=exception;
 223   resample_filter->image=ReferenceImage((Image *) image);
 224   resample_filter->view=AcquireCacheView(resample_filter->image);
 225
 226   resample_filter->debug=IsEventLogging();
 227   resample_filter->signature=MagickSignature;
 228
 229   resample_filter->image_area=(ssize_t) (image->columns*image->rows);
 230   resample_filter->average_defined = MagickFalse;
 231
 232   /* initialise the resampling filter settings */
 233   SetResampleFilter(resample_filter, image->filter, image->blur);
 234   (void) SetResampleFilterInterpolateMethod(resample_filter,
 235     image->interpolate);
 236   (void) SetResampleFilterVirtualPixelMethod(resample_filter,
 237     GetImageVirtualPixelMethod(image));
 238
 239   return(resample_filter);
 240 }
 241 \f
 242 /*
 243 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 244 %                                                                             %
 245 %                                                                             %
 246 %                                                                             %
 247 %   D e s t r o y R e s a m p l e I n f o                                     %
 248 %                                                                             %
 249 %                                                                             %
 250 %                                                                             %
 251 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 252 %
 253 %  DestroyResampleFilter() finalizes and cleans up the resampling
 254 %  resample_filter as returned by AcquireResampleFilter(), freeing any memory
 255 %  or other information as needed.
 256 %
 257 %  The format of the DestroyResampleFilter method is:
 258 %
 259 %      ResampleFilter *DestroyResampleFilter(ResampleFilter *resample_filter)
 260 %
 261 %  A description of each parameter follows:
 262 %
 263 %    o resample_filter: resampling information structure
 264 %
 265 */
 266 MagickExport ResampleFilter *DestroyResampleFilter(
 267   ResampleFilter *resample_filter)
 268 {
 269   assert(resample_filter != (ResampleFilter *) NULL);
 270   assert(resample_filter->signature == MagickSignature);
 271   assert(resample_filter->image != (Image *) NULL);
 272   if (resample_filter->debug != MagickFalse)
 273     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
 274       resample_filter->image->filename);
 275   resample_filter->view=DestroyCacheView(resample_filter->view);
 276   resample_filter->image=DestroyImage(resample_filter->image);
 277 #if ! FILTER_LUT
 278   resample_filter->filter_def=DestroyResizeFilter(resample_filter->filter_def);
 279 #endif
 280   resample_filter->signature=(~MagickSignature);
 281   resample_filter=(ResampleFilter *) RelinquishMagickMemory(resample_filter);
 282   return(resample_filter);
 283 }
 284 \f
 285 /*
 286 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 287 %                                                                             %
 288 %                                                                             %
 289 %                                                                             %
 290 %   R e s a m p l e P i x e l C o l o r                                       %
 291 %                                                                             %
 292 %                                                                             %
 293 %                                                                             %
 294 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 295 %
 296 %  ResamplePixelColor() samples the pixel values surrounding the location
 297 %  given using an elliptical weighted average, at the scale previously
 298 %  calculated, and in the most efficent manner possible for the
 299 %  VirtualPixelMethod setting.
 300 %
 301 %  The format of the ResamplePixelColor method is:
 302 %
 303 %     MagickBooleanType ResamplePixelColor(ResampleFilter *resample_filter,
 304 %       const double u0,const double v0,MagickPixelPacket *pixel)
 305 %
 306 %  A description of each parameter follows:
 307 %
 308 %    o resample_filter: the resample filter.
 309 %
 310 %    o u0,v0: A double representing the center of the area to resample,
 311 %        The distortion transformed transformed x,y coordinate.
 312 %
 313 %    o pixel: the resampled pixel is returned here.
 314 %
 315 */
 316 MagickExport MagickBooleanType ResamplePixelColor(
 317   ResampleFilter *resample_filter,const double u0,const double v0,
 318   MagickPixelPacket *pixel)
 319 {
 320   MagickBooleanType
 321     status;
 322
 323   ssize_t u,v, v1, v2, uw, hit;
 324   double u1;
 325   double U,V,Q,DQ,DDQ;
 326   double divisor_c,divisor_m;
 327   register double weight;
 328   register const PixelPacket *pixels;
 329   register const IndexPacket *indexes;
 330   assert(resample_filter != (ResampleFilter *) NULL);
 331   assert(resample_filter->signature == MagickSignature);
 332
 333   status=MagickTrue;
 334   /* GetMagickPixelPacket(resample_filter->image,pixel); */
 335   if ( resample_filter->do_interpolate ) {
 336     status=InterpolateMagickPixelPacket(resample_filter->image,
 337       resample_filter->view,resample_filter->interpolate,u0,v0,pixel,
 338       resample_filter->exception);
 339     return(status);
 340   }
 341
 342 #if DEBUG_ELLIPSE
 343   (void) FormatLocaleFile(stderr, "u0=%lf; v0=%lf;\n", u0, v0);
 344 #endif
 345
 346   /*
 347     Does resample area Miss the image?
 348     And is that area a simple solid color - then return that color
 349   */
 350   hit = 0;
 351   switch ( resample_filter->virtual_pixel ) {
 352     case BackgroundVirtualPixelMethod:
 353     case ConstantVirtualPixelMethod:
 354     case TransparentVirtualPixelMethod:
 355     case BlackVirtualPixelMethod:
 356     case GrayVirtualPixelMethod:
 357     case WhiteVirtualPixelMethod:
 358     case MaskVirtualPixelMethod:
 359       if ( resample_filter->limit_reached
 360            || u0 + resample_filter->Ulimit < 0.0
 361            || u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 362            || v0 + resample_filter->Vlimit < 0.0
 363            || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows
 364            )
 365         hit++;
 366       break;
 367
 368     case UndefinedVirtualPixelMethod:
 369     case EdgeVirtualPixelMethod:
 370       if (    ( u0 + resample_filter->Ulimit < 0.0 && v0 + resample_filter->Vlimit < 0.0 )
 371            || ( u0 + resample_filter->Ulimit < 0.0
 372                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows )
 373            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 374                 && v0 + resample_filter->Vlimit < 0.0 )
 375            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 376                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows )
 377            )
 378         hit++;
 379       break;
 380     case HorizontalTileVirtualPixelMethod:
 381       if (    v0 + resample_filter->Vlimit < 0.0
 382            || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows
 383            )
 384         hit++;  /* outside the horizontally tiled images. */
 385       break;
 386     case VerticalTileVirtualPixelMethod:
 387       if (    u0 + resample_filter->Ulimit < 0.0
 388            || u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
 389            )
 390         hit++;  /* outside the vertically tiled images. */
 391       break;
 392     case DitherVirtualPixelMethod:
 393       if (    ( u0 + resample_filter->Ulimit < -32.0 && v0 + resample_filter->Vlimit < -32.0 )
 394            || ( u0 + resample_filter->Ulimit < -32.0
 395                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+32.0 )
 396            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+32.0
 397                 && v0 + resample_filter->Vlimit < -32.0 )
 398            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+32.0
 399                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+32.0 )
 400            )
 401         hit++;
 402       break;
 403     case TileVirtualPixelMethod:
 404     case MirrorVirtualPixelMethod:
 405     case RandomVirtualPixelMethod:
 406     case HorizontalTileEdgeVirtualPixelMethod:
 407     case VerticalTileEdgeVirtualPixelMethod:
 408     case CheckerTileVirtualPixelMethod:
 409       /* resampling of area is always needed - no VP limits */
 410       break;
 411   }
 412   if ( hit ) {
 413     /* whole area is a solid color -- just return that color */
 414     status=InterpolateMagickPixelPacket(resample_filter->image,
 415       resample_filter->view,IntegerInterpolatePixel,u0,v0,pixel,
 416       resample_filter->exception);
 417     return(status);
 418   }
 419
 420   /*
 421     Scaling limits reached, return an 'averaged' result.
 422   */
 423   if ( resample_filter->limit_reached ) {
 424     switch ( resample_filter->virtual_pixel ) {
 425       /*  This is always handled by the above, so no need.
 426         case BackgroundVirtualPixelMethod:
 427         case ConstantVirtualPixelMethod:
 428         case TransparentVirtualPixelMethod:
 429         case GrayVirtualPixelMethod,
 430         case WhiteVirtualPixelMethod
 431         case MaskVirtualPixelMethod:
 432       */
 433       case UndefinedVirtualPixelMethod:
 434       case EdgeVirtualPixelMethod:
 435       case DitherVirtualPixelMethod:
 436       case HorizontalTileEdgeVirtualPixelMethod:
 437       case VerticalTileEdgeVirtualPixelMethod:
 438         /* We need an average edge pixel, from the correct edge!
 439            How should I calculate an average edge color?
 440            Just returning an averaged neighbourhood,
 441            works well in general, but falls down for TileEdge methods.
 442            This needs to be done properly!!!!!!
 443         */
 444         status=InterpolateMagickPixelPacket(resample_filter->image,
 445           resample_filter->view,AverageInterpolatePixel,u0,v0,pixel,
 446           resample_filter->exception);
 447         break;
 448       case HorizontalTileVirtualPixelMethod:
 449       case VerticalTileVirtualPixelMethod:
 450         /* just return the background pixel - Is there more direct way? */
 451         status=InterpolateMagickPixelPacket(resample_filter->image,
 452           resample_filter->view,IntegerInterpolatePixel,-1.0,-1.0,pixel,
 453           resample_filter->exception);
 454         break;
 455       case TileVirtualPixelMethod:
 456       case MirrorVirtualPixelMethod:
 457       case RandomVirtualPixelMethod:
 458       case CheckerTileVirtualPixelMethod:
 459       default:
 460         /* generate a average color of the WHOLE image */
 461         if ( resample_filter->average_defined == MagickFalse ) {
 462           Image
 463             *average_image;
 464
 465           CacheView
 466             *average_view;
 467
 468           GetMagickPixelPacket(resample_filter->image,(MagickPixelPacket *)
 469             &resample_filter->average_pixel);
 470           resample_filter->average_defined=MagickTrue;
 471
 472           /* Try to get an averaged pixel color of whole image */
 473           average_image=ResizeImage(resample_filter->image,1,1,BoxFilter,1.0,
 474             resample_filter->exception);
 475           if (average_image == (Image *) NULL)
 476             {
 477               *pixel=resample_filter->average_pixel; /* FAILED */
 478               break;
 479             }
 480           average_view=AcquireCacheView(average_image);
 481           pixels=(PixelPacket *)GetCacheViewVirtualPixels(average_view,0,0,1,1,
 482             resample_filter->exception);
 483           if (pixels == (const PixelPacket *) NULL) {
 484             average_view=DestroyCacheView(average_view);
 485             average_image=DestroyImage(average_image);
 486             *pixel=resample_filter->average_pixel; /* FAILED */
 487             break;
 488           }
 489           indexes=(IndexPacket *) GetCacheViewAuthenticIndexQueue(average_view);
 490           SetMagickPixelPacket(resample_filter->image,pixels,indexes,
 491             &(resample_filter->average_pixel));
 492           average_view=DestroyCacheView(average_view);
 493           average_image=DestroyImage(average_image);
 494
 495           if ( resample_filter->virtual_pixel == CheckerTileVirtualPixelMethod )
 496             {
 497               /* CheckerTile is avergae of image average half background */
 498               /* FUTURE: replace with a 50% blend of both pixels */
 499
 500               weight = QuantumScale*((MagickRealType)(QuantumRange-
 501                           resample_filter->average_pixel.opacity));
 502               resample_filter->average_pixel.red *= weight;
 503               resample_filter->average_pixel.green *= weight;
 504               resample_filter->average_pixel.blue *= weight;
 505               divisor_c = weight;
 506
 507               weight = QuantumScale*((MagickRealType)(QuantumRange-
 508                           resample_filter->image->background_color.opacity));
 509               resample_filter->average_pixel.red +=
 510                       weight*resample_filter->image->background_color.red;
 511               resample_filter->average_pixel.green +=
 512                       weight*resample_filter->image->background_color.green;
 513               resample_filter->average_pixel.blue +=
 514                       weight*resample_filter->image->background_color.blue;
 515               resample_filter->average_pixel.opacity +=
 516                       resample_filter->image->background_color.opacity;
 517               divisor_c += weight;
 518
 519               resample_filter->average_pixel.red /= divisor_c;
 520               resample_filter->average_pixel.green /= divisor_c;
 521               resample_filter->average_pixel.blue /= divisor_c;
 522               resample_filter->average_pixel.opacity /= 2;
 523
 524             }
 525         }
 526         *pixel=resample_filter->average_pixel;
 527         break;
 528     }
 529     return(status);
 530   }
 531
 532   /*
 533     Initialize weighted average data collection
 534   */
 535   hit = 0;
 536   divisor_c = 0.0;
 537   divisor_m = 0.0;
 538   pixel->red = pixel->green = pixel->blue = 0.0;
 539   if (pixel->matte != MagickFalse) pixel->opacity = 0.0;
 540   if (pixel->colorspace == CMYKColorspace) pixel->index = 0.0;
 541
 542   /*
 543     Determine the parellelogram bounding box fitted to the ellipse
 544     centered at u0,v0.  This area is bounding by the lines...
 545   */
 546   v1 = (ssize_t)ceil(v0 - resample_filter->Vlimit);  /* range of scan lines */
 547   v2 = (ssize_t)floor(v0 + resample_filter->Vlimit);
 548
 549   /* scan line start and width accross the parallelogram */
 550   u1 = u0 + (v1-v0)*resample_filter->slope - resample_filter->Uwidth;
 551   uw = (ssize_t)(2.0*resample_filter->Uwidth)+1;
 552
 553 #if DEBUG_ELLIPSE
 554   (void) FormatLocaleFile(stderr, "v1=%ld; v2=%ld\n", (long)v1, (long)v2);
 555   (void) FormatLocaleFile(stderr, "u1=%ld; uw=%ld\n", (long)u1, (long)uw);
 556 #else
 557 # define DEBUG_HIT_MISS 0 /* only valid if DEBUG_ELLIPSE is enabled */
 558 #endif
 559
 560   /*
 561     Do weighted resampling of all pixels,  within the scaled ellipse,
 562     bound by a Parellelogram fitted to the ellipse.
 563   */
 564   DDQ = 2*resample_filter->A;
 565   for( v=v1; v<=v2;  v++ ) {
 566 #if DEBUG_HIT_MISS
 567     long uu = ceil(u1);   /* actual pixel location (for debug only) */
 568     (void) FormatLocaleFile(stderr, "# scan line from pixel %ld, %ld\n", (long)uu, (long)v);
 569 #endif
 570     u = (ssize_t)ceil(u1);        /* first pixel in scanline */
 571     u1 += resample_filter->slope; /* start of next scan line */
 572
 573
 574     /* location of this first pixel, relative to u0,v0 */
 575     U = (double)u-u0;
 576     V = (double)v-v0;
 577
 578     /* Q = ellipse quotent ( if Q<F then pixel is inside ellipse) */
 579     Q = (resample_filter->A*U + resample_filter->B*V)*U + resample_filter->C*V*V;
 580     DQ = resample_filter->A*(2.0*U+1) + resample_filter->B*V;
 581
 582     /* get the scanline of pixels for this v */
 583     pixels=GetCacheViewVirtualPixels(resample_filter->view,u,v,(size_t) uw,
 584       1,resample_filter->exception);
 585     if (pixels == (const PixelPacket *) NULL)
 586       return(MagickFalse);
 587     indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
 588
 589     /* count up the weighted pixel colors */
 590     for( u=0; u<uw; u++ ) {
 591 #if FILTER_LUT
 592       /* Note that the ellipse has been pre-scaled so F = WLUT_WIDTH */
 593       if ( Q < (double)WLUT_WIDTH ) {
 594         weight = resample_filter->filter_lut[(int)Q];
 595 #else
 596       /* Note that the ellipse has been pre-scaled so F = support^2 */
 597       if ( Q < (double)resample_filter->F ) {
 598         weight = GetResizeFilterWeight(resample_filter->filter_def,
 599              sqrt(Q));    /* a SquareRoot!  Arrggghhhhh... */
 600 #endif
 601
 602         pixel->opacity  += weight*pixels->opacity;
 603         divisor_m += weight;
 604
 605         if (pixel->matte != MagickFalse)
 606           weight *= QuantumScale*((MagickRealType)(QuantumRange-pixels->opacity));
 607         pixel->red   += weight*pixels->red;
 608         pixel->green += weight*pixels->green;
 609         pixel->blue  += weight*pixels->blue;
 610         if (pixel->colorspace == CMYKColorspace)
 611           pixel->index += weight*(*indexes);
 612         divisor_c += weight;
 613
 614         hit++;
 615 #if DEBUG_HIT_MISS
 616         /* mark the pixel according to hit/miss of the ellipse */
 617         (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
 618                      (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
 619         (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
 620                      (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
 621       } else {
 622         (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
 623                      (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
 624         (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
 625                      (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
 626       }
 627       uu++;
 628 #else
 629       }
 630 #endif
 631       pixels++;
 632       indexes++;
 633       Q += DQ;
 634       DQ += DDQ;
 635     }
 636   }
 637 #if DEBUG_ELLIPSE
 638   (void) FormatLocaleFile(stderr, "Hit=%ld;  Total=%ld;\n", (long)hit, (long)uw*(v2-v1) );
 639 #endif
 640
 641   /*
 642     Result sanity check -- this should NOT happen
 643   */
 644   if ( hit == 0 ) {
 645     /* not enough pixels in resampling, resort to direct interpolation */
 646 #if DEBUG_NO_PIXEL_HIT
 647     pixel->opacity = pixel->red = pixel->green = pixel->blue = 0;
 648     pixel->red = QuantumRange; /* show pixels for which EWA fails */
 649 #else
 650     status=InterpolateMagickPixelPacket(resample_filter->image,
 651       resample_filter->view,resample_filter->interpolate,u0,v0,pixel,
 652       resample_filter->exception);
 653 #endif
 654     return status;
 655   }
 656
 657   /*
 658     Finialize results of resampling
 659   */
 660   divisor_m = 1.0/divisor_m;
 661   pixel->opacity = (MagickRealType) ClampToQuantum(divisor_m*pixel->opacity);
 662   divisor_c = 1.0/divisor_c;
 663   pixel->red   = (MagickRealType) ClampToQuantum(divisor_c*pixel->red);
 664   pixel->green = (MagickRealType) ClampToQuantum(divisor_c*pixel->green);
 665   pixel->blue  = (MagickRealType) ClampToQuantum(divisor_c*pixel->blue);
 666   if (pixel->colorspace == CMYKColorspace)
 667     pixel->index = (MagickRealType) ClampToQuantum(divisor_c*pixel->index);
 668   return(MagickTrue);
 669 }
 670 \f
 671 #if EWA && EWA_CLAMP
 672 /*
 673 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 674 %                                                                             %
 675 %                                                                             %
 676 %                                                                             %
 677 -   C l a m p U p A x e s                                                     %
 678 %                                                                             %
 679 %                                                                             %
 680 %                                                                             %
 681 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 682 %
 683 % ClampUpAxes() function converts the input vectors into a major and
 684 % minor axis unit vectors, and their magnitude.  This allows us to
 685 % ensure that the ellipse generated is never smaller than the unit
 686 % circle and thus never too small for use in EWA resampling.
 687 %
 688 % This purely mathematical 'magic' was provided by Professor Nicolas
 689 % Robidoux and his Masters student Chantal Racette.
 690 %
 691 % Reference: "We Recommend Singular Value Decomposition", David Austin
 692 %   http://www.ams.org/samplings/feature-column/fcarc-svd
 693 %
 694 % By generating major and minor axis vectors, we can actually use the
 695 % ellipse in its "canonical form", by remapping the dx,dy of the
 696 % sampled point into distances along the major and minor axis unit
 697 % vectors.
 698 %
 699 % Reference: http://en.wikipedia.org/wiki/Ellipse#Canonical_form
 700 */
 701 static inline void ClampUpAxes(const double dux,
 702                                const double dvx,
 703                                const double duy,
 704                                const double dvy,
 705                                double *major_mag,
 706                                double *minor_mag,
 707                                double *major_unit_x,
 708                                double *major_unit_y,
 709                                double *minor_unit_x,
 710                                double *minor_unit_y)
 711 {
 712   /*
 713    * ClampUpAxes takes an input 2x2 matrix
 714    *
 715    * [ a b ] = [ dux duy ]
 716    * [ c d ] = [ dvx dvy ]
 717    *
 718    * and computes from it the major and minor axis vectors [major_x,
 719    * major_y] and [minor_x,minor_y] of the smallest ellipse containing
 720    * both the unit disk and the ellipse which is the image of the unit
 721    * disk by the linear transformation
 722    *
 723    * [ dux duy ] [S] = [s]
 724    * [ dvx dvy ] [T] = [t]
 725    *
 726    * (The vector [S,T] is the difference between a position in output
 727    * space and [X,Y]; the vector [s,t] is the difference between a
 728    * position in input space and [x,y].)
 729    */
 730   /*
 731    * Output:
 732    *
 733    * major_mag is the half-length of the major axis of the "new"
 734    * ellipse.
 735    *
 736    * minor_mag is the half-length of the minor axis of the "new"
 737    * ellipse.
 738    *
 739    * major_unit_x is the x-coordinate of the major axis direction vector
 740    * of both the "old" and "new" ellipses.
 741    *
 742    * major_unit_y is the y-coordinate of the major axis direction vector.
 743    *
 744    * minor_unit_x is the x-coordinate of the minor axis direction vector.
 745    *
 746    * minor_unit_y is the y-coordinate of the minor axis direction vector.
 747    *
 748    * Unit vectors are useful for computing projections, in particular,
 749    * to compute the distance between a point in output space and the
 750    * center of a unit disk in output space, using the position of the
 751    * corresponding point [s,t] in input space. Following the clamping,
 752    * the square of this distance is
 753    *
 754    * ( ( s * major_unit_x + t * major_unit_y ) / major_mag )^2
 755    * +
 756    * ( ( s * minor_unit_x + t * minor_unit_y ) / minor_mag )^2
 757    *
 758    * If such distances will be computed for many [s,t]'s, it makes
 759    * sense to actually compute the reciprocal of major_mag and
 760    * minor_mag and multiply them by the above unit lengths.
 761    *
 762    * Now, if you want to modify the input pair of tangent vectors so
 763    * that it defines the modified ellipse, all you have to do is set
 764    *
 765    * newdux = major_mag * major_unit_x
 766    * newdvx = major_mag * major_unit_y
 767    * newduy = minor_mag * minor_unit_x = minor_mag * -major_unit_y
 768    * newdvy = minor_mag * minor_unit_y = minor_mag *  major_unit_x
 769    *
 770    * and use these tangent vectors as if they were the original ones.
 771    * Usually, this is a drastic change in the tangent vectors even if
 772    * the singular values are not clamped; for example, the minor axis
 773    * vector always points in a direction which is 90 degrees
 774    * counterclockwise from the direction of the major axis vector.
 775    */
 776   /*
 777    * Discussion:
 778    *
 779    * GOAL: Fix things so that the pullback, in input space, of a disk
 780    * of radius r in output space is an ellipse which contains, at
 781    * least, a disc of radius r. (Make this hold for any r>0.)
 782    *
 783    * ESSENCE OF THE METHOD: Compute the product of the first two
 784    * factors of an SVD of the linear transformation defining the
 785    * ellipse and make sure that both its columns have norm at least 1.
 786    * Because rotations and reflexions map disks to themselves, it is
 787    * not necessary to compute the third (rightmost) factor of the SVD.
 788    *
 789    * DETAILS: Find the singular values and (unit) left singular
 790    * vectors of Jinv, clampling up the singular values to 1, and
 791    * multiply the unit left singular vectors by the new singular
 792    * values in order to get the minor and major ellipse axis vectors.
 793    *
 794    * Image resampling context:
 795    *
 796    * The Jacobian matrix of the transformation at the output point
 797    * under consideration is defined as follows:
 798    *
 799    * Consider the transformation (x,y) -> (X,Y) from input locations
 800    * to output locations. (Anthony Thyssen, elsewhere in resample.c,
 801    * uses the notation (u,v) -> (x,y).)
 802    *
 803    * The Jacobian matrix of the transformation at (x,y) is equal to
 804    *
 805    *   J = [ A, B ] = [ dX/dx, dX/dy ]
 806    *       [ C, D ]   [ dY/dx, dY/dy ]
 807    *
 808    * that is, the vector [A,C] is the tangent vector corresponding to
 809    * input changes in the horizontal direction, and the vector [B,D]
 810    * is the tangent vector corresponding to input changes in the
 811    * vertical direction.
 812    *
 813    * In the context of resampling, it is natural to use the inverse
 814    * Jacobian matrix Jinv because resampling is generally performed by
 815    * pulling pixel locations in the output image back to locations in
 816    * the input image. Jinv is
 817    *
 818    *   Jinv = [ a, b ] = [ dx/dX, dx/dY ]
 819    *          [ c, d ]   [ dy/dX, dy/dY ]
 820    *
 821    * Note: Jinv can be computed from J with the following matrix
 822    * formula:
 823    *
 824    *   Jinv = 1/(A*D-B*C) [  D, -B ]
 825    *                      [ -C,  A ]
 826    *
 827    * What we do is modify Jinv so that it generates an ellipse which
 828    * is as close as possible to the original but which contains the
 829    * unit disk. This can be accomplished as follows:
 830    *
 831    * Let
 832    *
 833    *   Jinv = U Sigma V^T
 834    *
 835    * be an SVD decomposition of Jinv. (The SVD is not unique, but the
 836    * final ellipse does not depend on the particular SVD.)
 837    *
 838    * We could clamp up the entries of the diagonal matrix Sigma so
 839    * that they are at least 1, and then set
 840    *
 841    *   Jinv = U newSigma V^T.
 842    *
 843    * However, we do not need to compute V for the following reason:
 844    * V^T is an orthogonal matrix (that is, it represents a combination
 845    * of rotations and reflexions) so that it maps the unit circle to
 846    * itself. For this reason, the exact value of V does not affect the
 847    * final ellipse, and we can choose V to be the identity
 848    * matrix. This gives
 849    *
 850    *   Jinv = U newSigma.
 851    *
 852    * In the end, we return the two diagonal entries of newSigma
 853    * together with the two columns of U.
 854    */
 855   /*
 856    * ClampUpAxes was written by Nicolas Robidoux and Chantal Racette
 857    * of Laurentian University with insightful suggestions from Anthony
 858    * Thyssen and funding from the National Science and Engineering
 859    * Research Council of Canada. It is distinguished from its
 860    * predecessors by its efficient handling of degenerate cases.
 861    *
 862    * The idea of clamping up the EWA ellipse's major and minor axes so
 863    * that the result contains the reconstruction kernel filter support
 864    * is taken from Andreas Gustaffson's Masters thesis "Interactive
 865    * Image Warping", Helsinki University of Technology, Faculty of
 866    * Information Technology, 59 pages, 1993 (see Section 3.6).
 867    *
 868    * The use of the SVD to clamp up the singular values of the
 869    * Jacobian matrix of the pullback transformation for EWA resampling
 870    * is taken from the astrophysicist Craig DeForest.  It is
 871    * implemented in his PDL::Transform code (PDL = Perl Data
 872    * Language).
 873    */
 874   const double a = dux;
 875   const double b = duy;
 876   const double c = dvx;
 877   const double d = dvy;
 878   /*
 879    * n is the matrix Jinv * transpose(Jinv). Eigenvalues of n are the
 880    * squares of the singular values of Jinv.
 881    */
 882   const double aa = a*a;
 883   const double bb = b*b;
 884   const double cc = c*c;
 885   const double dd = d*d;
 886   /*
 887    * Eigenvectors of n are left singular vectors of Jinv.
 888    */
 889   const double n11 = aa+bb;
 890   const double n12 = a*c+b*d;
 891   const double n21 = n12;
 892   const double n22 = cc+dd;
 893   const double det = a*d-b*c;
 894   const double twice_det = det+det;
 895   const double frobenius_squared = n11+n22;
 896   const double discriminant =
 897     (frobenius_squared+twice_det)*(frobenius_squared-twice_det);
 898   const double sqrt_discriminant = sqrt(discriminant);
 899   /*
 900    * s1 is the largest singular value of the inverse Jacobian
 901    * matrix. In other words, its reciprocal is the smallest singular
 902    * value of the Jacobian matrix itself.
 903    * If s1 = 0, both singular values are 0, and any orthogonal pair of
 904    * left and right factors produces a singular decomposition of Jinv.
 905    */
 906   /*
 907    * Initially, we only compute the squares of the singular values.
 908    */
 909   const double s1s1 = 0.5*(frobenius_squared+sqrt_discriminant);
 910   /*
 911    * s2 the smallest singular value of the inverse Jacobian
 912    * matrix. Its reciprocal is the largest singular value of the
 913    * Jacobian matrix itself.
 914    */
 915   const double s2s2 = 0.5*(frobenius_squared-sqrt_discriminant);
 916   const double s1s1minusn11 = s1s1-n11;
 917   const double s1s1minusn22 = s1s1-n22;
 918   /*
 919    * u1, the first column of the U factor of a singular decomposition
 920    * of Jinv, is a (non-normalized) left singular vector corresponding
 921    * to s1. It has entries u11 and u21. We compute u1 from the fact
 922    * that it is an eigenvector of n corresponding to the eigenvalue
 923    * s1^2.
 924    */
 925   const double s1s1minusn11_squared = s1s1minusn11*s1s1minusn11;
 926   const double s1s1minusn22_squared = s1s1minusn22*s1s1minusn22;
 927   /*
 928    * The following selects the largest row of n-s1^2 I as the one
 929    * which is used to find the eigenvector. If both s1^2-n11 and
 930    * s1^2-n22 are zero, n-s1^2 I is the zero matrix.  In that case,
 931    * any vector is an eigenvector; in addition, norm below is equal to
 932    * zero, and, in exact arithmetic, this is the only case in which
 933    * norm = 0. So, setting u1 to the simple but arbitrary vector [1,0]
 934    * if norm = 0 safely takes care of all cases.
 935    */
 936   const double temp_u11 =
 937     ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? n12 : s1s1minusn22 );
 938   const double temp_u21 =
 939     ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? s1s1minusn11 : n21 );
 940   const double norm = sqrt(temp_u11*temp_u11+temp_u21*temp_u21);
 941   /*
 942    * Finalize the entries of first left singular vector (associated
 943    * with the largest singular value).
 944    */
 945   const double u11 = ( (norm>0.0) ? temp_u11/norm : 1.0 );
 946   const double u21 = ( (norm>0.0) ? temp_u21/norm : 0.0 );
 947   /*
 948    * Clamp the singular values up to 1.
 949    */
 950   *major_mag = ( (s1s1<=1.0) ? 1.0 : sqrt(s1s1) );
 951   *minor_mag = ( (s2s2<=1.0) ? 1.0 : sqrt(s2s2) );
 952   /*
 953    * Return the unit major and minor axis direction vectors.
 954    */
 955   *major_unit_x = u11;
 956   *major_unit_y = u21;
 957   *minor_unit_x = -u21;
 958   *minor_unit_y = u11;
 959 }
 960 \f
 961 #endif
 962 /*
 963 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 964 %                                                                             %
 965 %                                                                             %
 966 %                                                                             %
 967 %   S c a l e R e s a m p l e F i l t e r                                     %
 968 %                                                                             %
 969 %                                                                             %
 970 %                                                                             %
 971 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 972 %
 973 %  ScaleResampleFilter() does all the calculations needed to resample an image
 974 %  at a specific scale, defined by two scaling vectors.  This not using
 975 %  a orthogonal scaling, but two distorted scaling vectors, to allow the
 976 %  generation of a angled ellipse.
 977 %
 978 %  As only two deritive scaling vectors are used the center of the ellipse
 979 %  must be the center of the lookup.  That is any curvature that the
 980 %  distortion may produce is discounted.
 981 %
 982 %  The input vectors are produced by either finding the derivitives of the
 983 %  distortion function, or the partial derivitives from a distortion mapping.
 984 %  They do not need to be the orthogonal dx,dy scaling vectors, but can be
 985 %  calculated from other derivatives.  For example you could use  dr,da/r
 986 %  polar coordinate vector scaling vectors
 987 %
 988 %  If   u,v =  DistortEquation(x,y)   OR   u = Fu(x,y); v = Fv(x,y)
 989 %  Then the scaling vectors are determined from the deritives...
 990 %      du/dx, dv/dx     and    du/dy, dv/dy
 991 %  If the resulting scaling vectors is othogonally aligned then...
 992 %      dv/dx = 0   and   du/dy  =  0
 993 %  Producing an othogonally alligned ellipse in source space for the area to
 994 %  be resampled.
 995 %
 996 %  Note that scaling vectors are different to argument order.  Argument order
 997 %  is the general order the deritives are extracted from the distortion
 998 %  equations, and not the scaling vectors. As such the middle two vaules
 999 %  may be swapped from what you expect.  Caution is advised.
1000 %
1001 %  WARNING: It is assumed that any SetResampleFilter() method call will
1002 %  always be performed before the ScaleResampleFilter() method, so that the
1003 %  size of the ellipse will match the support for the resampling filter being
1004 %  used.
1005 %
1006 %  The format of the ScaleResampleFilter method is:
1007 %
1008 %     void ScaleResampleFilter(const ResampleFilter *resample_filter,
1009 %       const double dux,const double duy,const double dvx,const double dvy)
1010 %
1011 %  A description of each parameter follows:
1012 %
1013 %    o resample_filter: the resampling resample_filterrmation defining the
1014 %      image being resampled
1015 %
1016 %    o dux,duy,dvx,dvy:
1017 %         The deritives or scaling vectors defining the EWA ellipse.
1018 %         NOTE: watch the order, which is based on the order deritives
1019 %         are usally determined from distortion equations (see above).
1020 %         The middle two values may need to be swapped if you are thinking
1021 %         in terms of scaling vectors.
1022 %
1023 */
1024 MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
1025   const double dux,const double duy,const double dvx,const double dvy)
1026 {
1027   double A,B,C,F;
1028
1029   assert(resample_filter != (ResampleFilter *) NULL);
1030   assert(resample_filter->signature == MagickSignature);
1031
1032   resample_filter->limit_reached = MagickFalse;
1033
1034   /* A 'point' filter forces use of interpolation instead of area sampling */
1035   if ( resample_filter->filter == PointFilter )
1036     return; /* EWA turned off - nothing to do */
1037
1038 #if DEBUG_ELLIPSE
1039   (void) FormatLocaleFile(stderr, "# -----\n" );
1040   (void) FormatLocaleFile(stderr, "dux=%lf; dvx=%lf;   duy=%lf; dvy=%lf;\n",
1041        dux, dvx, duy, dvy);
1042 #endif
1043
1044   /* Find Ellipse Coefficents such that
1045         A*u^2 + B*u*v + C*v^2 = F
1046      With u,v relative to point around which we are resampling.
1047      And the given scaling dx,dy vectors in u,v space
1048          du/dx,dv/dx   and  du/dy,dv/dy
1049   */
1050 #if EWA
1051   /* Direct conversion of derivatives into elliptical coefficients
1052      However when magnifying images, the scaling vectors will be small
1053      resulting in a ellipse that is too small to sample properly.
1054      As such we need to clamp the major/minor axis to a minumum of 1.0
1055      to prevent it getting too small.
1056   */
1057 #if EWA_CLAMP
1058   { double major_mag,
1059            minor_mag,
1060            major_x,
1061            major_y,
1062            minor_x,
1063            minor_y;
1064
1065   ClampUpAxes(dux,dvx,duy,dvy, &major_mag, &minor_mag,
1066                 &major_x, &major_y, &minor_x, &minor_y);
1067   major_x *= major_mag;  major_y *= major_mag;
1068   minor_x *= minor_mag;  minor_y *= minor_mag;
1069 #if DEBUG_ELLIPSE
1070   (void) FormatLocaleFile(stderr, "major_x=%lf; major_y=%lf;  minor_x=%lf; minor_y=%lf;\n",
1071         major_x, major_y, minor_x, minor_y);
1072 #endif
1073   A = major_y*major_y+minor_y*minor_y;
1074   B = -2.0*(major_x*major_y+minor_x*minor_y);
1075   C = major_x*major_x+minor_x*minor_x;
1076   F = major_mag*minor_mag;
1077   F *= F; /* square it */
1078   }
1079 #else /* raw unclamped EWA */
1080   A = dvx*dvx+dvy*dvy;
1081   B = -2.0*(dux*dvx+duy*dvy);
1082   C = dux*dux+duy*duy;
1083   F = dux*dvy-duy*dvx;
1084   F *= F; /* square it */
1085 #endif /* EWA_CLAMP */
1086
1087 #else /* HQ_EWA */
1088   /*
1089     This Paul Heckbert's "Higher Quality EWA" formula, from page 60 in his
1090     thesis, which adds a unit circle to the elliptical area so as to do both
1091     Reconstruction and Prefiltering of the pixels in the resampling.  It also
1092     means it is always likely to have at least 4 pixels within the area of the
1093     ellipse, for weighted averaging.  No scaling will result with F == 4.0 and
1094     a circle of radius 2.0, and F smaller than this means magnification is
1095     being used.
1096
1097     NOTE: This method produces a very blury result at near unity scale while
1098     producing perfect results for strong minitification and magnifications.
1099
1100     However filter support is fixed to 2.0 (no good for Windowed Sinc filters)
1101   */
1102   A = dvx*dvx+dvy*dvy+1;
1103   B = -2.0*(dux*dvx+duy*dvy);
1104   C = dux*dux+duy*duy+1;
1105   F = A*C - B*B/4;
1106 #endif
1107
1108 #if DEBUG_ELLIPSE
1109   (void) FormatLocaleFile(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
1110
1111   /* Figure out the various information directly about the ellipse.
1112      This information currently not needed at this time, but may be
1113      needed later for better limit determination.
1114
1115      It is also good to have as a record for future debugging
1116   */
1117   { double alpha, beta, gamma, Major, Minor;
1118     double Eccentricity, Ellipse_Area, Ellipse_Angle;
1119
1120     alpha = A+C;
1121     beta  = A-C;
1122     gamma = sqrt(beta*beta + B*B );
1123
1124     if ( alpha - gamma <= MagickEpsilon )
1125       Major = MagickHuge;
1126     else
1127       Major = sqrt(2*F/(alpha - gamma));
1128     Minor = sqrt(2*F/(alpha + gamma));
1129
1130     (void) FormatLocaleFile(stderr, "# Major=%lf; Minor=%lf\n", Major, Minor );
1131
1132     /* other information about ellipse include... */
1133     Eccentricity = Major/Minor;
1134     Ellipse_Area = MagickPI*Major*Minor;
1135     Ellipse_Angle = atan2(B, A-C);
1136
1137     (void) FormatLocaleFile(stderr, "# Angle=%lf   Area=%lf\n",
1138          RadiansToDegrees(Ellipse_Angle), Ellipse_Area);
1139   }
1140 #endif
1141
1142   /* If one or both of the scaling vectors is impossibly large
1143      (producing a very large raw F value), we may as well not bother
1144      doing any form of resampling since resampled area is very large.
1145      In this case some alternative means of pixel sampling, such as
1146      the average of the whole image is needed to get a reasonable
1147      result. Calculate only as needed.
1148   */
1149   if ( (4*A*C - B*B) > MagickHuge ) {
1150     resample_filter->limit_reached = MagickTrue;
1151     return;
1152   }
1153
1154   /* Scale ellipse to match the filters support
1155      (that is, multiply F by the square of the support).
1156   */
1157   F *= resample_filter->support;
1158   F *= resample_filter->support;
1159
1160   /* Orthogonal bounds of the ellipse */
1161   resample_filter->Ulimit = sqrt(C*F/(A*C-0.25*B*B));
1162   resample_filter->Vlimit = sqrt(A*F/(A*C-0.25*B*B));
1163
1164   /* Horizontally aligned parallelogram fitted to Ellipse */
1165   resample_filter->Uwidth = sqrt(F/A); /* Half of the parallelogram width */
1166   resample_filter->slope = -B/(2.0*A); /* Reciprocal slope of the parallelogram */
1167
1168 #if DEBUG_ELLIPSE
1169   (void) FormatLocaleFile(stderr, "Ulimit=%lf; Vlimit=%lf; UWidth=%lf; Slope=%lf;\n",
1170            resample_filter->Ulimit, resample_filter->Vlimit,
1171            resample_filter->Uwidth, resample_filter->slope );
1172 #endif
1173
1174   /* Check the absolute area of the parallelogram involved.
1175    * This limit needs more work, as it is too slow for larger images
1176    * with tiled views of the horizon.
1177   */
1178   if ( (resample_filter->Uwidth * resample_filter->Vlimit)
1179          > (4.0*resample_filter->image_area)) {
1180     resample_filter->limit_reached = MagickTrue;
1181     return;
1182   }
1183
1184   /* Scale ellipse formula to directly index the Filter Lookup Table */
1185   { register double scale;
1186 #if FILTER_LUT
1187     /* scale so that F = WLUT_WIDTH; -- hardcoded */
1188     scale = (double)WLUT_WIDTH/F;
1189 #else
1190     /* scale so that F = resample_filter->F (support^2) */
1191     scale = resample_filter->F/F;
1192 #endif
1193     resample_filter->A = A*scale;
1194     resample_filter->B = B*scale;
1195     resample_filter->C = C*scale;
1196   }
1197 }
1198 \f
1199 /*
1200 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1201 %                                                                             %
1202 %                                                                             %
1203 %                                                                             %
1204 %   S e t R e s a m p l e F i l t e r                                         %
1205 %                                                                             %
1206 %                                                                             %
1207 %                                                                             %
1208 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1209 %
1210 %  SetResampleFilter() set the resampling filter lookup table based on a
1211 %  specific filter.  Note that the filter is used as a radial filter not as a
1212 %  two pass othogonally aligned resampling filter.
1213 %
1214 %  The default Filter, is Gaussian, which is the standard filter used by the
1215 %  original paper on the Elliptical Weighted Everage Algorithm. However other
1216 %  filters can also be used.
1217 %
1218 %  The format of the SetResampleFilter method is:
1219 %
1220 %    void SetResampleFilter(ResampleFilter *resample_filter,
1221 %      const FilterTypes filter,const double blur)
1222 %
1223 %  A description of each parameter follows:
1224 %
1225 %    o resample_filter: resampling resample_filterrmation structure
1226 %
1227 %    o filter: the resize filter for elliptical weighting LUT
1228 %
1229 %    o blur: filter blur factor (radial scaling) for elliptical weighting LUT
1230 %
1231 */
1232 MagickExport void SetResampleFilter(ResampleFilter *resample_filter,
1233   const FilterTypes filter,const double blur)
1234 {
1235   ResizeFilter
1236      *resize_filter;
1237
1238   assert(resample_filter != (ResampleFilter *) NULL);
1239   assert(resample_filter->signature == MagickSignature);
1240
1241   resample_filter->do_interpolate = MagickFalse;
1242   resample_filter->filter = filter;
1243
1244   if ( filter == PointFilter )
1245     {
1246       resample_filter->do_interpolate = MagickTrue;
1247       return;  /* EWA turned off - nothing more to do */
1248     }
1249
1250   /* Set a default cylindrical filter of a 'low blur' Jinc windowed Jinc */
1251   if ( filter == UndefinedFilter )
1252     resample_filter->filter = RobidouxFilter;
1253
1254   resize_filter = AcquireResizeFilter(resample_filter->image,
1255        resample_filter->filter,blur,MagickTrue,resample_filter->exception);
1256   if (resize_filter == (ResizeFilter *) NULL)
1257     {
1258       (void) ThrowMagickException(resample_filter->exception,GetMagickModule(),
1259            ModuleError, "UnableToSetFilteringValue",
1260            "Fall back to default EWA gaussian filter");
1261       resample_filter->filter = PointFilter;
1262     }
1263
1264   /* Get the practical working support for the filter,
1265    * after any API call blur factors have been accoded for.
1266    */
1267 #if EWA
1268   resample_filter->support = GetResizeFilterSupport(resize_filter);
1269 #else
1270   resample_filter->support = 2.0;  /* fixed support size for HQ-EWA */
1271 #endif
1272
1273 #if FILTER_LUT
1274   /* Fill the LUT with the weights from the selected filter function */
1275   { register int
1276        Q;
1277     double
1278        r_scale;
1279     /* Scale radius so the filter LUT covers the full support range */
1280     r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
1281     for(Q=0; Q<WLUT_WIDTH; Q++)
1282       resample_filter->filter_lut[Q] = (double)
1283            GetResizeFilterWeight(resize_filter,sqrt((double)Q)*r_scale);
1284
1285     /* finished with the resize filter */
1286     resize_filter = DestroyResizeFilter(resize_filter);
1287   }
1288 #else
1289   /* save the filter and the scaled ellipse bounds needed for filter */
1290   resample_filter->filter_def = resize_filter;
1291   resample_filter->F = resample_filter->support*resample_filter->support;
1292 #endif
1293
1294   /*
1295     Adjust the scaling of the default unit circle
1296     This assumes that any real scaling changes will always
1297     take place AFTER the filter method has been initialized.
1298   */
1299   ScaleResampleFilter(resample_filter, 1.0, 0.0, 0.0, 1.0);
1300
1301 #if 0
1302   /* This is old code kept as a reference only.  It is very wrong,
1303      and I don't understand exactly what it was attempting to do.
1304   */
1305   /*
1306     Create Normal Gaussian 2D Filter Weighted Lookup Table.
1307     A normal EWA guassual lookup would use   exp(Q*ALPHA)
1308     where  Q = distance squared from 0.0 (center) to 1.0 (edge)
1309     and    ALPHA = -4.0*ln(2.0)  ==>  -2.77258872223978123767
1310     The table is of length 1024, and equates to support radius of 2.0
1311     thus needs to be scaled by  ALPHA*4/1024 and any blur factor squared
1312
1313     The above came from some reference code provided by Fred Weinhaus
1314     and seems to have been a guess that was appropriate for its use
1315     in a 3d perspective landscape mapping program.
1316   */
1317   r_scale = -2.77258872223978123767/(WLUT_WIDTH*blur*blur);
1318   for(Q=0; Q<WLUT_WIDTH; Q++)
1319     resample_filter->filter_lut[Q] = exp((double)Q*r_scale);
1320   resample_filter->support = WLUT_WIDTH;
1321   break;
1322 #endif
1323
1324 #if FILTER_LUT
1325 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1326   #pragma omp single
1327 #endif
1328   { register int
1329        Q;
1330     double
1331        r_scale;
1332
1333     /* Scale radius so the filter LUT covers the full support range */
1334     r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
1335     if (IsMagickTrue(GetImageArtifact(resample_filter->image,"resample:verbose")) )
1336       {
1337         /* Debug output of the filter weighting LUT
1338           Gnuplot the LUT with hoizontal adjusted to 'r' using...
1339             plot [0:2][-.2:1] "lut.dat" using (sqrt($0/1024)*2):1 with lines
1340           The filter values is normalized for comparision
1341         */
1342         printf("#\n");
1343         printf("# Resampling Filter LUT (%d values)\n", WLUT_WIDTH);
1344         printf("#\n");
1345         printf("# Note: values in table are using a squared radius lookup.\n");
1346         printf("# And the whole table represents the filters support.\n");
1347         printf("\n"); /* generates a 'break' in gnuplot if multiple outputs */
1348         for(Q=0; Q<WLUT_WIDTH; Q++)
1349           printf("%8.*g %.*g\n",
1350                GetMagickPrecision(),sqrt((double)Q)*r_scale,
1351                GetMagickPrecision(),resample_filter->filter_lut[Q] );
1352       }
1353     /* output the above once only for each image, and each setting */
1354     (void) DeleteImageArtifact(resample_filter->image,"resample:verbose");
1355   }
1356 #endif /* FILTER_LUT */
1357   return;
1358 }
1359 \f
1360 /*
1361 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1362 %                                                                             %
1363 %                                                                             %
1364 %                                                                             %
1365 %   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       %
1366 %                                                                             %
1367 %                                                                             %
1368 %                                                                             %
1369 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1370 %
1371 %  SetResampleFilterInterpolateMethod() sets the resample filter interpolation
1372 %  method.
1373 %
1374 %  The format of the SetResampleFilterInterpolateMethod method is:
1375 %
1376 %      MagickBooleanType SetResampleFilterInterpolateMethod(
1377 %        ResampleFilter *resample_filter,const InterpolateMethod method)
1378 %
1379 %  A description of each parameter follows:
1380 %
1381 %    o resample_filter: the resample filter.
1382 %
1383 %    o method: the interpolation method.
1384 %
1385 */
1386 MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
1387   ResampleFilter *resample_filter,const InterpolatePixelMethod method)
1388 {
1389   assert(resample_filter != (ResampleFilter *) NULL);
1390   assert(resample_filter->signature == MagickSignature);
1391   assert(resample_filter->image != (Image *) NULL);
1392   if (resample_filter->debug != MagickFalse)
1393     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1394       resample_filter->image->filename);
1395   resample_filter->interpolate=method;
1396   return(MagickTrue);
1397 }
1398 \f
1399 /*
1400 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1401 %                                                                             %
1402 %                                                                             %
1403 %                                                                             %
1404 %   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     %
1405 %                                                                             %
1406 %                                                                             %
1407 %                                                                             %
1408 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1409 %
1410 %  SetResampleFilterVirtualPixelMethod() changes the virtual pixel method
1411 %  associated with the specified resample filter.
1412 %
1413 %  The format of the SetResampleFilterVirtualPixelMethod method is:
1414 %
1415 %      MagickBooleanType SetResampleFilterVirtualPixelMethod(
1416 %        ResampleFilter *resample_filter,const VirtualPixelMethod method)
1417 %
1418 %  A description of each parameter follows:
1419 %
1420 %    o resample_filter: the resample filter.
1421 %
1422 %    o method: the virtual pixel method.
1423 %
1424 */
1425 MagickExport MagickBooleanType SetResampleFilterVirtualPixelMethod(
1426   ResampleFilter *resample_filter,const VirtualPixelMethod method)
1427 {
1428   assert(resample_filter != (ResampleFilter *) NULL);
1429   assert(resample_filter->signature == MagickSignature);
1430   assert(resample_filter->image != (Image *) NULL);
1431   if (resample_filter->debug != MagickFalse)
1432     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
1433       resample_filter->image->filename);
1434   resample_filter->virtual_pixel=method;
1435   if (method != UndefinedVirtualPixelMethod)
1436     (void) SetCacheViewVirtualPixelMethod(resample_filter->view,method);
1437   return(MagickTrue);
1438 }