Inverse FFT for 2D histogram question

Hello all,

I must first give a shoutout to the folks directly involved with the development of ROOT. Thank you!

So here is the issue that I am presenting.

I have created a CT simulator in Geant4 that will use root to calculate the FFT and iFFT of the projection data to ultimately get the CT slice image and attenuation data. I have successfully installed and tested the FFTW and I am pretty comfortable with the 1D array discrete FFT and iFFT. So all installation issues can be ruled out.

Now to the meat of the problem. The way that I am calculating my FFT and iFFT follows this procedure (as very similar to the tutorial FFT.C)

1D Histogrammed data with 128 bins - This Works
//Define the histogram (I am attaching a file known as CT1D2.C which is the whole thing)
Then implement the following snippets:

   TH1 *hm =0;
   TVirtualFFT::SetTransform(0);
   hm = Test1D->FFT(hm, "MAG");
   TH1 *hp = 0;
   hp = Test1D->FFT(hp, "PH");
   hp->SetTitle("Phase of the 1st transform");
   int n=128;
   int holly=128;
   Double_t *re_full = new Double_t[n];
   Double_t *im_full = new Double_t[n];
   TVirtualFFT *fft = TVirtualFFT::GetCurrentTransform();
   fft->GetPointsComplex(re_full,im_full);
   
   //Now let's make a backward transform:
   TVirtualFFT *fft_back = TVirtualFFT::FFT(1, &holly, "C2R M K");
   fft_back->SetPointsComplex(re_full,im_full);
   fft_back->Transform();
   TH1 *hb = 0;
   //Let's look at the output
   hb = TH1::TransformHisto(fft_back,hb,"Re");
   hb->SetTitle("The backward transform result");
   hb->Draw();

The result is that we have a few different histograms where the Test1D is the original histogram, the hm is the magnitude FFT, the hp is the phase FFT, and hb is the iFFT. This works great!

2D Histogram FFT and iFFT - not working

Now we move on to the trouble with 2D (I will attach a code known as CTcanvas2.C) This should also be similar for other 2D histograms as well. Basically I just want to transform the image with FFT, and then get the resulting 2D iFFT in ROOT. (There will be additional processing in the middle, including interpolation and so forth, but I’m just trying to return the original 2D image for now).

I’ve spent about two weeks crashing ROOT trying to get the iFFT to work, but to no avail. I’m just not sure how to command the FFT function to do the iFFT for 2D.

Here is the culprit/code.

   //fourier stuff
      
   TH1 *hm =0;
   TVirtualFFT::SetTransform(0);
   hm = CT histogram Test->FFT(hm, "MAG");
   TH1 *hp = 0;
   hp = CT histogram Test->FFT(hp, "PH");
   hp->SetTitle("Phase of the 1st transform");
   int n=360;
   int resx=360;
   int resy=360;
   Double_t *re_full = new Double_t[resx][resy];
   Double_t *im_full = new Double_t[resx][resy];
   TVirtualFFT *fft = TVirtualFFT::GetCurrentTransform();
   fft->GetPointsComplex(re_full,im_full);
   
   //Now let's make a backward transform:
   TVirtualFFT *fft_back = TVirtualFFT::FFT(2, &n, "C2R M K");
   fft_back->SetPointsComplex(re_full,im_full);
   fft_back->Transform();
   //TH1 *hb = 0;
   //Let's look at the output
   //hb = TH1::TransformHisto(fft_back,hb,"Re");
   //hb->SetTitle("The backward transform result");
   //hb->Draw();

Thank you so much for any help, I am eagerly waiting to get this thing working!

John Perry

CTcanvas2.C (173 KB)
CT1D2.C (5.22 KB)

Hi everyone,

Another update (I believe I have solved the issue and am posting for posterity’s sake).

The trick is for the FFT(ndims, array of bins, "flags"),
ndims=2,
array of bins must be defined this way Int_t n[2]={128,360}

therefore FFT(2,n,"C2R M K"), I am posting the solution so far.

If anyone disagrees please let me know.

Also, moving on to something that is just as tricky, does anyone have a good interpolation program for placing rays of projection data onto a 2D histogram? Each ray will have an angle and so on.

Thank you again,
John Perry

CTcanvas2.C (173 KB)