| TerraLib 4.1

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files
- Examples

- File List
- File Members

E:/Projetos_Primeiro_Semestre_2012/TerraView/terralib/src/terralib/stat/TeKernelFunctions.h
Go to the documentation of this file.

    00001 /************************************************************************************
    00002 Exploring and analysis of geographical data using TerraLib and TerraView
    00003 Copyright � 2003-2007 INPE and LESTE/UFMG.
    00004
    00005 Partially funded by CNPq - Project SAUDAVEL, under grant no. 552044/2002-4,
    00006 SENASP-MJ and INPE
    00007
    00008 This program is free software; you can redistribute it and/or
    00009 modify it under the terms of the GNU Lesser General Public
    00010 License as published by the Free Software Foundation; either
    00011 version 2.1 of the License, or (at your option) any later version.
    00012
    00013 This program is distributed hoping it will be useful, however, WITHOUT ANY
    00014 WARRANTIES; neither to the implicit warranty of MERCHANTABILITY OR SUITABILITY FOR
    00015 AN SPECIFIC FINALITY. Consult the GNU General Public License for more details.
    00016
    00017 You must have received a copy of the GNU General Public License with this program.
    00018 In negative case, write to the Free Software Foundation, Inc. in the following
    00019 address: 59 Temple Street, Suite 330, Boston, MA 02111-1307 USA.
    00020 ***********************************************************************************/
    00025 #ifndef __TERRALIB_INTERNAL_KERNELFUNCTIONS_H
    00026 #define __TERRALIB_INTERNAL_KERNELFUNCTIONS_H
    00027
    00028 #include "../kernel/TeException.h"
    00029 #include "TeKernelParams.h"
    00030 #include "TeStatDefines.h"
    00031 #include "../kernel/TeGeometryAlgorithms.h"
    00032 #include "../kernel/TeProgress.h"
    00033
    00034
    00035 #define   INVALID_KERNEL -TeMAXFLOAT
    00036
    00037 #define IDX_ATTR            0
    00038 #define IDX_KERNEL          0
    00039 #define IDX_NUMERATOR       1
    00040 #define IDX_KERNEL1         0
    00041 #define IDX_KERNEL2         1
    00042
    00043 enum TeKernelError {
    00044  KERNEL_NO_CENTROID,
    00045    KERNEL_ALL_NULL,
    00046    KERNEL_INVALID_AREA,
    00047    KERNEL_INVALID_VALUE,
    00048    KERNEL_NO_VALUE
    00049 };
    00050
    00051 class STAT_DLL TeKernelException  {
    00052
    00053   TeKernelError error_;
    00054
    00055 public:
    00056
    00057   TeKernelException(TeKernelError e) {
    00058     error_ = e;
    00059   }
    00060
    00061   TeKernelError getErrorCode() {
    00062     return error_;
    00063   }
    00064 };
    00065
    00066
    00067
    00074 STAT_DLL double TeKernelQuartic(double tau, double distance, double intensity);
    00075 STAT_DLL double TeKernelNormal(double tau, double distance, double intensity);
    00076 STAT_DLL double TeKernelUniform(double tau, double distance, double intensity);
    00077 STAT_DLL double TeKernelTriangular(double tau, double distance, double intensity);
    00078 STAT_DLL double TeKernelNegExponential(double tau, double distance, double intensity);
    00079
    00080
    00086 template <typename ItReg>
    00087 double TeKernelGeometricMean(ItReg& beginReg, ItReg& endReg, int idx)
    00088 {
    00089    double mantissa, MediaE, MediaM;
    00090    int expoente, Cont;
    00091    int MediaETmp = 0;
    00092    double MediaMTmp = 0;
    00093
    00094    double log_two = log(2.0);
    00095    Cont = 0;
    00096    ItReg it = beginReg;
    00097    double val;
    00098    while (it != endReg) {
    00099          if ((*it).getDoubleProperty(idx, val) && (val > 0))
    00100          {
    00101             mantissa = frexp(val,&expoente);
    00102             MediaMTmp = MediaMTmp+log(mantissa);
    00103             MediaETmp = MediaETmp+expoente;
    00104             Cont++;
    00105          }
    00106          ++it;
    00107    }
    00108
    00109    MediaE = MediaETmp;
    00110    MediaM = (MediaMTmp+(MediaE*log_two))/Cont;
    00111    MediaM = exp(MediaM);
    00112
    00113    return MediaM;
    00114 }
    00115
    00116
    00117
    00123 template <typename ItEvt>
    00124 double KernelNormalizationFactor(ItEvt itBegin, ItEvt itEnd) {
    00125
    00126   double normFactor = 0;
    00127   double intensity;
    00128
    00129   ItEvt evtIter = itBegin;
    00130
    00131    while (evtIter != itEnd) {
    00132
    00133      (*evtIter).getDoubleProperty(IDX_ATTR, intensity);
    00134      normFactor += intensity;
    00135       ++evtIter;
    00136     }
    00137
    00138   return normFactor;
    00139
    00140 }
    00141
    00142
    00153 template <typename ItEvt>
    00154 double TeKernelValue (TeCoord2D& center, ItEvt& begin, ItEvt& end,
    00155                                                          double radius, TeKernelFunctionType kfunc)
    00156 {
    00157
    00158   string    intstr;
    00159
    00160            TeCoord2D location;              //Event spatial location
    00161   double    intensity;           //Event intensity
    00162
    00163   //Final kernel value
    00164            double kernel = 0;
    00165
    00166   //Iterates through events, getting their location and intensity values
    00167   //and updating kernel value
    00168            ItEvt it = begin;
    00169            while (it != end)
    00170            {
    00171                       (*it).centroid(location);
    00172
    00173                       //If there is no properties, assume intensity of one
    00174                if (!(*it).getDoubleProperty(IDX_ATTR, intensity)) {
    00175                                             intensity = 1.0;
    00176                       }
    00177
    00178                       //Evaluates kernel value for one event
    00179                       double distance = TeDistance (location,center);
    00180                       double localK;
    00181
    00182                       switch(kfunc) {
    00183                       case TeKQuartic:
    00184                                  localK = TeKernelQuartic(radius,distance,intensity);
    00185                                  break;
    00186                       case TeKNormal:
    00187                                  localK = TeKernelNormal(radius,distance,intensity);
    00188                                  break;
    00189                       case TeKTriangular:
    00190                                  localK = TeKernelTriangular(radius,distance,intensity);
    00191                                  break;
    00192                       case TeKNegExponential:
    00193                                  localK = TeKernelNegExponential(radius,distance,intensity);
    00194                                  break;
    00195                       case TeKUniform:
    00196                                  localK = TeKernelUniform(radius,distance,intensity);
    00197                                  break;
    00198                       }
    00199                       kernel += localK;
    00200                ++it;
    00201
    00202            }
    00203            return kernel;
    00204 }
    00205
    00215 template<typename EventSet, typename ItRegionSet>
    00216 bool TeKernelNormalize(EventSet& events,
    00217                    ItRegionSet& regionsBegin,
    00218                    ItRegionSet& regionsEnd,
    00219                    TeKernelComputeType  kType,
    00220                    double totKernel,
    00221                    int idxProp)
    00222 {
    00223
    00224   double normKernel, kernel, area;
    00225   try {
    00226
    00227   //Evaluates sum of attributes for normalization
    00228     double normFactor =
    00229       KernelNormalizationFactor(events.begin(), events.end());
    00230
    00231     //Reiterates through regions defining normalized kernel value,
    00232     //given the computation type
    00233     ItRegionSet regIter = regionsBegin;
    00234     while (regIter != regionsEnd)
    00235     {
    00236
    00237       (*regIter).getDoubleProperty(idxProp, kernel);
    00238       if (kernel < 0) {
    00239          throw TeKernelException(KERNEL_INVALID_VALUE);
    00240       }
    00241
    00242       //Evaluates value based on computation type
    00243       switch(kType) {
    00244       case TeKMovingAverage:
    00245         normKernel = (kernel * normFactor)/totKernel;
    00246         break;
    00247       case TeKDensity:
    00248         //If region does not have an area, assumes one
    00249         (*regIter).area(area);
    00250         if (area <= 0) {
    00251           throw TeKernelException(KERNEL_INVALID_AREA);
    00252         }
    00253         normKernel = ((kernel * normFactor)/totKernel)/area;
    00254         break;
    00255       case TeKProbability:
    00256         normKernel = kernel/totKernel;
    00257         break;
    00258       }
    00259
    00260       //Store final value
    00261       (*regIter).setDoubleProperty(idxProp, normKernel);
    00262       ++regIter;
    00263     }
    00264   } catch (TeKernelException /* e */) {
    00265 //    int a = e.getErrorCode();
    00266     return false;
    00267   }
    00268   return true;
    00269 }
    00270
    00271
    00282 template<typename EventSet, typename ItRegionSet>
    00283 bool TeStatIntKernel(EventSet& events,
    00284                    ItRegionSet& regionsBegin,
    00285                    ItRegionSet& regionsEnd,
    00286                    TeKernelFunctionType kfunc,
    00287                    TeKernelComputeType  kType,
    00288                    double radius,
    00289                    int idxProp)
    00290 {
    00291
    00292
    00293   //Support region attributes
    00294   TeCoord2D location;
    00295   double    kernel;
    00296
    00297   //Regions iterator
    00298   ItRegionSet regIter = regionsBegin;
    00299   // Sum of kernels
    00300   double totKernel = 0;
    00301
    00302   //Treats unexpected expections returning false
    00303   try {
    00304
    00305            int dt = CLOCKS_PER_SEC/4;
    00306            int dt2 = CLOCKS_PER_SEC * 5;
    00307            clock_t    t0, t1, t2;
    00308            t2 = clock();
    00309            t0 = t1 = t2;
    00310            int nsteps = 0;
    00311     //Iterates through regions
    00312     while (regIter != regionsEnd)
    00313     {
    00314       (*regIter).centroid(location);
    00315
    00316       //Define iterators for event sets -- may optimize
    00317       typename EventSet::iterator subSetBegin = events.begin(location,radius);
    00318       typename EventSet::iterator subSetEnd = events.end(location,radius);
    00319       if (kfunc == TeKNormal) {
    00320         subSetBegin = events.begin();
    00321         subSetEnd = events.end();
    00322       }
    00323
    00324       //Evaluates kernel value for region
    00325       kernel = TeKernelValue(location, subSetBegin, subSetEnd,
    00326         radius, kfunc);
    00327
    00328       totKernel += kernel;
    00329
    00330       (*regIter).setDoubleProperty(idxProp, kernel);
    00331       ++regIter;
    00332              t2 = clock();
    00333              ++nsteps;
    00334              if (TeProgress::instance())
    00335              {
    00336                        if (int(t2-t1) > dt)
    00337                        {
    00338                                  t1 = t2;
    00339                                  if (TeProgress::instance()->wasCancelled())
    00340                                             break;
    00341                                  if((int)(t2-t0) > dt2)
    00342                                             TeProgress::instance()->setProgress(nsteps);
    00343                        }
    00344              }
    00345            }
    00346            if (TeProgress::instance())
    00347                       TeProgress::instance()->reset();
    00348
    00349     if (totKernel <= 0) {
    00350       throw TeKernelException(KERNEL_ALL_NULL);
    00351     }
    00352
    00353     return TeKernelNormalize(events, regionsBegin, regionsEnd, kType, totKernel,
    00354       idxProp);
    00355
    00356   } catch (TeKernelException /* e */) {
    00357 //    int a = e.getErrorCode();
    00358     return false;
    00359   }
    00360
    00361   return true;
    00362 }
    00363
    00364
    00376 template<typename EventSet, typename ItRegionSet>
    00377 bool TeStatAdaptiveGeoMeanIntKernel(EventSet& events,
    00378                                     ItRegionSet& regionsBegin,
    00379                                     ItRegionSet& regionsEnd,
    00380                                     TeKernelFunctionType kfunc,
    00381                                     TeKernelComputeType ktype,
    00382                                     int /* numReg */,
    00383                                     double totalArea,
    00384                                     int idxProp)
    00385 {
    00386
    00387   try {
    00388
    00389     //Evaluate kernel with a fixed radius, based on formula, storing in vecValues
    00390     double radius = 0.68*pow((double)events.numObjects(),-0.2)*sqrt(totalArea);
    00391     double sqArea = sqrt(totalArea);
    00392     if (!TeStatIntKernel(events, regionsBegin, regionsEnd,
    00393       kfunc, ktype, radius, idxProp))
    00394       return false;
    00395
    00396     //Evaluate geometric mean of kernel values, to adjust radius
    00397     double meanKernel = TeKernelGeometricMean(regionsBegin, regionsEnd, idxProp);
    00398
    00399     if (meanKernel <= 0) {
    00400       throw TeKernelException(TeKernelException(KERNEL_ALL_NULL));
    00401     }
    00402
    00403     //Now, reassign radius, evaluating final value for kernel
    00404
    00405     ItRegionSet regIter = regionsBegin;
    00406     double    prevKernel;
    00407     double    newKernel;
    00408     double    newRadius;
    00409     double    totKernel = 0;
    00410     TeCoord2D location;
    00411
    00412            int dt = CLOCKS_PER_SEC/4;
    00413            int dt2 = CLOCKS_PER_SEC * 5;
    00414            clock_t    t0, t1, t2;
    00415            t2 = clock();
    00416            t0 = t1 = t2;
    00417            int nsteps = 0;
    00418
    00419     //Iterates through regions
    00420     while (regIter != regionsEnd)
    00421     {
    00422       (*regIter).centroid(location);
    00423
    00424       //Define iterators for event sets -- may optimize
    00425       typename EventSet::iterator subSetBegin = events.begin(location,radius);
    00426       typename EventSet::iterator subSetEnd = events.end(location,radius);
    00427       if (kfunc == TeKNormal) {
    00428         subSetBegin = events.begin();
    00429         subSetEnd = events.end();
    00430       }
    00431
    00432       if (!((*regIter).getDoubleProperty(idxProp, prevKernel)) || (prevKernel < 0)) {
    00433          throw TeKernelException(KERNEL_INVALID_VALUE);
    00434       }
    00435
    00436       if (prevKernel > 0) {
    00437         newRadius = radius*pow((meanKernel/prevKernel),0.5);
    00438         if (newRadius > sqArea/4.0)
    00439           newRadius = sqArea/4.0; //Para funcionar quando o valor eh muito pequeno
    00440         newKernel = TeKernelValue(location, subSetBegin, subSetEnd,  newRadius, kfunc);
    00441       }
    00442       else {
    00443         newKernel = 0.0;
    00444       }
    00445
    00446       totKernel += newKernel;
    00447       (*regIter).setDoubleProperty(idxProp, newKernel);
    00448       ++regIter;
    00449              t2 = clock();
    00450              ++nsteps;
    00451              if (TeProgress::instance())
    00452              {
    00453                        if (int(t2-t1) > dt)
    00454                        {
    00455                                  t1 = t2;
    00456                                  if (TeProgress::instance()->wasCancelled())
    00457                                             break;
    00458                                  if((int)(t2-t0) > dt2)
    00459                                             TeProgress::instance()->setProgress(nsteps);
    00460                        }
    00461              }
    00462     }
    00463            if (TeProgress::instance())
    00464                       TeProgress::instance()->reset();
    00465
    00466     if (totKernel <= 0) {
    00467       throw TeKernelException(KERNEL_ALL_NULL);
    00468     }
    00469
    00470     return TeKernelNormalize(events, regionsBegin, regionsEnd, ktype, totKernel,
    00471       idxProp);
    00472   } catch (TeKernelException /* e */) {
    00473 //      int a = e.getErrorCode();
    00474       return false;
    00475   }
    00476
    00477   return true;
    00478 }
    00479
    00480
    00481
    00485 class STAT_DLL TeStatKernel {
    00486
    00487 public:
    00488
    00490   TeKernelFunctionType kfunc_;
    00492   TeKernelComputeType  ktype_;
    00494   double               radius_;
    00496   double               totalArea_;
    00498   int                  numReg_;
    00499
    00500   TeStatKernel() {
    00501     kfunc_ = TeKQuartic;
    00502     ktype_ = TeKDensity;
    00503     radius_ = 0;
    00504     totalArea_ = 0;
    00505     numReg_ = 0;
    00506   }
    00507
    00515   template<typename EventSet, typename ItRegionSet>
    00516   bool apply(EventSet& events,
    00517     ItRegionSet regionsBegin,
    00518     ItRegionSet regionsEnd) {
    00519     if (totalArea_ <= 0) {
    00520       return false;
    00521     }
    00522
    00523     if (numReg_ <= 0) {
    00524       return false;
    00525     }
    00526
    00527     bool result;
    00528     //Radius defined, kernel with a fixed radius
    00529     if (radius_ > 0) {
    00530       result = TeStatIntKernel(events,regionsBegin, regionsEnd,
    00531         kfunc_,ktype_, radius_, IDX_KERNEL);
    00532
    00533     }
    00534     //Radius undefined, adaptive kernel
    00535     else {
    00536       result = TeStatAdaptiveGeoMeanIntKernel(events,  regionsBegin, regionsEnd,
    00537         kfunc_, ktype_, numReg_, totalArea_, IDX_KERNEL);
    00538
    00539     }
    00540     return result;
    00541   }
    00542
    00543 };
    00544
    00545
    00546
    00547 class STAT_DLL TeStatKernelRatio
    00548 {
    00549   public:
    00550   int                                                             numReg_;
    00551
    00552   public:
    00553
    00554   TeKernelCombinationType        kComb_;
    00555   TeKernelComputeType              ktype_;
    00556   TeKernelFunctionType                      kfunc1_;
    00557   double                                                       radius1_;
    00558   TeKernelFunctionType                      kfunc2_;
    00559   double                                               radius2_;
    00560   double                                               totalArea_;
    00561
    00562
    00563   TeStatKernelRatio() {
    00564     kComb_ = TeKRatio;
    00565     kfunc1_ =  kfunc2_ = TeKQuartic;
    00566     ktype_ = TeKDensity;
    00567     radius1_ = radius2_ = 0;
    00568     totalArea_ = 0;
    00569     numReg_ = 0;
    00570   }
    00571
    00572   template <typename EventSet1, typename EventSet2, typename ItRegionSet>
    00573   bool apply (EventSet1& ev1,
    00574     EventSet2& ev2,
    00575     ItRegionSet regBegin,
    00576     ItRegionSet regEnd) {
    00577
    00578
    00579     if (totalArea_ <= 0)
    00580       return false;
    00581
    00582     if (numReg_ <= 0)
    00583       return false;
    00584
    00585     //Calcula kernel para primeiro conjunto
    00586
    00587     bool result;
    00588     //Radius defined, kernel with a fixed radius
    00589     if (radius1_ > 0) {
    00590       result = TeStatIntKernel(ev1,regBegin, regEnd,
    00591         kfunc1_,ktype_, radius1_, IDX_KERNEL1);
    00592
    00593     }
    00594     //Radius undefined, adaptive kernel
    00595     else {
    00596       result = TeStatAdaptiveGeoMeanIntKernel(ev1,  regBegin, regEnd,
    00597         kfunc1_, ktype_, numReg_, totalArea_, IDX_KERNEL1);
    00598
    00599     }
    00600
    00601     if (!result)
    00602       return result;
    00603
    00604
    00605     //Calcula kernel para segundo conjunto
    00606      //Radius defined, kernel with a fixed radius
    00607     if (radius2_ > 0) {
    00608       result = TeStatIntKernel(ev2,regBegin, regEnd,
    00609         kfunc2_,ktype_, radius2_, IDX_KERNEL2);
    00610
    00611     }
    00612     //Radius undefined, adaptive kernel
    00613     else {
    00614       result = TeStatAdaptiveGeoMeanIntKernel(ev2,  regBegin, regEnd,
    00615         kfunc2_, ktype_, numReg_, totalArea_, IDX_KERNEL2);
    00616
    00617     }
    00618
    00619     if (!result)
    00620       return result;
    00621
    00622     //Itera calculando o tipo de combinacao e armazenando no valor da regiao
    00623     //Percorre todas as regioes
    00624 //    int Cont = 0;
    00625     double kernel;
    00626     double k1;
    00627     double k2;
    00628     double area;
    00629
    00630     ItRegionSet regIter = regBegin;
    00631
    00632     while (regIter != regEnd)
    00633     {
    00634
    00635       (*regIter).area(area);
    00636
    00637       (*regIter).getDoubleProperty(IDX_KERNEL1, k1);
    00638       (*regIter).getDoubleProperty(IDX_KERNEL2, k2);
    00639
    00640       switch(kComb_)
    00641       {
    00642       case TeKRatio:
    00643         if (k2 == 0.0)
    00644           kernel = 0.0;//INVALID_KERNEL;
    00645         else
    00646           kernel = k1/k2;
    00647         break;
    00648       case TeKLogRatio:
    00649         if (k2 == 0.0)
    00650           kernel = INVALID_KERNEL;
    00651         else
    00652           kernel = log(k1/k2);
    00653         break;
    00654       case TeAbsDifference:
    00655         kernel = k1 - k2;
    00656         break;
    00657       case TeRelDifference:
    00658         kernel = (k1 - k2) * area;
    00659         break;
    00660       case TeAbsSum:
    00661         kernel = k1+k2;
    00662         break;
    00663       case TeRelSum:
    00664         kernel = (k1 + k2) * area;
    00665         break;
    00666       }
    00667              (*regIter).setDoubleProperty(IDX_KERNEL, kernel);
    00668              ++regIter;
    00669     }
    00670
    00671     return true;
    00672   }
    00673
    00674 };
    00675
    00676
    00677 #endif
    00678
    00679
    00680
    00681
    00682

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines