/home/pclz-02/oschulz/Data/Science/Projects/Gerda/Software/blitzwave/src/arrayTools.h Source File

00001 // Copyright (C) 2003-2008 Oliver Schulz <oliver.schulz@tu-dortmund.de>
00002 
00003 // This program is free software; you can redistribute it and/or modify
00004 // it under the terms of the GNU General Public License as published by
00005 // the Free Software Foundation; either version 2 of the License, or
00006 // (at your option) any later version.
00007 //
00008 // This program 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
00011 // GNU General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU General Public License
00014 // along with this program; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
00016 
00017 
00018 #ifndef ARRAYTOOLS_H
00019 #define ARRAYTOOLS_H
00020 
00021 #include <cassert>
00022 #include <iostream>
00023 #include <blitz/array.h>
00024 
00025 namespace bwave {
00026 
00027 
00038 
00039 enum ExtensionMode {
00040         ZERO_EXT=0,
00041         CONSTANT_EXT=1,
00042         SYMMETRIC_EXT=2,
00043         SYMMETRIC2_EXT=3,
00044         CYCLIC_EXT=2
00045 };
00046 
00047 
00048 template<class tp_Type, int tp_rank>
00049         blitz::Array<tp_Type,1> slice(blitz::Array<tp_Type, tp_rank> &array,
00050         int dimension, blitz::TinyVector<int, tp_rank> position);
00051 
00052 template<class tp_Type>
00053         void fill(blitz::Array<tp_Type,1> &trg, const blitz::Array<tp_Type,1>
00054                 &src, int origin=0, ExtensionMode em=ZERO_EXT);
00055 
00056 
00057 template<class tp_Type> static tp_Type rollR(tp_Type x, tp_Type i) { assert(false); return 0; }
00058 template<> inline int rollR(int x, int i) { return x >> i; }
00059 template<> inline short rollR(short x, short i) { return x >> i; }
00060 template<> inline char rollR(char x, char i) { return x >> i; }
00061 
00062 template<class tp_Type> static tp_Type div2BS(tp_Type divisor) { return -1; }
00063 template<> inline int div2BS(int divisor) { return (divisor>0) && !(divisor & (divisor-1)) ? int( log((double)divisor) / log(2.0) ) : -1; }
00064 template<> inline short div2BS(short divisor) { return (divisor>0) && !(divisor & (divisor-1)) ? short( log((double)divisor) / log(2.0) ) : -1; }
00065 template<> inline char div2BS(char divisor) { return (divisor>0) && !(divisor & (divisor-1)) ? char( log((double)divisor) / log(2.0) ) : -1; }
00066 
00067 
00068 template<class tp_Type, int tp_size> class GenFilter {
00069 protected:
00070         tp_Type m_coeffs[tp_size];
00071         int m_origin;
00072         tp_Type m_divisor;
00073 
00074         template<class tp_Type2, int tp_size2> static tp_Type2 filterSP
00075                 (const GenFilter<tp_Type2, tp_size2> &coeffs, const blitz::Array<tp_Type2,1> &in, int i)
00076         {
00077                 tp_Type2 result=0;
00078                 for (int j=0; j<tp_size2; ++j) result += coeffs(j)*in(i+j);
00079                 return result;
00080         }
00081         
00082         template<class tp_Type2> static tp_Type2 filterSP
00083                 (const GenFilter<tp_Type2, 1> &coeffs, const blitz::Array<tp_Type2,1> &in, int i)
00084         { return coeffs(0)*in(i+0); }
00085 
00086         template<class tp_Type2> static tp_Type2 filterSP
00087                 (const GenFilter<tp_Type2, 2> &coeffs, const blitz::Array<tp_Type2,1> &in, int i)
00088         { return coeffs(0)*in(i+0) + coeffs(1)*in(i+1); }
00089 
00090         template<class tp_Type2> static tp_Type2 filterSP
00091                 (const GenFilter<tp_Type2, 3> &coeffs, const blitz::Array<tp_Type2,1> &in, int i)
00092         { return coeffs(0)*in(i+0) + coeffs(1)*in(i+1) + coeffs(2)*in(i+2); }
00093 
00094         template<class tp_Type2> static tp_Type2 filterSP
00095                 (const GenFilter<tp_Type2, 4> &coeffs, const blitz::Array<tp_Type2,1> &in, int i)
00096         { return coeffs(0)*in(i+0) + coeffs(1)*in(i+1) + coeffs(2)*in(i+2) + coeffs(3)*in(i+3); }
00097 
00098         template<class tp_Type2> static tp_Type2 filterSP
00099                 (const GenFilter<tp_Type2, 5> &coeffs, const blitz::Array<tp_Type2,1> &in, int i)
00100         { return coeffs(0)*in(i+0) + coeffs(1)*in(i+1) + coeffs(2)*in(i+2) + coeffs(3)*in(i+3) + coeffs(4)*in(i+4); }
00101 
00102         template<class tp_Type2> static tp_Type2 filterSP
00103                 (const GenFilter<tp_Type2, 6> &coeffs, const blitz::Array<tp_Type2,1> &in, int i)
00104         { return coeffs(0)*in(i+0) + coeffs(1)*in(i+1) + coeffs(2)*in(i+2) + coeffs(3)*in(i+3) + coeffs(4)*in(i+4) + coeffs(5)*in(i+5); }
00105 
00106 public:
00107         static void set(tp_Type &target, const tp_Type &value) { target = value; }
00108         static void inc(tp_Type &target, const tp_Type &value) { target += value; }
00109         static void dec(tp_Type &target, const tp_Type &value) { target -= value; }
00110         
00111         tp_Type operator()(unsigned i) const { return m_coeffs[i]; }
00112         tp_Type& operator()(unsigned i) { return m_coeffs[i]; }
00113 
00114         tp_Type divisor() const { return m_divisor; }
00115         tp_Type divisor(int newDiv) { m_divisor=newDiv; return divisor(); }
00116 
00117         template<void op(tp_Type &target, const tp_Type &value)>
00118                 void apply(const blitz::Array<tp_Type,1> &in, blitz::Array<tp_Type,1> &out, ExtensionMode be) const
00119         {
00120                 using namespace std;
00121                 assert( in.lbound()(0)==0 );
00122                 assert( out.lbound()(0)==0 );
00123                 
00124                 tp_Type bitShift = div2BS(m_divisor);
00125                 
00126                 const int nI=in.rows();
00127                 const int nO=out.rows();
00128                 const int fb = m_origin;
00129                 const int fe = m_origin+tp_size-1;
00130                 const int mainB = max(0, -fb), mainE = min(nI, nO)+min(0, -fe);
00131 
00132                 blitz::Array<tp_Type, 1> dummy(max(1, tp_size-1 + max(mainB, nO-mainE) ));
00133 
00134                 if (mainB>0) {
00135                         fill(dummy, in, mainB, be);
00136 
00137                         int inI=0;
00138                         if (bitShift>=0) for (int outI=0; outI<mainB; ++outI) {
00139                                 op( out(outI), rollR(filterSP(*this, dummy, inI), bitShift) );
00140                                 ++inI;
00141                         } else for (int outI=0; outI<mainB; ++outI) {
00142                                 op( out(outI), filterSP(*this, dummy, inI) / m_divisor );
00143                                 ++inI;
00144                         }
00145                 }
00146                 
00147                 int inI=mainB+fb;
00148                 if (bitShift>=0) for (int outI=mainB; outI<mainE; ++outI) {
00149                         op( out(outI), rollR(filterSP(*this, in, inI), bitShift) );
00150                         ++inI;
00151                 } else for (int outI=mainB; outI<mainE; ++outI) {
00152                         op( out(outI), filterSP(*this, in, inI) / m_divisor );
00153                         ++inI;
00154                 }
00155 
00156                 if (mainE<nO) {
00157                         fill(dummy, in, tp_size-1-nI, be);
00158 
00159                         int inI=0;
00160                         if (bitShift>=0) for (int outI=mainE; outI<nO; ++outI) {
00161                                 op( out(outI), rollR(filterSP(*this, dummy, inI), bitShift) );
00162                                 ++inI;
00163                         } else for (int outI=mainE; outI<nO; ++outI) {
00164                                 op( out(outI), filterSP(*this, dummy, inI) / m_divisor );
00165                                 ++inI;
00166                         }
00167                 }
00168         }
00169 
00170         void apply(const blitz::Array<tp_Type,1> &in, blitz::Array<tp_Type,1> &out, ExtensionMode be=ZERO_EXT) const
00171                 { apply<GenFilter<tp_Type, tp_size>::set>(in, out, be); }
00172 
00173         void applyAdd(const blitz::Array<tp_Type,1> &in, blitz::Array<tp_Type,1> &out, ExtensionMode be=ZERO_EXT) const
00174                 { apply<GenFilter<tp_Type, tp_size>::inc>(in, out, be); }
00175 
00176         void applySub(const blitz::Array<tp_Type,1> &in, blitz::Array<tp_Type,1> &out, ExtensionMode be=ZERO_EXT) const
00177                 { apply<GenFilter<tp_Type, tp_size>::dec>(in, out, be); }
00178 
00179         GenFilter(int origin) : m_origin(origin) {}
00180 
00181         GenFilter(int origin, const tp_Type *x, tp_Type divis)
00182                 : m_origin(origin), m_divisor(divis)
00183         {
00184                 for (int i=0; i<tp_size; ++i) m_coeffs[i] = x[i];
00185         }
00186         
00187         ~GenFilter() {}
00188 };
00189 
00190 
00191 double runTime();
00192 
00193 void readPNM(const std::string &filename, blitz::Array<unsigned char, 3> &target);
00194 
00195 void writePNM(const std::string &filename, blitz::Array<unsigned char, 3> &source);
00196 
00197 
00198 
00199 // -- template implementations: ---------------------------------------
00200 
00201 template<class tp_Type, int tp_rank>
00202         blitz::Array<tp_Type,1> slice(blitz::Array<tp_Type, tp_rank> &array,
00203         int dimension, blitz::TinyVector<int, tp_rank> position)
00204 {
00205         using namespace blitz;
00206 
00207         GeneralArrayStorage<1> newStorage;
00208         newStorage.ordering()(0) = firstDim;
00209         newStorage.ascendingFlag()(0) = array.stride(dimension)>=0;
00210         newStorage.base()(0) = array.base(dimension);
00211 
00212         position(dimension) = newStorage.ascendingFlag()(0) ?
00213                 array.lbound(dimension) : array.ubound(dimension);
00214         tp_Type *data = &(array(position));
00215 
00216         TinyVector<int, 1> newShape; newShape(0) = array.shape()(dimension);
00217         TinyVector<int, 1> newStride; newStride(0) = array.stride(dimension);
00218 
00219         return Array<tp_Type, 1>(data, newShape, newStride, neverDeleteData, newStorage);
00220 }
00221 
00222 
00223 template<class tp_Type>
00224         void fill(blitz::Array<tp_Type,1> &trg, const blitz::Array<tp_Type,1>
00225                 &src, int origin, ExtensionMode em)
00226 {
00227         using namespace blitz;
00228 
00229         // Symmetric and cyclic extension not implemented yet
00230         assert ( em!= SYMMETRIC_EXT );
00231         assert ( em!= SYMMETRIC2_EXT );
00232         assert ( em!= CYCLIC_EXT );
00233 
00234         //const int nSrc = src.rows();
00235         //const int nTrg = trg.rows();
00236         const int lSrc = src.lbound()(0);
00237         const int lTrg = trg.lbound()(0);
00238         const int uSrc = src.ubound()(0);
00239         const int uTrg = trg.ubound()(0);
00240 
00241         const int stop = uTrg+1;
00242         const int start = lTrg;
00243         const int sSrc = max(start-origin, 0)+lSrc;
00244         const int stage1 = min(max(origin, start), stop);
00245         const int stage2 = min(stage1+max(uSrc-sSrc+1,0), stop);
00246 
00247         tp_Type fillConst = (em==CONSTANT_EXT ? src(lSrc) : 0);
00248         for (int i=start; i<stage1; ++i) trg(i) = fillConst;
00249 
00250         const int shift = sSrc - stage1;
00251         for (int i=stage1; i<stage2; ++i) trg(i) = src(i+shift);
00252 
00253         fillConst = (em==CONSTANT_EXT ? src(uSrc) : 0);
00254         for (int i=stage2; i<stop; ++i) trg(i) = fillConst;
00255 }
00256 
00257 
00258 } // namespace bwave
00259 
00260 #endif