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