src/kdtree.cc
branchpyrit
changeset 9 3239f749e394
parent 7 bf17f9f84c91
child 10 f9fad94cd0cc
equal deleted inserted replaced
8:e6567b740c5e 9:3239f749e394
    18 		if (shapebb.R.y > bbox.R.y)  bbox.R.y = shapebb.R.y;
    18 		if (shapebb.R.y > bbox.R.y)  bbox.R.y = shapebb.R.y;
    19 		if (shapebb.R.z > bbox.R.z)  bbox.R.z = shapebb.R.z;
    19 		if (shapebb.R.z > bbox.R.z)  bbox.R.z = shapebb.R.z;
    20 	}
    20 	}
    21 };
    21 };
    22 
    22 
    23 void KdNode::subdivide(BBox bbox, int depth, int count)
    23 void KdNode::subdivide(BBox bbox, int depth)
    24 {
    24 {
    25 	int axis;
       
    26 	float pos;
       
    27 
       
    28 	//if (stopcriterionmet()) return
       
    29 
       
    30 	// choose split axis
    25 	// choose split axis
    31 	axis = 0;
    26 	axis = 0;
    32 	if (bbox.R.y - bbox.L.y > bbox.R.x - bbox.L.x)
    27 	if (bbox.R.y - bbox.L.y > bbox.R.x - bbox.L.x)
    33 		axis = 1;
    28 		axis = 1;
    34 	if (bbox.R.z - bbox.L.z > bbox.R.y - bbox.L.y)
    29 	if (bbox.R.z - bbox.L.z > bbox.R.y - bbox.L.y)
    47 		splitlist.push_back(SplitPos(sh->bbox.R.cell[axis]));
    42 		splitlist.push_back(SplitPos(sh->bbox.R.cell[axis]));
    48 	}
    43 	}
    49 	sort(splitlist.begin(), splitlist.end());
    44 	sort(splitlist.begin(), splitlist.end());
    50 
    45 
    51 	// find all posible splits
    46 	// find all posible splits
    52 	SplitList::iterator split;
    47 	SplitList::iterator spl;
    53 	int rest;
    48 	int rest;
    54 	for (split = splitlist.begin(); split != splitlist.end(); split++)
    49 	for (spl = splitlist.begin(); spl != splitlist.end(); spl++)
    55 	{
    50 	{
    56 		for (sh = sslist.begin(), rest = sslist.size(); sh != sslist.end(); sh++, rest--)
    51 		for (sh = sslist.begin(), rest = sslist.size(); sh != sslist.end(); sh++, rest--)
    57 		{
    52 		{
    58 			if (sh->hasMark())
    53 			if (sh->hasMark())
    59 				continue;
    54 				continue;
       
    55 			// if shape is completely contained in split plane
       
    56 			if (spl->pos == sh->bbox.L.cell[axis] == sh->bbox.R.cell[axis])
       
    57 			{
       
    58 				if (spl->pos - bbox.L.cell[axis] < bbox.R.cell[axis] - spl->pos)
       
    59 				{
       
    60 					// left subcell is smaller -> if not empty, put shape here
       
    61 					if (spl->lnum)
       
    62 						spl->lnum++;
       
    63 					else
       
    64 						spl->rnum++;
       
    65 				} else {
       
    66 					// right subcell is smaller
       
    67 					if (spl->rnum)
       
    68 						spl->rnum++;
       
    69 					else
       
    70 						spl->lnum++;
       
    71 				}
       
    72 			} else
    60 			// if shape is on left side of split plane
    73 			// if shape is on left side of split plane
    61 			if (sh->bbox.R.cell[axis] <= split->pos)
    74 			if (sh->bbox.R.cell[axis] <= spl->pos)
    62 			{
    75 			{
    63 				sh->setMark();
    76 				sh->setMark();
    64 				split->lnum++;
    77 				spl->lnum++;
    65 			}
    78 			} else
    66 			// if shape is completely contained in split plane
    79 			// if shape occupies both sides of split plane
    67 			if (sh->bbox.L.cell[axis] == sh->bbox.R.cell[axis] == split->pos)
    80 			if (sh->bbox.L.cell[axis] < spl->pos && sh->bbox.R.cell[axis] > spl->pos)
    68 			{
    81 			{
    69 				if (split->pos - bbox.L.cell[axis] < bbox.R.cell[axis] - split->pos)
    82 				spl->lnum++;
    70 				{
    83 				spl->rnum++;
    71 					// left subcell is smaller -> if not empty, put shape here
    84 			} else
    72 					if (split->lnum)
       
    73 						split->lnum++;
       
    74 					else
       
    75 						split->rnum++;
       
    76 				} else {
       
    77 					// right subcell is smaller
       
    78 					if (split->rnum)
       
    79 						split->rnum++;
       
    80 					else
       
    81 						split->lnum++;
       
    82 				}
       
    83 			}
       
    84 			// if shape is on right side of split plane
    85 			// if shape is on right side of split plane
    85 			if (sh->bbox.L.cell[axis] >= split->pos)
    86 			if (sh->bbox.L.cell[axis] >= spl->pos)
    86 			{
    87 			{
    87 				split->rnum += rest;
    88 				spl->rnum += rest;
    88 				break;
    89 				break;
    89 			}
       
    90 			// if shape occupies both sides of split plane
       
    91 			if (sh->bbox.L.cell[axis] < split->pos && sh->bbox.R.cell[axis] > split->pos)
       
    92 			{
       
    93 				split->lnum++;
       
    94 				split->rnum++;
       
    95 			}
    90 			}
    96 		}
    91 		}
    97 	}
    92 	}
    98 
    93 
    99 	// choose best split pos
    94 	// choose best split pos
   100 	// ... //
    95 	const float K = 1.4; // constant, K = cost of traversal / cost of ray-triangle intersection
       
    96 	float SAV = 2*(bbox.w()*bbox.h() + bbox.w()*bbox.d() + bbox.h()*bbox.d()); // surface area of node
       
    97 	float cost = SAV * (K + shapes.size()); // initial cost = non-split cost
       
    98 	SplitPos *splitpos = NULL;
       
    99 	leaf = true;
       
   100 	BBox lbb = bbox;
       
   101 	BBox rbb = bbox;
       
   102 	for (spl = splitlist.begin(); spl != splitlist.end(); spl++)
       
   103 	{
       
   104 		// calculate SAH cost of this split
       
   105 		lbb.R.cell[axis] = spl->pos;
       
   106 		rbb.L.cell[axis] = spl->pos;
       
   107 		float SAL = 2*(lbb.w()*lbb.h() + lbb.w()*lbb.d() + lbb.h()*lbb.d());
       
   108         float SAR = 2*(rbb.w()*rbb.h() + rbb.w()*rbb.d() + rbb.h()*rbb.d());
       
   109         float splitcost = K + SAL/SAV*(K+spl->lnum) + SAR/SAV*(K+spl->rnum);
   101 
   110 
   102 	KdNode *children = new KdNode[2];
   111 		if (splitcost < cost)
   103 
   112 		{
   104 	ShapeList::iterator shape;
   113 			leaf = false;
   105 	for (shape = shapes.begin(); shape != shapes.end(); shape++)
   114 			cost = splitcost;
   106 	{
   115 			splitpos = &*spl;
   107 		if (((Triangle*)*shape)->A.cell[axis-1] < pos
   116 		}
   108 		|| ((Triangle*)*shape)->B.cell[axis-1] < pos
       
   109 		|| ((Triangle*)*shape)->C.cell[axis-1] < pos)
       
   110 			children[0].addShape(*shape);
       
   111 		if (((Triangle*)*shape)->A.cell[axis-1] >= pos
       
   112 		|| ((Triangle*)*shape)->B.cell[axis-1] >= pos
       
   113 		|| ((Triangle*)*shape)->C.cell[axis-1] >= pos)
       
   114 			children[1].addShape(*shape);
       
   115 	}
   117 	}
   116 
   118 
   117 	bbox.R[axis-1] = pos;
   119 	if (leaf)
   118 	children[0].subdivide(bbox, depth+1, children[0].shapes.size());
   120 		return;
   119 
   121 
   120 	bbox.L[axis-1] = pos;
   122 	split = splitpos->pos;
   121 	children[1].subdivide(bbox, depth+1, children[1].shapes.size());
   123 	float lnum = splitpos->lnum;
       
   124 	float rnum = splitpos->rnum;
       
   125 
       
   126 	// split this node
       
   127 	children = new KdNode[2];
       
   128 	int state = 0;
       
   129 	for (sh = sslist.begin(); sh != sslist.end(); sh++)
       
   130 	{
       
   131 		// if shape is on left side of split plane
       
   132 		if (state == 1)
       
   133 		{ // only right
       
   134 			children[1].addShape(sh->shape);
       
   135 			continue;
       
   136 		}
       
   137 		if (state == 0)
       
   138 		{
       
   139 			if (sh->bbox.R.cell[axis] < split)
       
   140 			{ // left
       
   141 				children[0].addShape(sh->shape);
       
   142 			} else
       
   143 			if (sh->bbox.R.cell[axis] > split)
       
   144 			{
       
   145 				if (sh->bbox.L.cell[axis] < split)
       
   146 				{ // both
       
   147 					children[0].addShape(sh->shape);
       
   148 					children[1].addShape(sh->shape);
       
   149 				} else
       
   150 				{ // right
       
   151 					children[1].addShape(sh->shape);
       
   152 					state = 1;
       
   153 				}
       
   154 			} else
       
   155 			{ // R == split
       
   156 				if (sh->bbox.L.cell[axis] < split)
       
   157 				{ // left
       
   158 					children[0].addShape(sh->shape);
       
   159 				} else
       
   160 				{ // contained
       
   161 					if (split - bbox.L.cell[axis] < bbox.R.cell[axis] - split)
       
   162 					{
       
   163 						// left subcell is smaller -> if not empty, put shape here
       
   164 						if (lnum)
       
   165 							children[0].addShape(sh->shape);
       
   166 						else
       
   167 							children[1].addShape(sh->shape);
       
   168 					} else {
       
   169 						// right subcell is smaller
       
   170 						if (rnum)
       
   171 							children[1].addShape(sh->shape);
       
   172 						else
       
   173 							children[0].addShape(sh->shape);
       
   174 					}
       
   175 				}
       
   176 			}
       
   177 		}
       
   178 	}
       
   179 
       
   180 	lbb.R.cell[axis] = split;
       
   181 	rbb.L.cell[axis] = split;
       
   182 
       
   183 	children[0].subdivide(lbb, depth+1);
       
   184 	children[1].subdivide(rbb, depth+1);
   122 }
   185 }
   123 
   186 
   124 void KdTree::build()
   187 void KdTree::build()
   125 {
   188 {
   126 	root = new KdNode();
   189 	root = new KdNode();
   127 	root->shapes = shapes;
   190 	root->shapes = shapes;
   128 	root->subdivide(bbox, 0, shapes.size());
   191 	root->subdivide(bbox, 0);
   129 }
   192 }