arbritary slicing of a volume...

Hello,

I would like to compute an arbritary slice through a volume.

Say the volume is 64x64x64, and that I have the plane equation & normal for a user oriented plane.

Hence I would like to find the 2d image this plane represents/slices through my volume.

Here is what i currently do:

I traverse the whole volume and check wether the current index is infront or behind my plane.

This is identified by a flag, 1 for infront, 2 for behind.

Then i check if the current point’s flag is the same as and it’s +/- y neighbours and +/- x neighbours.

If the current point is different from any of it’s neighbours, then i know it is part of the slice.

With this i can identify which indices in my volume are part of the slice.

  1. But i can’t figure out how to move them from just flagged indices of my volume to indices of a 2d slice. Im stuck :frowning:

  2. This seems a bit brute force, maybe even incorrect, so are there better ways to do this?

I look forward to suggestions, and take it slowly with the maths, clearly I’m just learning :wink:

Thanks,

UT.

hi, how about checking out the thread
“How to draw practically any poly from verts???”

ignore my above post.

why don’t you just interpolate the slice plane using the 4 points found on the outer edge of the volume (the points formed at the intersecting planes of your cube), this does lead to the question, at what resolution do i use. I’m assuming you are using an array store to hold your data however.

Hello again,

Well I found a somewhat cough ugly cough solution to my problem…

Here is what i do:

Say I have the following,

A volume which is 256x256x256.
An XY axis aligned plane.
A 2d image array.

What I do is allow the user to fully orient the plane, and keep the rotation matix. This plane is the arbritary slicing plane.

So now what we want to do is retrieve the image this plane represents through the volume.

So I find the inverse of the rotation matix, via a transpose.

Then I traverse my 2d image array, and treat the current index as a point (x_i,y_i,0).

I then multiply the point by my inverted rotation matrix and this gives me the position of my index in the volume.

Lastly I just look for the closest voxel to the transformed index, and fill in the image array.

Long, ugly and probably better ways to do this…
Any suggestions on how to improve are always welcome.

UT.

It would propably easier if you use a 3D texture. As a special extra, you’ll get interpolation for free :wink:

Hello,

yes I tried this too… I actually use 3d texturing to render the volume.

What I was doing is,

enable 3d texturing,
draw the single user oriented plane,
disable 3d texturing.

And as i had the auto gen of texture coords, this would indeed display correctly the image of interest.

But what i need is to be able to read back that image back to the system. To my understanding the only way i found to do this was via a render to texture. Which meant I had to re-project the whole scene so that I get only the image of interest, and not the whole scene…

Is there any other way to get a slice of a 3d texture back to the system?

thanks,

UT.

if your images are large, and interpolation is good for you, then the hardware way is a nice one. maybe you could use the texure matrix stack for slicing the vollume, and render just a simple screen aligned quad with glOrtho projection. rendering to a texture with for example FBO’s is quite fast and lacks only antialiasing, but no prob when mipmapping is used. So just do it this simple way.

if your vollume gets smaller, and no interpolation is needed, then probably cpu work is better, this is mainly because of the relatively slow download of rendered data from grafics board to ram.

One efficient algorithmus used to draw lines, planes etc. or map them with textures is called “Bresenham”, or generalisations thereof. (which i think are also used inside the GPU for all kinds of work) I googled around but found no actual code-near description… mabye you could do better.
(Something like 3d texturing bresenham)

hope this helps
Paul

Hi UT666 , all

hope this helps

This is a simple object slicer CUBE’s, Pyramids etc, it works ONLY in triangle based surfaces and NON complex objects, but is fairly easy to convert.

I pump in for pobj->Vpos the 3dtexture coords of a cube and directly render the returned poly coords to both the Vertices (Scaled) and the 3dTexCoords without the requirement for autogen.

How to use

pobj->Vpos[]=Vertices for object (3d cube texture coords)

pobj->Indice_Pool[surfaceNo*3 + surface indice order (only triangle surfaces)]=pobj->Vpos[^data]

Bind object to the slicer
Bind_Volume(^object)

loop

Set the Slicer Normal
	Set_Volume_Slice_Normal(_array3 Vnrm)
 
Extract Polygon
	Extract_Cap(GLfloat plane_dist) (Between Max/Min distances)

till finished with obj

Free Slicer
Free_Slicer(true)

////C++ Code
</font><blockquote><font size=“1” face=“Verdana, Arial”>code:</font><hr /><pre style=“font-size:x-small; font-family: monospace;”>
#define FreePtr(pt)
if (pt!=NULL) free(pt);
pt=NULL;

typedef GLfloat _array3[3];

typedef struct tag_SurfaceCut
{
int IndiceLinks[2][2];
bool used;
} SurfaceCut;

template<typename T>T* T_MemAlloc(unsigned int scount)
{
T *Temp;
Temp=(T *) malloc(sizeof(T)*scount);
if (Temp==NULL) return NULL;
memset(Temp,0,sizeof(T)*scount); // zero mem
return Temp;
}

class Volume_Slicer
{
public:

	Volume_Slicer()
		{
		pobj=NULL;
		Vertex_Distance=NULL;
		Poly_Cap=NULL;
		Edge_Links=NULL;
		EdgeToEdge=NULL;
		}
	~Volume_Slicer()
		{
		Free_Slicer(true);
		}

	object_pool *pobj;  //bind object

	_array3 SliceNormal;           //Surface Cut Plane Normal
	GLfloat *Vertex_Distance;      //ASSIGNED T_Mem
	GLfloat Max_Dist,Min_Dist;     //Returned Max/Min Vertex Distance to Surface Plane, based on Surface Normal

	_array3 *Poly_Cap;             //Returned Drawing Polygon Slice
	int			No_PolyVertices;       //No of Vertices in Polygon

	SurfaceCut *Edge_Links; //pointer to indice
	int No_EdgeLinks;       //etc

	int *EdgeToEdge;        //Edges found for cutting
	int NoEdgeToEdge;       //etc

	void Free_Slicer(bool finished)  //TRUE FINISHED WITH ALL DATA, FALSE STILL USING POBJ
		{
		if(finished)
			{
			FreePtr(Vertex_Distance);
			FreePtr(Edge_Links);
			pobj=NULL;
			}
		FreePtr(Poly_Cap);
		FreePtr(EdgeToEdge);
		}

	bool Bind_Volume(object_pool *pobjp) //Bind object to Slicer (Only using VERTS/SURFACE INDICES)
		{
		if (pobjp==NULL) return false;

		if (pobj!=NULL) Free_Slicer(true);	
		pobj=pobjp; 

		//create MaxPoly(vertices)buffer

		Vertex_Distance=T_MemAlloc&lt;GLfloat&gt;(pobj-&gt;NoVertices);
		if (Vertex_Distance==NULL) 
			{
			pobj=NULL;				
			return false;
			}

		Edge_Links=T_MemAlloc&lt;SurfaceCut&gt;(pobj-&gt;NoSurfaces);
		if (Edge_Links==NULL) 
			{
			pobj=NULL;				
			free(Vertex_Distance);
			return false;
			}

		return true;
		} //error invalid pobjp, not enough mem verts or MaxPoly

	int Set_Volume_Slice_Normal(_array3 Vnrm)  //Set Slice Normal and dot Vertices with normal 
		{
		if (pobj==NULL) return 0;

		Vertex(SliceNormal,Vnrm);      

		Max_Dist=-3.4E38;
		Min_Dist=3.4E38;
		for (int i=0;i&lt;pobj-&gt;NoVertices;i++)				
			{			
			Vertex_Distance[i]=DotProduct(pobj-&gt;Vpos[i],SliceNormal);
			Max_Dist=get_max(Max_Dist,Vertex_Distance[i]);
			Min_Dist=get_min(Min_Dist,Vertex_Distance[i]);
			}			
		return i; //error i=0, no data or invalid pointer pobj
		}


	bool Extract_Cap(GLfloat plane_dist) //Creates Polygon with Min&lt;=plane_dist&lt;=Max
		{
		if (pobj==NULL) return false;
		Free_Slicer(false);
		
		GLfloat PlaneEq;
	
		//compile edges cut

		int si,vi,indxA,indxB;

		int rcEL,rcSC;
		rcEL=0;

		for (si=0;si&lt;pobj-&gt;NoSurfaces;si++)
			{
			rcSC=0;
			for(vi=0;vi&lt;3;vi++)
				{
				indxA=pobj-&gt;Indice_Pool[si*3+vi];         
				indxB=pobj-&gt;Indice_Pool[si*3+((vi+1)%3)];

				PlaneEq=(Vertex_Distance[indxA]-plane_dist)*(Vertex_Distance[indxB]-plane_dist);
				if (PlaneEq&lt;=0) //Plane Cuts
					{
					Edge_Links[rcEL].IndiceLinks[rcSC][0]=indxA;  //store vertices for edge
					Edge_Links[rcEL].IndiceLinks[rcSC][1]=indxB;
					Edge_Links[rcEL].used=false;
					rcSC++;
					}
				}
			if (rcSC&gt;0) rcEL++;
			}
		
		//create linked list of indices then compute actual cuts
		NoEdgeToEdge=rcEL*2;
		EdgeToEdge=T_MemAlloc&lt;int&gt;(NoEdgeToEdge);
		if (EdgeToEdge==NULL) return false;
		EdgeToEdge[0]=0;  //n*2 + which surface edge
		EdgeToEdge[1]=1;

		int inspecting,cmpA,cmpB,cmpA2,cmpB2;
  
		inspecting=1; //rcEL 0 rcSC 1
		cmpA=Edge_Links[0].IndiceLinks[1][0];
		cmpB=Edge_Links[0].IndiceLinks[1][1];
		Edge_Links[0].used=true;

		int usedfound=1,chksum=0,EtoE=2;
		bool cmpfound;
		si=0;
		while ((usedfound&lt;rcEL)&&(chksum&lt;rcEL))
			{
			if (Edge_Links[si].used==false)
				{					
				cmpfound=false;
				for(vi=0;vi&lt;2;vi++)
					{
					cmpA2=Edge_Links[si].IndiceLinks[vi][0];
					cmpB2=Edge_Links[si].IndiceLinks[vi][1];
					if (((cmpA==cmpA2)&&(cmpB==cmpB2))

Dont know why this is truncated but heres the rest from si=

</font><blockquote><font size=“1” face=“Verdana, Arial”>code:</font><hr /><pre style=“font-size:x-small; font-family: monospace;”>

si=0;
		while ((usedfound&lt;rcEL)&&(chksum&lt;rcEL))
			{
			if (Edge_Links[si].used==false)
				{					
				cmpfound=false;
				for(vi=0;vi&lt;2;vi++)
					{
					cmpA2=Edge_Links[si].IndiceLinks[vi][0];
					cmpB2=Edge_Links[si].IndiceLinks[vi][1];
					if (((cmpA==cmpA2)&&(cmpB==cmpB2))

ARRGHH , it appears to be ‘or’ truncating

heres it is as a txt file

http://myweb.tiscali.co.uk/reapers/slicercode.txt

if there is any missing info let me know, thx

Hope this helps, any question …

regards

ps, would any one please comment on my code as i’ve never had any feed back, thx in advance.