NFFT 3.5.3alpha
mri.c
1/*
2 * Copyright (c) 2002, 2017 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#include "config.h"
20
21#include <string.h>
22#include <math.h>
23#ifdef HAVE_COMPLEX_H
24#include <complex.h>
25#endif
26#include "nfft3.h"
27#include "infft.h"
28
34typedef struct window_funct_plan_ {
35 int d;
36 int m;
37 int n[1];
38 double sigma[1];
39 double *b;
41
45static void window_funct_init(window_funct_plan* ths, int m, int n, double sigma) {
46 ths->d=1;
47 ths->m=m;
48 ths->n[0]=n;
49 ths->sigma[0]=sigma;
50 WINDOW_HELP_INIT
51}
52
53/*
54 * mri_inh_2d1d
55 */
56
57void mri_inh_2d1d_trafo(mri_inh_2d1d_plan *that) {
58 int l,j;
59 double _Complex *f = (double _Complex*) nfft_malloc(that->M_total*sizeof(double _Complex));
60 double _Complex *f_hat = (double _Complex*) nfft_malloc(that->N_total*sizeof(double _Complex));
61
63 window_funct_init(ths,that->plan.m,that->N3,that->sigma3);
64
65 /* the pointers that->f and that->f_hat have been modified by the solver */
66 that->plan.f = that->f;
67 that->plan.f_hat = that->f_hat;
68
69
70 memset(f,0,that->M_total*sizeof(double _Complex));
71 for(j=0;j<that->N_total;j++)
72 {
73 f_hat[j]=that->f_hat[j];
74 }
75
76 for(l=-ths->n[0]/2;l<=ths->n[0]/2;l++) {
77 for(j=0;j<that->N_total;j++)
78 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);
79 nfft_trafo(&that->plan);
80 for(j=0;j<that->M_total;j++){
81 /* PHI has compact support */
82 if(fabs(that->t[j]-((double)l)/((double)ths->n[0]))<that->plan.m/((double)ths->n[0]))
83 {
84 double phi_val = PHI(ths->n[0],that->t[j]-((double)l)/((double)ths->n[0]),0);
85 f[j]+=that->f[j]*phi_val;
86// the line below causes internal compiler error for gcc 4.7.1
87// f[j]+=that->f[j]*PHI(ths->n[0],that->t[j]-((double)l)/((double)ths->n[0]),0);
88 }
89 }
90 for(j=0;j<that->N_total;j++)
91 that->f_hat[j]=f_hat[j];
92 }
93
94 nfft_free(that->plan.f);
95 that->f=f;
96 that->plan.f = that->f;
97
98 nfft_free(f_hat);
99
100 WINDOW_HELP_FINALIZE
101 nfft_free(ths);
102}
103
104void mri_inh_2d1d_adjoint(mri_inh_2d1d_plan *that) {
105 int l,j;
106 double _Complex *f = (double _Complex*) nfft_malloc(that->M_total*sizeof(double _Complex));
107 double _Complex *f_hat = (double _Complex*) nfft_malloc(that->N_total*sizeof(double _Complex));
108
110 window_funct_init(ths,that->plan.m,that->N3,that->sigma3);
111
112 memset(f_hat,0,that->N_total*sizeof(double _Complex));
113
114 /* the pointers that->f and that->f_hat have been modified by the solver */
115 that->plan.f = that->f;
116 that->plan.f_hat = that->f_hat;
117
118 for(j=0;j<that->M_total;j++)
119 {
120 f[j]=that->f[j];
121 }
122
123
124
125 for(l=-ths->n[0]/2;l<=ths->n[0]/2;l++) {
126
127 for(j=0;j<that->M_total;j++) {
128 /* PHI has compact support */
129 if(fabs(that->t[j]-((double)l)/((double)ths->n[0]))<that->plan.m/((double)ths->n[0]))
130 that->f[j]*=PHI(ths->n[0],that->t[j]-((double)l)/((double)ths->n[0]),0);
131 else
132 that->f[j]=0.0;
133 }
134 nfft_adjoint(&that->plan);
135 for(j=0;j<that->N_total;j++)
136 f_hat[j]+=that->f_hat[j]*cexp(2*KPI*_Complex_I*that->w[j]*((double)l));
137 for(j=0;j<that->M_total;j++)
138 that->f[j]=f[j];
139 }
140
141 for(j=0;j<that->N_total;j++)
142 {
143 f_hat[j] /= PHI_HUT(ths->n[0],ths->n[0]*that->w[j],0);
144 }
145
146 nfft_free(that->plan.f_hat);
147 that->f_hat=f_hat;
148 that->plan.f_hat = that->f_hat;
149
150 nfft_free(f);
151
152 WINDOW_HELP_FINALIZE
153 nfft_free(ths);
154}
155
156void mri_inh_2d1d_init_guru(mri_inh_2d1d_plan *ths, int *N, int M, int *n,
157 int m, double sigma, unsigned nfft_flags, unsigned fftw_flags) {
158
159 nfft_init_guru(&ths->plan,2,N,M,n,m,nfft_flags,fftw_flags);
160 ths->N3=N[2];
161 ths->sigma3=sigma;
162 ths->N_total = ths->plan.N_total;
163 ths->M_total = ths->plan.M_total;
164 ths->f = ths->plan.f;
165 ths->f_hat = ths->plan.f_hat;
166
167 ths->t = (double*) nfft_malloc(ths->M_total*sizeof(double));
168 ths->w = (double*) nfft_malloc(ths->N_total*sizeof(double));
169
170 ths->mv_trafo = (void (*) (void* ))mri_inh_2d1d_trafo;
171 ths->mv_adjoint = (void (*) (void* ))mri_inh_2d1d_adjoint;
172}
173
174void mri_inh_2d1d_finalize(mri_inh_2d1d_plan *ths) {
175 nfft_free(ths->t);
176 nfft_free(ths->w);
177
178 /* the pointers ths->f and ths->f_hat have been modified by the solver */
179 ths->plan.f = ths->f;
180 ths->plan.f_hat = ths->f_hat;
181
182 nfft_finalize(&ths->plan);
183}
184
185/*
186 * mri_inh_3d
187 */
188
189void mri_inh_3d_trafo(mri_inh_3d_plan *that) {
190 int l,j;
192 window_funct_init(ths,that->plan.m,that->N3,that->sigma3);
193
194 /* the pointers that->f has been modified by the solver */
195 that->plan.f =that->f ;
196
197
198
199 for(j=0;j<that->N_total;j++) {
200 for(l=-ths->n[0]/2;l<ths->n[0]/2;l++)
201 {
202 /* PHI has compact support */
203 if(fabs(that->w[j]-((double)l)/((double)ths->n[0]))<ths->m/((double)ths->n[0]))
204 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);
205 else
206 that->plan.f_hat[j*ths->n[0]+(l+ths->n[0]/2)]=0.0;
207 }
208 }
209
210 nfft_trafo(&that->plan);
211
212 for(j=0;j<that->M_total;j++)
213 {
214 that->f[j] /= PHI_HUT(ths->n[0],ths->n[0]*that->plan.x[3*j+2],0);
215 }
216
217 WINDOW_HELP_FINALIZE
218 nfft_free(ths);
219}
220
221void mri_inh_3d_adjoint(mri_inh_3d_plan *that) {
222 int l,j;
224 window_funct_init(ths,that->plan.m,that->N3,that->sigma3);
225
226 /* the pointers that->f has been modified by the solver */
227 that->plan.f =that->f ;
228
229 for(j=0;j<that->M_total;j++)
230 {
231 that->f[j] /= PHI_HUT(ths->n[0],ths->n[0]*that->plan.x[3*j+2],0);
232 }
233
234 nfft_adjoint(&that->plan);
235
236 for(j=0;j<that->N_total;j++) {
237 that->f_hat[j]=0.0;
238 for(l=-ths->n[0]/2;l<ths->n[0]/2;l++)
239 {
240 /* PHI has compact support */
241 if(fabs(that->w[j]-((double)l)/((double)ths->n[0]))<ths->m/((double)ths->n[0]))
242 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);
243 }
244 }
245
246
247 WINDOW_HELP_FINALIZE
248 nfft_free(ths);
249}
250
251void mri_inh_3d_init_guru(mri_inh_3d_plan *ths, int *N, int M, int *n,
252 int m, double sigma, unsigned nfft_flags, unsigned fftw_flags) {
253 ths->N3=N[2];
254 ths->sigma3=sigma;
255 nfft_init_guru(&ths->plan,3,N,M,n,m,nfft_flags,fftw_flags);
256 ths->N_total = N[0]*N[1];
257 ths->M_total = ths->plan.M_total;
258 ths->f = ths->plan.f;
259 ths->f_hat = (double _Complex*) nfft_malloc(ths->N_total*sizeof(double _Complex));
260 ths->w = (double*) nfft_malloc(ths->N_total*sizeof(double));
261
262 ths->mv_trafo = (void (*) (void* ))mri_inh_3d_trafo;
263 ths->mv_adjoint = (void (*) (void* ))mri_inh_3d_adjoint;
264}
265
266void mri_inh_3d_finalize(mri_inh_3d_plan *ths) {
267 nfft_free(ths->w);
268 nfft_free(ths->f_hat);
269 nfft_finalize(&ths->plan);
270}
Internal header file for auxiliary definitions and functions.
Header file for the nfft3 library.
void * nfft_malloc(size_t n)
void nfft_free(void *p)
NFFT_INT N_total
Total number of Fourier coefficients.
Definition nfft3.h:532
NFFT_INT M_total
Total number of samples.
Definition nfft3.h:532
void(* mv_adjoint)(void *)
Adjoint transform.
Definition nfft3.h:532
fftw_complex * f
Samples.
Definition nfft3.h:532
void(* mv_trafo)(void *)
Transform.
Definition nfft3.h:532
fftw_complex * f_hat
Fourier coefficients.
Definition nfft3.h:532
void(* mv_trafo)(void *)
Transform.
Definition nfft3.h:532
NFFT_INT M_total
Total number of samples.
Definition nfft3.h:532
fftw_complex * f_hat
Fourier coefficients.
Definition nfft3.h:532
fftw_complex * f
Samples.
Definition nfft3.h:532
void(* mv_adjoint)(void *)
Adjoint transform.
Definition nfft3.h:532
NFFT_INT N_total
Total number of Fourier coefficients.
Definition nfft3.h:532
window_funct_plan is a plan to use the window functions independent of the nfft
Definition mri.c:34