octree.cpp

用c++实现的八叉树建立

#include "stdafx.h"
#include <iostream>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "octree.h"

using namespace std;
/* ------------------------------------------------------------------------ */
/* ------------------------------------------------------------------------ */

Vec3 makeVec3( double x, double y, double z)
{
  Vec3 v3 = (Vec3) malloc(3 * sizeof(double));
  v3[0] = x; v3[1] = y; v3[2] = z;
  return v3;
}

Vec3 copyVec3( Vec3 src )
{
  Vec3 v3 = (Vec3) malloc(3 * sizeof(double));
  v3[0] = src[0]; v3[1] = src[1]; v3[2] = src[2];
  return v3;
}

/* ------------------------------------------------------------------------ */

Octree* make_octree( Vec3 min, Vec3 max )
{
  //Octree* o = (Octree*) malloc(sizeof(Octree));
	Octree* o = new Octree;
  o->min = copyVec3(min);
  o->max = copyVec3(max);
  o->children = 0;
  o->at_max_depth = 0;

  /*
  printf("creating octree %.3lf,%.3lf,%.3lf ... %.3lf,%.3lf,%.3lf\n", 
	 o->min[2], o->min[1], o->min[0],
	 o->max[2], o->max[1], o->max[0] );
  */

  return o;
}
void subpoint( Octree* o,int oc,Vec3 min, Vec3 max)
{
	pvex_nor *m_p1,*m_p2;
	POSITION pos1,pos2;
	for(pos1=o->vex.GetHeadPosition(),pos2=o->normal.GetHeadPosition();pos1!=NULL;)
	{
		//pos1=o->vex.FindIndex(i);//pos2=o->normal.FindIndex(i);
		m_p1=(pvex_nor*)o->vex.GetNext(pos1);m_p2=(pvex_nor*)o->normal.GetNext(pos2);
		if((m_p1->x>min[0]&&m_p1->x<max[0])&&(m_p1->y>min[1]&&m_p1->y<max[1])
			&&(m_p1->z>min[2]&&m_p1->z<max[2]))
		{
			o->children[oc]->vex.AddHead(new pvex_nor(m_p1->x,m_p1->y,m_p1->z));
			o->children[oc]->normal.AddHead(new pvex_nor(m_p2->x,m_p2->y,m_p2->z));
		}
	}

}
void split_octree( Octree* o )
{

  double oc_min[3];
  double oc_max[3];
  
  Vec3 mid = makeVec3( (o->min[0] + o->max[0]) * 0.5, 
		       (o->min[1] + o->max[1]) * 0.5, 
		       (o->min[2] + o->max[2]) * 0.5 );
  int xp, yp, zp;
  int oc = 0;

  //o->children = (Octree**) malloc( 8 * sizeof(Octree*));
  o->children = new Octree*;
  for(zp=0; zp < 2; zp++)
  {
    if(zp == 0)
    { 
      oc_min[2] = o->min[2];
      oc_max[2] = mid[2];
    }
    else
    {
      oc_min[2] = mid[2];
      oc_max[2] = o->max[2];
    }
    
    for(yp=0; yp < 2; yp++)
    {
      if(yp == 0)
      { 
	oc_min[1] = o->min[1];
	oc_max[1] = mid[1];
      }
      else
      {
	oc_min[1] = mid[1];
	oc_max[1] = o->max[1];
      }
      
      for(xp=0; xp < 2; xp++)
      {
	if(xp == 0)
	{ 
	  oc_min[0] = o->min[0];
	  oc_max[0] = mid[0];
	}
	else
	{
	  oc_min[0] = mid[0];
	  oc_max[0] = o->max[0];
	}
	
	o->children[ (zp*4) + (yp*2) + xp ] = make_octree( oc_min, oc_max );
    subpoint( o,(zp*4) + (yp*2) + xp,oc_min,oc_max);
      }
    }
  }

}

/* ------------------------------------------------------------------------ */


int recursively_evaluate_octree( int min_depth, int max_depth, int current_depth, Octree* o )
{
  int deepest_child = current_depth;
  //int xp, yp, zp;
  int oc;
  int cd;
  //pvex_nor *m_p;
  //POSITION pos;
  /*Vec3 point = makeVec3(0,0,0);
  int p = 0;
  for(zp=0; zp < 2; zp++)
  {
    point[2] = (zp == 0) ? o->min[2] : o->max[2];

    for(yp=0; yp < 2; yp++)
    { 
      point[1] = (yp == 0) ? o->min[1] : o->max[1];

      for(xp=0; xp < 2; xp++)
      { 
	   point[0] = (xp == 0) ? o->min[0] : o->max[0];
	
	   o->value[ p++] = evaluate_point( point,o);
      }
    }
  }
  */
  o->density=evaluate1_point(o);

  current_depth++;
  
  o->not_fully_divided = (char) octree_needs_to_be_split( o );

  if( current_depth <= max_depth )
  {
    if(( current_depth <= min_depth) || ( o->not_fully_divided ))
    {
      //if(deepest_child==current_depth||deepest_child==0)
	  //{
	      
	      split_octree( o );
	 // }
      for(oc = 0; oc < 8; oc++)
      {
	     /*Vex.RemoveAll();
         for( pos = vex[oc].GetHeadPosition(); pos != NULL; )
		 {
                m_p=(pvex_nor*)vex[oc].GetNext( pos );
	        	Vex.AddHead(new pvex_nor(m_p->x,m_p->y,m_p->z));
		 }*/
         //if(deepest_child==current_depth||deepest_child==0)
		 //{
	     cd = recursively_evaluate_octree( min_depth, max_depth, current_depth, o->children[ oc ] );
		 //}
    
     	if(cd > deepest_child)
	    deepest_child = cd;
      }
    }
  }
  else
  {
    o->at_max_depth = 1;
  }
  
  return deepest_child;
}

/* ------------------------------------------------------------------------ */

int subdivide_octree( int min_depth, int max_depth, Octree* o )
{
  return recursively_evaluate_octree(min_depth, max_depth, 0, o );
}

double demo1( Vec3 pos )
{
  /* demo 1: the surface is a sphere of radius 0.75 centered at 0,0,0 
     
     function returns 1.0 if point inside sphere, 0.0 otherwise 
  */
  
  double dist_sq = (pos[0] * pos[0]) + (pos[1] * pos[1]) + (pos[2] * pos[2]);
  
  return ( dist_sq < 0.5625 ) ? 1.0 : 0.0;
}

double demo2( Vec3 pos )
{
  /* demo 2: the surface is two spheres, 
           A: radius 0.25 centered at -.25,-.5,0 
       and B: radius 0.5  centered at -0.5,0,0 

     function returns 1.0 if point inside sphere A, 2.0 for sphere B, 0.0 for neither
  */

  double dist_sq_a = ((pos[0]+.25) * (pos[0]+.25)) + ((pos[1]+.5) * (pos[1]+.5)) + (pos[2] * pos[2]);
  double dist_sq_b = ((pos[0]+.8) * (pos[0]+.8)) + (pos[1] * pos[1]) + (pos[2] * pos[2]);

  if( dist_sq_a <= .0625 )
    return 1.0;


  if( dist_sq_b <= .25 )
    return 2.0;

  return 0.0;
}

double demo3( Vec3 pos )
{
  /* demo 3: the surface is tiny sphere, radius 0.1 centered at -.5,.5,0 

     function returns 1.0 if point inside sphere A, 0.0 otherwise
  */


  double dist_sq = ((pos[0]+.5) * (pos[0]+.5)) + ((pos[1]-.5) * (pos[1]-.5)) + (pos[2] * pos[2]);

  return ( dist_sq < 0.01 ) ? 1.0 : 0.0;

}

double demo4( Vec3 pos )
{
  /* demo 4: wavey surface
     
     function returns 1.0 if point 'near' surface , 0.0 otherwise
     
  */
  double surface_height = sin( (pos[0] * 3.0) ) * cos ( (pos[1] * 3.0) );
  
  double distance_sq = (pos[2] - surface_height) * (pos[2] - surface_height);

  return ( distance_sq < 0.01 ) ? 1.0 : 0.0;
}

double demo5( Vec3 pos )
{
  /* demo 5: hemisphere, center 0,0,0 radius 0.5, cut by plane at z=0

   */

  double abs_dist_sq = ((pos[0]) * (pos[0])) + ((pos[1]) * (pos[1])) + (pos[2] * pos[2]);
  
  double surf_dist_sq = abs_dist_sq - 0.5625;
  if(surf_dist_sq  < 0)
    surf_dist_sq = -surf_dist_sq;

  if( (pos[2] > 0) && (surf_dist_sq < 0.1 ))
    return 1.0;
  else
    return .0;
}

double demo6( Vec3 pos )
{
  /* demo 6: another wavey surface
     
     function returns 1.0 if point 'near' surface , 0.0 otherwise
     
  */
  double surface_height = sin( (pos[0] * 2.0) ) + sin ( (pos[1] * 2.0) );
  
  double distance_sq = (pos[2] - surface_height) * (pos[2] - surface_height);

  return ( distance_sq < 0.01 ) ? 1.0 : 0.0;
}


double demo7( Vec3 pos )
{
  /* demo 7: a cylinder
     
     function returns 1.0 if point 'near' surface , 0.0 otherwise