1 /* 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 3 % % 4 % % 5 % % 6 % FFFFF OOO U U RRRR IIIII EEEEE RRRR % 7 % F O O U U R R I E R R % 8 % FFF O O U U RRRR I EEE RRRR % 9 % F O O U U R R I E R R % 10 % F OOO UUU R R IIIII EEEEE R R % 11 % % 12 % % 13 % MagickCore Discrete Fourier Transform Methods % 14 % % 15 % Software Design % 16 % Sean Burke % 17 % Fred Weinhaus % 18 % Cristy % 19 % July 2009 % 20 % % 21 % % 22 % Copyright 1999-2016 ImageMagick Studio LLC, a non-profit organization % 23 % dedicated to making software imaging solutions freely available. % 24 % % 25 % You may not use this file except in compliance with the License. You may % 26 % obtain a copy of the License at % 27 % % 28 % http://www.imagemagick.org/script/license.php % 29 % % 30 % Unless required by applicable law or agreed to in writing, software % 31 % distributed under the License is distributed on an "AS IS" BASIS, % 32 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. % 33 % See the License for the specific language governing permissions and % 34 % limitations under the License. % 35 % % 36 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 37 % 38 % 39 % 40 */ 41 42 /* 44 Include declarations. 45 */ 46 #include "MagickCore/studio.h" 47 #include "MagickCore/artifact.h" 48 #include "MagickCore/attribute.h" 49 #include "MagickCore/blob.h" 50 #include "MagickCore/cache.h" 51 #include "MagickCore/image.h" 52 #include "MagickCore/image-private.h" 53 #include "MagickCore/list.h" 54 #include "MagickCore/fourier.h" 55 #include "MagickCore/log.h" 56 #include "MagickCore/memory_.h" 57 #include "MagickCore/monitor.h" 58 #include "MagickCore/monitor-private.h" 59 #include "MagickCore/pixel-accessor.h" 60 #include "MagickCore/pixel-private.h" 61 #include "MagickCore/property.h" 62 #include "MagickCore/quantum-private.h" 63 #include "MagickCore/resource_.h" 64 #include "MagickCore/string-private.h" 65 #include "MagickCore/thread-private.h" 66 #if defined(MAGICKCORE_FFTW_DELEGATE) 67 #if defined(MAGICKCORE_HAVE_COMPLEX_H) 68 #include <complex.h> 69 #endif 70 #include <fftw3.h> 71 #if !defined(MAGICKCORE_HAVE_CABS) 72 #define cabs(z) (sqrt(z[0]*z[0]+z[1]*z[1])) 73 #endif 74 #if !defined(MAGICKCORE_HAVE_CARG) 75 #define carg(z) (atan2(cimag(z),creal(z))) 76 #endif 77 #if !defined(MAGICKCORE_HAVE_CIMAG) 78 #define cimag(z) (z[1]) 79 #endif 80 #if !defined(MAGICKCORE_HAVE_CREAL) 81 #define creal(z) (z[0]) 82 #endif 83 #endif 84 85 /* 87 Typedef declarations. 88 */ 89 typedef struct _FourierInfo 90 { 91 PixelChannel 92 channel; 93 94 MagickBooleanType 95 modulus; 96 97 size_t 98 width, 99 height; 100 101 ssize_t 102 center; 103 } FourierInfo; 104 105 /* 107 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 108 % % 109 % % 110 % % 111 % C o m p l e x I m a g e s % 112 % % 113 % % 114 % % 115 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 116 % 117 % ComplexImages() performs complex mathematics on an image sequence. 118 % 119 % The format of the ComplexImages method is: 120 % 121 % MagickBooleanType ComplexImages(Image *images,const ComplexOperator op, 122 % ExceptionInfo *exception) 123 % 124 % A description of each parameter follows: 125 % 126 % o image: the image. 127 % 128 % o op: A complex operator. 129 % 130 % o exception: return any errors or warnings in this structure. 131 % 132 */ 133 MagickExport Image *ComplexImages(const Image *images,const ComplexOperator op, 134 ExceptionInfo *exception) 135 { 136 #define ComplexImageTag "Complex/Image" 137 138 CacheView 139 *Ai_view, 140 *Ar_view, 141 *Bi_view, 142 *Br_view, 143 *Ci_view, 144 *Cr_view; 145 146 const char 147 *artifact; 148 149 const Image 150 *Ai_image, 151 *Ar_image, 152 *Bi_image, 153 *Br_image; 154 155 double 156 snr; 157 158 Image 159 *Ci_image, 160 *complex_images, 161 *Cr_image, 162 *image; 163 164 MagickBooleanType 165 status; 166 167 MagickOffsetType 168 progress; 169 170 ssize_t 171 y; 172 173 assert(images != (Image *) NULL); 174 assert(images->signature == MagickCoreSignature); 175 if (images->debug != MagickFalse) 176 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",images->filename); 177 assert(exception != (ExceptionInfo *) NULL); 178 assert(exception->signature == MagickCoreSignature); 179 if (images->next == (Image *) NULL) 180 { 181 (void) ThrowMagickException(exception,GetMagickModule(),ImageError, 182 "ImageSequenceRequired","`%s'",images->filename); 183 return((Image *) NULL); 184 } 185 image=CloneImage(images,images->columns,images->rows,MagickTrue,exception); 186 if (image == (Image *) NULL) 187 return((Image *) NULL); 188 if (SetImageStorageClass(image,DirectClass,exception) == MagickFalse) 189 { 190 image=DestroyImageList(image); 191 return(image); 192 } 193 image->depth=32UL; 194 complex_images=NewImageList(); 195 AppendImageToList(&complex_images,image); 196 image=CloneImage(images,images->columns,images->rows,MagickTrue,exception); 197 if (image == (Image *) NULL) 198 { 199 complex_images=DestroyImageList(complex_images); 200 return(complex_images); 201 } 202 AppendImageToList(&complex_images,image); 203 /* 204 Apply complex mathematics to image pixels. 205 */ 206 artifact=GetImageArtifact(image,"complex:snr"); 207 snr=0.0; 208 if (artifact != (const char *) NULL) 209 snr=StringToDouble(artifact,(char **) NULL); 210 Ar_image=images; 211 Ai_image=images->next; 212 Br_image=images; 213 Bi_image=images->next; 214 if ((images->next->next != (Image *) NULL) && 215 (images->next->next->next != (Image *) NULL)) 216 { 217 Br_image=images->next->next; 218 Bi_image=images->next->next->next; 219 } 220 Cr_image=complex_images; 221 Ci_image=complex_images->next; 222 Ar_view=AcquireVirtualCacheView(Ar_image,exception); 223 Ai_view=AcquireVirtualCacheView(Ai_image,exception); 224 Br_view=AcquireVirtualCacheView(Br_image,exception); 225 Bi_view=AcquireVirtualCacheView(Bi_image,exception); 226 Cr_view=AcquireAuthenticCacheView(Cr_image,exception); 227 Ci_view=AcquireAuthenticCacheView(Ci_image,exception); 228 status=MagickTrue; 229 progress=0; 230 #if defined(MAGICKCORE_OPENMP_SUPPORT) 231 #pragma omp parallel for schedule(static,4) shared(progress,status) \ 232 magick_threads(images,complex_images,images->rows,1L) 233 #endif 234 for (y=0; y < (ssize_t) images->rows; y++) 235 { 236 register const Quantum 237 *magick_restrict Ai, 238 *magick_restrict Ar, 239 *magick_restrict Bi, 240 *magick_restrict Br; 241 242 register Quantum 243 *magick_restrict Ci, 244 *magick_restrict Cr; 245 246 register ssize_t 247 x; 248 249 if (status == MagickFalse) 250 continue; 251 Ar=GetCacheViewVirtualPixels(Ar_view,0,y,Ar_image->columns,1,exception); 252 Ai=GetCacheViewVirtualPixels(Ai_view,0,y,Ai_image->columns,1,exception); 253 Br=GetCacheViewVirtualPixels(Br_view,0,y,Br_image->columns,1,exception); 254 Bi=GetCacheViewVirtualPixels(Bi_view,0,y,Bi_image->columns,1,exception); 255 Cr=QueueCacheViewAuthenticPixels(Cr_view,0,y,Cr_image->columns,1,exception); 256 Ci=QueueCacheViewAuthenticPixels(Ci_view,0,y,Ci_image->columns,1,exception); 257 if ((Ar == (const Quantum *) NULL) || (Ai == (const Quantum *) NULL) || 258 (Br == (const Quantum *) NULL) || (Bi == (const Quantum *) NULL) || 259 (Cr == (Quantum *) NULL) || (Ci == (Quantum *) NULL)) 260 { 261 status=MagickFalse; 262 continue; 263 } 264 for (x=0; x < (ssize_t) images->columns; x++) 265 { 266 register ssize_t 267 i; 268 269 for (i=0; i < (ssize_t) GetPixelChannels(images); i++) 270 { 271 switch (op) 272 { 273 case AddComplexOperator: 274 { 275 Cr[i]=Ar[i]+Br[i]; 276 Ci[i]=Ai[i]+Bi[i]; 277 break; 278 } 279 case ConjugateComplexOperator: 280 default: 281 { 282 Cr[i]=Ar[i]; 283 Ci[i]=(-Bi[i]); 284 break; 285 } 286 case DivideComplexOperator: 287 { 288 double 289 gamma; 290 291 gamma=PerceptibleReciprocal(Br[i]*Br[i]+Bi[i]*Bi[i]+snr); 292 Cr[i]=gamma*(Ar[i]*Br[i]+Ai[i]*Bi[i]); 293 Ci[i]=gamma*(Ai[i]*Br[i]-Ar[i]*Bi[i]); 294 break; 295 } 296 case MagnitudePhaseComplexOperator: 297 { 298 Cr[i]=sqrt(Ar[i]*Ar[i]+Ai[i]*Ai[i]); 299 Ci[i]=atan2(Ai[i],Ar[i])/(2.0*MagickPI)+0.5; 300 break; 301 } 302 case MultiplyComplexOperator: 303 { 304 Cr[i]=QuantumScale*(Ar[i]*Br[i]-Ai[i]*Bi[i]); 305 Ci[i]=QuantumScale*(Ai[i]*Br[i]+Ar[i]*Bi[i]); 306 break; 307 } 308 case RealImaginaryComplexOperator: 309 { 310 Cr[i]=Ar[i]*cos(2.0*MagickPI*(Ai[i]-0.5)); 311 Ci[i]=Ar[i]*sin(2.0*MagickPI*(Ai[i]-0.5)); 312 break; 313 } 314 case SubtractComplexOperator: 315 { 316 Cr[i]=Ar[i]-Br[i]; 317 Ci[i]=Ai[i]-Bi[i]; 318 break; 319 } 320 } 321 } 322 Ar+=GetPixelChannels(Ar_image); 323 Ai+=GetPixelChannels(Ai_image); 324 Br+=GetPixelChannels(Br_image); 325 Bi+=GetPixelChannels(Bi_image); 326 Cr+=GetPixelChannels(Cr_image); 327 Ci+=GetPixelChannels(Ci_image); 328 } 329 if (SyncCacheViewAuthenticPixels(Ci_view,exception) == MagickFalse) 330 status=MagickFalse; 331 if (SyncCacheViewAuthenticPixels(Cr_view,exception) == MagickFalse) 332 status=MagickFalse; 333 if (images->progress_monitor != (MagickProgressMonitor) NULL) 334 { 335 MagickBooleanType 336 proceed; 337 338 #if defined(MAGICKCORE_OPENMP_SUPPORT) 339 #pragma omp critical (MagickCore_ComplexImages) 340 #endif 341 proceed=SetImageProgress(images,ComplexImageTag,progress++, 342 images->rows); 343 if (proceed == MagickFalse) 344 status=MagickFalse; 345 } 346 } 347 Cr_view=DestroyCacheView(Cr_view); 348 Ci_view=DestroyCacheView(Ci_view); 349 Br_view=DestroyCacheView(Br_view); 350 Bi_view=DestroyCacheView(Bi_view); 351 Ar_view=DestroyCacheView(Ar_view); 352 Ai_view=DestroyCacheView(Ai_view); 353 if (status == MagickFalse) 354 complex_images=DestroyImageList(complex_images); 355 return(complex_images); 356 } 357 358 /* 360 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 361 % % 362 % % 363 % % 364 % F o r w a r d F o u r i e r T r a n s f o r m I m a g e % 365 % % 366 % % 367 % % 368 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 369 % 370 % ForwardFourierTransformImage() implements the discrete Fourier transform 371 % (DFT) of the image either as a magnitude / phase or real / imaginary image 372 % pair. 373 % 374 % The format of the ForwadFourierTransformImage method is: 375 % 376 % Image *ForwardFourierTransformImage(const Image *image, 377 % const MagickBooleanType modulus,ExceptionInfo *exception) 378 % 379 % A description of each parameter follows: 380 % 381 % o image: the image. 382 % 383 % o modulus: if true, return as transform as a magnitude / phase pair 384 % otherwise a real / imaginary image pair. 385 % 386 % o exception: return any errors or warnings in this structure. 387 % 388 */ 389 390 #if defined(MAGICKCORE_FFTW_DELEGATE) 391 392 static MagickBooleanType RollFourier(const size_t width,const size_t height, 393 const ssize_t x_offset,const ssize_t y_offset,double *roll_pixels) 394 { 395 double 396 *source_pixels; 397 398 MemoryInfo 399 *source_info; 400 401 register ssize_t 402 i, 403 x; 404 405 ssize_t 406 u, 407 v, 408 y; 409 410 /* 411 Move zero frequency (DC, average color) from (0,0) to (width/2,height/2). 412 */ 413 source_info=AcquireVirtualMemory(width,height*sizeof(*source_pixels)); 414 if (source_info == (MemoryInfo *) NULL) 415 return(MagickFalse); 416 source_pixels=(double *) GetVirtualMemoryBlob(source_info); 417 i=0L; 418 for (y=0L; y < (ssize_t) height; y++) 419 { 420 if (y_offset < 0L) 421 v=((y+y_offset) < 0L) ? y+y_offset+(ssize_t) height : y+y_offset; 422 else 423 v=((y+y_offset) > ((ssize_t) height-1L)) ? y+y_offset-(ssize_t) height : 424 y+y_offset; 425 for (x=0L; x < (ssize_t) width; x++) 426 { 427 if (x_offset < 0L) 428 u=((x+x_offset) < 0L) ? x+x_offset+(ssize_t) width : x+x_offset; 429 else 430 u=((x+x_offset) > ((ssize_t) width-1L)) ? x+x_offset-(ssize_t) width : 431 x+x_offset; 432 source_pixels[v*width+u]=roll_pixels[i++]; 433 } 434 } 435 (void) CopyMagickMemory(roll_pixels,source_pixels,height*width* 436 sizeof(*source_pixels)); 437 source_info=RelinquishVirtualMemory(source_info); 438 return(MagickTrue); 439 } 440 441 static MagickBooleanType ForwardQuadrantSwap(const size_t width, 442 const size_t height,double *source_pixels,double *forward_pixels) 443 { 444 MagickBooleanType 445 status; 446 447 register ssize_t 448 x; 449 450 ssize_t 451 center, 452 y; 453 454 /* 455 Swap quadrants. 456 */ 457 center=(ssize_t) (width/2L)+1L; 458 status=RollFourier((size_t) center,height,0L,(ssize_t) height/2L, 459 source_pixels); 460 if (status == MagickFalse) 461 return(MagickFalse); 462 for (y=0L; y < (ssize_t) height; y++) 463 for (x=0L; x < (ssize_t) (width/2L); x++) 464 forward_pixels[y*width+x+width/2L]=source_pixels[y*center+x]; 465 for (y=1; y < (ssize_t) height; y++) 466 for (x=0L; x < (ssize_t) (width/2L); x++) 467 forward_pixels[(height-y)*width+width/2L-x-1L]= 468 source_pixels[y*center+x+1L]; 469 for (x=0L; x < (ssize_t) (width/2L); x++) 470 forward_pixels[width/2L-x-1L]=source_pixels[x+1L]; 471 return(MagickTrue); 472 } 473 474 static void CorrectPhaseLHS(const size_t width,const size_t height, 475 double *fourier_pixels) 476 { 477 register ssize_t 478 x; 479 480 ssize_t 481 y; 482 483 for (y=0L; y < (ssize_t) height; y++) 484 for (x=0L; x < (ssize_t) (width/2L); x++) 485 fourier_pixels[y*width+x]*=(-1.0); 486 } 487 488 static MagickBooleanType ForwardFourier(const FourierInfo *fourier_info, 489 Image *image,double *magnitude,double *phase,ExceptionInfo *exception) 490 { 491 CacheView 492 *magnitude_view, 493 *phase_view; 494 495 double 496 *magnitude_pixels, 497 *phase_pixels; 498 499 Image 500 *magnitude_image, 501 *phase_image; 502 503 MagickBooleanType 504 status; 505 506 MemoryInfo 507 *magnitude_info, 508 *phase_info; 509 510 register Quantum 511 *q; 512 513 register ssize_t 514 x; 515 516 ssize_t 517 i, 518 y; 519 520 magnitude_image=GetFirstImageInList(image); 521 phase_image=GetNextImageInList(image); 522 if (phase_image == (Image *) NULL) 523 { 524 (void) ThrowMagickException(exception,GetMagickModule(),ImageError, 525 "ImageSequenceRequired","`%s'",image->filename); 526 return(MagickFalse); 527 } 528 /* 529 Create "Fourier Transform" image from constituent arrays. 530 */ 531 magnitude_info=AcquireVirtualMemory((size_t) fourier_info->width, 532 fourier_info->height*sizeof(*magnitude_pixels)); 533 phase_info=AcquireVirtualMemory((size_t) fourier_info->width, 534 fourier_info->height*sizeof(*phase_pixels)); 535 if ((magnitude_info == (MemoryInfo *) NULL) || 536 (phase_info == (MemoryInfo *) NULL)) 537 { 538 if (phase_info != (MemoryInfo *) NULL) 539 phase_info=RelinquishVirtualMemory(phase_info); 540 if (magnitude_info != (MemoryInfo *) NULL) 541 magnitude_info=RelinquishVirtualMemory(magnitude_info); 542 (void) ThrowMagickException(exception,GetMagickModule(), 543 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 544 return(MagickFalse); 545 } 546 magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info); 547 (void) ResetMagickMemory(magnitude_pixels,0,fourier_info->width* 548 fourier_info->height*sizeof(*magnitude_pixels)); 549 phase_pixels=(double *) GetVirtualMemoryBlob(phase_info); 550 (void) ResetMagickMemory(phase_pixels,0,fourier_info->width* 551 fourier_info->height*sizeof(*phase_pixels)); 552 status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height, 553 magnitude,magnitude_pixels); 554 if (status != MagickFalse) 555 status=ForwardQuadrantSwap(fourier_info->width,fourier_info->height,phase, 556 phase_pixels); 557 CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels); 558 if (fourier_info->modulus != MagickFalse) 559 { 560 i=0L; 561 for (y=0L; y < (ssize_t) fourier_info->height; y++) 562 for (x=0L; x < (ssize_t) fourier_info->width; x++) 563 { 564 phase_pixels[i]/=(2.0*MagickPI); 565 phase_pixels[i]+=0.5; 566 i++; 567 } 568 } 569 magnitude_view=AcquireAuthenticCacheView(magnitude_image,exception); 570 i=0L; 571 for (y=0L; y < (ssize_t) fourier_info->height; y++) 572 { 573 q=GetCacheViewAuthenticPixels(magnitude_view,0L,y,fourier_info->width,1UL, 574 exception); 575 if (q == (Quantum *) NULL) 576 break; 577 for (x=0L; x < (ssize_t) fourier_info->width; x++) 578 { 579 switch (fourier_info->channel) 580 { 581 case RedPixelChannel: 582 default: 583 { 584 SetPixelRed(magnitude_image,ClampToQuantum(QuantumRange* 585 magnitude_pixels[i]),q); 586 break; 587 } 588 case GreenPixelChannel: 589 { 590 SetPixelGreen(magnitude_image,ClampToQuantum(QuantumRange* 591 magnitude_pixels[i]),q); 592 break; 593 } 594 case BluePixelChannel: 595 { 596 SetPixelBlue(magnitude_image,ClampToQuantum(QuantumRange* 597 magnitude_pixels[i]),q); 598 break; 599 } 600 case BlackPixelChannel: 601 { 602 SetPixelBlack(magnitude_image,ClampToQuantum(QuantumRange* 603 magnitude_pixels[i]),q); 604 break; 605 } 606 case AlphaPixelChannel: 607 { 608 SetPixelAlpha(magnitude_image,ClampToQuantum(QuantumRange* 609 magnitude_pixels[i]),q); 610 break; 611 } 612 } 613 i++; 614 q+=GetPixelChannels(magnitude_image); 615 } 616 status=SyncCacheViewAuthenticPixels(magnitude_view,exception); 617 if (status == MagickFalse) 618 break; 619 } 620 magnitude_view=DestroyCacheView(magnitude_view); 621 i=0L; 622 phase_view=AcquireAuthenticCacheView(phase_image,exception); 623 for (y=0L; y < (ssize_t) fourier_info->height; y++) 624 { 625 q=GetCacheViewAuthenticPixels(phase_view,0L,y,fourier_info->width,1UL, 626 exception); 627 if (q == (Quantum *) NULL) 628 break; 629 for (x=0L; x < (ssize_t) fourier_info->width; x++) 630 { 631 switch (fourier_info->channel) 632 { 633 case RedPixelChannel: 634 default: 635 { 636 SetPixelRed(phase_image,ClampToQuantum(QuantumRange* 637 phase_pixels[i]),q); 638 break; 639 } 640 case GreenPixelChannel: 641 { 642 SetPixelGreen(phase_image,ClampToQuantum(QuantumRange* 643 phase_pixels[i]),q); 644 break; 645 } 646 case BluePixelChannel: 647 { 648 SetPixelBlue(phase_image,ClampToQuantum(QuantumRange* 649 phase_pixels[i]),q); 650 break; 651 } 652 case BlackPixelChannel: 653 { 654 SetPixelBlack(phase_image,ClampToQuantum(QuantumRange* 655 phase_pixels[i]),q); 656 break; 657 } 658 case AlphaPixelChannel: 659 { 660 SetPixelAlpha(phase_image,ClampToQuantum(QuantumRange* 661 phase_pixels[i]),q); 662 break; 663 } 664 } 665 i++; 666 q+=GetPixelChannels(phase_image); 667 } 668 status=SyncCacheViewAuthenticPixels(phase_view,exception); 669 if (status == MagickFalse) 670 break; 671 } 672 phase_view=DestroyCacheView(phase_view); 673 phase_info=RelinquishVirtualMemory(phase_info); 674 magnitude_info=RelinquishVirtualMemory(magnitude_info); 675 return(status); 676 } 677 678 static MagickBooleanType ForwardFourierTransform(FourierInfo *fourier_info, 679 const Image *image,double *magnitude_pixels,double *phase_pixels, 680 ExceptionInfo *exception) 681 { 682 CacheView 683 *image_view; 684 685 const char 686 *value; 687 688 double 689 *source_pixels; 690 691 fftw_complex 692 *forward_pixels; 693 694 fftw_plan 695 fftw_r2c_plan; 696 697 MemoryInfo 698 *forward_info, 699 *source_info; 700 701 register const Quantum 702 *p; 703 704 register ssize_t 705 i, 706 x; 707 708 ssize_t 709 y; 710 711 /* 712 Generate the forward Fourier transform. 713 */ 714 source_info=AcquireVirtualMemory((size_t) fourier_info->width, 715 fourier_info->height*sizeof(*source_pixels)); 716 if (source_info == (MemoryInfo *) NULL) 717 { 718 (void) ThrowMagickException(exception,GetMagickModule(), 719 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 720 return(MagickFalse); 721 } 722 source_pixels=(double *) GetVirtualMemoryBlob(source_info); 723 ResetMagickMemory(source_pixels,0,fourier_info->width*fourier_info->height* 724 sizeof(*source_pixels)); 725 i=0L; 726 image_view=AcquireVirtualCacheView(image,exception); 727 for (y=0L; y < (ssize_t) fourier_info->height; y++) 728 { 729 p=GetCacheViewVirtualPixels(image_view,0L,y,fourier_info->width,1UL, 730 exception); 731 if (p == (const Quantum *) NULL) 732 break; 733 for (x=0L; x < (ssize_t) fourier_info->width; x++) 734 { 735 switch (fourier_info->channel) 736 { 737 case RedPixelChannel: 738 default: 739 { 740 source_pixels[i]=QuantumScale*GetPixelRed(image,p); 741 break; 742 } 743 case GreenPixelChannel: 744 { 745 source_pixels[i]=QuantumScale*GetPixelGreen(image,p); 746 break; 747 } 748 case BluePixelChannel: 749 { 750 source_pixels[i]=QuantumScale*GetPixelBlue(image,p); 751 break; 752 } 753 case BlackPixelChannel: 754 { 755 source_pixels[i]=QuantumScale*GetPixelBlack(image,p); 756 break; 757 } 758 case AlphaPixelChannel: 759 { 760 source_pixels[i]=QuantumScale*GetPixelAlpha(image,p); 761 break; 762 } 763 } 764 i++; 765 p+=GetPixelChannels(image); 766 } 767 } 768 image_view=DestroyCacheView(image_view); 769 forward_info=AcquireVirtualMemory((size_t) fourier_info->width, 770 (fourier_info->height/2+1)*sizeof(*forward_pixels)); 771 if (forward_info == (MemoryInfo *) NULL) 772 { 773 (void) ThrowMagickException(exception,GetMagickModule(), 774 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 775 source_info=(MemoryInfo *) RelinquishVirtualMemory(source_info); 776 return(MagickFalse); 777 } 778 forward_pixels=(fftw_complex *) GetVirtualMemoryBlob(forward_info); 779 #if defined(MAGICKCORE_OPENMP_SUPPORT) 780 #pragma omp critical (MagickCore_ForwardFourierTransform) 781 #endif 782 fftw_r2c_plan=fftw_plan_dft_r2c_2d(fourier_info->width,fourier_info->height, 783 source_pixels,forward_pixels,FFTW_ESTIMATE); 784 fftw_execute_dft_r2c(fftw_r2c_plan,source_pixels,forward_pixels); 785 fftw_destroy_plan(fftw_r2c_plan); 786 source_info=(MemoryInfo *) RelinquishVirtualMemory(source_info); 787 value=GetImageArtifact(image,"fourier:normalize"); 788 if ((value == (const char *) NULL) || (LocaleCompare(value,"forward") == 0)) 789 { 790 double 791 gamma; 792 793 /* 794 Normalize fourier transform. 795 */ 796 i=0L; 797 gamma=PerceptibleReciprocal((double) fourier_info->width* 798 fourier_info->height); 799 for (y=0L; y < (ssize_t) fourier_info->height; y++) 800 for (x=0L; x < (ssize_t) fourier_info->center; x++) 801 { 802 #if defined(MAGICKCORE_HAVE_COMPLEX_H) 803 forward_pixels[i]*=gamma; 804 #else 805 forward_pixels[i][0]*=gamma; 806 forward_pixels[i][1]*=gamma; 807 #endif 808 i++; 809 } 810 } 811 /* 812 Generate magnitude and phase (or real and imaginary). 813 */ 814 i=0L; 815 if (fourier_info->modulus != MagickFalse) 816 for (y=0L; y < (ssize_t) fourier_info->height; y++) 817 for (x=0L; x < (ssize_t) fourier_info->center; x++) 818 { 819 magnitude_pixels[i]=cabs(forward_pixels[i]); 820 phase_pixels[i]=carg(forward_pixels[i]); 821 i++; 822 } 823 else 824 for (y=0L; y < (ssize_t) fourier_info->height; y++) 825 for (x=0L; x < (ssize_t) fourier_info->center; x++) 826 { 827 magnitude_pixels[i]=creal(forward_pixels[i]); 828 phase_pixels[i]=cimag(forward_pixels[i]); 829 i++; 830 } 831 forward_info=(MemoryInfo *) RelinquishVirtualMemory(forward_info); 832 return(MagickTrue); 833 } 834 835 static MagickBooleanType ForwardFourierTransformChannel(const Image *image, 836 const PixelChannel channel,const MagickBooleanType modulus, 837 Image *fourier_image,ExceptionInfo *exception) 838 { 839 double 840 *magnitude_pixels, 841 *phase_pixels; 842 843 FourierInfo 844 fourier_info; 845 846 MagickBooleanType 847 status; 848 849 MemoryInfo 850 *magnitude_info, 851 *phase_info; 852 853 fourier_info.width=image->columns; 854 fourier_info.height=image->rows; 855 if ((image->columns != image->rows) || ((image->columns % 2) != 0) || 856 ((image->rows % 2) != 0)) 857 { 858 size_t extent=image->columns < image->rows ? image->rows : image->columns; 859 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent; 860 } 861 fourier_info.height=fourier_info.width; 862 fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L; 863 fourier_info.channel=channel; 864 fourier_info.modulus=modulus; 865 magnitude_info=AcquireVirtualMemory((size_t) fourier_info.width, 866 (fourier_info.height/2+1)*sizeof(*magnitude_pixels)); 867 phase_info=AcquireVirtualMemory((size_t) fourier_info.width, 868 (fourier_info.height/2+1)*sizeof(*phase_pixels)); 869 if ((magnitude_info == (MemoryInfo *) NULL) || 870 (phase_info == (MemoryInfo *) NULL)) 871 { 872 if (phase_info != (MemoryInfo *) NULL) 873 phase_info=RelinquishVirtualMemory(phase_info); 874 if (magnitude_info == (MemoryInfo *) NULL) 875 magnitude_info=RelinquishVirtualMemory(magnitude_info); 876 (void) ThrowMagickException(exception,GetMagickModule(), 877 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 878 return(MagickFalse); 879 } 880 magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info); 881 phase_pixels=(double *) GetVirtualMemoryBlob(phase_info); 882 status=ForwardFourierTransform(&fourier_info,image,magnitude_pixels, 883 phase_pixels,exception); 884 if (status != MagickFalse) 885 status=ForwardFourier(&fourier_info,fourier_image,magnitude_pixels, 886 phase_pixels,exception); 887 phase_info=RelinquishVirtualMemory(phase_info); 888 magnitude_info=RelinquishVirtualMemory(magnitude_info); 889 return(status); 890 } 891 #endif 892 893 MagickExport Image *ForwardFourierTransformImage(const Image *image, 894 const MagickBooleanType modulus,ExceptionInfo *exception) 895 { 896 Image 897 *fourier_image; 898 899 fourier_image=NewImageList(); 900 #if !defined(MAGICKCORE_FFTW_DELEGATE) 901 (void) modulus; 902 (void) ThrowMagickException(exception,GetMagickModule(), 903 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)", 904 image->filename); 905 #else 906 { 907 Image 908 *magnitude_image; 909 910 size_t 911 height, 912 width; 913 914 width=image->columns; 915 height=image->rows; 916 if ((image->columns != image->rows) || ((image->columns % 2) != 0) || 917 ((image->rows % 2) != 0)) 918 { 919 size_t extent=image->columns < image->rows ? image->rows : 920 image->columns; 921 width=(extent & 0x01) == 1 ? extent+1UL : extent; 922 } 923 height=width; 924 magnitude_image=CloneImage(image,width,height,MagickTrue,exception); 925 if (magnitude_image != (Image *) NULL) 926 { 927 Image 928 *phase_image; 929 930 magnitude_image->storage_class=DirectClass; 931 magnitude_image->depth=32UL; 932 phase_image=CloneImage(image,width,height,MagickTrue,exception); 933 if (phase_image == (Image *) NULL) 934 magnitude_image=DestroyImage(magnitude_image); 935 else 936 { 937 MagickBooleanType 938 is_gray, 939 status; 940 941 phase_image->storage_class=DirectClass; 942 phase_image->depth=32UL; 943 AppendImageToList(&fourier_image,magnitude_image); 944 AppendImageToList(&fourier_image,phase_image); 945 status=MagickTrue; 946 is_gray=IsImageGray(image); 947 #if defined(MAGICKCORE_OPENMP_SUPPORT) 948 #pragma omp parallel sections 949 #endif 950 { 951 #if defined(MAGICKCORE_OPENMP_SUPPORT) 952 #pragma omp section 953 #endif 954 { 955 MagickBooleanType 956 thread_status; 957 958 if (is_gray != MagickFalse) 959 thread_status=ForwardFourierTransformChannel(image, 960 GrayPixelChannel,modulus,fourier_image,exception); 961 else 962 thread_status=ForwardFourierTransformChannel(image, 963 RedPixelChannel,modulus,fourier_image,exception); 964 if (thread_status == MagickFalse) 965 status=thread_status; 966 } 967 #if defined(MAGICKCORE_OPENMP_SUPPORT) 968 #pragma omp section 969 #endif 970 { 971 MagickBooleanType 972 thread_status; 973 974 thread_status=MagickTrue; 975 if (is_gray == MagickFalse) 976 thread_status=ForwardFourierTransformChannel(image, 977 GreenPixelChannel,modulus,fourier_image,exception); 978 if (thread_status == MagickFalse) 979 status=thread_status; 980 } 981 #if defined(MAGICKCORE_OPENMP_SUPPORT) 982 #pragma omp section 983 #endif 984 { 985 MagickBooleanType 986 thread_status; 987 988 thread_status=MagickTrue; 989 if (is_gray == MagickFalse) 990 thread_status=ForwardFourierTransformChannel(image, 991 BluePixelChannel,modulus,fourier_image,exception); 992 if (thread_status == MagickFalse) 993 status=thread_status; 994 } 995 #if defined(MAGICKCORE_OPENMP_SUPPORT) 996 #pragma omp section 997 #endif 998 { 999 MagickBooleanType 1000 thread_status; 1001 1002 thread_status=MagickTrue; 1003 if (image->colorspace == CMYKColorspace) 1004 thread_status=ForwardFourierTransformChannel(image, 1005 BlackPixelChannel,modulus,fourier_image,exception); 1006 if (thread_status == MagickFalse) 1007 status=thread_status; 1008 } 1009 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1010 #pragma omp section 1011 #endif 1012 { 1013 MagickBooleanType 1014 thread_status; 1015 1016 thread_status=MagickTrue; 1017 if (image->alpha_trait != UndefinedPixelTrait) 1018 thread_status=ForwardFourierTransformChannel(image, 1019 AlphaPixelChannel,modulus,fourier_image,exception); 1020 if (thread_status == MagickFalse) 1021 status=thread_status; 1022 } 1023 } 1024 if (status == MagickFalse) 1025 fourier_image=DestroyImageList(fourier_image); 1026 fftw_cleanup(); 1027 } 1028 } 1029 } 1030 #endif 1031 return(fourier_image); 1032 } 1033 1034 /* 1036 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1037 % % 1038 % % 1039 % % 1040 % I n v e r s e F o u r i e r T r a n s f o r m I m a g e % 1041 % % 1042 % % 1043 % % 1044 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 1045 % 1046 % InverseFourierTransformImage() implements the inverse discrete Fourier 1047 % transform (DFT) of the image either as a magnitude / phase or real / 1048 % imaginary image pair. 1049 % 1050 % The format of the InverseFourierTransformImage method is: 1051 % 1052 % Image *InverseFourierTransformImage(const Image *magnitude_image, 1053 % const Image *phase_image,const MagickBooleanType modulus, 1054 % ExceptionInfo *exception) 1055 % 1056 % A description of each parameter follows: 1057 % 1058 % o magnitude_image: the magnitude or real image. 1059 % 1060 % o phase_image: the phase or imaginary image. 1061 % 1062 % o modulus: if true, return transform as a magnitude / phase pair 1063 % otherwise a real / imaginary image pair. 1064 % 1065 % o exception: return any errors or warnings in this structure. 1066 % 1067 */ 1068 1069 #if defined(MAGICKCORE_FFTW_DELEGATE) 1070 static MagickBooleanType InverseQuadrantSwap(const size_t width, 1071 const size_t height,const double *source,double *destination) 1072 { 1073 register ssize_t 1074 x; 1075 1076 ssize_t 1077 center, 1078 y; 1079 1080 /* 1081 Swap quadrants. 1082 */ 1083 center=(ssize_t) (width/2L)+1L; 1084 for (y=1L; y < (ssize_t) height; y++) 1085 for (x=0L; x < (ssize_t) (width/2L+1L); x++) 1086 destination[(height-y)*center-x+width/2L]=source[y*width+x]; 1087 for (y=0L; y < (ssize_t) height; y++) 1088 destination[y*center]=source[y*width+width/2L]; 1089 for (x=0L; x < center; x++) 1090 destination[x]=source[center-x-1L]; 1091 return(RollFourier(center,height,0L,(ssize_t) height/-2L,destination)); 1092 } 1093 1094 static MagickBooleanType InverseFourier(FourierInfo *fourier_info, 1095 const Image *magnitude_image,const Image *phase_image, 1096 fftw_complex *fourier_pixels,ExceptionInfo *exception) 1097 { 1098 CacheView 1099 *magnitude_view, 1100 *phase_view; 1101 1102 double 1103 *inverse_pixels, 1104 *magnitude_pixels, 1105 *phase_pixels; 1106 1107 MagickBooleanType 1108 status; 1109 1110 MemoryInfo 1111 *inverse_info, 1112 *magnitude_info, 1113 *phase_info; 1114 1115 register const Quantum 1116 *p; 1117 1118 register ssize_t 1119 i, 1120 x; 1121 1122 ssize_t 1123 y; 1124 1125 /* 1126 Inverse fourier - read image and break down into a double array. 1127 */ 1128 magnitude_info=AcquireVirtualMemory((size_t) fourier_info->width, 1129 fourier_info->height*sizeof(*magnitude_pixels)); 1130 phase_info=AcquireVirtualMemory((size_t) fourier_info->width, 1131 fourier_info->height*sizeof(*phase_pixels)); 1132 inverse_info=AcquireVirtualMemory((size_t) fourier_info->width, 1133 (fourier_info->height/2+1)*sizeof(*inverse_pixels)); 1134 if ((magnitude_info == (MemoryInfo *) NULL) || 1135 (phase_info == (MemoryInfo *) NULL) || 1136 (inverse_info == (MemoryInfo *) NULL)) 1137 { 1138 if (magnitude_info != (MemoryInfo *) NULL) 1139 magnitude_info=RelinquishVirtualMemory(magnitude_info); 1140 if (phase_info != (MemoryInfo *) NULL) 1141 phase_info=RelinquishVirtualMemory(phase_info); 1142 if (inverse_info != (MemoryInfo *) NULL) 1143 inverse_info=RelinquishVirtualMemory(inverse_info); 1144 (void) ThrowMagickException(exception,GetMagickModule(), 1145 ResourceLimitError,"MemoryAllocationFailed","`%s'", 1146 magnitude_image->filename); 1147 return(MagickFalse); 1148 } 1149 magnitude_pixels=(double *) GetVirtualMemoryBlob(magnitude_info); 1150 phase_pixels=(double *) GetVirtualMemoryBlob(phase_info); 1151 inverse_pixels=(double *) GetVirtualMemoryBlob(inverse_info); 1152 i=0L; 1153 magnitude_view=AcquireVirtualCacheView(magnitude_image,exception); 1154 for (y=0L; y < (ssize_t) fourier_info->height; y++) 1155 { 1156 p=GetCacheViewVirtualPixels(magnitude_view,0L,y,fourier_info->width,1UL, 1157 exception); 1158 if (p == (const Quantum *) NULL) 1159 break; 1160 for (x=0L; x < (ssize_t) fourier_info->width; x++) 1161 { 1162 switch (fourier_info->channel) 1163 { 1164 case RedPixelChannel: 1165 default: 1166 { 1167 magnitude_pixels[i]=QuantumScale*GetPixelRed(magnitude_image,p); 1168 break; 1169 } 1170 case GreenPixelChannel: 1171 { 1172 magnitude_pixels[i]=QuantumScale*GetPixelGreen(magnitude_image,p); 1173 break; 1174 } 1175 case BluePixelChannel: 1176 { 1177 magnitude_pixels[i]=QuantumScale*GetPixelBlue(magnitude_image,p); 1178 break; 1179 } 1180 case BlackPixelChannel: 1181 { 1182 magnitude_pixels[i]=QuantumScale*GetPixelBlack(magnitude_image,p); 1183 break; 1184 } 1185 case AlphaPixelChannel: 1186 { 1187 magnitude_pixels[i]=QuantumScale*GetPixelAlpha(magnitude_image,p); 1188 break; 1189 } 1190 } 1191 i++; 1192 p+=GetPixelChannels(magnitude_image); 1193 } 1194 } 1195 magnitude_view=DestroyCacheView(magnitude_view); 1196 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height, 1197 magnitude_pixels,inverse_pixels); 1198 (void) CopyMagickMemory(magnitude_pixels,inverse_pixels,fourier_info->height* 1199 fourier_info->center*sizeof(*magnitude_pixels)); 1200 i=0L; 1201 phase_view=AcquireVirtualCacheView(phase_image,exception); 1202 for (y=0L; y < (ssize_t) fourier_info->height; y++) 1203 { 1204 p=GetCacheViewVirtualPixels(phase_view,0,y,fourier_info->width,1, 1205 exception); 1206 if (p == (const Quantum *) NULL) 1207 break; 1208 for (x=0L; x < (ssize_t) fourier_info->width; x++) 1209 { 1210 switch (fourier_info->channel) 1211 { 1212 case RedPixelChannel: 1213 default: 1214 { 1215 phase_pixels[i]=QuantumScale*GetPixelRed(phase_image,p); 1216 break; 1217 } 1218 case GreenPixelChannel: 1219 { 1220 phase_pixels[i]=QuantumScale*GetPixelGreen(phase_image,p); 1221 break; 1222 } 1223 case BluePixelChannel: 1224 { 1225 phase_pixels[i]=QuantumScale*GetPixelBlue(phase_image,p); 1226 break; 1227 } 1228 case BlackPixelChannel: 1229 { 1230 phase_pixels[i]=QuantumScale*GetPixelBlack(phase_image,p); 1231 break; 1232 } 1233 case AlphaPixelChannel: 1234 { 1235 phase_pixels[i]=QuantumScale*GetPixelAlpha(phase_image,p); 1236 break; 1237 } 1238 } 1239 i++; 1240 p+=GetPixelChannels(phase_image); 1241 } 1242 } 1243 if (fourier_info->modulus != MagickFalse) 1244 { 1245 i=0L; 1246 for (y=0L; y < (ssize_t) fourier_info->height; y++) 1247 for (x=0L; x < (ssize_t) fourier_info->width; x++) 1248 { 1249 phase_pixels[i]-=0.5; 1250 phase_pixels[i]*=(2.0*MagickPI); 1251 i++; 1252 } 1253 } 1254 phase_view=DestroyCacheView(phase_view); 1255 CorrectPhaseLHS(fourier_info->width,fourier_info->height,phase_pixels); 1256 if (status != MagickFalse) 1257 status=InverseQuadrantSwap(fourier_info->width,fourier_info->height, 1258 phase_pixels,inverse_pixels); 1259 (void) CopyMagickMemory(phase_pixels,inverse_pixels,fourier_info->height* 1260 fourier_info->center*sizeof(*phase_pixels)); 1261 inverse_info=RelinquishVirtualMemory(inverse_info); 1262 /* 1263 Merge two sets. 1264 */ 1265 i=0L; 1266 if (fourier_info->modulus != MagickFalse) 1267 for (y=0L; y < (ssize_t) fourier_info->height; y++) 1268 for (x=0L; x < (ssize_t) fourier_info->center; x++) 1269 { 1270 #if defined(MAGICKCORE_HAVE_COMPLEX_H) 1271 fourier_pixels[i]=magnitude_pixels[i]*cos(phase_pixels[i])+I* 1272 magnitude_pixels[i]*sin(phase_pixels[i]); 1273 #else 1274 fourier_pixels[i][0]=magnitude_pixels[i]*cos(phase_pixels[i]); 1275 fourier_pixels[i][1]=magnitude_pixels[i]*sin(phase_pixels[i]); 1276 #endif 1277 i++; 1278 } 1279 else 1280 for (y=0L; y < (ssize_t) fourier_info->height; y++) 1281 for (x=0L; x < (ssize_t) fourier_info->center; x++) 1282 { 1283 #if defined(MAGICKCORE_HAVE_COMPLEX_H) 1284 fourier_pixels[i]=magnitude_pixels[i]+I*phase_pixels[i]; 1285 #else 1286 fourier_pixels[i][0]=magnitude_pixels[i]; 1287 fourier_pixels[i][1]=phase_pixels[i]; 1288 #endif 1289 i++; 1290 } 1291 magnitude_info=RelinquishVirtualMemory(magnitude_info); 1292 phase_info=RelinquishVirtualMemory(phase_info); 1293 return(status); 1294 } 1295 1296 static MagickBooleanType InverseFourierTransform(FourierInfo *fourier_info, 1297 fftw_complex *fourier_pixels,Image *image,ExceptionInfo *exception) 1298 { 1299 CacheView 1300 *image_view; 1301 1302 const char 1303 *value; 1304 1305 double 1306 *source_pixels; 1307 1308 fftw_plan 1309 fftw_c2r_plan; 1310 1311 MemoryInfo 1312 *source_info; 1313 1314 register Quantum 1315 *q; 1316 1317 register ssize_t 1318 i, 1319 x; 1320 1321 ssize_t 1322 y; 1323 1324 source_info=AcquireVirtualMemory((size_t) fourier_info->width, 1325 fourier_info->height*sizeof(*source_pixels)); 1326 if (source_info == (MemoryInfo *) NULL) 1327 { 1328 (void) ThrowMagickException(exception,GetMagickModule(), 1329 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename); 1330 return(MagickFalse); 1331 } 1332 source_pixels=(double *) GetVirtualMemoryBlob(source_info); 1333 value=GetImageArtifact(image,"fourier:normalize"); 1334 if (LocaleCompare(value,"inverse") == 0) 1335 { 1336 double 1337 gamma; 1338 1339 /* 1340 Normalize inverse transform. 1341 */ 1342 i=0L; 1343 gamma=PerceptibleReciprocal((double) fourier_info->width* 1344 fourier_info->height); 1345 for (y=0L; y < (ssize_t) fourier_info->height; y++) 1346 for (x=0L; x < (ssize_t) fourier_info->center; x++) 1347 { 1348 #if defined(MAGICKCORE_HAVE_COMPLEX_H) 1349 fourier_pixels[i]*=gamma; 1350 #else 1351 fourier_pixels[i][0]*=gamma; 1352 fourier_pixels[i][1]*=gamma; 1353 #endif 1354 i++; 1355 } 1356 } 1357 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1358 #pragma omp critical (MagickCore_InverseFourierTransform) 1359 #endif 1360 fftw_c2r_plan=fftw_plan_dft_c2r_2d(fourier_info->width,fourier_info->height, 1361 fourier_pixels,source_pixels,FFTW_ESTIMATE); 1362 fftw_execute_dft_c2r(fftw_c2r_plan,fourier_pixels,source_pixels); 1363 fftw_destroy_plan(fftw_c2r_plan); 1364 i=0L; 1365 image_view=AcquireAuthenticCacheView(image,exception); 1366 for (y=0L; y < (ssize_t) fourier_info->height; y++) 1367 { 1368 if (y >= (ssize_t) image->rows) 1369 break; 1370 q=GetCacheViewAuthenticPixels(image_view,0L,y,fourier_info->width > 1371 image->columns ? image->columns : fourier_info->width,1UL,exception); 1372 if (q == (Quantum *) NULL) 1373 break; 1374 for (x=0L; x < (ssize_t) fourier_info->width; x++) 1375 { 1376 if (x < (ssize_t) image->columns) 1377 switch (fourier_info->channel) 1378 { 1379 case RedPixelChannel: 1380 default: 1381 { 1382 SetPixelRed(image,ClampToQuantum(QuantumRange*source_pixels[i]),q); 1383 break; 1384 } 1385 case GreenPixelChannel: 1386 { 1387 SetPixelGreen(image,ClampToQuantum(QuantumRange*source_pixels[i]), 1388 q); 1389 break; 1390 } 1391 case BluePixelChannel: 1392 { 1393 SetPixelBlue(image,ClampToQuantum(QuantumRange*source_pixels[i]), 1394 q); 1395 break; 1396 } 1397 case BlackPixelChannel: 1398 { 1399 SetPixelBlack(image,ClampToQuantum(QuantumRange*source_pixels[i]), 1400 q); 1401 break; 1402 } 1403 case AlphaPixelChannel: 1404 { 1405 SetPixelAlpha(image,ClampToQuantum(QuantumRange*source_pixels[i]), 1406 q); 1407 break; 1408 } 1409 } 1410 i++; 1411 q+=GetPixelChannels(image); 1412 } 1413 if (SyncCacheViewAuthenticPixels(image_view,exception) == MagickFalse) 1414 break; 1415 } 1416 image_view=DestroyCacheView(image_view); 1417 source_info=RelinquishVirtualMemory(source_info); 1418 return(MagickTrue); 1419 } 1420 1421 static MagickBooleanType InverseFourierTransformChannel( 1422 const Image *magnitude_image,const Image *phase_image, 1423 const PixelChannel channel,const MagickBooleanType modulus, 1424 Image *fourier_image,ExceptionInfo *exception) 1425 { 1426 fftw_complex 1427 *inverse_pixels; 1428 1429 FourierInfo 1430 fourier_info; 1431 1432 MagickBooleanType 1433 status; 1434 1435 MemoryInfo 1436 *inverse_info; 1437 1438 fourier_info.width=magnitude_image->columns; 1439 fourier_info.height=magnitude_image->rows; 1440 if ((magnitude_image->columns != magnitude_image->rows) || 1441 ((magnitude_image->columns % 2) != 0) || 1442 ((magnitude_image->rows % 2) != 0)) 1443 { 1444 size_t extent=magnitude_image->columns < magnitude_image->rows ? 1445 magnitude_image->rows : magnitude_image->columns; 1446 fourier_info.width=(extent & 0x01) == 1 ? extent+1UL : extent; 1447 } 1448 fourier_info.height=fourier_info.width; 1449 fourier_info.center=(ssize_t) (fourier_info.width/2L)+1L; 1450 fourier_info.channel=channel; 1451 fourier_info.modulus=modulus; 1452 inverse_info=AcquireVirtualMemory((size_t) fourier_info.width, 1453 (fourier_info.height/2+1)*sizeof(*inverse_pixels)); 1454 if (inverse_info == (MemoryInfo *) NULL) 1455 { 1456 (void) ThrowMagickException(exception,GetMagickModule(), 1457 ResourceLimitError,"MemoryAllocationFailed","`%s'", 1458 magnitude_image->filename); 1459 return(MagickFalse); 1460 } 1461 inverse_pixels=(fftw_complex *) GetVirtualMemoryBlob(inverse_info); 1462 status=InverseFourier(&fourier_info,magnitude_image,phase_image, 1463 inverse_pixels,exception); 1464 if (status != MagickFalse) 1465 status=InverseFourierTransform(&fourier_info,inverse_pixels,fourier_image, 1466 exception); 1467 inverse_info=RelinquishVirtualMemory(inverse_info); 1468 return(status); 1469 } 1470 #endif 1471 1472 MagickExport Image *InverseFourierTransformImage(const Image *magnitude_image, 1473 const Image *phase_image,const MagickBooleanType modulus, 1474 ExceptionInfo *exception) 1475 { 1476 Image 1477 *fourier_image; 1478 1479 assert(magnitude_image != (Image *) NULL); 1480 assert(magnitude_image->signature == MagickCoreSignature); 1481 if (magnitude_image->debug != MagickFalse) 1482 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s", 1483 magnitude_image->filename); 1484 if (phase_image == (Image *) NULL) 1485 { 1486 (void) ThrowMagickException(exception,GetMagickModule(),ImageError, 1487 "ImageSequenceRequired","`%s'",magnitude_image->filename); 1488 return((Image *) NULL); 1489 } 1490 #if !defined(MAGICKCORE_FFTW_DELEGATE) 1491 fourier_image=(Image *) NULL; 1492 (void) modulus; 1493 (void) ThrowMagickException(exception,GetMagickModule(), 1494 MissingDelegateWarning,"DelegateLibrarySupportNotBuiltIn","`%s' (FFTW)", 1495 magnitude_image->filename); 1496 #else 1497 { 1498 fourier_image=CloneImage(magnitude_image,magnitude_image->columns, 1499 magnitude_image->rows,MagickTrue,exception); 1500 if (fourier_image != (Image *) NULL) 1501 { 1502 MagickBooleanType 1503 is_gray, 1504 status; 1505 1506 status=MagickTrue; 1507 is_gray=IsImageGray(magnitude_image); 1508 if (is_gray != MagickFalse) 1509 is_gray=IsImageGray(phase_image); 1510 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1511 #pragma omp parallel sections 1512 #endif 1513 { 1514 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1515 #pragma omp section 1516 #endif 1517 { 1518 MagickBooleanType 1519 thread_status; 1520 1521 if (is_gray != MagickFalse) 1522 thread_status=InverseFourierTransformChannel(magnitude_image, 1523 phase_image,GrayPixelChannel,modulus,fourier_image,exception); 1524 else 1525 thread_status=InverseFourierTransformChannel(magnitude_image, 1526 phase_image,RedPixelChannel,modulus,fourier_image,exception); 1527 if (thread_status == MagickFalse) 1528 status=thread_status; 1529 } 1530 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1531 #pragma omp section 1532 #endif 1533 { 1534 MagickBooleanType 1535 thread_status; 1536 1537 thread_status=MagickTrue; 1538 if (is_gray == MagickFalse) 1539 thread_status=InverseFourierTransformChannel(magnitude_image, 1540 phase_image,GreenPixelChannel,modulus,fourier_image,exception); 1541 if (thread_status == MagickFalse) 1542 status=thread_status; 1543 } 1544 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1545 #pragma omp section 1546 #endif 1547 { 1548 MagickBooleanType 1549 thread_status; 1550 1551 thread_status=MagickTrue; 1552 if (is_gray == MagickFalse) 1553 thread_status=InverseFourierTransformChannel(magnitude_image, 1554 phase_image,BluePixelChannel,modulus,fourier_image,exception); 1555 if (thread_status == MagickFalse) 1556 status=thread_status; 1557 } 1558 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1559 #pragma omp section 1560 #endif 1561 { 1562 MagickBooleanType 1563 thread_status; 1564 1565 thread_status=MagickTrue; 1566 if (magnitude_image->colorspace == CMYKColorspace) 1567 thread_status=InverseFourierTransformChannel(magnitude_image, 1568 phase_image,BlackPixelChannel,modulus,fourier_image,exception); 1569 if (thread_status == MagickFalse) 1570 status=thread_status; 1571 } 1572 #if defined(MAGICKCORE_OPENMP_SUPPORT) 1573 #pragma omp section 1574 #endif 1575 { 1576 MagickBooleanType 1577 thread_status; 1578 1579 thread_status=MagickTrue; 1580 if (magnitude_image->alpha_trait != UndefinedPixelTrait) 1581 thread_status=InverseFourierTransformChannel(magnitude_image, 1582 phase_image,AlphaPixelChannel,modulus,fourier_image,exception); 1583 if (thread_status == MagickFalse) 1584 status=thread_status; 1585 } 1586 } 1587 if (status == MagickFalse) 1588 fourier_image=DestroyImage(fourier_image); 1589 } 1590 fftw_cleanup(); 1591 } 1592 #endif 1593 return(fourier_image); 1594 } 1595
