00001
00022 #ifndef DGVNOISE_H
00023 #define DGVNOISE_H
00024
00025 #include "DGVAliases.h"
00026
00028 #include <random/discrete-uniform.h>
00029 #include <random/normal.h>
00030 #include <random/exponential.h>
00031
00032 using namespace ranlib;
00033
00035 enum noiseType {deltaFunction = 0, uniformRandom, gaussian, poisson};
00036
00052 template <typename type, int rank>
00053 class DGVNoise
00054 {
00055 public:
00061 DGVNoise();
00065 virtual ~DGVNoise();
00067
00073 inline void setMax(type maxVal)
00074 {
00075 max = maxVal;
00076 }
00080 inline void setSignalToNoiseRatio(type SNR)
00081 {
00082 snr = SNR;
00083 }
00087 inline void setData(Array<type,rank> &data)
00088 {
00089 noise = data;
00090 }
00094 inline void setOffset(coordinate offsetValue)
00095 {
00096 offset = offsetValue;
00097 }
00101 inline void toDeltaFunction()
00102 {
00103 form = deltaFunction;
00104 }
00108 inline void toUniformRandom()
00109 {
00110 form = uniformRandom;
00111 }
00115 inline void toGaussianForm()
00116 {
00117 form = gaussian;
00118 }
00122 inline void toPoissonForm()
00123 {
00124 form = poisson;
00125 }
00127
00133 void generateNoiseFor(Array<type,rank> &data);
00137 void apply(Array<type,rank> &data);
00141 void applyAsIntegers(Array<type,rank> &data);
00145 inline void remove(Array<type,rank> &data)
00146 {
00147 data -= noise;
00148 }
00152 inline Array<type,rank>& data()
00153 {
00154 return noise;
00155 }
00159 void meanErrors(Array<type,rank> &actual, Array<type,rank> &result, Array<type,rank> &errors, type &mse, type &rmse, type &psnr);
00161
00167 type randomNumber_Normal(type mean, type stdev);
00171 long randomNumber_Poisson(type mean);
00175 long randomNumber_Uniform(long max);
00179 type randomNumber_Exponential(type mean);
00181
00182 protected:
00183
00184
00185 Array<type,rank> noise;
00186 type max;
00187 type snr;
00188 coordinate offset;
00189 noiseType form;
00190 };
00191
00192 template <typename type, int rank>
00193 DGVNoise<type,rank>::DGVNoise()
00194 {
00195 max = 1;
00196 snr = 0.1;
00197 offset = coordinate(0,0);
00198 form = deltaFunction;
00199 }
00200
00201 template <typename type, int rank>
00202 DGVNoise<type,rank>::~DGVNoise()
00203 {
00204
00205 }
00206
00207 template <typename type, int rank>
00208 void DGVNoise<type,rank>::generateNoiseFor(Array<type,rank> &data)
00209 {
00210 const int xSize = data.rows(), ySize = data.cols(), zSize = data.depth();;
00211 type noiseLevel = max*(1.0-snr);
00212 int noiseRow = 0, noiseCol = 0, noiseHeight = 0;
00213
00214 noise.resize(data.shape());
00215 noise = 0;
00216 if (rank == 1)
00217 {
00218 switch (form)
00219 {
00220 case deltaFunction:
00221 noise(offset.imag()) += noiseLevel;
00222 break;
00223 case uniformRandom:
00224 for (int j = 0; j < xSize/2; j ++)
00225 {
00226 noiseCol = randomNumber_Uniform(xSize);
00227 noise(noiseCol) = randomNumber_Normal(noiseLevel,sqrt(noiseLevel)) - noiseLevel;
00228 }
00229 break;
00230 case gaussian:
00231 for(int j = 0; j < xSize; j ++)
00232 noise(j) = randomNumber_Normal(noiseLevel,sqrt(noiseLevel)) - noiseLevel;
00233 break;
00234 case poisson:
00235 for(int j = 0; j < xSize; j ++)
00236 noise(j) = randomNumber_Poisson(noiseLevel) - noiseLevel;
00237 break;
00238 }
00239 }
00240 else if (rank == 2)
00241 {
00242 switch (form)
00243 {
00244 case deltaFunction:
00245 noise(offset.real(),offset.imag()) += noiseLevel;
00246 break;
00247 case uniformRandom:
00248 for (int j = 0; j < 2*xSize; j ++)
00249 {
00250 noiseRow = randomNumber_Uniform(xSize);
00251 noiseCol = randomNumber_Uniform(ySize);
00252 noise(noiseRow,noiseCol) = randomNumber_Normal(noiseLevel,sqrt(noiseLevel)) - noiseLevel;
00253 }
00254 break;
00255 case gaussian:
00256 for(int j = 0; j < xSize; j ++)
00257 for(int k = 0; k < ySize; k ++)
00258 noise(j,k) = randomNumber_Normal(noiseLevel,sqrt(noiseLevel)) - noiseLevel;
00259 break;
00260 case poisson:
00261 for(int j = 0; j < xSize; j ++)
00262 for(int k = 0; k < ySize; k ++)
00263 noise(j,k) = randomNumber_Poisson(noiseLevel) - noiseLevel;
00264 break;
00265 }
00266 }
00267 else if (rank == 3)
00268 {
00269 switch (form)
00270 {
00271 case deltaFunction:
00272 noise(offset.real(),offset.imag(),0) += noiseLevel;
00273 break;
00274 case uniformRandom:
00275 for (int j = 0; j < xSize*ySize; j ++)
00276 {
00277 noiseRow = randomNumber_Uniform(xSize);
00278 noiseCol = randomNumber_Uniform(ySize);
00279 noiseHeight = randomNumber_Uniform(zSize);
00280 noise(noiseRow,noiseCol,noiseHeight) = randomNumber_Normal(noiseLevel,sqrt(noiseLevel)) - noiseLevel;
00281 }
00282 break;
00283 case gaussian:
00284 for(int j = 0; j < xSize; j ++)
00285 for(int k = 0; k < ySize; k ++)
00286 for(int l = 0; l < zSize; l ++)
00287 noise(j,k,l) = randomNumber_Normal(noiseLevel,sqrt(noiseLevel)) - noiseLevel;
00288 break;
00289 case poisson:
00290 for(int j = 0; j < xSize; j ++)
00291 for(int k = 0; k < ySize; k ++)
00292 for(int l = 0; l < zSize; l ++)
00293 noise(j,k,l) = randomNumber_Poisson(noiseLevel) - noiseLevel;
00294 break;
00295 }
00296 }
00297 }
00298
00299 template <typename type, int rank>
00300 void DGVNoise<type,rank>::apply(Array<type,rank> &data)
00301 {
00302 int xSize = noise.rows(), ySize = noise.cols();
00303
00304 if(xSize > data.rows())
00305 xSize = data.rows();
00306 if(ySize > data.cols())
00307 ySize = data.cols();
00308
00309 for(int j = 0; j < xSize; j ++)
00310 for(int k = 0; k < ySize; k ++)
00311 data(j,k) += noise(j,k);
00312 }
00313
00314 template <typename type, int rank>
00315 void DGVNoise<type,rank>::applyAsIntegers(Array<type,rank> &data)
00316 {
00317 int xSize = noise.rows(), ySize = noise.cols();
00318
00319 if(xSize > data.rows())
00320 xSize = data.rows();
00321 if(ySize > data.cols())
00322 ySize = data.cols();
00323
00324 for(int j = 0; j < xSize; j ++)
00325 for(int k = 0; k < ySize; k ++)
00326 data(j,k) += static_cast<long>(noise(j,k) + 0.5);
00327 }
00328
00329 template <typename type, int rank>
00330 void DGVNoise<type,rank>::meanErrors(Array<type,rank> &actual, Array<type,rank> &result, Array<type,rank> &errors, type &mse, type &rmse, type &psnr)
00331 {
00332 const int xSize = actual.rows(), ySize = actual.cols();
00333
00334 errors = actual - result;
00335 mse = 0;
00336 for(int j = 0; j < xSize; j ++)
00337 for(int k = 0; k < ySize; k ++)
00338 mse += errors(j,k)*errors(j,k);
00339 mse /= xSize*ySize;
00340 rmse = sqrt(mse);
00341 psnr = 10*log10(255*255/mse);
00342 }
00343
00344 template <typename type, int rank>
00345 type DGVNoise<type,rank>::randomNumber_Normal(type mean, type stdev)
00346 {
00347 Normal<type> valueRNG(mean,stdev);
00348
00349 return valueRNG.random();
00350 }
00351
00352 template <typename type, int rank>
00353 long DGVNoise<type,rank>::randomNumber_Poisson(type mean)
00354 {
00355 Exponential<type> valueRNG(1.0);
00356 double sum = 0.0;
00357 long x = 0;
00358
00359 while (sum < mean)
00360 {
00361 sum += valueRNG.random();
00362 x ++;
00363 }
00364
00365 return x - 1;
00366 }
00367
00368 template <typename type, int rank>
00369 long DGVNoise<type,rank>::randomNumber_Uniform(long max)
00370 {
00371 DiscreteUniform<long> valueRNG(max);
00372
00373 return valueRNG.random()%max;
00374 }
00375
00376 template <typename type, int rank>
00377 type DGVNoise<type,rank>::randomNumber_Exponential(type mean)
00378 {
00379 Exponential<type> valueRNG(mean);
00380
00381 return valueRNG.random();
00382 }
00383
00384 #endif // DGVNOISE_H
