Changeset 814
- Timestamp:
- May 28, 2004, 3:11:01 PM (22 years ago)
- Location:
- trunk/psLib/src
- Files:
-
- 6 edited
-
dataManip/psFFT.c (modified) (25 diffs)
-
dataManip/psFFT.h (modified) (2 diffs)
-
dataManip/psVectorFFT.c (modified) (25 diffs)
-
dataManip/psVectorFFT.h (modified) (2 diffs)
-
fft/psVectorFFT.c (modified) (25 diffs)
-
fft/psVectorFFT.h (modified) (2 diffs)
Legend:
- Unmodified
- Added
- Removed
-
trunk/psLib/src/dataManip/psFFT.c
r809 r814 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $8 * @date $Date: 2004-05-2 8 21:07:17$7 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-05-29 01:11:01 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 25 25 static bool p_fftwWisdomImported = false; 26 26 27 psImage* psImageFFT( const psImage *image, int flags)27 psImage* psImageFFT(psImage* out, const psImage* in, int flags) 28 28 { 29 29 static fftwf_plan lastForwardPlan = NULL; … … 39 39 psElemType type; 40 40 fftwf_plan plan; 41 psImage* out = NULL;42 41 43 42 /* got good image data? */ 44 if (i mage==NULL) {43 if (in==NULL) { 45 44 return NULL; 46 45 } … … 51 50 } 52 51 53 type = i mage->type.type;52 type = in->type.type; 54 53 55 54 if ( (type != PS_TYPE_F32) && (type != PS_TYPE_C32) ) { … … 72 71 } 73 72 74 numRows = i mage->numRows;75 numCols = (image->numCols-1)*2;73 numRows = in->numRows; 74 numCols = in->numCols; 76 75 // n.b., numCols needs adjustment for c->r, see below. 77 76 … … 80 79 // Assuming image is from a r2c transform. 81 80 // (FFTW only uses nx/2+1 for r->c transform results) 82 numCols = ( image->numCols-1)*2;83 84 out = psImage Alloc(numCols,numRows,PS_TYPE_F32);81 numCols = (numCols-1)*2; 82 83 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 85 84 86 85 if (lastReversePlanRows != numRows || lastReversePlanCols != numCols … … 88 87 || lastReversePlan == NULL) { 89 88 lastReversePlan = fftwf_plan_dft_c2r_2d(numCols, numRows, 90 (fftwf_complex*)image->data.C32[0],out->data.F32[0], 89 (fftwf_complex*)in->data.C32[0], 90 out->data.F32[0], 91 91 P_FFTW_PLAN_RIGOR); 92 92 } 93 93 } else { // complex result desired 94 out = psImage Alloc(numCols,numRows,PS_TYPE_C32);94 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 95 95 if (lastReversePlanRows != numRows || lastReversePlanCols != numCols 96 96 || lastReversePlanResultType != PS_TYPE_C32 97 97 || lastReversePlan == NULL) { 98 98 lastReversePlan = fftwf_plan_dft_2d(numCols, numRows, 99 (fftwf_complex*)image->data.C32[0], 100 (fftwf_complex*)out->data.C32[0],FFTW_BACKWARD, 99 (fftwf_complex*)in->data.C32[0], 100 (fftwf_complex*)out->data.C32[0], 101 FFTW_BACKWARD, 101 102 P_FFTW_PLAN_RIGOR); 102 103 } … … 106 107 if (type == PS_TYPE_F32) { // real data input 107 108 // (FFTW only uses nx/2+1 for r->c transform results) 108 out = psImage Alloc(numCols/2+1,numRows,PS_TYPE_C32);109 out = psImageRecycle(out,numCols/2+1,numRows,PS_TYPE_C32); 109 110 110 111 if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols … … 112 113 /* need new plan created */ 113 114 lastForwardPlan = fftwf_plan_dft_r2c_2d(numCols, numRows, 114 image->data.F32[0],(fftwf_complex*)out->data.C32[0], 115 in->data.F32[0], 116 (fftwf_complex*)out->data.C32[0], 115 117 P_FFTW_PLAN_RIGOR); 116 118 } 117 119 } else if (type == PS_TYPE_C32) { // complex data input 118 out = psImage Alloc(numCols,numRows,PS_TYPE_C32);120 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 119 121 120 122 if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols 121 123 || lastForwardPlanType != type || lastForwardPlan == NULL) { 122 124 lastForwardPlan = fftwf_plan_dft_2d(numCols, numRows, 123 (fftwf_complex*)i mage->data.C32[0],125 (fftwf_complex*)in->data.C32[0], 124 126 (fftwf_complex*)out->data.C32[0], 125 FFTW_FORWARD,P_FFTW_PLAN_RIGOR); 127 FFTW_FORWARD, 128 P_FFTW_PLAN_RIGOR); 126 129 } 127 130 } … … 171 174 psC32* inRow; 172 175 173 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F32);176 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 174 177 for (unsigned int row=0;row<numRows;row++) { 175 178 outRow = out->data.F32[row]; … … 184 187 psC64* inRow; 185 188 186 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F64);189 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64); 187 190 for (unsigned int row=0;row<numRows;row++) { 188 191 outRow = out->data.F64[row]; … … 223 226 psLogMsg(__func__,PS_LOG_WARN,"Imaginary portion of a non-Complex type called for. " 224 227 "A zero image was returned."); 225 out = psImageRe alloc(out,numCols,numRows,type);228 out = psImageRecycle(out,numCols,numRows,type); 226 229 memset(out->data.v[0],0,PSELEMTYPE_SIZEOF(type)*numCols*numRows); 227 230 return out; … … 232 235 psC32* inRow; 233 236 234 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F32);237 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 235 238 for (unsigned int row=0;row<numRows;row++) { 236 239 outRow = out->data.F32[row]; … … 245 248 psC64* inRow; 246 249 247 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F64);250 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64); 248 251 for (unsigned int row=0;row<numRows;row++) { 249 252 outRow = out->data.F64[row]; … … 264 267 } 265 268 266 psImage *psImageComplex(psImage *real, const psImage *imag)269 psImage *psImageComplex(psImage* out, psImage *real, const psImage *imag) 267 270 { 268 271 psElemType type; 269 272 unsigned int numCols; 270 273 unsigned int numRows; 271 psImage* out = NULL;272 274 273 275 … … 282 284 if (imag->type.type != type) { 283 285 psError(__func__,"The inputs to psImageComplex must be the same type."); 286 return NULL; 287 } 288 289 if (imag->numCols != numCols || 290 imag->numRows != numRows) { 291 psError(__func__,"The inputs to psImageComplex must be the same dimensions."); 284 292 return NULL; 285 293 } … … 300 308 psF32* imagRow; 301 309 302 out = psImage Alloc(numCols,numRows,PS_TYPE_C32);310 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 303 311 304 312 for (unsigned int row=0;row<numRows;row++) { … … 311 319 } 312 320 } 313 } else if (type == PS_TYPE_ C64) {321 } else if (type == PS_TYPE_F64) { 314 322 psC64* outRow; 315 323 psF64* realRow; 316 324 psF64* imagRow; 317 325 318 out = psImage Alloc(numCols,numRows,PS_TYPE_C64);326 out = psImageRecycle(out, numCols,numRows,PS_TYPE_C64); 319 327 for (unsigned int row=0;row<numRows;row++) { 320 328 outRow = out->data.C64[row]; … … 363 371 psC32* inRow; 364 372 365 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_C32);373 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 366 374 for (unsigned int row=0;row<numRows;row++) { 367 375 outRow = out->data.C32[row]; … … 376 384 psC64* inRow; 377 385 378 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_C64);386 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C64); 379 387 for (unsigned int row=0;row<numRows;row++) { 380 388 outRow = out->data.C64[row]; … … 395 403 } 396 404 397 psImage *psImagePowerSpectrum( const psImage* in)405 psImage *psImagePowerSpectrum(psImage* out, const psImage* in) 398 406 { 399 407 psElemType type; … … 401 409 unsigned int numRows; 402 410 int numElementsSquared; 403 psImage* out = NULL;404 405 411 406 412 if (in == NULL) { … … 426 432 427 433 428 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F32);434 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 429 435 for (unsigned int row=0;row<numRows;row++) { 430 436 outRow = out->data.F32[row]; … … 444 450 445 451 446 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F64);452 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64); 447 453 for (unsigned int row=0;row<numRows;row++) { 448 454 outRow = out->data.F64[row]; … … 465 471 466 472 } 473 474 /************************************** Vector Functions ***************************************/ 475 476 psVector* psVectorFFT(psVector* out, const psVector* in, int flags) 477 { 478 static fftwf_plan lastForwardPlan = NULL; 479 static int lastForwardPlanType = 0; 480 static int lastForwardPlanElements = 0; 481 static fftwf_plan lastReversePlan = NULL; 482 static int lastReversePlanResultType = 0; 483 static int lastReversePlanElements = 0; 484 unsigned int numElements; 485 psElemType type; 486 fftwf_plan plan; 487 488 /* got good image data? */ 489 if (in==NULL) { 490 return NULL; 491 } 492 493 if ( ((flags & PS_FFT_REVERSE) == 0) && ((flags & PS_FFT_REAL_RESULT) != 0) ) { 494 psError(__func__,"Can not perform a forward FFT to real result."); 495 return NULL; 496 } 497 498 type = in->type.type; 499 500 if ( (type != PS_TYPE_F32) && (type != PS_TYPE_C32) ) { 501 psError(__func__,"Input image must be a 32-bit float or complex image (type=%d)", 502 type); 503 return NULL; 504 } 505 506 if ( (type != PS_TYPE_C32) && ((flags & PS_FFT_REVERSE) != 0) ) { 507 psError(__func__,"Input image must be complex image for reverse FFT (type=%d).", 508 type); 509 return NULL; 510 511 } 512 513 /* make sure the system-level wisdom information is imported. */ 514 if (! p_fftwWisdomImported) { 515 fftwf_import_system_wisdom(); 516 p_fftwWisdomImported = true; 517 } 518 519 numElements = in->n; 520 // n.b., numElements needs adjustment for c->r, see below. 521 522 if (flags & PS_FFT_REVERSE) { // reverse transform 523 if (flags & PS_FFT_REAL_RESULT) { // real result desired (i.e., c2r) 524 // Assuming image is from a r2c transform. 525 // (FFTW only uses nx/2+1 for r->c transform results) 526 numElements = (numElements-1)*2; 527 528 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 529 530 if (lastReversePlanElements != numElements || 531 lastReversePlanResultType != PS_TYPE_F32 || 532 lastReversePlan == NULL) { 533 lastReversePlan = fftwf_plan_dft_c2r_1d(numElements, 534 (fftwf_complex*)in->vec.cf, 535 out->vec.f, 536 P_FFTW_PLAN_RIGOR); 537 } 538 } else { // complex result desired 539 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 540 if (lastReversePlanElements != numElements || 541 lastReversePlanResultType != PS_TYPE_C32 || 542 lastReversePlan == NULL) { 543 lastReversePlan = fftwf_plan_dft_1d(numElements, 544 (fftwf_complex*)in->vec.cf, 545 (fftwf_complex*)out->vec.cf, 546 FFTW_BACKWARD, 547 P_FFTW_PLAN_RIGOR); 548 } 549 } 550 plan = lastReversePlan; 551 } else { // forward transform 552 if (type == PS_TYPE_F32) { // real data input 553 // (FFTW only uses nx/2+1 for r->c transform results) 554 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 555 556 if (lastForwardPlanElements != numElements || 557 lastForwardPlanType != type || 558 lastForwardPlan == NULL) { 559 /* need new plan created */ 560 lastForwardPlan = fftwf_plan_dft_r2c_1d(numElements, 561 in->vec.f, 562 (fftwf_complex*)out->vec.cf, 563 P_FFTW_PLAN_RIGOR); 564 } 565 } else if (type == PS_TYPE_C32) { // complex data input 566 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 567 568 if (lastForwardPlanElements != numElements || 569 lastForwardPlanType != type || 570 lastForwardPlan == NULL) { 571 lastForwardPlan = fftwf_plan_dft_1d(numElements, 572 (fftwf_complex*)in->vec.cf, 573 (fftwf_complex*)out->vec.cf, 574 FFTW_FORWARD, 575 P_FFTW_PLAN_RIGOR); 576 } 577 } 578 plan = lastForwardPlan; 579 } 580 581 /* check if a plan exists now*/ 582 if (plan == NULL) { 583 psError(__func__,"Failed to create FFTW plan."); 584 psVectorFree(out); 585 return NULL; 586 } 587 588 /* finally, call FFTW with the plan made above */ 589 fftwf_execute(plan); 590 591 return out; 592 593 } 594 595 596 psVector *psVectorReal(psVector *out, const psVector* in) 597 { 598 psElemType type; 599 unsigned int numElements; 600 601 if (in == NULL) { 602 return NULL; 603 } 604 605 type = in->type.type; 606 numElements = in->n; 607 608 /* if not a complex number, this is logically just a copy */ 609 if (! PS_IS_PSELEMTYPE_COMPLEX(type)) { 610 // Warn user, as this is probably not expected 611 psLogMsg(__func__,PS_LOG_WARN,"Real portion of a non-Complex type called called for. " 612 "Just a vector copy was performed."); 613 out = psVectorRecycle(out,type,numElements); 614 memcpy(out->vec.v,in->vec.v,numElements*PSELEMTYPE_SIZEOF(type)); 615 return out; 616 } 617 618 if (type == PS_TYPE_C32) { 619 psF32* outVec; 620 psC32* inVec = in->vec.cf; 621 622 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 623 outVec = out->vec.f; 624 625 for (unsigned int i=0;i<numElements;i++) { 626 outVec[i] = crealf(inVec[i]); 627 } 628 } else { 629 psError(__func__,"Can not extract real component from given vector type (%d).", 630 type); 631 psVectorFree(out); 632 return NULL; 633 } 634 635 return out; 636 } 637 638 psVector *psVectorImaginary(psVector *out, const psVector* in) 639 { 640 psElemType type; 641 unsigned int numElements; 642 643 644 if (in == NULL) { 645 return NULL; 646 } 647 648 type = in->type.type; 649 numElements = in->n; 650 651 /* if not a complex number, this is logically just zeroed image of same size */ 652 if (! PS_IS_PSELEMTYPE_COMPLEX(type)) { 653 // Warn user, as this is probably not expected 654 psLogMsg(__func__,PS_LOG_WARN,"Imaginary portion of a non-Complex type called for. " 655 "A zeroed vector was returned."); 656 out = psVectorRecycle(out,numElements,type); 657 memset(out->vec.v,0,PSELEMTYPE_SIZEOF(type)*numElements); 658 return out; 659 } 660 661 if (type == PS_TYPE_C32) { 662 psF32* outVec; 663 psC32* inVec = in->vec.cf; 664 665 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 666 outVec = out->vec.f; 667 668 for (unsigned int i=0;i<numElements;i++) { 669 outVec[i] = cimagf(inVec[i]); 670 } 671 } else { 672 psError(__func__,"Can not extract imaginary component from given vector type (%d).", 673 type); 674 psVectorFree(out); 675 return NULL; 676 } 677 678 return out; 679 } 680 681 psVector *psVectorComplex(psVector* out, psVector *real, const psVector *imag) 682 { 683 psElemType type; 684 unsigned int numElements; 685 686 687 if (real == NULL || imag == NULL) { 688 return NULL; 689 } 690 691 type = real->type.type; 692 if (real->n >= imag->n) { 693 numElements = real->n; 694 } else { 695 numElements = imag->n; 696 } 697 698 if (imag->type.type != type) { 699 psError(__func__,"The inputs to psVectorComplex must be the same type."); 700 return NULL; 701 } 702 703 if (PS_IS_PSELEMTYPE_COMPLEX(type)) { 704 psError(__func__,"The inputs to psVectorComplex can not be complex."); 705 return NULL; 706 } 707 708 if (type == PS_TYPE_F32) { 709 psC32* outVec; 710 psF32* realVec = real->vec.f; 711 psF32* imagVec = imag->vec.f; 712 713 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 714 outVec = out->vec.cf; 715 716 for (unsigned int i=0;i<numElements;i++) { 717 outVec[i] = realVec[i] + I*imagVec[i]; 718 } 719 } else { 720 psError(__func__,"Can not merge real and imaginary portions for given vector type (%d).", 721 type); 722 psVectorFree(out); 723 return NULL; 724 } 725 726 return out; 727 } 728 729 psVector *psVectorConjugate(psVector *out, const psVector *in) 730 { 731 psElemType type; 732 unsigned int numElements; 733 734 735 if (in == NULL) { 736 return NULL; 737 } 738 739 type = in->type.type; 740 numElements = in->n; 741 742 /* if not a complex number, this is logically just a image copy */ 743 if (! PS_IS_PSELEMTYPE_COMPLEX(type)) { 744 // Warn user, as this is probably not expected 745 psLogMsg(__func__,PS_LOG_WARN,"Complex Conjugate of a non-Complex type called for. " 746 "Vector copy was performed instead."); 747 748 out = psVectorRecycle(out,type,numElements); 749 memcpy(out->vec.v,in->vec.v,PSELEMTYPE_SIZEOF(type)*numElements); 750 return out; 751 } 752 753 if (type == PS_TYPE_C32) { 754 psC32* outVec; 755 psC32* inVec = in->vec.cf; 756 757 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 758 outVec = out->vec.cf; 759 760 for (unsigned int i=0;i<numElements;i++) { 761 outVec[i] = crealf(inVec[i]) - I*cimagf(inVec[i]); 762 } 763 } else { 764 psError(__func__,"Can not compute complex conjugate for given vector type (%d).", 765 type); 766 psVectorFree(out); 767 return NULL; 768 } 769 770 return out; 771 } 772 773 psVector *psVectorPowerSpectrum(psVector* out, const psVector* in) 774 { 775 psElemType type; 776 unsigned int numElements; 777 int numElementsSquared; 778 779 if (in == NULL) { 780 return NULL; 781 } 782 783 type = in->type.type; 784 numElements = in->n; 785 numElementsSquared = numElements*numElements; 786 787 /* if not a complex number, this is not implemented */ 788 if (! PS_IS_PSELEMTYPE_COMPLEX(type)) { 789 psError(__func__,"Power Spectrum for non-complex inputs is not implemented."); 790 return NULL; 791 } 792 793 if (type == PS_TYPE_C32) { 794 psF32* outVec; 795 psC32* inVec = in->vec.cf; 796 psF32 real; 797 psF32 imag; 798 799 800 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 801 outVec = out->vec.f; 802 803 for (unsigned int i=0;i<numElements;i++) { 804 real = crealf(inVec[i]); 805 imag = cimagf(inVec[i]); 806 outVec[i] = (real*real+imag*imag)/numElementsSquared; 807 } 808 } else { 809 psError(__func__,"Can not power spectrum for given vector type (%d).", 810 type); 811 psVectorFree(out); 812 return NULL; 813 } 814 815 return out; 816 817 } -
trunk/psLib/src/dataManip/psFFT.h
r801 r814 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 1$ $Name: not supported by cvs2svn $8 * @date $Date: 2004-05-2 8 03:15:40$7 * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-05-29 01:11:01 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 23 23 24 24 typedef enum { 25 PS_FFT_FORWARD = 0, ///< psImageFFT/psVectorFFT should perform a forward FFT. 26 PS_FFT_REVERSE = 1, ///< psImageFFT/psVectorFFT should perform a reverse FFT. 25 /// psImageFFT/psVectorFFT should perform a forward FFT. 26 PS_FFT_FORWARD = 0, 27 28 ///< psImageFFT/psVectorFFT should perform a reverse FFT. 29 PS_FFT_REVERSE = 1, 30 31 /** psImageFFT/psVectorFFT should return a real image. This is valid for 32 * only reverse FFT, i.e., the psImageFFT/psVectorFFT flag parameter 33 * should appear as PS_FFT_REVERSE+PS_FFT_REAL_RESULT. 34 */ 27 35 PS_FFT_REAL_RESULT = 2 28 /**< psImageFFT/psVectorFFT should return a real image. This is valid for29 * only reverse FFT, i.e., the psImageFFT/psVectorFFT flag parameter30 * should appear as PS_FFT_REVERSE+PS_FFT_REAL_RESULT.31 */32 36 } psFftFlags; 33 37 34 psImage *psImageFFT(const psImage *image, int flags); 35 psImage *psImageReal(psImage *out, const psImage* in); 36 psImage *psImageImaginary(psImage *out, const psImage* in); 37 psImage *psImageComplex(psImage *real, const psImage *imag); 38 psImage *psImageConjugate(psImage *out, const psImage *in); 39 psImage *psImagePowerSpectrum(const psImage* in); 38 psImage* psImageFFT(psImage* out, const psImage* in, int flags); 39 psImage* psImageReal(psImage *out, const psImage* in); 40 psImage* psImageImaginary(psImage *out, const psImage* in); 41 psImage* psImageComplex(psImage* out, psImage *real, const psImage *imag); 42 psImage* psImageConjugate(psImage *out, const psImage *in); 43 psImage* psImagePowerSpectrum(psImage* out, const psImage* in); 44 45 psVector* psVectorFFT(psVector* out, const psVector* in, int flags); 46 psVector* psVectorReal(psVector* out, const psVector* in); 47 psVector* psVectorImaginary(psVector* out, const psVector* in); 48 psVector* psVectorComplex(psVector* out, psVector* real, const psVector* imag); 49 psVector* psVectorConjugate(psVector* out, const psVector* in); 50 psVector* psVectorPowerSpectrum(psVector* out, const psVector* in); 40 51 41 52 -
trunk/psLib/src/dataManip/psVectorFFT.c
r809 r814 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $8 * @date $Date: 2004-05-2 8 21:07:17$7 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-05-29 01:11:01 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 25 25 static bool p_fftwWisdomImported = false; 26 26 27 psImage* psImageFFT( const psImage *image, int flags)27 psImage* psImageFFT(psImage* out, const psImage* in, int flags) 28 28 { 29 29 static fftwf_plan lastForwardPlan = NULL; … … 39 39 psElemType type; 40 40 fftwf_plan plan; 41 psImage* out = NULL;42 41 43 42 /* got good image data? */ 44 if (i mage==NULL) {43 if (in==NULL) { 45 44 return NULL; 46 45 } … … 51 50 } 52 51 53 type = i mage->type.type;52 type = in->type.type; 54 53 55 54 if ( (type != PS_TYPE_F32) && (type != PS_TYPE_C32) ) { … … 72 71 } 73 72 74 numRows = i mage->numRows;75 numCols = (image->numCols-1)*2;73 numRows = in->numRows; 74 numCols = in->numCols; 76 75 // n.b., numCols needs adjustment for c->r, see below. 77 76 … … 80 79 // Assuming image is from a r2c transform. 81 80 // (FFTW only uses nx/2+1 for r->c transform results) 82 numCols = ( image->numCols-1)*2;83 84 out = psImage Alloc(numCols,numRows,PS_TYPE_F32);81 numCols = (numCols-1)*2; 82 83 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 85 84 86 85 if (lastReversePlanRows != numRows || lastReversePlanCols != numCols … … 88 87 || lastReversePlan == NULL) { 89 88 lastReversePlan = fftwf_plan_dft_c2r_2d(numCols, numRows, 90 (fftwf_complex*)image->data.C32[0],out->data.F32[0], 89 (fftwf_complex*)in->data.C32[0], 90 out->data.F32[0], 91 91 P_FFTW_PLAN_RIGOR); 92 92 } 93 93 } else { // complex result desired 94 out = psImage Alloc(numCols,numRows,PS_TYPE_C32);94 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 95 95 if (lastReversePlanRows != numRows || lastReversePlanCols != numCols 96 96 || lastReversePlanResultType != PS_TYPE_C32 97 97 || lastReversePlan == NULL) { 98 98 lastReversePlan = fftwf_plan_dft_2d(numCols, numRows, 99 (fftwf_complex*)image->data.C32[0], 100 (fftwf_complex*)out->data.C32[0],FFTW_BACKWARD, 99 (fftwf_complex*)in->data.C32[0], 100 (fftwf_complex*)out->data.C32[0], 101 FFTW_BACKWARD, 101 102 P_FFTW_PLAN_RIGOR); 102 103 } … … 106 107 if (type == PS_TYPE_F32) { // real data input 107 108 // (FFTW only uses nx/2+1 for r->c transform results) 108 out = psImage Alloc(numCols/2+1,numRows,PS_TYPE_C32);109 out = psImageRecycle(out,numCols/2+1,numRows,PS_TYPE_C32); 109 110 110 111 if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols … … 112 113 /* need new plan created */ 113 114 lastForwardPlan = fftwf_plan_dft_r2c_2d(numCols, numRows, 114 image->data.F32[0],(fftwf_complex*)out->data.C32[0], 115 in->data.F32[0], 116 (fftwf_complex*)out->data.C32[0], 115 117 P_FFTW_PLAN_RIGOR); 116 118 } 117 119 } else if (type == PS_TYPE_C32) { // complex data input 118 out = psImage Alloc(numCols,numRows,PS_TYPE_C32);120 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 119 121 120 122 if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols 121 123 || lastForwardPlanType != type || lastForwardPlan == NULL) { 122 124 lastForwardPlan = fftwf_plan_dft_2d(numCols, numRows, 123 (fftwf_complex*)i mage->data.C32[0],125 (fftwf_complex*)in->data.C32[0], 124 126 (fftwf_complex*)out->data.C32[0], 125 FFTW_FORWARD,P_FFTW_PLAN_RIGOR); 127 FFTW_FORWARD, 128 P_FFTW_PLAN_RIGOR); 126 129 } 127 130 } … … 171 174 psC32* inRow; 172 175 173 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F32);176 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 174 177 for (unsigned int row=0;row<numRows;row++) { 175 178 outRow = out->data.F32[row]; … … 184 187 psC64* inRow; 185 188 186 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F64);189 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64); 187 190 for (unsigned int row=0;row<numRows;row++) { 188 191 outRow = out->data.F64[row]; … … 223 226 psLogMsg(__func__,PS_LOG_WARN,"Imaginary portion of a non-Complex type called for. " 224 227 "A zero image was returned."); 225 out = psImageRe alloc(out,numCols,numRows,type);228 out = psImageRecycle(out,numCols,numRows,type); 226 229 memset(out->data.v[0],0,PSELEMTYPE_SIZEOF(type)*numCols*numRows); 227 230 return out; … … 232 235 psC32* inRow; 233 236 234 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F32);237 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 235 238 for (unsigned int row=0;row<numRows;row++) { 236 239 outRow = out->data.F32[row]; … … 245 248 psC64* inRow; 246 249 247 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F64);250 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64); 248 251 for (unsigned int row=0;row<numRows;row++) { 249 252 outRow = out->data.F64[row]; … … 264 267 } 265 268 266 psImage *psImageComplex(psImage *real, const psImage *imag)269 psImage *psImageComplex(psImage* out, psImage *real, const psImage *imag) 267 270 { 268 271 psElemType type; 269 272 unsigned int numCols; 270 273 unsigned int numRows; 271 psImage* out = NULL;272 274 273 275 … … 282 284 if (imag->type.type != type) { 283 285 psError(__func__,"The inputs to psImageComplex must be the same type."); 286 return NULL; 287 } 288 289 if (imag->numCols != numCols || 290 imag->numRows != numRows) { 291 psError(__func__,"The inputs to psImageComplex must be the same dimensions."); 284 292 return NULL; 285 293 } … … 300 308 psF32* imagRow; 301 309 302 out = psImage Alloc(numCols,numRows,PS_TYPE_C32);310 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 303 311 304 312 for (unsigned int row=0;row<numRows;row++) { … … 311 319 } 312 320 } 313 } else if (type == PS_TYPE_ C64) {321 } else if (type == PS_TYPE_F64) { 314 322 psC64* outRow; 315 323 psF64* realRow; 316 324 psF64* imagRow; 317 325 318 out = psImage Alloc(numCols,numRows,PS_TYPE_C64);326 out = psImageRecycle(out, numCols,numRows,PS_TYPE_C64); 319 327 for (unsigned int row=0;row<numRows;row++) { 320 328 outRow = out->data.C64[row]; … … 363 371 psC32* inRow; 364 372 365 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_C32);373 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 366 374 for (unsigned int row=0;row<numRows;row++) { 367 375 outRow = out->data.C32[row]; … … 376 384 psC64* inRow; 377 385 378 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_C64);386 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C64); 379 387 for (unsigned int row=0;row<numRows;row++) { 380 388 outRow = out->data.C64[row]; … … 395 403 } 396 404 397 psImage *psImagePowerSpectrum( const psImage* in)405 psImage *psImagePowerSpectrum(psImage* out, const psImage* in) 398 406 { 399 407 psElemType type; … … 401 409 unsigned int numRows; 402 410 int numElementsSquared; 403 psImage* out = NULL;404 405 411 406 412 if (in == NULL) { … … 426 432 427 433 428 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F32);434 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 429 435 for (unsigned int row=0;row<numRows;row++) { 430 436 outRow = out->data.F32[row]; … … 444 450 445 451 446 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F64);452 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64); 447 453 for (unsigned int row=0;row<numRows;row++) { 448 454 outRow = out->data.F64[row]; … … 465 471 466 472 } 473 474 /************************************** Vector Functions ***************************************/ 475 476 psVector* psVectorFFT(psVector* out, const psVector* in, int flags) 477 { 478 static fftwf_plan lastForwardPlan = NULL; 479 static int lastForwardPlanType = 0; 480 static int lastForwardPlanElements = 0; 481 static fftwf_plan lastReversePlan = NULL; 482 static int lastReversePlanResultType = 0; 483 static int lastReversePlanElements = 0; 484 unsigned int numElements; 485 psElemType type; 486 fftwf_plan plan; 487 488 /* got good image data? */ 489 if (in==NULL) { 490 return NULL; 491 } 492 493 if ( ((flags & PS_FFT_REVERSE) == 0) && ((flags & PS_FFT_REAL_RESULT) != 0) ) { 494 psError(__func__,"Can not perform a forward FFT to real result."); 495 return NULL; 496 } 497 498 type = in->type.type; 499 500 if ( (type != PS_TYPE_F32) && (type != PS_TYPE_C32) ) { 501 psError(__func__,"Input image must be a 32-bit float or complex image (type=%d)", 502 type); 503 return NULL; 504 } 505 506 if ( (type != PS_TYPE_C32) && ((flags & PS_FFT_REVERSE) != 0) ) { 507 psError(__func__,"Input image must be complex image for reverse FFT (type=%d).", 508 type); 509 return NULL; 510 511 } 512 513 /* make sure the system-level wisdom information is imported. */ 514 if (! p_fftwWisdomImported) { 515 fftwf_import_system_wisdom(); 516 p_fftwWisdomImported = true; 517 } 518 519 numElements = in->n; 520 // n.b., numElements needs adjustment for c->r, see below. 521 522 if (flags & PS_FFT_REVERSE) { // reverse transform 523 if (flags & PS_FFT_REAL_RESULT) { // real result desired (i.e., c2r) 524 // Assuming image is from a r2c transform. 525 // (FFTW only uses nx/2+1 for r->c transform results) 526 numElements = (numElements-1)*2; 527 528 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 529 530 if (lastReversePlanElements != numElements || 531 lastReversePlanResultType != PS_TYPE_F32 || 532 lastReversePlan == NULL) { 533 lastReversePlan = fftwf_plan_dft_c2r_1d(numElements, 534 (fftwf_complex*)in->vec.cf, 535 out->vec.f, 536 P_FFTW_PLAN_RIGOR); 537 } 538 } else { // complex result desired 539 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 540 if (lastReversePlanElements != numElements || 541 lastReversePlanResultType != PS_TYPE_C32 || 542 lastReversePlan == NULL) { 543 lastReversePlan = fftwf_plan_dft_1d(numElements, 544 (fftwf_complex*)in->vec.cf, 545 (fftwf_complex*)out->vec.cf, 546 FFTW_BACKWARD, 547 P_FFTW_PLAN_RIGOR); 548 } 549 } 550 plan = lastReversePlan; 551 } else { // forward transform 552 if (type == PS_TYPE_F32) { // real data input 553 // (FFTW only uses nx/2+1 for r->c transform results) 554 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 555 556 if (lastForwardPlanElements != numElements || 557 lastForwardPlanType != type || 558 lastForwardPlan == NULL) { 559 /* need new plan created */ 560 lastForwardPlan = fftwf_plan_dft_r2c_1d(numElements, 561 in->vec.f, 562 (fftwf_complex*)out->vec.cf, 563 P_FFTW_PLAN_RIGOR); 564 } 565 } else if (type == PS_TYPE_C32) { // complex data input 566 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 567 568 if (lastForwardPlanElements != numElements || 569 lastForwardPlanType != type || 570 lastForwardPlan == NULL) { 571 lastForwardPlan = fftwf_plan_dft_1d(numElements, 572 (fftwf_complex*)in->vec.cf, 573 (fftwf_complex*)out->vec.cf, 574 FFTW_FORWARD, 575 P_FFTW_PLAN_RIGOR); 576 } 577 } 578 plan = lastForwardPlan; 579 } 580 581 /* check if a plan exists now*/ 582 if (plan == NULL) { 583 psError(__func__,"Failed to create FFTW plan."); 584 psVectorFree(out); 585 return NULL; 586 } 587 588 /* finally, call FFTW with the plan made above */ 589 fftwf_execute(plan); 590 591 return out; 592 593 } 594 595 596 psVector *psVectorReal(psVector *out, const psVector* in) 597 { 598 psElemType type; 599 unsigned int numElements; 600 601 if (in == NULL) { 602 return NULL; 603 } 604 605 type = in->type.type; 606 numElements = in->n; 607 608 /* if not a complex number, this is logically just a copy */ 609 if (! PS_IS_PSELEMTYPE_COMPLEX(type)) { 610 // Warn user, as this is probably not expected 611 psLogMsg(__func__,PS_LOG_WARN,"Real portion of a non-Complex type called called for. " 612 "Just a vector copy was performed."); 613 out = psVectorRecycle(out,type,numElements); 614 memcpy(out->vec.v,in->vec.v,numElements*PSELEMTYPE_SIZEOF(type)); 615 return out; 616 } 617 618 if (type == PS_TYPE_C32) { 619 psF32* outVec; 620 psC32* inVec = in->vec.cf; 621 622 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 623 outVec = out->vec.f; 624 625 for (unsigned int i=0;i<numElements;i++) { 626 outVec[i] = crealf(inVec[i]); 627 } 628 } else { 629 psError(__func__,"Can not extract real component from given vector type (%d).", 630 type); 631 psVectorFree(out); 632 return NULL; 633 } 634 635 return out; 636 } 637 638 psVector *psVectorImaginary(psVector *out, const psVector* in) 639 { 640 psElemType type; 641 unsigned int numElements; 642 643 644 if (in == NULL) { 645 return NULL; 646 } 647 648 type = in->type.type; 649 numElements = in->n; 650 651 /* if not a complex number, this is logically just zeroed image of same size */ 652 if (! PS_IS_PSELEMTYPE_COMPLEX(type)) { 653 // Warn user, as this is probably not expected 654 psLogMsg(__func__,PS_LOG_WARN,"Imaginary portion of a non-Complex type called for. " 655 "A zeroed vector was returned."); 656 out = psVectorRecycle(out,numElements,type); 657 memset(out->vec.v,0,PSELEMTYPE_SIZEOF(type)*numElements); 658 return out; 659 } 660 661 if (type == PS_TYPE_C32) { 662 psF32* outVec; 663 psC32* inVec = in->vec.cf; 664 665 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 666 outVec = out->vec.f; 667 668 for (unsigned int i=0;i<numElements;i++) { 669 outVec[i] = cimagf(inVec[i]); 670 } 671 } else { 672 psError(__func__,"Can not extract imaginary component from given vector type (%d).", 673 type); 674 psVectorFree(out); 675 return NULL; 676 } 677 678 return out; 679 } 680 681 psVector *psVectorComplex(psVector* out, psVector *real, const psVector *imag) 682 { 683 psElemType type; 684 unsigned int numElements; 685 686 687 if (real == NULL || imag == NULL) { 688 return NULL; 689 } 690 691 type = real->type.type; 692 if (real->n >= imag->n) { 693 numElements = real->n; 694 } else { 695 numElements = imag->n; 696 } 697 698 if (imag->type.type != type) { 699 psError(__func__,"The inputs to psVectorComplex must be the same type."); 700 return NULL; 701 } 702 703 if (PS_IS_PSELEMTYPE_COMPLEX(type)) { 704 psError(__func__,"The inputs to psVectorComplex can not be complex."); 705 return NULL; 706 } 707 708 if (type == PS_TYPE_F32) { 709 psC32* outVec; 710 psF32* realVec = real->vec.f; 711 psF32* imagVec = imag->vec.f; 712 713 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 714 outVec = out->vec.cf; 715 716 for (unsigned int i=0;i<numElements;i++) { 717 outVec[i] = realVec[i] + I*imagVec[i]; 718 } 719 } else { 720 psError(__func__,"Can not merge real and imaginary portions for given vector type (%d).", 721 type); 722 psVectorFree(out); 723 return NULL; 724 } 725 726 return out; 727 } 728 729 psVector *psVectorConjugate(psVector *out, const psVector *in) 730 { 731 psElemType type; 732 unsigned int numElements; 733 734 735 if (in == NULL) { 736 return NULL; 737 } 738 739 type = in->type.type; 740 numElements = in->n; 741 742 /* if not a complex number, this is logically just a image copy */ 743 if (! PS_IS_PSELEMTYPE_COMPLEX(type)) { 744 // Warn user, as this is probably not expected 745 psLogMsg(__func__,PS_LOG_WARN,"Complex Conjugate of a non-Complex type called for. " 746 "Vector copy was performed instead."); 747 748 out = psVectorRecycle(out,type,numElements); 749 memcpy(out->vec.v,in->vec.v,PSELEMTYPE_SIZEOF(type)*numElements); 750 return out; 751 } 752 753 if (type == PS_TYPE_C32) { 754 psC32* outVec; 755 psC32* inVec = in->vec.cf; 756 757 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 758 outVec = out->vec.cf; 759 760 for (unsigned int i=0;i<numElements;i++) { 761 outVec[i] = crealf(inVec[i]) - I*cimagf(inVec[i]); 762 } 763 } else { 764 psError(__func__,"Can not compute complex conjugate for given vector type (%d).", 765 type); 766 psVectorFree(out); 767 return NULL; 768 } 769 770 return out; 771 } 772 773 psVector *psVectorPowerSpectrum(psVector* out, const psVector* in) 774 { 775 psElemType type; 776 unsigned int numElements; 777 int numElementsSquared; 778 779 if (in == NULL) { 780 return NULL; 781 } 782 783 type = in->type.type; 784 numElements = in->n; 785 numElementsSquared = numElements*numElements; 786 787 /* if not a complex number, this is not implemented */ 788 if (! PS_IS_PSELEMTYPE_COMPLEX(type)) { 789 psError(__func__,"Power Spectrum for non-complex inputs is not implemented."); 790 return NULL; 791 } 792 793 if (type == PS_TYPE_C32) { 794 psF32* outVec; 795 psC32* inVec = in->vec.cf; 796 psF32 real; 797 psF32 imag; 798 799 800 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 801 outVec = out->vec.f; 802 803 for (unsigned int i=0;i<numElements;i++) { 804 real = crealf(inVec[i]); 805 imag = cimagf(inVec[i]); 806 outVec[i] = (real*real+imag*imag)/numElementsSquared; 807 } 808 } else { 809 psError(__func__,"Can not power spectrum for given vector type (%d).", 810 type); 811 psVectorFree(out); 812 return NULL; 813 } 814 815 return out; 816 817 } -
trunk/psLib/src/dataManip/psVectorFFT.h
r801 r814 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 1$ $Name: not supported by cvs2svn $8 * @date $Date: 2004-05-2 8 03:15:40$7 * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-05-29 01:11:01 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 23 23 24 24 typedef enum { 25 PS_FFT_FORWARD = 0, ///< psImageFFT/psVectorFFT should perform a forward FFT. 26 PS_FFT_REVERSE = 1, ///< psImageFFT/psVectorFFT should perform a reverse FFT. 25 /// psImageFFT/psVectorFFT should perform a forward FFT. 26 PS_FFT_FORWARD = 0, 27 28 ///< psImageFFT/psVectorFFT should perform a reverse FFT. 29 PS_FFT_REVERSE = 1, 30 31 /** psImageFFT/psVectorFFT should return a real image. This is valid for 32 * only reverse FFT, i.e., the psImageFFT/psVectorFFT flag parameter 33 * should appear as PS_FFT_REVERSE+PS_FFT_REAL_RESULT. 34 */ 27 35 PS_FFT_REAL_RESULT = 2 28 /**< psImageFFT/psVectorFFT should return a real image. This is valid for29 * only reverse FFT, i.e., the psImageFFT/psVectorFFT flag parameter30 * should appear as PS_FFT_REVERSE+PS_FFT_REAL_RESULT.31 */32 36 } psFftFlags; 33 37 34 psImage *psImageFFT(const psImage *image, int flags); 35 psImage *psImageReal(psImage *out, const psImage* in); 36 psImage *psImageImaginary(psImage *out, const psImage* in); 37 psImage *psImageComplex(psImage *real, const psImage *imag); 38 psImage *psImageConjugate(psImage *out, const psImage *in); 39 psImage *psImagePowerSpectrum(const psImage* in); 38 psImage* psImageFFT(psImage* out, const psImage* in, int flags); 39 psImage* psImageReal(psImage *out, const psImage* in); 40 psImage* psImageImaginary(psImage *out, const psImage* in); 41 psImage* psImageComplex(psImage* out, psImage *real, const psImage *imag); 42 psImage* psImageConjugate(psImage *out, const psImage *in); 43 psImage* psImagePowerSpectrum(psImage* out, const psImage* in); 44 45 psVector* psVectorFFT(psVector* out, const psVector* in, int flags); 46 psVector* psVectorReal(psVector* out, const psVector* in); 47 psVector* psVectorImaginary(psVector* out, const psVector* in); 48 psVector* psVectorComplex(psVector* out, psVector* real, const psVector* imag); 49 psVector* psVectorConjugate(psVector* out, const psVector* in); 50 psVector* psVectorPowerSpectrum(psVector* out, const psVector* in); 40 51 41 52 -
trunk/psLib/src/fft/psVectorFFT.c
r809 r814 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 3$ $Name: not supported by cvs2svn $8 * @date $Date: 2004-05-2 8 21:07:17$7 * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-05-29 01:11:01 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 25 25 static bool p_fftwWisdomImported = false; 26 26 27 psImage* psImageFFT( const psImage *image, int flags)27 psImage* psImageFFT(psImage* out, const psImage* in, int flags) 28 28 { 29 29 static fftwf_plan lastForwardPlan = NULL; … … 39 39 psElemType type; 40 40 fftwf_plan plan; 41 psImage* out = NULL;42 41 43 42 /* got good image data? */ 44 if (i mage==NULL) {43 if (in==NULL) { 45 44 return NULL; 46 45 } … … 51 50 } 52 51 53 type = i mage->type.type;52 type = in->type.type; 54 53 55 54 if ( (type != PS_TYPE_F32) && (type != PS_TYPE_C32) ) { … … 72 71 } 73 72 74 numRows = i mage->numRows;75 numCols = (image->numCols-1)*2;73 numRows = in->numRows; 74 numCols = in->numCols; 76 75 // n.b., numCols needs adjustment for c->r, see below. 77 76 … … 80 79 // Assuming image is from a r2c transform. 81 80 // (FFTW only uses nx/2+1 for r->c transform results) 82 numCols = ( image->numCols-1)*2;83 84 out = psImage Alloc(numCols,numRows,PS_TYPE_F32);81 numCols = (numCols-1)*2; 82 83 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 85 84 86 85 if (lastReversePlanRows != numRows || lastReversePlanCols != numCols … … 88 87 || lastReversePlan == NULL) { 89 88 lastReversePlan = fftwf_plan_dft_c2r_2d(numCols, numRows, 90 (fftwf_complex*)image->data.C32[0],out->data.F32[0], 89 (fftwf_complex*)in->data.C32[0], 90 out->data.F32[0], 91 91 P_FFTW_PLAN_RIGOR); 92 92 } 93 93 } else { // complex result desired 94 out = psImage Alloc(numCols,numRows,PS_TYPE_C32);94 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 95 95 if (lastReversePlanRows != numRows || lastReversePlanCols != numCols 96 96 || lastReversePlanResultType != PS_TYPE_C32 97 97 || lastReversePlan == NULL) { 98 98 lastReversePlan = fftwf_plan_dft_2d(numCols, numRows, 99 (fftwf_complex*)image->data.C32[0], 100 (fftwf_complex*)out->data.C32[0],FFTW_BACKWARD, 99 (fftwf_complex*)in->data.C32[0], 100 (fftwf_complex*)out->data.C32[0], 101 FFTW_BACKWARD, 101 102 P_FFTW_PLAN_RIGOR); 102 103 } … … 106 107 if (type == PS_TYPE_F32) { // real data input 107 108 // (FFTW only uses nx/2+1 for r->c transform results) 108 out = psImage Alloc(numCols/2+1,numRows,PS_TYPE_C32);109 out = psImageRecycle(out,numCols/2+1,numRows,PS_TYPE_C32); 109 110 110 111 if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols … … 112 113 /* need new plan created */ 113 114 lastForwardPlan = fftwf_plan_dft_r2c_2d(numCols, numRows, 114 image->data.F32[0],(fftwf_complex*)out->data.C32[0], 115 in->data.F32[0], 116 (fftwf_complex*)out->data.C32[0], 115 117 P_FFTW_PLAN_RIGOR); 116 118 } 117 119 } else if (type == PS_TYPE_C32) { // complex data input 118 out = psImage Alloc(numCols,numRows,PS_TYPE_C32);120 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 119 121 120 122 if (lastForwardPlanRows != numRows || lastForwardPlanCols != numCols 121 123 || lastForwardPlanType != type || lastForwardPlan == NULL) { 122 124 lastForwardPlan = fftwf_plan_dft_2d(numCols, numRows, 123 (fftwf_complex*)i mage->data.C32[0],125 (fftwf_complex*)in->data.C32[0], 124 126 (fftwf_complex*)out->data.C32[0], 125 FFTW_FORWARD,P_FFTW_PLAN_RIGOR); 127 FFTW_FORWARD, 128 P_FFTW_PLAN_RIGOR); 126 129 } 127 130 } … … 171 174 psC32* inRow; 172 175 173 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F32);176 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 174 177 for (unsigned int row=0;row<numRows;row++) { 175 178 outRow = out->data.F32[row]; … … 184 187 psC64* inRow; 185 188 186 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F64);189 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64); 187 190 for (unsigned int row=0;row<numRows;row++) { 188 191 outRow = out->data.F64[row]; … … 223 226 psLogMsg(__func__,PS_LOG_WARN,"Imaginary portion of a non-Complex type called for. " 224 227 "A zero image was returned."); 225 out = psImageRe alloc(out,numCols,numRows,type);228 out = psImageRecycle(out,numCols,numRows,type); 226 229 memset(out->data.v[0],0,PSELEMTYPE_SIZEOF(type)*numCols*numRows); 227 230 return out; … … 232 235 psC32* inRow; 233 236 234 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F32);237 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 235 238 for (unsigned int row=0;row<numRows;row++) { 236 239 outRow = out->data.F32[row]; … … 245 248 psC64* inRow; 246 249 247 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F64);250 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64); 248 251 for (unsigned int row=0;row<numRows;row++) { 249 252 outRow = out->data.F64[row]; … … 264 267 } 265 268 266 psImage *psImageComplex(psImage *real, const psImage *imag)269 psImage *psImageComplex(psImage* out, psImage *real, const psImage *imag) 267 270 { 268 271 psElemType type; 269 272 unsigned int numCols; 270 273 unsigned int numRows; 271 psImage* out = NULL;272 274 273 275 … … 282 284 if (imag->type.type != type) { 283 285 psError(__func__,"The inputs to psImageComplex must be the same type."); 286 return NULL; 287 } 288 289 if (imag->numCols != numCols || 290 imag->numRows != numRows) { 291 psError(__func__,"The inputs to psImageComplex must be the same dimensions."); 284 292 return NULL; 285 293 } … … 300 308 psF32* imagRow; 301 309 302 out = psImage Alloc(numCols,numRows,PS_TYPE_C32);310 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 303 311 304 312 for (unsigned int row=0;row<numRows;row++) { … … 311 319 } 312 320 } 313 } else if (type == PS_TYPE_ C64) {321 } else if (type == PS_TYPE_F64) { 314 322 psC64* outRow; 315 323 psF64* realRow; 316 324 psF64* imagRow; 317 325 318 out = psImage Alloc(numCols,numRows,PS_TYPE_C64);326 out = psImageRecycle(out, numCols,numRows,PS_TYPE_C64); 319 327 for (unsigned int row=0;row<numRows;row++) { 320 328 outRow = out->data.C64[row]; … … 363 371 psC32* inRow; 364 372 365 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_C32);373 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C32); 366 374 for (unsigned int row=0;row<numRows;row++) { 367 375 outRow = out->data.C32[row]; … … 376 384 psC64* inRow; 377 385 378 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_C64);386 out = psImageRecycle(out,numCols,numRows,PS_TYPE_C64); 379 387 for (unsigned int row=0;row<numRows;row++) { 380 388 outRow = out->data.C64[row]; … … 395 403 } 396 404 397 psImage *psImagePowerSpectrum( const psImage* in)405 psImage *psImagePowerSpectrum(psImage* out, const psImage* in) 398 406 { 399 407 psElemType type; … … 401 409 unsigned int numRows; 402 410 int numElementsSquared; 403 psImage* out = NULL;404 405 411 406 412 if (in == NULL) { … … 426 432 427 433 428 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F32);434 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F32); 429 435 for (unsigned int row=0;row<numRows;row++) { 430 436 outRow = out->data.F32[row]; … … 444 450 445 451 446 out = psImageRe alloc(out,numCols,numRows,PS_TYPE_F64);452 out = psImageRecycle(out,numCols,numRows,PS_TYPE_F64); 447 453 for (unsigned int row=0;row<numRows;row++) { 448 454 outRow = out->data.F64[row]; … … 465 471 466 472 } 473 474 /************************************** Vector Functions ***************************************/ 475 476 psVector* psVectorFFT(psVector* out, const psVector* in, int flags) 477 { 478 static fftwf_plan lastForwardPlan = NULL; 479 static int lastForwardPlanType = 0; 480 static int lastForwardPlanElements = 0; 481 static fftwf_plan lastReversePlan = NULL; 482 static int lastReversePlanResultType = 0; 483 static int lastReversePlanElements = 0; 484 unsigned int numElements; 485 psElemType type; 486 fftwf_plan plan; 487 488 /* got good image data? */ 489 if (in==NULL) { 490 return NULL; 491 } 492 493 if ( ((flags & PS_FFT_REVERSE) == 0) && ((flags & PS_FFT_REAL_RESULT) != 0) ) { 494 psError(__func__,"Can not perform a forward FFT to real result."); 495 return NULL; 496 } 497 498 type = in->type.type; 499 500 if ( (type != PS_TYPE_F32) && (type != PS_TYPE_C32) ) { 501 psError(__func__,"Input image must be a 32-bit float or complex image (type=%d)", 502 type); 503 return NULL; 504 } 505 506 if ( (type != PS_TYPE_C32) && ((flags & PS_FFT_REVERSE) != 0) ) { 507 psError(__func__,"Input image must be complex image for reverse FFT (type=%d).", 508 type); 509 return NULL; 510 511 } 512 513 /* make sure the system-level wisdom information is imported. */ 514 if (! p_fftwWisdomImported) { 515 fftwf_import_system_wisdom(); 516 p_fftwWisdomImported = true; 517 } 518 519 numElements = in->n; 520 // n.b., numElements needs adjustment for c->r, see below. 521 522 if (flags & PS_FFT_REVERSE) { // reverse transform 523 if (flags & PS_FFT_REAL_RESULT) { // real result desired (i.e., c2r) 524 // Assuming image is from a r2c transform. 525 // (FFTW only uses nx/2+1 for r->c transform results) 526 numElements = (numElements-1)*2; 527 528 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 529 530 if (lastReversePlanElements != numElements || 531 lastReversePlanResultType != PS_TYPE_F32 || 532 lastReversePlan == NULL) { 533 lastReversePlan = fftwf_plan_dft_c2r_1d(numElements, 534 (fftwf_complex*)in->vec.cf, 535 out->vec.f, 536 P_FFTW_PLAN_RIGOR); 537 } 538 } else { // complex result desired 539 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 540 if (lastReversePlanElements != numElements || 541 lastReversePlanResultType != PS_TYPE_C32 || 542 lastReversePlan == NULL) { 543 lastReversePlan = fftwf_plan_dft_1d(numElements, 544 (fftwf_complex*)in->vec.cf, 545 (fftwf_complex*)out->vec.cf, 546 FFTW_BACKWARD, 547 P_FFTW_PLAN_RIGOR); 548 } 549 } 550 plan = lastReversePlan; 551 } else { // forward transform 552 if (type == PS_TYPE_F32) { // real data input 553 // (FFTW only uses nx/2+1 for r->c transform results) 554 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 555 556 if (lastForwardPlanElements != numElements || 557 lastForwardPlanType != type || 558 lastForwardPlan == NULL) { 559 /* need new plan created */ 560 lastForwardPlan = fftwf_plan_dft_r2c_1d(numElements, 561 in->vec.f, 562 (fftwf_complex*)out->vec.cf, 563 P_FFTW_PLAN_RIGOR); 564 } 565 } else if (type == PS_TYPE_C32) { // complex data input 566 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 567 568 if (lastForwardPlanElements != numElements || 569 lastForwardPlanType != type || 570 lastForwardPlan == NULL) { 571 lastForwardPlan = fftwf_plan_dft_1d(numElements, 572 (fftwf_complex*)in->vec.cf, 573 (fftwf_complex*)out->vec.cf, 574 FFTW_FORWARD, 575 P_FFTW_PLAN_RIGOR); 576 } 577 } 578 plan = lastForwardPlan; 579 } 580 581 /* check if a plan exists now*/ 582 if (plan == NULL) { 583 psError(__func__,"Failed to create FFTW plan."); 584 psVectorFree(out); 585 return NULL; 586 } 587 588 /* finally, call FFTW with the plan made above */ 589 fftwf_execute(plan); 590 591 return out; 592 593 } 594 595 596 psVector *psVectorReal(psVector *out, const psVector* in) 597 { 598 psElemType type; 599 unsigned int numElements; 600 601 if (in == NULL) { 602 return NULL; 603 } 604 605 type = in->type.type; 606 numElements = in->n; 607 608 /* if not a complex number, this is logically just a copy */ 609 if (! PS_IS_PSELEMTYPE_COMPLEX(type)) { 610 // Warn user, as this is probably not expected 611 psLogMsg(__func__,PS_LOG_WARN,"Real portion of a non-Complex type called called for. " 612 "Just a vector copy was performed."); 613 out = psVectorRecycle(out,type,numElements); 614 memcpy(out->vec.v,in->vec.v,numElements*PSELEMTYPE_SIZEOF(type)); 615 return out; 616 } 617 618 if (type == PS_TYPE_C32) { 619 psF32* outVec; 620 psC32* inVec = in->vec.cf; 621 622 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 623 outVec = out->vec.f; 624 625 for (unsigned int i=0;i<numElements;i++) { 626 outVec[i] = crealf(inVec[i]); 627 } 628 } else { 629 psError(__func__,"Can not extract real component from given vector type (%d).", 630 type); 631 psVectorFree(out); 632 return NULL; 633 } 634 635 return out; 636 } 637 638 psVector *psVectorImaginary(psVector *out, const psVector* in) 639 { 640 psElemType type; 641 unsigned int numElements; 642 643 644 if (in == NULL) { 645 return NULL; 646 } 647 648 type = in->type.type; 649 numElements = in->n; 650 651 /* if not a complex number, this is logically just zeroed image of same size */ 652 if (! PS_IS_PSELEMTYPE_COMPLEX(type)) { 653 // Warn user, as this is probably not expected 654 psLogMsg(__func__,PS_LOG_WARN,"Imaginary portion of a non-Complex type called for. " 655 "A zeroed vector was returned."); 656 out = psVectorRecycle(out,numElements,type); 657 memset(out->vec.v,0,PSELEMTYPE_SIZEOF(type)*numElements); 658 return out; 659 } 660 661 if (type == PS_TYPE_C32) { 662 psF32* outVec; 663 psC32* inVec = in->vec.cf; 664 665 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 666 outVec = out->vec.f; 667 668 for (unsigned int i=0;i<numElements;i++) { 669 outVec[i] = cimagf(inVec[i]); 670 } 671 } else { 672 psError(__func__,"Can not extract imaginary component from given vector type (%d).", 673 type); 674 psVectorFree(out); 675 return NULL; 676 } 677 678 return out; 679 } 680 681 psVector *psVectorComplex(psVector* out, psVector *real, const psVector *imag) 682 { 683 psElemType type; 684 unsigned int numElements; 685 686 687 if (real == NULL || imag == NULL) { 688 return NULL; 689 } 690 691 type = real->type.type; 692 if (real->n >= imag->n) { 693 numElements = real->n; 694 } else { 695 numElements = imag->n; 696 } 697 698 if (imag->type.type != type) { 699 psError(__func__,"The inputs to psVectorComplex must be the same type."); 700 return NULL; 701 } 702 703 if (PS_IS_PSELEMTYPE_COMPLEX(type)) { 704 psError(__func__,"The inputs to psVectorComplex can not be complex."); 705 return NULL; 706 } 707 708 if (type == PS_TYPE_F32) { 709 psC32* outVec; 710 psF32* realVec = real->vec.f; 711 psF32* imagVec = imag->vec.f; 712 713 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 714 outVec = out->vec.cf; 715 716 for (unsigned int i=0;i<numElements;i++) { 717 outVec[i] = realVec[i] + I*imagVec[i]; 718 } 719 } else { 720 psError(__func__,"Can not merge real and imaginary portions for given vector type (%d).", 721 type); 722 psVectorFree(out); 723 return NULL; 724 } 725 726 return out; 727 } 728 729 psVector *psVectorConjugate(psVector *out, const psVector *in) 730 { 731 psElemType type; 732 unsigned int numElements; 733 734 735 if (in == NULL) { 736 return NULL; 737 } 738 739 type = in->type.type; 740 numElements = in->n; 741 742 /* if not a complex number, this is logically just a image copy */ 743 if (! PS_IS_PSELEMTYPE_COMPLEX(type)) { 744 // Warn user, as this is probably not expected 745 psLogMsg(__func__,PS_LOG_WARN,"Complex Conjugate of a non-Complex type called for. " 746 "Vector copy was performed instead."); 747 748 out = psVectorRecycle(out,type,numElements); 749 memcpy(out->vec.v,in->vec.v,PSELEMTYPE_SIZEOF(type)*numElements); 750 return out; 751 } 752 753 if (type == PS_TYPE_C32) { 754 psC32* outVec; 755 psC32* inVec = in->vec.cf; 756 757 out = psVectorRecycle(out,PS_TYPE_C32,numElements); 758 outVec = out->vec.cf; 759 760 for (unsigned int i=0;i<numElements;i++) { 761 outVec[i] = crealf(inVec[i]) - I*cimagf(inVec[i]); 762 } 763 } else { 764 psError(__func__,"Can not compute complex conjugate for given vector type (%d).", 765 type); 766 psVectorFree(out); 767 return NULL; 768 } 769 770 return out; 771 } 772 773 psVector *psVectorPowerSpectrum(psVector* out, const psVector* in) 774 { 775 psElemType type; 776 unsigned int numElements; 777 int numElementsSquared; 778 779 if (in == NULL) { 780 return NULL; 781 } 782 783 type = in->type.type; 784 numElements = in->n; 785 numElementsSquared = numElements*numElements; 786 787 /* if not a complex number, this is not implemented */ 788 if (! PS_IS_PSELEMTYPE_COMPLEX(type)) { 789 psError(__func__,"Power Spectrum for non-complex inputs is not implemented."); 790 return NULL; 791 } 792 793 if (type == PS_TYPE_C32) { 794 psF32* outVec; 795 psC32* inVec = in->vec.cf; 796 psF32 real; 797 psF32 imag; 798 799 800 out = psVectorRecycle(out,PS_TYPE_F32,numElements); 801 outVec = out->vec.f; 802 803 for (unsigned int i=0;i<numElements;i++) { 804 real = crealf(inVec[i]); 805 imag = cimagf(inVec[i]); 806 outVec[i] = (real*real+imag*imag)/numElementsSquared; 807 } 808 } else { 809 psError(__func__,"Can not power spectrum for given vector type (%d).", 810 type); 811 psVectorFree(out); 812 return NULL; 813 } 814 815 return out; 816 817 } -
trunk/psLib/src/fft/psVectorFFT.h
r801 r814 5 5 * @author Robert DeSonia, MHPCC 6 6 * 7 * @version $Revision: 1. 1$ $Name: not supported by cvs2svn $8 * @date $Date: 2004-05-2 8 03:15:40$7 * @version $Revision: 1.2 $ $Name: not supported by cvs2svn $ 8 * @date $Date: 2004-05-29 01:11:01 $ 9 9 * 10 10 * Copyright 2004 Maui High Performance Computing Center, University of Hawaii … … 23 23 24 24 typedef enum { 25 PS_FFT_FORWARD = 0, ///< psImageFFT/psVectorFFT should perform a forward FFT. 26 PS_FFT_REVERSE = 1, ///< psImageFFT/psVectorFFT should perform a reverse FFT. 25 /// psImageFFT/psVectorFFT should perform a forward FFT. 26 PS_FFT_FORWARD = 0, 27 28 ///< psImageFFT/psVectorFFT should perform a reverse FFT. 29 PS_FFT_REVERSE = 1, 30 31 /** psImageFFT/psVectorFFT should return a real image. This is valid for 32 * only reverse FFT, i.e., the psImageFFT/psVectorFFT flag parameter 33 * should appear as PS_FFT_REVERSE+PS_FFT_REAL_RESULT. 34 */ 27 35 PS_FFT_REAL_RESULT = 2 28 /**< psImageFFT/psVectorFFT should return a real image. This is valid for29 * only reverse FFT, i.e., the psImageFFT/psVectorFFT flag parameter30 * should appear as PS_FFT_REVERSE+PS_FFT_REAL_RESULT.31 */32 36 } psFftFlags; 33 37 34 psImage *psImageFFT(const psImage *image, int flags); 35 psImage *psImageReal(psImage *out, const psImage* in); 36 psImage *psImageImaginary(psImage *out, const psImage* in); 37 psImage *psImageComplex(psImage *real, const psImage *imag); 38 psImage *psImageConjugate(psImage *out, const psImage *in); 39 psImage *psImagePowerSpectrum(const psImage* in); 38 psImage* psImageFFT(psImage* out, const psImage* in, int flags); 39 psImage* psImageReal(psImage *out, const psImage* in); 40 psImage* psImageImaginary(psImage *out, const psImage* in); 41 psImage* psImageComplex(psImage* out, psImage *real, const psImage *imag); 42 psImage* psImageConjugate(psImage *out, const psImage *in); 43 psImage* psImagePowerSpectrum(psImage* out, const psImage* in); 44 45 psVector* psVectorFFT(psVector* out, const psVector* in, int flags); 46 psVector* psVectorReal(psVector* out, const psVector* in); 47 psVector* psVectorImaginary(psVector* out, const psVector* in); 48 psVector* psVectorComplex(psVector* out, psVector* real, const psVector* imag); 49 psVector* psVectorConjugate(psVector* out, const psVector* in); 50 psVector* psVectorPowerSpectrum(psVector* out, const psVector* in); 40 51 41 52
Note:
See TracChangeset
for help on using the changeset viewer.
