fix fftw failure when using MKL based conda stack

lsst · Oct 2, 2021 · 59802e7 · 59802e7
1 parent 93a5c75
commit 59802e7
Showing 1 changed file with 8 additions and 6 deletions.
diff --git a/src/SincCoeffs.cc b/src/SincCoeffs.cc
@@ -423,12 +423,14 @@ std::shared_ptr<afw::image::Image<PixelT>> calcImageKSpaceReal(double const rad1
     int xcen = wid / 2, ycen = wid / 2;
     FftShifter fftshift(wid);
 
-    boost::shared_array<double> cimg(new double[wid * wid]);
-    double* c = cimg.get();
+    boost::shared_array<std::complex<double>> cimg(new std::complex<double>[wid * wid]);
+    std::complex<double>* c = cimg.get();
     // fftplan args: nx, ny, *in, *out, kindx, kindy, flags
     // - done in-situ if *in == *out
-    fftw_plan plan = fftw_plan_r2r_2d(wid, wid, c, c, FFTW_R2HC, FFTW_R2HC, FFTW_ESTIMATE);
 
+    //fftw_plan plan = fftw_plan_r2r_2d(wid, wid, c, c, FFTW_R2HC, FFTW_R2HC, FFTW_ESTIMATE);
+    fftw_plan plan = fftw_plan_dft_2d(wid, wid, reinterpret_cast<fftw_complex*>(c),
+                                      reinterpret_cast<fftw_complex*>(c), FFTW_BACKWARD, FFTW_ESTIMATE);
     // compute the k-space values and put them in the cimg array
     double const twoPiRad1 = geom::TWOPI * rad1;
     double const twoPiRad2 = geom::TWOPI * rad2;
@@ -447,11 +449,11 @@ std::shared_ptr<afw::image::Image<PixelT>> calcImageKSpaceReal(double const rad1
             double const airy1 = (rad1 > 0 ? rad1 * J1(twoPiRad1 * k) : 0.0) / k;
             double const airy2 = rad2 * J1(twoPiRad2 * k) / k;
             double const airy = airy2 - airy1;
-            c[fY * wid + fX] = airy;
+            c[fY * wid + fX] = {airy, 0.};
         }
     }
     int fxy = fftshift.shift(wid / 2);
-    c[fxy * wid + fxy] = geom::PI * (rad2 * rad2 - rad1 * rad1);
+    c[fxy * wid + fxy] = {geom::PI * (rad2 * rad2 - rad1 * rad1), 0.};
 
     // perform the fft and clean up after ourselves
     fftw_execute(plan);
@@ -468,7 +470,7 @@ std::shared_ptr<afw::image::Image<PixelT>> calcImageKSpaceReal(double const rad1
             // now need to reflect the quadrant we solved to the other three
             int fX = iX < hwid ? hwid - iX - 1 : iX - hwid + 1;
             int fY = iY < hwid ? hwid - iY - 1 : iY - hwid + 1;
-            *ptr = static_cast<PixelT>(c[fY * wid + fX] / (wid * wid));
+            *ptr = static_cast<PixelT>(c[fY * wid + fX].real() / (wid * wid));
             iX++;
         }
     }