00001
00022 #ifndef DGVRADONTRANSFORM_H
00023 #define DGVRADONTRANSFORM_H
00024
00025 #include <typeinfo>
00026
00027 #include "DGVAliases.h"
00028 #include "DGVFourierTransform.h"
00029
00030 extern "C"
00031 {
00032 #include "radon.h"
00033 }
00034
00043 template <typename type>
00044 class DGVRadonTransform
00045 {
00046 public:
00052 DGVRadonTransform();
00056 virtual ~DGVRadonTransform();
00058
00066 void DiscreteRT(Array<type,2> &field, Array<type,2> &bins);
00073 void iDiscreteRT(Array<type,2> &bins, Array<type,2> &result, bool norm = true);
00074
00080 void FastRT(Array<type,2> &field, Array<type,2> &bins);
00087 void iFastRT(Array<type,2> &bins, Array<type,2> &result, bool norm = true);
00088
00092 void extractSlices(Array<complex<type>,2> &field, Array<complex<type>,2> &slices);
00093 void extractSlices(Array<type,2> &field, Array<type,2> &slices);
00097 void getSlice(int m, Array<complex<type>,2> &field, Array<complex<type>,1> &slice);
00098 void getSlice(int m, Array<type,2> &field, Array<type,1> &slice);
00102 void setSlice(int m, Array<complex<type>,2> &field, Array<complex<type>,1> &slice);
00103 void setSlice(int m, Array<type,2> &field, Array<type,1> &slice);
00105
00111 inline void init()
00112 { Isum = initializer; }
00116 inline void setIsum(type sum)
00117 { Isum = sum; }
00121 inline type getIsum()
00122 { return Isum; }
00127 void ISumFromBins(Array<type,2> &bins);
00132 inline void setInitializer(type value)
00133 { initializer = value; }
00138 inline type getInitializer()
00139 { return initializer; }
00143 bool isTypeInteger();
00145
00146 protected:
00147 type Isum;
00148 type initializer;
00149 };
00150
00151 template <typename type>
00152 DGVRadonTransform<type>::DGVRadonTransform()
00153 {
00154 initializer = blitz::zero(initializer);
00155 }
00156
00157 template <typename type>
00158 DGVRadonTransform<type>::~DGVRadonTransform()
00159 {
00160
00161 }
00162
00163 template <typename type>
00164 void DGVRadonTransform<type>::DiscreteRT(Array<type,2> &field, Array<type,2> &bins)
00165 {
00166 const int p = field.cols(), translate_x = 1;
00167 int translate_y, y, k;
00168 bins.resizeAndPreserve(p+1,p);
00169 Isum = initializer;
00170 bins = initializer;
00171
00172
00173 for(int j = 0; j < p; j ++)
00174 {
00175 y = 0;
00176 for(int m = 0; m < p; m++)
00177 {
00178 y += j;
00179 k = (m + translate_x)%p;
00180 for(translate_y = 0; translate_y < p; translate_y ++)
00181 bins(j,translate_y) += field(k, (translate_y + y)%p);
00182 }
00183 }
00184
00185 for (int j = 0; j < p; j ++)
00186 for(k = 0; k < p; k ++)
00187 {
00188 bins(p,j) += field(j,k);
00189 Isum += field(j,k);
00190 }
00191 }
00192
00193 template <typename type>
00194 void DGVRadonTransform<type>::iDiscreteRT(Array<type,2> &bins, Array<type,2> &result, bool norm)
00195 {
00196 const int p = bins.cols(), translate_x = 1;
00197 int translate_y, y, k;
00198 result.resizeAndPreserve(p,p);
00199 result = initializer;
00200
00201 if(Isum == 0)
00202 ISumFromBins(bins);
00203
00204
00205 for(int j = 0; j < p; j ++)
00206 {
00207 y = 0;
00208 for(int m = 0; m < p; m++)
00209 {
00210 y += p-j;
00211 k = (m + translate_x)%p;
00212 for(translate_y = 0; translate_y < p; translate_y ++)
00213 result(j,translate_y) += bins(k, (translate_y + y)%p);
00214 }
00215 }
00216
00217 for(int j = 0; j < p; j ++)
00218 for(k = 0; k < p; k++)
00219 {
00220 result(j,k) += bins(p,j);
00221 result(j,k) = static_cast<type>(result(j,k) - Isum);
00222 if(norm)
00223 result(j,k) /= p;
00224 }
00225 }
00226
00227 template <typename type>
00228 void DGVRadonTransform<type>::FastRT(Array<type,2> &field, Array<type,2> &bins)
00229 {
00230 const int n = field.cols();
00231
00232 Isum = initializer;
00233
00234 if(n%2 == 1)
00235 {
00236 bins.resizeAndPreserve(n+1,n);
00237 bins = initializer;
00238
00239 if( isTypeInteger() )
00240 frt( reinterpret_cast<nttw_integer *>(field.data()), reinterpret_cast<nttw_integer *>(bins.data()), n );
00241 else
00242 frt_double( reinterpret_cast<double *>(field.data()), reinterpret_cast<double *>(bins.data()), n );
00243 }
00244 else
00245 {
00246 bins.resizeAndPreserve(n+n/2,n);
00247 bins = initializer;
00248
00249 if( isTypeInteger() )
00250 frt_dyadic( reinterpret_cast<nttw_integer *>(field.data()), reinterpret_cast<nttw_integer *>(bins.data()), n );
00251 else
00252 frt_dyadic_double( reinterpret_cast<double *>(field.data()), reinterpret_cast<double *>(bins.data()), n );
00253 }
00254 }
00255
00256 template <typename type>
00257 void DGVRadonTransform<type>::iFastRT(Array<type,2> &bins, Array<type,2> &result, bool norm)
00258 {
00259 const int n = bins.cols();
00260
00261 if(Isum == 0)
00262 ISumFromBins(bins);
00263
00264 result.resizeAndPreserve(n,n);
00265 result = initializer;
00266
00267 if(n == bins.rows()-1)
00268 {
00269 if( isTypeInteger() )
00270 ifrt( reinterpret_cast<nttw_integer *>(bins.data()), reinterpret_cast<nttw_integer *>(result.data()), n, norm );
00271 else
00272 ifrt_double( reinterpret_cast<double *>(bins.data()), reinterpret_cast<double *>(result.data()), n, norm );
00273 }
00274 else
00275 {
00276 if( isTypeInteger() )
00277 ifrt_dyadic( reinterpret_cast<nttw_integer *>(bins.data()), reinterpret_cast<nttw_integer *>(result.data()), n, norm );
00278 else
00279 ifrt_dyadic_double( reinterpret_cast<double *>(bins.data()), reinterpret_cast<double *>(result.data()), n, norm );
00280 }
00281 }
00282
00283 template <typename type>
00284 void DGVRadonTransform<type>::extractSlices(Array<complex<type>,2> &field, Array<complex<type>,2> &slices)
00285 {
00286 const int p = field.cols();
00287 DGVFourierTransform<type,1> Fourier1D;
00288 Array<complex<type>,1> slice(p), projection(p);
00289
00291 slices.resize(p+1,p);
00292 for(int m = 0; m < p+1; m ++)
00293 {
00294 getSlice(m,field,slice);
00295 Fourier1D.iFFT(slice,projection);
00296 Fourier1D.normalise(projection);
00297 slices(m,Range::all()) = projection;
00298 }
00299 }
00300
00301 template <typename type>
00302 void DGVRadonTransform<type>::extractSlices(Array<type,2> &field, Array<type,2> &slices)
00303 {
00304 const int p = field.cols();
00305 Array<type,1> slice(p);
00306
00308 slices.resize(p+1,p);
00309 for(int m = 0; m < p+1; m ++)
00310 {
00311 getSlice(m,field,slice);
00312 slices(m,Range::all()) = slice;
00313 }
00314 }
00315
00316 template <typename type>
00317 void DGVRadonTransform<type>::getSlice(int m, Array<complex<type>,2> &field, Array<complex<type>,1> &slice)
00318 {
00319 const int p = slice.size();
00320 int index = 0;
00321
00322 if(m < p && m >= 0)
00323 {
00324 slice(0) = field(0,0);
00325 for(int k = 1; k < p; k ++)
00326 {
00327 index = ( p-(k*m)%p )%p;
00328 slice(k) = field(index,k);
00329 }
00330 }
00331 else if(m == p)
00332 slice = field(Range::all(),0);
00333 }
00334
00335 template <typename type>
00336 void DGVRadonTransform<type>::getSlice(int m, Array<type,2> &field, Array<type,1> &slice)
00337 {
00338 const int p = slice.size();
00339 int index = 0;
00340
00341 if(m < p && m >= 0)
00342 {
00343 slice(0) = field(0,0);
00344 for(int k = 1; k < p; k ++)
00345 {
00346 index = ( p-(k*m)%p )%p;
00347 slice(k) = field(index,k);
00348 }
00349 }
00350 else if(m == p)
00351 slice = field(Range::all(),0);
00352 }
00353
00354 template <typename type>
00355 void DGVRadonTransform<type>::setSlice(int m, Array<complex<type>,2> &field, Array<complex<type>,1> &slice)
00356 {
00357 const int p = slice.size();
00358 int index = 0;
00359
00360 if(m < p && m >= 0)
00361 {
00362 field(0,0) += slice(0);
00363 for(int k = 1; k < p; k ++)
00364 {
00365 index = ( p-(k*m)%p )%p;
00366 field(index,k) = slice(k);
00367 }
00368 }
00369 else if(m == p)
00370 field(Range::all(),0) = slice;
00371 }
00372
00373 template <typename type>
00374 void DGVRadonTransform<type>::setSlice(int m, Array<type,2> &field, Array<type,1> &slice)
00375 {
00376 const int p = slice.size();
00377 int index = 0;
00378
00379 if(m < p && m >= 0)
00380 {
00381 field(0,0) += slice(0);
00382 for(int k = 1; k < p; k ++)
00383 {
00384 index = ( p-(k*m)%p )%p;
00385 field(index,k) = slice(k);
00386 }
00387 }
00388 else if(m == p)
00389 field(Range::all(),0) = slice;
00390 }
00391
00392 template <typename type>
00393 void DGVRadonTransform<type>::ISumFromBins(Array<type,2> &bins)
00394 {
00395 Isum = initializer;
00396 for(int k = 0; k < bins.cols(); k ++)
00397 Isum += bins(0,k);
00398 }
00399
00400 template <typename type>
00401 bool DGVRadonTransform<type>::isTypeInteger()
00402 {
00403 if( typeid(type).name() == "int" || typeid(type).name() == "long" || typeid(type).name() == "short" || typeid(type).name() == "char" ||
00404 typeid(type).name() == "long int" || typeid(type).name() == "long long" || typeid(type).name() == "long short" || typeid(type).name() == "unsigned char" ||
00405 typeid(type).name() == "unsigned int" || typeid(type).name() == "unsigned long" || typeid(type).name() == "unsigned short" || typeid(type).name() == "unsigned long long")
00406 return true;
00407 else
00408 return false;
00409 }
00410
00411 #endif // DGVRADONTRANSFORM_H