| [932] | 1 | // This may look like C code, but it is really -*- C++ -*- | 
|---|
|  | 2 | //-------------------------------------------------------------------------- | 
|---|
|  | 3 | // File and Version Information: | 
|---|
|  | 4 | //        Simple Operations on Maps | 
|---|
|  | 5 | // Description: | 
|---|
|  | 6 | // | 
|---|
|  | 7 | // History (add to end): | 
|---|
|  | 8 | //      Sophie   April, 2000  - creation | 
|---|
|  | 9 | //        code from Reza | 
|---|
|  | 10 | //------------------------------------------------------------------------ | 
|---|
|  | 11 | #include "machdefs.h" | 
|---|
|  | 12 | #include <math.h> | 
|---|
|  | 13 | #include "mapoperation.h" | 
|---|
|  | 14 | /*! | 
|---|
| [1371] | 15 | * \class SOPHYA::MapOperations | 
|---|
| [932] | 16 | * This class is for simple operations on maps | 
|---|
|  | 17 | */ | 
|---|
|  | 18 |  | 
|---|
|  | 19 | template <class T> | 
|---|
|  | 20 | void MapOperations<T>::sphTosph(PixelMap<T>& in, PixelMap<T>& out) | 
|---|
|  | 21 | // Cette fonction remplit la sphere out a partir de la sphere in | 
|---|
|  | 22 | // Les spheres peuvent etre de type et de pixelisations differentes | 
|---|
|  | 23 | { | 
|---|
|  | 24 | int kin,kout; | 
|---|
|  | 25 | double teta,phi; | 
|---|
|  | 26 |  | 
|---|
|  | 27 | int* cnt = new int[out.NbPixels()+1]; | 
|---|
|  | 28 | for(kout=0; kout<out.NbPixels(); kout++) | 
|---|
|  | 29 | { cnt[kout] = 0; out(kout) = 0.; } | 
|---|
|  | 30 |  | 
|---|
|  | 31 | for(kin=0; kin<in.NbPixels(); kin++) { | 
|---|
|  | 32 | in.PixThetaPhi(kin, teta, phi); | 
|---|
|  | 33 | kout = out.PixIndexSph(teta, phi); | 
|---|
| [1950] | 34 | if (kout < 0) continue; | 
|---|
| [932] | 35 | out(kout) += in(kin); | 
|---|
|  | 36 | cnt[kout] ++; | 
|---|
|  | 37 | } | 
|---|
|  | 38 |  | 
|---|
|  | 39 | double moy, sig, dcn; | 
|---|
|  | 40 | moy = 0.; sig = 0.; | 
|---|
|  | 41 | for(kout=0; kout<out.NbPixels(); kout++) { | 
|---|
|  | 42 | dcn = cnt[kout];  moy += dcn;        sig += (dcn*dcn); | 
|---|
|  | 43 | if (cnt[kout] > 0)  out(kout) /= dcn; | 
|---|
|  | 44 | else { | 
|---|
|  | 45 | out.PixThetaPhi(kout, teta, phi); | 
|---|
|  | 46 | int pixel = in.PixIndexSph(teta,phi); | 
|---|
|  | 47 | out(kout) = in.PixVal(pixel); | 
|---|
|  | 48 | } | 
|---|
|  | 49 | } | 
|---|
|  | 50 |  | 
|---|
|  | 51 | moy /= out.NbPixels(); | 
|---|
|  | 52 | sig = sig/out.NbPixels() - moy*moy; | 
|---|
|  | 53 | if (sig >= 0.)        sig = sqrt(sig); | 
|---|
|  | 54 |  | 
|---|
|  | 55 | delete[] cnt; | 
|---|
|  | 56 |  | 
|---|
|  | 57 | } | 
|---|
|  | 58 |  | 
|---|
|  | 59 | template <class T> | 
|---|
|  | 60 | void  MapOperations<T>::meanSig(NDataBlock<T> const & dbl, double& gmoy, double& gsig) | 
|---|
|  | 61 | { | 
|---|
|  | 62 | gmoy=0.; | 
|---|
|  | 63 | gsig = 0.; | 
|---|
|  | 64 | double valok; | 
|---|
|  | 65 | for(int k=0; k<(int)dbl.Size(); k++) { | 
|---|
|  | 66 | valok = dbl(k); | 
|---|
|  | 67 | gmoy += valok;  gsig += valok*valok; | 
|---|
|  | 68 | } | 
|---|
|  | 69 | gmoy /= (double)dbl.Size(); | 
|---|
|  | 70 | gsig = gsig/(double)dbl.Size() - gmoy*gmoy; | 
|---|
|  | 71 | if (gsig >= 0.) gsig = sqrt(gsig); | 
|---|
|  | 72 | } | 
|---|
|  | 73 |  | 
|---|
|  | 74 | /////////////////////////////////////////////////////////////// | 
|---|
|  | 75 | #ifdef __CXX_PRAGMA_TEMPLATES__ | 
|---|
|  | 76 | #pragma define_template MapOperations<r_4> | 
|---|
|  | 77 | #pragma define_template MapOperations<r_8> | 
|---|
|  | 78 | #endif | 
|---|
|  | 79 |  | 
|---|
|  | 80 | #if defined(ANSI_TEMPLATES) || defined(GNU_TEMPLATES) | 
|---|
|  | 81 | template class MapOperations<r_4>; | 
|---|
| [940] | 82 | template class MapOperations<r_8>; | 
|---|
| [932] | 83 | #endif | 
|---|
|  | 84 |  | 
|---|