C++ BVHNode类代码示例

/**
 * Recursively intersects Ray Charles with his node and returns the distance to
 * the closest intersection and a stores the id of the sphere in sphereId.
 */
inline float Intersect(const Ray charles, float t, 
                       const BVHNode& node, const vector<BVHNode>& nodes, 
                       const vector<Sphere> spheres, unsigned int &sphereId) {

    if (node.GetType() == BVHNode::LEAF) {
        // Intersect leaf
        return Exhaustive(charles, t, node, spheres, sphereId);
    } else {
        // Traverse further
        const BVHNode left = nodes[node.GetLeftChild()];
        float tLeft;
        if (!left.aabb.ClosestIntersection(charles, tLeft)) tLeft = 1e32f;
        
        const BVHNode right = nodes[node.GetRightChild()];
        float tRight;
        if (!right.aabb.ClosestIntersection(charles, tRight)) tRight = 1e32f;
    
        if (tLeft < tRight) { // Intersect left first
            if (tLeft < t) t = Intersect(charles, t, left, nodes, spheres, sphereId);
            if (tRight < t) t = Intersect(charles, t, right, nodes, spheres, sphereId);
        } else { // Intersect right first
            if (tRight < t) t = Intersect(charles, t, right, nodes, spheres, sphereId);
            if (tLeft < t) t = Intersect(charles, t, left, nodes, spheres, sphereId);
        }

        return t;
    }
}

void BVH::build(Progress& progress)
{
	progress.set_substatus("Building BVH");

	/* build nodes */
	BVHBuild bvh_build(objects,
	                   pack.prim_type,
	                   pack.prim_index,
	                   pack.prim_object,
	                   params,
	                   progress);
	BVHNode *root = bvh_build.run();

	if(progress.get_cancel()) {
		if(root) root->deleteSubtree();
		return;
	}

	/* pack triangles */
	progress.set_substatus("Packing BVH triangles and strands");
	pack_primitives();

	if(progress.get_cancel()) {
		root->deleteSubtree();
		return;
	}

	/* pack nodes */
	progress.set_substatus("Packing BVH nodes");
	pack_nodes(root);

	/* free build nodes */
	root->deleteSubtree();
}

BVHNode* BVHBuild::run()
{
	BVHRange root;

	/* add references */
	add_references(root);

	if(progress.get_cancel())
		return NULL;

	/* init spatial splits */
	if(params.top_level) /* todo: get rid of this */
		params.use_spatial_split = false;

	spatial_min_overlap = root.bounds().safe_area() * params.spatial_split_alpha;
	spatial_right_bounds.clear();
	spatial_right_bounds.resize(max(root.size(), (int)BVHParams::NUM_SPATIAL_BINS) - 1);

	/* init progress updates */
	progress_start_time = time_dt();
	progress_count = 0;
	progress_total = references.size();
	progress_original_total = progress_total;

	prim_segment.resize(references.size());
	prim_index.resize(references.size());
	prim_object.resize(references.size());

	/* build recursively */
	BVHNode *rootnode;

	if(params.use_spatial_split) {
		/* singlethreaded spatial split build */
		rootnode = build_node(root, 0);
	}
	else {
		/* multithreaded binning build */
		BVHObjectBinning rootbin(root, (references.size())? &references[0]: NULL);
		rootnode = build_node(rootbin, 0);
		task_pool.wait_work();
	}

	/* delete if we canceled */
	if(rootnode) {
		if(progress.get_cancel()) {
			rootnode->deleteSubtree();
			rootnode = NULL;
		}
		else if(!params.use_spatial_split) {
			/*rotate(rootnode, 4, 5);*/
			rootnode->update_visibility();
		}
	}

	return rootnode;
}

//--------------------------------------------------------------
// return count of descendents
int BVHNode::nodeCount()
{
  int count = 1;
  for (int i = 0; i < m_children.length(); i++)
  {
    BVHNode* child = (BVHNode*) m_children[i];
    count += child->nodeCount();
  }
  return count;
}

float BVHNode::computeSubtreeSAHCost(const BVHParams& p, float probability) const
{
	float SAH = probability * p.cost(num_children(), num_triangles());

	for(int i = 0; i < num_children(); i++) {
		BVHNode *child = get_child(i);
		SAH += child->computeSubtreeSAHCost(p, probability * child->bounds.safe_area()/bounds.safe_area());
	}

	return SAH;
}

void PhotonAccelerator::build(const std::vector<Hitpoint*>& hitpoints)
{
	hps = hitpoints;
	nodes.reserve(2 * hitpoints.size() - 1);
	BVHNode root;
	AABB& world_box = root.getAABB();
	std::for_each(hps.begin(), hps.end(), [&] (Hitpoint* hp) {
		AABB aabb(hp->is.mPosition);
		aabb.grow(hp->radius);
		world_box.include(aabb);
	});
	nodes.push_back(root);
	build_recursive(0, hps.size(), &nodes[0], 0);
}

void IKTree::reset(int frame)
{
  MT_Quaternion q;
  BVHNode *n;

  for (int i=0; i<numBones; i++) {
    bone[i].pos = origin;
    bone[i].lRot = identity;
    bone[i].gRot = identity;
    n = bone[i].node;
    Rotation rot=n->frameData(frame).rotation();
    Position pos=n->frameData(frame).position();

    for (int k=0; k<n->numChannels; k++) {  // rotate each axis in order
      q = identity;
      switch (n->channelType[k]) {
      case BVH_XROT:
        q.setRotation(xAxis, rot.x * M_PI / 180);
      break;
      case BVH_YROT:
        q.setRotation(yAxis, rot.y * M_PI / 180);
      break;
      case BVH_ZROT:
        q.setRotation(zAxis, rot.z * M_PI / 180);
      break;
      case BVH_XPOS: bone[i].pos[0] = pos.x; break;
      case BVH_YPOS: bone[i].pos[1] = pos.y; break;
      case BVH_ZPOS: bone[i].pos[2] = pos.z; break;
      }
      bone[i].lRot = q * bone[i].lRot;
    }
/*
    for (int k=0; k<3; k++) {  // rotate each axis in order
      rad = n->frame[frame][k] * M_PI / 180;
      q = identity;
      switch (n->channelType[k]) {
      case BVH_XROT: q.setRotation(xAxis, rad); break;
      case BVH_YROT: q.setRotation(yAxis, rad); break;
      case BVH_ZROT: q.setRotation(zAxis, rad); break;
      case BVH_XPOS: bone[i].pos[0] = n->frame[frame][k]; break;
      case BVH_YPOS: bone[i].pos[1] = n->frame[frame][k]; break;
      case BVH_ZPOS: bone[i].pos[2] = n->frame[frame][k]; break;
      }
      bone[i].lRot = q * bone[i].lRot;
    }
*/
  }
  updateBones(0);
}

inline float Exhaustive(const Ray charles, float t, const BVHNode& node, 
                        const vector<Sphere> spheres, unsigned int &sphereId) {

    exhaustives += node.GetPrimitiveRange();

    for (unsigned int p = node.GetFirstPrimitive(); 
         p < node.GetFirstPrimitive() + node.GetPrimitiveRange(); ++p) {

        const Sphere sphere = spheres[p];
        const float tSphere = sphere.Intersect(charles);
        if (0 < tSphere && tSphere < t) {
            sphereId = p;
            t = tSphere;
        }
    }

    return t;
}

/**
 * Output
 */
std::string BVHNode::toString() const {
	using std::stringstream;
	using std::endl;

	stringstream ss;
	ss << "\t" << id << " [width=.5,height=1,style=filled,color=\".5 .5 .5\",shape=box,label=\"" << id << ",\\n" << box.min << ",\\n" << box.max << "\"];" << std::endl;
	if (leaf) {
		for(uint i = 0; i < 2; ++i) {
			Object* o = childs[i].leaf;
			if (o) ss << "\t" << id << " -> " << "box" << o->id << std::endl << o->toString();
		}
		ss << std::endl;
	} else {
		for(uint i = 0; i < 2; ++i) {
			BVHNode* o = childs[i].node;
			ss << "\t" << id << " -> " << o->id << std::endl << o->toString();
		}
		ss << std::endl;
	}
	return ss.str();
}

/**
 * Iteratively intersect rays with the given nodes and returns the distance to
 * the closest intersection and stores the id of the intersected sphere in
 * sphereId.
 */
inline float ItrIntersect(const Ray charles, const vector<BVHNode>& nodes, 
                          const vector<Sphere> spheres, unsigned int &sphereId) {
    
    // cout << "=== ItrIntersect:" << charles.ToString() << " ===" << endl;
    
    sphereId = -1;
    float t = 1e30f;
    vector<std::pair<int, float> > stack(60);
    int i = 0;
    stack[i] = std::pair<int, float>(0, 0.0f);

    do {
        // cout << "\n----" << endl;
        // for (int j = i; j >= 0; --j)
        //     cout << "| " << j << ": " << stack[j].second << " - " << nodes[stack[j].first].ToString() << endl;
        // cout << "----" << endl;

        std::pair<int, float> next = stack[i]; --i;
        if (t < next.second) {
            // cout << i << "discard: " << nodes[next.first].ToString() << endl;
            continue;
        }
        
        BVHNode currentNode = nodes[next.first];

        if (currentNode.GetType() == BVHNode::LEAF) // Intersect leaf
            t = Exhaustive(charles, t, currentNode, spheres, sphereId);
        else {
            rayAABBCheck += 2;
        
            const BVHNode& left = nodes[currentNode.GetLeftChild()];
            float tLeft;
            if (!left.aabb.ClosestIntersection(charles, tLeft)) tLeft = 1e32f;
            
            const BVHNode& right = nodes[currentNode.GetRightChild()];
            float tRight;
            if (!right.aabb.ClosestIntersection(charles, tRight)) tRight = 1e32f;

            if (tLeft < tRight) { // Intersect left first
                if (tRight < t) stack[++i] = std::pair<int, float>(currentNode.GetRightChild(), tRight);
                if (tLeft < t) stack[++i] = std::pair<int, float>(currentNode.GetLeftChild(), tLeft);
            } else { // Intersect right first
                if (tLeft < t) stack[++i] = std::pair<int, float>(currentNode.GetLeftChild(), tLeft);
                if (tRight < t) stack[++i] = std::pair<int, float>(currentNode.GetRightChild(), tRight);
            }
        }
    } while (i >= 0);

    return t;
}

BVHNode * BVHTree::buildSubTree(std::vector<Shape *> shapes, int depth) {
    int size = (int) shapes.size();
    int axis = depth % dims; // 0-2

    // Construct leaf BVHNode containing shapes
    if (size == shapesPerLeaf) {
        BVHNode *leaf = new BVHNode();
        leaf->addShape(shapes[0]); // Add shape and init bbox
        leaf->leaf = true;
        return leaf;
    } else if (size == 0) {
        return NULL;
    }
    // Otherwise split the shapes into two subsets, and divide them amongst
    // the left and right child nodes

    // Sort shapes based on the current axis
    std::sort(shapes.begin(), shapes.end(), comparators[axis]);

    // Find the median 
    int median = size / 2;
    std::vector<Shape *>::iterator mid = shapes.begin() + median;

    // Construct tree BVHNode
    BVHNode *treeNode = new BVHNode();

    // Construct left child BVHNode
    std::vector<Shape *> lShapes(shapes.begin(), mid);
    treeNode->left = buildSubTree(lShapes, depth + 1);

    // Construct right child BVHNode
    std::vector<Shape *> rShapes(mid, shapes.end());
    treeNode->right = buildSubTree(rShapes, depth + 1);

    // Store the bbox for the treeNode
    treeNode->bbox = treeNode->left->bbox + treeNode->right->bbox;
    return treeNode;
}

void PhotonAccelerator::accumulate(const Intersection& is, Color flux)
{
	BVHNode* currentNode = &nodes[0];
	AABB& box = currentNode->getAABB();
	Point3D point = is.mPosition;
	if (!box.inside(point))
		return;
	std::stack<BVHNode> intersect_stack;
	while (true) {
		if (currentNode->isLeaf()) {
			Vector3D v(point);
			for (unsigned int i = currentNode->getIndex(); i < currentNode->getIndex() + currentNode->getNObjs(); ++i) {
				Hitpoint& hp = *hps[i];
				if ((v - hp.is.mPosition).length() <= hp.radius && hp.is.mNormal * is.mNormal >= 0) {
					hp.newPhotonCount += 1;
					hp.totalFlux += flux * hp.is.mMaterial->evalBRDF(hp.is, -is.mRay.dir);
				}
			}
		} else {
			BVHNode& left_node = nodes[currentNode->getIndex()];
			BVHNode& right_node = nodes[currentNode->getIndex() + 1];
			AABB& left_box = left_node.getAABB();
			AABB& right_box = right_node.getAABB();

			bool leftHit, rightHit;
			leftHit = left_box.inside(point);
			rightHit = right_box.inside(point);
			if (leftHit && rightHit) {
				currentNode = &left_node;
				intersect_stack.push(right_node);
				continue;
			} else if (leftHit) {
				currentNode = &left_node;
				continue;
			} else if (rightHit) {
				currentNode = &right_node;
				continue;
			}
		}
		if (intersect_stack.empty())
			return;
		currentNode = &intersect_stack.top();
		intersect_stack.pop();
	}
}

BVHNode* BVH::RecursiveBuild(uint32_t start, uint32_t end, uint32_t depth)
{
    maxDepth = fmaxf(depth, maxDepth);
    totalNodes++;
    BVHNode* node = new BVHNode;

    //compute bounds of all primitives in BVH node
    BBox bbox;
    for(auto i = start; i < end; ++i)
        bbox = Union(bbox, workList[i].bounds);

    uint32_t nPrims = end - start;
    //if number of primitives are less than threshold, create leaf node
    if(nPrims <= MAX_LEAF_PRIM_NUM)
    {
        uint32_t firstPrimOffset = orderedPrims.size();
        for(auto i = start; i < end; ++i)
        {
            auto pIdx = workList[i].pIdx;
            orderedPrims.push_back(mesh.faces[pIdx]);
        }
        node->InitLeaf(firstPrimOffset, nPrims, bbox);
    }
    else
    {
        //compute bound of primitive centroids, choose split dimension
        BBox centroidBounds;
        for(auto i = start; i < end; ++i)
            centroidBounds = Union(centroidBounds, workList[i].bounds.bcenter);

        // split along max span axis
        int dim = centroidBounds.MaxExtent();

        //partition primitives into two sets and build children
        uint32_t mid = (end + start) / 2;
        // if max span axis is too small, create a leaf node
        if((centroidBounds.bmax[dim] - centroidBounds.bmin[dim]) < 1e-4)
        {
            uint32_t firstPrimOffset = orderedPrims.size();
            for(auto i = start; i < end; ++i)
            {
                auto pIdx = workList[i].pIdx;
                orderedPrims.push_back(mesh.faces[pIdx]);
            }
            node->InitLeaf(firstPrimOffset, nPrims, bbox);

            return node;
        }

        //partition primitives based on SAH
        std::vector<BucketInfo> buckets(nBuckets);
        float extent = centroidBounds.bmax[dim] - centroidBounds.bmin[dim];
        for(auto i = start; i < end; ++i)
        {
            uint32_t b = nBuckets * ((workList[i].bounds.bcenter[dim] - centroidBounds.bmin[dim]) / extent);
            if(b == nBuckets) b -= 1;
            buckets[b].count++;
            buckets[b].bounds = Union(buckets[b].bounds, workList[i].bounds);
        }

        //compute costs for splitting after each bucket
        float cost[nBuckets - 1];
        for(auto i = 0; i < nBuckets - 1; ++i)
        {
            BBox b0, b1;
            int count0 = 0, count1 = 0;

            for(auto j = 0; j <= i; ++j)
            {
                b0 = Union(b0, buckets[j].bounds);
                count0 += buckets[j].count;
            }
            for(auto j = i + 1; j < nBuckets; ++j)
            {
                b1 = Union(b1, buckets[j].bounds);
                count1 += buckets[j].count;
            }

            cost[i] = (count0 * b0.SurfaceArea() + count1 * b1.SurfaceArea()) / bbox.SurfaceArea();
        }

        //find best split
        float minCost = cost[0];
        uint32_t bestSplit = 0;
        for(auto i = 1; i < nBuckets - 1; ++i)
        {
            if(cost[i] < minCost)
            {
                minCost = cost[i];
                bestSplit = i;
            }
        }

        //either create leaf or split at selected SAH bucket
        if(nPrims > MAX_LEAF_PRIM_NUM || minCost < nPrims)
        {
            auto compare = [&](BVHPrimitiveInfo& p) {
                auto b = nBuckets * ((p.bounds.bcenter[dim] - centroidBounds.bmin[dim]) / extent);
                b = (b == nBuckets) ? (b - 1) : b;
                return b <= bestSplit;
//.........这里部分代码省略.........

void Blender::EvaluateRelativeLimbWeights(QList<TimelineTrail*>* trails, int trailsCount)       //TODO: into below method?
{
  int minPosIndex = 999999999;
  int maxPosIndex = -1;
  TrailItem** currentItems = new TrailItem*[trailsCount];

  //Find first and last occupied time-line position. And initialize currentItems with first non-shadows
  for(int i=0; i<trailsCount; i++)
  {
    TrailItem* firstItem = trails->at(i)->firstItem();
    if(firstItem==NULL)
    {
      currentItems[i] = NULL;
      continue;
    }

    //Shadow items don't matter much because the nearest/furthest item is always non-shadow
    if(firstItem->beginIndex() < minPosIndex)
      minPosIndex = firstItem->beginIndex();

    TrailItem* lastItem = trails->at(i)->lastItem();
    if(lastItem->endIndex() > maxPosIndex)
      maxPosIndex = trails->at(i)->lastItem()->endIndex();

    while(firstItem->isShadow())
      firstItem = firstItem->nextItem();
    currentItems[i] = firstItem;
  }
  if(maxPosIndex == -1)                       //Time-line is empty
    return;

  int curPosIndex = minPosIndex;
  //Position pseudo-node is not included, but it's OK as it can't be highlighted
  QStringList boneNames = BVH::getValidNodeNames();
  while(curPosIndex<=maxPosIndex)
  {
    foreach(QString bName, boneNames)
    {
      int sumWeight = 0;
      int weightsUsed = 0;
      bool emptyPosition = true;

      //Temporary helper structure. Value is array of 2 integers: frame index inside
      //respective animation and its frame weight
      QMap<BVHNode*, QVector<int> >* usedData = new QMap<BVHNode*, QVector<int> >();

      for(int i=0; i<trailsCount; i++)
      {
        TrailItem* item = currentItems[i];

        if(item == NULL)                                 //No more items on i-th trail
          continue;
        if(item->beginIndex() > curPosIndex)             //We're before first item of this trail
          continue;

        if(item->endIndex() < curPosIndex)               //Move to next item
        {
          currentItems[i] = item->nextItem();
          item = currentItems[i];
        }
        if(item == NULL || item->isShadow())             //Still here?
          continue;
        else
          emptyPosition = false;

        if(item->beginIndex() <= curPosIndex && item->endIndex() >= curPosIndex)
        {
          int frameIndex = curPosIndex - item->beginIndex();
          int frameWeight = item->getWeight(frameIndex);

          if(!item->getAnimation()->bones()->contains(bName))
            Announcer::Exception(NULL, "Exception: can't evaluate relative weight for limb " +bName);


//          BVHNode* limb = item->getAnimation()->bones()->value(bName);      NOT WORKING. REALLY NEED TO KNOW WHY
          //desperate DEBUG
          BVHNode* limb = item->getAnimation()->getNodeByName(bName);

          QVector<int> tempData;
          tempData << frameIndex << frameWeight;
          usedData->insert(limb, tempData);
          sumWeight += frameWeight * limb->frameData(frameIndex).weight();
          weightsUsed++;
        }
//        else Announcer::Exception(NULL, "Lame programmer exception :(");    //DEBUG
      }

      if(emptyPosition)                         //There were no valid item at this position
        break;                                  //so jump to the next one



      QMapIterator<BVHNode*, QVector<int> > iter(*usedData);
      while(iter.hasNext())
      {
        iter.next();
        BVHNode* bone = iter.key();
        QVector<int> tempData = iter.value();
        int frame = tempData[0];
        int frameW = tempData[1];
//.........这里部分代码省略.........

BVHNode *BVHParser::parseNode( Buffer & buf, BVHNode * parent )
{
	BVHNode * node = new BVHNode();
	node->parent = parent;
	
	_linearNodes.push_back(node);
	
	buf.getToken(node->name);
	
	std::string str;
	
	buf.getToken(str);
	if( str != "{" )
	{
		debugPrint("Error didnt find opening bracket\n");
		delete node;
		return 0;
	}
	
	// parse offset
	buf.getToken(str);
	if(str!="OFFSET")
	{
		debugPrint("Error didnt find OFFSET\n");
		delete node;
		return 0;
	}
	
	buf.readLine(str);
	
	sscanf(str.c_str(),"%f %f %f",&node->offset.x,&node->offset.y,&node->offset.z);
	
	//buf.getToken(str);
	if(!buf.testToken("CHANNELS"))//, <#bool caseSensitive#>)str!="CHANNELS")
	{
		debugPrint("Error didnt find CHANNELS\n");
		delete node;
		return 0;
	}

	// parse channels 
	
	int nChannels = 0;
	
	buf.readInt(&nChannels);
	
	for( int i = 0; i < nChannels; i++ )
	{
		buf.getToken(str);
		BVH_CHANNEL chan = parseChannel(str);
		if(chan == CHANNEL_UNKNOWN) 
		{
			debugPrint("Error undefined channel %s\n",str.c_str());
			delete node;
			return 0;
		}
		node->channels.push_back(chan);
	}
	
	buf.getToken(str);
	
	// parse children ( if any... )
	while( str!="}" )
	{
		if( str=="JOINT" )
		{
			BVHNode * child = parseNode(buf,node);
			if(child==0)
			{
				debugPrint("Could not add child\n");
				delete node;
				return 0;
			}
			
			node->addChild(child);
			//buf.getToken(str);
		}
		else if( str=="End" )
		{
			buf.getToken(str); // skip site
			
			BVHNode * child = parseEndSite(buf,node);
			if(child==0)
			{
				debugPrint("Could not add child\n");
				delete node;
				return 0;
			}
			
			node->addChild(child);
		}
		
		buf.getToken(str);
	}
	
	// Hacky compute bone length...
	if( node->children.size() )
	{
		BVHNode * child = node->children[0];
	}
//.........这里部分代码省略.........