#include "testDB.h"
#include "polyPairTest.h"
#include "matrix.h"
#include "collIntMat.h"

static ConNode *objObjIntersection(ObjNode *, ObjNode *, ConNode *, Status *);
static void findMinMax(Transform *, Tupple, Tupple, Tupple *, Tupple *);
static Boolean findSubMatrix(Tupple, Tupple, Tupple, Tupple, Tupple, int *, 
					int *);
static int listPolys(SEADS *, Tupple, Tupple, int **);
static PairNode *addPairNodeList(int *, int *, int, int, PairNode *);
static Boolean globalIntersect(Tupple, Tupple, Tupple, Tupple);

#ifdef STATS
/* statistical variables */
extern Ulong boundBoxTestCount;
extern Ulong boundBoxTestHit;
extern Ulong voxelTestCount;
extern Ulong polyTestCount;
#endif

/*
 * Given the database of objects, tests for intersections etc.  Returns the
 * status of the database (objects intersected, objects just in contact,
 * nothing) for the given time period, and if appropriate the list of contact
 * points/areas, or points of intersection.  Each object has a transformation 
 * matrix that enables the points in it's local reference frame to be 
 * transformed to the global frame at this point in time.
 *
 * Params	objdb		object database
 *			s			return the status of the database
 *
 * Globals	none
 *
 * Returns	the list of contact or intersection information (if any), and 
 *			the status of the database in "s"
 */
ConNode *test4Intersection(objdb, s)
ObjDB objdb;
Status *s;
{
	ObjNode *obja, *objb;
	ConNode *cn;
	Status stemp = Nothing;
	Ulong ui, uj, ulimit;

	cn = NULL;
	*s = Nothing;

	if (objdb.cim == NULL) {
		/* collision interest matrix isn't being used so test each object
		 * against all others
		 */
		
		for (obja=objdb.objects; obja != NULL; obja=obja->next) {
			if (!(obja->valid)) continue;
			for (objb=obja->next; objb != NULL; objb=objb->next) {
				if (!(objb->valid)) continue;
				
				if (globalIntersect(obja->gmin, obja->gmax, objb->gmin, 
					objb->gmax)) {

					#ifdef STATS
					boundBoxTestCount++;
					#endif
					
					/* global bounding boxes intersect, so test these objects */
					cn = objObjIntersection(obja, objb, cn, &stemp);

					#ifdef STATS
					if (stemp == Intersection || stemp == Contact)
						boundBoxTestHit++;
					#endif

					if (stemp == Intersection) {
						*s = Intersection;
					} else if (stemp == Contact && *s != Intersection) {
						*s = Contact;
					}
				}
			}
		}
	} else {

		/* using the collision interest matrix */
		ulimit = objdb.objects->id;
		for (ui=ulimit,obja=objdb.objects; ui > 0; ui--,obja=obja->next) {
			if (!(obja->valid)) continue;
		for (uj=ui-1,objb=obja->next; objb != NULL; uj--,objb=objb->next) {
			if (!(objb->valid)) continue;
			if (queryCollIntMatrix(ui, uj, objdb.cim)) {
			/* objects have to be tested*/
				if (globalIntersect(obja->gmin, obja->gmax, objb->gmin, 
					objb->gmax)) {

					#ifdef STATS
					boundBoxTestCount++;
					#endif

					/* global bounding boxes intersect, so test these objects */
					cn = objObjIntersection(obja, objb, cn, &stemp);

					#ifdef STATS
					if (stemp == Intersection || stemp == Contact)
						boundBoxTestHit++;
					#endif

					if (stemp == Intersection) {
						*s = Intersection;
					} else if (stemp == Contact && *s != Intersection) {
						*s = Contact;
					}
				}
					
			}
		}
		}
	}

	return cn;
}

/* 
 * given an object, calculates its bounding box in the global reference frame
 * by applying the transformation matrix to each of the corners of its local
 * reference box, and finding the maximum and minimum values.  Stores the
 * result in obj->gmin and obj->gmax.
 *
 * Params	obj		object
 *
 * Globals	none
 *
 * Returns	nothing
 */
void calcBoundingBox(obj)
ObjNode *obj;
{
	Tupple cnr[8], temp;
	int i;

	cnr[0] = obj->min;
	cnr[1].i = obj->min.i; cnr[1].j = obj->min.j; cnr[1].k = obj->max.k;
	cnr[2].i = obj->min.i; cnr[2].j = obj->max.j; cnr[2].k = obj->min.k;
	cnr[3].i = obj->max.i; cnr[3].j = obj->min.j; cnr[3].k = obj->min.k;
	cnr[4].i = obj->max.i; cnr[4].j = obj->max.j; cnr[4].k = obj->min.k;
	cnr[5].i = obj->max.i; cnr[5].j = obj->min.j; cnr[5].k = obj->max.k;
	cnr[6].i = obj->min.i; cnr[6].j = obj->max.j; cnr[6].k = obj->max.k;
	cnr[7] = obj->max;

	applyTransform(obj->loc2glob, cnr[0], &temp);
	obj->gmin = temp;
	obj->gmax = temp;

	for (i=1; i<8 ;i++) {
		applyTransform(obj->loc2glob, cnr[i], &temp);
		if (obj->gmin.i > temp.i) obj->gmin.i = temp.i;
		if (obj->gmin.j > temp.j) obj->gmin.j = temp.j;
		if (obj->gmin.k > temp.k) obj->gmin.k = temp.k;
		if (obj->gmax.i < temp.i) obj->gmax.i = temp.i;
		if (obj->gmax.j < temp.j) obj->gmax.j = temp.j;
		if (obj->gmax.k < temp.k) obj->gmax.k = temp.k;
	}
}

/*
 * Given two objects, test for interpenetration or contact.  First determines
 * which subset of voxels of SEADS of each object is involved.  If one of the
 * subsets (of object A or object B) is empty then return "nothing", else
 * transform each voxel in B's subset into A's reference frame and find out 
 * which voxels of A it interferes with and then form polygon pairs to test.
 * Remove any duplicates, then test the polygons.
 *
 * The objects must have loc2glob transforms -- if they don't have the inverse
 * transforms, they are calculated.
 *
 * Params	obja		object A
 *			objb		object B
 *			cn			ConNode pointer (where to store the contact info)
 *			s			return the status
 *
 * Globals	none
 *
 * Returns	the contact list, and the status of the objects using "s"i
 */
static ConNode *objObjIntersection(obja, objb, cn, s)
ObjNode *obja, *objb;
ConNode *cn;
Status *s;
{
	Transform tab, tba;
	int i, j, k, dims[3], Bindices[6], Aindices[6], temp1, temp2;
	int *alist, total, index;
	Tupple min, max, tmin, tmax;
	PairNode *plist;
	Boolean occupied;

	#ifdef STATS
	int voxCount = 0;
	int tCount;
	#endif

	plist = NULL;
	alist = NULL;

	/* Calculate the inverse transforms of the objects if needed.
	 * loc2glob transforms are needed for all objects, but inverse transforms
	 * aren't, so calculate them when needed. That is, in this function
	 * because the objects are being tested.
	 */
	if (obja->glob2loc == NULL) {
		obja->glob2loc = (Transform *)malloc(sizeof(Transform));
		if (obja->glob2loc == NULL) {
			fprintf(stderr, "Error: couldn't malloc Transform\n");
			exit(1);
		}
		if (!matrixInvert(obja->loc2glob, obja->glob2loc)) {
			fprintf(stderr, "Error: couldn't invert transformation matrix of ");
			fprintf(stderr, "object %u\n", obja->id);
			exit(2);
		}
	}

	if (objb->glob2loc == NULL) {
		objb->glob2loc = (Transform *)malloc(sizeof(Transform));
		if (objb->glob2loc == NULL) {
			fprintf(stderr, "Error: couldn't malloc Transform\n");
			exit(1);
		}
		if (!matrixInvert(objb->loc2glob, objb->glob2loc)) {
			fprintf(stderr, "Error: couldn't invert transformation matrix of ");
			fprintf(stderr, "object %u\n", objb->id);
			exit(2);
		}
	}

	/* calculate Trans BA = Trans AG^-1 x Trans BG */
	matrixMult(obja->glob2loc, objb->loc2glob, &tba);
	/* find the minimum and maximum values of B in A's frame */
	findMinMax(&tba, objb->min, objb->max, &min, &max);

	/* have minimum and maximum values of obj B in A's frame, so work out
	 * the interference submatrix of A, ie. the range of indices involved
	 */
	dims[0] = obja->seads->dimi;
	dims[1] = obja->seads->dimj;
	dims[2] = obja->seads->dimk;
	if (!findSubMatrix(min, max, obja->min, obja->max, 
		obja->seads->voxSizes, dims, Aindices)) {

		/* no voxels overlapping so no intersection can take place */
		*s = Nothing;
		return cn;
	}

	#ifdef STATS
	voxCount = 1;
	for (i=5; i < 1; i=i-2) {
		tCount = Aindices[i] - Aindices[i-1];
		if (tCount <= 0) tCount = 1;
		voxCount *= tCount;
	}
	voxelTestCount += voxCount;
	#endif

	/* check to see if the interference submatrix contains any polygons */
	occupied = False;
	for (k=Aindices[4]; k <= Aindices[5]; k++) {
		temp2 = k * dims[0] * dims[1];
		for (j=Aindices[2]; j <= Aindices[3]; j++) {
			temp1 = j * dims[0];
			for (i=Aindices[0]; i <= Aindices[1]; i++) {
				if (obja->seads->matrix[i+temp1+temp2].total > 0) {
					occupied = True;
					break;
				}
			}
			if (occupied) break;
		}
		if (occupied) break;
	}

	if (!occupied) {
		/* no polygons in the interference sub matrix so no collision
		 * possible
		 */

		*s = Nothing;
		return cn;
	}

	/* calculate Trans AB = Trans BG^-1 x Trans AG */
	/* matrixMult(obja->loc2glob, objb->glob2loc, &tab); */
	matrixMult(objb->glob2loc, obja->loc2glob, &tab);
	/* find min and max of A in B's reference frame */
	findMinMax(&tab, obja->min, obja->max, &min, &max);

	/* have minimum and maximum values of obj A in B's frame, so work out
	 * the interference submatrix of B, ie. the range of indices involved
	 */
	dims[0] = objb->seads->dimi;
	dims[1] = objb->seads->dimj;
	dims[2] = objb->seads->dimk;

	if (!findSubMatrix(min, max, objb->min, objb->max, 
		objb->seads->voxSizes, dims, Bindices)) {

		/* no voxels overlapping so no intersection can take place */
		*s = Nothing;
		return cn;
	}


	#ifdef STATS
	voxCount = 1;
	for (i=5; i < 1; i=i-2) {
		tCount = Bindices[i] - Bindices[i-1];
		if (tCount <= 0) tCount = 1;
		voxCount *= tCount;
	}
	voxelTestCount += voxCount;
	#endif

	/* Now we have B's interference sub matrix, so go through each voxel.
	 * If it contains polygons, transform it and find which voxels in A's
	 * frame it interferes with, and then form pairs of polygons that need
	 * to be tested.
	 */
	for (k=Bindices[4]; k <= Bindices[5]; k++) {
		temp2 = k * dims[0] * dims[1];
		for (j=Bindices[2]; j <= Bindices[3]; j++) {
			temp1 = j * dims[0];
			for (i=Bindices[0]; i <= Bindices[1]; i++) {
				if (objb->seads->matrix[(index=i+temp1+temp2)].total > 0) {

					/* calculate corners of voxel */
					tmin.i = i * objb->seads->voxSizes.i;
					tmin.j = j * objb->seads->voxSizes.j;
					tmin.k = k * objb->seads->voxSizes.k;
					tmax.i = tmin.i + objb->seads->voxSizes.i;
					tmax.j = tmin.j + objb->seads->voxSizes.j;
					tmax.k = tmin.k + objb->seads->voxSizes.k;

					findMinMax(&tba, tmin, tmax, &min, &max);

					total = listPolys(obja->seads, min, max, &alist);

					/* have two arrays of polygon indices, create pairs
					 * and insert into list
					 */
					plist = addPairNodeList(alist, 
								objb->seads->matrix[index].a, total, 
								objb->seads->matrix[index].total, plist);
					if (alist != NULL)
						free(alist);
					alist = NULL;
				}
			}
		}
	}

	*s = Nothing;
	/* now we have a list of polygon pairs that need to be checked for
	 * intersection, so check them.
	 */
	return polyPairTest(plist, obja, objb, &tab, &tba, cn, s);
}

/*
 * Finds the minimum and maximum values of a box after it undergoes
 * transformation from obj B's reference frame into obj A's reference frame.
 *
 * Params	trans		transformation (B to A)
 *			Bmin		box's min in B frame
 *			Bmax		box's max in B frame
 *			Amin		box's min in A frame (not simply Bmin transformed)
 *			Amax		box's max in A frame (not simply Bmax transformed)
 *
 * Globals	none
 *
 * Returns	the min and max values of obj B in A's ref frame
 */
static void findMinMax(trans, Bmin, Bmax, Amin, Amax)
Transform *trans;
Tupple Bmin, Bmax, *Amin, *Amax;
{
	int i, j, k;
	Tupple temp1, temp2;
	Real error;

	/* transform each corner of B's box into A's frame and find min and
	 * max values.
	 */
	applyTransform(trans, Bmin, Amin); /* setup initial values */
	*Amax = *Amin;

	for (i=0; i < 2; i++) {
		if (i == 0) temp1.i = Bmin.i;
		else temp1.i = Bmax.i;
		for (j=0; j < 2; j++) {
			if (j == 0) temp1.j = Bmin.j;
			else temp1.j = Bmax.j;
			for (k=0; k < 2; k++) {
				if (k == 0) temp1.k = Bmin.k;
				else temp1.k = Bmax.k;
				applyTransform(trans, temp1, &temp2);

				/* decide if the transformed values are minimums or maximums
				 * in A's ref frame
				 */
				if (temp2.i < Amin->i) Amin->i = temp2.i;
				if (temp2.i > Amax->i) Amax->i = temp2.i;
				if (temp2.j < Amin->j) Amin->j = temp2.j;
				if (temp2.j > Amax->j) Amax->j = temp2.j;
				if (temp2.k < Amin->k) Amin->k = temp2.k;
				if (temp2.k > Amax->k) Amax->k = temp2.k;
			}
		}
	}

}

/*
 * finds the interference submatrix (the range of indices) given the minimums
 * and maximums of the 2 objects in the same reference frame.  Because of
 * floating point arithmetic errors, the minimum and maximum values are
 * "fuzzed" so that we don't miss out on voxels that should be part of the
 * interference submatrix.
 *
 * Params	Amin		minimums of obj A
 *			Amax		maximums of obj A
 *			Bmin		minimums of obj B
 *			Bmax		maximums of obj B
 *			voxSizes	sizes of A's voxels
 *			indices		used to return the range of indices
 *
 * Globals	none
 *
 * Returns	True if the given values have an intersection set (ie form a box)
 *				(also gives the range of indices in "indices")
 *			False if they don't
 */
static Boolean findSubMatrix(Amin, Amax, Bmin, Bmax, voxSizes, dims, indices)
Tupple Amin, Amax, Bmin, Bmax, voxSizes;
int *dims, *indices;
{
	int i;
	Real min1, min2, max1, max2, *fi, *fa;
	Real iFuz, jFuz, kFuz;
	Tupple setMin, setMax;

	iFuz = voxSizes.i * ETA;
	jFuz = voxSizes.j * ETA;
	kFuz = voxSizes.k * ETA;

	for (i=0; i < 3; i++) {
		switch (i) {
			case 0: min1 = Amin.i - iFuz;
					max1 = Amax.i + iFuz;
					min2 = Bmin.i - iFuz;
					max2 = Bmax.i + iFuz;
					fi = &setMin.i; fa = &setMax.i;
					break;
			case 1: min1 = Amin.j - jFuz;
					max1 = Amax.j + jFuz;
					min2 = Bmin.j - jFuz;
					max2 = Bmax.j + jFuz;
					fi = &setMin.j; fa = &setMax.j;
					break;
			case 2: min1 = Amin.k - kFuz;
					max1 = Amax.k + kFuz;
					min2 = Bmin.k - kFuz;
					max2 = Bmax.k + kFuz;
					fi = &setMin.k; fa = &setMax.k;
					break;
		}

		if (min1 <= min2 && min2 <= max1) {
			*fi = min2;
			*fa = MIN(max1, max2);
		} else if (min2 <= min1 && min1 <= max2) {
			*fi = min1;
			*fa = MIN(max1, max2);
		} else {
			/* there is no overlap between the two ranges along this
			 * axis, so there can be no voxels in common */
			return False;
		}
	}

	/* have found the intersection sets, translate these into voxel indices */
	if (voxSizes.i == 0.0) {
		indices[0] = 0; indices[1] = 0;
	} else {

		indices[0] = (int) floor(setMin.i / voxSizes.i);
		if (indices[0] >= dims[0]) indices[0] = dims[0] - 1;
		else if (indices[0] < 0) indices[0] = 0;

		indices[1] = (int) floor(setMax.i / voxSizes.i);
		if (indices[1] >= dims[0]) indices[1] = dims[0] - 1;
		else if (indices[1] < 0) indices[1] = 0;
	}
	if (voxSizes.j == 0.0) {
		indices[2] = 0; indices[3] = 0;
	} else {
		indices[2] = (int) floor(setMin.j / voxSizes.j);
		if (indices[2] >= dims[1]) indices[2] = dims[1] - 1;
		else if (indices[2] < 0) indices[2] = 0;

		indices[3] = (int) floor(setMax.j / voxSizes.j);
		if (indices[3] >= dims[1]) indices[3] = dims[1] - 1;
		else if (indices[3] < 0) indices[3] = 0;
	}
	if (voxSizes.k == 0.0) {
		indices[4] = 0; indices[5] = 0;
	} else {
		indices[4] = (int) floor(setMin.k / voxSizes.k);
		if (indices[4] >= dims[2]) indices[4] = dims[2] - 1;
		else if (indices[4] < 0) indices[4] = 0;

		indices[5] = (int) floor(setMax.k / voxSizes.k);
		if (indices[5] >= dims[2]) indices[5] = dims[2] - 1;
		else if (indices[5] < 0) indices[5] = 0;
	}

	return True;
}

/*
 * returns a list of polygon indices that belong to the group of voxels
 * lying within the min and max bounds, and also the total number found.
 *
 * Params	s		SEADS structure
 *			min		minimum corner of box in question
 *			max		maximum corner of box in question
 *			ap		pointer to pointer where the polygon indices will be 
 *						returned
 *
 * Globals	none
 *
 * Returns	the number of polygons in the array, and the array
 */
static int listPolys(s, min, max, ap)
SEADS *s;
Tupple min, max;
int **ap;
{
	int total, i, j, k, q, r, temp1, temp2, index, start[3], stop[3];
	int *ta, *a, type;
	Real *dp, error;
	Boolean found = False;

	total = 0;

	if (*ap != NULL)
		free(*ap);

	/* decrease minimums and increase maximums slightly */
	for (i=0; i < 6; i++) {
		switch (i) {
			case 0: dp = &min.i; type = 1; break;
			case 1: dp = &min.j; type = 1; break;
			case 2: dp = &min.k; type = 1; break;
			case 3: dp = &max.i; type = 2; break;
			case 4: dp = &max.j; type = 2; break;
			case 5: dp = &max.k; type = 2; break;
		}

		if (*dp >= -ETA && *dp <= ETA)
			error = ETA;
		else
			error = fabs(ETA * *dp);
		
		if (type == 1)
			*dp = *dp - error;
		else
			*dp = *dp + error;
	}

	/* set index bounds of voxels to look at */
	if (s->voxSizes.i == 0.0) {
		start[0] = 0; stop[0] = 0;
	} else {
		start[0] = (int) floor(min.i / s->voxSizes.i);
		if (start[0] >= s->dimi) start[0] = s->dimi - 1;
		else if (start[0] < 0) start[0] = 0;

		stop[0] = (int) floor(max.i / s->voxSizes.i);
		if (stop[0] >= s->dimi) stop[0] = s->dimi - 1;
		else if (stop[0] < 0) stop[0] = 0;
	}
	if (s->voxSizes.j == 0.0) {
		start[1] = 0; stop[1] = 0;
	} else {
		start[1] = (int) floor(min.j / s->voxSizes.j);
		if (start[1] >= s->dimj) start[1] = s->dimj - 1;
		else if (start[1] < 0) start[1] = 0;

		stop[1] = (int) floor(max.j / s->voxSizes.j);
		if (stop[1] >= s->dimj) stop[1] = s->dimj - 1;
		else if (stop[1] < 0) stop[1] = 0;
	}
	if (s->voxSizes.k == 0.0) {
		start[2] = 0; stop[2] = 0;
	} else {
		start[2] = (int) floor(min.k / s->voxSizes.k);
		if (start[2] >= s->dimk) start[2] = s->dimk - 1;
		else if (start[2] < 0) start[2] = 0;

		stop[2] = (int) floor(max.k / s->voxSizes.k);
		if (stop[2] >= s->dimk) stop[2] = s->dimk - 1;
		else if (stop[2] < 0) stop[2] = 0;
	}

	a = NULL;
	/* go through all the voxels in question */
	for(k=start[2]; k<=stop[2]; k++) {
		temp2 = k * s->dimi * s->dimj;
		for(j=start[1]; j<=stop[1]; j++) {
			temp1 = j * s->dimi;
			for(i=start[0]; i<=stop[0]; i++) {
				ta = s->matrix[(index=i+temp1+temp2)].a;
				/* go through the list of polygons in the current voxel */
				for (q=0; q < s->matrix[index].total; q++) {
					/* and see if it is already in the list of polygons */
					found = False;
					for (r=0; r < total; r++) {
						if (a[r] == ta[q]) {
							found = True;
							break;
						}
					}
					if (found) continue;

					/* not already in list so add it */
					total++;
					if (total == 1)
						a = (int *)malloc(sizeof(int));
					else
						a = (int *)realloc(a, sizeof(int) * total);
					if (a == NULL) {
						fprintf(stderr, "Error: couldn't realloc int array\n");
						exit(1);
					}
					a[total-1] = ta[q];
				}
			}
		}
	}

	*ap = a;
	return total;
}

/*
 * create/add to PairNode list.  Takes two arrays of integers, and forms pairs 
 * and inserts them into the list if they aren't already in the list
 *
 * Params	a, b	arrays of indices (each array has no duplicates)
 *			ta, tb	number of elements in each array
 *			p		PairNode list to add to (sorted Real linked list)
 *
 * Globals	none
 *
 * Returns	pointer to the head of the PairNode list
 */
static PairNode *addPairNodeList(a, b, ta, tb, p)
int *a, *b, ta, tb;
PairNode *p;
{
	int i, j;
	PairNode *cp, *tp;

	for (i=0; i < ta; i++) {
		for (j=0; j < tb; j++) {
			tp = (PairNode *)malloc(sizeof(PairNode));
			if (tp == NULL) {
				fprintf(stderr, "Error: couldn't malloc PairNode\n");
				exit(1);
			}
			tp->a = a[i]; tp->b = b[j];
			tp->next = NULL;
			tp->prev = NULL;

			if (p == NULL) {
				p = tp;
				continue;
			}

			/* find point to insert this pair in sorted doubly linked list */
			cp = p;
			while (True) {
				if (a[i] == cp->a && b[j] == cp->b) {
					/* this pair of indices already exists */
					free(tp);
					break;
				} else if ((a[i] < cp->a) || (a[i] == cp->a && b[j] < cp->b)) {
					/* insert here */
					if (cp->prev != NULL) cp->prev->next = tp;
					tp->prev = cp->prev;
					cp->prev = tp;
					tp->next = cp;
					break;
				} else if (cp->next == NULL) {
					cp->next = tp;
					tp->prev = cp;
					break;
				} else {
					cp = cp->next;
				}
			}
		}
	}

	return p;
}

/*
 * Given the global minimuns and maximums of two objects, tests to see if the
 * two bounding boxes intersect or not (taking into account some error in the
 * bounding values).
 *
 * Params	mina		minimums of obj A
 *			maxa		maximums of obj A
 *			minb		minimums of obj B
 *			maxb		maximums of obj B
 *
 * Globals	none
 *
 * Returns	True if the global bounding boxes overlap, otherwise false
 */
static Boolean globalIntersect(mina, maxa, minb, maxb)
Tupple mina, maxa, minb, maxb;
{
	int i, type;
	Real *dp, error;

	for (i=0; i < 12; i++) {
		switch (i) {
			case 0: dp = &mina.i; type = 1; break;
			case 1: dp = &mina.j; type = 1; break;
			case 2: dp = &mina.k; type = 1; break;
			case 3: dp = &maxa.i; type = 2; break;
			case 4: dp = &maxa.j; type = 2; break;
			case 5: dp = &maxa.k; type = 2; break;
			case 6: dp = &minb.i; type = 1; break;
			case 7: dp = &minb.j; type = 1; break;
			case 8: dp = &minb.k; type = 1; break;
			case 9: dp = &maxb.i; type = 2; break;
			case 10: dp = &maxb.j; type = 2; break;
			case 11: dp = &maxb.k; type = 2; break;
		}

		if (*dp >= - ETA && *dp <= ETA)
			error = ETA;
		else
			error = fabs(*dp * ETA);

		if (type == 1)
			*dp = *dp - error;
		else
			*dp = *dp + error;
	}

	if ((maxa.i >= minb.i && maxb.i >= mina.i) &&
		(maxa.j >= minb.j && maxb.j >= mina.j) &&
		(maxa.k >= minb.k && maxb.k >= mina.k))
		return True;
	else
		return False;
	
}