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-2009 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 "magick/studio.h"
00044 #include "magick/artifact.h"
00045 #include "magick/color-private.h"
00046 #include "magick/cache.h"
00047 #include "magick/draw.h"
00048 #include "magick/exception-private.h"
00049 #include "magick/gem.h"
00050 #include "magick/image.h"
00051 #include "magick/image-private.h"
00052 #include "magick/log.h"
00053 #include "magick/memory_.h"
00054 #include "magick/pixel-private.h"
00055 #include "magick/quantum.h"
00056 #include "magick/random_.h"
00057 #include "magick/resample.h"
00058 #include "magick/resize.h"
00059 #include "magick/resize-private.h"
00060 #include "magick/transform.h"
00061 #include "magick/signature-private.h"
00062 /*
00063   Typedef declarations.
00064 */
00065 #define WLUT_WIDTH 1024
00066 struct _ResampleFilter
00067 {
00068   Image
00069     *image;
00070 
00071   CacheView
00072     *view;
00073 
00074   ExceptionInfo
00075     *exception;
00076 
00077   MagickBooleanType
00078     debug;
00079 
00080   /* Information about image being resampled */
00081   long
00082     image_area;
00083 
00084   InterpolatePixelMethod
00085     interpolate;
00086 
00087   VirtualPixelMethod
00088     virtual_pixel;
00089 
00090   FilterTypes
00091     filter;
00092 
00093   /* processing settings needed */
00094   MagickBooleanType
00095     limit_reached,
00096     do_interpolate,
00097     average_defined;
00098 
00099   MagickPixelPacket
00100     average_pixel;
00101 
00102   /* current ellipitical area being resampled around center point */
00103   double
00104     A, B, C,
00105     sqrtA, sqrtC, sqrtU, slope;
00106 
00107   /* LUT of weights for filtered average in elliptical area */
00108   double
00109     filter_lut[WLUT_WIDTH],
00110     support;
00111 
00112   unsigned long
00113     signature;
00114 };
00115 
00116 /*
00117 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00118 %                                                                             %
00119 %                                                                             %
00120 %                                                                             %
00121 %   A c q u i r e R e s a m p l e I n f o                                     %
00122 %                                                                             %
00123 %                                                                             %
00124 %                                                                             %
00125 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00126 %
00127 %  AcquireResampleFilter() initializes the information resample needs do to a
00128 %  scaled lookup of a color from an image, using area sampling.
00129 %
00130 %  The algorithm is based on a Elliptical Weighted Average, where the pixels
00131 %  found in a large elliptical area is averaged together according to a
00132 %  weighting (filter) function.  For more details see "Fundamentals of Texture
00133 %  Mapping and Image Warping" a master's thesis by Paul.S.Heckbert, June 17,
00134 %  1989.  Available for free from, http://www.cs.cmu.edu/~ph/
00135 %
00136 %  As EWA resampling (or any sort of resampling) can require a lot of
00137 %  calculations to produce a distorted scaling of the source image for each
00138 %  output pixel, the ResampleFilter structure generated holds that information
00139 %  between individual image resampling.
00140 %
00141 %  This function will make the appropriate AcquireCacheView() calls
00142 %  to view the image, calling functions do not need to open a cache view.
00143 %
00144 %  Usage Example...
00145 %      resample_filter=AcquireResampleFilter(image,exception);
00146 %      for (y=0; y < (long) image->rows; y++) {
00147 %        for (x=0; x < (long) image->columns; x++) {
00148 %          X= ....;   Y= ....;
00149 %          ScaleResampleFilter(resample_filter, ... scaling vectors ...);
00150 %          (void) ResamplePixelColor(resample_filter,X,Y,&pixel);
00151 %          ... assign resampled pixel value ...
00152 %        }
00153 %      }
00154 %      DestroyResampleFilter(resample_filter);
00155 %
00156 %  The format of the AcquireResampleFilter method is:
00157 %
00158 %     ResampleFilter *AcquireResampleFilter(const Image *image,
00159 %       ExceptionInfo *exception)
00160 %
00161 %  A description of each parameter follows:
00162 %
00163 %    o image: the image.
00164 %
00165 %    o exception: return any errors or warnings in this structure.
00166 %
00167 */
00168 MagickExport ResampleFilter *AcquireResampleFilter(const Image *image,
00169   ExceptionInfo *exception)
00170 {
00171   register ResampleFilter
00172     *resample_filter;
00173 
00174   assert(image != (Image *) NULL);
00175   assert(image->signature == MagickSignature);
00176   if (image->debug != MagickFalse)
00177     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
00178   assert(exception != (ExceptionInfo *) NULL);
00179   assert(exception->signature == MagickSignature);
00180 
00181   resample_filter=(ResampleFilter *) AcquireMagickMemory(
00182     sizeof(*resample_filter));
00183   if (resample_filter == (ResampleFilter *) NULL)
00184     ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
00185   (void) ResetMagickMemory(resample_filter,0,sizeof(*resample_filter));
00186 
00187   resample_filter->image=ReferenceImage((Image *) image);
00188   resample_filter->view=AcquireCacheView(resample_filter->image);
00189   resample_filter->exception=exception;
00190 
00191   resample_filter->debug=IsEventLogging();
00192   resample_filter->signature=MagickSignature;
00193 
00194   resample_filter->image_area = (long) resample_filter->image->columns *
00195     resample_filter->image->rows;
00196   resample_filter->average_defined = MagickFalse;
00197 
00198   /* initialise the resampling filter settings */
00199   SetResampleFilter(resample_filter, resample_filter->image->filter,
00200     resample_filter->image->blur);
00201   resample_filter->interpolate = resample_filter->image->interpolate;
00202   resample_filter->virtual_pixel=GetImageVirtualPixelMethod(image);
00203 
00204   /* init scale to a default of a unit circle */
00205   ScaleResampleFilter(resample_filter, 1.0, 0.0, 0.0, 1.0);
00206 
00207   return(resample_filter);
00208 }
00209 
00210 /*
00211 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00212 %                                                                             %
00213 %                                                                             %
00214 %                                                                             %
00215 %   D e s t r o y R e s a m p l e I n f o                                     %
00216 %                                                                             %
00217 %                                                                             %
00218 %                                                                             %
00219 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00220 %
00221 %  DestroyResampleFilter() finalizes and cleans up the resampling
00222 %  resample_filter as returned by AcquireResampleFilter(), freeing any memory
00223 %  or other information as needed.
00224 %
00225 %  The format of the DestroyResampleFilter method is:
00226 %
00227 %      ResampleFilter *DestroyResampleFilter(ResampleFilter *resample_filter)
00228 %
00229 %  A description of each parameter follows:
00230 %
00231 %    o resample_filter: resampling information structure
00232 %
00233 */
00234 MagickExport ResampleFilter *DestroyResampleFilter(
00235   ResampleFilter *resample_filter)
00236 {
00237   assert(resample_filter != (ResampleFilter *) NULL);
00238   assert(resample_filter->signature == MagickSignature);
00239   assert(resample_filter->image != (Image *) NULL);
00240   if (resample_filter->debug != MagickFalse)
00241     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
00242       resample_filter->image->filename);
00243   resample_filter->view=DestroyCacheView(resample_filter->view);
00244   resample_filter->image=DestroyImage(resample_filter->image);
00245   resample_filter->signature=(~MagickSignature);
00246   resample_filter=(ResampleFilter *) RelinquishMagickMemory(resample_filter);
00247   return(resample_filter);
00248 }
00249 
00250 /*
00251 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00252 %                                                                             %
00253 %                                                                             %
00254 %                                                                             %
00255 %   I n t e r p o l a t e R e s a m p l e F i l t e r                         %
00256 %                                                                             %
00257 %                                                                             %
00258 %                                                                             %
00259 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00260 %
00261 %  InterpolateResampleFilter() applies bi-linear or tri-linear interpolation
00262 %  between a floating point coordinate and the pixels surrounding that
00263 %  coordinate.  No pixel area resampling, or scaling of the result is
00264 %  performed.
00265 %
00266 %  The format of the InterpolateResampleFilter method is:
00267 %
00268 %      MagickBooleanType InterpolateResampleFilter(
00269 %        ResampleInfo *resample_filter,const InterpolatePixelMethod method,
00270 %        const double x,const double y,MagickPixelPacket *pixel)
00271 %
00272 %  A description of each parameter follows:
00273 %
00274 %    o resample_filter: the resample filter.
00275 %
00276 %    o method: the pixel clor interpolation method.
00277 %
00278 %    o x,y: A double representing the current (x,y) position of the pixel.
00279 %
00280 %    o pixel: return the interpolated pixel here.
00281 %
00282 */
00283 
00284 static inline double MagickMax(const double x,const double y)
00285 {
00286   if (x > y)
00287     return(x);
00288   return(y);
00289 }
00290 
00291 static void BicubicInterpolate(const MagickPixelPacket *pixels,const double dx,
00292   MagickPixelPacket *pixel)
00293 {
00294   MagickRealType
00295     dx2,
00296     p,
00297     q,
00298     r,
00299     s;
00300 
00301   dx2=dx*dx;
00302   p=(pixels[3].red-pixels[2].red)-(pixels[0].red-pixels[1].red);
00303   q=(pixels[0].red-pixels[1].red)-p;
00304   r=pixels[2].red-pixels[0].red;
00305   s=pixels[1].red;
00306   pixel->red=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
00307   p=(pixels[3].green-pixels[2].green)-(pixels[0].green-pixels[1].green);
00308   q=(pixels[0].green-pixels[1].green)-p;
00309   r=pixels[2].green-pixels[0].green;
00310   s=pixels[1].green;
00311   pixel->green=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
00312   p=(pixels[3].blue-pixels[2].blue)-(pixels[0].blue-pixels[1].blue);
00313   q=(pixels[0].blue-pixels[1].blue)-p;
00314   r=pixels[2].blue-pixels[0].blue;
00315   s=pixels[1].blue;
00316   pixel->blue=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
00317   p=(pixels[3].opacity-pixels[2].opacity)-(pixels[0].opacity-pixels[1].opacity);
00318   q=(pixels[0].opacity-pixels[1].opacity)-p;
00319   r=pixels[2].opacity-pixels[0].opacity;
00320   s=pixels[1].opacity;
00321   pixel->opacity=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
00322   if (pixel->colorspace == CMYKColorspace)
00323     {
00324       p=(pixels[3].index-pixels[2].index)-(pixels[0].index-pixels[1].index);
00325       q=(pixels[0].index-pixels[1].index)-p;
00326       r=pixels[2].index-pixels[0].index;
00327       s=pixels[1].index;
00328       pixel->index=(dx*dx2*p)+(dx2*q)+(dx*r)+s;
00329     }
00330 }
00331 
00332 static inline MagickRealType CubicWeightingFunction(const MagickRealType x)
00333 {
00334   MagickRealType
00335     alpha,
00336     gamma;
00337 
00338   alpha=MagickMax(x+2.0,0.0);
00339   gamma=1.0*alpha*alpha*alpha;
00340   alpha=MagickMax(x+1.0,0.0);
00341   gamma-=4.0*alpha*alpha*alpha;
00342   alpha=MagickMax(x+0.0,0.0);
00343   gamma+=6.0*alpha*alpha*alpha;
00344   alpha=MagickMax(x-1.0,0.0);
00345   gamma-=4.0*alpha*alpha*alpha;
00346   return(gamma/6.0);
00347 }
00348 
00349 static inline double MeshInterpolate(const PointInfo *delta,const double p,
00350   const double x,const double y)
00351 {
00352   return(delta->x*x+delta->y*y+(1.0-delta->x-delta->y)*p);
00353 }
00354 
00355 static inline long NearestNeighbor(MagickRealType x)
00356 {
00357   if (x >= 0.0)
00358     return((long) (x+0.5));
00359   return((long) (x-0.5));
00360 }
00361 
00362 static MagickBooleanType InterpolateResampleFilter(
00363   ResampleFilter *resample_filter,const InterpolatePixelMethod method,
00364   const double x,const double y,MagickPixelPacket *pixel)
00365 {
00366   MagickBooleanType
00367     status;
00368 
00369   register const IndexPacket
00370     *indexes;
00371 
00372   register const PixelPacket
00373     *p;
00374 
00375   register long
00376     i;
00377 
00378   assert(resample_filter != (ResampleFilter *) NULL);
00379   assert(resample_filter->signature == MagickSignature);
00380   status=MagickTrue;
00381   switch (method)
00382   {
00383     case AverageInterpolatePixel:
00384     {
00385       MagickPixelPacket
00386         pixels[16];
00387 
00388       MagickRealType
00389         alpha[16],
00390         gamma;
00391 
00392       p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x)-1,(long)
00393         floor(y)-1,4,4,resample_filter->exception);
00394       if (p == (const PixelPacket *) NULL)
00395         {
00396           status=MagickFalse;
00397           break;
00398         }
00399       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
00400       for (i=0; i < 16L; i++)
00401       {
00402         GetMagickPixelPacket(resample_filter->image,pixels+i);
00403         SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
00404         alpha[i]=1.0;
00405         if (resample_filter->image->matte != MagickFalse)
00406           {
00407             alpha[i]=QuantumScale*((MagickRealType) QuantumRange-p->opacity);
00408             pixels[i].red*=alpha[i];
00409             pixels[i].green*=alpha[i];
00410             pixels[i].blue*=alpha[i];
00411             if (resample_filter->image->colorspace == CMYKColorspace)
00412               pixels[i].index*=alpha[i];
00413           }
00414         gamma=alpha[i];
00415         gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
00416         pixel->red+=gamma*0.0625*pixels[i].red;
00417         pixel->green+=gamma*0.0625*pixels[i].green;
00418         pixel->blue+=gamma*0.0625*pixels[i].blue;
00419         pixel->opacity+=0.0625*pixels[i].opacity;
00420         if (resample_filter->image->colorspace == CMYKColorspace)
00421           pixel->index+=gamma*0.0625*pixels[i].index;
00422         p++;
00423       }
00424       break;
00425     }
00426     case BicubicInterpolatePixel:
00427     {
00428       MagickPixelPacket
00429         pixels[16],
00430         u[4];
00431 
00432       MagickRealType
00433         alpha[16];
00434 
00435       PointInfo
00436         delta;
00437 
00438       p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x)-1,(long)
00439         floor(y)-1,4,4,resample_filter->exception);
00440       if (p == (const PixelPacket *) NULL)
00441         {
00442           status=MagickFalse;
00443           break;
00444         }
00445       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
00446       for (i=0; i < 16L; i++)
00447       {
00448         GetMagickPixelPacket(resample_filter->image,pixels+i);
00449         SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
00450         alpha[i]=1.0;
00451         if (resample_filter->image->matte != MagickFalse)
00452           {
00453             alpha[i]=QuantumScale*((MagickRealType) QuantumRange-p->opacity);
00454             pixels[i].red*=alpha[i];
00455             pixels[i].green*=alpha[i];
00456             pixels[i].blue*=alpha[i];
00457             if (resample_filter->image->colorspace == CMYKColorspace)
00458               pixels[i].index*=alpha[i];
00459           }
00460         p++;
00461       }
00462       delta.x=x-floor(x);
00463       for (i=0; i < 4L; i++)
00464         BicubicInterpolate(pixels+4*i,delta.x,u+i);
00465       delta.y=y-floor(y);
00466       BicubicInterpolate(u,delta.y,pixel);
00467       break;
00468     }
00469     case BilinearInterpolatePixel:
00470     default:
00471     {
00472       MagickPixelPacket
00473         pixels[4];
00474 
00475       MagickRealType
00476         alpha[4],
00477         gamma;
00478 
00479       PointInfo
00480         delta,
00481         epsilon;
00482 
00483       p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x),(long)
00484         floor(y),2,2,resample_filter->exception);
00485       if (p == (const PixelPacket *) NULL)
00486         {
00487           status=MagickFalse;
00488           break;
00489         }
00490       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
00491       for (i=0; i < 4L; i++)
00492       {
00493         pixels[i].red=(MagickRealType) p[i].red;
00494         pixels[i].green=(MagickRealType) p[i].green;
00495         pixels[i].blue=(MagickRealType) p[i].blue;
00496         pixels[i].opacity=(MagickRealType) p[i].opacity;
00497         alpha[i]=1.0;
00498       }
00499       if (resample_filter->image->matte != MagickFalse)
00500         for (i=0; i < 4L; i++)
00501         {
00502           alpha[i]=QuantumScale*((MagickRealType) QuantumRange-p[i].opacity);
00503           pixels[i].red*=alpha[i];
00504           pixels[i].green*=alpha[i];
00505           pixels[i].blue*=alpha[i];
00506         }
00507       if (indexes != (IndexPacket *) NULL)
00508         for (i=0; i < 4L; i++)
00509         {
00510           pixels[i].index=(MagickRealType) indexes[i];
00511           if (resample_filter->image->colorspace == CMYKColorspace)
00512             pixels[i].index*=alpha[i];
00513         }
00514       delta.x=x-floor(x);
00515       delta.y=y-floor(y);
00516       epsilon.x=1.0-delta.x;
00517       epsilon.y=1.0-delta.y;
00518       gamma=((epsilon.y*(epsilon.x*alpha[0]+delta.x*alpha[1])+delta.y*
00519         (epsilon.x*alpha[2]+delta.x*alpha[3])));
00520       gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
00521       pixel->red=gamma*(epsilon.y*(epsilon.x*pixels[0].red+delta.x*
00522         pixels[1].red)+delta.y*(epsilon.x*pixels[2].red+delta.x*pixels[3].red));
00523       pixel->green=gamma*(epsilon.y*(epsilon.x*pixels[0].green+delta.x*
00524         pixels[1].green)+delta.y*(epsilon.x*pixels[2].green+delta.x*
00525         pixels[3].green));
00526       pixel->blue=gamma*(epsilon.y*(epsilon.x*pixels[0].blue+delta.x*
00527         pixels[1].blue)+delta.y*(epsilon.x*pixels[2].blue+delta.x*
00528         pixels[3].blue));
00529       pixel->opacity=(epsilon.y*(epsilon.x*pixels[0].opacity+delta.x*
00530         pixels[1].opacity)+delta.y*(epsilon.x*pixels[2].opacity+delta.x*
00531         pixels[3].opacity));
00532       if (resample_filter->image->colorspace == CMYKColorspace)
00533         pixel->index=gamma*(epsilon.y*(epsilon.x*pixels[0].index+delta.x*
00534           pixels[1].index)+delta.y*(epsilon.x*pixels[2].index+delta.x*
00535           pixels[3].index));
00536       break;
00537     }
00538     case FilterInterpolatePixel:
00539     {
00540       Image
00541         *excerpt_image,
00542         *filter_image;
00543 
00544       MagickPixelPacket
00545         pixels[1];
00546 
00547       RectangleInfo
00548         geometry;
00549 
00550       CacheView
00551         *filter_view;
00552 
00553       geometry.width=4L;
00554       geometry.height=4L;
00555       geometry.x=(long) floor(x)-1L;
00556       geometry.y=(long) floor(y)-1L;
00557       excerpt_image=ExcerptImage(resample_filter->image,&geometry,
00558         resample_filter->exception);
00559       if (excerpt_image == (Image *) NULL)
00560         {
00561           status=MagickFalse;
00562           break;
00563         }
00564       filter_image=ResizeImage(excerpt_image,1,1,resample_filter->image->filter,
00565         resample_filter->image->blur,resample_filter->exception);
00566       excerpt_image=DestroyImage(excerpt_image);
00567       if (filter_image == (Image *) NULL)
00568         break;
00569       filter_view=AcquireCacheView(filter_image);
00570       p=GetCacheViewVirtualPixels(filter_view,0,0,1,1,
00571         resample_filter->exception);
00572       if (p != (const PixelPacket *) NULL)
00573         {
00574           indexes=GetVirtualIndexQueue(filter_image);
00575           GetMagickPixelPacket(resample_filter->image,pixels);
00576           SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
00577         }
00578       filter_view=DestroyCacheView(filter_view);
00579       filter_image=DestroyImage(filter_image);
00580       break;
00581     }
00582     case IntegerInterpolatePixel:
00583     {
00584       MagickPixelPacket
00585         pixels[1];
00586 
00587       p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x),(long)
00588         floor(y),1,1,resample_filter->exception);
00589       if (p == (const PixelPacket *) NULL)
00590         {
00591           status=MagickFalse;
00592           break;
00593         }
00594       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
00595       GetMagickPixelPacket(resample_filter->image,pixels);
00596       SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
00597       break;
00598     }
00599     case MeshInterpolatePixel:
00600     {
00601       MagickPixelPacket
00602         pixels[4];
00603 
00604       MagickRealType
00605         alpha[4],
00606         gamma;
00607 
00608       PointInfo
00609         delta,
00610         luminance;
00611 
00612       p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x),(long)
00613         floor(y),2,2,resample_filter->exception);
00614       if (p == (const PixelPacket *) NULL)
00615         {
00616           status=MagickFalse;
00617           break;
00618         }
00619       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
00620       for (i=0; i < 4L; i++)
00621       {
00622         GetMagickPixelPacket(resample_filter->image,pixels+i);
00623         SetMagickPixelPacket(resample_filter->image,p,indexes+i,pixels+i);
00624         alpha[i]=1.0;
00625         if (resample_filter->image->matte != MagickFalse)
00626           {
00627             alpha[i]=QuantumScale*((MagickRealType) QuantumRange-p->opacity);
00628             pixels[i].red*=alpha[i];
00629             pixels[i].green*=alpha[i];
00630             pixels[i].blue*=alpha[i];
00631             if (resample_filter->image->colorspace == CMYKColorspace)
00632               pixels[i].index*=alpha[i];
00633           }
00634         p++;
00635       }
00636       delta.x=x-floor(x);
00637       delta.y=y-floor(y);
00638       luminance.x=MagickPixelLuminance(pixels+0)-MagickPixelLuminance(pixels+3);
00639       luminance.y=MagickPixelLuminance(pixels+1)-MagickPixelLuminance(pixels+2);
00640       if (fabs(luminance.x) < fabs(luminance.y))
00641         {
00642           /*
00643             Diagonal 0-3 NW-SE.
00644           */
00645           if (delta.x <= delta.y)
00646             {
00647               /*
00648                 Bottom-left triangle  (pixel:2, diagonal: 0-3).
00649               */
00650               delta.y=1.0-delta.y;
00651               gamma=MeshInterpolate(&delta,alpha[2],alpha[3],alpha[0]);
00652               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
00653               pixel->red=gamma*MeshInterpolate(&delta,pixels[2].red,
00654                 pixels[3].red,pixels[0].red);
00655               pixel->green=gamma*MeshInterpolate(&delta,pixels[2].green,
00656                 pixels[3].green,pixels[0].green);
00657               pixel->blue=gamma*MeshInterpolate(&delta,pixels[2].blue,
00658                 pixels[3].blue,pixels[0].blue);
00659               pixel->opacity=gamma*MeshInterpolate(&delta,pixels[2].opacity,
00660                 pixels[3].opacity,pixels[0].opacity);
00661               if (resample_filter->image->colorspace == CMYKColorspace)
00662                 pixel->index=gamma*MeshInterpolate(&delta,pixels[2].index,
00663                   pixels[3].index,pixels[0].index);
00664             }
00665           else
00666             {
00667               /*
00668                 Top-right triangle (pixel:1, diagonal: 0-3).
00669               */
00670               delta.x=1.0-delta.x;
00671               gamma=MeshInterpolate(&delta,alpha[1],alpha[0],alpha[3]);
00672               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
00673               pixel->red=gamma*MeshInterpolate(&delta,pixels[1].red,
00674                 pixels[0].red,pixels[3].red);
00675               pixel->green=gamma*MeshInterpolate(&delta,pixels[1].green,
00676                 pixels[0].green,pixels[3].green);
00677               pixel->blue=gamma*MeshInterpolate(&delta,pixels[1].blue,
00678                 pixels[0].blue,pixels[3].blue);
00679               pixel->opacity=gamma*MeshInterpolate(&delta,pixels[1].opacity,
00680                 pixels[0].opacity,pixels[3].opacity);
00681               if (resample_filter->image->colorspace == CMYKColorspace)
00682                 pixel->index=gamma*MeshInterpolate(&delta,pixels[1].index,
00683                   pixels[0].index,pixels[3].index);
00684             }
00685         }
00686       else
00687         {
00688           /*
00689             Diagonal 1-2 NE-SW.
00690           */
00691           if (delta.x <= (1.0-delta.y))
00692             {
00693               /*
00694                 Top-left triangle (pixel 0, diagonal: 1-2).
00695               */
00696               gamma=MeshInterpolate(&delta,alpha[0],alpha[1],alpha[2]);
00697               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
00698               pixel->red=gamma*MeshInterpolate(&delta,pixels[0].red,
00699                 pixels[1].red,pixels[2].red);
00700               pixel->green=gamma*MeshInterpolate(&delta,pixels[0].green,
00701                 pixels[1].green,pixels[2].green);
00702               pixel->blue=gamma*MeshInterpolate(&delta,pixels[0].blue,
00703                 pixels[1].blue,pixels[2].blue);
00704               pixel->opacity=gamma*MeshInterpolate(&delta,pixels[0].opacity,
00705                 pixels[1].opacity,pixels[2].opacity);
00706               if (resample_filter->image->colorspace == CMYKColorspace)
00707                 pixel->index=gamma*MeshInterpolate(&delta,pixels[0].index,
00708                   pixels[1].index,pixels[2].index);
00709             }
00710           else
00711             {
00712               /*
00713                 Bottom-right triangle (pixel: 3, diagonal: 1-2).
00714               */
00715               delta.x=1.0-delta.x;
00716               delta.y=1.0-delta.y;
00717               gamma=MeshInterpolate(&delta,alpha[3],alpha[2],alpha[1]);
00718               gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
00719               pixel->red=gamma*MeshInterpolate(&delta,pixels[3].red,
00720                 pixels[2].red,pixels[1].red);
00721               pixel->green=gamma*MeshInterpolate(&delta,pixels[3].green,
00722                 pixels[2].green,pixels[1].green);
00723               pixel->blue=gamma*MeshInterpolate(&delta,pixels[3].blue,
00724                 pixels[2].blue,pixels[1].blue);
00725               pixel->opacity=gamma*MeshInterpolate(&delta,pixels[3].opacity,
00726                 pixels[2].opacity,pixels[1].opacity);
00727               if (resample_filter->image->colorspace == CMYKColorspace)
00728                 pixel->index=gamma*MeshInterpolate(&delta,pixels[3].index,
00729                   pixels[2].index,pixels[1].index);
00730             }
00731         }
00732       break;
00733     }
00734     case NearestNeighborInterpolatePixel:
00735     {
00736       MagickPixelPacket
00737         pixels[1];
00738 
00739       p=GetCacheViewVirtualPixels(resample_filter->view,NearestNeighbor(x),
00740         NearestNeighbor(y),1,1,resample_filter->exception);
00741       if (p == (const PixelPacket *) NULL)
00742         {
00743           status=MagickFalse;
00744           break;
00745         }
00746       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
00747       GetMagickPixelPacket(resample_filter->image,pixels);
00748       SetMagickPixelPacket(resample_filter->image,p,indexes,pixel);
00749       break;
00750     }
00751     case SplineInterpolatePixel:
00752     {
00753       long
00754         j,
00755         n;
00756 
00757       MagickPixelPacket
00758         pixels[16];
00759 
00760       MagickRealType
00761         alpha[16],
00762         dx,
00763         dy,
00764         gamma;
00765 
00766       PointInfo
00767         delta;
00768 
00769       p=GetCacheViewVirtualPixels(resample_filter->view,(long) floor(x)-1,(long)
00770         floor(y)-1,4,4,resample_filter->exception);
00771       if (p == (const PixelPacket *) NULL)
00772         {
00773           status=MagickFalse;
00774           break;
00775         }
00776       indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
00777       n=0;
00778       delta.x=x-floor(x);
00779       delta.y=y-floor(y);
00780       for (i=(-1); i < 3L; i++)
00781       {
00782         dy=CubicWeightingFunction((MagickRealType) i-delta.y);
00783         for (j=(-1); j < 3L; j++)
00784         {
00785           GetMagickPixelPacket(resample_filter->image,pixels+n);
00786           SetMagickPixelPacket(resample_filter->image,p,indexes+n,pixels+n);
00787           alpha[n]=1.0;
00788           if (resample_filter->image->matte != MagickFalse)
00789             {
00790               alpha[n]=QuantumScale*((MagickRealType) QuantumRange-p->opacity);
00791               pixels[n].red*=alpha[n];
00792               pixels[n].green*=alpha[n];
00793               pixels[n].blue*=alpha[n];
00794               if (resample_filter->image->colorspace == CMYKColorspace)
00795                 pixels[n].index*=alpha[n];
00796             }
00797           dx=CubicWeightingFunction(delta.x-(MagickRealType) j);
00798           gamma=alpha[n];
00799           gamma=1.0/(fabs((double) gamma) <= MagickEpsilon ? 1.0 : gamma);
00800           pixel->red+=gamma*dx*dy*pixels[n].red;
00801           pixel->green+=gamma*dx*dy*pixels[n].green;
00802           pixel->blue+=gamma*dx*dy*pixels[n].blue;
00803           if (resample_filter->image->matte != MagickFalse)
00804             pixel->opacity+=dx*dy*pixels[n].opacity;
00805           if (resample_filter->image->colorspace == CMYKColorspace)
00806             pixel->index+=gamma*dx*dy*pixels[n].index;
00807           n++;
00808           p++;
00809         }
00810       }
00811       break;
00812     }
00813   }
00814   return(status);
00815 }
00816 
00817 /*
00818 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00819 %                                                                             %
00820 %                                                                             %
00821 %                                                                             %
00822 %   R e s a m p l e P i x e l C o l o r                                       %
00823 %                                                                             %
00824 %                                                                             %
00825 %                                                                             %
00826 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
00827 %
00828 %  ResamplePixelColor() samples the pixel values surrounding the location
00829 %  given using an elliptical weighted average, at the scale previously
00830 %  calculated, and in the most efficent manner possible for the
00831 %  VirtualPixelMethod setting.
00832 %
00833 %  The format of the ResamplePixelColor method is:
00834 %
00835 %     MagickBooleanType ResamplePixelColor(ResampleFilter *resample_filter,
00836 %       const double u0,const double v0,MagickPixelPacket *pixel)
00837 %
00838 %  A description of each parameter follows:
00839 %
00840 %    o resample_filter: the resample filter.
00841 %
00842 %    o u0,v0: A double representing the center of the area to resample,
00843 %        The distortion transformed transformed x,y coordinate.
00844 %
00845 %    o pixel: the resampled pixel is returned here.
00846 %
00847 */
00848 MagickExport MagickBooleanType ResamplePixelColor(
00849   ResampleFilter *resample_filter,const double u0,const double v0,
00850   MagickPixelPacket *pixel)
00851 {
00852   MagickBooleanType
00853     status;
00854 
00855   long u,v, uw,v1,v2, hit;
00856   double u1;
00857   double U,V,Q,DQ,DDQ;
00858   double divisor_c,divisor_m;
00859   register double weight;
00860   register const PixelPacket *pixels;
00861   register const IndexPacket *indexes;
00862   assert(resample_filter != (ResampleFilter *) NULL);
00863   assert(resample_filter->signature == MagickSignature);
00864 
00865   status=MagickTrue;
00866   GetMagickPixelPacket(resample_filter->image,pixel);
00867   if ( resample_filter->do_interpolate ) {
00868     status=InterpolateResampleFilter(resample_filter,
00869       resample_filter->interpolate,u0,v0,pixel);
00870     return(status);
00871   }
00872 
00873   /*
00874     Does resample area Miss the image?
00875     And is that area a simple solid color - then return that color
00876   */
00877   hit = 0;
00878   switch ( resample_filter->virtual_pixel ) {
00879     case BackgroundVirtualPixelMethod:
00880     case ConstantVirtualPixelMethod:
00881     case TransparentVirtualPixelMethod:
00882     case BlackVirtualPixelMethod:
00883     case GrayVirtualPixelMethod:
00884     case WhiteVirtualPixelMethod:
00885     case MaskVirtualPixelMethod:
00886       if ( resample_filter->limit_reached
00887            || u0 + resample_filter->sqrtC < 0.0
00888            || u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
00889            || v0 + resample_filter->sqrtA < 0.0
00890            || v0 - resample_filter->sqrtA > (double) resample_filter->image->rows
00891            )
00892         hit++;
00893       break;
00894 
00895     case UndefinedVirtualPixelMethod:
00896     case EdgeVirtualPixelMethod:
00897       if (    ( u0 + resample_filter->sqrtC < 0.0 && v0 + resample_filter->sqrtA < 0.0 )
00898            || ( u0 + resample_filter->sqrtC < 0.0
00899                 && v0 - resample_filter->sqrtA > (double) resample_filter->image->rows )
00900            || ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
00901                 && v0 + resample_filter->sqrtA < 0.0 )
00902            || ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
00903                 && v0 - resample_filter->sqrtA > (double) resample_filter->image->rows )
00904            )
00905         hit++;
00906       break;
00907     case HorizontalTileVirtualPixelMethod:
00908       if (    v0 + resample_filter->sqrtA < 0.0
00909            || v0 - resample_filter->sqrtA > (double) resample_filter->image->rows
00910            )
00911         hit++;  /* outside the horizontally tiled images. */
00912       break;
00913     case VerticalTileVirtualPixelMethod:
00914       if (    u0 + resample_filter->sqrtC < 0.0
00915            || u0 - resample_filter->sqrtC > (double) resample_filter->image->columns
00916            )
00917         hit++;  /* outside the vertically tiled images. */
00918       break;
00919     case DitherVirtualPixelMethod:
00920       if (    ( u0 + resample_filter->sqrtC < -32.0 && v0 + resample_filter->sqrtA < -32.0 )
00921            || ( u0 + resample_filter->sqrtC < -32.0
00922                 && v0 - resample_filter->sqrtA > (double) resample_filter->image->rows+32.0 )
00923            || ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns+32.0
00924                 && v0 + resample_filter->sqrtA < -32.0 )
00925            || ( u0 - resample_filter->sqrtC > (double) resample_filter->image->columns+32.0
00926                 && v0 - resample_filter->sqrtA > (double) resample_filter->image->rows+32.0 )
00927            )
00928         hit++;
00929       break;
00930     case TileVirtualPixelMethod:
00931     case MirrorVirtualPixelMethod:
00932     case RandomVirtualPixelMethod:
00933     case HorizontalTileEdgeVirtualPixelMethod:
00934     case VerticalTileEdgeVirtualPixelMethod:
00935     case CheckerTileVirtualPixelMethod:
00936       /* resampling of area is always needed - no VP limits */
00937       break;
00938   }
00939   if ( hit ) {
00940     /* whole area is a solid color -- just return that color */
00941     status=InterpolateResampleFilter(resample_filter,IntegerInterpolatePixel,
00942       u0,v0,pixel);
00943     return(status);
00944   }
00945 
00946   /*
00947     Scaling limits reached, return an 'averaged' result.
00948   */
00949   if ( resample_filter->limit_reached ) {
00950     switch ( resample_filter->virtual_pixel ) {
00951       /*  This is always handled by the above, so no need.
00952         case BackgroundVirtualPixelMethod:
00953         case ConstantVirtualPixelMethod:
00954         case TransparentVirtualPixelMethod:
00955         case GrayVirtualPixelMethod,
00956         case WhiteVirtualPixelMethod
00957         case MaskVirtualPixelMethod:
00958       */
00959       case UndefinedVirtualPixelMethod:
00960       case EdgeVirtualPixelMethod:
00961       case DitherVirtualPixelMethod:
00962       case HorizontalTileEdgeVirtualPixelMethod:
00963       case VerticalTileEdgeVirtualPixelMethod:
00964         /* We need an average edge pixel, for the right edge!
00965            How should I calculate an average edge color?
00966            Just returning an averaged neighbourhood,
00967            works well in general, but falls down for TileEdge methods.
00968            This needs to be done properly!!!!!!
00969         */
00970         status=InterpolateResampleFilter(resample_filter,
00971           AverageInterpolatePixel,u0,v0,pixel);
00972         break;
00973       case HorizontalTileVirtualPixelMethod:
00974       case VerticalTileVirtualPixelMethod:
00975         /* just return the background pixel - Is there more direct way? */
00976         status=InterpolateResampleFilter(resample_filter,
00977            IntegerInterpolatePixel,(double)-1,(double)-1,pixel);
00978         break;
00979       case TileVirtualPixelMethod:
00980       case MirrorVirtualPixelMethod:
00981       case RandomVirtualPixelMethod:
00982       case CheckerTileVirtualPixelMethod:
00983       default:
00984         /* generate a average color of the WHOLE image */
00985         if ( resample_filter->average_defined == MagickFalse ) {
00986           Image
00987             *average_image;
00988 
00989           CacheView
00990             *average_view;
00991 
00992           GetMagickPixelPacket(resample_filter->image,
00993                 (MagickPixelPacket *)&(resample_filter->average_pixel));
00994           resample_filter->average_defined = MagickTrue;
00995 
00996           /* Try to get an averaged pixel color of whole image */
00997           average_image=ResizeImage(resample_filter->image,1,1,BoxFilter,1.0,
00998            resample_filter->exception);
00999           if (average_image == (Image *) NULL)
01000             {
01001               *pixel=resample_filter->average_pixel; /* FAILED */
01002               break;
01003             }
01004           average_view=AcquireCacheView(average_image);
01005           pixels=(PixelPacket *)GetCacheViewVirtualPixels(average_view,0,0,1,1,
01006             resample_filter->exception);
01007           if (pixels == (const PixelPacket *) NULL) {
01008             average_view=DestroyCacheView(average_view);
01009             average_image=DestroyImage(average_image);
01010             *pixel=resample_filter->average_pixel; /* FAILED */
01011             break;
01012           }
01013           indexes=(IndexPacket *) GetCacheViewAuthenticIndexQueue(average_view);
01014           SetMagickPixelPacket(resample_filter->image,pixels,indexes,
01015             &(resample_filter->average_pixel));
01016           average_view=DestroyCacheView(average_view);
01017           average_image=DestroyImage(average_image);
01018 #if 0
01019           /* CheckerTile should average the image with background color */
01020           //if ( resample_filter->virtual_pixel == CheckerTileVirtualPixelMethod ) {
01021 #if 0
01022             resample_filter->average_pixel.red =
01023                       ( resample_filter->average_pixel.red +
01024                           resample_filter->image->background_color.red ) /2;
01025             resample_filter->average_pixel.green =
01026                       ( resample_filter->average_pixel.green +
01027                           resample_filter->image->background_color.green ) /2;
01028             resample_filter->average_pixel.blue =
01029                       ( resample_filter->average_pixel.blue +
01030                           resample_filter->image->background_color.blue ) /2;
01031             resample_filter->average_pixel.matte =
01032                       ( resample_filter->average_pixel.matte +
01033                           resample_filter->image->background_color.matte ) /2;
01034             resample_filter->average_pixel.black =
01035                       ( resample_filter->average_pixel.black +
01036                           resample_filter->image->background_color.black ) /2;
01037 #else
01038             resample_filter->average_pixel =
01039                           resample_filter->image->background_color;
01040 #endif
01041           }
01042 #endif
01043         }
01044         *pixel=resample_filter->average_pixel;
01045         break;
01046     }
01047     return(status);
01048   }
01049 
01050   /*
01051     Initialize weighted average data collection
01052   */
01053   hit = 0;
01054   divisor_c = 0.0;
01055   divisor_m = 0.0;
01056   pixel->red = pixel->green = pixel->blue = 0.0;
01057   if (resample_filter->image->matte != MagickFalse) pixel->opacity = 0.0;
01058   if (resample_filter->image->colorspace == CMYKColorspace) pixel->index = 0.0;
01059 
01060   /*
01061     Determine the parellelogram bounding box fitted to the ellipse
01062     centered at u0,v0.  This area is bounding by the lines...
01063         v = +/- sqrt(A)
01064         u = -By/2A  +/- sqrt(F/A)
01065     Which has been pre-calculated above.
01066   */
01067   v1 = (long)(v0 - resample_filter->sqrtA);               /* range of scan lines */
01068   v2 = (long)(v0 + resample_filter->sqrtA + 1);
01069 
01070   u1 = u0 + (v1-v0)*resample_filter->slope - resample_filter->sqrtU; /* start of scanline for v=v1 */
01071   uw = (long)(2*resample_filter->sqrtU)+1;       /* width of parallelogram */
01072 
01073   /*
01074     Do weighted resampling of all pixels,  within the scaled ellipse,
01075     bound by a Parellelogram fitted to the ellipse.
01076   */
01077   DDQ = 2*resample_filter->A;
01078   for( v=v1; v<=v2;  v++, u1+=resample_filter->slope ) {
01079     u = (long)u1;       /* first pixel in scanline  ( floor(u1) ) */
01080     U = (double)u-u0;   /* location of that pixel, relative to u0,v0 */
01081     V = (double)v-v0;
01082 
01083     /* Q = ellipse quotent ( if Q<F then pixel is inside ellipse) */
01084     Q = U*(resample_filter->A*U + resample_filter->B*V) + resample_filter->C*V*V;
01085     DQ = resample_filter->A*(2.0*U+1) + resample_filter->B*V;
01086 
01087     /* get the scanline of pixels for this v */
01088     pixels=GetCacheViewVirtualPixels(resample_filter->view,u,v,(unsigned long) uw,
01089       1,resample_filter->exception);
01090     if (pixels == (const PixelPacket *) NULL)
01091       return(MagickFalse);
01092     indexes=GetCacheViewVirtualIndexQueue(resample_filter->view);
01093 
01094     /* count up the weighted pixel colors */
01095     for( u=0; u<uw; u++ ) {
01096       /* Note that the ellipse has been pre-scaled so F = WLUT_WIDTH */
01097       if ( Q < (double)WLUT_WIDTH ) {
01098         weight = resample_filter->filter_lut[(int)Q];
01099 
01100         pixel->opacity  += weight*pixels->opacity;
01101         divisor_m += weight;
01102 
01103         if (resample_filter->image->matte != MagickFalse)
01104           weight *= QuantumScale*((MagickRealType)(QuantumRange-pixels->opacity));
01105         pixel->red   += weight*pixels->red;
01106         pixel->green += weight*pixels->green;
01107         pixel->blue  += weight*pixels->blue;
01108         if (resample_filter->image->colorspace == CMYKColorspace)
01109           pixel->index += weight*(*indexes);
01110         divisor_c += weight;
01111 
01112         hit++;
01113       }
01114       pixels++;
01115       indexes++;
01116       Q += DQ;
01117       DQ += DDQ;
01118     }
01119   }
01120 
01121   /*
01122     Result sanity check -- this should NOT happen
01123   */
01124   if ( hit < 4 || divisor_c < 1.0 ) {
01125     /* not enough pixels in resampling, resort to direct interpolation */
01126     status=InterpolateResampleFilter(resample_filter,
01127       resample_filter->interpolate,u0,v0,pixel);
01128     return status;
01129   }
01130 
01131   /*
01132     Finialize results of resampling
01133   */
01134   divisor_m = 1.0/divisor_m;
01135   pixel->opacity = (MagickRealType) RoundToQuantum(divisor_m*pixel->opacity);
01136   divisor_c = 1.0/divisor_c;
01137   pixel->red   = (MagickRealType) RoundToQuantum(divisor_c*pixel->red);
01138   pixel->green = (MagickRealType) RoundToQuantum(divisor_c*pixel->green);
01139   pixel->blue  = (MagickRealType) RoundToQuantum(divisor_c*pixel->blue);
01140   if (resample_filter->image->colorspace == CMYKColorspace)
01141     pixel->index = (MagickRealType) RoundToQuantum(divisor_c*pixel->index);
01142   return(MagickTrue);
01143 }
01144 
01145 /*
01146 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01147 %                                                                             %
01148 %                                                                             %
01149 %                                                                             %
01150 %   S c a l e R e s a m p l e F i l t e r                                     %
01151 %                                                                             %
01152 %                                                                             %
01153 %                                                                             %
01154 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01155 %
01156 %  ScaleResampleFilter() does all the calculations needed to resample an image
01157 %  at a specific scale, defined by two scaling vectors.  This not using
01158 %  a orthogonal scaling, but two distorted scaling vectors, to allow the
01159 %  generation of a angled ellipse.
01160 %
01161 %  As only two deritive scaling vectors are used the center of the ellipse
01162 %  must be the center of the lookup.  That is any curvature that the
01163 %  distortion may produce is discounted.
01164 %
01165 %  The input vectors are produced by either finding the derivitives of the
01166 %  distortion function, or the partial derivitives from a distortion mapping.
01167 %  They do not need to be the orthogonal dx,dy scaling vectors, but can be
01168 %  calculated from other derivatives.  For example you could use  dr,da/r
01169 %  polar coordinate vector scaling vectors
01170 %
01171 %  If   u,v =  DistortEquation(x,y)
01172 %  Then the scaling vectors dx,dy (in u,v space) are the derivitives...
01173 %      du/dx, dv/dx     and    du/dy, dv/dy
01174 %  If the scaling is only othogonally aligned then...
01175 %      dv/dx = 0   and   du/dy  =  0
01176 %  Producing an othogonally alligned ellipse for the area to be resampled.
01177 %
01178 %  Note that scaling vectors are different to argument order.  Argument order
01179 %  is the general order the deritives are extracted from the distortion
01180 %  equations, EG: U(x,y), V(x,y).  Caution is advised if you are trying to
01181 %  define the ellipse directly from scaling vectors.
01182 %
01183 %  The format of the ScaleResampleFilter method is:
01184 %
01185 %     void ScaleResampleFilter(const ResampleFilter *resample_filter,
01186 %       const double dux,const double duy,const double dvx,const double dvy)
01187 %
01188 %  A description of each parameter follows:
01189 %
01190 %    o resample_filter: the resampling resample_filterrmation defining the
01191 %      image being resampled
01192 %
01193 %    o dux,duy,dvx,dvy:
01194 %         The partial derivitives or scaling vectors for resampling.
01195 %           dx = du/dx, dv/dx    and  dy = du/dy, dv/dy
01196 %
01197 %         The values are used to define the size and angle of the
01198 %         elliptical resampling area, centered on the lookup point.
01199 %
01200 */
01201 MagickExport void ScaleResampleFilter(ResampleFilter *resample_filter,
01202   const double dux,const double duy,const double dvx,const double dvy)
01203 {
01204   double A,B,C,F, area;
01205 
01206   assert(resample_filter != (ResampleFilter *) NULL);
01207   assert(resample_filter->signature == MagickSignature);
01208 
01209   resample_filter->limit_reached = MagickFalse;
01210   resample_filter->do_interpolate = MagickFalse;
01211 
01212   /* A 'point' filter forces use of interpolation instead of area sampling */
01213   if ( resample_filter->filter == PointFilter ) {
01214     resample_filter->do_interpolate = MagickTrue;
01215     return;
01216   }
01217 
01218   /* Find Ellipse Coefficents such that
01219         A*u^2 + B*u*v + C*v^2 = F
01220      With u,v relative to point around which we are resampling.
01221      And the given scaling dx,dy vectors in u,v space
01222          du/dx,dv/dx   and  du/dy,dv/dy
01223   */
01224 #if 0
01225   /* Direct conversions of derivatives to elliptical coefficients
01226      No scaling will result in F == 1.0 and a unit circle.
01227   */
01228   A = dvx*dvx+dvy*dvy;
01229   B = (dux*dvx+duy*dvy)*-2.0;
01230   C = dux*dux+duy*duy;
01231   F = dux*dvy+duy*dvx;
01232   F *= F;
01233 #define F_UNITY 1.0
01234 #else
01235   /* This Paul Heckbert's recomended "Higher Quality EWA" formula, from page
01236      60 in his thesis, which adds a unit circle to the elliptical area so are
01237      to do both Reconstruction and Prefiltering of the pixels in the
01238      resampling.  It also means it is likely to have at least 4 pixels within
01239      the area of the ellipse, for weighted averaging.
01240      No scaling will result if F == 4.0 and a circle of radius 2.0
01241   */
01242   A = dvx*dvx+dvy*dvy+1;
01243   B = (dux*dvx+duy*dvy)*-2.0;
01244   C = dux*dux+duy*duy+1;
01245   F = A*C - B*B/4;
01246 #define F_UNITY 4.0
01247 #endif
01248 
01249 /* DEBUGGING OUTPUT */
01250 #if 0
01251   fprintf(stderr, "dux=%lf; dvx=%lf;   duy=%lf; dvy%lf;\n",
01252        dux, dvx, duy, dvy);
01253   fprintf(stderr, "A=%lf; B=%lf; C=%lf; F=%lf\n", A,B,C,F);
01254 #endif
01255 
01256 #if 0
01257   /* Figure out the Ellipses Major and Minor Axis, and other info.
01258      This information currently not needed at this time, but may be
01259      needed later for better limit determination.
01260   */
01261   { double alpha, beta, gamma, Major, Minor;
01262     double Eccentricity, Ellipse_Area, Ellipse_angle;
01263     double max_horizontal_cross_section, max_vertical_cross_section;
01264     alpha = A+C;
01265     beta  = A-C;
01266     gamma = sqrt(beta*beta + B*B );
01267 
01268     if ( alpha - gamma <= MagickEpsilon )
01269       Major = MagickHuge;
01270     else
01271       Major = sqrt(2*F/(alpha - gamma));
01272     Minor = sqrt(2*F/(alpha + gamma));
01273 
01274     fprintf(stderr, "\tMajor=%lf; Minor=%lf\n",
01275          Major, Minor );
01276 
01277     /* other information about ellipse include... */
01278     Eccentricity = Major/Minor;
01279     Ellipse_Area = MagickPI*Major*Minor;
01280     Ellipse_angle =  atan2(B, A-C);
01281 
01282     fprintf(stderr, "\tAngle=%lf Area=%lf\n",
01283          RadiansToDegrees(Ellipse_angle), Ellipse_Area );
01284 
01285     /* Ellipse limits */
01286 
01287     /* orthogonal rectangle - improved ellipse */
01288     max_horizontal_orthogonal = sqrt(A); /* = sqrt(4*A*F/(4*A*C-B*B)) */
01289     max_vertical_orthogonal   = sqrt(C); /* = sqrt(4*C*F/(4*A*C-B*B)) */
01290 
01291     /* parallelogram bounds -- what we are using */
01292     max_horizontal_cross_section = sqrt(F/A);
01293     max_vertical_cross_section   = sqrt(F/C);
01294   }
01295 #endif
01296 
01297   /* Is default elliptical area, too small? Image being magnified?
01298      Switch to doing pure 'point' interpolation of the pixel.
01299      That is turn off  EWA Resampling.
01300   */
01301   if ( F <= F_UNITY ) {
01302     resample_filter->do_interpolate = MagickTrue;
01303     return;
01304   }
01305 
01306 
01307   /* If F is impossibly large, we may as well not bother doing any
01308    * form of resampling, as you risk an infinite resampled area.
01309   */
01310   if ( F > MagickHuge ) {
01311     resample_filter->limit_reached = MagickTrue;
01312     return;
01313   }
01314 
01315   /* Othogonal bounds of the ellipse */
01316   resample_filter->sqrtA = sqrt(A)+1.0;     /* Vertical Orthogonal Limit */
01317   resample_filter->sqrtC = sqrt(C)+1.0;     /* Horizontal Orthogonal Limit */
01318 
01319   /* Horizontally aligned Parallelogram fitted to ellipse */
01320   resample_filter->sqrtU = sqrt(F/A)+1.0;   /* Parallelogram Width */
01321   resample_filter->slope = -B/(2*A);        /* Slope of the parallelogram */
01322 
01323   /* The size of the area of the parallelogram we will be sampling */
01324   area = 4 * resample_filter->sqrtA * resample_filter->sqrtU;
01325 
01326   /* Absolute limit on the area to be resampled
01327    * This limit needs more work, as it gets too slow for
01328    * larger images involved with tiled views of the horizon. */
01329   if ( area > 20.0*resample_filter->image_area ) {
01330     resample_filter->limit_reached = MagickTrue;
01331     return;
01332   }
01333 
01334   /* Scale ellipse formula to directly fit the Filter Lookup Table */
01335   { register double scale;
01336     scale = (double)WLUT_WIDTH/F;
01337     resample_filter->A = A*scale;
01338     resample_filter->B = B*scale;
01339     resample_filter->C = C*scale;
01340     /* ..ple_filter->F = WLUT_WIDTH; -- hardcoded */
01341   }
01342 }
01343 
01344 /*
01345 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01346 %                                                                             %
01347 %                                                                             %
01348 %                                                                             %
01349 %   S e t R e s a m p l e F i l t e r                                         %
01350 %                                                                             %
01351 %                                                                             %
01352 %                                                                             %
01353 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01354 %
01355 %  SetResampleFilter() set the resampling filter lookup table based on a
01356 %  specific filter.  Note that the filter is used as a radial filter not as a
01357 %  two pass othogonally aligned resampling filter.
01358 %
01359 %  The default Filter, is Gaussian, which is the standard filter used by the
01360 %  original paper on the Elliptical Weighted Everage Algorithm. However other
01361 %  filters can also be used.
01362 %
01363 %  The format of the SetResampleFilter method is:
01364 %
01365 %    void SetResampleFilter(ResampleFilter *resample_filter,
01366 %      const FilterTypes filter,const double blur)
01367 %
01368 %  A description of each parameter follows:
01369 %
01370 %    o resample_filter: resampling resample_filterrmation structure
01371 %
01372 %    o filter: the resize filter for elliptical weighting LUT
01373 %
01374 %    o blur: filter blur factor (radial scaling) for elliptical weighting LUT
01375 %
01376 */
01377 MagickExport void SetResampleFilter(ResampleFilter *resample_filter,
01378   const FilterTypes filter,const double blur)
01379 {
01380   register int
01381      Q;
01382 
01383   double
01384      r_scale;
01385 
01386   ResizeFilter
01387      *resize_filter;
01388 
01389   assert(resample_filter != (ResampleFilter *) NULL);
01390   assert(resample_filter->signature == MagickSignature);
01391 
01392   resample_filter->filter = filter;
01393 
01394   /* Scale radius so it equals 1.0, at edge of ellipse when a
01395      default blurring factor of 1.0 is used.
01396 
01397      Note that these filters are being used as a radial filter, not as
01398      an othoginally alligned filter. How this effects results is still
01399      being worked out.
01400 
01401      Future: Direct use of teh resize filters in "resize.c" to set the lookup
01402      table, based on the filters working support window.
01403   */
01404   r_scale = sqrt(1.0/(double)WLUT_WIDTH)/blur;
01405   r_scale *= 2; /* for 2 pixel radius of Improved Elliptical Formula */
01406 
01407   switch ( filter ) {
01408   case PointFilter:
01409     /* This equivelent to turning off the EWA algroithm.
01410        Only Interpolated lookup will be used.  */
01411     break;
01412   default:
01413     /*
01414       Fill the LUT with a 1D resize filter function
01415       But make the Sinc/Bessel tapered window 2.0
01416       I also normalize the result so the filter is 1.0
01417     */
01418     resize_filter = AcquireResizeFilter(resample_filter->image,filter,
01419          (MagickRealType)1.0,MagickTrue,resample_filter->exception);
01420     if (resize_filter != (ResizeFilter *) NULL) {
01421       resample_filter->support = GetResizeFilterSupport(resize_filter);
01422       resample_filter->support /= blur; /* taken into account above */
01423       resample_filter->support *= resample_filter->support;
01424       resample_filter->support *= (double)WLUT_WIDTH/4;
01425       if ( resample_filter->support >= (double)WLUT_WIDTH )
01426            resample_filter->support = (double)WLUT_WIDTH;  /* hack */
01427       for(Q=0; Q<WLUT_WIDTH; Q++)
01428         if ( (double) Q < resample_filter->support )
01429           resample_filter->filter_lut[Q] = (double)
01430                GetResizeFilterWeight(resize_filter,sqrt((double)Q)*r_scale);
01431         else
01432           resample_filter->filter_lut[Q] = 0.0;
01433       resize_filter = DestroyResizeFilter(resize_filter);
01434       break;
01435     }
01436     else {
01437       (void) ThrowMagickException(resample_filter->exception,GetMagickModule(),
01438            ModuleError, "UnableToSetFilteringValue",
01439            "Fall back to default EWA gaussian filter");
01440     }
01441     /* FALLTHRU - on exception */
01442   /*case GaussianFilter:*/
01443   case UndefinedFilter:
01444     /*
01445       Create Normal Gaussian 2D Filter Weighted Lookup Table.
01446       A normal EWA guassual lookup would use   exp(Q*ALPHA)
01447       where  Q = distantce squared from 0.0 (center) to 1.0 (edge)
01448       and    ALPHA = -4.0*ln(2.0)  ==>  -2.77258872223978123767
01449       However the table is of length 1024, and equates to a radius of 2px
01450       thus needs to be scaled by  ALPHA*4/1024 and any blur factor squared
01451     */
01452     /*r_scale = -2.77258872223978123767*4/WLUT_WIDTH/blur/blur;*/
01453     r_scale = -2.77258872223978123767/WLUT_WIDTH/blur/blur;
01454     for(Q=0; Q<WLUT_WIDTH; Q++)
01455       resample_filter->filter_lut[Q] = exp((double)Q*r_scale);
01456     resample_filter->support = WLUT_WIDTH;
01457     break;
01458   }
01459   if (GetImageArtifact(resample_filter->image,"resample:verbose")
01460         != (const char *) NULL)
01461     /* Debug output of the filter weighting LUT
01462       Gnuplot the LUT with hoizontal adjusted to 'r' using...
01463         plot [0:2][-.2:1] "lut.dat" using (sqrt($0/1024)*2):1 with lines
01464       THe filter values is normalized for comparision
01465     */
01466     for(Q=0; Q<WLUT_WIDTH; Q++)
01467       printf("%lf\n", resample_filter->filter_lut[Q]
01468                         /resample_filter->filter_lut[0] );
01469   return;
01470 }
01471 
01472 /*
01473 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01474 %                                                                             %
01475 %                                                                             %
01476 %                                                                             %
01477 %   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       %
01478 %                                                                             %
01479 %                                                                             %
01480 %                                                                             %
01481 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01482 %
01483 %  SetResampleFilterInterpolateMethod() changes the interpolation method
01484 %  associated with the specified resample filter.
01485 %
01486 %  The format of the SetResampleFilterInterpolateMethod method is:
01487 %
01488 %      MagickBooleanType SetResampleFilterInterpolateMethod(
01489 %        ResampleFilter *resample_filter,const InterpolateMethod method)
01490 %
01491 %  A description of each parameter follows:
01492 %
01493 %    o resample_filter: the resample filter.
01494 %
01495 %    o method: the interpolation method.
01496 %
01497 */
01498 MagickExport MagickBooleanType SetResampleFilterInterpolateMethod(
01499   ResampleFilter *resample_filter,const InterpolatePixelMethod method)
01500 {
01501   assert(resample_filter != (ResampleFilter *) NULL);
01502   assert(resample_filter->signature == MagickSignature);
01503   assert(resample_filter->image != (Image *) NULL);
01504   if (resample_filter->debug != MagickFalse)
01505     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
01506       resample_filter->image->filename);
01507   resample_filter->interpolate=method;
01508   return(MagickTrue);
01509 }
01510 
01511 /*
01512 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01513 %                                                                             %
01514 %                                                                             %
01515 %                                                                             %
01516 %   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     %
01517 %                                                                             %
01518 %                                                                             %
01519 %                                                                             %
01520 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
01521 %
01522 %  SetResampleFilterVirtualPixelMethod() changes the virtual pixel method
01523 %  associated with the specified resample filter.
01524 %
01525 %  The format of the SetResampleFilterVirtualPixelMethod method is:
01526 %
01527 %      MagickBooleanType SetResampleFilterVirtualPixelMethod(
01528 %        ResampleFilter *resample_filter,const VirtualPixelMethod method)
01529 %
01530 %  A description of each parameter follows:
01531 %
01532 %    o resample_filter: the resample filter.
01533 %
01534 %    o method: the virtual pixel method.
01535 %
01536 */
01537 MagickExport MagickBooleanType SetResampleFilterVirtualPixelMethod(
01538   ResampleFilter *resample_filter,const VirtualPixelMethod method)
01539 {
01540   assert(resample_filter != (ResampleFilter *) NULL);
01541   assert(resample_filter->signature == MagickSignature);
01542   assert(resample_filter->image != (Image *) NULL);
01543   if (resample_filter->debug != MagickFalse)
01544     (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",
01545       resample_filter->image->filename);
01546   resample_filter->virtual_pixel=method;
01547   (void) SetCacheViewVirtualPixelMethod(resample_filter->view,method);
01548   return(MagickTrue);
01549 }

Generated on 19 Nov 2009 for MagickCore by  doxygen 1.6.1