0001: //////////////////////////////////////////////////
0002: //
0003: //                      ppa-r
0004: //  point process analysis system: pair correlation
0005: //
0006: //      Copyright  October 2001 by TAKENAKA, A.
0007: //           rev.  November 2004
0008: 
0009: #include <iostream>
0010: 
0011: #include "plotarea.h"
0012: 
0013: //////////////////////////////
0014: //
0015: //  Constructor  # Does nothing.
0016: 
0017: PlotArea::PlotArea(void)
0018: {}
0019: 
0020: //////////////////////////////
0021: //
0022: 
0023: bool PlotArea::loadPoints(const char* file1, const char* file2)
0024: {
0025:     Points1.erase(Points1.begin(), Points1.end());
0026:     Points2.erase(Points2.begin(), Points2.end());
0027: 
0028:     Points1.load(file1);
0029:     if (Points1.size() == 0) {  //  No data.
0030:         return false;
0031:     }
0032: 
0033:     if (file2 != NULL) {  // If two files are specified...
0034:         Points2.load(file2);
0035:         if (Points2.size() == 0) {  //  No data.
0036:             return false;
0037:         }
0038:     }
0039: 
0040:   //  The area of the first data set is stored.
0041: 
0042:     AreaSize  = Points1.getAreaSize();
0043:     Area      = AreaSize.X * AreaSize.Y;
0044:     Perimeter = 2.0 * (AreaSize.X + AreaSize.Y);
0045: 
0046:     Density1  = Points1.size() / Area;
0047:     MeanDensity = Density1;
0048: 
0049:     if (Points2.size() > 0) {
0050:         Density2 = Points2.size() / Area;
0051:         MeanDensity = sqrt(Density1 * Density2);
0052:     }
0053: 
0054:     if ((Points2.size()) > 0 && !(AreaSize == Points2.getAreaSize()) ) {
0055:         std::cerr << "! WARNING: Area of two data sets does not match."
0056:                   << std::endl;
0057:     }
0058: 
0059:     Delta = 0.15 / sqrt(MeanDensity);  //  width of Epanechnikov kernel.
0060: 
0061:     std::cout << "Delta\t" << Delta << std::endl;
0062: 
0063:     return true;
0064: }
0065: 
0066: //////////////////////////////
0067: //
0068: 
0069: bool  PlotArea::calc_cor(double r_unit, double r_upto)
0070: {
0071:     if (r_unit <= 0 || r_unit > r_upto) {  //  Invalid r specification.
0072:         std::cerr << "Invalid r specification." << std::endl;
0073:         return false;
0074:     }
0075: 
0076:     R_Unit = r_unit;
0077:     R_Upto = r_upto;
0078: 
0079:     //  R is kept shorter than the width of the area.
0080: 
0081:     double min_side = (AreaSize.X < AreaSize.Y) ? AreaSize.X : AreaSize.Y;
0082:     if (R_Upto > min_side) {
0083:         R_Upto = min_side;
0084:     }
0085: 
0086:     //  r = R_Unit, 2 * R_Unit, 3 * R_Unit, ...(MaxR_Index - 1) * R_Unit.
0087: 
0088:     MaxR_Index = int (R_Upto / R_Unit);
0089: 
0090:     initArray();  //  Initialize array to store the results.
0091: 
0092:     PointList::iterator itr1, itr2;
0093:     PointList::iterator itr2_begin, itr2_end;
0094: 
0095:     if (Points2.size() > 0) {  //  between two sets of points.
0096:         itr2_begin = Points2.begin();
0097:         itr2_end   = Points2.end();
0098:     }
0099:     else {                     //  within a set of points.
0100:         itr2_begin = Points1.begin();
0101:         itr2_end   = Points1.end();
0102:     }
0103: 
0104:     for (itr1 = Points1.begin(); itr1 != Points1.end(); ++itr1) {
0105:         for (itr2 = itr2_begin; itr2 != itr2_end; ++itr2) {
0106: 
0107:             if (itr1 == itr2) {
0108:                 continue;
0109:             }
0110: 
0111:             double r = (*itr1).distance(*itr2);
0112: 
0113:             for (int i_r = 0; i_r < MaxR_Index; ++i_r) {
0114: 
0115:                 double d = (i_r + 0.5) * R_Unit;
0116:                 double w = get_w(r - d);  //  Epanechnikov kernel
0117: 
0118:                 if (w > 0) {
0119:                     double s = get_s(d);  // Ohser's edge correction factor.
0120:                     ResultsArray[i_r].second += w / s / d;
0121:                 }
0122:             }
0123:         }
0124:     }
0125: 
0126:     //  Normalize for point density.
0127: 
0128:     for (int i = 0; i < MaxR_Index; ++i) {
0129:         ResultsArray[i].second /= (MeanDensity * MeanDensity * 2.0 * M_PI);
0130:     }
0131: 
0132:     return true;
0133: }
0134: 
0135: //////////////////////////////
0136: //
0137: 
0138: const r_X_Array& PlotArea::getResults(void)
0139: {
0140:     return ResultsArray;
0141: }
0142: 
0143: //////////////////////////////
0144: //
0145: //  Initialize the result array.
0146: 
0147: void PlotArea::initArray(void)
0148: {
0149:     ResultsArray.erase(ResultsArray.begin(), ResultsArray.end());
0150: 
0151:     for (int i = 0; i < MaxR_Index; ++i) {
0152:         double r = (i + 0.5) * R_Unit;
0153:         ResultsArray.push_back(r_X_Pair(r, 0.0));
0154:     }
0155: }
0156: 
0157: 
0158: //////////////////////////////
0159: //
0160: //   Ohser's edge correction
0161: 
0162: double PlotArea::get_s(double x)
0163: {
0164:     return Area - x * (Perimeter - x) / M_PI;
0165: }
0166: 
0167: //////////////////////////////
0168: //
0169: //  Epanechnikov kernel
0170: 
0171: double PlotArea::get_w(double d)
0172: {
0173:     if (d < 0) {
0174:         d = -d;
0175:     }
0176: 
0177:     if (d >= Delta) {
0178:         return 0;
0179:     }
0180:     else {
0181:         return  0.75 / Delta * (1.0 - d * d / (Delta * Delta));
0182:     }
0183: }