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