Source

Pleasant3D_OpenSource / ToolPlugins / Slice / sliceVertices.cl


//
//  sliceVertices.cl
//  Slice
//
//  Created by Eberhard Rensch on 12.08.09.
//  Copyright 2009 Pleasant Software. All rights reserved.
//

#include "clSliceTypes.h"
#include "clEdgeIndexTypes.h"

typedef struct tagSTLFacetDiscreet {
	float	normalx;
	float	normaly;
	float	normalz;
	float	p0x;
	float	p0y;
	float	p0z;
	float	p1x;
	float	p1y;
	float	p1z;
	float	p2x;
	float	p2y;
	float	p2z;
	ushort	attrib;
}	STLFacetDiscreet;

const size_t kSlicedEdgeSize=sizeof(SlicedEdge);
const size_t kSTLFacetSize=12*sizeof(float)+2;//sizeof(ushort); // STLData is not aligned!
const size_t kLineOffsetSize=sizeof(LineOffset);

const float k2Pi = (2.f*(float)M_PI);

// This function calculates the needed memory for slicing a specific facet
__kernel void calcLineOffsets(__global LineOffset *lineOffsets, __constant const char *facets, __constant const float zStart, __constant const float layerHeight)
{		
	const uint facetIndex = get_global_id(0);
	STLFacetDiscreet* facet = (STLFacetDiscreet*)(facets+facetIndex*kSTLFacetSize);
		
	float minZ = facet->p0z;
	float maxZ = facet->p0z;	
	minZ=min(minZ, facet->p1z);
	maxZ=max(maxZ, facet->p1z);
	minZ=min(minZ, facet->p2z);
	maxZ=max(maxZ, facet->p2z);
	
//	// Calculate the min layer and the max layer affecting by this facet
	lineOffsets[facetIndex].minLayer = (uint)((minZ-zStart)/layerHeight);
	lineOffsets[facetIndex].maxLayer = max((uint)((maxZ-zStart)/layerHeight), lineOffsets[facetIndex].minLayer)+1;
	
	// Calculate the number of layers affecting this layer
	// This number of sliced lines needs to be reserved for this facet
	// Slice the facet. Only try slicing between the min/max layers calculated in calcLineOffsets
	
	lineOffsets[facetIndex].offset = 0;
	for(uint layerIndex = lineOffsets[facetIndex].minLayer; layerIndex<lineOffsets[facetIndex].maxLayer; layerIndex++)
	{
		float facetPX, facetPY, facetPZ;
		float facetQX, facetQY, facetQZ;
		
		float z=zStart+(float)layerIndex*layerHeight+layerHeight/2.; // The height for this slice
		uint edge1;
		for(edge1 = 0; edge1<3; edge1++)
		{
			switch(edge1)
			{
				case 0: facetPX=facet->p0x; facetPY=facet->p0y; facetPZ=facet->p0z;
						facetQX=facet->p1x; facetQY=facet->p1y; facetQZ=facet->p1z;
						break;
				case 1: facetPX=facet->p0x; facetPY=facet->p0y; facetPZ=facet->p0z;
						facetQX=facet->p2x; facetQY=facet->p2y; facetQZ=facet->p2z;
						break;
				case 2: facetPX=facet->p1x; facetPY=facet->p1y; facetPZ=facet->p1z;
						facetQX=facet->p2x; facetQY=facet->p2y; facetQZ=facet->p2z;
						break;
			}
			
			if((facetPZ<=z && facetQZ>z) || (facetPZ>z && facetQZ<=z)) // This edge is cut by the layer
				break;
		}

		if(edge1<2) // There's still a chance for a second edge!
		{
			for(uint edge2 = edge1+1; edge2<3; edge2++)
			{
				switch(edge2)
				{
					case 0: facetPX=facet->p0x; facetPY=facet->p0y; facetPZ=facet->p0z;
							facetQX=facet->p1x; facetQY=facet->p1y; facetQZ=facet->p1z;
							break;
					case 1: facetPX=facet->p0x; facetPY=facet->p0y; facetPZ=facet->p0z;
							facetQX=facet->p2x; facetQY=facet->p2y; facetQZ=facet->p2z;
							break;
					case 2: facetPX=facet->p1x; facetPY=facet->p1y; facetPZ=facet->p1z;
							facetQX=facet->p2x; facetQY=facet->p2y; facetQZ=facet->p2z;
							break;
				}
				if((facetPZ<=z && facetQZ>z) || (facetPZ>z && facetQZ<=z)) // This edge is cut by the layer, we have a winner!
				{
					lineOffsets[facetIndex].offset++;
					break;
				}
			}
		}
	}
}

__kernel void sliceTriangles(__global SlicedEdge* slicedEdges, __global SlicedEdgeLayerReference* slicedEdgeLayerReference, __constant const char *facets, __constant const LineOffset* lineOffsets, __constant const float zStart, __constant const float layerHeight)
{
	const unsigned int facetIndex = get_global_id(0);
	STLFacetDiscreet* facet = (STLFacetDiscreet*)(facets+facetIndex*kSTLFacetSize);
	
	if(lineOffsets[facetIndex].offset != (unsigned int)-1) // else there's nothing to slice!
	{
		// the start index for the sliced lines in the slicedEdges array for this line is calculated in the offset field
		unsigned int slicedEdgeIndex = lineOffsets[facetIndex].offset;

		// Slice the facet. Only try slicing between the min/max layers calculated in calcLineOffsets
		for(unsigned int layerIndex = lineOffsets[facetIndex].minLayer; layerIndex<lineOffsets[facetIndex].maxLayer; layerIndex++)
		{
			float facetPX, facetPY, facetPZ;
			float facetQX, facetQY, facetQZ;
			float z=zStart+(float)layerIndex*layerHeight+layerHeight/2.; // The height for this slice
			
			uint	edge1, 
					edge2=3;
			float2 edgeStart;
			
			// Look for an edge which is cut by the current z layer line
			for(edge1 = 0; edge1<3; edge1++)
			{
				switch(edge1)
				{
					case 0: facetPX=facet->p0x; facetPY=facet->p0y; facetPZ=facet->p0z;
							facetQX=facet->p1x; facetQY=facet->p1y; facetQZ=facet->p1z;
							break;
					case 1: facetPX=facet->p0x; facetPY=facet->p0y; facetPZ=facet->p0z;
							facetQX=facet->p2x; facetQY=facet->p2y; facetQZ=facet->p2z;
							break;
					case 2: facetPX=facet->p1x; facetPY=facet->p1y; facetPZ=facet->p1z;
							facetQX=facet->p2x; facetQY=facet->p2y; facetQZ=facet->p2z;
							break;
				}
				
				if((facetPZ<=z && facetQZ>z) || (facetPZ>z && facetQZ<=z)) // This edge is cut by the layer
				{
					float2 p = (float2)(facetPX, facetPY);
					float2 q = (float2)(facetQX, facetQY);
					float pZz = z-facetPZ;
					float faktor = pZz/(facetQZ-facetPZ);
					
					edgeStart=p+(q-p)*faktor;
					break;
				}
			}
			
			if(edge1<2) // There's still a chance for a second edge!
			{
				// Look for the second edge edge which is cut by the current z layer line
				for(edge2 = edge1+1; edge2<3; edge2++)
				{
					switch(edge2)
					{
						case 0: facetPX=facet->p0x; facetPY=facet->p0y; facetPZ=facet->p0z;
								facetQX=facet->p1x; facetQY=facet->p1y; facetQZ=facet->p1z;
								break;
						case 1: facetPX=facet->p0x; facetPY=facet->p0y; facetPZ=facet->p0z;
								facetQX=facet->p2x; facetQY=facet->p2y; facetQZ=facet->p2z;
								break;
						case 2: facetPX=facet->p1x; facetPY=facet->p1y; facetPZ=facet->p1z;
								facetQX=facet->p2x; facetQY=facet->p2y; facetQZ=facet->p2z;
								break;
					}
					if((facetPZ<=z && facetQZ>z) || (facetPZ>z && facetQZ<=z)) // This edge is cut by the layer
					{
						float2 p = (float2)(facetPX, facetPY);
						float2 q = (float2)(facetQX, facetQY);
						float pZz = z-facetPZ;
						float faktor = pZz/(facetQZ-facetPZ);
						
						// Payload
						slicedEdges[slicedEdgeIndex].endPoint = p+(q-p)*faktor;
						slicedEdges[slicedEdgeIndex].startPoint = edgeStart;

						// Meta-Data
						slicedEdges[slicedEdgeIndex].parentFacet = facetIndex;				
						slicedEdges[slicedEdgeIndex].parentStartEdge = edge1;				
						slicedEdges[slicedEdgeIndex].parentEndEdge = edge2;	
						slicedEdges[slicedEdgeIndex].normal = normalize((float2)(facet->normalx,facet->normaly));	
						slicedEdges[slicedEdgeIndex].connectsTo = (uint)-1;
						slicedEdges[slicedEdgeIndex].layerIndex = layerIndex; // This also validates the edge
						slicedEdges[slicedEdgeIndex].optimizedIndex = (uint)-1;

						slicedEdgeLayerReference[slicedEdgeIndex].layerIndex = layerIndex;
						slicedEdgeLayerReference[slicedEdgeIndex].slicedEdgeIndex = slicedEdgeIndex;
						slicedEdgeIndex++;
						break;						
					}
				}
			}
		}
	}
}

// The optimization needs 2 steps (2 separate kernels). Otherwise we might run into problems with parallel processing
__kernel void optimizeCornerPoints(__global uint* outOptimizedConnections, __global SlicedEdge *inSlicedEdges)
{
	const uint edgeIdex = get_global_id(0);
	SlicedEdge thisEdge = inSlicedEdges[edgeIdex];
	outOptimizedConnections[edgeIdex]=(uint)-1;
	if(thisEdge.connectsTo!=(uint)-1)
	{
		SlicedEdge nextEdge= inSlicedEdges[thisEdge.connectsTo];
		if(nextEdge.connectsTo!=(uint)-1)
		{
			float m1 = atan2(thisEdge.endPoint.y-thisEdge.startPoint.y, thisEdge.endPoint.x-thisEdge.startPoint.x);
			float m2 = atan2(nextEdge.endPoint.y-nextEdge.startPoint.y, nextEdge.endPoint.x-nextEdge.startPoint.x);
			float diff = fabs(m1-m2);
			int candidate = (diff<FLT_EPSILON);
			if(candidate) // Same Slope -> The next edge is unnecessary
			{
				outOptimizedConnections[edgeIdex]=nextEdge.connectsTo;
				//printf("Kandidat: %d\n", edgeIdex);
			}
		}
	}
}

__kernel void optimizeConnections(__global uint* outOptimizedConnections, __global SlicedEdge *inOutSlicedEdges)
{
	const uint edgeIdex = get_global_id(0);
//	printf("edgeIndex=%d, inOutSlicedEdges[edgeIdex].connectsTo=%d outOptimizedConnections[edgeIdex]=%d\n", edgeIdex, inOutSlicedEdges[edgeIdex].connectsTo, outOptimizedConnections[edgeIdex]);
//	__global SlicedEdge* thisEdge = &inOutSlicedEdges[edgeIdex];
//	if(outOptimizedConnections[edgeIdex]!=(uint)-1 && thisEdge->connectsTo != outOptimizedConnections[edgeIdex])
//	{
//		if(thisEdge->connectsTo!=(uint)-1)
//		{
//			__global SlicedEdge* nextEdge= &inOutSlicedEdges[thisEdge->connectsTo];
//			nextEdge->connectsTo = (uint)-1;
//			thisEdge->endPoint=nextEdge->endPoint;
//		}
//		thisEdge->connectsTo = outOptimizedConnections[edgeIdex];
//	}
}

__kernel void insetLoop(__global InsetLoopCorner* outInsetLoopCorners, __constant SlicedEdge *inSlicedEdges, __constant const float inset)
{
	const uint edgeIdex = get_global_id(0);
//	printf("sizeof(InsetLoopCorner)=%d, sizeof(SlicedEdge)=%d\n", sizeof(InsetLoopCorner), sizeof(SlicedEdge));
	
	SlicedEdge thisEdge = inSlicedEdges[edgeIdex];
		
	if(thisEdge.connectsTo!=(uint)-1 && thisEdge.optimizedIndex!=(uint)-1)
	{
//		printf("Handling Edge #%d Layer %d Facet %d (%f|%f)[E%d] -> (%f|%f)[E%d] N(%f|%f) ConnectsTo:%d LayerIndex:%d OptimizedIndex:%d\n",
//		edgeIdex, thisEdge.layerIndex, thisEdge.parentFacet, thisEdge.startPoint.x,thisEdge.startPoint.y, thisEdge.parentStartEdge, thisEdge.endPoint.x, thisEdge.endPoint.y, thisEdge.parentEndEdge, thisEdge.normal.x, thisEdge.normal.y, thisEdge.connectsTo, thisEdge.optimizedIndex);
		SlicedEdge nextEdge= inSlicedEdges[thisEdge.connectsTo];

//		float2 g1v = thisEdge.normal*inset;
//		float2 g1p1 = thisEdge.startPoint-g1v;
//		float2 g1p2 = thisEdge.endPoint-g1v;
//				
//		float2 g2v = nextEdge.normal*inset;
//		float2 g2p1 = nextEdge.startPoint-g2v;
//		float2 g2p2 = nextEdge.endPoint-g2v;
//		
//		// Slope of g1
//		float m1 = INFINITY;
//		if(g1p2.x!=g1p1.x)
//			m1 = (g1p2.y-g1p1.y)/(g1p2.x-g1p1.x);
//			
//		// Slope of g2
//		float m2 = INFINITY;
//		if(g2p2.x!=g2p1.x)
//			m2 = (g2p2.y-g2p1.y)/(g2p2.x-g2p1.x);
//
//		if(m1==m2)
//		{
//			// TODO: This point is not necessary and should be deleted from the loop!
//			outInsetLoopCorners[edgeIdex].point = g1p2;
//			outInsetLoopCorners[edgeIdex].sort = 0;
//		}
//		else if(m1==INFINITY)
//		{
//			float t2 = g2p1.y-m2*g2p1.x;
//			outInsetLoopCorners[edgeIdex].point = (float2)(g1p2.x, m2*g1p2.x+t2);
//			outInsetLoopCorners[edgeIdex].sort = 1;
//		}
//		else if(m2==INFINITY)
//		{
//			float t1 = g1p1.y-m1*g1p1.x;
//			outInsetLoopCorners[edgeIdex].point = (float2)(g2p2.x, m1*g2p2.x+t1);
//			outInsetLoopCorners[edgeIdex].sort = 2;
//		}	
//        else
//		{
//			float t1 = g1p1.y-m1*g1p1.x;
//			float t2 = g2p1.y-m2*g2p1.x;
//			
//			outInsetLoopCorners[edgeIdex].point = (float2)((t2-t1)/(m1-m2), (m1*t2-m2*t1)/(m1-m2));
//			outInsetLoopCorners[edgeIdex].sort = 3;
//		}
//		
		// Calculate the cornerpoint's normal
		
		float rad1 = atan2(thisEdge.normal.y, thisEdge.normal.x);
		float rad2 = atan2(nextEdge.normal.y, nextEdge.normal.x);
		
		float drad1=rad1-rad2;
		if(drad1<0.) drad1+=k2Pi;
		
		float drad2=rad2-rad1;
		if(drad2<0.) drad2+=k2Pi;		
		
		float prad;
		if(fabs(drad1)>fabs(drad2))
			prad = rad1+drad2/2.f;
		else
			prad = rad2+drad1/2.f;
		
		outInsetLoopCorners[thisEdge.optimizedIndex].normal = (float2)(cos(prad), sin(prad));		
		outInsetLoopCorners[thisEdge.optimizedIndex].point = thisEdge.endPoint;//-outInsetLoopCorners[edgeIdex].normal*inset;
	}	
}