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