NFFT  3.3.0
mri.c
1 /*
2  * Copyright (c) 2002, 2015 Jens Keiner, Stefan Kunis, Daniel Potts
3  *
4  * This program is free software; you can redistribute it and/or modify it under
5  * the terms of the GNU General Public License as published by the Free Software
6  * Foundation; either version 2 of the License, or (at your option) any later
7  * version.
8  *
9  * This program is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
11  * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
12  * details.
13  *
14  * You should have received a copy of the GNU General Public License along with
15  * this program; if not, write to the Free Software Foundation, Inc., 51
16  * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
17  */
18 
19 /* $Id$ */
20 
21 #include "config.h"
22 
23 #include <string.h>
24 #include <math.h>
25 #ifdef HAVE_COMPLEX_H
26 #include <complex.h>
27 #endif
28 #include "nfft3.h"
29 #include "infft.h"
30 
36 typedef struct window_funct_plan_ {
37  int d;
38  int m;
39  int n[1];
40  double sigma[1];
41  double *b;
42  double *spline_coeffs;
45 
49 static void window_funct_init(window_funct_plan* ths, int m, int n, double sigma) {
50  ths->d=1;
51  ths->m=m;
52  ths->n[0]=n;
53  ths->sigma[0]=sigma;
54  WINDOW_HELP_INIT
55 }
56 
57 /*
58  * mri_inh_2d1d
59  */
60 
61 void mri_inh_2d1d_trafo(mri_inh_2d1d_plan *that) {
62  int l,j;
63  double _Complex *f = (double _Complex*) nfft_malloc(that->M_total*sizeof(double _Complex));
64  double _Complex *f_hat = (double _Complex*) nfft_malloc(that->N_total*sizeof(double _Complex));
65 
67  window_funct_init(ths,that->plan.m,that->N3,that->sigma3);
68 
69  /* the pointers that->f and that->f_hat have been modified by the solver */
70  that->plan.f = that->f;
71  that->plan.f_hat = that->f_hat;
72 
73 
74  memset(f,0,that->M_total*sizeof(double _Complex));
75  for(j=0;j<that->N_total;j++)
76  {
77  f_hat[j]=that->f_hat[j];
78  }
79 
80  for(l=-ths->n[0]/2;l<=ths->n[0]/2;l++) {
81  for(j=0;j<that->N_total;j++)
82  that->f_hat[j]*=cexp(-2*KPI*_Complex_I*that->w[j]*((double)l))/PHI_HUT(ths->n[0], ths->n[0]*that->w[j],0);
83  nfft_trafo(&that->plan);
84  for(j=0;j<that->M_total;j++){
85  /* PHI has compact support */
86  if(fabs(that->t[j]-((double)l)/((double)ths->n[0]))<that->plan.m/((double)ths->n[0]))
87  {
88  double phi_val = PHI(ths->n[0],that->t[j]-((double)l)/((double)ths->n[0]),0);
89  f[j]+=that->f[j]*phi_val;
90 // the line below causes internal compiler error for gcc 4.7.1
91 // f[j]+=that->f[j]*PHI(ths->n[0],that->t[j]-((double)l)/((double)ths->n[0]),0);
92  }
93  }
94  for(j=0;j<that->N_total;j++)
95  that->f_hat[j]=f_hat[j];
96  }
97 
98  nfft_free(that->plan.f);
99  that->f=f;
100  that->plan.f = that->f;
101 
102  nfft_free(f_hat);
103 
104  WINDOW_HELP_FINALIZE
105  nfft_free(ths);
106 }
107 
108 void mri_inh_2d1d_adjoint(mri_inh_2d1d_plan *that) {
109  int l,j;
110  double _Complex *f = (double _Complex*) nfft_malloc(that->M_total*sizeof(double _Complex));
111  double _Complex *f_hat = (double _Complex*) nfft_malloc(that->N_total*sizeof(double _Complex));
112 
114  window_funct_init(ths,that->plan.m,that->N3,that->sigma3);
115 
116  memset(f_hat,0,that->N_total*sizeof(double _Complex));
117 
118  /* the pointers that->f and that->f_hat have been modified by the solver */
119  that->plan.f = that->f;
120  that->plan.f_hat = that->f_hat;
121 
122  for(j=0;j<that->M_total;j++)
123  {
124  f[j]=that->f[j];
125  }
126 
127 
128 
129  for(l=-ths->n[0]/2;l<=ths->n[0]/2;l++) {
130 
131  for(j=0;j<that->M_total;j++) {
132  /* PHI has compact support */
133  if(fabs(that->t[j]-((double)l)/((double)ths->n[0]))<that->plan.m/((double)ths->n[0]))
134  that->f[j]*=PHI(ths->n[0],that->t[j]-((double)l)/((double)ths->n[0]),0);
135  else
136  that->f[j]=0.0;
137  }
138  nfft_adjoint(&that->plan);
139  for(j=0;j<that->N_total;j++)
140  f_hat[j]+=that->f_hat[j]*cexp(2*KPI*_Complex_I*that->w[j]*((double)l));
141  for(j=0;j<that->M_total;j++)
142  that->f[j]=f[j];
143  }
144 
145  for(j=0;j<that->N_total;j++)
146  {
147  f_hat[j] /= PHI_HUT(ths->n[0],ths->n[0]*that->w[j],0);
148  }
149 
150  nfft_free(that->plan.f_hat);
151  that->f_hat=f_hat;
152  that->plan.f_hat = that->f_hat;
153 
154  nfft_free(f);
155 
156  WINDOW_HELP_FINALIZE
157  nfft_free(ths);
158 }
159 
160 void mri_inh_2d1d_init_guru(mri_inh_2d1d_plan *ths, int *N, int M, int *n,
161  int m, double sigma, unsigned nfft_flags, unsigned fftw_flags) {
162 
163  nfft_init_guru(&ths->plan,2,N,M,n,m,nfft_flags,fftw_flags);
164  ths->N3=N[2];
165  ths->sigma3=sigma;
166  ths->N_total = ths->plan.N_total;
167  ths->M_total = ths->plan.M_total;
168  ths->f = ths->plan.f;
169  ths->f_hat = ths->plan.f_hat;
170 
171  ths->t = (double*) nfft_malloc(ths->M_total*sizeof(double));
172  ths->w = (double*) nfft_malloc(ths->N_total*sizeof(double));
173 
174  ths->mv_trafo = (void (*) (void* ))mri_inh_2d1d_trafo;
175  ths->mv_adjoint = (void (*) (void* ))mri_inh_2d1d_adjoint;
176 }
177 
178 void mri_inh_2d1d_finalize(mri_inh_2d1d_plan *ths) {
179  nfft_free(ths->t);
180  nfft_free(ths->w);
181 
182  /* the pointers ths->f and ths->f_hat have been modified by the solver */
183  ths->plan.f = ths->f;
184  ths->plan.f_hat = ths->f_hat;
185 
186  nfft_finalize(&ths->plan);
187 }
188 
189 /*
190  * mri_inh_3d
191  */
192 
193 void mri_inh_3d_trafo(mri_inh_3d_plan *that) {
194  int l,j;
196  window_funct_init(ths,that->plan.m,that->N3,that->sigma3);
197 
198  /* the pointers that->f has been modified by the solver */
199  that->plan.f =that->f ;
200 
201 
202 
203  for(j=0;j<that->N_total;j++) {
204  for(l=-ths->n[0]/2;l<ths->n[0]/2;l++)
205  {
206  /* PHI has compact support */
207  if(fabs(that->w[j]-((double)l)/((double)ths->n[0]))<ths->m/((double)ths->n[0]))
208  that->plan.f_hat[j*ths->n[0]+(l+ths->n[0]/2)]= that->f_hat[j]*PHI(ths->n[0],that->w[j]-((double)l)/((double)ths->n[0]),0);
209  else
210  that->plan.f_hat[j*ths->n[0]+(l+ths->n[0]/2)]=0.0;
211  }
212  }
213 
214  nfft_trafo(&that->plan);
215 
216  for(j=0;j<that->M_total;j++)
217  {
218  that->f[j] /= PHI_HUT(ths->n[0],ths->n[0]*that->plan.x[3*j+2],0);
219  }
220 
221  WINDOW_HELP_FINALIZE
222  nfft_free(ths);
223 }
224 
225 void mri_inh_3d_adjoint(mri_inh_3d_plan *that) {
226  int l,j;
228  window_funct_init(ths,that->plan.m,that->N3,that->sigma3);
229 
230  /* the pointers that->f has been modified by the solver */
231  that->plan.f =that->f ;
232 
233  for(j=0;j<that->M_total;j++)
234  {
235  that->f[j] /= PHI_HUT(ths->n[0],ths->n[0]*that->plan.x[3*j+2],0);
236  }
237 
238  nfft_adjoint(&that->plan);
239 
240  for(j=0;j<that->N_total;j++) {
241  that->f_hat[j]=0.0;
242  for(l=-ths->n[0]/2;l<ths->n[0]/2;l++)
243  {
244  /* PHI has compact support */
245  if(fabs(that->w[j]-((double)l)/((double)ths->n[0]))<ths->m/((double)ths->n[0]))
246  that->f_hat[j]+= that->plan.f_hat[j*ths->n[0]+(l+ths->n[0]/2)]*PHI(ths->n[0],that->w[j]-((double)l)/((double)ths->n[0]),0);
247  }
248  }
249 
250 
251  WINDOW_HELP_FINALIZE
252  nfft_free(ths);
253 }
254 
255 void mri_inh_3d_init_guru(mri_inh_3d_plan *ths, int *N, int M, int *n,
256  int m, double sigma, unsigned nfft_flags, unsigned fftw_flags) {
257  ths->N3=N[2];
258  ths->sigma3=sigma;
259  nfft_init_guru(&ths->plan,3,N,M,n,m,nfft_flags,fftw_flags);
260  ths->N_total = N[0]*N[1];
261  ths->M_total = ths->plan.M_total;
262  ths->f = ths->plan.f;
263  ths->f_hat = (double _Complex*) nfft_malloc(ths->N_total*sizeof(double _Complex));
264  ths->w = (double*) nfft_malloc(ths->N_total*sizeof(double));
265 
266  ths->mv_trafo = (void (*) (void* ))mri_inh_3d_trafo;
267  ths->mv_adjoint = (void (*) (void* ))mri_inh_3d_adjoint;
268 }
269 
270 void mri_inh_3d_finalize(mri_inh_3d_plan *ths) {
271  nfft_free(ths->w);
272  nfft_free(ths->f_hat);
273  nfft_finalize(&ths->plan);
274 }
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:527
void(* mv_trafo)(void *)
Transform.
Definition: nfft3.h:527
void(* mv_trafo)(void *)
Transform.
Definition: nfft3.h:527
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:194
window_funct_plan is a plan to use the window functions independent of the nfft
Definition: mri.c:36
void(* mv_adjoint)(void *)
Adjoint transform.
Definition: nfft3.h:527
fftw_complex * f
Samples.
Definition: nfft3.h:527
fftw_complex * f_hat
Fourier coefficients.
Definition: nfft3.h:527
fftw_complex * f
Samples.
Definition: nfft3.h:194
fftw_complex * f
Samples.
Definition: nfft3.h:527
NFFT_INT N_total
Total number of Fourier coefficients.
Definition: nfft3.h:527
void nfft_free(void *p)
NFFT_INT N_total
Total number of Fourier coefficients.
Definition: nfft3.h:194
NFFT_INT M_total
Total number of samples.
Definition: nfft3.h:527
NFFT_INT M_total
Total number of samples.
Definition: nfft3.h:527
NFFT_INT M_total
Total number of samples.
Definition: nfft3.h:194
double * spline_coeffs
input for de Boor algorithm, if B_SPLINE or SINC_2m is defined
Definition: mri.c:42
void(* mv_adjoint)(void *)
Adjoint transform.
Definition: nfft3.h:527
void * nfft_malloc(size_t n)
double * x
Nodes in time/spatial domain, size is doubles.
Definition: nfft3.h:194
NFFT_INT N_total
Total number of Fourier coefficients.
Definition: nfft3.h:527
NFFT_INT m
Cut-off parameter for window function.
Definition: nfft3.h:194