[299] | 1 | // ----------------------------------------------------------------------- |
---|
| 2 | // filter.cc: Defining a filter function for wavelet reconstruction. |
---|
| 3 | // ----------------------------------------------------------------------- |
---|
| 4 | // Copyright (C) 2006, Matthew Whiting, ATNF |
---|
| 5 | // |
---|
| 6 | // This program is free software; you can redistribute it and/or modify it |
---|
| 7 | // under the terms of the GNU General Public License as published by the |
---|
| 8 | // Free Software Foundation; either version 2 of the License, or (at your |
---|
| 9 | // option) any later version. |
---|
| 10 | // |
---|
| 11 | // Duchamp is distributed in the hope that it will be useful, but WITHOUT |
---|
| 12 | // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
---|
| 13 | // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
---|
| 14 | // for more details. |
---|
| 15 | // |
---|
| 16 | // You should have received a copy of the GNU General Public License |
---|
| 17 | // along with Duchamp; if not, write to the Free Software Foundation, |
---|
| 18 | // Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA |
---|
| 19 | // |
---|
| 20 | // Correspondence concerning Duchamp may be directed to: |
---|
| 21 | // Internet email: Matthew.Whiting [at] atnf.csiro.au |
---|
| 22 | // Postal address: Dr. Matthew Whiting |
---|
| 23 | // Australia Telescope National Facility, CSIRO |
---|
| 24 | // PO Box 76 |
---|
| 25 | // Epping NSW 1710 |
---|
| 26 | // AUSTRALIA |
---|
| 27 | // ----------------------------------------------------------------------- |
---|
[139] | 28 | #include <iostream> |
---|
| 29 | #include <sstream> |
---|
[393] | 30 | #include <duchamp/duchamp.hh> |
---|
| 31 | #include <duchamp/ATrous/filter.hh> |
---|
[3] | 32 | #include <math.h> |
---|
[258] | 33 | #include <vector> |
---|
[3] | 34 | |
---|
[378] | 35 | namespace duchamp |
---|
[3] | 36 | { |
---|
[365] | 37 | |
---|
[378] | 38 | Filter::Filter() |
---|
| 39 | { |
---|
| 40 | this->sigmaFactors.resize(3); |
---|
| 41 | for(int i=0;i<3;i++) this->sigmaFactors[i] = new std::vector<double>(20); |
---|
| 42 | this->loadSpline(); |
---|
| 43 | } |
---|
[365] | 44 | |
---|
[378] | 45 | Filter::Filter(const Filter& f) |
---|
| 46 | { |
---|
| 47 | operator=(f); |
---|
| 48 | } |
---|
[365] | 49 | |
---|
[378] | 50 | Filter& Filter::operator=(const Filter& f) |
---|
| 51 | { |
---|
| 52 | if(this==&f) return *this; |
---|
| 53 | this->name = f.name; |
---|
| 54 | this->filter1D = f.filter1D; |
---|
| 55 | this->maxNumScales = f.maxNumScales; |
---|
| 56 | this->sigmaFactors = f.sigmaFactors; |
---|
| 57 | return *this; |
---|
| 58 | } |
---|
[3] | 59 | |
---|
[378] | 60 | //----------------------------------------------------------------------- |
---|
[3] | 61 | |
---|
[378] | 62 | Filter::~Filter() |
---|
| 63 | { |
---|
[468] | 64 | this->filter1D.clear(); |
---|
| 65 | this->maxNumScales.clear(); |
---|
[378] | 66 | //for(int i=0;i<3;i++) delete this->sigmaFactors[i]; |
---|
[468] | 67 | //for(int i=0;i<3;i++) this->sigmaFactors[i]->clear(); |
---|
| 68 | // for(int i=0;i<3;i++){ |
---|
| 69 | // std::vector<double>* p = this->sigmaFactors[i]; |
---|
| 70 | // delete p; |
---|
| 71 | // } |
---|
| 72 | this->sigmaFactors.clear(); |
---|
[378] | 73 | } |
---|
| 74 | //----------------------------------------------------------------------- |
---|
| 75 | |
---|
| 76 | void Filter::define(int filtercode) |
---|
| 77 | { |
---|
| 78 | switch(filtercode) |
---|
| 79 | { |
---|
| 80 | case 2: |
---|
| 81 | this->loadTriangle(); |
---|
| 82 | break; |
---|
| 83 | case 3: |
---|
| 84 | this->loadHaar(); |
---|
| 85 | break; |
---|
| 86 | case 4: |
---|
| 87 | // this->loadTopHat(); |
---|
| 88 | case 1: |
---|
| 89 | default: |
---|
| 90 | if(filtercode!=1){ |
---|
| 91 | // if(filtercode==4) { |
---|
| 92 | // std::stringstream errmsg; |
---|
| 93 | // errmsg << "TopHat Wavelet not being used currently." |
---|
| 94 | // << "Using B3 spline instead.\n"; |
---|
| 95 | // duchampWarning("Filter::define", errmsg.str()); |
---|
| 96 | // } |
---|
| 97 | // else { |
---|
[913] | 98 | DUCHAMPWARN("Wavelet Filter", "Filter code " << filtercode << " undefined. Using B3 spline."); |
---|
[378] | 99 | // } |
---|
| 100 | } |
---|
| 101 | this->loadSpline(); |
---|
| 102 | break; |
---|
[3] | 103 | } |
---|
[378] | 104 | |
---|
| 105 | } |
---|
| 106 | //----------------------------------------------------------------------- |
---|
| 107 | |
---|
[846] | 108 | unsigned int Filter::getNumScales(long length) |
---|
[378] | 109 | { |
---|
[846] | 110 | unsigned int num; |
---|
[378] | 111 | switch(this->filter1D.size()){ |
---|
| 112 | case 5: |
---|
[855] | 113 | num = int(log(double(length-1))/M_LN2) - 1; |
---|
[3] | 114 | break; |
---|
[378] | 115 | case 3: |
---|
[855] | 116 | num = int(log(double(length-1))/M_LN2); |
---|
[378] | 117 | break; |
---|
| 118 | default: |
---|
[855] | 119 | num = 1 + int(log(double(length-1)/double(this->filter1D.size()-1))/M_LN2); |
---|
[378] | 120 | break; |
---|
[3] | 121 | } |
---|
[846] | 122 | return num; |
---|
[220] | 123 | } |
---|
[378] | 124 | //----------------------------------------------------------------------- |
---|
[220] | 125 | |
---|
[846] | 126 | unsigned int Filter::getMaxSize(int scale) |
---|
[378] | 127 | { |
---|
| 128 | switch(this->filter1D.size()){ |
---|
| 129 | case 5: |
---|
| 130 | return int(pow(2,scale+1)) + 1; |
---|
| 131 | break; |
---|
| 132 | case 3: |
---|
| 133 | return int(pow(2,scale)) + 1; |
---|
| 134 | break; |
---|
| 135 | default: |
---|
| 136 | return int(pow(2,scale-1))*(this->filter1D.size()-1) + 1; |
---|
| 137 | break; |
---|
| 138 | } |
---|
[220] | 139 | } |
---|
[378] | 140 | //----------------------------------------------------------------------- |
---|
[220] | 141 | |
---|
[378] | 142 | void Filter::loadSpline() |
---|
| 143 | { |
---|
| 144 | double filter[5] = {0.0625, 0.25, 0.375, 0.25, 0.0625}; |
---|
| 145 | this->name = "B3 spline function"; |
---|
| 146 | this->filter1D.resize(5); |
---|
| 147 | for(int i=0;i<5;i++) this->filter1D[i] = filter[i]; |
---|
| 148 | this->sigmaFactors.resize(3); |
---|
| 149 | this->maxNumScales.resize(3); |
---|
[3] | 150 | |
---|
[378] | 151 | this->maxNumScales[0] = 18; |
---|
| 152 | double sigmaFactors1D[19] = {1.00000000000,7.23489806e-1,2.85450405e-1,1.77947535e-1, |
---|
| 153 | 1.22223156e-1,8.58113122e-2,6.05703043e-2,4.28107206e-2, |
---|
| 154 | 3.02684024e-2,2.14024008e-2,1.51336781e-2,1.07011079e-2, |
---|
| 155 | 7.56682272e-3,5.35055108e-3,3.78341085e-3,2.67527545e-3, |
---|
| 156 | 1.89170541e-3,1.33763772e-3,9.45852704e-4}; |
---|
| 157 | this->sigmaFactors[0]->resize(19); |
---|
| 158 | for(int i=0;i<19;i++)(*this->sigmaFactors[0])[i] = sigmaFactors1D[i]; |
---|
[3] | 159 | |
---|
[859] | 160 | this->maxNumScales[1] = 12; |
---|
| 161 | double sigmaFactors2D[13] = {1.00000000000,8.90796310e-1,2.00663851e-1,8.55075048e-2, |
---|
[378] | 162 | 4.12474444e-2,2.04249666e-2,1.01897592e-2,5.09204670e-3, |
---|
[859] | 163 | 2.54566946e-3,1.27279050e-3,6.36389722e-4,3.18194170e-4, |
---|
| 164 | 1.59096999e-4}; |
---|
| 165 | this->sigmaFactors[1]->resize(13); |
---|
| 166 | for(int i=0;i<13;i++)(*this->sigmaFactors[1])[i] = sigmaFactors2D[i]; |
---|
[3] | 167 | |
---|
[378] | 168 | this->maxNumScales[2] = 7; |
---|
| 169 | double sigmaFactors3D[8] = {1.00000000000,9.56543592e-1,1.20336499e-1,3.49500154e-2, |
---|
| 170 | 1.18164242e-2,4.13233507e-3,1.45703714e-3,5.14791120e-4}; |
---|
| 171 | this->sigmaFactors[2]->resize(8); |
---|
| 172 | for(int i=0;i<12;i++)(*this->sigmaFactors[2])[i] = sigmaFactors3D[i]; |
---|
| 173 | } |
---|
| 174 | //----------------------------------------------------------------------- |
---|
[3] | 175 | |
---|
[378] | 176 | void Filter::loadTriangle() |
---|
| 177 | { |
---|
| 178 | double filter[3] = {1./4., 1./2., 1./4.}; |
---|
| 179 | this->filter1D.resize(3); |
---|
| 180 | for(int i=0;i<3;i++) this->filter1D[i] = filter[i]; |
---|
| 181 | this->name = "Triangle function"; |
---|
| 182 | this->sigmaFactors.resize(3); |
---|
| 183 | this->maxNumScales.resize(3); |
---|
[3] | 184 | |
---|
[378] | 185 | this->maxNumScales[0] = 18; |
---|
| 186 | double sigmaFactors1D[19] = {1.00000000000,6.12372436e-1,3.30718914e-1,2.11947812e-1, |
---|
| 187 | 1.45740298e-1,1.02310944e-1,7.22128185e-2,5.10388224e-2, |
---|
| 188 | 3.60857673e-2,2.55157615e-2,1.80422389e-2,1.27577667e-2, |
---|
| 189 | 9.02109930e-3,6.37887978e-3,4.51054902e-3,3.18942978e-3, |
---|
| 190 | 2.25527449e-3,1.59471988e-3,1.12763724e-4}; |
---|
| 191 | this->sigmaFactors[0]->resize(19); |
---|
| 192 | for(int i=0;i<19;i++)(*this->sigmaFactors[0])[i] = sigmaFactors1D[i]; |
---|
[3] | 193 | |
---|
[378] | 194 | this->maxNumScales[1] = 12; |
---|
| 195 | double sigmaFactors2D[13] = {1.00000000000,8.00390530e-1,2.72878894e-1,1.19779282e-1, |
---|
| 196 | 5.77664785e-2,2.86163283e-2,1.42747506e-2,7.13319703e-3, |
---|
| 197 | 3.56607618e-3,1.78297280e-3,8.91478237e-4,4.45738098e-4, |
---|
| 198 | 2.22868922e-4}; |
---|
| 199 | this->sigmaFactors[1]->resize(13); |
---|
| 200 | for(int i=0;i<12;i++)(*this->sigmaFactors[1])[i] = sigmaFactors2D[i]; |
---|
[3] | 201 | |
---|
[378] | 202 | this->maxNumScales[2] = 8; |
---|
| 203 | double sigmaFactors3D[9] = {1.00000000000,8.959544490e-1,1.92033014e-1,5.76484078e-2, |
---|
| 204 | 1.94912393e-2,6.812783870e-3,2.40175885e-3,8.48538128e-4, |
---|
| 205 | 2.99949455e-4}; |
---|
| 206 | this->sigmaFactors[2]->resize(9); |
---|
| 207 | for(int i=0;i<12;i++)(*this->sigmaFactors[2])[i] = sigmaFactors3D[i]; |
---|
| 208 | } |
---|
| 209 | //----------------------------------------------------------------------- |
---|
[3] | 210 | |
---|
[378] | 211 | void Filter::loadHaar() |
---|
| 212 | { |
---|
| 213 | double filter[3] = {0., 1./2., 1./2.}; |
---|
| 214 | this->name = "Haar wavelet"; |
---|
| 215 | this->filter1D.resize(3); |
---|
| 216 | for(int i=0;i<3;i++) this->filter1D[i] = filter[i]; |
---|
| 217 | this->sigmaFactors.resize(3); |
---|
| 218 | this->maxNumScales.resize(3); |
---|
[3] | 219 | |
---|
[378] | 220 | this->maxNumScales[0] = 6; |
---|
| 221 | double sigmaFactors1D[7] = {1.00000000000,7.07167810e-1,5.00000000e-1,3.53553391e-1, |
---|
| 222 | 2.50000000e-1,1.76776695e-1,1.25000000e-1}; |
---|
| 223 | this->sigmaFactors[0]->resize(7); |
---|
| 224 | for(int i=0;i<19;i++)(*this->sigmaFactors[0])[i] = sigmaFactors1D[i]; |
---|
[3] | 225 | |
---|
[378] | 226 | this->maxNumScales[1] = 6; |
---|
| 227 | double sigmaFactors2D[7] = {1.00000000000,4.33012702e-1,2.16506351e-1,1.08253175e-1, |
---|
| 228 | 5.41265877e-2,2.70632939e-2,1.35316469e-2}; |
---|
| 229 | this->sigmaFactors[1]->resize(7); |
---|
| 230 | for(int i=0;i<12;i++)(*this->sigmaFactors[1])[i] = sigmaFactors2D[i]; |
---|
[3] | 231 | |
---|
[378] | 232 | this->maxNumScales[2] = 8; |
---|
| 233 | double sigmaFactors3D[9] = {1.00000000000,9.35414347e-1,3.30718914e-1,1.16926793e-1, |
---|
| 234 | 4.13398642e-2,1.46158492e-2,5.16748303e-3,1.82698115e-3, |
---|
| 235 | 6.45935379e-4}; |
---|
| 236 | this->sigmaFactors[2]->resize(9); |
---|
| 237 | for(int i=0;i<12;i++)(*this->sigmaFactors[2])[i] = sigmaFactors3D[i]; |
---|
| 238 | } |
---|
| 239 | |
---|
[3] | 240 | } |
---|