![]() |
NFFT
3.3.1
|
00001 /* 00002 * Copyright (c) 2002, 2016 Jens Keiner, Stefan Kunis, Daniel Potts 00003 * 00004 * This program is free software; you can redistribute it and/or modify it under 00005 * the terms of the GNU General Public License as published by the Free Software 00006 * Foundation; either version 2 of the License, or (at your option) any later 00007 * version. 00008 * 00009 * This program is distributed in the hope that it will be useful, but WITHOUT 00010 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 00011 * FOR A PARTICULAR PURPOSE. See the GNU General Public License for more 00012 * details. 00013 * 00014 * You should have received a copy of the GNU General Public License along with 00015 * this program; if not, write to the Free Software Foundation, Inc., 51 00016 * Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. 00017 */ 00018 00019 #include "config.h" 00020 00021 #include <string.h> 00022 #include <math.h> 00023 #ifdef HAVE_COMPLEX_H 00024 #include <complex.h> 00025 #endif 00026 #include "nfft3.h" 00027 #include "infft.h" 00028 00034 typedef struct window_funct_plan_ { 00035 int d; 00036 int m; 00037 int n[1]; 00038 double sigma[1]; 00039 double *b; 00040 } window_funct_plan; 00041 00045 static void window_funct_init(window_funct_plan* ths, int m, int n, double sigma) { 00046 ths->d=1; 00047 ths->m=m; 00048 ths->n[0]=n; 00049 ths->sigma[0]=sigma; 00050 WINDOW_HELP_INIT 00051 } 00052 00053 /* 00054 * mri_inh_2d1d 00055 */ 00056 00057 void mri_inh_2d1d_trafo(mri_inh_2d1d_plan *that) { 00058 int l,j; 00059 double _Complex *f = (double _Complex*) nfft_malloc(that->M_total*sizeof(double _Complex)); 00060 double _Complex *f_hat = (double _Complex*) nfft_malloc(that->N_total*sizeof(double _Complex)); 00061 00062 window_funct_plan *ths = (window_funct_plan*) nfft_malloc(sizeof(window_funct_plan)); 00063 window_funct_init(ths,that->plan.m,that->N3,that->sigma3); 00064 00065 /* the pointers that->f and that->f_hat have been modified by the solver */ 00066 that->plan.f = that->f; 00067 that->plan.f_hat = that->f_hat; 00068 00069 00070 memset(f,0,that->M_total*sizeof(double _Complex)); 00071 for(j=0;j<that->N_total;j++) 00072 { 00073 f_hat[j]=that->f_hat[j]; 00074 } 00075 00076 for(l=-ths->n[0]/2;l<=ths->n[0]/2;l++) { 00077 for(j=0;j<that->N_total;j++) 00078 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); 00079 nfft_trafo(&that->plan); 00080 for(j=0;j<that->M_total;j++){ 00081 /* PHI has compact support */ 00082 if(fabs(that->t[j]-((double)l)/((double)ths->n[0]))<that->plan.m/((double)ths->n[0])) 00083 { 00084 double phi_val = PHI(ths->n[0],that->t[j]-((double)l)/((double)ths->n[0]),0); 00085 f[j]+=that->f[j]*phi_val; 00086 // the line below causes internal compiler error for gcc 4.7.1 00087 // f[j]+=that->f[j]*PHI(ths->n[0],that->t[j]-((double)l)/((double)ths->n[0]),0); 00088 } 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 /* the pointers that->f and that->f_hat have been modified by the solver */ 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 /* PHI has compact support */ 00129 if(fabs(that->t[j]-((double)l)/((double)ths->n[0]))<that->plan.m/((double)ths->n[0])) 00130 that->f[j]*=PHI(ths->n[0],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*KPI*_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],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 /* the pointers ths->f and ths->f_hat have been modified by the solver */ 00179 ths->plan.f = ths->f; 00180 ths->plan.f_hat = ths->f_hat; 00181 00182 nfft_finalize(&ths->plan); 00183 } 00184 00185 /* 00186 * mri_inh_3d 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 /* the pointers that->f has been modified by the solver */ 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 /* PHI has compact support */ 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(ths->n[0],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],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 /* the pointers that->f has been modified by the solver */ 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],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 /* PHI has compact support */ 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(ths->n[0],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 }