00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "config.h"
00022
00023 #include <string.h>
00024 #include <math.h>
00025 #ifdef HAVE_COMPLEX_H
00026 #include <complex.h>
00027 #endif
00028 #include "nfft3util.h"
00029 #include "nfft3.h"
00030 #include "infft.h"
00031
00037 typedef struct window_funct_plan_ {
00038 int d;
00039 int m;
00040 int n[1];
00041 double sigma[1];
00042 double *b;
00043 double *spline_coeffs;
00045 } window_funct_plan;
00046
00050 static void window_funct_init(window_funct_plan* ths, int m, int n, double sigma) {
00051 ths->d=1;
00052 ths->m=m;
00053 ths->n[0]=n;
00054 ths->sigma[0]=sigma;
00055 WINDOW_HELP_INIT
00056 }
00057
00058
00059
00060
00061
00062 void mri_inh_2d1d_trafo(mri_inh_2d1d_plan *that) {
00063 int l,j;
00064 double _Complex *f = (double _Complex*) nfft_malloc(that->M_total*sizeof(double _Complex));
00065 double _Complex *f_hat = (double _Complex*) nfft_malloc(that->N_total*sizeof(double _Complex));
00066
00067 window_funct_plan *ths = (window_funct_plan*) nfft_malloc(sizeof(window_funct_plan));
00068 window_funct_init(ths,that->plan.m,that->N3,that->sigma3);
00069
00070
00071 that->plan.f = that->f;
00072 that->plan.f_hat = that->f_hat;
00073
00074
00075 memset(f,0,that->M_total*sizeof(double _Complex));
00076 for(j=0;j<that->N_total;j++)
00077 {
00078 f_hat[j]=that->f_hat[j];
00079 }
00080
00081 for(l=-ths->n[0]/2;l<=ths->n[0]/2;l++) {
00082 for(j=0;j<that->N_total;j++)
00083 that->f_hat[j]*=cexp(-2*PI*_Complex_I*that->w[j]*((double)l))/PHI_HUT(ths->n[0]*that->w[j],0);
00084 nfft_trafo(&that->plan);
00085 for(j=0;j<that->M_total;j++){
00086
00087 if(fabs(that->t[j]-((double)l)/((double)ths->n[0]))<that->plan.m/((double)ths->n[0]))
00088 f[j]+=that->f[j]*PHI(that->t[j]-((double)l)/((double)ths->n[0]),0);
00089 }
00090 for(j=0;j<that->N_total;j++)
00091 that->f_hat[j]=f_hat[j];
00092 }
00093
00094 nfft_free(that->plan.f);
00095 that->f=f;
00096 that->plan.f = that->f;
00097
00098 nfft_free(f_hat);
00099
00100 WINDOW_HELP_FINALIZE
00101 nfft_free(ths);
00102 }
00103
00104 void mri_inh_2d1d_adjoint(mri_inh_2d1d_plan *that) {
00105 int l,j;
00106 double _Complex *f = (double _Complex*) nfft_malloc(that->M_total*sizeof(double _Complex));
00107 double _Complex *f_hat = (double _Complex*) nfft_malloc(that->N_total*sizeof(double _Complex));
00108
00109 window_funct_plan *ths = (window_funct_plan*) nfft_malloc(sizeof(window_funct_plan));
00110 window_funct_init(ths,that->plan.m,that->N3,that->sigma3);
00111
00112 memset(f_hat,0,that->N_total*sizeof(double _Complex));
00113
00114
00115 that->plan.f = that->f;
00116 that->plan.f_hat = that->f_hat;
00117
00118 for(j=0;j<that->M_total;j++)
00119 {
00120 f[j]=that->f[j];
00121 }
00122
00123
00124
00125 for(l=-ths->n[0]/2;l<=ths->n[0]/2;l++) {
00126
00127 for(j=0;j<that->M_total;j++) {
00128
00129 if(fabs(that->t[j]-((double)l)/((double)ths->n[0]))<that->plan.m/((double)ths->n[0]))
00130 that->f[j]*=PHI(that->t[j]-((double)l)/((double)ths->n[0]),0);
00131 else
00132 that->f[j]=0.0;
00133 }
00134 nfft_adjoint(&that->plan);
00135 for(j=0;j<that->N_total;j++)
00136 f_hat[j]+=that->f_hat[j]*cexp(2*PI*_Complex_I*that->w[j]*((double)l));
00137 for(j=0;j<that->M_total;j++)
00138 that->f[j]=f[j];
00139 }
00140
00141 for(j=0;j<that->N_total;j++)
00142 {
00143 f_hat[j] /= PHI_HUT(ths->n[0]*that->w[j],0);
00144 }
00145
00146 nfft_free(that->plan.f_hat);
00147 that->f_hat=f_hat;
00148 that->plan.f_hat = that->f_hat;
00149
00150 nfft_free(f);
00151
00152 WINDOW_HELP_FINALIZE
00153 nfft_free(ths);
00154 }
00155
00156 void mri_inh_2d1d_init_guru(mri_inh_2d1d_plan *ths, int *N, int M, int *n,
00157 int m, double sigma, unsigned nfft_flags, unsigned fftw_flags) {
00158
00159 nfft_init_guru(&ths->plan,2,N,M,n,m,nfft_flags,fftw_flags);
00160 ths->N3=N[2];
00161 ths->sigma3=sigma;
00162 ths->N_total = ths->plan.N_total;
00163 ths->M_total = ths->plan.M_total;
00164 ths->f = ths->plan.f;
00165 ths->f_hat = ths->plan.f_hat;
00166
00167 ths->t = (double*) nfft_malloc(ths->M_total*sizeof(double));
00168 ths->w = (double*) nfft_malloc(ths->N_total*sizeof(double));
00169
00170 ths->mv_trafo = (void (*) (void* ))mri_inh_2d1d_trafo;
00171 ths->mv_adjoint = (void (*) (void* ))mri_inh_2d1d_adjoint;
00172 }
00173
00174 void mri_inh_2d1d_finalize(mri_inh_2d1d_plan *ths) {
00175 nfft_free(ths->t);
00176 nfft_free(ths->w);
00177
00178
00179 ths->plan.f = ths->f;
00180 ths->plan.f_hat = ths->f_hat;
00181
00182 nfft_finalize(&ths->plan);
00183 }
00184
00185
00186
00187
00188
00189 void mri_inh_3d_trafo(mri_inh_3d_plan *that) {
00190 int l,j;
00191 window_funct_plan *ths = (window_funct_plan*) nfft_malloc(sizeof(window_funct_plan));
00192 window_funct_init(ths,that->plan.m,that->N3,that->sigma3);
00193
00194
00195 that->plan.f =that->f ;
00196
00197
00198
00199 for(j=0;j<that->N_total;j++) {
00200 for(l=-ths->n[0]/2;l<ths->n[0]/2;l++)
00201 {
00202
00203 if(fabs(that->w[j]-((double)l)/((double)ths->n[0]))<ths->m/((double)ths->n[0]))
00204 that->plan.f_hat[j*ths->n[0]+(l+ths->n[0]/2)]= that->f_hat[j]*PHI(that->w[j]-((double)l)/((double)ths->n[0]),0);
00205 else
00206 that->plan.f_hat[j*ths->n[0]+(l+ths->n[0]/2)]=0.0;
00207 }
00208 }
00209
00210 nfft_trafo(&that->plan);
00211
00212 for(j=0;j<that->M_total;j++)
00213 {
00214 that->f[j] /= PHI_HUT(ths->n[0]*that->plan.x[3*j+2],0);
00215 }
00216
00217 WINDOW_HELP_FINALIZE
00218 nfft_free(ths);
00219 }
00220
00221 void mri_inh_3d_adjoint(mri_inh_3d_plan *that) {
00222 int l,j;
00223 window_funct_plan *ths = (window_funct_plan*) nfft_malloc(sizeof(window_funct_plan));
00224 window_funct_init(ths,that->plan.m,that->N3,that->sigma3);
00225
00226
00227 that->plan.f =that->f ;
00228
00229 for(j=0;j<that->M_total;j++)
00230 {
00231 that->f[j] /= PHI_HUT(ths->n[0]*that->plan.x[3*j+2],0);
00232 }
00233
00234 nfft_adjoint(&that->plan);
00235
00236 for(j=0;j<that->N_total;j++) {
00237 that->f_hat[j]=0.0;
00238 for(l=-ths->n[0]/2;l<ths->n[0]/2;l++)
00239 {
00240
00241 if(fabs(that->w[j]-((double)l)/((double)ths->n[0]))<ths->m/((double)ths->n[0]))
00242 that->f_hat[j]+= that->plan.f_hat[j*ths->n[0]+(l+ths->n[0]/2)]*PHI(that->w[j]-((double)l)/((double)ths->n[0]),0);
00243 }
00244 }
00245
00246
00247 WINDOW_HELP_FINALIZE
00248 nfft_free(ths);
00249 }
00250
00251 void mri_inh_3d_init_guru(mri_inh_3d_plan *ths, int *N, int M, int *n,
00252 int m, double sigma, unsigned nfft_flags, unsigned fftw_flags) {
00253 ths->N3=N[2];
00254 ths->sigma3=sigma;
00255 nfft_init_guru(&ths->plan,3,N,M,n,m,nfft_flags,fftw_flags);
00256 ths->N_total = N[0]*N[1];
00257 ths->M_total = ths->plan.M_total;
00258 ths->f = ths->plan.f;
00259 ths->f_hat = (double _Complex*) nfft_malloc(ths->N_total*sizeof(double _Complex));
00260 ths->w = (double*) nfft_malloc(ths->N_total*sizeof(double));
00261
00262 ths->mv_trafo = (void (*) (void* ))mri_inh_3d_trafo;
00263 ths->mv_adjoint = (void (*) (void* ))mri_inh_3d_adjoint;
00264 }
00265
00266 void mri_inh_3d_finalize(mri_inh_3d_plan *ths) {
00267 nfft_free(ths->w);
00268 nfft_free(ths->f_hat);
00269 nfft_finalize(&ths->plan);
00270 }