in layer1/Basis.cpp [2854:3599]
int BasisMakeMap(CBasis * I, int *vert2prim, CPrimitive * prim, int n_prim,
float *volume,
int group_id, int block_base,
int perspective, float front, float size_hint)
{
float *v;
float ll;
CPrimitive *prm;
int i;
int *tempRef = NULL;
int n = 0, h, q, x, y, z, j, k, l, e;
int extra_vert = 0;
float p[3], dd[3], *d1, *d2, vd[3], cx[3], cy[3];
float *tempVertex = NULL;
float xs, ys;
int remapMode = true; /* remap mode means that some objects will span more
* than one voxel, so we have to worry about populating
* those voxels and also about eliminating duplicates
* when traversing the neighbor lists */
float min[3], max[3], extent[6];
float sep;
float diagonal[3];
float l1, l2;
float bh, ch;
int n_voxel;
int ok = true;
const float _0 = 0.0;
const float _p5 = 0.5;
PRINTFD(I->G, FB_Ray)
" BasisMakeMap: I->NVertex %d [(%8.3f, %8.3f, %8.3f),...]\n", I->NVertex,
I->Vertex[0], I->Vertex[1], I->Vertex[2]
ENDFD;
sep = I->MinVoxel;
if(sep == _0) {
remapMode = false;
sep = I->MaxRadius; /* also will imply no remapping of vertices */
}
/* we need to get a sense of the actual size in order to avoid sep being too small */
v = I->Vertex;
min[0] = max[0] = v[0];
min[1] = max[1] = v[1];
min[2] = max[2] = v[2];
v += 3;
{
int a;
for(a = 1; a < I->NVertex; a++) {
if(min[0] > v[0])
min[0] = v[0];
if(max[0] < v[0])
max[0] = v[0];
if(min[1] > v[1])
min[1] = v[1];
if(max[1] < v[1])
max[1] = v[1];
if(min[2] > v[2])
min[2] = v[2];
if(max[2] < v[2])
max[2] = v[2];
v += 3;
}
}
if(volume) {
/*
if(min[0] > volume[0]) min[0] = volume[0];
if(max[0] < volume[1]) max[0] = volume[1];
if(min[1] > volume[2]) min[1] = volume[2];
if(max[1] < volume[3]) max[1] = volume[3];
if(min[2] > (-volume[5])) min[2] = (-volume[5]);
if(max[2] < (-volume[4])) max[2] = (-volume[4]);
*/
if(min[0] < volume[0])
min[0] = volume[0];
if(max[0] > volume[1])
max[0] = volume[1];
if(min[1] < volume[2])
min[1] = volume[2];
if(max[1] > volume[3])
max[1] = volume[3];
if(min[2] > (-volume[5]))
min[2] = (-volume[5]);
if(max[2] < (-volume[4]))
max[2] = (-volume[4]);
PRINTFB(I->G, FB_Ray, FB_Debugging)
" BasisMakeMap: (%8.3f,%8.3f),(%8.3f,%8.3f),(%8.3f,%8.3f)\n",
volume[0], volume[1], volume[2], volume[3], volume[4], volume[5]
ENDFB(I->G);
}
/* don't break up space unnecessarily if we only have a few vertices... */
if(I->NVertex) {
l1 = (float) fabs(max[0] - min[0]);
l2 = (float) fabs(max[1] - min[1]);
if(l1 < l2)
l1 = l2;
l2 = (float) fabs(max[2] - min[2]);
if(l1 < l2)
l1 = l2;
if(l1 < kR_SMALL4)
l1 = 100.0;
if(I->NVertex < (l1 / sep))
sep = (l1 / I->NVertex);
}
if(Feedback(I->G, FB_Ray, FB_Debugging)) {
dump3f(min, " BasisMakeMap: min");
dump3f(max, " BasisMakeMap: max");
dump3f(I->Vertex, " BasisMakeMap: I->Vertex");
fflush(stdout);
}
sep = MapGetSeparation(I->G, sep, max, min, diagonal); /* this needs to be a minimum
* estimate of the actual value */
/* here we have to carry out a complicated work-around in order to
* efficiently encode our primitives into the map in a way that doesn't
* require expanding the map cutoff to the size of the largest object*/
if(remapMode) {
int a, b, c;
if(sep < size_hint) /* this keeps us from wasting time & memory on unnecessary subdivision */
sep = size_hint;
{
int *vert2prim_a = vert2prim;
for(a = 0; a < I->NVertex; a++) {
prm = prim + *(vert2prim_a++);
switch (prm->type) {
case cPrimTriangle:
case cPrimCharacter:
if(a == prm->vert) { /* only do this calculation for one of the three vertices */
l1 = (float) length3f(I->Precomp + I->Vert2Normal[a] * 3);
l2 = (float) length3f(I->Precomp + I->Vert2Normal[a] * 3 + 3);
if((l1 >= sep) || (l2 >= sep)) {
b = (int) ceil(l1 / sep) + 1;
c = (int) ceil(l2 / sep) + 1;
extra_vert += 4 * b * c;
}
}
break;
case cPrimCone:
case cPrimCylinder:
case cPrimSausage:
if((prm->l1 + 2 * prm->r1) >= sep) {
q = ((int) (2 * (floor(prm->r1 / sep) + 1))) + 1;
q = q * q * ((int) ceil((prm->l1 + 2 * prm->r1) / sep) + 2);
extra_vert += q;
}
break;
case cPrimEllipsoid:
case cPrimSphere:
if(prm->r1 >= sep) {
b = (int) (2 * floor(prm->r1 / sep) + 1);
extra_vert += (b * b * b);
}
break;
}
} /* for */
}
extra_vert += I->NVertex;
if (ok)
tempVertex =
CacheAlloc(I->G, float, extra_vert * 3, group_id, cCache_basis_tempVertex);
CHECKOK(ok, tempVertex);
if (ok)
tempRef = CacheAlloc(I->G, int, extra_vert, group_id, cCache_basis_tempRef);
CHECKOK(ok, tempRef);
ErrChkPtr(I->G, tempVertex); /* can happen if extra vert is unreasonable */
ErrChkPtr(I->G, tempRef);
/* lower indexes->flags, top is ref->lower index */
ok &= !I->G->Interrupt;
if (ok){
float *vv, *d;
int *vert2prim_a = vert2prim;
n = I->NVertex;
v = tempVertex;
vv = I->Vertex;
memcpy(v, vv, n * sizeof(float) * 3);
vv += n * 3;
v += n * 3;
for(a = 0; ok && a < I->NVertex; a++) {
prm = prim + *(vert2prim_a++);
switch (prm->type) {
case cPrimTriangle:
case cPrimCharacter:
if(a == prm->vert) {
/* only do this calculation for one of the three vertices */
d1 = I->Precomp + I->Vert2Normal[a] * 3;
d2 = I->Precomp + I->Vert2Normal[a] * 3 + 3;
vv = I->Vertex + a * 3;
l1 = (float) length3f(d1);
l2 = (float) length3f(d2);
if((l1 >= sep) || (l2 >= sep)) {
b = (int) floor(l1 / sep) + 1;
c = (int) floor(l2 / sep) + 1;
extra_vert += b * c;
bh = (float) (b / 2) + 1;
ch = (float) (c / 2) + 1;
for(x = 0; x < bh; x++) {
const float xb = (float) x / b;
for(y = 0; y < ch; y++) {
const float yc = (float) y / c;
*(v++) = vv[0] + (d1[0] * xb) + (d2[0] * yc);
*(v++) = vv[1] + (d1[1] * xb) + (d2[1] * yc);
*(v++) = vv[2] + (d1[2] * xb) + (d2[2] * yc);
tempRef[n++] = a;
}
}
for(x = 0; x < bh; x++) {
const float xb = (float) x / b;
for(y = 0; y < ch; y++) {
const float yc = (float) y / c;
if((xb + yc) < _p5) {
*(v++) = vv[0] + d1[0] * (_p5 + xb) + (d2[0] * yc);
*(v++) = vv[1] + d1[1] * (_p5 + xb) + (d2[1] * yc);
*(v++) = vv[2] + d1[2] * (_p5 + xb) + (d2[2] * yc);
tempRef[n++] = a;
*(v++) = vv[0] + (d1[0] * xb) + d2[0] * (_p5 + yc);
*(v++) = vv[1] + (d1[1] * xb) + d2[1] * (_p5 + yc);
*(v++) = vv[2] + (d1[2] * xb) + d2[2] * (_p5 + yc);
tempRef[n++] = a;
}
}
}
}
} /* if */
break;
case cPrimCone:
case cPrimCylinder:
case cPrimSausage:
ll = prm->l1 + 2 * prm->r1;
if(ll >= sep) {
d = I->Normal + 3 * I->Vert2Normal[a];
vv = I->Vertex + a * 3;
get_system1f3f(d, cx, cy); /* creates an orthogonal system about d */
p[0] = vv[0] - d[0] * prm->r1;
p[1] = vv[1] - d[1] * prm->r1;
p[2] = vv[2] - d[2] * prm->r1;
dd[0] = d[0] * sep;
dd[1] = d[1] * sep;
dd[2] = d[2] * sep;
q = (int) floor(prm->r1 / sep) + 1;
while(1) {
vd[0] = (p[0] += dd[0]);
vd[1] = (p[1] += dd[1]);
vd[2] = (p[2] += dd[2]);
for(x = -q; x <= q; x++) {
for(y = -q; y <= q; y++) {
xs = x * sep;
ys = y * sep;
*(v++) = vd[0] + xs * cx[0] + ys * cy[0];
*(v++) = vd[1] + xs * cx[1] + ys * cy[1];
*(v++) = vd[2] + xs * cx[2] + ys * cy[2];
tempRef[n++] = a;
}
}
if(ll <= _0)
break;
ll -= sep;
}
}
break;
case cPrimEllipsoid:
case cPrimSphere:
if(prm->r1 >= sep) {
q = (int) floor(prm->r1 / sep);
vv = I->Vertex + a * 3;
for(x = -q; x <= q; x++) {
for(y = -q; y <= q; y++) {
for(z = -q; z <= q; z++) {
*(v++) = vv[0] + x * sep;
*(v++) = vv[1] + y * sep;
*(v++) = vv[2] + z * sep;
tempRef[n++] = a;
}
}
}
}
break;
} /* end of switch */
ok &= !I->G->Interrupt;
}
}
if(n > extra_vert) {
printf("BasisMakeMap: %d>%d\n", n, extra_vert);
ErrFatal(I->G, "BasisMakeMap",
"used too many extra vertices (this indicates a bug)...\n");
}
PRINTFB(I->G, FB_Ray, FB_Blather)
" BasisMakeMap: %d total vertices\n", n ENDFB(I->G);
if (ok){
if(volume) {
v = tempVertex;
min[0] = max[0] = v[0];
min[1] = max[1] = v[1];
min[2] = max[2] = v[2];
v += 3;
if(Feedback(I->G, FB_Ray, FB_Debugging)) {
dump3f(min, " BasisMakeMap: remapped min");
dump3f(max, " BasisMakeMap: remapped max");
fflush(stdout);
}
for(a = 1; a < n; a++) {
if(min[0] > v[0])
min[0] = v[0];
if(max[0] < v[0])
max[0] = v[0];
if(min[1] > v[1])
min[1] = v[1];
if(max[1] < v[1])
max[1] = v[1];
if(min[2] > v[2])
min[2] = v[2];
if(max[2] < v[2])
max[2] = v[2];
v += 3;
}
if(Feedback(I->G, FB_Ray, FB_Debugging)) {
dump3f(min, " BasisMakeMap: remapped min");
dump3f(max, " BasisMakeMap: remapped max");
fflush(stdout);
}
if(min[0] < volume[0])
min[0] = volume[0];
if(max[0] > volume[1])
max[0] = volume[1];
if(min[1] < volume[2])
min[1] = volume[2];
if(max[1] > volume[3])
max[1] = volume[3];
if(min[2] < (-volume[5]))
min[2] = (-volume[5]);
if(max[2] > (-volume[4]))
max[2] = (-volume[4]);
extent[0] = min[0];
extent[1] = max[0];
extent[2] = min[1];
extent[3] = max[1];
extent[4] = min[2];
extent[5] = max[2];
PRINTFB(I->G, FB_Ray, FB_Blather)
" BasisMakeMap: Extent [%8.2f %8.2f] [%8.2f %8.2f] [%8.2f %8.2f]\n",
extent[0], extent[1], extent[2], extent[3], extent[4], extent[5]
ENDFB(I->G);
I->Map = MapNewCached(I->G, -sep, tempVertex, n, extent, group_id, block_base);
CHECKOK(ok, I->Map);
} else {
I->Map = MapNewCached(I->G, sep, tempVertex, n, NULL, group_id, block_base);
CHECKOK(ok, I->Map);
}
}
if (ok)
n_voxel = I->Map->Dim[0] * I->Map->Dim[1] * I->Map->Dim[2];
if (!ok){
} else if(perspective) {
/* this is a new optimization which prevents primitives
contained entirely within a single voxel from spilling over
into neighboring voxels for performance of intersection
checks (NOTE: this currently only helps with triangles, and
is only useful in optimizing between Z layers). The main
impact of this optimization is to reduce the size of the
express list (I->Map->Elist) array by about 3-fold */
int *prm_spanner = NULL;
int *spanner = NULL;
float *v = tempVertex;
int j, k, l, jj, kk, ll;
int nVertex = I->NVertex;
int prm_index;
MapType *map = I->Map;
prm_spanner = Calloc(int, n_prim);
CHECKOK(ok, prm_spanner);
if (ok)
spanner = Calloc(int, n);
CHECKOK(ok, spanner);
/* figure out which primitives span more than one voxel */
for(a = 0; ok && a < n; a++) {
if(a < nVertex)
prm_index = vert2prim[a];
else
prm_index = vert2prim[tempRef[a]];
if(!prm_spanner[prm_index]) {
prm = prim + prm_index;
{
float *vv = prm->vert * 3 + I->Vertex;
switch (prm->type) {
case cPrimTriangle:
case cPrimCharacter:
MapLocus(map, v, &j, &k, &l);
MapLocus(map, vv, &jj, &kk, &ll);
if((j != jj) || (k != kk) || (l != ll)) {
prm_spanner[prm_index] = 1;
}
break;
default: /* currently we aren't optimizing other primitives */
prm_spanner[prm_index] = 1;
break;
}
}
}
v += 3;
ok &= !I->G->Interrupt;
}
/* now flag associated vertices */
for(a = 0; ok && a < n; a++) {
if(a < nVertex)
prm_index = vert2prim[a];
else
prm_index = vert2prim[tempRef[a]];
spanner[a] = prm_spanner[prm_index];
ok &= !I->G->Interrupt;
}
/* and do the optimized expansion */
ok &= MapSetupExpressPerp(I->Map, tempVertex, front, n, true, spanner);
FreeP(spanner);
FreeP(prm_spanner);
} else if(n_voxel < (3 * n)) {
ok &= MapSetupExpressXY(I->Map, n, true);
} else {
ok &= MapSetupExpressXYVert(I->Map, tempVertex, n, true);
}
if (ok) {
MapType *map = I->Map;
int *sp, *ip, *ip0, ii;
int *elist = map->EList, *ehead = map->EHead;
int *elist_new = elist, *ehead_new = ehead;
int newelem = 0, neelem = -map->NEElem;
int i_nVertex = I->NVertex;
const int iMin0 = map->iMin[0];
const int iMin1 = map->iMin[1];
const int iMin2 = map->iMin[2];
const int iMax0 = map->iMax[0];
const int iMax1 = map->iMax[1];
const int iMax2 = map->iMax[2];
/* now do a filter-reassignment pass to remap fake vertices
to the original line vertex while deleting duplicate entries */
memset(tempRef, 0, sizeof(int) * i_nVertex);
if(n_voxel < (3 * n)) { /* faster to traverse the entire map */
int *start;
for(a = iMin0; a <= iMax0; a++) {
for(b = iMin1; b <= iMax1; b++) {
for(c = iMin2; c <= iMax2; c++) {
start = MapEStart(map, a, b, c);
h = *start;
if((h < 0) && (h >= neelem)) {
sp = elist - h;
*(start) = -h; /* flip sign */
i = *(sp++);
if(ehead_new != ehead) {
ehead_new[(start - ehead) - 1] = newelem;
ip = ip0 = (elist_new + newelem);
} else {
ip = ip0 = (sp - 1);
}
while(i >= 0) {
int ii = *(sp++);
if(i >= i_nVertex) /* reference -- remap */
i = tempRef[i];
if(!tempRef[i]) { /*eliminate duplicates */
tempRef[i] = 1;
*(ip++) = i;
}
i = ii;
}
*(ip++) = -1; /* terminate list */
newelem += (ip - ip0);
/* now reset flags efficiently */
i = *(ip0++);
while(i >= 0) { /* unroll */
ii = *(ip0++);
tempRef[i] = 0;
if((i = ii) >= 0) {
ii = *(ip0++);
tempRef[i] = 0;
if((i = ii) >= 0) {
ii = *(ip0++);
tempRef[i] = 0;
i = ii;
}
}
}
}
} /* for c */
} /* for b */
} /* for a */
} else {
int *start;
int jm1, jp1, km1, kp1, lm1, lp1;
MapType *map = I->Map;
v = tempVertex;
for(e = 0; e < n; e++) {
MapLocus(map, v, &j, &k, &l);
jm1 = j - 1;
jp1 = j + 1;
km1 = k - 1;
kp1 = k + 1;
lm1 = l - 1;
lp1 = l + 1;
if(perspective) { /* for normal maps */
int *iPtr1 =
map->EHead + ((j - 1) * map->D1D2) + ((k - 1) * map->Dim[2]) + (l - 1);
for(a = jm1; a <= jp1; a++) {
int *iPtr2 = iPtr1;
if((a >= iMin0) && (a <= iMax0)) {
for(b = km1; b <= kp1; b++) {
int *iPtr3 = iPtr2;
if((b >= iMin1) && (b <= iMax1)) {
for(c = lm1; c <= lp1; c++) {
if((c >= iMin2) && (c <= iMax2)) {
start = iPtr3;
/*start = MapEStart(map,a,b,c); */
h = *start;
if((h < 0) && (h >= neelem)) {
sp = elist - h;
(*start) = -h; /* no repeat visits */
i = *(sp++);
if(ehead_new != ehead) {
ehead_new[(start - ehead)] = newelem;
ip = ip0 = (elist_new + newelem);
} else {
ip = ip0 = (sp - 1);
}
while(i >= 0) {
int ii = *(sp++);
if(i >= i_nVertex)
i = tempRef[i];
if(!tempRef[i]) { /*eliminate duplicates */
tempRef[i] = 1;
*(ip++) = i;
}
i = ii;
}
*(ip++) = -1; /* terminate list */
newelem += (ip - ip0);
/* now reset flags efficiently */
i = *(ip0++);
while(i >= 0) { /* unroll */
ii = *(ip0++);
tempRef[i] = 0;
if((i = ii) >= 0) {
ii = *(ip0++);
tempRef[i] = 0;
if((i = ii) >= 0) {
ii = *(ip0++);
tempRef[i] = 0;
i = ii;
}
}
}
} /* h > 0 */
}
iPtr3++;
}
}
iPtr2 += map->Dim[2];
}
}
iPtr1 += map->D1D2;
}
} else { /* for XY maps... */
c = l;
{
int *iPtr1 =
ehead + ((j - 1) * map->D1D2) + ((k - 1) * map->Dim[2]) + c;
if((c >= iMin2) && (c <= iMax2)) {
for(a = jm1; a <= jp1; a++) {
int *iPtr2 = iPtr1;
if((a >= iMin0) && (a <= iMax0)) {
for(b = km1; b <= kp1; b++) {
if((b >= iMin1) && (b <= iMax1)) {
start = iPtr2;
/*start = MapEStart(map,a,b,c); */
h = *start;
if((h < 0) && (h >= neelem)) {
sp = elist - h;
(*start) = -h; /* no repeat visits */
i = *(sp++);
if(ehead_new != ehead) {
ehead_new[(start - ehead)] = newelem;
ip = ip0 = (elist_new + newelem);
} else {
ip = ip0 = sp - 1;
}
while(i >= 0) {
int ii = *(sp++);
if(i >= i_nVertex)
i = tempRef[i];
if(!tempRef[i]) { /*eliminate duplicates */
tempRef[i] = 1;
*(ip++) = i;
}
i = ii;
}
*(ip++) = -1; /* terminate list */
newelem += (ip - ip0);
/* now reset flags efficiently */
i = *(ip0++);
while(i >= 0) { /* unroll */
ii = *(ip0++);
tempRef[i] = 0;
if((i = ii) >= 0) {
ii = *(ip0++);
tempRef[i] = 0;
if((i = ii) >= 0) {
ii = *(ip0++);
tempRef[i] = 0;
i = ii;
}
}
}
} /* h > 0 */
}
iPtr2 += map->Dim[2];
} /* for b */
}
iPtr1 += map->D1D2;
} /* for a */
}
}
}
v += 3; /* happens for EVERY e! */
} /* for e */
}
if(ehead != ehead_new) {
CacheFreeP(I->G, map->EHead, group_id, block_base + cCache_map_ehead_offset,
false);
VLACacheFreeP(I->G, map->EList, group_id, block_base + cCache_map_elist_offset,
false);
map->EList = elist_new;
map->EHead = ehead_new;
map->NEElem = newelem;
MemoryCacheReplaceBlock(I->G, group_id, block_base + cCache_map_ehead_new_offset,
block_base + cCache_map_ehead_offset);
MemoryCacheReplaceBlock(I->G, group_id, block_base + cCache_map_elist_new_offset,
block_base + cCache_map_elist_offset);
VLACacheSize(G, map->EList, int, map->NEElem, group_id,
block_base + cCache_map_elist_offset);
}
}
CacheFreeP(I->G, tempVertex, group_id, cCache_basis_tempVertex, false);
CacheFreeP(I->G, tempRef, group_id, cCache_basis_tempRef, false);
} else {
/* simple sphere mode */
I->Map = MapNewCached(I->G, -sep, I->Vertex, I->NVertex, NULL, group_id, block_base);
CHECKOK(ok, I->Map);
if (ok){
if(perspective) {
ok &= MapSetupExpressPerp(I->Map, I->Vertex, front, I->NVertex, false, NULL);
} else {
ok &= MapSetupExpressXYVert(I->Map, I->Vertex, I->NVertex, false);
}
}
}
return ok;
}