G4KDTree.cc

Go to the documentation of this file.

//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
/*
 * G4KDTree.cc
 *
 *  Created on: 22 oct. 2013
 *      Author: kara
 */
 
#include "globals.hh"
#include <stdio.h>
#include <stdlib.h>
#include <cmath>
#include "G4KDTree.hh"
#include "G4KDMap.hh"
#include "G4KDNode.hh"
#include "G4KDTreeResult.hh"
#include <list>
#include <iostream>
 
using namespace std;
 
G4Allocator<G4KDTree>*& G4KDTree::fgAllocator()
{
    G4ThreadLocalStatic G4Allocator<G4KDTree>* _instance = nullptr;
    return _instance;
}
 
//______________________________________________________________________
// KDTree methods
G4KDTree::G4KDTree(size_t k) :
    fKDMap(new G4KDMap(k))
{
  fDim = k;
  fRoot = 0;
  fRect = 0;
  fNbNodes = 0;
  fNbActiveNodes = 0;
}
 
G4KDTree::~G4KDTree()
{
  if (fRoot) __Clear_Rec(fRoot);
  fRoot = 0;
 
  if (fRect)
  {
    delete fRect;
    fRect = 0;
  }
 
  if (fKDMap) delete fKDMap;
}
 
void* G4KDTree::operator new(size_t)
{
  if (!fgAllocator()) fgAllocator() = new G4Allocator<G4KDTree>;
  return (void *) fgAllocator()->MallocSingle();
}
 
void G4KDTree::operator delete(void *aNode)
{
  fgAllocator()->FreeSingle((G4KDTree*) aNode);
}
 
void G4KDTree::Print(std::ostream& out) const
{
  if (fRoot) fRoot->Print(out);
}
 
void G4KDTree::Clear()
{
  __Clear_Rec(fRoot);
  fRoot = 0;
  fNbNodes = 0;
 
  if (fRect)
  {
    delete fRect;
    fRect = 0;
  }
}
 
void G4KDTree::__Clear_Rec(G4KDNode_Base* node)
{
  if (!node) return;
 
  if (node->GetLeft()) __Clear_Rec(node->GetLeft());
  if (node->GetRight()) __Clear_Rec(node->GetRight());
 
  delete node;
}
 
void G4KDTree::__InsertMap(G4KDNode_Base *node)
{
  fKDMap->Insert(node);
}
 
void G4KDTree::Build()
{
  size_t Nnodes = fKDMap->GetSize();
 
  G4cout << "********************" << G4endl;
  G4cout << "template<typename PointT> G4KDTree<PointT>::Build" << G4endl;
  G4cout << "Map size = " << Nnodes << G4endl;
 
  G4KDNode_Base* root = fKDMap->PopOutMiddle(0);
 
  if(root == 0) return;
 
  fRoot = root;
  fNbActiveNodes++;
  fRect = new HyperRect(fDim);
  fRect->SetMinMax(*fRoot, *fRoot);
 
  Nnodes--;
 
  G4KDNode_Base* parent = fRoot;
 
  for (size_t n = 0; n < Nnodes; n += fDim)
  {
    for (size_t dim = 0; dim < fDim; dim++)
    {
      G4KDNode_Base* node = fKDMap->PopOutMiddle(dim);
      if (node)
      {
        parent->Insert(node);
        fNbActiveNodes++;
        fRect->Extend(*node);
        parent = node;
      }
    }
  }
}
 
G4KDTreeResultHandle G4KDTree::Nearest(G4KDNode_Base* node)
{
  //    G4cout << "Nearest(node)" << G4endl ;
  if (!fRect)
  {
    G4cout << "Tree empty" << G4endl;
    return 0;
  }
 
  std::vector<G4KDNode_Base* > result;
  double dist_sq = DBL_MAX;
 
  /* Duplicate the bounding hyperrectangle, we will work on the copy */
  HyperRect *newrect = new HyperRect(*fRect);
 
  /* Search for the nearest neighbour recursively */
  int nbresult = 0;
 
  __NearestToNode(node, fRoot, *node, result, &dist_sq, newrect, nbresult);
 
  /* Free the copy of the hyperrect */
  delete newrect;
 
  /* Store the result */
  if (!result.empty())
  {
    G4KDTreeResultHandle rset(new G4KDTreeResult(this));
    int j = 0;
    while (j<nbresult)
    {
      rset->Insert(dist_sq, result[j]);
      j++;
    }
    rset->Rewind();
 
    return rset;
  }
  else
  {
    return 0;
  }
}
 
G4KDTreeResultHandle G4KDTree::NearestInRange(G4KDNode_Base* node,
                                              const double& range)
{
  if (!node) return 0;
  int ret(-1);
 
  G4KDTreeResult *rset = new G4KDTreeResult(this);
 
  const double range_sq = sqr(range);
 
  if ((ret = __NearestInRange(fRoot, *node, range_sq, range, *rset, 0, node)) == -1)
  {
    delete rset;
    return 0;
  }
  rset->Sort();
  rset->Rewind();
  return rset;
}