/* * fractalobject.c * * Peter Janssen * */ #include "geom.h" #include "fractalobject.h" #include "bounds.h" #include #include #include FractalClosure *Fractalclosure; PointPool *Pointpool; TrianglePool *Trianglepool; EntityPool *Entitypool; static Methods *iFractalObjMethods = NULL; static char fractalobjName[] = "fractalobject"; static unsigned long raynumber = 1; static FirstPoint *SideList; unsigned long FractalObjTests, FractalObjHits; static Float Scale; /* * Random Generator */ #define rand_m (unsigned long)2147483647 #define rand_q (unsigned long)127773 #define rand_a (unsigned int)16807 #define rand_r (unsigned int)2836 /* * F(z) = (az)%m * = az-m(az/m) * * F(z) = G(z)+mT(z) * G(z) = a(z%q)- r(z/q) * T(z) = (z/q) - (az/m) * * F(z) = a(z%q)- rz/q+ m((z/q) - a(z/m)) * = a(z%q)- rz/q+ m(z/q) - az */ static long rand_seed; static Float rand_number() { long lo, hi, test; static Float inv_rand_m; inv_rand_m = 1.0 / rand_m; hi = rand_seed/rand_q; lo = rand_seed%rand_q; test = rand_a*lo - rand_r*hi; if (test > 0) rand_seed = test; else rand_seed = test + rand_m; return (rand_seed * inv_rand_m * 2) - 1; } /* * Calculate the boundingbox and the size of a voxel * for every subobject in the fractalobject */ static void CalculateBoundingBoxesAndVoxelSizes(fractalobj) FractalObj *fractalobj; { int i, j, CurrPointNumber; SubObject *CurrSubobject; FractalPoint *Points; Points = Fractalclosure->Pointpool->Points; for(CurrSubobject = fractalobj->SubObjects; CurrSubobject; CurrSubobject = CurrSubobject->Next) { for (j = 0; j <= 2; j++) CurrSubobject->Bounds[LOW][j] = MAXVALUE; for (j = 0; j <= 2; j++) CurrSubobject->Bounds[HIGH][j] = -MAXVALUE; for (i = 0; i < CurrSubobject->Trianglepool->NumberOfTriangles; i++) { for (j = 0; j <= 2; j++) { CurrPointNumber = CurrSubobject->Trianglepool->Triangles[i].Point[j]; CurrSubobject->Bounds[LOW][X] = min(CurrSubobject->Bounds[LOW][X], Points[CurrPointNumber].x); CurrSubobject->Bounds[HIGH][X] = max(CurrSubobject->Bounds[HIGH][X], Points[CurrPointNumber].x); CurrSubobject->Bounds[LOW][Y] = min(CurrSubobject->Bounds[LOW][Y], Points[CurrPointNumber].y); CurrSubobject->Bounds[HIGH][Y] = max(CurrSubobject->Bounds[HIGH][Y], Points[CurrPointNumber].y); CurrSubobject->Bounds[LOW][Z] = min(CurrSubobject->Bounds[LOW][Z], Points[CurrPointNumber].z); CurrSubobject->Bounds[HIGH][Z] = max(CurrSubobject->Bounds[HIGH][Z], Points[CurrPointNumber].z); } } CurrSubobject->VoxelSize[X] = (CurrSubobject->Bounds[HIGH][X] - CurrSubobject->Bounds[LOW][X]) / CurrSubobject->XSize; CurrSubobject->VoxelSize[Y] = (CurrSubobject->Bounds[HIGH][Y] - CurrSubobject->Bounds[LOW][Y]) / CurrSubobject->YSize; CurrSubobject->VoxelSize[Z] = (CurrSubobject->Bounds[HIGH][Z] - CurrSubobject->Bounds[LOW][Z]) / CurrSubobject->ZSize; } } /* * Calculate new points on the midpoints of the * three edges of the triangle */ static void CalculatePoints(fractalobj, NewPoints, FTriangle, Delta) FractalObj *fractalobj; int NewPoints[3]; FractalTriangle *FTriangle; DeltaInfo Delta[3]; { int i, j, fpoint, lpoint, NPoints, Pointsallocated; FirstPoint *Currfp, *chasingfp, *Newfp; LastPoint *Currlp, *chasinglp, *Newlp; FractalPoint *Points, *TPoints; PointPool *Pointpool; Float displacement, size1, size2; Vector n1, n2; Pointpool = Fractalclosure->Pointpool; NPoints = Pointpool->NumberOfPoints; Pointsallocated = Pointpool->PointsAllocated; Points = Pointpool->Points; for(i=0; i<=2; i++) { fpoint = min(FTriangle->Point[(i + 1)%3], FTriangle->Point[i]); lpoint = max(FTriangle->Point[(i + 1)%3], FTriangle->Point[i]); for(Currfp = SideList; fpoint > Currfp->Point; chasingfp = Currfp, Currfp = Currfp->Next); if (Currfp->Point == fpoint) { /* * sidelist has an entry with firstpoint */ chasinglp = NULL; for(Currlp = Currfp->Info; Currlp && (Currlp->Point < lpoint); chasinglp = Currlp, Currlp = Currlp->Next); if (Currlp && (Currlp->Point == lpoint)) { /* * sidelist has an entry for pointpair fpoint, lpoint */ NewPoints[i] = Currlp->NewPoint; continue; } } /* * Calculate New Point */ displacement = Delta[i].distance * pow((fractalobj->Excentricity * Delta[i].distance * Scale), fractalobj->Dimension); /* n1 */ if (Delta[i].z != 0) { n1.x = -Delta[i].z; n1.y = -Delta[i].z; n1.z = Delta[i].x + Delta[i].y; } else if (Delta[i].y != 0) { n1.x = -Delta[i].y; n1.y = Delta[i].x + Delta[i].z; n1.z = -Delta[i].y; } else { n1.x = Delta[i].y + Delta[i].z; n1.y = -Delta[i].x; n1.z = -Delta[i].x; } VecNormalize(&n1); /* n2 */ VecCross(&Delta[i], &n1, &n2); VecNormalize(&n2); size1 = rand_number(); size2 = rand_number(); CheckPoints(Points, NPoints, Pointsallocated); Points[NPoints].x = (Points[fpoint].x + Points[lpoint].x) / 2 + displacement * size1 * n1.x + displacement * size2 * n2.x; Points[NPoints].y = (Points[fpoint].y + Points[lpoint].y) / 2 + displacement * size1 * n1.y + displacement * size2 * n2.y; Points[NPoints].z = (Points[fpoint].z + Points[lpoint].z) / 2 + displacement * size1 * n1.z + displacement * size2 * n2.z; NewPoints[i] = NPoints; /* * Currfp points to position equal to or just bigger than fpoint */ Newlp = (LastPoint *)share_malloc(sizeof(LastPoint)); Newlp->Point = lpoint; if (Currfp->Point == fpoint) { /* * Currlp points tot position equal to or just bigger than * lpoint or is equal to NULL * Insert between chasinglp and Currlp */ if (chasinglp == NULL) /* * the first in line */ Currfp->Info = Newlp; else chasinglp->Next = Newlp; Newlp->Next = Currlp; } else { Newfp = (FirstPoint *)share_malloc(sizeof(FirstPoint)); chasingfp->Next = Newfp; Newfp->Next = Currfp; Newfp->Info = Newlp; Newfp->Point = fpoint; Newlp->Next = NULL; } Newlp->NewPoint = NPoints++; } Pointpool->NumberOfPoints = NPoints; Pointpool->PointsAllocated = Pointsallocated; Pointpool->Points = Points; } /* * Take triangulated object and distort it */ static void BuildFractalObject(fractalobj) FractalObj *fractalobj; { Float Maxsize, Delta_x, Delta_y, Delta_z, distance; int Changed, Poolchanged, ONTriangles, NTriangles, Trianglesallocated; FractalPoint *FPoints; SubObject *CurrSubobject; TrianglePool *Trianglepool; FractalTriangle *FTriangles; DeltaInfo Delta[3]; int i, j, k, NewPoints[3], t1, t2; FirstPoint *Currfp, *chasingfp; LastPoint *Currlp, *chasinglp; Maxsize = fractalobj->MaximumSize; Changed = TRUE; Scale = 0; FPoints = Pointpool->Points; for (CurrSubobject = fractalobj->SubObjects; CurrSubobject; CurrSubobject = CurrSubobject->Next) { FTriangles = CurrSubobject->Trianglepool->Triangles; for (i=0; i < CurrSubobject->Trianglepool->NumberOfTriangles; i++) { for(j = 0; j <= 2; j++) { t1 = FTriangles[i].Point[j]; t2 = FTriangles[i].Point[(j+1)%3]; Delta_x = FPoints[t1].x - FPoints[t2].x; Delta_y = FPoints[t1].y - FPoints[t2].y; Delta_z = FPoints[t1].z - FPoints[t2].z; distance = sqrt(Delta_x * Delta_x + Delta_y * Delta_y + Delta_z * Delta_z); Scale = max(Scale, distance); } } } Scale = 1.0 / Scale; while (Changed) { Changed = FALSE; for (CurrSubobject = fractalobj->SubObjects; CurrSubobject; CurrSubobject = CurrSubobject->Next) { Trianglepool = CurrSubobject->Trianglepool; if (Trianglepool->PoolChanged) { Poolchanged = FALSE; ONTriangles = NTriangles = Trianglepool->NumberOfTriangles; Trianglesallocated = Trianglepool->TrianglesAllocated; FTriangles = Trianglepool->Triangles; for (i=0; i < ONTriangles; i++) { FPoints = Pointpool->Points; for(j = 0; j <= 2; j++) { t1 = FTriangles[i].Point[j]; t2 = FTriangles[i].Point[(j+1)%3]; Delta[j].x = FPoints[t1].x - FPoints[t2].x; Delta[j].y = FPoints[t1].y - FPoints[t2].y; Delta[j].z = FPoints[t1].z - FPoints[t2].z; Delta[j].distance = sqrt(Delta[j].x * Delta[j].x + Delta[j].y * Delta[j].y + Delta[j].z * Delta[j].z); } if ((Delta[0].distance > Maxsize) || (Delta[1].distance > Maxsize) || (Delta[2].distance > Maxsize)) { Poolchanged = Changed = TRUE; if (NTriangles + 3 >= Trianglesallocated) { Trianglesallocated += TrianglesToAllocate; FTriangles = (FractalTriangle *)realloc((char *)FTriangles, Trianglesallocated * sizeof(FractalTriangle)); } NTriangles += 3; CalculatePoints(fractalobj, NewPoints, &FTriangles[i], Delta); FTriangles[NTriangles - 1].Point[0] = NewPoints[0]; FTriangles[NTriangles - 1].Point[1] = NewPoints[1]; FTriangles[NTriangles - 1].Point[2] = NewPoints[2]; FTriangles[NTriangles - 3].Point[0] = NewPoints[0]; FTriangles[NTriangles - 3].Point[1] = FTriangles[i].Point[1]; FTriangles[NTriangles - 3].Point[2] = NewPoints[1]; FTriangles[NTriangles - 2].Point[0] = NewPoints[2]; FTriangles[NTriangles - 2].Point[1] = NewPoints[1]; FTriangles[NTriangles - 2].Point[2] = FTriangles[i].Point[2]; FTriangles[i].Point[1] = NewPoints[0]; FTriangles[i].Point[2] = NewPoints[2]; } } Trianglepool->PoolChanged = Poolchanged; Trianglepool->NumberOfTriangles = NTriangles; Trianglepool->TrianglesAllocated = Trianglesallocated; Trianglepool->Triangles = FTriangles; } } /* * release sidelist and setup new one */ for(Currfp = SideList->Next; Currfp->Point != MAXVALUE; chasingfp = Currfp, Currfp = Currfp->Next, free(chasingfp)) for(Currlp = Currfp->Info; Currlp; chasinglp = Currlp, Currlp = Currlp->Next, free(chasinglp)); SideList->Next = Currfp; } CalculateBoundingBoxesAndVoxelSizes(fractalobj); } /* * Convert worldcoordinates to gridcoordinates */ static void Pos2Grid(subobj, pos, index) SubObject *subobj; Float pos[3]; int index[3]; { index[X] = (int)(Fx2voxel(subobj, pos[0])); index[Y] = (int)(Fy2voxel(subobj, pos[1])); index[Z] = (int)(Fz2voxel(subobj, pos[2])); if (index[X] == subobj->XSize) index[X]--; if (index[Y] == subobj->YSize) index[Y]--; if (index[Z] == subobj->ZSize) index[Z]--; } /* * Put triangles of fractalobject in voxels * for efficient raytracing */ static void ConvertFractalObject(fractalobj) FractalObj *fractalobj; { int i, x, y, z, low[3], high[3]; Float bounds[2][3]; SubObject *CurrSubobject; FractalTriangle *Triangles, *CurrTriangle; TriangleInfo *Triangleinfo; TriangleClosure *Triangleclosure; FractalPoint *Points; int p0, p1, p2; Points = Fractalclosure->Pointpool->Points; for(CurrSubobject = fractalobj->SubObjects; CurrSubobject; CurrSubobject = CurrSubobject->Next) { Triangles = CurrSubobject->Trianglepool->Triangles; CurrSubobject->Cells = (TriangleClosure ****)share_malloc(CurrSubobject->XSize * sizeof(TriangleClosure ***)); for (x = 0; x < CurrSubobject->XSize; x++) { CurrSubobject->Cells[x] = (TriangleClosure ***)share_malloc(CurrSubobject->YSize * sizeof(TriangleClosure **)); for (y = 0; y < CurrSubobject->YSize; y++) CurrSubobject->Cells[x][y] = (TriangleClosure **)share_calloc( CurrSubobject->ZSize , sizeof(TriangleClosure *)); } for (i = 0; i < CurrSubobject->Trianglepool->NumberOfTriangles; i++) { CurrTriangle = &Triangles[i]; p0 = CurrTriangle->Point[0]; p1 = CurrTriangle->Point[1]; p2 = CurrTriangle->Point[2]; bounds[LOW][X] = min3(Points[p0].x, Points[p1].x, Points[p2].x); bounds[HIGH][X] = max3(Points[p0].x, Points[p1].x, Points[p2].x); bounds[LOW][Y] = min3(Points[p0].y, Points[p1].y, Points[p2].y); bounds[HIGH][Y] = max3(Points[p0].y, Points[p1].y, Points[p2].y); bounds[LOW][Z] = min3(Points[p0].z, Points[p1].z, Points[p2].z); bounds[HIGH][Z] = max3(Points[p0].z, Points[p1].z, Points[p2].z); Pos2Grid(CurrSubobject, bounds[LOW], low); Pos2Grid(CurrSubobject, bounds[HIGH], high); Triangleinfo = (TriangleInfo *)share_malloc(sizeof(TriangleInfo)); Triangleinfo->Point[0] = p0; Triangleinfo->Point[1] = p1; Triangleinfo->Point[2] = p2; Triangleinfo->Counter = 0; for (x = low[X]; x <= high[X]; x++) for (y = low[Y]; y <= high[Y]; y++) for (z = low[Z]; z <= high[Z]; z++) { Triangleclosure = (TriangleClosure *)share_malloc(sizeof(TriangleClosure)); Triangleclosure->Triangle = Triangleinfo; Triangleclosure->Next = CurrSubobject->Cells[x][y][z]; CurrSubobject->Cells[x][y][z] = Triangleclosure; } } /* * release trianglepool */ free(CurrSubobject->Trianglepool->Triangles); free(CurrSubobject->Trianglepool); } } /* * Create FractalObject and return reference to it */ FractalObj *FractalObjCreate(Number, Dimension, MaxSize, Excentricity, fc) Float MaxSize, Dimension, Excentricity; int Number; FractalClosure *fc; { FractalObj *fractalobj; int *TList, TriangleNumber; FractalTriangle *InitialTriangles, *LocalTriangles, *TTriangles; FractalEntity *CurrEntity; int NPoints, NTriangles, NEntities, i, j, LocalAllocated, LocalNTriangles; SubObject *Subobjects, *TSubobject; TriangleList *CurrTriangle, *chasingTriangle; int CurrPoint; fractalobj = (FractalObj *)share_malloc(sizeof(FractalObj)); fractalobj->MaximumSize = MaxSize; fractalobj->Dimension = Dimension; fractalobj->Excentricity = Excentricity; fractalobj->Closure = fc; rand_seed = Number; Fractalclosure = fc; Pointpool = Fractalclosure->Pointpool; Trianglepool = Fractalclosure->Trianglepool; Entitypool = Fractalclosure->Entitypool; NPoints = Pointpool->NumberOfPoints; NTriangles = Trianglepool->NumberOfTriangles; NEntities = Entitypool->NumberOfEntities; if (NPoints <= 0) RLerror(RL_ABORT, "Illegal number of Points : %d", NPoints); if (NTriangles <= 0) RLerror(RL_ABORT, "Illegal number of Triangles : %d", NTriangles); if (NEntities <= 0) RLerror(RL_ABORT, "Illegal number of Entities : %d", NEntities); TList = (int *)Malloc(sizeof(int) * NTriangles); Subobjects = NULL; InitialTriangles = Trianglepool->Triangles; for (i=0; i < NEntities; i++) { CurrEntity = &Entitypool->Entities[i]; TSubobject = (SubObject *)share_malloc(sizeof(SubObject)); TSubobject->Trianglepool = (TrianglePool *)share_malloc(sizeof(TrianglePool)); LocalTriangles = (FractalTriangle *)share_malloc(TrianglesToAllocate * sizeof(FractalTriangle)); LocalAllocated = TrianglesToAllocate; TSubobject->XSize = CurrEntity->XSize; TSubobject->YSize = CurrEntity->YSize; TSubobject->ZSize = CurrEntity->ZSize; LocalNTriangles = 0; for (CurrTriangle = CurrEntity->Triangles; CurrTriangle; chasingTriangle = CurrTriangle, CurrTriangle = CurrTriangle->Next, free(chasingTriangle)) { TriangleNumber = CurrTriangle->TriangleNumber; if (TriangleNumber < 0 || TriangleNumber >= NTriangles) RLerror(RL_ABORT, "Illegal Triangle number : %d\n", TriangleNumber+1); CheckTriangles(LocalTriangles, LocalNTriangles, LocalAllocated); CurrPoint = InitialTriangles[TriangleNumber].Point[0]; if ((CurrPoint < 0) || (CurrPoint >= NPoints)) RLerror(RL_ABORT, "Illegal Point number : %d\n", CurrPoint+1); LocalTriangles[LocalNTriangles].Point[0] = CurrPoint; CurrPoint = InitialTriangles[TriangleNumber].Point[1]; if ((CurrPoint < 0) || (CurrPoint >= NPoints)) RLerror(RL_ABORT, "Illegal Point number : %d\n", CurrPoint+1); LocalTriangles[LocalNTriangles].Point[1] = CurrPoint; CurrPoint = InitialTriangles[TriangleNumber].Point[2]; if ((CurrPoint < 0) || (CurrPoint >= NPoints)) RLerror(RL_ABORT, "Illegal Point number : %d\n", CurrPoint+1); LocalTriangles[LocalNTriangles].Point[2] = CurrPoint; LocalNTriangles++; TList[TriangleNumber] = TRUE; } TSubobject->Trianglepool->NumberOfTriangles = LocalNTriangles; TSubobject->Trianglepool->TrianglesAllocated = LocalAllocated; TSubobject->Trianglepool->Triangles = LocalTriangles; TSubobject->Trianglepool->PoolChanged = TRUE; TSubobject->Next = Subobjects; Subobjects = TSubobject; } fractalobj->SubObjects = Subobjects; /* * Check if every triangle is placed in a subobject */ for (i=0; i < NTriangles; i++) if (!TList[i]) RLerror(RL_ADVISE, "Triangle %d not placed in entity.\n", i+1); free(TList); /* * release triangles uit trianglepool */ free(Trianglepool->Triangles); free(Trianglepool); /* * release entitypool */ free(Entitypool->Entities); free(Entitypool); /* * Build empty sidelist */ SideList = (FirstPoint *)share_malloc(sizeof(FirstPoint)); SideList->Point = -1; SideList->Next = (FirstPoint *)share_malloc(sizeof(FirstPoint)); SideList->Next->Point = MAXVALUE; /* * Distort object */ BuildFractalObject(fractalobj); /* * free sidelist */ free(SideList->Next); free(SideList); /* * Convert to voxels */ ConvertFractalObject(fractalobj); return fractalobj; } /* * Intersect ray with fractalobject. */ int FractalObjIntersect(fractalobj, ray, mindist, maxdist) FractalObj *fractalobj; Ray *ray; Float mindist, *maxdist; { int hit; Float offset; Vector curpos; unsigned long counter; SubObject *CurrSubobject; FractalObjTests++; hit = FALSE; /* * If outside of the bounding box, check to see if we * hit it. */ VecAddScaled(ray->pos, mindist, ray->dir, &curpos); offset = *maxdist; if (OutOfBounds(&curpos, fractalobj->Bounds) && !BoundsIntersect(ray, fractalobj->Bounds, mindist, &offset)) /* * Ray never hit grid space. */ return hit; counter = raynumber++; for (CurrSubobject = fractalobj->SubObjects; CurrSubobject; CurrSubobject = CurrSubobject->Next) if (SubObjectIntersect(fractalobj, CurrSubobject, ray, mindist, maxdist, counter)) hit = TRUE; return hit; } /* * Test for intersection of ray with every of the subobjects * * (based on the code of grid.c) */ int SubObjectIntersect(fractalobj, subobject, ray, mindist, maxdist, counter) FractalObj *fractalobj; SubObject *subobject; Ray *ray; Float mindist; Float *maxdist; unsigned long counter; { int hit; Float offset, tMaxX, tMaxY, tMaxZ; Float tDeltaX, tDeltaY, tDeltaZ; int stepX, stepY, stepZ, outX, outY, outZ, x, y, z; Vector curpos; TriangleClosure *list; hit = FALSE; /* * If outside of the bounding box, check to see if we * hit it. */ VecAddScaled(ray->pos, mindist, ray->dir, &curpos); if (OutOfBounds(&curpos, subobject->Bounds)) { offset = *maxdist; if (!BoundsIntersect(ray, subobject->Bounds, mindist, &offset)) /* * Ray never hit grid space. */ return hit; VecAddScaled(ray->pos, offset, ray->dir, &curpos); } else offset = mindist; /* * tMaxX is the absolute distance from the ray origin we must move * to get to the next voxel in the X * direction. It is incrementally updated * by DDA as we move from voxel-to-voxel. * tDeltaX is the relative amount along the ray we move to * get to the next voxel in the X direction. Thus, * when we decide to move in the X direction, * we increment tMaxX by tDeltaX. */ x = Fx2voxel(subobject, curpos.x); if (x == subobject->XSize) x--; if (fabs(ray->dir.x) < EPSILON) { tMaxX = FAR_AWAY; tDeltaX = 0.; } else if (ray->dir.x < 0.) { tMaxX = offset + (Fvoxel2x(subobject, x) - curpos.x) / ray->dir.x; tDeltaX = subobject->VoxelSize[X] / - ray->dir.x; stepX = outX = -1; } else { tMaxX = offset + (Fvoxel2x(subobject, x+1) - curpos.x) / ray->dir.x; tDeltaX = subobject->VoxelSize[X] / ray->dir.x; stepX = 1; outX = subobject->XSize; } y = Fy2voxel(subobject, curpos.y); if (y == subobject->YSize) y--; if (fabs(ray->dir.y) < EPSILON) { tMaxY = FAR_AWAY; tDeltaY = 0.; } else if (ray->dir.y < 0.) { tMaxY = offset + (Fvoxel2y(subobject, y) - curpos.y) / ray->dir.y; tDeltaY = subobject->VoxelSize[Y] / - ray->dir.y; stepY = outY = -1; } else { tMaxY = offset + (Fvoxel2y(subobject, y+1) - curpos.y) / ray->dir.y; tDeltaY = subobject->VoxelSize[Y] / ray->dir.y; stepY = 1; outY = subobject->YSize; } z = Fz2voxel(subobject, curpos.z); if (z == subobject->ZSize) z--; if (fabs(ray->dir.z) < EPSILON) { tMaxZ = FAR_AWAY; tDeltaZ = 0.; } else if (ray->dir.z < 0.) { tMaxZ = offset + (Fvoxel2z(subobject, z) - curpos.z) / ray->dir.z; tDeltaZ = subobject->VoxelSize[Z] / - ray->dir.z; stepZ = outZ = -1; } else { tMaxZ = offset + (Fvoxel2z(subobject, z+1) - curpos.z) / ray->dir.z; tDeltaZ = subobject->VoxelSize[Z] / ray->dir.z; stepZ = 1; outZ = subobject->ZSize; } /* * Do DDA-algoritm to traverse the grid */ while (TRUE) { list = subobject->Cells[x][y][z]; if (tMaxX < tMaxY && tMaxX < tMaxZ) { if (list) if (ListIntersect(fractalobj, list,ray,counter,offset,maxdist)) hit = TRUE; x += stepX; if (*maxdist < tMaxX || x == outX) break; tMaxX += tDeltaX; } else if (tMaxZ < tMaxY) { if (list) if (ListIntersect(fractalobj, list,ray,counter,offset,maxdist)) hit = TRUE; z += stepZ; if (*maxdist < tMaxZ || z == outZ) break; tMaxZ += tDeltaZ; } else { if (list) if (ListIntersect(fractalobj, list,ray,counter,offset,maxdist)) hit = TRUE; y += stepY; if (*maxdist < tMaxY || y == outY) break; tMaxY += tDeltaY; } } return hit; } /* * Do intersection test with every triangle in the voxel */ static int ListIntersect(fractalobj, list, ray, counter, mindist, maxdist) FractalObj *fractalobj; TriangleClosure *list; Ray *ray; unsigned long counter; Float mindist, *maxdist; { int hit; hit = FALSE; for(;list; list = list->Next) { /* * If object's counter is greater than or equal to the * number associated with the current grid, * don't bother checking again. */ if (list->Triangle->Counter < counter ) { list->Triangle->Counter = counter; if (TriangleIntersect(fractalobj, list->Triangle->Point, ray, mindist, maxdist)) hit = TRUE; } } return hit; } /* * Do intersection test with triangle */ static int TriangleIntersect(fractalobj, point, ray, mindist, maxdist) FractalObj *fractalobj; int point[3]; Ray *ray; Float mindist, *maxdist; { Vector edge1, edge2, edge3; Float k, s; Vector normal, absnorm; Vector pos, dir; Float qi1, qi2, b0, b1, b2; FractalPoint *Points; FractalPoint p0, p1, p2; Points = fractalobj->Closure->Pointpool->Points; p0 = Points[point[0]]; p1 = Points[point[1]]; p2 = Points[point[2]]; VecSub(p1, p0, &edge1); VecSub(p2, p1, &edge2); VecSub(p0, p2, &edge3); /* * Find plane normal. */ VecCross(&edge1, &edge2, &normal); pos = ray->pos; dir = ray->dir; /* * Plane intersection. */ k = dotp(&normal, &dir); if (fabs(k) < EPSILON) return FALSE; s = (dotp(&normal, &p0) - dotp(&normal, &pos)) / k; if (s < mindist || s > *maxdist) return FALSE; /* * Find "dominant" part of normal vector. */ absnorm.x = fabs(normal.x); absnorm.y = fabs(normal.y); absnorm.z = fabs(normal.z); if (absnorm.x > absnorm.y && absnorm.x > absnorm.z) { qi1 = pos.y + s * dir.y; qi2 = pos.z + s * dir.z; b0 = (edge2.y * (qi2 - p1.z) - edge2.z * (qi1 - p1.y)) / normal.x; if (b0 < 0. || b0 > 1.) return FALSE; b1 = (edge3.y * (qi2 - p2.z) - edge3.z * (qi1 - p2.y)) / normal.x; if (b1 < 0. || b1 > 1.) return FALSE; b2 = (edge1.y * (qi2 - p0.z) - edge1.z * (qi1 - p0.y)) / normal.x; if (b2 < 0. || b2 > 1.) return FALSE; } else if (absnorm.y > absnorm.z) { qi1 = pos.x + s * dir.x; qi2 = pos.z + s * dir.z; b0 = (edge2.z * (qi1 - p1.x) - edge2.x * (qi2 - p1.z)) / normal.y; if (b0 < 0. || b0 > 1.) return FALSE; b1 = (edge3.z * (qi1 - p2.x) - edge3.x * (qi2 - p2.z)) / normal.y; if (b1 < 0. || b1 > 1.) return FALSE; b2 = (edge1.z * (qi1 - p0.x) - edge1.x * (qi2 - p0.z)) / normal.y; if (b2 < 0. || b2 > 1.) return FALSE; } else { qi1 = pos.x + s * dir.x; qi2 = pos.y + s * dir.y; b0 = (edge2.x * (qi2 - p1.y) - edge2.y * (qi1 - p1.x)) / normal.z; if (b0 < 0. || b0 > 1.) return FALSE; b1 = (edge3.x * (qi2 - p2.y) - edge3.y * (qi1 - p2.x)) / normal.z; if (b1 < 0. || b1 > 1.) return FALSE; b2 = (edge1.x * (qi2 - p0.y) - edge1.y * (qi1 - p0.x)) / normal.z; if (b2 < 0. || b2 > 1.) return FALSE; } FractalObjHits++; *maxdist = s; fractalobj->Normal = normal; VecAddScaled(ray->pos, s, ray->dir, &fractalobj->IntersectPos); return TRUE; } #define VecEqual(a,b) ((equal((a).x, (b).x)) && (equal((a).y, (b).y)) && (equal((a).z, (b).z))) /* * Calculate Normal on triangle at pos */ int FractalObjNormal(fractalobj, pos, nrm, gnrm) FractalObj *fractalobj; Vector *pos, *nrm, *gnrm; { if (!VecEqual(*pos, fractalobj->IntersectPos)) RLerror(RL_ABORT, "Can't compute fractalobj normal in point other than the last intersection.\n"); *gnrm = fractalobj->Normal; VecNormalize(gnrm); *nrm = *gnrm; return FALSE; } char *FractalObjName() { return fractalobjName; } void FractalObjStats(tests, hits) unsigned long *tests, *hits; { *tests = FractalObjTests; *hits = FractalObjHits; } void FractalObjBounds(fractalobj, bounds) FractalObj *fractalobj; Float bounds[2][3]; { int i, j; SubObject *CurrSubobject; for (i = 0; i <= 2; i++) { bounds[LOW][i] = MAXVALUE; bounds[HIGH][i] = -MAXVALUE; } for(CurrSubobject = fractalobj->SubObjects; CurrSubobject; CurrSubobject = CurrSubobject->Next) for (i = 0; i <= 2; i++) { bounds[LOW][i] = min(bounds[LOW][i], CurrSubobject->Bounds[LOW][i]); bounds[HIGH][i] = max(bounds[HIGH][i], CurrSubobject->Bounds[HIGH][i]); } for (i = 0; i <= 2; i++) for (j = 0; j <= 1; j++) fractalobj->Bounds[j][i] = bounds[j][i]; } void FractalObjMethodRegister(meth) UserMethodType meth; { if (iFractalObjMethods) iFractalObjMethods->user = meth; } Methods *FractalObjMethods() { if (iFractalObjMethods == (Methods *)NULL) { iFractalObjMethods = MethodsCreate(); iFractalObjMethods->create = (GeomCreateFunc *)FractalObjCreate; iFractalObjMethods->methods = FractalObjMethods; iFractalObjMethods->name = FractalObjName; iFractalObjMethods->intersect = FractalObjIntersect; iFractalObjMethods->normal = FractalObjNormal; iFractalObjMethods->bounds = FractalObjBounds; iFractalObjMethods->stats = FractalObjStats; iFractalObjMethods->checkbounds = FALSE; iFractalObjMethods->closed = FALSE; } return iFractalObjMethods; }