plotarea.cpp

//////////////////////////////////////////////////
//
//  This file is to be used for point process analysis
//
//      Copyright  October 2001 by TAKENAKA, A.

#include <iostream>

#include "plotarea.h"

#define M_PI 3.14159265358979323846

//////////////////////////////
//
//  Constructor  # Does nothing.

PlotArea::PlotArea(void)
{}

//////////////////////////////
//

bool PlotArea::loadPoints(const char* file1, const char* file2)
{
    Points1.erase(Points1.begin(), Points1.end());
    Points2.erase(Points2.begin(), Points2.end());

    Points1.load(file1);
    if (Points1.size() == 0) {  //  No data.
        return false;
    }

    if (file2 != NULL) {  // If two files are specified...
        Points2.load(file2);
        if (Points2.size() == 0) {  //  No data.
            return false;
        }
    }

  //  The area of the first data set is stored.

    AreaSize  = Points1.getAreaSize();
    Area      = AreaSize.X * AreaSize.Y;
    Perimeter = 2.0 * (AreaSize.X + AreaSize.Y);

    Density1  = Points1.size() / Area;
    if (Points2.size() > 0) {
        Density2 = Points2.size() / Area;
    }

    if ((Points2.size()) > 0 && !(AreaSize == Points2.getAreaSize()) ) {
        std::cerr << "! WARNING: Area of two data sets does not match."
                  << std::endl;
    }

    return true;
}

//////////////////////////////
//
//  Calculate L function.
//
//  r_unit : increment of the distance, r.
//  r_upto : the furthest distance to calculate L.

bool PlotArea::calc_L(double r_unit, double r_upto)
{
    if (r_unit <= 0 || r_unit > r_upto) {  //  Invalid r specification.
        std::cerr << "Invalid r specification." << std::endl;
        return false;
    }

    R_Unit = r_unit;
    R_Upto = r_upto;

    //  R is kept shorter than the width of the area.

    double min_side = (AreaSize.X < AreaSize.Y) ? AreaSize.X : AreaSize.Y;
    if (R_Upto > min_side) {
        R_Upto = min_side;
    }

    //  r = R_Unit, 2 * R_Unit, 3 * R_Unit, ...(MaxR_Index - 1) * R_Unit.

    MaxR_Index = R_Upto / R_Unit;

    initArray();  //  Initialize array to store the results.
    calc_K();
    fromKtoL();

    return true;
}

//////////////////////////////
//
//  Calculate K function.

void   PlotArea::calc_K(void)
{
    PointList::iterator itr1, itr2;
    PointList::iterator itr2_begin, itr2_end;

    if (Points2.size() > 0) {  //  cross K between two sets of points.
        itr2_begin = Points2.begin();
        itr2_end   = Points2.end();
    }
    else {                     //  K within a set of points.
        itr2_begin = Points1.begin();
        itr2_end   = Points1.end();
    }

    for (itr1 = Points1.begin(); itr1 != Points1.end(); ++itr1) {
        for (itr2 = itr2_begin; itr2 != itr2_end; ++itr2) {

            if (itr1 == itr2) {  //  Does not happen in case of cross K.
                continue;
            }

            double r = (*itr1).distance(*itr2);
            int    r_index = r_to_index(r);

            if (r_index < 0) {
                continue;  //  r_to_index(r) returns -1 if points are too far.
            }

  //  if the point (*itr2) is within a dinstance n *R_Unit from the
  //  point (*itr1), the point (*itr2) is also within a distance
  //  (n+1) * R_Unit, (n+2) * R_Unit ,... (MaxR_Index - 1) * R_Unit.

            double s = get_s(r);  // Ohser's edge correction factor.

            for (int i = r_index; i < MaxR_Index; ++i) {
                ResultsArray[i].second += 1.0 / s;
            }
        }
    }

    //  Normalize for point density.

    double dd = (Points2.size() == 0) ? (Density1 * Density1)
                                     : (Density1 * Density2);

    for (int i = 0; i < MaxR_Index; ++i) {
        ResultsArray[i].second /= dd;
    }
}

//////////////////////////////
//

const r_X_Array& PlotArea::getResults(void)
{
    return ResultsArray;
}

//////////////////////////////
//
//  Initialize the result array.

void PlotArea::initArray(void)
{
    ResultsArray.erase(ResultsArray.begin(), ResultsArray.end());

    for (int i = 0; i < MaxR_Index; ++i) {
        double r = (i + 1) * R_Unit;
        ResultsArray.push_back(r_X_Pair(r, 0.0));
    }
}

//////////////////////////////
//
//  Returns index of the element in the result array (ResultsArray)
//  corresponding to the given distance, r.
//
//  For too large r (index > MaxR_Index), -1 will be returned.


int PlotArea::r_to_index(double r)
{
    int index = (r / R_Unit);

  // If r = n * R_Unit (n is integer), the index should be n-1, not n.
    if (index * R_Unit == r) {
        index--;
    }

    if (index >= MaxR_Index) {   //  Too far.
        return -1;
    }

    return index;
}

//////////////////////////////
//
//   Ohser's edge correction

double PlotArea::get_s(double x)
{
    return Area - x * (Perimeter - x) / M_PI;
}

//////////////////////////////
//
//  Replace K function values in the result array by L values.

void PlotArea::fromKtoL(void)
{
    double r, k, l;

    for (int i = 0; i < MaxR_Index; ++i) {

        r = ResultsArray[i].first;
        k = ResultsArray[i].second;
        l = sqrt(k / M_PI) - r;

        ResultsArray[i].second  = l;  //  replace K by L.
    }
}