DGV: include/DGVRadonTransform.h Source File

Go to the documentation of this file.
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" //FRTW Library
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     //Get Projections where j denotes projection number
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     //Get the last projection of rows
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     //Get Projections where j denotes projection number
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