// ================================================================ // $Id: lsds_kdtree.cpp,v 1.18 2005/04/16 09:34:21 bittner Exp $ // **************************************************************** /** The KD tree based LSDS */ // Initial coding by /** @author Jiri Bittner */ // Standard headers #include #include #include #include #include #include "VssTree.h" #include "Environment.h" #include "VssRay.h" #include "Intersectable.h" #include "Ray.h" #define DEBUG_DIR_SPLIT 0 // Static variables int VssTreeLeaf::mailID = 0; inline void AddObject2Pvs(Intersectable *object, const int side, int &pvsBack, int &pvsFront) { if (!object) return; if (side <= 0) { if (!object->Mailed() && !object->Mailed(2)) { pvsBack++; if (object->Mailed(1)) object->Mail(2); else object->Mail(); } } if (side >= 0) { if (!object->Mailed(1) && !object->Mailed(2)) { pvsFront++; if (object->Mailed()) object->Mail(2); else object->Mail(1); } } } // Constructor VssTree::VssTree() { environment->GetIntValue("VssTree.maxDepth", termMaxDepth); environment->GetIntValue("VssTree.minPvs", termMinPvs); environment->GetIntValue("VssTree.minRays", termMinRays); environment->GetFloatValue("VssTree.maxRayContribution", termMaxRayContribution); environment->GetFloatValue("VssTree.maxCostRatio", termMaxCostRatio); environment->GetFloatValue("VssTree.minSize", termMinSize); termMinSize = sqr(termMinSize); environment->GetFloatValue("VssTree.refDirBoxMaxSize", refDirBoxMaxSize); refDirBoxMaxSize = sqr(refDirBoxMaxSize); environment->GetFloatValue("VssTree.epsilon", epsilon); environment->GetFloatValue("VssTree.ct_div_ci", ct_div_ci); environment->GetFloatValue("VssTree.maxTotalMemory", maxTotalMemory); environment->GetFloatValue("VssTree.maxStaticMemory", maxStaticMemory); float refDirAngle; environment->GetFloatValue("VssTree.refDirAngle", refDirAngle); environment->GetIntValue("VssTree.accessTimeThreshold", accessTimeThreshold); //= 1000; environment->GetIntValue("VssTree.minCollapseDepth", minCollapseDepth); // int minCollapseDepth = 4; // pRefDirThresh = cos(0.5*M_PI - M_PI*refDirAngle/180.0); // cosRefDir = cos(M_PI*refDirAngle/180.0); // sinRefDir = sin(M_PI*refDirAngle/180.0); // split type char sname[128]; environment->GetStringValue("VssTree.splitType", sname); string name(sname); if (name.compare("regular") == 0) splitType = ESplitRegular; else if (name.compare("heuristic") == 0) splitType = ESplitHeuristic; else { cerr<<"Invalid VssTree split type "<GetBoolValue("VssTree.randomize", randomize); root = NULL; splitCandidates = new vector; } VssTree::~VssTree() { if (root) delete root; } void VssStatistics::Print(ostream &app) const { app << "===== VssTree statistics ===============\n"; app << "#N_RAYS ( Number of rays )\n" << rays <IsActive()) { Intersectable *object; #if BIDIRECTIONAL_RAY object = (*ri).mRay->mOriginObject; if (object && !object->Mailed()) { pvsSize++; object->Mail(); } #endif object = (*ri).mRay->mTerminationObject; if (object && !object->Mailed()) { pvsSize++; object->Mail(); } } mPvsSize = pvsSize; mValidPvs = true; } } void VssTree::Construct( VssRayContainer &rays, AxisAlignedBox3 *forcedBoundingBox ) { stat.Start(); maxMemory = maxStaticMemory; if (root) delete root; root = new VssTreeLeaf(NULL, rays.size()); // first construct a leaf that will get subdivide VssTreeLeaf *leaf = (VssTreeLeaf *) root; stat.nodes = 1; bbox.Initialize(); dirBBox.Initialize(); for(VssRayContainer::const_iterator ri = rays.begin(); ri != rays.end(); ri++) { leaf->AddRay(VssTreeNode::RayInfo(*ri)); bbox.Include((*ri)->GetOrigin()); bbox.Include((*ri)->GetTermination()); dirBBox.Include(Vector3( (*ri)->GetDirParametrization(0), (*ri)->GetDirParametrization(1), 0 ) ); } if ( forcedBoundingBox ) bbox = *forcedBoundingBox; cout<<"Bbox = "< maxMemory ) { // count statistics on unprocessed leafs while (!tStack.empty()) { EvaluateLeafStats(tStack.top()); tStack.pop(); } break; } TraversalData data = tStack.top(); tStack.pop(); VssTreeNode *node = SubdivideNode((VssTreeLeaf *) data.node, data.bbox, backBox, frontBox ); if (result == NULL) result = node; if (!node->IsLeaf()) { VssTreeInterior *interior = (VssTreeInterior *) node; // push the children on the stack tStack.push(TraversalData(interior->back, backBox, data.depth+1)); tStack.push(TraversalData(interior->front, frontBox, data.depth+1)); } else { EvaluateLeafStats(data); } } return result; } // returns selected plane for subdivision int VssTree::SelectPlane( VssTreeLeaf *leaf, const AxisAlignedBox3 &box, float &position, int &raysBack, int &raysFront, int &pvsBack, int &pvsFront ) { int minDirDepth = 6; int axis; float costRatio; if (splitType == ESplitRegular) { costRatio = BestCostRatioRegular(leaf, axis, position, raysBack, raysFront, pvsBack, pvsFront ); } else { if (splitType == ESplitHeuristic) costRatio = BestCostRatioHeuristic(leaf, axis, position, raysBack, raysFront, pvsBack, pvsFront ); else { cerr<<"VssTree: Unknown split heuristics\n"; exit(1); } } if (costRatio > termMaxCostRatio) { // cout<<"Too big cost ratio "<rays.size(); float oldCost = pvsSize; float ratio = newCost/oldCost; // cout<<"ratio="<rays.size(); int pl=0, pr = leaf->GetPvsSize(); float minBox = box.Min(axis); float maxBox = box.Max(axis); float sizeBox = maxBox - minBox; float minBand = minBox + 0.1*(maxBox - minBox); float maxBand = minBox + 0.9*(maxBox - minBox); float sum = rr*sizeBox; float minSum = 1e20; Intersectable::NewMail(); // set all object as belonging to the fron pvs for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin(); ri != leaf->rays.end(); ri++) if ((*ri).mRay->IsActive()) { Intersectable *object = (*ri).mRay->mTerminationObject; if (object) if (!object->Mailed()) { object->Mail(); object->mCounter = 1; } else object->mCounter++; } Intersectable::NewMail(); for(vector::const_iterator ci = splitCandidates->begin(); ci < splitCandidates->end(); ci++) { VssRay *ray; switch ((*ci).type) { case SortableEntry::ERayMin: { rl++; ray = (VssRay *) (*ci).data; Intersectable *object = ray->mTerminationObject; if (object && !object->Mailed()) { object->Mail(); pl++; } break; } case SortableEntry::ERayMax: { rr--; ray = (VssRay *) (*ci).data; Intersectable *object = ray->mTerminationObject; if (object) { if (--object->mCounter == 0) pr--; } break; } } if ((*ci).value > minBand && (*ci).value < maxBand) { sum = pl*((*ci).value - minBox) + pr*(maxBox - (*ci).value); // cout<<"pos="<<(*ci).value<<"\t q=("<clear(); int requestedSize = 2*(node->rays.size()); // creates a sorted split candidates array if (splitCandidates->capacity() > 500000 && requestedSize < (int)(splitCandidates->capacity()/10) ) { delete splitCandidates; splitCandidates = new vector; } splitCandidates->reserve(requestedSize); // insert all queries for(VssTreeNode::RayInfoContainer::const_iterator ri = node->rays.begin(); ri < node->rays.end(); ri++) { bool positive = (*ri).mRay->HasPosDir(axis); splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMin : SortableEntry::ERayMax, (*ri).ExtrapOrigin(axis), (void *)&*ri) ); splitCandidates->push_back(SortableEntry(positive ? SortableEntry::ERayMax : SortableEntry::ERayMin, (*ri).ExtrapTermination(axis), (void *)&*ri) ); } stable_sort(splitCandidates->begin(), splitCandidates->end()); } void VssTree::EvaluateLeafStats(const TraversalData &data) { // the node became a leaf -> evaluate stats for leafs VssTreeLeaf *leaf = (VssTreeLeaf *)data.node; if (data.depth >= termMaxDepth) stat.maxDepthNodes++; // if ( (int)(leaf->rays.size()) < termMinCost) // stat.minCostNodes++; if ( leaf->GetPvsSize() < termMinPvs) stat.minPvsNodes++; if ( leaf->GetPvsSize() < termMinRays) stat.minRaysNodes++; if (0 && leaf->GetAvgRayContribution() > termMaxRayContribution ) stat.maxRayContribNodes++; if (SqrMagnitude(data.bbox.Size()) <= termMinSize) { stat.minSizeNodes++; } if ( (int)(leaf->rays.size()) > stat.maxRayRefs) stat.maxRayRefs = leaf->rays.size(); } VssTreeNode * VssTree::SubdivideNode( VssTreeLeaf *leaf, const AxisAlignedBox3 &box, AxisAlignedBox3 &backBBox, AxisAlignedBox3 &frontBBox ) { if ( (leaf->GetPvsSize() < termMinPvs) || (leaf->rays.size() < termMinRays) || // (leaf->GetAvgRayContribution() > termMaxRayContribution ) || (leaf->depth >= termMaxDepth) || SqrMagnitude(box.Size()) <= termMinSize ) { #if 0 if (leaf->depth >= termMaxDepth) { cout<<"Warning: max depth reached depth="<<(int)leaf->depth<<" rays="<rays.size()<rays.size() > (unsigned)termMinCost && (leaf->GetPvsSize() >= termMinPvs) && SqrMagnitude(leafBBox.Size()) > sizeThreshold) { // memory check and realese... if (GetMemUsage() > maxTotalMemory) { ReleaseMemory( pass ); } AxisAlignedBox3 backBBox, frontBBox; // subdivide the node node = SubdivideNode(leaf, leafBBox, backBBox, frontBBox ); } return node; } void VssTree::UpdateRays(VssRayContainer &remove, VssRayContainer &add ) { VssTreeLeaf::NewMail(); // schedule rays for removal for(VssRayContainer::const_iterator ri = remove.begin(); ri != remove.end(); ri++) { (*ri)->ScheduleForRemoval(); } int inactive=0; for(VssRayContainer::const_iterator ri = remove.begin(); ri != remove.end(); ri++) { if ((*ri)->ScheduledForRemoval()) // RemoveRay(*ri, NULL, false); // !!! BUG - with true it does not work correctly - aggreated delete RemoveRay(*ri, NULL, true); else inactive++; } // cout<<"all/inactive"< *affectedLeaves, const bool removeAllScheduledRays ) { stack tstack; tstack.push(RayTraversalData(root, VssTreeNode::RayInfo(ray))); RayTraversalData data; // cout<<"Number of ray refs = "<RefCount()<IsLeaf()) { // split the set of rays in two groups intersecting the // two subtrees TraverseInternalNode(data, tstack); } else { // remove the ray from the leaf // find the ray in the leaf and swap it with the last ray... VssTreeLeaf *leaf = (VssTreeLeaf *)data.node; if (!leaf->Mailed()) { leaf->Mail(); if (affectedLeaves) affectedLeaves->push_back(leaf); if (removeAllScheduledRays) { int tail = leaf->rays.size()-1; for (int i=0; i < (int)(leaf->rays.size()); i++) { if (leaf->rays[i].mRay->ScheduledForRemoval()) { // find a ray to replace it with while (tail >= i && leaf->rays[tail].mRay->ScheduledForRemoval()) { stat.removedRayRefs++; leaf->rays[tail].mRay->Unref(); leaf->rays.pop_back(); tail--; } if (tail < i) break; stat.removedRayRefs++; leaf->rays[i].mRay->Unref(); leaf->rays[i] = leaf->rays[tail]; leaf->rays.pop_back(); tail--; } } } } if (!removeAllScheduledRays) for (int i=0; i < (int)leaf->rays.size(); i++) { if (leaf->rays[i].mRay == ray) { stat.removedRayRefs++; ray->Unref(); leaf->rays[i] = leaf->rays[leaf->rays.size()-1]; leaf->rays.pop_back(); // check this ray again break; } } } } if (ray->RefCount() != 0) { cerr<<"Error: Number of remaining refs = "<RefCount()< tstack; tstack.push(RayTraversalData(root, VssTreeNode::RayInfo(ray))); RayTraversalData data; while (!tstack.empty()) { data = tstack.top(); tstack.pop(); if (!data.node->IsLeaf()) { TraverseInternalNode(data, tstack); } else { // remove the ray from the leaf // find the ray in the leaf and swap it with the last ray... VssTreeLeaf *leaf = (VssTreeLeaf *)data.node; leaf->AddRay(data.rayData); stat.addedRayRefs++; } } } void VssTree::TraverseInternalNode( RayTraversalData &data, stack &tstack) { VssTreeInterior *in = (VssTreeInterior *) data.node; if (in->axis <= VssTreeNode::SPLIT_Z) { // determine the side of this ray with respect to the plane int side = in->ComputeRayIntersection(data.rayData, data.rayData.mRay->mT); if (side == 0) { if (data.rayData.mRay->HasPosDir(in->axis)) { tstack.push(RayTraversalData(in->back, VssTreeNode::RayInfo(data.rayData.mRay, data.rayData.mMinT, data.rayData.mRay->mT)) ); tstack.push(RayTraversalData(in->front, VssTreeNode::RayInfo(data.rayData.mRay, data.rayData.mRay->mT, data.rayData.mMaxT )) ); } else { tstack.push(RayTraversalData(in->back, VssTreeNode::RayInfo(data.rayData.mRay, data.rayData.mRay->mT, data.rayData.mMaxT )) ); tstack.push(RayTraversalData(in->front, VssTreeNode::RayInfo(data.rayData.mRay, data.rayData.mMinT, data.rayData.mRay->mT)) ); } } else if (side == 1) tstack.push(RayTraversalData(in->front, data.rayData)); else tstack.push(RayTraversalData(in->back, data.rayData)); } else { // directional split if (data.rayData.mRay->GetDirParametrization(in->axis - 3) > in->position) tstack.push(RayTraversalData(in->front, data.rayData)); else tstack.push(RayTraversalData(in->back, data.rayData)); } } int VssTree::CollapseSubtree(VssTreeNode *sroot, const int time) { // first count all rays in the subtree // use mail 1 for this purpose stack tstack; int rayCount = 0; int totalRayCount = 0; int collapsedNodes = 0; #if DEBUG_COLLAPSE cout<<"Collapsing subtree"< tstack; tstack.push(root); Intersectable::NewMail(); int pvsSize = 0; while (!tstack.empty()) { VssTreeNode *node = tstack.top(); tstack.pop(); if (node->IsLeaf()) { VssTreeLeaf *leaf = (VssTreeLeaf *)node; for(VssTreeNode::RayInfoContainer::iterator ri = leaf->rays.begin(); ri != leaf->rays.end(); ri++) if ((*ri).mRay->IsActive()) { Intersectable *object; #if BIDIRECTIONAL_RAY object = (*ri).mRay->mOriginObject; if (object && !object->Mailed()) { pvsSize++; object->Mail(); } #endif object = (*ri).mRay->mTerminationObject; if (object && !object->Mailed()) { pvsSize++; object->Mail(); } } } else { VssTreeInterior *in = (VssTreeInterior *)node; if (in->axis < 3) { if (box.Max(in->axis) >= in->position ) tstack.push(in->front); if (box.Min(in->axis) <= in->position ) tstack.push(in->back); } else { // both nodes for directional splits tstack.push(in->front); tstack.push(in->back); } } } return pvsSize; } void VssTree::GetRayContributionStatistics( float &minRayContribution, float &maxRayContribution, float &avgRayContribution ) { stack tstack; tstack.push(root); minRayContribution = 1.0f; maxRayContribution = 0.0f; float sumRayContribution = 0.0f; int leaves = 0; while (!tstack.empty()) { VssTreeNode *node = tstack.top(); tstack.pop(); if (node->IsLeaf()) { leaves++; VssTreeLeaf *leaf = (VssTreeLeaf *)node; float c = leaf->GetAvgRayContribution(); if (c > maxRayContribution) maxRayContribution = c; if (c < minRayContribution) minRayContribution = c; sumRayContribution += c; } else { VssTreeInterior *in = (VssTreeInterior *)node; // both nodes for directional splits tstack.push(in->front); tstack.push(in->back); } } cout<<"sum="<