MagickCore  6.7.5
resample.c
Go to the documentation of this file.
00001 /*
00002 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00003 %                                                                             %
00004 %                                                                             %
00005 %                                                                             %
00006 %           RRRR    EEEEE   SSSSS   AAA   M   M  PPPP   L      EEEEE          %
00007 %           R   R   E       SS     A   A  MM MM  P   P  L      E              %
00008 %           RRRR    EEE      SSS   AAAAA  M M M  PPPP   L      EEE            %
00009 %           R R     E          SS  A   A  M   M  P      L      E              %
00010 %           R  R    EEEEE   SSSSS  A   A  M   M  P      LLLLL  EEEEE          %
00011 %                                                                             %
00012 %                                                                             %
00013 %                      MagickCore Pixel Resampling Methods                    %
00014 %                                                                             %
00015 %                              Software Design                                %
00016 %                                John Cristy                                  %
00017 %                              Anthony Thyssen                                %
00018 %                                August 2007                                  %
00019 %                                                                             %
00020 %                                                                             %
00021 %  Copyright 1999-2012 ImageMagick Studio LLC, a non-profit organization      %
00022 %  dedicated to making software imaging solutions freely available.           %
00023 %                                                                             %
00024 %  You may not use this file except in compliance with the License.  You may  %
00025 %  obtain a copy of the License at                                            %
00026 %                                                                             %
00027 %    http://www.imagemagick.org/script/license.php                            %
00028 %                                                                             %
00029 %  Unless required by applicable law or agreed to in writing, software        %
00030 %  distributed under the License is distributed on an "AS IS" BASIS,          %
00031 %  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.   %
00032 %  See the License for the specific language governing permissions and        %
00033 %  limitations under the License.                                             %
00034 %                                                                             %
00035 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00036 %
00037 %
00038 */
00039 
00040 /*
00041   Include declarations.
00042 */
00043 #include "MagickCore/studio.h"
00044 #include "MagickCore/artifact.h"
00045 #include "MagickCore/color-private.h"
00046 #include "MagickCore/cache.h"
00047 #include "MagickCore/draw.h"
00048 #include "MagickCore/exception-private.h"
00049 #include "MagickCore/gem.h"
00050 #include "MagickCore/image.h"
00051 #include "MagickCore/image-private.h"
00052 #include "MagickCore/log.h"
00053 #include "MagickCore/magick.h"
00054 #include "MagickCore/memory_.h"
00055 #include "MagickCore/pixel.h"
00056 #include "MagickCore/pixel-accessor.h"
00057 #include "MagickCore/quantum.h"
00058 #include "MagickCore/random_.h"
00059 #include "MagickCore/resample.h"
00060 #include "MagickCore/resize.h"
00061 #include "MagickCore/resize-private.h"
00062 #include "MagickCore/transform.h"
00063 #include "MagickCore/signature-private.h"
00064 #include "MagickCore/utility.h"
00065 #include "MagickCore/utility-private.h"
00066 /*
00067   EWA Resampling Options
00068 */
00069 
00070 /* select ONE resampling method */
00071 #define EWA 1                 /* Normal EWA handling - raw or clamped */
00072                               /* if 0 then use "High Quality EWA" */
00073 #define EWA_CLAMP 1           /* EWA Clamping from Nicolas Robidoux */
00074 
00075 #define FILTER_LUT 1          /* Use a LUT rather then direct filter calls */
00076 
00077 /* output debugging information */
00078 #define DEBUG_ELLIPSE 0       /* output ellipse info for debug */
00079 #define DEBUG_HIT_MISS 0      /* output hit/miss pixels (as gnuplot commands) */
00080 #define DEBUG_NO_PIXEL_HIT 0  /* Make pixels that fail to hit anything - RED */
00081 
00082 #if ! FILTER_DIRECT
00083 #define WLUT_WIDTH 1024       /* size of the filter cache */
00084 #endif
00085 
00086 /*
00087   Typedef declarations.
00088 */
00089 struct _ResampleFilter
00090 {
00091   CacheView
00092     *view;
00093 
00094   Image
00095     *image;
00096 
00097   ExceptionInfo
00098     *exception;
00099 
00100   MagickBooleanType
00101     debug;
00102 
00103   /* Information about image being resampled */
00104   ssize_t
00105     image_area;
00106 
00107   PixelInterpolateMethod
00108     interpolate;
00109 
00110   VirtualPixelMethod
00111     virtual_pixel;
00112 
00113   FilterTypes
00114     filter;
00115 
00116   /* processing settings needed */
00117   MagickBooleanType
00118     limit_reached,
00119     do_interpolate,
00120     average_defined;
00121 
00122   PixelInfo
00123     average_pixel;
00124 
00125   /* current ellipitical area being resampled around center point */
00126   double
00127     A, B, C,
00128     Vlimit, Ulimit, Uwidth, slope;
00129 
00130 #if FILTER_LUT
00131   /* LUT of weights for filtered average in elliptical area */
00132   double
00133     filter_lut[WLUT_WIDTH];
00134 #else
00135   /* Use a Direct call to the filter functions */
00136   ResizeFilter
00137     *filter_def;
00138 
00139   double
00140     F;
00141 #endif
00142 
00143   /* the practical working support of the filter */
00144   double
00145     support;
00146 
00147   size_t
00148     signature;
00149 };
00150 
00151 /*
00152 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00153 %                                                                             %
00154 %                                                                             %
00155 %                                                                             %
00156 %   A c q u i r e R e s a m p l e I n f o                                     %
00157 %                                                                             %
00158 %                                                                             %
00159 %                                                                             %
00160 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00161 %
00162 %  AcquireResampleFilter() initializes the information resample needs do to a
00163 %  scaled lookup of a color from an image, using area sampling.
00164 %
00165 %  The algorithm is based on a Elliptical Weighted Average, where the pixels
00166 %  found in a large elliptical area is averaged together according to a
00167 %  weighting (filter) function.  For more details see "Fundamentals of Texture
00168 %  Mapping and Image Warping" a master's thesis by Paul.S.Heckbert, June 17,
00169 %  1989.  Available for free from, http://www.cs.cmu.edu/~ph/
00170 %
00171 %  As EWA resampling (or any sort of resampling) can require a lot of
00172 %  calculations to produce a distorted scaling of the source image for each
00173 %  output pixel, the ResampleFilter structure generated holds that information
00174 %  between individual image resampling.
00175 %
00176 %  This function will make the appropriate AcquireCacheView() calls
00177 %  to view the image, calling functions do not need to open a cache view.
00178 %
00179 %  Usage Example...
00180 %      resample_filter=AcquireResampleFilter(image,exception);
00181 %      SetResampleFilter(resample_filter, GaussianFilter, 1.0);
00182 %      for (y=0; y < (ssize_t) image->rows; y++) {
00183 %        for (x=0; x < (ssize_t) image->columns; x++) {
00184 %          u= ....;   v= ....;
00185 %          ScaleResampleFilter(resample_filter, ... scaling vectors ...);
00186 %          (void) ResamplePixelColor(resample_filter,u,v,&pixel);
00187 %          ... assign resampled pixel value ...
00188 %        }
00189 %      }
00190 %      DestroyResampleFilter(resample_filter);
00191 %
00192 %  The format of the AcquireResampleFilter method is:
00193 %
00194 %     ResampleFilter *AcquireResampleFilter(const Image *image,
00195 %       ExceptionInfo *exception)
00196 %
00197 %  A description of each parameter follows:
00198 %
00199 %    o image: the image.
00200 %
00201 %    o exception: return any errors or warnings in this structure.
00202 %
00203 */
00204 MagickExport ResampleFilter *AcquireResampleFilter(const Image *image,
00205   ExceptionInfo *exception)
00206 {
00207   register ResampleFilter
00208     *resample_filter;
00209 
00210   assert(image != (Image *) NULL);
00211   assert(image->signature == MagickSignature);
00212   if (image->debug != MagickFalse)
00213     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00214   assert(exception != (ExceptionInfo *) NULL);
00215   assert(exception->signature == MagickSignature);
00216 
00217   resample_filter=(ResampleFilter *) AcquireMagickMemory(
00218     sizeof(*resample_filter));
00219   if (resample_filter == (ResampleFilter *) NULL)
00220     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
00221   (void) ResetMagickMemory(resample_filter,0,sizeof(*resample_filter));
00222 
00223   resample_filter->exception=exception;
00224   resample_filter->image=ReferenceImage((Image *) image);
00225   resample_filter->view=AcquireCacheView(resample_filter->image);
00226 
00227   resample_filter->debug=IsEventLogging();
00228   resample_filter->signature=MagickSignature;
00229 
00230   resample_filter->image_area=(ssize_t) (image->columns*image->rows);
00231   resample_filter->average_defined = MagickFalse;
00232 
00233   /* initialise the resampling filter settings */
00234   SetResampleFilter(resample_filter, image->filter, image->blur);
00235   (void) SetResampleFilterInterpolateMethod(resample_filter,
00236     image->interpolate);
00237   (void) SetResampleFilterVirtualPixelMethod(resample_filter,
00238     GetImageVirtualPixelMethod(image));
00239 
00240   return(resample_filter);
00241 }
00242 
00243 /*
00244 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00245 %                                                                             %
00246 %                                                                             %
00247 %                                                                             %
00248 %   D e s t r o y R e s a m p l e I n f o                                     %
00249 %                                                                             %
00250 %                                                                             %
00251 %                                                                             %
00252 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00253 %
00254 %  DestroyResampleFilter() finalizes and cleans up the resampling
00255 %  resample_filter as returned by AcquireResampleFilter(), freeing any memory
00256 %  or other information as needed.
00257 %
00258 %  The format of the DestroyResampleFilter method is:
00259 %
00260 %      ResampleFilter *DestroyResampleFilter(ResampleFilter *resample_filter)
00261 %
00262 %  A description of each parameter follows:
00263 %
00264 %    o resample_filter: resampling information structure
00265 %
00266 */
00267 MagickExport ResampleFilter *DestroyResampleFilter(
00268   ResampleFilter *resample_filter)
00269 {
00270   assert(resample_filter != (ResampleFilter *) NULL);
00271   assert(resample_filter->signature == MagickSignature);
00272   assert(resample_filter->image != (Image *) NULL);
00273   if (resample_filter->debug != MagickFalse)
00274     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
00275       resample_filter->image->filename);
00276   resample_filter->view=DestroyCacheView(resample_filter->view);
00277   resample_filter->image=DestroyImage(resample_filter->image);
00278 #if ! FILTER_LUT
00279   resample_filter->filter_def=DestroyResizeFilter(resample_filter->filter_def);
00280 #endif
00281   resample_filter->signature=(~MagickSignature);
00282   resample_filter=(ResampleFilter *) RelinquishMagickMemory(resample_filter);
00283   return(resample_filter);
00284 }
00285 
00286 /*
00287 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00288 %                                                                             %
00289 %                                                                             %
00290 %                                                                             %
00291 %   R e s a m p l e P i x e l C o l o r                                       %
00292 %                                                                             %
00293 %                                                                             %
00294 %                                                                             %
00295 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00296 %
00297 %  ResamplePixelColor() samples the pixel values surrounding the location
00298 %  given using an elliptical weighted average, at the scale previously
00299 %  calculated, and in the most efficent manner possible for the
00300 %  VirtualPixelMethod setting.
00301 %
00302 %  The format of the ResamplePixelColor method is:
00303 %
00304 %     MagickBooleanType ResamplePixelColor(ResampleFilter *resample_filter,
00305 %       const double u0,const double v0,PixelInfo *pixel)
00306 %
00307 %  A description of each parameter follows:
00308 %
00309 %    o resample_filter: the resample filter.
00310 %
00311 %    o u0,v0: A double representing the center of the area to resample,
00312 %        The distortion transformed transformed x,y coordinate.
00313 %
00314 %    o pixel: the resampled pixel is returned here.
00315 %
00316 */
00317 MagickExport MagickBooleanType ResamplePixelColor(
00318   ResampleFilter *resample_filter,const double u0,const double v0,
00319   PixelInfo *pixel)
00320 {
00321   MagickBooleanType
00322     status;
00323 
00324   ssize_t u,v, v1, v2, uw, hit;
00325   double u1;
00326   double U,V,Q,DQ,DDQ;
00327   double divisor_c,divisor_m;
00328   register double weight;
00329   register const Quantum *pixels;
00330   assert(resample_filter != (ResampleFilter *) NULL);
00331   assert(resample_filter->signature == MagickSignature);
00332 
00333   status=MagickTrue;
00334   /* GetPixelInfo(resample_filter->image,pixel); */
00335   if ( resample_filter->do_interpolate ) {
00336     status=InterpolatePixelInfo(resample_filter->image,resample_filter->view,
00337       resample_filter->interpolate,u0,v0,pixel,resample_filter->exception);
00338     return(status);
00339   }
00340 
00341 #if DEBUG_ELLIPSE
00342   (void) FormatLocaleFile(stderr, "u0=%lf; v0=%lf;\n", u0, v0);
00343 #endif
00344 
00345   /*
00346     Does resample area Miss the image?
00347     And is that area a simple solid color - then return that color
00348   */
00349   hit = 0;
00350   switch ( resample_filter->virtual_pixel ) {
00351     case BackgroundVirtualPixelMethod:
00352     case TransparentVirtualPixelMethod:
00353     case BlackVirtualPixelMethod:
00354     case GrayVirtualPixelMethod:
00355     case WhiteVirtualPixelMethod:
00356     case MaskVirtualPixelMethod:
00357       if ( resample_filter->limit_reached
00358            || u0 + resample_filter->Ulimit < 0.0
00359            || u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
00360            || v0 + resample_filter->Vlimit < 0.0
00361            || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows
00362            )
00363         hit++;
00364       break;
00365 
00366     case UndefinedVirtualPixelMethod:
00367     case EdgeVirtualPixelMethod:
00368       if (    ( u0 + resample_filter->Ulimit < 0.0 && v0 + resample_filter->Vlimit < 0.0 )
00369            || ( u0 + resample_filter->Ulimit < 0.0
00370                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows )
00371            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
00372                 && v0 + resample_filter->Vlimit < 0.0 )
00373            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
00374                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows )
00375            )
00376         hit++;
00377       break;
00378     case HorizontalTileVirtualPixelMethod:
00379       if (    v0 + resample_filter->Vlimit < 0.0
00380            || v0 - resample_filter->Vlimit > (double) resample_filter->image->rows
00381            )
00382         hit++;  /* outside the horizontally tiled images. */
00383       break;
00384     case VerticalTileVirtualPixelMethod:
00385       if (    u0 + resample_filter->Ulimit < 0.0
00386            || u0 - resample_filter->Ulimit > (double) resample_filter->image->columns
00387            )
00388         hit++;  /* outside the vertically tiled images. */
00389       break;
00390     case DitherVirtualPixelMethod:
00391       if (    ( u0 + resample_filter->Ulimit < -32.0 && v0 + resample_filter->Vlimit < -32.0 )
00392            || ( u0 + resample_filter->Ulimit < -32.0
00393                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+32.0 )
00394            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+32.0
00395                 && v0 + resample_filter->Vlimit < -32.0 )
00396            || ( u0 - resample_filter->Ulimit > (double) resample_filter->image->columns+32.0
00397                 && v0 - resample_filter->Vlimit > (double) resample_filter->image->rows+32.0 )
00398            )
00399         hit++;
00400       break;
00401     case TileVirtualPixelMethod:
00402     case MirrorVirtualPixelMethod:
00403     case RandomVirtualPixelMethod:
00404     case HorizontalTileEdgeVirtualPixelMethod:
00405     case VerticalTileEdgeVirtualPixelMethod:
00406     case CheckerTileVirtualPixelMethod:
00407       /* resampling of area is always needed - no VP limits */
00408       break;
00409   }
00410   if ( hit ) {
00411     /* whole area is a solid color -- just return that color */
00412     status=InterpolatePixelInfo(resample_filter->image,
00413       resample_filter->view,IntegerInterpolatePixel,u0,v0,pixel,
00414       resample_filter->exception);
00415     return(status);
00416   }
00417 
00418   /*
00419     Scaling limits reached, return an 'averaged' result.
00420   */
00421   if ( resample_filter->limit_reached ) {
00422     switch ( resample_filter->virtual_pixel ) {
00423       /*  This is always handled by the above, so no need.
00424         case BackgroundVirtualPixelMethod:
00425         case ConstantVirtualPixelMethod:
00426         case TransparentVirtualPixelMethod:
00427         case GrayVirtualPixelMethod,
00428         case WhiteVirtualPixelMethod
00429         case MaskVirtualPixelMethod:
00430       */
00431       case UndefinedVirtualPixelMethod:
00432       case EdgeVirtualPixelMethod:
00433       case DitherVirtualPixelMethod:
00434       case HorizontalTileEdgeVirtualPixelMethod:
00435       case VerticalTileEdgeVirtualPixelMethod:
00436         /* We need an average edge pixel, from the correct edge!
00437            How should I calculate an average edge color?
00438            Just returning an averaged neighbourhood,
00439            works well in general, but falls down for TileEdge methods.
00440            This needs to be done properly!!!!!!
00441         */
00442         status=InterpolatePixelInfo(resample_filter->image,
00443           resample_filter->view,AverageInterpolatePixel,u0,v0,pixel,
00444           resample_filter->exception);
00445         break;
00446       case HorizontalTileVirtualPixelMethod:
00447       case VerticalTileVirtualPixelMethod:
00448         /* just return the background pixel - Is there more direct way? */
00449         status=InterpolatePixelInfo(resample_filter->image,
00450           resample_filter->view,IntegerInterpolatePixel,-1.0,-1.0,pixel,
00451           resample_filter->exception);
00452         break;
00453       case TileVirtualPixelMethod:
00454       case MirrorVirtualPixelMethod:
00455       case RandomVirtualPixelMethod:
00456       case CheckerTileVirtualPixelMethod:
00457       default:
00458         /* generate a average color of the WHOLE image */
00459         if ( resample_filter->average_defined == MagickFalse ) {
00460           Image
00461             *average_image;
00462 
00463           CacheView
00464             *average_view;
00465 
00466           GetPixelInfo(resample_filter->image,(PixelInfo *)
00467             &resample_filter->average_pixel);
00468           resample_filter->average_defined=MagickTrue;
00469 
00470           /* Try to get an averaged pixel color of whole image */
00471           average_image=ResizeImage(resample_filter->image,1,1,BoxFilter,1.0,
00472             resample_filter->exception);
00473           if (average_image == (Image *) NULL)
00474             {
00475               *pixel=resample_filter->average_pixel; /* FAILED */
00476               break;
00477             }
00478           average_view=AcquireCacheView(average_image);
00479           pixels=GetCacheViewVirtualPixels(average_view,0,0,1,1,
00480             resample_filter->exception);
00481           if (pixels == (const Quantum *) NULL) {
00482             average_view=DestroyCacheView(average_view);
00483             average_image=DestroyImage(average_image);
00484             *pixel=resample_filter->average_pixel; /* FAILED */
00485             break;
00486           }
00487           GetPixelInfoPixel(resample_filter->image,pixels,
00488             &(resample_filter->average_pixel));
00489           average_view=DestroyCacheView(average_view);
00490           average_image=DestroyImage(average_image);
00491 
00492           if ( resample_filter->virtual_pixel == CheckerTileVirtualPixelMethod )
00493             {
00494               /* CheckerTile is avergae of image average half background */
00495               /* FUTURE: replace with a 50% blend of both pixels */
00496 
00497               weight = QuantumScale*((MagickRealType)
00498                 resample_filter->average_pixel.alpha);
00499               resample_filter->average_pixel.red *= weight;
00500               resample_filter->average_pixel.green *= weight;
00501               resample_filter->average_pixel.blue *= weight;
00502               divisor_c = weight;
00503 
00504               weight = QuantumScale*((MagickRealType)
00505                 resample_filter->image->background_color.alpha);
00506               resample_filter->average_pixel.red +=
00507                       weight*resample_filter->image->background_color.red;
00508               resample_filter->average_pixel.green +=
00509                       weight*resample_filter->image->background_color.green;
00510               resample_filter->average_pixel.blue +=
00511                       weight*resample_filter->image->background_color.blue;
00512               resample_filter->average_pixel.alpha +=
00513                       resample_filter->image->background_color.alpha;
00514               divisor_c += weight;
00515 
00516               resample_filter->average_pixel.red /= divisor_c;
00517               resample_filter->average_pixel.green /= divisor_c;
00518               resample_filter->average_pixel.blue /= divisor_c;
00519               resample_filter->average_pixel.alpha /= 2;
00520 
00521             }
00522         }
00523         *pixel=resample_filter->average_pixel;
00524         break;
00525     }
00526     return(status);
00527   }
00528 
00529   /*
00530     Initialize weighted average data collection
00531   */
00532   hit = 0;
00533   divisor_c = 0.0;
00534   divisor_m = 0.0;
00535   pixel->red = pixel->green = pixel->blue = 0.0;
00536   if (pixel->colorspace == CMYKColorspace)
00537     pixel->black = 0.0;
00538   if (pixel->matte != MagickFalse)
00539     pixel->alpha = 0.0;
00540 
00541   /*
00542     Determine the parellelogram bounding box fitted to the ellipse
00543     centered at u0,v0.  This area is bounding by the lines...
00544   */
00545   v1 = (ssize_t)ceil(v0 - resample_filter->Vlimit);  /* range of scan lines */
00546   v2 = (ssize_t)floor(v0 + resample_filter->Vlimit);
00547 
00548   /* scan line start and width accross the parallelogram */
00549   u1 = u0 + (v1-v0)*resample_filter->slope - resample_filter->Uwidth;
00550   uw = (ssize_t)(2.0*resample_filter->Uwidth)+1;
00551 
00552 #if DEBUG_ELLIPSE
00553   (void) FormatLocaleFile(stderr, "v1=%ld; v2=%ld\n", (long)v1, (long)v2);
00554   (void) FormatLocaleFile(stderr, "u1=%ld; uw=%ld\n", (long)u1, (long)uw);
00555 #else
00556 # define DEBUG_HIT_MISS 0 /* only valid if DEBUG_ELLIPSE is enabled */
00557 #endif
00558 
00559   /*
00560     Do weighted resampling of all pixels,  within the scaled ellipse,
00561     bound by a Parellelogram fitted to the ellipse.
00562   */
00563   DDQ = 2*resample_filter->A;
00564   for( v=v1; v<=v2;  v++ ) {
00565 #if DEBUG_HIT_MISS
00566     long uu = ceil(u1);   /* actual pixel location (for debug only) */
00567     (void) FormatLocaleFile(stderr, "# scan line from pixel %ld, %ld\n", (long)uu, (long)v);
00568 #endif
00569     u = (ssize_t)ceil(u1);        /* first pixel in scanline */
00570     u1 += resample_filter->slope; /* start of next scan line */
00571 
00572 
00573     /* location of this first pixel, relative to u0,v0 */
00574     U = (double)u-u0;
00575     V = (double)v-v0;
00576 
00577     /* Q = ellipse quotent ( if Q<F then pixel is inside ellipse) */
00578     Q = (resample_filter->A*U + resample_filter->B*V)*U + resample_filter->C*V*V;
00579     DQ = resample_filter->A*(2.0*U+1) + resample_filter->B*V;
00580 
00581     /* get the scanline of pixels for this v */
00582     pixels=GetCacheViewVirtualPixels(resample_filter->view,u,v,(size_t) uw,
00583       1,resample_filter->exception);
00584     if (pixels == (const Quantum *) NULL)
00585       return(MagickFalse);
00586 
00587     /* count up the weighted pixel colors */
00588     for( u=0; u<uw; u++ ) {
00589 #if FILTER_LUT
00590       /* Note that the ellipse has been pre-scaled so F = WLUT_WIDTH */
00591       if ( Q < (double)WLUT_WIDTH ) {
00592         weight = resample_filter->filter_lut[(int)Q];
00593 #else
00594       /* Note that the ellipse has been pre-scaled so F = support^2 */
00595       if ( Q < (double)resample_filter->F ) {
00596         weight = GetResizeFilterWeight(resample_filter->filter_def,
00597              sqrt(Q));    /* a SquareRoot!  Arrggghhhhh... */
00598 #endif
00599 
00600         pixel->alpha  += weight*GetPixelAlpha(resample_filter->image,pixels);
00601         divisor_m += weight;
00602 
00603         if (pixel->matte != MagickFalse)
00604           weight *= QuantumScale*((MagickRealType) GetPixelAlpha(resample_filter->image,pixels));
00605         pixel->red   += weight*GetPixelRed(resample_filter->image,pixels);
00606         pixel->green += weight*GetPixelGreen(resample_filter->image,pixels);
00607         pixel->blue  += weight*GetPixelBlue(resample_filter->image,pixels);
00608         if (pixel->colorspace == CMYKColorspace)
00609           pixel->black += weight*GetPixelBlack(resample_filter->image,pixels);
00610         divisor_c += weight;
00611 
00612         hit++;
00613 #if DEBUG_HIT_MISS
00614         /* mark the pixel according to hit/miss of the ellipse */
00615         (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
00616                      (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
00617         (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 3\n",
00618                      (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
00619       } else {
00620         (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
00621                      (long)uu-.1,(double)v-.1,(long)uu+.1,(long)v+.1);
00622         (void) FormatLocaleFile(stderr, "set arrow from %lf,%lf to %lf,%lf nohead ls 1\n",
00623                      (long)uu+.1,(double)v-.1,(long)uu-.1,(long)v+.1);
00624       }
00625       uu++;
00626 #else
00627       }
00628 #endif
00629       pixels+=GetPixelChannels(resample_filter->image);
00630       Q += DQ;
00631       DQ += DDQ;
00632     }
00633   }
00634 #if DEBUG_ELLIPSE
00635   (void) FormatLocaleFile(stderr, "Hit=%ld;  Total=%ld;\n", (long)hit, (long)uw*(v2-v1) );
00636 #endif
00637 
00638   /*
00639     Result sanity check -- this should NOT happen
00640   */
00641   if ( hit == 0 ) {
00642     /* not enough pixels in resampling, resort to direct interpolation */
00643 #if DEBUG_NO_PIXEL_HIT
00644     pixel->alpha = pixel->red = pixel->green = pixel->blue = 0;
00645     pixel->red = QuantumRange; /* show pixels for which EWA fails */
00646 #else
00647     status=InterpolatePixelInfo(resample_filter->image,
00648       resample_filter->view,resample_filter->interpolate,u0,v0,pixel,
00649       resample_filter->exception);
00650 #endif
00651     return status;
00652   }
00653 
00654   /*
00655     Finialize results of resampling
00656   */
00657   divisor_m = 1.0/divisor_m;
00658   pixel->alpha = (MagickRealType) ClampToQuantum(divisor_m*pixel->alpha);
00659   divisor_c = 1.0/divisor_c;
00660   pixel->red   = (MagickRealType) ClampToQuantum(divisor_c*pixel->red);
00661   pixel->green = (MagickRealType) ClampToQuantum(divisor_c*pixel->green);
00662   pixel->blue  = (MagickRealType) ClampToQuantum(divisor_c*pixel->blue);
00663   if (pixel->colorspace == CMYKColorspace)
00664     pixel->black = (MagickRealType) ClampToQuantum(divisor_c*pixel->black);
00665   return(MagickTrue);
00666 }
00667 
00668 #if EWA && EWA_CLAMP
00669 /*
00670 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00671 %                                                                             %
00672 %                                                                             %
00673 %                                                                             %
00674 -   C l a m p U p A x e s                                                     %
00675 %                                                                             %
00676 %                                                                             %
00677 %                                                                             %
00678 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00679 %
00680 % ClampUpAxes() function converts the input vectors into a major and
00681 % minor axis unit vectors, and their magnitude.  This allows us to
00682 % ensure that the ellipse generated is never smaller than the unit
00683 % circle and thus never too small for use in EWA resampling.
00684 %
00685 % This purely mathematical 'magic' was provided by Professor Nicolas
00686 % Robidoux and his Masters student Chantal Racette.
00687 %
00688 % Reference: "We Recommend Singular Value Decomposition", David Austin
00689 %   http://www.ams.org/samplings/feature-column/fcarc-svd
00690 %
00691 % By generating major and minor axis vectors, we can actually use the
00692 % ellipse in its "canonical form", by remapping the dx,dy of the
00693 % sampled point into distances along the major and minor axis unit
00694 % vectors.
00695 %
00696 % Reference: http://en.wikipedia.org/wiki/Ellipse#Canonical_form
00697 */
00698 static inline void ClampUpAxes(const double dux,
00699                                const double dvx,
00700                                const double duy,
00701                                const double dvy,
00702                                double *major_mag,
00703                                double *minor_mag,
00704                                double *major_unit_x,
00705                                double *major_unit_y,
00706                                double *minor_unit_x,
00707                                double *minor_unit_y)
00708 {
00709   /*
00710    * ClampUpAxes takes an input 2x2 matrix
00711    *
00712    * [ a b ] = [ dux duy ]
00713    * [ c d ] = [ dvx dvy ]
00714    *
00715    * and computes from it the major and minor axis vectors [major_x,
00716    * major_y] and [minor_x,minor_y] of the smallest ellipse containing
00717    * both the unit disk and the ellipse which is the image of the unit
00718    * disk by the linear transformation
00719    *
00720    * [ dux duy ] [S] = [s]
00721    * [ dvx dvy ] [T] = [t]
00722    *
00723    * (The vector [S,T] is the difference between a position in output
00724    * space and [X,Y]; the vector [s,t] is the difference between a
00725    * position in input space and [x,y].)
00726    */
00727   /*
00728    * Output:
00729    *
00730    * major_mag is the half-length of the major axis of the "new"
00731    * ellipse.
00732    *
00733    * minor_mag is the half-length of the minor axis of the "new"
00734    * ellipse.
00735    *
00736    * major_unit_x is the x-coordinate of the major axis direction vector
00737    * of both the "old" and "new" ellipses.
00738    *
00739    * major_unit_y is the y-coordinate of the major axis direction vector.
00740    *
00741    * minor_unit_x is the x-coordinate of the minor axis direction vector.
00742    *
00743    * minor_unit_y is the y-coordinate of the minor axis direction vector.
00744    *
00745    * Unit vectors are useful for computing projections, in particular,
00746    * to compute the distance between a point in output space and the
00747    * center of a unit disk in output space, using the position of the
00748    * corresponding point [s,t] in input space. Following the clamping,
00749    * the square of this distance is
00750    *
00751    * ( ( s * major_unit_x + t * major_unit_y ) / major_mag )^2
00752    * +
00753    * ( ( s * minor_unit_x + t * minor_unit_y ) / minor_mag )^2
00754    *
00755    * If such distances will be computed for many [s,t]'s, it makes
00756    * sense to actually compute the reciprocal of major_mag and
00757    * minor_mag and multiply them by the above unit lengths.
00758    *
00759    * Now, if you want to modify the input pair of tangent vectors so
00760    * that it defines the modified ellipse, all you have to do is set
00761    *
00762    * newdux = major_mag * major_unit_x
00763    * newdvx = major_mag * major_unit_y
00764    * newduy = minor_mag * minor_unit_x = minor_mag * -major_unit_y
00765    * newdvy = minor_mag * minor_unit_y = minor_mag *  major_unit_x
00766    *
00767    * and use these tangent vectors as if they were the original ones.
00768    * Usually, this is a drastic change in the tangent vectors even if
00769    * the singular values are not clamped; for example, the minor axis
00770    * vector always points in a direction which is 90 degrees
00771    * counterclockwise from the direction of the major axis vector.
00772    */
00773   /*
00774    * Discussion:
00775    *
00776    * GOAL: Fix things so that the pullback, in input space, of a disk
00777    * of radius r in output space is an ellipse which contains, at
00778    * least, a disc of radius r. (Make this hold for any r>0.)
00779    *
00780    * ESSENCE OF THE METHOD: Compute the product of the first two
00781    * factors of an SVD of the linear transformation defining the
00782    * ellipse and make sure that both its columns have norm at least 1.
00783    * Because rotations and reflexions map disks to themselves, it is
00784    * not necessary to compute the third (rightmost) factor of the SVD.
00785    *
00786    * DETAILS: Find the singular values and (unit) left singular
00787    * vectors of Jinv, clampling up the singular values to 1, and
00788    * multiply the unit left singular vectors by the new singular
00789    * values in order to get the minor and major ellipse axis vectors.
00790    *
00791    * Image resampling context:
00792    *
00793    * The Jacobian matrix of the transformation at the output point
00794    * under consideration is defined as follows:
00795    *
00796    * Consider the transformation (x,y) -> (X,Y) from input locations
00797    * to output locations. (Anthony Thyssen, elsewhere in resample.c,
00798    * uses the notation (u,v) -> (x,y).)
00799    *
00800    * The Jacobian matrix of the transformation at (x,y) is equal to
00801    *
00802    *   J = [ A, B ] = [ dX/dx, dX/dy ]
00803    *       [ C, D ]   [ dY/dx, dY/dy ]
00804    *
00805    * that is, the vector [A,C] is the tangent vector corresponding to
00806    * input changes in the horizontal direction, and the vector [B,D]
00807    * is the tangent vector corresponding to input changes in the
00808    * vertical direction.
00809    *
00810    * In the context of resampling, it is natural to use the inverse
00811    * Jacobian matrix Jinv because resampling is generally performed by
00812    * pulling pixel locations in the output image back to locations in
00813    * the input image. Jinv is
00814    *
00815    *   Jinv = [ a, b ] = [ dx/dX, dx/dY ]
00816    *          [ c, d ]   [ dy/dX, dy/dY ]
00817    *
00818    * Note: Jinv can be computed from J with the following matrix
00819    * formula:
00820    *
00821    *   Jinv = 1/(A*D-B*C) [  D, -B ]
00822    *                      [ -C,  A ]
00823    *
00824    * What we do is modify Jinv so that it generates an ellipse which
00825    * is as close as possible to the original but which contains the
00826    * unit disk. This can be accomplished as follows:
00827    *
00828    * Let
00829    *
00830    *   Jinv = U Sigma V^T
00831    *
00832    * be an SVD decomposition of Jinv. (The SVD is not unique, but the
00833    * final ellipse does not depend on the particular SVD.)
00834    *
00835    * We could clamp up the entries of the diagonal matrix Sigma so
00836    * that they are at least 1, and then set
00837    *
00838    *   Jinv = U newSigma V^T.
00839    *
00840    * However, we do not need to compute V for the following reason:
00841    * V^T is an orthogonal matrix (that is, it represents a combination
00842    * of rotations and reflexions) so that it maps the unit circle to
00843    * itself. For this reason, the exact value of V does not affect the
00844    * final ellipse, and we can choose V to be the identity
00845    * matrix. This gives
00846    *
00847    *   Jinv = U newSigma.
00848    *
00849    * In the end, we return the two diagonal entries of newSigma
00850    * together with the two columns of U.
00851    */
00852   /*
00853    * ClampUpAxes was written by Nicolas Robidoux and Chantal Racette
00854    * of Laurentian University with insightful suggestions from Anthony
00855    * Thyssen and funding from the National Science and Engineering
00856    * Research Council of Canada. It is distinguished from its
00857    * predecessors by its efficient handling of degenerate cases.
00858    *
00859    * The idea of clamping up the EWA ellipse's major and minor axes so
00860    * that the result contains the reconstruction kernel filter support
00861    * is taken from Andreas Gustaffson's Masters thesis "Interactive
00862    * Image Warping", Helsinki University of Technology, Faculty of
00863    * Information Technology, 59 pages, 1993 (see Section 3.6).
00864    *
00865    * The use of the SVD to clamp up the singular values of the
00866    * Jacobian matrix of the pullback transformation for EWA resampling
00867    * is taken from the astrophysicist Craig DeForest.  It is
00868    * implemented in his PDL::Transform code (PDL = Perl Data
00869    * Language).
00870    */
00871   const double a = dux;
00872   const double b = duy;
00873   const double c = dvx;
00874   const double d = dvy;
00875   /*
00876    * n is the matrix Jinv * transpose(Jinv). Eigenvalues of n are the
00877    * squares of the singular values of Jinv.
00878    */
00879   const double aa = a*a;
00880   const double bb = b*b;
00881   const double cc = c*c;
00882   const double dd = d*d;
00883   /*
00884    * Eigenvectors of n are left singular vectors of Jinv.
00885    */
00886   const double n11 = aa+bb;
00887   const double n12 = a*c+b*d;
00888   const double n21 = n12;
00889   const double n22 = cc+dd;
00890   const double det = a*d-b*c;
00891   const double twice_det = det+det;
00892   const double frobenius_squared = n11+n22;
00893   const double discriminant =
00894     (frobenius_squared+twice_det)*(frobenius_squared-twice_det);
00895   const double sqrt_discriminant = sqrt(discriminant);
00896   /*
00897    * s1 is the largest singular value of the inverse Jacobian
00898    * matrix. In other words, its reciprocal is the smallest singular
00899    * value of the Jacobian matrix itself.
00900    * If s1 = 0, both singular values are 0, and any orthogonal pair of
00901    * left and right factors produces a singular decomposition of Jinv.
00902    */
00903   /*
00904    * Initially, we only compute the squares of the singular values.
00905    */
00906   const double s1s1 = 0.5*(frobenius_squared+sqrt_discriminant);
00907   /*
00908    * s2 the smallest singular value of the inverse Jacobian
00909    * matrix. Its reciprocal is the largest singular value of the
00910    * Jacobian matrix itself.
00911    */
00912   const double s2s2 = 0.5*(frobenius_squared-sqrt_discriminant);
00913   const double s1s1minusn11 = s1s1-n11;
00914   const double s1s1minusn22 = s1s1-n22;
00915   /*
00916    * u1, the first column of the U factor of a singular decomposition
00917    * of Jinv, is a (non-normalized) left singular vector corresponding
00918    * to s1. It has entries u11 and u21. We compute u1 from the fact
00919    * that it is an eigenvector of n corresponding to the eigenvalue
00920    * s1^2.
00921    */
00922   const double s1s1minusn11_squared = s1s1minusn11*s1s1minusn11;
00923   const double s1s1minusn22_squared = s1s1minusn22*s1s1minusn22;
00924   /*
00925    * The following selects the largest row of n-s1^2 I as the one
00926    * which is used to find the eigenvector. If both s1^2-n11 and
00927    * s1^2-n22 are zero, n-s1^2 I is the zero matrix.  In that case,
00928    * any vector is an eigenvector; in addition, norm below is equal to
00929    * zero, and, in exact arithmetic, this is the only case in which
00930    * norm = 0. So, setting u1 to the simple but arbitrary vector [1,0]
00931    * if norm = 0 safely takes care of all cases.
00932    */
00933   const double temp_u11 =
00934     ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? n12 : s1s1minusn22 );
00935   const double temp_u21 =
00936     ( (s1s1minusn11_squared>=s1s1minusn22_squared) ? s1s1minusn11 : n21 );
00937   const double norm = sqrt(temp_u11*temp_u11+temp_u21*temp_u21);
00938   /*
00939    * Finalize the entries of first left singular vector (associated
00940    * with the largest singular value).
00941    */
00942   const double u11 = ( (norm>0.0) ? temp_u11/norm : 1.0 );
00943   const double u21 = ( (norm>0.0) ? temp_u21/norm : 0.0 );
00944   /*
00945    * Clamp the singular values up to 1.
00946    */
00947   *major_mag = ( (s1s1<=1.0) ? 1.0 : sqrt(s1s1) );
00948   *minor_mag = ( (s2s2<=1.0) ? 1.0 : sqrt(s2s2) );
00949   /*
00950    * Return the unit major and minor axis direction vectors.
00951    */
00952   *major_unit_x = u11;
00953   *major_unit_y = u21;
00954   *minor_unit_x = -u21;
00955   *minor_unit_y = u11;
00956 }
00957 
00958 #endif
00959 /*
00960 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00961 %                                                                             %
00962 %                                                                             %
00963 %                                                                             %
00964 %   S c a l e R e s a m p l e F i l t e r                                     %
00965 %                                                                             %
00966 %                                                                             %
00967 %                                                                             %
00968 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00969 %
00970 %  ScaleResampleFilter() does all the calculations needed to resample an image
00971 %  at a specific scale, defined by two scaling vectors.  This not using
00972 %  a orthogonal scaling, but two distorted scaling vectors, to allow the
00973 %  generation of a angled ellipse.
00974 %
00975 %  As only two deritive scaling vectors are used the center of the ellipse
00976 %  must be the center of the lookup.  That is any curvature that the
00977 %  distortion may produce is discounted.
00978 %
00979 %  The input vectors are produced by either finding the derivitives of the
00980 %  distortion function, or the partial derivitives from a distortion mapping.
00981 %  They do not need to be the orthogonal dx,dy scaling vectors, but can be
00982 %  calculated from other derivatives.  For example you could use  dr,da/r
00983 %  polar coordinate vector scaling vectors
00984 %
00985 %  If   u,v =  DistortEquation(x,y)   OR   u = Fu(x,y); v = Fv(x,y)
00986 %  Then the scaling vectors are determined from the deritives...
00987 %      du/dx, dv/dx     and    du/dy, dv/dy
00988 %  If the resulting scaling vectors is othogonally aligned then...
00989 %      dv/dx = 0   and   du/dy  =  0
00990 %  Producing an othogonally alligned ellipse in source space for the area to
00991 %  be resampled.
00992 %
00993 %  Note that scaling vectors are different to argument order.  Argument order
00994 %  is the general order the deritives are extracted from the distortion
00995 %  equations, and not the scaling vectors. As such the middle two vaules
00996 %  may be swapped from what you expect.  Caution is advised.
00997 %
00998 %  WARNING: It is assumed that any SetResampleFilter() method call will
00999 %  always be performed before the ScaleResampleFilter() method, so that the
01000 %  size of the ellipse will match the support for the resampling filter being
01001 %  used.
01002 %
01003 %  The format of the ScaleResampleFilter method is:
01004 %
01005 %     void ScaleResampleFilter(const ResampleFilter *resample_filter,
01006 %       const double dux,const double duy,const double dvx,const double dvy)
01007 %
01008 %  A description of each parameter follows:
01009 %
01010 %    o resample_filter: the resampling resample_filterrmation defining the
01011 %      image being resampled
01012 %
01013 %    o dux,duy,dvx,dvy:
01014 %         The deritives or scaling vectors defining the EWA ellipse.
01015 %         NOTE: watch the order, which is based on the order deritives
01016 %         are usally determined from distortion equations (see above).
01017 %         The middle two values may need to be swapped if you are thinking
01018 %         in terms of scaling vectors.
01019 %
01020 */
01021 MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
01022   const double dux,const double duy,const double dvx,const double dvy)
01023 {
01024   double A,B,C,F;
01025 
01026   assert(resample_filter != (ResampleFilter *) NULL);
01027   assert(resample_filter->signature == MagickSignature);
01028 
01029   resample_filter->limit_reached = MagickFalse;
01030 
01031   /* A 'point' filter forces use of interpolation instead of area sampling */
01032   if ( resample_filter->filter == PointFilter )
01033     return; /* EWA turned off - nothing to do */
01034 
01035 #if DEBUG_ELLIPSE
01036   (void) FormatLocaleFile(stderr, "# -----\n" );
01037   (void) FormatLocaleFile(stderr, "dux=%lf; dvx=%lf;   duy=%lf; dvy=%lf;\n",
01038        dux, dvx, duy, dvy);
01039 #endif
01040 
01041   /* Find Ellipse Coefficents such that
01042         A*u^2 + B*u*v + C*v^2 = F
01043      With u,v relative to point around which we are resampling.
01044      And the given scaling dx,dy vectors in u,v space
01045          du/dx,dv/dx   and  du/dy,dv/dy
01046   */
01047 #if EWA
01048   /* Direct conversion of derivatives into elliptical coefficients
01049      However when magnifying images, the scaling vectors will be small
01050      resulting in a ellipse that is too small to sample properly.
01051      As such we need to clamp the major/minor axis to a minumum of 1.0
01052      to prevent it getting too small.
01053   */
01054 #if EWA_CLAMP
01055   { double major_mag,
01056            minor_mag,
01057            major_x,
01058            major_y,
01059            minor_x,
01060            minor_y;
01061 
01062   ClampUpAxes(dux,dvx,duy,dvy, &major_mag, &minor_mag,
01063                 &major_x, &major_y, &minor_x, &minor_y);
01064   major_x *= major_mag;  major_y *= major_mag;
01065   minor_x *= minor_mag;  minor_y *= minor_mag;
01066 #if DEBUG_ELLIPSE
01067   (void) FormatLocaleFile(stderr, "major_x=%lf; major_y=%lf;  minor_x=%lf; minor_y=%lf;\n",
01068         major_x, major_y, minor_x, minor_y);
01069 #endif
01070   A = major_y*major_y+minor_y*minor_y;
01071   B = -2.0*(major_x*major_y+minor_x*minor_y);
01072   C = major_x*major_x+minor_x*minor_x;
01073   F = major_mag*minor_mag;
01074   F *= F; /* square it */
01075   }
01076 #else /* raw unclamped EWA */
01077   A = dvx*dvx+dvy*dvy;
01078   B = -2.0*(dux*dvx+duy*dvy);
01079   C = dux*dux+duy*duy;
01080   F = dux*dvy-duy*dvx;
01081   F *= F; /* square it */
01082 #endif /* EWA_CLAMP */
01083 
01084 #else /* HQ_EWA */
01085   /*
01086     This Paul Heckbert's "Higher Quality EWA" formula, from page 60 in his
01087     thesis, which adds a unit circle to the elliptical area so as to do both
01088     Reconstruction and Prefiltering of the pixels in the resampling.  It also
01089     means it is always likely to have at least 4 pixels within the area of the
01090     ellipse, for weighted averaging.  No scaling will result with F == 4.0 and
01091     a circle of radius 2.0, and F smaller than this means magnification is
01092     being used.
01093 
01094     NOTE: This method produces a very blury result at near unity scale while
01095     producing perfect results for strong minitification and magnifications.
01096 
01097     However filter support is fixed to 2.0 (no good for Windowed Sinc filters)
01098   */
01099   A = dvx*dvx+dvy*dvy+1;
01100   B = -2.0*(dux*dvx+duy*dvy);
01101   C = dux*dux+duy*duy+1;
01102   F = A*C - B*B/4;
01103 #endif
01104 
01105 #if DEBUG_ELLIPSE
01106   (void) FormatLocaleFile(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
01107 
01108   /* Figure out the various information directly about the ellipse.
01109      This information currently not needed at this time, but may be
01110      needed later for better limit determination.
01111 
01112      It is also good to have as a record for future debugging
01113   */
01114   { double alpha, beta, gamma, Major, Minor;
01115     double Eccentricity, Ellipse_Area, Ellipse_Angle;
01116 
01117     alpha = A+C;
01118     beta  = A-C;
01119     gamma = sqrt(beta*beta + B*B );
01120 
01121     if ( alpha - gamma <= MagickEpsilon )
01122       Major = MagickHuge;
01123     else
01124       Major = sqrt(2*F/(alpha - gamma));
01125     Minor = sqrt(2*F/(alpha + gamma));
01126 
01127     (void) FormatLocaleFile(stderr, "# Major=%lf; Minor=%lf\n", Major, Minor );
01128 
01129     /* other information about ellipse include... */
01130     Eccentricity = Major/Minor;
01131     Ellipse_Area = MagickPI*Major*Minor;
01132     Ellipse_Angle = atan2(B, A-C);
01133 
01134     (void) FormatLocaleFile(stderr, "# Angle=%lf   Area=%lf\n",
01135          RadiansToDegrees(Ellipse_Angle), Ellipse_Area);
01136   }
01137 #endif
01138 
01139   /* If one or both of the scaling vectors is impossibly large
01140      (producing a very large raw F value), we may as well not bother
01141      doing any form of resampling since resampled area is very large.
01142      In this case some alternative means of pixel sampling, such as
01143      the average of the whole image is needed to get a reasonable
01144      result. Calculate only as needed.
01145   */
01146   if ( (4*A*C - B*B) > MagickHuge ) {
01147     resample_filter->limit_reached = MagickTrue;
01148     return;
01149   }
01150 
01151   /* Scale ellipse to match the filters support
01152      (that is, multiply F by the square of the support).
01153   */
01154   F *= resample_filter->support;
01155   F *= resample_filter->support;
01156 
01157   /* Orthogonal bounds of the ellipse */
01158   resample_filter->Ulimit = sqrt(C*F/(A*C-0.25*B*B));
01159   resample_filter->Vlimit = sqrt(A*F/(A*C-0.25*B*B));
01160 
01161   /* Horizontally aligned parallelogram fitted to Ellipse */
01162   resample_filter->Uwidth = sqrt(F/A); /* Half of the parallelogram width */
01163   resample_filter->slope = -B/(2.0*A); /* Reciprocal slope of the parallelogram */
01164 
01165 #if DEBUG_ELLIPSE
01166   (void) FormatLocaleFile(stderr, "Ulimit=%lf; Vlimit=%lf; UWidth=%lf; Slope=%lf;\n",
01167            resample_filter->Ulimit, resample_filter->Vlimit,
01168            resample_filter->Uwidth, resample_filter->slope );
01169 #endif
01170 
01171   /* Check the absolute area of the parallelogram involved.
01172    * This limit needs more work, as it is too slow for larger images
01173    * with tiled views of the horizon.
01174   */
01175   if ( (resample_filter->Uwidth * resample_filter->Vlimit)
01176          > (4.0*resample_filter->image_area)) {
01177     resample_filter->limit_reached = MagickTrue;
01178     return;
01179   }
01180 
01181   /* Scale ellipse formula to directly index the Filter Lookup Table */
01182   { register double scale;
01183 #if FILTER_LUT
01184     /* scale so that F = WLUT_WIDTH; -- hardcoded */
01185     scale = (double)WLUT_WIDTH/F;
01186 #else
01187     /* scale so that F = resample_filter->F (support^2) */
01188     scale = resample_filter->F/F;
01189 #endif
01190     resample_filter->A = A*scale;
01191     resample_filter->B = B*scale;
01192     resample_filter->C = C*scale;
01193   }
01194 }
01195 
01196 /*
01197 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01198 %                                                                             %
01199 %                                                                             %
01200 %                                                                             %
01201 %   S e t R e s a m p l e F i l t e r                                         %
01202 %                                                                             %
01203 %                                                                             %
01204 %                                                                             %
01205 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01206 %
01207 %  SetResampleFilter() set the resampling filter lookup table based on a
01208 %  specific filter.  Note that the filter is used as a radial filter not as a
01209 %  two pass othogonally aligned resampling filter.
01210 %
01211 %  The default Filter, is Gaussian, which is the standard filter used by the
01212 %  original paper on the Elliptical Weighted Everage Algorithm. However other
01213 %  filters can also be used.
01214 %
01215 %  The format of the SetResampleFilter method is:
01216 %
01217 %    void SetResampleFilter(ResampleFilter *resample_filter,
01218 %      const FilterTypes filter,const double blur)
01219 %
01220 %  A description of each parameter follows:
01221 %
01222 %    o resample_filter: resampling resample_filterrmation structure
01223 %
01224 %    o filter: the resize filter for elliptical weighting LUT
01225 %
01226 %    o blur: filter blur factor (radial scaling) for elliptical weighting LUT
01227 %
01228 */
01229 MagickExport void SetResampleFilter(ResampleFilter *resample_filter,
01230   const FilterTypes filter,const double blur)
01231 {
01232   ResizeFilter
01233      *resize_filter;
01234 
01235   assert(resample_filter != (ResampleFilter *) NULL);
01236   assert(resample_filter->signature == MagickSignature);
01237 
01238   resample_filter->do_interpolate = MagickFalse;
01239   resample_filter->filter = filter;
01240 
01241   if ( filter == PointFilter )
01242     {
01243       resample_filter->do_interpolate = MagickTrue;
01244       return;  /* EWA turned off - nothing more to do */
01245     }
01246 
01247   /* Set a default cylindrical filter of a 'low blur' Jinc windowed Jinc */
01248   if ( filter == UndefinedFilter )
01249     resample_filter->filter = RobidouxFilter;
01250 
01251   resize_filter = AcquireResizeFilter(resample_filter->image,
01252        resample_filter->filter,blur,MagickTrue,resample_filter->exception);
01253   if (resize_filter == (ResizeFilter *) NULL)
01254     {
01255       (void) ThrowMagickException(resample_filter->exception,GetMagickModule(),
01256            ModuleError, "UnableToSetFilteringValue",
01257            "Fall back to default EWA gaussian filter");
01258       resample_filter->filter = PointFilter;
01259     }
01260 
01261   /* Get the practical working support for the filter,
01262    * after any API call blur factors have been accoded for.
01263    */
01264 #if EWA
01265   resample_filter->support = GetResizeFilterSupport(resize_filter);
01266 #else
01267   resample_filter->support = 2.0;  /* fixed support size for HQ-EWA */
01268 #endif
01269 
01270 #if FILTER_LUT
01271   /* Fill the LUT with the weights from the selected filter function */
01272   { register int
01273        Q;
01274     double
01275        r_scale;
01276     /* Scale radius so the filter LUT covers the full support range */
01277     r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
01278     for(Q=0; Q<WLUT_WIDTH; Q++)
01279       resample_filter->filter_lut[Q] = (double)
01280            GetResizeFilterWeight(resize_filter,sqrt((double)Q)*r_scale);
01281 
01282     /* finished with the resize filter */
01283     resize_filter = DestroyResizeFilter(resize_filter);
01284   }
01285 #else
01286   /* save the filter and the scaled ellipse bounds needed for filter */
01287   resample_filter->filter_def = resize_filter;
01288   resample_filter->F = resample_filter->support*resample_filter->support;
01289 #endif
01290 
01291   /*
01292     Adjust the scaling of the default unit circle
01293     This assumes that any real scaling changes will always
01294     take place AFTER the filter method has been initialized.
01295   */
01296   ScaleResampleFilter(resample_filter, 1.0, 0.0, 0.0, 1.0);
01297 
01298 #if 0
01299   /* This is old code kept as a reference only.  It is very wrong,
01300      and I don't understand exactly what it was attempting to do.
01301   */
01302   /*
01303     Create Normal Gaussian 2D Filter Weighted Lookup Table.
01304     A normal EWA guassual lookup would use   exp(Q*ALPHA)
01305     where  Q = distance squared from 0.0 (center) to 1.0 (edge)
01306     and    ALPHA = -4.0*ln(2.0)  ==>  -2.77258872223978123767
01307     The table is of length 1024, and equates to support radius of 2.0
01308     thus needs to be scaled by  ALPHA*4/1024 and any blur factor squared
01309 
01310     The above came from some reference code provided by Fred Weinhaus
01311     and seems to have been a guess that was appropriate for its use
01312     in a 3d perspective landscape mapping program.
01313   */
01314   r_scale = -2.77258872223978123767/(WLUT_WIDTH*blur*blur);
01315   for(Q=0; Q<WLUT_WIDTH; Q++)
01316     resample_filter->filter_lut[Q] = exp((double)Q*r_scale);
01317   resample_filter->support = WLUT_WIDTH;
01318   break;
01319 #endif
01320 
01321 #if FILTER_LUT
01322 #if defined(MAGICKCORE_OPENMP_SUPPORT)
01323   #pragma omp single
01324 #endif
01325   { register int
01326        Q;
01327     double
01328        r_scale;
01329 
01330     /* Scale radius so the filter LUT covers the full support range */
01331     r_scale = resample_filter->support*sqrt(1.0/(double)WLUT_WIDTH);
01332     if (IsMagickTrue(GetImageArtifact(resample_filter->image,"resample:verbose")) )
01333       {
01334         /* Debug output of the filter weighting LUT
01335           Gnuplot the LUT with hoizontal adjusted to 'r' using...
01336             plot [0:2][-.2:1] "lut.dat" using (sqrt($0/1024)*2):1 with lines
01337           The filter values is normalized for comparision
01338         */
01339         printf("#\n");
01340         printf("# Resampling Filter LUT (%d values)\n", WLUT_WIDTH);
01341         printf("#\n");
01342         printf("# Note: values in table are using a squared radius lookup.\n");
01343         printf("# And the whole table represents the filters support.\n");
01344         printf("\n"); /* generates a 'break' in gnuplot if multiple outputs */
01345         for(Q=0; Q<WLUT_WIDTH; Q++)
01346           printf("%8.*g %.*g\n",
01347                GetMagickPrecision(),sqrt((double)Q)*r_scale,
01348                GetMagickPrecision(),resample_filter->filter_lut[Q] );
01349       }
01350     /* output the above once only for each image, and each setting */
01351     (void) DeleteImageArtifact(resample_filter->image,"resample:verbose");
01352   }
01353 #endif /* FILTER_LUT */
01354   return;
01355 }
01356 
01357 /*
01358 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01359 %                                                                             %
01360 %                                                                             %
01361 %                                                                             %
01362 %   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       %
01363 %                                                                             %
01364 %                                                                             %
01365 %                                                                             %
01366 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01367 %
01368 %  SetResampleFilterInterpolateMethod() sets the resample filter interpolation
01369 %  method.
01370 %
01371 %  The format of the SetResampleFilterInterpolateMethod method is:
01372 %
01373 %      MagickBooleanType SetResampleFilterInterpolateMethod(
01374 %        ResampleFilter *resample_filter,const InterpolateMethod method)
01375 %
01376 %  A description of each parameter follows:
01377 %
01378 %    o resample_filter: the resample filter.
01379 %
01380 %    o method: the interpolation method.
01381 %
01382 */
01383 MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
01384   ResampleFilter *resample_filter,const PixelInterpolateMethod method)
01385 {
01386   assert(resample_filter != (ResampleFilter *) NULL);
01387   assert(resample_filter->signature == MagickSignature);
01388   assert(resample_filter->image != (Image *) NULL);
01389   if (resample_filter->debug != MagickFalse)
01390     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
01391       resample_filter->image->filename);
01392   resample_filter->interpolate=method;
01393   return(MagickTrue);
01394 }
01395 
01396 /*
01397 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01398 %                                                                             %
01399 %                                                                             %
01400 %                                                                             %
01401 %   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     %
01402 %                                                                             %
01403 %                                                                             %
01404 %                                                                             %
01405 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01406 %
01407 %  SetResampleFilterVirtualPixelMethod() changes the virtual pixel method
01408 %  associated with the specified resample filter.
01409 %
01410 %  The format of the SetResampleFilterVirtualPixelMethod method is:
01411 %
01412 %      MagickBooleanType SetResampleFilterVirtualPixelMethod(
01413 %        ResampleFilter *resample_filter,const VirtualPixelMethod method)
01414 %
01415 %  A description of each parameter follows:
01416 %
01417 %    o resample_filter: the resample filter.
01418 %
01419 %    o method: the virtual pixel method.
01420 %
01421 */
01422 MagickExport MagickBooleanType SetResampleFilterVirtualPixelMethod(
01423   ResampleFilter *resample_filter,const VirtualPixelMethod method)
01424 {
01425   assert(resample_filter != (ResampleFilter *) NULL);
01426   assert(resample_filter->signature == MagickSignature);
01427   assert(resample_filter->image != (Image *) NULL);
01428   if (resample_filter->debug != MagickFalse)
01429     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
01430       resample_filter->image->filename);
01431   resample_filter->virtual_pixel=method;
01432   if (method != UndefinedVirtualPixelMethod)
01433     (void) SetCacheViewVirtualPixelMethod(resample_filter->view,method);
01434   return(MagickTrue);
01435 }