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