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