From b6d08c51e1666814b156d42aef650f41d33a87ba Mon Sep 17 00:00:00 2001 From: anthony Date: Mon, 13 Sep 2010 01:17:04 +0000 Subject: [PATCH] Addition of LanczosChebyshev resize filter Fix of cylindrical Filter handling --- ChangeLog | 2 + magick/distort.c | 3 +- magick/layer.c | 125 ++++++++++++++-------------- magick/option.c | 1 + magick/resample.c | 6 ++ magick/resample.h | 1 + magick/resize.c | 203 +++++++++++++++++++++++++++++++++------------- 7 files changed, 222 insertions(+), 119 deletions(-) diff --git a/ChangeLog b/ChangeLog index 92c866acc..b09169de8 100644 --- a/ChangeLog +++ b/ChangeLog @@ -2,6 +2,8 @@ * The RGBO format is now listed as a supported format. 2010-09-07 6.6.4-1 Anthony Thyssen + * Added the Nicolas Robidoux and Chantal Racette Lanczos resize filter + function as "LanczosChebyshev" as faster alturnative to Lanczos. * Re-code Nicolas Robidoux and Chantal Racette Polynomial Approximation of the Sinc Trigonometric resize filter, as a proper filter to allow direct comparision and speed testing of the filter. diff --git a/magick/distort.c b/magick/distort.c index 9afc2deac..6fd2bad3b 100644 --- a/magick/distort.c +++ b/magick/distort.c @@ -1925,8 +1925,9 @@ MagickExport Image *DistortImage(const Image *image,DistortImageMethod method, exception); if (distort_image == (Image *) NULL) return((Image *) NULL); + /* if image is ColorMapped - change it to DirectClass */ if (SetImageStorageClass(distort_image,DirectClass) == MagickFalse) - { /* if image is ColorMapped - change it to DirectClass */ + { InheritException(exception,&distort_image->exception); distort_image=DestroyImage(distort_image); return((Image *) NULL); diff --git a/magick/layer.c b/magick/layer.c index e3ad302dd..db92fbd6b 100644 --- a/magick/layer.c +++ b/magick/layer.c @@ -388,7 +388,10 @@ MagickExport Image *DisposeImages(const Image *image,ExceptionInfo *exception) *dispose_images; register Image - *next; + *curr; + + RectangleInfo + bounds; /* Run the image through the animation sequence @@ -399,19 +402,19 @@ MagickExport Image *DisposeImages(const Image *image,ExceptionInfo *exception) (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename); assert(exception != (ExceptionInfo *) NULL); assert(exception->signature == MagickSignature); - next=GetFirstImageInList(image); - dispose_image=CloneImage(next,next->page.width,next->page.height,MagickTrue, + curr=GetFirstImageInList(image); + dispose_image=CloneImage(curr,curr->page.width,curr->page.height,MagickTrue, exception); if (dispose_image == (Image *) NULL) return((Image *) NULL); - dispose_image->page=next->page; + dispose_image->page=curr->page; dispose_image->page.x=0; dispose_image->page.y=0; dispose_image->dispose=NoneDispose; dispose_image->background_color.opacity=(Quantum) TransparentOpacity; (void) SetImageBackgroundColor(dispose_image); dispose_images=NewImageList(); - for ( ; next != (Image *) NULL; next=GetNextImageInList(next)) + for ( ; curr != (Image *) NULL; curr=GetNextImageInList(curr)) { Image *current_image; @@ -426,19 +429,17 @@ MagickExport Image *DisposeImages(const Image *image,ExceptionInfo *exception) dispose_image=DestroyImage(dispose_image); return((Image *) NULL); } - (void) CompositeImage(current_image,next->matte != MagickFalse ? - OverCompositeOp : CopyCompositeOp,next,next->page.x,next->page.y); + (void) CompositeImage(current_image,curr->matte != MagickFalse ? + OverCompositeOp : CopyCompositeOp,curr,curr->page.x,curr->page.y); + /* Handle Background dispose: image is displayed for the delay period. */ - if (next->dispose == BackgroundDispose) + if (curr->dispose == BackgroundDispose) { - RectangleInfo - bounds; - - bounds=next->page; - bounds.width=next->columns; - bounds.height=next->rows; + bounds=curr->page; + bounds.width=curr->columns; + bounds.height=curr->rows; if (bounds.x < 0) { bounds.width+=bounds.x; @@ -458,20 +459,21 @@ MagickExport Image *DisposeImages(const Image *image,ExceptionInfo *exception) /* Select the appropriate previous/disposed image. */ - if (next->dispose == PreviousDispose) + if (curr->dispose == PreviousDispose) current_image=DestroyImage(current_image); else { dispose_image=DestroyImage(dispose_image); dispose_image=current_image; + current_image=(Image *)NULL; } + /* + Save the dispose image just calculated for return. + */ { Image *dispose; - /* - Save the dispose image just calculated for return. - */ dispose=CloneImage(dispose_image,0,0,MagickTrue,exception); if (dispose == (Image *) NULL) { @@ -479,12 +481,12 @@ MagickExport Image *DisposeImages(const Image *image,ExceptionInfo *exception) dispose_image=DestroyImage(dispose_image); return((Image *) NULL); } - (void) CloneImageProfiles(dispose,next); - (void) CloneImageProperties(dispose,next); - (void) CloneImageArtifacts(dispose,next); + (void) CloneImageProfiles(dispose,curr); + (void) CloneImageProperties(dispose,curr); + (void) CloneImageArtifacts(dispose,curr); dispose->page.x=0; dispose->page.y=0; - dispose->dispose=next->dispose; + dispose->dispose=curr->dispose; AppendImageToList(&dispose_images,dispose); } } @@ -989,7 +991,7 @@ static Image *OptimizeLayerFrames(const Image *image, *disposals; register const Image - *next; + *curr; register ssize_t i; @@ -1011,10 +1013,10 @@ static Image *OptimizeLayerFrames(const Image *image, /* Ensure all the images are the same size */ - next=GetFirstImageInList(image); - for (; next != (Image *) NULL; next=GetNextImageInList(next)) + curr=GetFirstImageInList(image); + for (; curr != (Image *) NULL; curr=GetNextImageInList(curr)) { - if ((next->columns != image->columns) || (next->rows != image->rows)) + if ((curr->columns != image->columns) || (curr->rows != image->rows)) ThrowImageException(OptionError,"ImagesAreNotTheSameSize"); /* FUTURE: also check that image is also fully coalesced (full page) @@ -1024,9 +1026,9 @@ static Image *OptimizeLayerFrames(const Image *image, /* Allocate memory (times 2 if we allow the use of frame duplications) */ - next=GetFirstImageInList(image); + curr=GetFirstImageInList(image); bounds=(RectangleInfo *) AcquireQuantumMemory((size_t) - GetImageListLength(next),(add_frames != MagickFalse ? 2UL : 1UL)* + GetImageListLength(curr),(add_frames != MagickFalse ? 2UL : 1UL)* sizeof(*bounds)); if (bounds == (RectangleInfo *) NULL) ThrowImageException(ResourceLimitError,"MemoryAllocationFailed"); @@ -1041,7 +1043,7 @@ static Image *OptimizeLayerFrames(const Image *image, /* Initialise Previous Image as fully transparent */ - prev_image=CloneImage(next,next->page.width,next->page.height, + prev_image=CloneImage(curr,curr->page.width,curr->page.height, MagickTrue,exception); if (prev_image == (Image *) NULL) { @@ -1049,7 +1051,7 @@ static Image *OptimizeLayerFrames(const Image *image, disposals=(DisposeType *) RelinquishMagickMemory(disposals); return((Image *) NULL); } - prev_image->page=next->page; /* ERROR: <-- should not be need, but is! */ + prev_image->page=curr->page; /* ERROR: <-- should not be need, but is! */ prev_image->page.x=0; prev_image->page.y=0; prev_image->dispose=NoneDispose; @@ -1065,7 +1067,7 @@ static Image *OptimizeLayerFrames(const Image *image, fprintf(stderr, "frame %.20g :-\n", (double) i); #endif disposals[0]=NoneDispose; - bounds[0]=CompareImageBounds(prev_image,next,CompareAnyLayer,exception); + bounds[0]=CompareImageBounds(prev_image,curr,CompareAnyLayer,exception); #if DEBUG_OPT_FRAME fprintf(stderr, "overlay: %.20gx%.20g%+.20g%+.20g\n\n", (double) bounds[i].width,(double) bounds[i].height, @@ -1081,8 +1083,8 @@ static Image *OptimizeLayerFrames(const Image *image, dup_bounds.height=0; dup_bounds.x=0; dup_bounds.y=0; - next=GetNextImageInList(next); - for ( ; next != (const Image *) NULL; next=GetNextImageInList(next)) + curr=GetNextImageInList(curr); + for ( ; curr != (const Image *) NULL; curr=GetNextImageInList(curr)) { #if DEBUG_OPT_FRAME fprintf(stderr, "frame %.20g :-\n", (double) i); @@ -1090,8 +1092,8 @@ static Image *OptimizeLayerFrames(const Image *image, /* Assume none disposal is the best */ - bounds[i]=CompareImageBounds(next->previous,next,CompareAnyLayer,exception); - cleared=IsBoundsCleared(next->previous,next,&bounds[i],exception); + bounds[i]=CompareImageBounds(curr->previous,curr,CompareAnyLayer,exception); + cleared=IsBoundsCleared(curr->previous,curr,&bounds[i],exception); disposals[i-1]=NoneDispose; #if DEBUG_OPT_FRAME fprintf(stderr, "overlay: %.20gx%.20g%+.20g%+.20g%s%s\n", @@ -1120,8 +1122,8 @@ static Image *OptimizeLayerFrames(const Image *image, /* Compare a none disposal against a previous disposal */ - try_bounds=CompareImageBounds(prev_image,next,CompareAnyLayer,exception); - try_cleared=IsBoundsCleared(prev_image,next,&try_bounds,exception); + try_bounds=CompareImageBounds(prev_image,curr,CompareAnyLayer,exception); + try_cleared=IsBoundsCleared(prev_image,curr,&try_bounds,exception); #if DEBUG_OPT_FRAME fprintf(stderr, "test_prev: %.20gx%.20g%+.20g%+.20g%s\n", (double) try_bounds.width,(double) try_bounds.height, @@ -1150,8 +1152,8 @@ static Image *OptimizeLayerFrames(const Image *image, dup_bounds.width=dup_bounds.height=0; /* no dup, no pixel added */ if ( add_frames ) { - dup_image=CloneImage(next->previous,next->previous->page.width, - next->previous->page.height,MagickTrue,exception); + dup_image=CloneImage(curr->previous,curr->previous->page.width, + curr->previous->page.height,MagickTrue,exception); if (dup_image == (Image *) NULL) { bounds=(RectangleInfo *) RelinquishMagickMemory(bounds); @@ -1159,9 +1161,9 @@ static Image *OptimizeLayerFrames(const Image *image, prev_image=DestroyImage(prev_image); return((Image *) NULL); } - dup_bounds=CompareImageBounds(dup_image,next,CompareClearLayer,exception); + dup_bounds=CompareImageBounds(dup_image,curr,CompareClearLayer,exception); ClearBounds(dup_image,&dup_bounds); - try_bounds=CompareImageBounds(dup_image,next,CompareAnyLayer,exception); + try_bounds=CompareImageBounds(dup_image,curr,CompareAnyLayer,exception); if ( cleared || dup_bounds.width*dup_bounds.height +try_bounds.width*try_bounds.height @@ -1178,8 +1180,8 @@ static Image *OptimizeLayerFrames(const Image *image, /* Now compare against a simple background disposal */ - bgnd_image=CloneImage(next->previous,next->previous->page.width, - next->previous->page.height,MagickTrue,exception); + bgnd_image=CloneImage(curr->previous,curr->previous->page.width, + curr->previous->page.height,MagickTrue,exception); if (bgnd_image == (Image *) NULL) { bounds=(RectangleInfo *) RelinquishMagickMemory(bounds); @@ -1191,8 +1193,8 @@ static Image *OptimizeLayerFrames(const Image *image, } bgnd_bounds=bounds[i-1]; /* interum bounds of the previous image */ ClearBounds(bgnd_image,&bgnd_bounds); - try_bounds=CompareImageBounds(bgnd_image,next,CompareAnyLayer,exception); - try_cleared=IsBoundsCleared(bgnd_image,next,&try_bounds,exception); + try_bounds=CompareImageBounds(bgnd_image,curr,CompareAnyLayer,exception); + try_cleared=IsBoundsCleared(bgnd_image,curr,&try_bounds,exception); #if DEBUG_OPT_FRAME fprintf(stderr, "background: %s\n", try_cleared?"(pixels cleared)":""); @@ -1205,8 +1207,7 @@ static Image *OptimizeLayerFrames(const Image *image, include the pixels that are cleared. This guaranteed to work, though may not be the most optimized solution. */ - // ERROR HERE???? next->previous??? - try_bounds=CompareImageBounds(next->previous,next,CompareClearLayer,exception); + try_bounds=CompareImageBounds(curr->previous,curr,CompareClearLayer,exception); #if DEBUG_OPT_FRAME fprintf(stderr, "expand_clear: %.20gx%.20g%+.20g%+.20g%s\n", (double) try_bounds.width,(double) try_bounds.height, @@ -1263,20 +1264,20 @@ static Image *OptimizeLayerFrames(const Image *image, * to see, or writet he image at this point it is hard to tell what is wrong! * Only CompareOverlay seemed to return something sensible. */ - try_bounds=CompareImageBounds(bgnd_image,next,CompareClearLayer,exception); + try_bounds=CompareImageBounds(bgnd_image,curr,CompareClearLayer,exception); fprintf(stderr, "expand_ctst: %.20gx%.20g%+.20g%+.20g\n", (double) try_bounds.width,(double) try_bounds.height, (double) try_bounds.x,(double) try_bounds.y ); - try_bounds=CompareImageBounds(bgnd_image,next,CompareAnyLayer,exception); - try_cleared=IsBoundsCleared(bgnd_image,next,&try_bounds,exception); + try_bounds=CompareImageBounds(bgnd_image,curr,CompareAnyLayer,exception); + try_cleared=IsBoundsCleared(bgnd_image,curr,&try_bounds,exception); fprintf(stderr, "expand_any : %.20gx%.20g%+.20g%+.20g%s\n", (double) try_bounds.width,(double) try_bounds.height, (double) try_bounds.x,(double) try_bounds.y, try_cleared?" (pixels cleared)":""); #endif - try_bounds=CompareImageBounds(bgnd_image,next,CompareOverlayLayer,exception); + try_bounds=CompareImageBounds(bgnd_image,curr,CompareOverlayLayer,exception); #if DEBUG_OPT_FRAME - try_cleared=IsBoundsCleared(bgnd_image,next,&try_bounds,exception); + try_cleared=IsBoundsCleared(bgnd_image,curr,&try_bounds,exception); fprintf(stderr, "expand_test: %.20gx%.20g%+.20g%+.20g%s\n", (double) try_bounds.width,(double) try_bounds.height, (double) try_bounds.x,(double) try_bounds.y, @@ -1335,8 +1336,8 @@ static Image *OptimizeLayerFrames(const Image *image, bgnd_image=DestroyImage(bgnd_image); if ( disposals[i-1] == NoneDispose ) { - prev_image=CloneImage(next->previous,next->previous->page.width, - next->previous->page.height,MagickTrue,exception); + prev_image=CloneImage(curr->previous,curr->previous->page.width, + curr->previous->page.height,MagickTrue,exception); if (prev_image == (Image *) NULL) { bounds=(RectangleInfo *) RelinquishMagickMemory(bounds); @@ -1371,21 +1372,21 @@ static Image *OptimizeLayerFrames(const Image *image, */ sans_exception=AcquireExceptionInfo(); i=0; - next=GetFirstImageInList(image); + curr=GetFirstImageInList(image); optimized_image=NewImageList(); - while ( next != (const Image *) NULL ) + while ( curr != (const Image *) NULL ) { - prev_image=CloneImage(next,0,0,MagickTrue,exception); + prev_image=CloneImage(curr,0,0,MagickTrue,exception); if (prev_image == (Image *) NULL) break; if ( disposals[i] == DelDispose ) { size_t time = 0; while ( disposals[i] == DelDispose ) { - time += next->delay*1000/next->ticks_per_second; - next=GetNextImageInList(next); + time += curr->delay*1000/curr->ticks_per_second; + curr=GetNextImageInList(curr); i++; } - time += next->delay*1000/next->ticks_per_second; + time += curr->delay*1000/curr->ticks_per_second; prev_image->ticks_per_second = 100L; prev_image->delay = time*prev_image->ticks_per_second/1000; } @@ -1399,14 +1400,14 @@ static Image *OptimizeLayerFrames(const Image *image, bgnd_image->dispose=NoneDispose; } else - next=GetNextImageInList(next); + curr=GetNextImageInList(curr); AppendImageToList(&optimized_image,bgnd_image); i++; } sans_exception=DestroyExceptionInfo(sans_exception); bounds=(RectangleInfo *) RelinquishMagickMemory(bounds); disposals=(DisposeType *) RelinquishMagickMemory(disposals); - if (next != (Image *) NULL) + if (curr != (Image *) NULL) { optimized_image=DestroyImageList(optimized_image); return((Image *) NULL); diff --git a/magick/option.c b/magick/option.c index d484350f1..5cebcc04c 100644 --- a/magick/option.c +++ b/magick/option.c @@ -947,6 +947,7 @@ static const OptionInfo { "Kaiser", (ssize_t) KaiserFilter, MagickFalse }, { "Lagrange", (ssize_t) LagrangeFilter, MagickFalse }, { "Lanczos", (ssize_t) LanczosFilter, MagickFalse }, + { "LanczosChebyshev", (ssize_t) LanczosChebyshevFilter, MagickFalse }, { "Mitchell", (ssize_t) MitchellFilter, MagickFalse }, { "Parzen", (ssize_t) ParzenFilter, MagickFalse }, { "Point", (ssize_t) PointFilter, MagickFalse }, diff --git a/magick/resample.c b/magick/resample.c index 296964501..b4c8611e8 100644 --- a/magick/resample.c +++ b/magick/resample.c @@ -1464,6 +1464,12 @@ MagickExport void SetResampleFilter(ResampleFilter *resample_filter, plot [0:2][-.2:1] "lut.dat" using (sqrt($0/1024)*2):1 with lines The filter values is normalized for comparision */ + printf("# Resampling Filter LUT (%d values)\n", WLUT_WIDTH); + printf("#\n"); + printf("# Plot using the gnuplot command\n"); + printf("# plot [0:2][-.2:1] \"lut.dat\" using "); + printf("sqrt($0/%d)*2):1 with lines\n", WLUT_WIDTH); + printf("#\n"); for(Q=0; Qfilter_lut[Q] /resample_filter->filter_lut[0] ); diff --git a/magick/resample.h b/magick/resample.h index cd8b6e62e..06794c0b6 100644 --- a/magick/resample.h +++ b/magick/resample.h @@ -55,6 +55,7 @@ typedef enum BohmanFilter, BartlettFilter, SincPolynomialFilter, + LanczosChebyshevFilter, SentinelFilter /* a count of all the filters, not a real filter */ } FilterTypes; diff --git a/magick/resize.c b/magick/resize.c index 0a7a0451e..a3dc23d9f 100644 --- a/magick/resize.c +++ b/magick/resize.c @@ -96,7 +96,8 @@ struct _ResizeFilter */ static MagickRealType I0(MagickRealType x), - BesselOrderOne(MagickRealType); + BesselOrderOne(MagickRealType), + Sinc(const MagickRealType,const ResizeFilter *); /* %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% @@ -144,8 +145,8 @@ static MagickRealType Bessel(const MagickRealType x, http://www.ph.ed.ac.uk/%7ewjh/teaching/mo/slides/lens/lens.pdf. */ if (x == 0.0) - return((MagickRealType) (0.25*MagickPIL)); - return(BesselOrderOne((MagickRealType) (MagickPIL*x))/(x+x)); + return(0.25*MagickPIL); + return(BesselOrderOne(MagickPIL*x)/(x+x)); } static MagickRealType Blackman(const MagickRealType x, @@ -157,7 +158,7 @@ static MagickRealType Blackman(const MagickRealType x, Refactored by Chantal Racette and Nicolas Robidoux to one trig call and five flops. */ - const double pix = (double) (MagickPIL*x); + const double pix = MagickPIL*x; const MagickRealType cospix = cos(pix); return(0.34+cospix*(0.5+cospix*0.16)); } @@ -169,8 +170,8 @@ static MagickRealType Bohman(const MagickRealType x, Bohman: 2rd Order cosine windowing function: (1-x) cos(pi x) + sin(pi x) / pi. */ - const double pix = (double) (MagickPIL*x); - return((MagickRealType) ((1-x)*cos(pix)+(1.0/MagickPIL)*sin(pix))); + const double pix = MagickPIL*x; + return((1-x)*cos(pix)+(1.0/MagickPIL)*sin(pix)); } static MagickRealType Box(const MagickRealType magick_unused(x), @@ -227,18 +228,11 @@ static MagickRealType Gaussian(const MagickRealType x, const ResizeFilter *magick_unused(resize_filter)) { /* - Normalized Gaussian with variance 1/2 (by default): - 1/sqrt(2 pi sigma^2) exp(-x^2/(2 sigma^2)) + 1D Gaussian with sigm=1/2 + exp(-2 x^2)/sqrt(pi/2)) */ - #define MagickGAUSSIANSIGMAL 0.5L - /* - Change the value of MagickGAUSSIANSIGMAL if you want to override - the default. - */ - const MagickRealType sigma2 = MagickGAUSSIANSIGMAL*MagickGAUSSIANSIGMAL; - const MagickRealType alpha = -1.0/(2.0*sigma2); - const MagickRealType normalizer = sqrt((double) (1.0/(2.0*MagickPIL*sigma2))); - return(normalizer*exp((double) (alpha*x*x))); + const MagickRealType alpha = 2.0/MagickSQ2PI; + return(exp(-(double)(2.0*x*x))*alpha); } static MagickRealType Hanning(const MagickRealType x, @@ -247,7 +241,7 @@ static MagickRealType Hanning(const MagickRealType x, /* Cosine window function: .5 + .5 cos(pi x). */ - const double pix = (double) (MagickPIL*x); + const double pix = MagickPIL*x; const MagickRealType cospix = cos(pix); return(0.5+0.5*cospix); } @@ -258,7 +252,7 @@ static MagickRealType Hamming(const MagickRealType x, /* Offset cosine window function: .54 + .46 cos(pi x). */ - const double pix = (double) (MagickPIL*x); + const double pix = MagickPIL*x; const MagickRealType cospix = cos(pix); return(0.54+0.46*cospix); } @@ -312,6 +306,93 @@ static MagickRealType Lagrange(const MagickRealType x, return(value); } +static MagickRealType LanczosChebyshev(const MagickRealType x, + const ResizeFilter *resize_filter) +{ + /* + Lanczos is normally a Sinc() windowed Sinc filter, but that requires + 2 calls to a trigonometric function, whcih can slow it down. + + This version is computed with only one call to a trigonometric function + with a recursive formula, based on the Chebyshev method for the + computation of sines and cosines of multiples of an angle, discovered by + Nicolas Robidoux (pending the discovery of an earlier discoverer) with the + assistance of Chantal Racette. + + Chebyshev method of sines of multiple angles + http://en.wikipedia.org/wiki/List_of_trigonometric_identities#Chebyshev_method + + The filter only works with an interger support to a limit of 10.0, at + which point it will fallback to a Sinc windowed Sinc (after all the work + was done anyway). + + FUTURE: AcquireResizeFilter() should override this function if support + greater than 10, or a non-integer support is requested. + */ + const double support = resize_filter->support; + const MagickRealType x2 = x*x; + if (x2 == 0) + return(1.0); + + { + const MagickRealType pi2 = MagickPIL*MagickPIL; + const MagickRealType c = cos((double) ((MagickPIL/support)*x)); + const MagickRealType s2 = 1 - c * c; + const MagickRealType s2c = s2 * c; + const MagickRealType ss2 = s2c + s2c; + /* + Really, we want to use Lanczos 2 if the support is <= 2, + but doing things as follows allows for inaccuracy in the + support computation. + */ + if (support<2.5) + return((2/pi2)/x2*ss2); + { + const MagickRealType ss3 = s2*(3+-4*s2); + if (support<3.5) + return((3/pi2)/x2*ss3); + { + const MagickRealType t4 = c*ss3; + const MagickRealType ss4 = t4-ss2+t4; + if (support<4.5) + return((4/pi2)/x2*ss4); + { + const MagickRealType t5 = c*ss4; + const MagickRealType ss5 = t5-ss3+t5; + if (support<5.5) + return((5/pi2)/x2*ss5); + { + const MagickRealType t6 = c*ss5; + const MagickRealType ss6 = t6-ss4+t6; + if (support<6.5) + return((6/pi2)/x2*ss6); + { + const MagickRealType t7 = c*ss6; + const MagickRealType ss7 = t7-ss5+t7; + if (support<7.5) + return((7/pi2)/x2*ss7); + { + const MagickRealType t8 = c*ss7; + const MagickRealType ss8 = t8-ss6+t8; + if (support<8.5) + return((8/pi2)/x2*ss8); + { + const MagickRealType t9 = c*ss8; + const MagickRealType ss9 = t9-ss7+t9; + if (support<9.5) + return((9/pi2)/x2*ss9); + if (support<10.5) + { + const MagickRealType t10 = c*ss9; + const MagickRealType ss10 = t10-ss8+t10; + return((10/pi2)/x2*ss10); + } + return( + Sinc((double) x,resize_filter) * + Sinc((double) (x/support),resize_filter)); + } } } } } } } } +} + static MagickRealType Quadratic(const MagickRealType x, const ResizeFilter *magick_unused(resize_filter)) { @@ -334,7 +415,7 @@ static MagickRealType Sinc(const MagickRealType x, if (x == 0.0) return(1.0); { - const MagickRealType pix = (MagickRealType) (MagickPIL*x); + const MagickRealType pix = MagickPIL*x; const MagickRealType sinpix = sin((double) pix); return(sinpix/pix); } @@ -352,7 +433,7 @@ static MagickRealType SincPolynomial(const MagickRealType x, const MagickRealType xx = x*x; if (xx > 16.0) { - const MagickRealType pix = (MagickRealType) (MagickPIL*x); + const MagickRealType pix = MagickPIL*x; const MagickRealType sinpix = sin((double) pix); return(sinpix/pix); } @@ -421,7 +502,7 @@ static MagickRealType SincPolynomial(const MagickRealType x, const MagickRealType c14 = 0.374841980075726557899013574367932640586e-25L; const MagickRealType c15 = -0.138632329047117683500928913798808544919e-27L; const MagickRealType p = c0+xx*(c1+xx*(c2+xx*(c3+xx*(c4+xx*(c5+xx*(c6+xx* - (c7+xx*(c8+xx*(c9+xx*(c10+xx*(c11+xx*(c12+xx*(c13+xx*(c14+xx*c15 + (c7+xx*(c8+xx*(c9+xx*(c10+xx*(c11+xx*(c12+xx*(c13+xx*(c14+xx*15 )))))))))))))); #endif return((xx-1.0)*(xx-4.0)*(xx-9.0)*(xx-16.0)*p); @@ -477,6 +558,9 @@ static MagickRealType Welsh(const MagickRealType x, % Blackman Hanning Hamming % Kaiser Lanczos (Sinc) % +% Polynomial Approximations (high precision fast versions) +% SincPolynomial +% % FIR filters are used as is, and are limited by that filters support window % (unless over-ridden). 'Gaussian' while classed as an IIR filter, is also % simply clipped by its support size (1.5). @@ -492,8 +576,8 @@ static MagickRealType Welsh(const MagickRealType x, % recommended as it removes the correct filter selection for different % filtering image operations. Selecting a window filtering method is better. % -% Lanczos is purely special case of a Sinc windowed Sinc, but defaulting to -% a 3 lobe support, rather that the default 4 lobe support. +% Lanczos is a special case of a Sinc windowed Sinc, but defaulting to +% a 3 lobe support, rather that the default 4 lobe support of the others. % % Special options can be used to override specific, or all the filter % settings. However doing so is not advisible unless you have expert @@ -523,25 +607,27 @@ static MagickRealType Welsh(const MagickRealType x, % used for simple filters like FIR filters, and the Gaussian Filter. % This will override any 'filter:lobes' option. % +% "filter:win-support" Scale windowing function to this size instead. +% This causes the windowing (or self-windowing Lagrange filter) to act +% is if the support window it much much larger than what is actually +% supplied to the calling operator. The filter however is still +% clipped to the real support size given, by the support range suppiled +% to the caller. If unset this will equal the normal filter support +% size. +% % "filter:blur" Scale the filter and support window by this amount. % A value >1 will generally result in a more burred image with % more ringing effects, while a value <1 will sharpen the % resulting image with more aliasing and Morie effects. % -% "filter:win-support" Scale windowing function to this size instead. -% This causes the windowing (or self-windowing Lagrange filter) -% to act is if the support winodw it much much larger than what -% is actually supplied to the calling operator. The filter however -% is still clipped to the real support size given. If unset this -% will equal the normal filter support size. -% % "filter:b" % "filter:c" Override the preset B,C values for a Cubic type of filter % If only one of these are given it is assumes to be a 'Keys' % type of filter such that B+2C=1, where Keys 'alpha' value = C % -% "filter:verbose" Output verbose plotting data for graphing the -% resulting filter over the whole support range (with blur effect). +% "filter:verbose" Output the exact results of the filter selections +% made, as well as plotting data for graphing the resulting filter +% over support range (blur adjusted). % % Set a true un-windowed Sinc filter with 10 lobes (very slow) % -set option:filter:filter Sinc @@ -562,9 +648,10 @@ static MagickRealType Welsh(const MagickRealType x, % o image: the image. % % o filter: the filter type, defining a preset filter, window and support. +% The artifact settings listed above will override those selections. % % o blur: blur the filter by this amount, use 1.0 if unknown. Image -% artifact "filter:blur" will override this old usage +% artifact "filter:blur" will override this internal usage. % % o radial: 1D orthogonal filter (Sinc) or 2D radial filter (Bessel) % @@ -634,7 +721,8 @@ MagickExport ResizeFilter *AcquireResizeFilter(const Image *image, { LagrangeFilter, BoxFilter }, /* Lagrange self-windowing filter */ { SincFilter, BohmanFilter }, /* Bohman -- 2*Cosine-Sinc */ { SincFilter, TriangleFilter }, /* Bartlett -- Triangle-Sinc */ - { SincPolynomialFilter, BlackmanFilter } /* Polynomial Approximated Sinc */ + { SincPolynomialFilter, BlackmanFilter },/* Polynomial Approximated Sinc */ + { LanczosChebyshevFilter, BoxFilter } /* Lanczos using Chebyshev Opt */ }; /* Table maping the filter/window function from the above table to the actual @@ -676,7 +764,8 @@ MagickExport ResizeFilter *AcquireResizeFilter(const Image *image, { Lagrange, 2.0f, 1.0f, 0.0f, 0.0f }, /* Lagrangian Filter */ { Bohman, 1.0f, 1.0f, 0.0f, 0.0f }, /* Bohman, 2*Cosine windowing */ { Triangle, 1.0f, 1.0f, 0.0f, 0.0f }, /* Bartlett, Triangle windowing */ - { SincPolynomial, 4.0f, 1.0f, 0.0f, 0.0f } /* Poly Approx Sinc */ + { SincPolynomial, 4.f, 1.f, 0.f, 0.f }, /* Polynomial Approximate Sinc */ + { LanczosChebyshev, 3.f, 1.f, 0.f, 0.f } /* Lanczos using Chevbyshev Opt */ }; /* The known zero crossings of the Bessel() or the Jinc(x*PI) function found @@ -732,15 +821,17 @@ MagickExport ResizeFilter *AcquireResizeFilter(const Image *image, resize_filter->blur=StringToDouble(artifact); if (resize_filter->blur < MagickEpsilon) resize_filter->blur=(MagickRealType) MagickEpsilon; - if ((cylindrical != MagickFalse) && (filter != SincFilter)) + if (cylindrical != MagickFalse) switch (filter_type) { case SincFilter: { /* Promote 1D Sinc Filter to a 2D Bessel filter. + As long as the user did not directly request a 'Sinc' filter */ - filter_type=BesselFilter; + if ( filter != SincFilter ) + filter_type=BesselFilter; break; } case LanczosFilter: @@ -759,18 +850,18 @@ MagickExport ResizeFilter *AcquireResizeFilter(const Image *image, Paul Heckbert's paper on EWA resampling. FUTURE: to be reviewed. */ - resize_filter->blur*=2.0*log(2.0)/sqrt((double) (2.0/MagickPI)); - break; - } - case BesselFilter: - { - /* - Filters with a 1.0 zero root crossing by the first bessel zero. - */ - resize_filter->blur*=bessel_zeros[0]; + /*resize_filter->blur*=2.0*log(2.0)/sqrt(2.0/MagickPI);*/ + /* It seems the formula above is wrong as both 1D and 2D gaussian + * is the same function just with different normalization weights + */ break; } default: + /* What about other filters to make them 'cylindrical frendly? + * For example Mitchel is actually quite close to a cyldrical + * Lanczos (Bessel-Bessel) support 2, though not quite there. + * Are there other well known 'cylindrical' specific filters? + */ break; } artifact=GetImageArtifact(image,"filter:filter"); @@ -864,7 +955,7 @@ MagickExport ResizeFilter *AcquireResizeFilter(const Image *image, if (artifact != (const char *) NULL) resize_filter->window_support=fabs(StringToDouble(artifact)); /* - Adjust window function X scaling to fit + Adjust window function X scaling to fit Avoids a division on every filter call. */ resize_filter->scale /= resize_filter->window_support; @@ -937,22 +1028,22 @@ MagickExport ResizeFilter *AcquireResizeFilter(const Image *image, /* Report Filter Details */ - support=GetResizeFilterSupport(resize_filter); /* support range */ + support = GetResizeFilterSupport(resize_filter); /* support range */ (void) fprintf(stdout,"#\n# Resize Filter (for graphing)\n#\n"); (void) fprintf(stdout,"# filter = %s\n", - MagickOptionToMnemonic(MagickFilterOptions,filter_type)); + MagickOptionToMnemonic(MagickFilterOptions, filter_type) ); (void) fprintf(stdout,"# window = %s\n", - MagickOptionToMnemonic(MagickFilterOptions,window_type)); + MagickOptionToMnemonic(MagickFilterOptions, window_type) ); (void) fprintf(stdout,"# support = %.*g\n", - GetMagickPrecision(),(double) resize_filter->support); + GetMagickPrecision(),resize_filter->support ); (void) fprintf(stdout,"# win-support = %.*g\n", - GetMagickPrecision(),(double) resize_filter->window_support); + GetMagickPrecision(),resize_filter->window_support ); (void) fprintf(stdout,"# blur = %.*g\n", - GetMagickPrecision(),(double) resize_filter->blur); + GetMagickPrecision(),resize_filter->blur ); (void) fprintf(stdout,"# blurred_support = %.*g\n", - GetMagickPrecision(),(double) support); + GetMagickPrecision(),support); (void) fprintf(stdout,"# B,C = %.*g,%.*g\n", - GetMagickPrecision(),B,GetMagickPrecision(),(double) C); + GetMagickPrecision(),B, GetMagickPrecision(),C); (void) fprintf(stdout,"#\n"); /* Output values of resulting filter graph -- for graphing filter result. -- 2.50.1