olena_ dwt.hh source file.txt

DWT变换源代码

dwt.hh


00001 // Copyright (C) 2001, 2002, 2003, 2004  EPITA Research and Development Laboratory
00002 //
00003 // This file is part of the Olena Library.  This library is free
00004 // software; you can redistribute it and/or modify it under the terms
00005 // of the GNU General Public License version 2 as published by the
00006 // Free Software Foundation.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU General Public License
00014 // along with this library; see the file COPYING.  If not, write to
00015 // the Free Software Foundation, 59 Temple Place - Suite 330, Boston,
00016 // MA 02111-1307, USA.
00017 //
00018 // As a special exception, you may use this file as part of a free
00019 // software library without restriction.  Specifically, if other files
00020 // instantiate templates or use macros or inline functions from this
00021 // file, or you compile this file and link it with other files to
00022 // produce an executable, this file does not by itself cause the
00023 // resulting executable to be covered by the GNU General Public
00024 // License.  This exception does not however invalidate any other
00025 // reasons why the executable file might be covered by the GNU General
00026 // Public License.
00027 
00028 #ifndef OLENA_TRANSFORMS_DWT_HH
00029 # define OLENA_TRANSFORMS_DWT_HH
00030 
00031 # include <oln/config/system.hh>
00032 # include <oln/basics.hh>
00033 
00034 # include <ntg/basics.hh>
00035 # include <ntg/float.hh>
00036 # include <ntg/utils/cast.hh>
00037 # include <mlc/array/1d.hh>
00038 
00039 // FIXME: this file has not been adjusted to the coding style since it
00040 // will be completely rewritten in next release.
00041 
00042 namespace oln {
00043 
00047   namespace transforms {
00048 
00049     // macros used to define all the wavelets coefficients
00050 
00051 # define Wavelet_coeffs_definition(Name, Type, Size) \
00052     struct Name : public oln::internal::wavelet_coeffs_<Type, Size, Name > \
00053     { \
00054       Name()
00055 
00056 # define Wavelet_coeffs_begin \
00057       h = (wc_start =
00058 
00059 # define Wavelet_coeffs_end \
00060         , end); \
00061       init(); \
00062     }
00063 
00067     typedef enum {
00068       dwt_std,
00069       dwt_non_std
00070     } dwt_transform_type;
00071 
00072   } // end of namespace transforms
00073 
00074   namespace internal
00075   {
00077     static const ntg::float_d ln_2_ = 0.6931471805599453092;
00078 
00086     template <class T, unsigned N, class Self>
00087     struct wavelet_coeffs_
00088     {
00089       typedef T         value_t; 
00090       typedef Self      self_t; 
00091 
00092     public:
00093 
00094       // accessors
00096       const value_t getG(unsigned i) const { return g[i]; }
00098       const value_t getInvG(unsigned i) const { return ig[i]; }
00100       const value_t getH(unsigned i) const { return h[i]; }
00102       const value_t getInvH(unsigned i) const { return ih[i]; }
00104       const unsigned size() const { return size_; }
00105 
00106     protected:
00112       void init()
00113       {
00114         for (unsigned i = 0; i < size_; i += 2) {
00115           g[i] = -h[size_ - 1 - i];
00116           g[i + 1] = h[size_ - 2 - i];
00117           ig[size_ - 1 - i] = g[i + 1];
00118           ig[size_ - 2 - i] = h[i + 1];
00119           ih[size_ - 1 - i] = g[i];
00120           ih[size_ - 2 - i] = h[i];
00121         }
00122       }
00123 
00127       wavelet_coeffs_(){}
00128 
00134       ~wavelet_coeffs_()
00135       {
00136         mlc::is_false<N % 2>::ensure();
00137       }
00138 
00139       mlc::array1d< mlc::array1d_info<N>, value_t>      h; 
00140       mlc::internal::array1d_start_<value_t>            wc_start; 
00141 
00142     private:
00143 
00144       value_t   g[N]; 
00145       value_t   ih[N]; 
00146       value_t   ig[N]; 
00147       enum {
00148         size_ = N
00149       };
00150     };
00151 
00164     template <class I, class K>
00165     void dwt_transform_step_(abstract::image<I>& im,
00166                              const oln_point_type(I)& p_,
00167                              const unsigned d,
00168                              const unsigned n,
00169                              const K& coeffs)
00170     {
00171       assertion(n >= coeffs.size());
00172 
00173       const unsigned    half = n >> 1;
00174       unsigned          lim = n + 1 - coeffs.size();
00175       ntg::float_d*     tmp = new ntg::float_d[n];
00176       oln_point_type(I) p(p_);
00177       unsigned          i, j, k;
00178 
00179       for (i = 0, j = 0; j < lim; j += 2, i++) {
00180         tmp[i] = 0;
00181         tmp[i + half] = 0;
00182         for (k = 0; k < coeffs.size(); k++) {
00183           p.nth(d) = j + k;
00184           tmp[i] += im[p] * coeffs.getH(k);
00185           tmp[i + half] += im[p] * coeffs.getG(k);
00186         }
00187       }
00188       for (; i < half; j += 2, i++) {
00189         tmp[i] = 0;
00190         tmp[i + half] = 0;
00191         for (unsigned k = 0; k < coeffs.size(); k++) {
00192           p.nth(d) = (j + k) % n;
00193           tmp[i] += im[p] * coeffs.getH(k);
00194           tmp[i + half] += im[p] * coeffs.getG(k);
00195         }
00196       }
00197 
00198       for (i = 0; i < n; i++) {
00199         p.nth(d) = i;
00200         im[p] = tmp[i];
00201       }
00202 
00203       delete[] tmp;
00204     }
00205 
00218     template <class I, class K>
00219     void dwt_transform_inv_step_(abstract::image<I>& im,
00220                                  const oln_point_type(I)& p_,
00221                                  const unsigned d,
00222                                  const unsigned n,
00223                                  const K& coeffs)
00224     {
00225       assertion(n >= coeffs.size());
00226 
00227       const unsigned half = n >> 1;
00228       unsigned lim = coeffs.size() - 2;
00229       ntg::float_d* tmp = new ntg::float_d[n];
00230      oln_point_type(I) p(p_), q(p_);
00231       unsigned i, j, k, l;
00232 
00233       for (i = half - lim / 2, j = 0; j < lim; i++, j += 2) {
00234         tmp[j] = 0;
00235         tmp[j + 1] = 0;
00236         for (k = 0, l = 0; l < coeffs.size(); k++, l += 2) {
00237           p.nth(d) = (i + k) % half;
00238           q.nth(d) = p.nth(d) + half;
00239           tmp[j] += im[p] * coeffs.getInvH(l) + im[q] * coeffs.getInvH(l + 1);
00240           tmp[j + 1] += im[p] * coeffs.getInvG(l) + im[q] * coeffs.getInvG(l + 1);
00241         }
00242       }
00243       lim = n - 1;
00244       for (i = 0; j < lim; i++, j += 2) {
00245         tmp[j] = 0;
00246         tmp[j + 1] = 0;
00247         for (k = 0, l = 0; l < coeffs.size(); k++, l += 2) {
00248           p.nth(d) = i + k;
00249           q.nth(d) = p.nth(d) + half;
00250           tmp[j] += im[p] * coeffs.getInvH(l) + im[q] * coeffs.getInvH(l + 1);
00251           tmp[j + 1] += im[p] * coeffs.getInvG(l) + im[q] * coeffs.getInvG(l + 1);
00252         }
00253       }
00254 
00255       for (i = 0; i < n; i++) {
00256         p.nth(d) = i;
00257         im[p] = tmp[i];
00258       }
00259 
00260       delete[] tmp;
00261     }
00262 
00263     // Functions used to iterate over all dimensions except one
00267     typedef enum {
00268       dwt_fwd,
00269       dwt_bwd
00270     } dwt_transform_dir_;
00271 
00279     template <unsigned dim, unsigned skip,
00280               unsigned current = dim>
00281     struct dim_skip_iterate_rec_
00282     {
00286       template <class I, class K>
00287       static void doit(abstract::image<I>& im,
00288                       oln_point_type(I)& p,
00289                        const unsigned l1,
00290                        const unsigned l2,
00291                        const K& coeffs,
00292                        dwt_transform_dir_ d)
00293       {
00294         if (current != skip) {
00295           for (unsigned i = 0; i < l2; i++) {
00296             p.nth(current - 1) = i;
00297             dim_skip_iterate_rec_<dim, skip, current - 1>::
00298 	      doit(im, p, l1, l2, coeffs, d);
00299           }
00300           p.nth(current - 1) = 0;
00301         }
00302         else
00303           dim_skip_iterate_rec_<dim, skip, current - 1>::
00304 	    doit(im, p, l1, l2, coeffs, d);
00305       }