quickhull.jai

#scope_file
#import "Basic";
#import "Sort";
btAssert :: assert;

operator == :: (using a: $T, b: T) -> bool 	{
	return (x == b.x) && (y == b.y) && (z == b.z);
}

operator != :: (using a: $T, b: T) -> bool	{
	return (x != b.x) || (y != b.y) || (z != b.z);
}

operator+ :: (using a: $T, b: T) -> T {
    result: T;
	result.x = x + b.x;
	result.y = y + b.y;
	result.z = z + b.z;
    return result;
}

operator- :: (using a: $T, b: T) -> T {
    result: T;
    result.x = x - b.x;
    result.y = y - b.y;
    result.z = z - b.z;
	return result;
}

isZero :: (this: $T) -> bool {
	return (this.x == 0) && (this.y == 0) && (this.z == 0);
}

dot :: (a: $T, b: T) -> type_of(T.x) {
	return a.x * b.x + a.y * b.y + a.z * b.z;
}

dot :: (a: $T, b: $U, $RETURN_T: Type) -> RETURN_T {

    result: RETURN_T;
    result = a.x * b.x;
    result += a.y * b.y;
    result += a.z * b.z;
	return result;
}

// cross :: (a: $T, b: T ) -> T {
//     result: T;
//     result.x = (cast,force(s64)a.y) * b.z - (cast,force(s64)a.z) * b.y;
//     result.y = (cast,force(s64)a.z) * b.x - (cast,force(s64)a.x) * b.z;
//     result.z = (cast,force(s64)a.x) * b.y - (cast,force(s64)a.y) * b.x;
// 	return result;
// }
cross :: (a: $T, b: $U , $RETURN_T: Type) -> RETURN_T {
    result: RETURN_T;
    result.x = (cast(s64)a.y) * b.z - (cast(s64)a.z) * b.y;
    result.y = (cast(s64)a.z) * b.x - (cast(s64)a.x) * b.z;
    result.z = (cast(s64)a.x) * b.y - (cast(s64)a.y) * b.x;
	return result;
}

Point64 :: struct	{
		x: s64;
		y: s64;
		z: s64;
};

Point32 :: struct {
	x: s32;
	y: s32;
	z: s32;
	index: s32 = -1;
};

Int128 :: struct
{
	low: u64;
	high: u64;
};
	get_Int128 :: (low: u64, high: u64)  -> Int128 {
		a: Int128;
		a.low = low;
		a.high = high;
		return a;
	}

	get_Int128 :: (low: u64)-> Int128{
		a: Int128;
		a.low=low;
		return a;
	}

	get_Int128 :: (value: s64) -> Int128  {
        the_val : u64 : 0xffffffffffffffff;
	    return Int128.{low  = cast,no_check(u64)value,
        high = ifx value >= 0 then 0 else the_val };
        // a: Int128;
        // a.low = value;
        // a.high = ifx value >= 0 then 0 else cast,force(u64)(-1);
        // return a;
	}


	operator- :: (using this: Int128) -> Int128
	{
		return get_Int128(cast(u64)( - cast(s64)low), ~high + cast(u64)(low == 0));
	}

	operator+ :: (using this: Int128, b: Int128) -> Int128
	{
		lo :u64= low + b.low;
		return get_Int128(lo, high + b.high + cast,no_check(u64)(lo < low));
	}

    operator_plus_plus :: (a: *u64) {
	    a.* += 1;
	}

	operator_plus_plus :: (a: *Int128) {
	    a.low += 1;
		if (a.low == 0) {
		      a.high += 1;
		}
	}

	operator- :: (this: Int128, b: Int128) -> Int128
	{
		return this + -b;
	}

	operator+= :: (this: Int128, b: Int128) -> Int128
	{
        result: Int128 = this;
		lo :u64= result.low + b.low;
		if (lo < result.low)
		{
			result.high += 1;
		}
		result.low = lo;
		result.high += b.high;
		return result;
	}


	toScalar :: (using this: Int128) -> float
	{
	  hexfloat:float:0x100000000;
	  return ifx (cast,no_check(s64)high >= 0) then cast(float)(high) * (hexfloat * hexfloat) + cast(float)(low) else -(toScalar(-this));
	}

	getSign :: (using this: Int128) -> s32
	{
		return ifx (cast,no_check(s64)high < 0)
		       then cast(s32)-1
		       else ifx (high != 0 || low != 0)
		            then cast(s32)1
		            else cast(s32)0;
	}

	operator< :: (a:Int128, b: Int128) -> bool
	{
		return (a.high < b.high) || ((a.high == b.high) && (a.low < b.low));
	}

	ucmp ::(using a:Int128, b: Int128) -> s32
	{
		if (high < b.high)			{
			return -1;
		}
		if (high > b.high)			{
			return 1;
		}
		if (low < b.low)			{
			return -1;
		}
		if (low > b.low)
		{
			return 1;
		}
		return 0;
	}

get_Rational64 :: (numerator: s64, denominator: s64) -> Rational64 {
    using rat: Rational64;
	if (numerator > 0) {
		sign = 1;
		m_numerator = cast(u64)numerator;
	}
	else if (numerator < 0) {
		sign = -1;
		m_numerator = cast(u64)-numerator;
	} else {
		sign = 0;
		m_numerator = 0;
	}
	if (denominator > 0) {
		m_denominator = cast(u64)denominator;
	} else if (denominator < 0) {
		sign = -sign;
		m_denominator = cast(u64)-denominator;
	} else {
		m_denominator = 0;
	}
	return rat;
}

isNegativeInfinity :: (using rat: Rational64) -> bool {
	return (sign < 0) && (m_denominator == 0);
}

isNaN :: (using rat: Rational64) -> bool {
	return (sign == 0) && (m_denominator == 0);
}

toScalar :: (using rat: Rational64) -> float {
	return sign * (ifx (m_denominator == 0) then FLOAT32_INFINITY else cast(float)m_numerator / m_denominator);
}
toScalar :: (using this: Rational128) -> float {
	return sign * (ifx (getSign(denominator) == 0)  then FLOAT32_INFINITY else toScalar(numerator) / toScalar(denominator));
}

Rational64 :: struct {
	m_numerator: u64;
	m_denominator: u64;
	sign: s32;

};

Rational128 :: struct{
	numerator: Int128;
	denominator: Int128;
	sign: s32;
	isInt64: bool;

};
get_Rational128 :: ( value: s64 ) -> Rational128
{
    result: Rational128;
	if (value > 0) {
		result.sign = 1;
		result.numerator = get_Int128(value);
	}
	else if (value < 0) {
		result.sign = -1;
		result.numerator = get_Int128(-value);
	}
	else
	{
		result.sign = 0;
		result.numerator = get_Int128(cast(u64)0);
	}
	result.denominator = get_Int128(cast(u64)1);
	result.isInt64 = true;
	return result;
}


get_Rational128 :: (numerator: Int128, denominator: Int128) -> Rational128 {
    result: Rational128;
	result.sign = getSign(numerator);
	if (result.sign >= 0) {
		result.numerator = numerator;
	} else {
		result.numerator = -numerator;
	}
	dsign :s32= getSign(denominator);
	if (dsign >= 0)	{
		result.denominator = denominator;
	} else {
		result.sign = -result.sign;
		result.denominator = -denominator;
	}
	result.isInt64 = false;
	return result;
}



xvalue :: (using this: PointR128) -> float {
	return toScalar(x) / toScalar(denominator);
}

yvalue :: (using this: PointR128) -> float {
	return toScalar(y) / toScalar(denominator);
}

zvalue :: (using this: PointR128) -> float {
	return toScalar(z) / toScalar(denominator);
}
PointR128 :: struct
{
	x: Int128;
	y: Int128;
	z: Int128;
	denominator: Int128;


};


operator- :: ( a: Vertex, b: Vertex) -> Point32 {
	return a.point - b.point;
}

dot :: (v: Vertex, b: Point64) -> Rational128 {
	return ifx (v.point.index >= 0) then get_Rational128(dot(v.point, b,s64))
								else get_Rational128(v.point128.x * b.x + v.point128.y * b.y + v.point128.z * b.z, v.point128.denominator);
}

xvalue :: (v: Vertex)  -> float {
	return ifx (v.point.index >= 0) then cast(float)(v.point.x) else xvalue(v.point128);
}

yvalue :: (v: Vertex)  -> float	{
	return ifx (v.point.index >= 0) then cast(float)(v.point.y) else yvalue(v.point128);
}

zvalue :: (v: Vertex)  -> float	{
	return ifx (v.point.index >= 0) then cast(float)(v.point.z) else zvalue(v.point128);
}

receiveNearbyFaces :: (using this: *Vertex, src: *Vertex)
{
	if (lastNearbyFace)		{
		lastNearbyFace.nextWithSameNearbyVertex = src.firstNearbyFace;
	} else {
		firstNearbyFace = src.firstNearbyFace;
	}
	if (src.lastNearbyFace) {
		lastNearbyFace = src.lastNearbyFace;
	}
    // for (Face* f = src.firstNearbyFace; f; f = f.nextWithSameNearbyVertex) {
    f := src.firstNearbyFace;
	while f {
	    defer f = f.nextWithSameNearbyVertex;
		btAssert(f.nearbyVertex == src);
		f.nearbyVertex = this;
	}
	src.firstNearbyFace = null;
	src.lastNearbyFace = null;
}

Vertex :: struct {
	next:            *Vertex;
	prev:            *Vertex;
	edges:           *Edge;
	firstNearbyFace: *Face;
	lastNearbyFace:  *Face;
	point128:        PointR128;
	point:           Point32;
	copy:            s32 = -1;
};

destroy_Edge :: (edge: *Edge) {
	edge.next = null;
	edge.prev = null;
	edge.reverse = null;
	edge.target = null;
	edge.face = null;
}

link :: (using this: *Edge, n: *Edge) {
	btAssert(reverse.target == n.reverse.target);
	next = n;
	n.prev = this;
}

Edge :: struct {
	next:   *Edge;
	prev:   *Edge;
	reverse:*Edge;
	target: *Vertex;
	face:   *Face;
	copy:   s32;
};

init :: (using face: *Face, a: *Vertex, b: *Vertex, c: *Vertex)
{
	nearbyVertex = a;
	origin = a.point;
	dir0 = b.* - a.*;
	dir1 = c.* - a.*;
	if (a.lastNearbyFace) {
		a.lastNearbyFace.nextWithSameNearbyVertex = face;
	} else {
		a.firstNearbyFace = face;
	}
	a.lastNearbyFace = face;
}

getNormal :: (using face: *Face) -> Point64 {
	return cross(dir0, dir1, Point64);
}

Face :: struct{
	next: *Face;
	nearbyVertex: *Vertex;
	nextWithSameNearbyVertex: *Face;
	origin: Point32 ;
	dir0: Point32 ;
	dir1: Point32 ;
};


DMul :: struct(UWord: Type, UHWord: Type)
{
	 high :: (value: u64 ) ->  u32 {
		return cast,no_check(u32)(value >> 32);
	}

	 low :: (value: u64 ) ->  u32 {
		return cast,no_check(u32)value;
	}

	 mul :: (a: u32 , b: u32) ->  u64 {
		return cast(u64)a * cast(u64)b;
	}

	 shlHalf :: ( value: *u64)  {
		value.* <<= 32;
	}

	 high :: (value: Int128 ) ->  u64 {
		return value.high;
	}

	 low :: (value: Int128 ) ->  u64 {
		return value.low;
	}

	 mul :: (a: u64 , b: u64) ->  Int128 {
		return mul(a, b);
	}

	 shlHalf :: ( value: *Int128) {
		value.high = value.low;
		value.low = 0;
	}

	mul :: (a: UWord , b: UWord ,  resLow: *UWord,  resHigh: *UWord)
	{
		p00: UWord = mul(low(a), low(b));
		p01: UWord = mul(low(a), high(b));
		p10: UWord = mul(high(a), low(b));
		p11: UWord = mul(high(a), high(b));
        // p0110: UWord = (low(p01)).(UWord) + (low(p10)).(UWord);
		#if UWord == u64 {
            p0110: UWord = low(p01) + low(p10);
    		p11 += high(p01);
    		p11 += high(p10);
    		p11 += high(p0110);

        } else {
            p0110: UWord = get_Int128(low(p01)) + get_Int128(low(p10));
    		p11 += get_Int128(high(p01));
    		p11 += get_Int128(high(p10));
    		p11 += get_Int128(high(p0110));
        }
		shlHalf(*p0110);
		p00 += p0110;
		if (p00 < p0110)
		{
			operator_plus_plus(*p11);
		}
		resLow.* = p00;
		resHigh.* = p11;
	}
};

IntermediateHull :: struct
{
	minXy: *Vertex;
	maxXy: *Vertex;
	minYx: *Vertex;
	maxYx: *Vertex;
};

Orientation :: enum
{
	NONE;
	CLOCKWISE;
	COUNTER_CLOCKWISE;
};

#scope_export
btEdge :: struct {
		next: int;
		reverse: int;
		targetVertex: int;
}

btConvexHullComputer :: struct {
	vertices_list: [..]Vector3;        // Vertices of the output hull
	original_vertex_index: [..]int; // The original vertex index in the input coords array
	edges_list: [..]btEdge;                // Edges of the output hull
	faces_list: [..]int;                 // Faces of the convex hull. Each entry is an index into the "edges" array pointing to an edge of the face. Faces are planar n-gons
}

btConvexHullInternal :: struct {
    // buffer: Flat_Pool;
    // vertexPool:  [..]Vertex;
    // edgePool:  [..]Edge;
    // facePool:  [..]Face;
	originalVertices: [..]*Vertex;
	vertexList: *Vertex;

	scaling: Vector3;
	center: Vector3;
	mergeStamp: s32;
	minAxis: int;
	medAxis: int;
	maxAxis: int;
	usedEdgePairs: s32;
	maxUsedEdgePairs: s32;

	removeEdgePair :: (using this: *btConvexHullInternal, edge: *Edge)
	{
		n := edge.next;
		r := edge.reverse;

		btAssert(edge.target && r.target);

		if (n != edge) {
			n.prev = edge.prev;
			edge.prev.next = n;
			r.target.edges = n;
		} else {
			r.target.edges = null;
		}

		n = r.next;

		if (n != r) {
			n.prev = r.prev;
			r.prev.next = n;
			edge.target.edges = n;
		} else {
			edge.target.edges = null;
		}
		//I jsut use an edgePool: [..]*Edges; as memory.
        // edgePool.freeObject(edge); Origial
        destroy_Edge(edge);//maybe I need to call this? because jai doesnt have destructors
        free(edge);
        // array_ordered_remove_by_index(*edgePool, edge-edgePool.data);
        // edgePool.freeObject(r); Origial
        destroy_Edge(r);//maybe I need to call this? because jai doesnt have destructors
        free(r);
        // array_ordered_remove_by_index(*edgePool, r-edgePool.data);
		usedEdgePairs-= 1;
	}
};

operator * ::(using this: Int128, b: s64) -> Int128
{
	negative : bool= cast(s64)high < 0;
	a := ifx negative then -this else this;
	if (b < 0)	{
		negative = !negative;
		b = -b;
	}
	result := mul(a.low, cast(u64)b);
	result.high += a.high * cast(u64)b;
	return ifx negative then -result else result;
}

mul :: (a: s64, b: s64) -> Int128 {
	result: Int128;

	negative: bool = a < 0;
	if negative{
		a = -a;
	}

	if b < 0
	{
		negative = !negative;
		b = -b;
	}
	DMul(u64, u32).mul(cast(u64)a, cast(u64)b, *result.low, *result.high);
	return ifx negative then -result else result;
}

mul :: (a: u64, b: u64) -> Int128
{
	result: Int128;
	DMul(u64, u32).mul(a, b, *result.low, *result.high);

	return result;
}

compare :: (using a: Rational64, b: Rational64) -> s32
{
	if (sign != b.sign) {
		return sign - b.sign;
	}
	else if (sign == 0) {
		return 0;
	}
	m1 :Int128 = mul(m_numerator, b.m_denominator);
	m2 :Int128 = mul(m_denominator, b.m_numerator);
	xd := ucmp(m1, m2);
    return sign * xd;

}

compare :: (using a: Rational128, b: Rational128) -> s32
{
	if (sign != b.sign) {
		return sign - b.sign;
	} else if (sign == 0) {
		return 0;
	}

	if (isInt64) {
		return -compare(b, sign * cast(s64)numerator.low);
	}

	nbdLow: Int128; nbdHigh: Int128; dbnLow: Int128; dbnHigh: Int128;
	DMul(Int128, u64).mul(numerator, b.denominator, *nbdLow, *nbdHigh);
	DMul(Int128, u64).mul(denominator, b.numerator, *dbnLow, *dbnHigh);

	cmp := ucmp(nbdHigh,dbnHigh);
	if (cmp) {
		return cmp * sign;
	}
	return ucmp(nbdLow,dbnLow) * sign;
}

compare :: (using this: Rational128, b: s64) -> s32
{
	if (isInt64) {
		a := sign * cast(s64)numerator.low;
		return ifx (a > b) then cast(s32)1 else ifx (a < b) then -1 else cast(s32)0;
	}
	if (b > 0) {
		if (sign <= 0) {
			return -1;
		}
	} else if (b < 0) {
		if (sign >= 0) {
			return 1;
		}
		b = -b;
	} else {
		return sign;
	}

	return ucmp(numerator, denominator * b) * sign;
}

newEdgePair :: (using this: *btConvexHullInternal, from: *Vertex, to: *Vertex) -> *Edge {
	btAssert(from && to);
	e := New(Edge);//edgePool.newObject();
	r := New(Edge);//edgePool.newObject();
	e.reverse = r;
	r.reverse = e;
	e.copy = mergeStamp;
	r.copy = mergeStamp;
	e.target = to;
	r.target = from;
	e.face = null;
	r.face = null;
	usedEdgePairs += 1;
	if (usedEdgePairs > maxUsedEdgePairs) {
		maxUsedEdgePairs = usedEdgePairs;
	}
	return e;
}

mergeProjection :: (h0: *IntermediateHull, h1: *IntermediateHull, c0: **Vertex, c1: **Vertex) -> bool
{
	v0 := h0.maxYx;
	v1 := h1.minYx;
	if ((v0.point.x == v1.point.x) && (v0.point.y == v1.point.y)) {
		btAssert(v0.point.z < v1.point.z);
		v1p := v1.prev;
		if (v1p == v1) {
			c0.* = v0;
			if (v1.edges) {
				btAssert(v1.edges.next == v1.edges);
				v1 = v1.edges.target;
				btAssert(v1.edges.next == v1.edges);
			}
			c1.* = v1;
			return false;
		}
		v1n := v1.next;
		v1p.next = v1n;
		v1n.prev = v1p;
		if (v1 == h1.minXy) {
			if (v1n.point.x < v1p.point.x) || ((v1n.point.x == v1p.point.x) && (v1n.point.y < v1p.point.y)) {
				h1.minXy = v1n;
			} else {
				h1.minXy = v1p;
			}
		}
		if (v1 == h1.maxXy) {
			if (v1n.point.x > v1p.point.x) || ((v1n.point.x == v1p.point.x) && (v1n.point.y > v1p.point.y)) {
				h1.maxXy = v1n;
			} else {
				h1.maxXy = v1p;
			}
		}
	}

	v0 = h0.maxXy;
	v1 = h1.maxXy;
	v00 :*Vertex;
	v10 :*Vertex;
	sign: s32 = 1;

    // for (s32 side = 0; side <= 1; side++)
	for side:  0..1
	{
	    if side == 1 btAssert(sign == -1);
		dx:s32 = (v1.point.x - v0.point.x) * sign;
		if (dx > 0) {
			while (true) {
				dy:s32 = v1.point.y - v0.point.y;

				w0 := ifx side then v0.next else v0.prev;
				if w0 != v0 {
					dx0:s32 = (w0.point.x - v0.point.x) * sign;
					dy0:s32 = w0.point.y - v0.point.y;
					if (dy0 <= 0) && ((dx0 == 0) || ((dx0 < 0) && (dy0 * dx <= dy * dx0))) {
						v0 = w0;
						dx = (v1.point.x - v0.point.x) * sign;
						continue;
					}
				}

				w1 := ifx side then v1.next else v1.prev;
				if (w1 != v1) {
					dx1:s32 = (w1.point.x - v1.point.x) * sign;
					dy1:s32 = w1.point.y - v1.point.y;
					dxn:s32 = (w1.point.x - v0.point.x) * sign;
					if ((dxn > 0) && (dy1 < 0) && ((dx1 == 0) || ((dx1 < 0) && (dy1 * dx < dy * dx1)))) {
						v1 = w1;
						dx = dxn;
						continue;
					}
				}

				break;
			}
		}
		else if (dx < 0) {
			while (true)
			{
				 dy: = v1.point.y - v0.point.y;

				w1 := ifx side  then v1.prev else v1.next;
				if (w1 != v1) {
					 dx1: = (w1.point.x - v1.point.x) * sign;
					 dy1: = w1.point.y - v1.point.y;
					if ((dy1 >= 0) && ((dx1 == 0) || ((dx1 < 0) && (dy1 * dx <= dy * dx1))))					{
						v1 = w1;
						dx = (v1.point.x - v0.point.x) * sign;
						continue;
					}
				}

				w0 := ifx side then v0.prev else v0.next;
				if (w0 != v0)				{
					 dx0: = (w0.point.x - v0.point.x) * sign;
					 dy0: = w0.point.y - v0.point.y;
					 dxn: = (v1.point.x - w0.point.x) * sign;
					if ((dxn < 0) && (dy0 > 0) && ((dx0 == 0) || ((dx0 < 0) && (dy0 * dx < dy * dx0))))					{
						v0 = w0;
						dx = dxn;
						continue;
					}
				}

				break;
			}
		} else {
			x: = v0.point.x;
			y0: = v0.point.y;
			w0 := v0;
			t: *Vertex;
			t = ifx side then w0.next else w0.prev;
        // while (((t = side ? w0->next : w0->prev) != v0) && (t->point.x == x) && (t->point.y <= y0))
			while ((t != v0) && (t.point.x == x) && (t.point.y <= y0))
			{
				w0 = t;
				y0 = t.point.y;
			    t = ifx side then w0.next else w0.prev;
			}
			v0 = w0;

			y1 := v1.point.y;
			w1 := v1;
			t = ifx side then w1.prev else w1.next;
        // while (((t = side ? w1->prev : w1->next) != v1) && (t->point.x == x) && (t->point.y >= y1))
			while ((t != v1) && (t.point.x == x) && (t.point.y >= y1))
			{
				w1 = t;
				y1 = t.point.y;
			    t = ifx side then w1.prev else w1.next;
			}
			v1 = w1;
		}

		if (side == 0) {
			v00 = v0;
			v10 = v1;

			v0 = h0.minXy;
			v1 = h1.minXy;
			sign = -1;
		}
	}

	v0.prev = v1;
	v1.next = v0;

	v00.next = v10;
	v10.prev = v00;

	if (h1.minXy.point.x < h0.minXy.point.x)
	{
		h0.minXy = h1.minXy;
	}
	if (h1.maxXy.point.x >= h0.maxXy.point.x)
	{
		h0.maxXy = h1.maxXy;
	}

	h0.maxYx = h1.maxYx;

	c0.* = v00;
	c1.* = v10;

	return true;
}
/*
bool btConvexHullInternal::mergeProjection(IntermediateHull& h0, IntermediateHull& h1, Vertex*& c0, Vertex*& c1)
{
	Vertex* v0 = h0.maxYx;
	Vertex* v1 = h1.minYx;
	if ((v0->point.x == v1->point.x) && (v0->point.y == v1->point.y))
	{
		btAssert(v0->point.z < v1->point.z);
		Vertex* v1p = v1->prev;
		if (v1p == v1)
		{
			c0 = v0;
			if (v1->edges)
			{
				btAssert(v1->edges->next == v1->edges);
				v1 = v1->edges->target;
				btAssert(v1->edges->next == v1->edges);
			}
			c1 = v1;
			return false;
		}
		Vertex* v1n = v1->next;
		v1p->next = v1n;
		v1n->prev = v1p;
		if (v1 == h1.minXy)
		{
			if ((v1n->point.x < v1p->point.x) || ((v1n->point.x == v1p->point.x) && (v1n->point.y < v1p->point.y)))
			{
				h1.minXy = v1n;
			}
			else
			{
				h1.minXy = v1p;
			}
		}
		if (v1 == h1.maxXy)
		{
			if ((v1n->point.x > v1p->point.x) || ((v1n->point.x == v1p->point.x) && (v1n->point.y > v1p->point.y)))
			{
				h1.maxXy = v1n;
			}
			else
			{
				h1.maxXy = v1p;
			}
		}
	}

	v0 = h0.maxXy;
	v1 = h1.maxXy;
	Vertex* v00 = NULL;
	Vertex* v10 = NULL;
	int32_t sign = 1;

	for (int side = 0; side <= 1; side++)
	{
		int32_t dx = (v1->point.x - v0->point.x) * sign;
		if (dx > 0)
		{
			while (true)
			{
				int32_t dy = v1->point.y - v0->point.y;

				Vertex* w0 = side ? v0->next : v0->prev;
				if (w0 != v0)
				{
					int32_t dx0 = (w0->point.x - v0->point.x) * sign;
					int32_t dy0 = w0->point.y - v0->point.y;
					if ((dy0 <= 0) && ((dx0 == 0) || ((dx0 < 0) && (dy0 * dx <= dy * dx0))))
					{
						v0 = w0;
						dx = (v1->point.x - v0->point.x) * sign;
						continue;
					}
				}

				Vertex* w1 = side ? v1->next : v1->prev;
				if (w1 != v1)
				{
					int32_t dx1 = (w1->point.x - v1->point.x) * sign;
					int32_t dy1 = w1->point.y - v1->point.y;
					int32_t dxn = (w1->point.x - v0->point.x) * sign;
					if ((dxn > 0) && (dy1 < 0) && ((dx1 == 0) || ((dx1 < 0) && (dy1 * dx < dy * dx1))))
					{
						v1 = w1;
						dx = dxn;
						continue;
					}
				}

				break;
			}
		}
		else if (dx < 0)
		{
			while (true)
			{
				int32_t dy = v1->point.y - v0->point.y;

				Vertex* w1 = side ? v1->prev : v1->next;
				if (w1 != v1)
				{
					int32_t dx1 = (w1->point.x - v1->point.x) * sign;
					int32_t dy1 = w1->point.y - v1->point.y;
					if ((dy1 >= 0) && ((dx1 == 0) || ((dx1 < 0) && (dy1 * dx <= dy * dx1))))
					{
						v1 = w1;
						dx = (v1->point.x - v0->point.x) * sign;
						continue;
					}
				}

				Vertex* w0 = side ? v0->prev : v0->next;
				if (w0 != v0)
				{
					int32_t dx0 = (w0->point.x - v0->point.x) * sign;
					int32_t dy0 = w0->point.y - v0->point.y;
					int32_t dxn = (v1->point.x - w0->point.x) * sign;
					if ((dxn < 0) && (dy0 > 0) && ((dx0 == 0) || ((dx0 < 0) && (dy0 * dx < dy * dx0))))
					{
						v0 = w0;
						dx = dxn;
						continue;
					}
				}

				break;
			}
		}
		else
		{
			int32_t x = v0->point.x;
			int32_t y0 = v0->point.y;
			Vertex* w0 = v0;
			Vertex* t;
			while (((t = side ? w0->next : w0->prev) != v0) && (t->point.x == x) && (t->point.y <= y0))
			{
				w0 = t;
				y0 = t->point.y;
			}
			v0 = w0;

			int32_t y1 = v1->point.y;
			Vertex* w1 = v1;
			while (((t = side ? w1->prev : w1->next) != v1) && (t->point.x == x) && (t->point.y >= y1))
			{
				w1 = t;
				y1 = t->point.y;
			}
			v1 = w1;
		}

		if (side == 0)
		{
			v00 = v0;
			v10 = v1;

			v0 = h0.minXy;
			v1 = h1.minXy;
			sign = -1;
		}
	}

	v0->prev = v1;
	v1->next = v0;

	v00->next = v10;
	v10->prev = v00;

	if (h1.minXy->point.x < h0.minXy->point.x)
	{
		h0.minXy = h1.minXy;
	}
	if (h1.maxXy->point.x >= h0.maxXy->point.x)
	{
		h0.maxXy = h1.maxXy;
	}

	h0.maxYx = h1.maxYx;

	c0 = v00;
	c1 = v10;

	return true;
}
*/
computeInternal ::(using this: *btConvexHullInternal, start: s32,  end: s32,  result: *IntermediateHull)
{
	n := end - start;
	if (n) ==
	{
		case 0;
			result.minXy = null;
			result.maxXy = null;
			result.minYx = null;
			result.maxYx = null;
			return;
		case 2;
		{
			v := originalVertices[start];
			w := originalVertices[start+1];
			if (v.point != w.point)
			{
				dx := v.point.x - w.point.x;
				dy := v.point.y - w.point.y;

				if ((dx == 0) && (dy == 0))
				{
					if (v.point.z > w.point.z)
					{
						t := w;
						w = v;
						v = t;
					}
					btAssert(v.point.z < w.point.z);
					v.next = v;
					v.prev = v;
					result.minXy = v;
					result.maxXy = v;
					result.minYx = v;
					result.maxYx = v;
				}
				else
				{
					v.next = w;
					v.prev = w;
					w.next = v;
					w.prev = v;

					if ((dx < 0) || ((dx == 0) && (dy < 0)))
					{
						result.minXy = v;
						result.maxXy = w;
					}
					else
					{
						result.minXy = w;
						result.maxXy = v;
					}

					if ((dy < 0) || ((dy == 0) && (dx < 0)))
					{
						result.minYx = v;
						result.maxYx = w;
					}
					else
					{
						result.minYx = w;
						result.maxYx = v;
					}
				}

				e := newEdgePair(this, v, w);
				link(e,e);
				v.edges = e;

				e = e.reverse;
				link(e,e);
				w.edges = e;

				return;
			}
			{
				 v := originalVertices[start];
				v.edges = null;
				v.next = v;
				v.prev = v;

				result.minXy = v;
				result.maxXy = v;
				result.minYx = v;
				result.maxYx = v;
			}

			return;
		}

		case 1;
		{
			v: = originalVertices[start];
			v.edges = null;
			v.next = v;
			v.prev = v;

			result.minXy = v;
			result.maxXy = v;
			result.minYx = v;
			result.maxYx = v;

			return;
		}
	}

	split0:s32 = start + n / 2;
	p: = originalVertices[split0 - 1].point;
	split1:s32 = split0;
	while ((split1 < end) && (originalVertices[split1].point == p))
	{
		split1+=1;
	}
	computeInternal(this, start, split0, result);
	temp_hull: IntermediateHull;
	computeInternal(this, split1, end, *temp_hull);
	merge(this, result, *temp_hull);
}


getOrientation :: (prev: *Edge, next: *Edge, s: Point32, t: Point32) -> Orientation
{
	btAssert(prev.reverse.target == next.reverse.target);
	if (prev.next == next)
	{
		if (prev.prev == next)
		{
			 n :Point64 = cross(t, s, Point64);
			 m :Point64 = cross((prev.target.* - next.reverse.target.*),next.target.* - next.reverse.target.*, Point64);
			btAssert(!isZero(m));
			component := dot(n, m);
			btAssert(component != 0);
			return ifx (component > 0) then .COUNTER_CLOCKWISE else .CLOCKWISE;
		}
		return .COUNTER_CLOCKWISE;
	}
	else if (prev.prev == next)
	{
		return .CLOCKWISE;
	}
	else
	{
		return .NONE;
	}
}

findMaxAngle :: (using this: *btConvexHullInternal, ccw: bool, start: *Vertex, s: Point32, rxs: Point64, sxrxs: Point64, minCot: *Rational64) -> *Edge {
	minEdge: *Edge;
	e := start.edges;
	if (e) {
		while true {
		    defer if (e == start.edges) break;

			if (e.copy > mergeStamp) {
				t := e.target.* - start.*;
				cot := get_Rational64(dot(t, sxrxs, s64), dot(t, rxs, s64));
				if (isNaN(cot)) {
					btAssert(ifx ccw then (dot(t,s, s64) < 0) else (dot(t,s, s64) > 0));
				} else {
					cmp: s32;
					if (minEdge == null) {
						minCot.* = cot;
						minEdge = e;
					} else {
					  cmp = compare(cot, minCot.*);
					  if cmp < 0 {
					   minCot.* = cot;
					   minEdge = e;
					  } else if ((cmp == 0) && (ccw == (getOrientation(minEdge, e, s, t) == .COUNTER_CLOCKWISE))) {
					   minEdge = e;
					  }
				  }
				}
			}
			e = e.next;
		};
	}
	return minEdge;
}

findEdgeForCoplanarFaces :: (using this: *btConvexHullInternal, c0: *Vertex, c1: *Vertex, e0: *Edge, e1: *Edge, stop0: *Vertex, stop1: *Vertex)
{
	start0: *Edge = e0;
	start1: *Edge = e1;
	et0:Point32 = ifx start0 then start0.target.point else c0.point;
	et1:Point32 = ifx start1 then start1.target.point else c1.point;
	s:Point32 = c1.point - c0.point;
    // Point64 normal = ((start0 ? start0 : start1)->target->point - c0->point).cross(s);
	normal:Point64 = cross(((ifx start0 then start0 else start1).target.point - c0.point),s, Point64);
	 dist:s64 = dot(c0.point,normal,s64);
	btAssert(!start1 || dot(start1.target.point,normal,s64) == dist);
	perp:Point64 = cross(s,normal, Point64);
	btAssert(!isZero(perp));


	maxDot0:s64 = dot(et0,perp, s64);
	if (e0)
	{
		while (e0.target != stop0)
		{
			e: = e0.reverse.prev;
			if (dot(e.target.point, normal, s64) < dist)
			{
				break;
			}
			btAssert(dot(e.target.point, normal, s64) == dist);
			if (e.copy == mergeStamp)
			{
				break;
			}
			component := dot(e.target.point, perp, s64);
			if (component <= maxDot0)
			{
				break;
			}
			maxDot0 = component;
			e0 = e;
			et0 = e.target.point;
		}
	}

	maxDot1: s64 = dot(et1, perp, s64);
	if (e1)
	{
		while (e1.target != stop1)
		{
			e: = e1.reverse.next;
			if (dot(e.target.point, normal, s64) < dist)
			{
				break;
			}
			btAssert(dot(e.target.point, normal, s64) == dist);
			if (e.copy == mergeStamp)
			{
				break;
			}
			component := dot(e.target.point, perp, s64);
			if (component <= maxDot1)
			{
				break;
			}
			maxDot1 = component;
			e1 = e;
			et1 = e.target.point;
		}
	}

	dx := maxDot1 - maxDot0;
	if (dx > 0)
	{
		while (true)
		{
			dy: = dot((et1 - et0),s);

			if (e0 && (e0.target != stop0))
			{
				 f0 := e0.next.reverse;
				if (f0.copy > mergeStamp)
				{
					dx0:s64 = dot((f0.target.point - et0), perp, s64);
					dy0:s64 = dot((f0.target.point - et0), s, s64);
					if (ifx (dx0 == 0) then (dy0 < 0) else ((dx0 < 0) && (compare(get_Rational64(dy0, dx0),get_Rational64(dy, dx)) >= 0)))
					{
						et0 = f0.target.point;
						dx = dot((et1 - et0),perp, s64);
						e0 = ifx (e0 == start0) then null else f0;
						continue;
					}
				}
			}

			if (e1 && (e1.target != stop1))
			{
				f1 := e1.reverse.next;
				if (f1.copy > mergeStamp)
				{
					d1 := f1.target.point - et1;
					if (dot(d1,normal, s64) == 0)
					{
						dx1:s64 = dot(d1,perp, s64);
						dy1:s64 = dot(d1,s, s64);
						dxn:s64 = dot((f1.target.point - et0),perp, s64);
                    // if (cond1 ? (dy1 < 0) : (c1 && (Rational64(dy1, dx1).compare(Rational64(dy, dx)) > 0))))
                    // if ((dxn > 0) && ((dx1 == 0) ? (dy1 < 0) : ((dx1 < 0) && (Rational64(dy1, dx1).compare(Rational64(dy, dx)) > 0))))

						if ((dxn > 0) && (ifx(dx1 == 0) then (dy1 < 0) else ((dx1 < 0) && (compare(get_Rational64(dy1, dx1),get_Rational64(dy, dx)) > 0))))
						{
							e1 = f1;
							et1 = e1.target.point;
							dx = dxn;
							continue;
						}
					}
					else
					{
						btAssert((e1 == start1) && (dot(d1, normal, s64) < 0));
					}
				}
			}

			break;
		}
	}
	else if (dx < 0)
	{
		while (true)
		{
			dy : s64= dot((et1 - et0),s, s64);

			if (e1 && (e1.target != stop1))
			{
				f1 := e1.prev.reverse;
				if (f1.copy > mergeStamp)
				{
					dx1: = dot((f1.target.point - et1),perp, s64);
					dy1: = dot((f1.target.point - et1),s, s64);
					if (ifx (dx1 == 0) then (dy1 > 0) else ((dx1 < 0) && (compare(get_Rational64(dy1, dx1), get_Rational64(dy, dx)) <= 0)))
					{
						et1 = f1.target.point;
						dx = dot((et1 - et0),perp, s64);
						e1 = ifx (e1 == start1)  then null else f1;
						continue;
					}
				}
			}

			if (e0 && (e0.target != stop0))
			{
				f0 := e0.reverse.prev;
				if (f0.copy > mergeStamp)
				{
					d0 := f0.target.point - et0;
					if (dot(d0, normal, s64) == 0)
					{
						dx0:s64 = dot(d0, perp, s64);
						dy0:s64 = dot(d0, s, s64);
						dxn:s64 = dot((et1 - f0.target.point), perp, s64);
						if ((dxn < 0) && (ifx (dx0 == 0) then (dy0 > 0) else ((dx0 < 0) && (compare(get_Rational64(dy0, dx0),get_Rational64(dy, dx)) < 0))))
						{
							e0 = f0;
							et0 = e0.target.point;
							dx = dxn;
							continue;
						}
					} else {
						btAssert((e0 == start0) && (dot(d0, normal, s64) < 0));
					}
				}
			}

			break;
		}
	}
}

merge :: (using this: *btConvexHullInternal, h0: *IntermediateHull, h1: *IntermediateHull) {
	if (!h1.maxXy)	{
		return;
	}
	if (!h0.maxXy)	{
		h0 = h1;
		return;
	}

	mergeStamp-=1;

	c0: *Vertex;
	toPrev0: *Edge;
	firstNew0: *Edge;
	pendingHead0: *Edge;
	pendingTail0: *Edge;
	c1: *Vertex;
	toPrev1: *Edge;
	firstNew1: *Edge;
	pendingHead1: *Edge;
	pendingTail1: *Edge;
	prevPoint: Point32;

	if (mergeProjection(h0, h1, *c0, *c1))
	{
		s: = c1.* - c0.*;
		normal :Point64 = cross(Point32.{0, 0, -1,-1},s,Point64);
		print("c0 = %\n", c0.point);
        print("c1 = %\n", c1.point);
        print("s  = %\n", s);
        print("normal = %\n", normal);
		t := cross(s,normal,Point64);
		btAssert(!isZero(t));

		e: = c0.edges;
		start0 : *Edge;
		if (e)
		{
			while true {
			      defer if (e == c0.edges) break;
				component := dot((e.target.* - c0.*), normal, s64);
				btAssert(component <= 0, "The expectation is that for every edge leaving c0, every neighboring vertex lies on or behind the supporting line/plane defined by normal");
				if ((component == 0) && (dot((e.target.* - c0.*),t,s64) > 0))
				{
					if (!start0 || (getOrientation(start0, e, s, Point32.{0, 0, -1,-1}) == .CLOCKWISE))
					{
						start0 = e;
					}
				}
				e = e.next;
			}
		}

		e = c1.edges;
		start1: *Edge;
		if (e)
		{
			while true { defer if (e == c1.edges) break;

				component :s64= dot((e.target.* - c1.*),normal, s64);
				btAssert(component <= 0);
				if ((component == 0) && (dot((e.target.* - c1.*),t, s64) > 0))
				{
					if (!start1 || (getOrientation(start1, e, s, Point32.{0, 0, -1,-1}) == .COUNTER_CLOCKWISE))
					{
						start1 = e;
					}
				}
				e = e.next;
			} ;
		}

		if (start0 || start1)
		{
			findEdgeForCoplanarFaces(this, c0, c1, start0, start1, null, null);
			if (start0)
			{
				c0 = start0.target;
			}
			if (start1)
			{
				c1 = start1.target;
			}
		}

		prevPoint = c1.point;
		prevPoint.z+=1;
	}
	else
	{
		prevPoint = c1.point;
		prevPoint.x+=1;
	}

	first0 := c0;
	first1 := c1;
	firstRun := true;

	while (true)
	{
		 s: = c1.* - c0.*;
		r:Point32 = prevPoint - c0.point;
		rxs:Point64 = cross(r,s,Point64);
		sxrxs:Point64 = cross(s,rxs,Point64);

		minCot0 := get_Rational64 (0, 0);
		min0 := findMaxAngle(this,false, c0, s, rxs, sxrxs, *minCot0);
		minCot1 := get_Rational64 (0, 0);
		min1 := findMaxAngle(this,true, c1, s, rxs, sxrxs, *minCot1);
		if (!min0 && !min1)
		{
			e: = newEdgePair(this, c0, c1);
			link(e,e);
			c0.edges = e;

			e = e.reverse;
			link(e,e);
			c1.edges = e;
			return;
		}
		else
		{
			cmp :s32= ifx !min0 then 1 else ifx !min1  then -1 else compare(minCot0, minCot1);
			if (firstRun || (ifx(cmp >= 0) then !isNegativeInfinity(minCot1) else !isNegativeInfinity(minCot0)))
			{
				e: = newEdgePair(this, c0, c1);
				if (pendingTail0)
				{
					pendingTail0.prev = e;
				}
				else
				{
					pendingHead0 = e;
				}
				e.next = pendingTail0;
				pendingTail0 = e;

				e = e.reverse;
				if (pendingTail1)
				{
					pendingTail1.next = e;
				}
				else
				{
					pendingHead1 = e;
				}
				e.prev = pendingTail1;
				pendingTail1 = e;
			}

			e0 := min0;
			e1 := min1;


			if (cmp == 0)
			{
				findEdgeForCoplanarFaces(this, c0, c1, e0, e1, null, null);
			}

			if ((cmp >= 0) && e1)
			{
				if (toPrev1)
				{
                    // for (Edge *e = toPrev1.next, *n = null; e != min1; e = n)
                    e := toPrev1.next;
                    n: *Edge;
					while e != min1 {
					   defer e = n;
						n = e.next;
						removeEdgePair(this, e);
					}
				}

				if (pendingTail1)
				{
					if (toPrev1)
					{
						link(toPrev1,pendingHead1);
					}
					else
					{
						link(min1.prev,pendingHead1);
						firstNew1 = pendingHead1;
					}
					link(pendingTail1,min1);
					pendingHead1 = null;
					pendingTail1 = null;
				}
				else if (!toPrev1)
				{
					firstNew1 = min1;
				}

				prevPoint = c1.point;
				c1 = e1.target;
				toPrev1 = e1.reverse;
			}

			if ((cmp <= 0) && e0)
			{
				if (toPrev0)
				{
                    // for (Edge *e = toPrev0.prev, *n = null; e != min0; e = n)
                    e := toPrev0.prev;
                    n: *Edge;
					while  e != min0 {
					    defer e = n;
						n = e.prev;
						removeEdgePair(this, e);
					}
				}

				if (pendingTail0)
				{
					if (toPrev0)
					{
						link(pendingHead0,toPrev0);
					}
					else
					{
						link(pendingHead0,min0.next);
						firstNew0 = pendingHead0;
					}
					link(min0,pendingTail0);
					pendingHead0 = null;
					pendingTail0 = null;
				}
				else if (!toPrev0)
				{
					firstNew0 = min0;
				}

				prevPoint = c0.point;
				c0 = e0.target;
				toPrev0 = e0.reverse;
			}
		}

		if ((c0 == first0) && (c1 == first1))
		{
			if (toPrev0 == null)
			{
				link(pendingHead0,pendingTail0);
				c0.edges = pendingTail0;
			}
			else
			{
                // for (Edge *e = toPrev0.prev, *n = null; e != firstNew0; e = n)
                e := toPrev0.prev;
                n: *Edge;
                while e != firstNew0
				{
				    defer e=n;
					n = e.prev;
					removeEdgePair(this, e);
				}
				if (pendingTail0)
				{
					link(pendingHead0,toPrev0);
					link(firstNew0,pendingTail0);
				}
			}

			if (toPrev1 == null)
			{
				link(pendingTail1,pendingHead1);
				c1.edges = pendingTail1;
			}
			else
			{
                // for (Edge *e = toPrev1.next, *n = null; e != firstNew1; e = n)
                e := toPrev1.next;
                n: *Edge;
                while e != firstNew1
				{
				    defer e=n;
					n = e.next;
					removeEdgePair(this, e);
				}
				if (pendingTail1)
				{
					link(toPrev1, pendingHead1);
					link(pendingTail1, firstNew1);
				}
			}

			return;
		}

		firstRun = false;
	}
}

pointCmp :: (a: Point32, b: Point32) -> int {
    if a.y != b.y return cast(s64)a.y - cast(s64)b.y;
    if a.x != b.x return cast(s64)a.x - cast(s64)b.x;
    if a.z != b.z return cast(s64)a.z - cast(s64)b.z;
    return 0;
}
// pointCmp :: (p: Point32, q: Point32) -> int {
// 	return cast(int)((p.y < q.y) || ((p.y == q.y) && ((p.x < q.x) || ((p.x == q.x) && (p.z < q.z)))));
// }

// pointCmp :: (p: Point32, q: Point32) -> int {
// 	return ifx ((p.y < q.y) || ((p.y == q.y) && ((p.x < q.x) || ((p.x == q.x) && (p.z < q.z))))) then -1 else 1;
// }
get_minAxis :: (a: Vector3) -> int {
	y_is_bigger_than_x := a.x < a.y;
	if y_is_bigger_than_x
		return (ifx a.x < a.z then 0 else 2);
	else
		return (ifx a.y < a.z then 1 else 2);
}

	/**@brief Return the axis with the largest value
   * Note return values are 0,1,2 for x, y, or z */
get_maxAxis :: (a: Vector3) -> int	{
	y_is_bigger_than_x := a.x < a.y;
	if y_is_bigger_than_x
		return (ifx a.y < a.z then 2 else 1);
	else
		return (ifx a.x < a.z then 2 else 0);
}

compute :: (using hull: *btConvexHullInternal, cloud_of_points: []Vector3)
{
	min_v3: Vector3= .{FLOAT32_INFINITY,FLOAT32_INFINITY,FLOAT32_INFINITY};
	max_v3: Vector3= .{ -FLOAT32_INFINITY, -FLOAT32_INFINITY, -FLOAT32_INFINITY};

	for cloud_of_points	{
		min_v3 = min(it, min_v3);
		max_v3 = max(it, max_v3);
	}
	center = (min_v3 + max_v3) * 0.5;
	s := max_v3 - min_v3;
	maxAxis = get_maxAxis(s);
	minAxis = get_minAxis(s);
	if (minAxis == maxAxis)	{
		minAxis = (maxAxis + 1) % 3;
	}
	medAxis = 3 - maxAxis - minAxis;
    assert(medAxis != maxAxis, "minmax same");
    assert(medAxis != minAxis, "minmax same");
    assert(maxAxis != minAxis, "minmax same");
    //above works checked 1 time.
	s /= 10216.0;
	if ((medAxis + 1) % 3) != maxAxis {
		s *= -1;
	}
	scaling = s;

	if s.component[0] != 0 {
		s.component[0] = 1.0/s.component[0];
	}
	if s.component[1] != 0 {
		s.component[1] = 1.0/s.component[1];
	}
	if s.component[2] != 0 {
		s.component[2] = 1.0/s.component[2];
	}


	points: [..]Point32;
	array_reserve(*points, cloud_of_points.count);
	points.count = cloud_of_points.count;
	{
		for cloud_of_points {
			p := it;
			p = (p - center) * s;
			points[it_index].x = cast(s32)p.component[medAxis];
			points[it_index].y = cast(s32)p.component[maxAxis];
			points[it_index].z = cast(s32)p.component[minAxis];
            // assert(p.y > p.x, "1");
            // assert(p.x > p.z, "2 % %", p.x,points[it_index]);
            // assert(p.y > p.z, "3");
			points[it_index].index = cast(s32)it_index;
		}
	}
	quick_sort(points, f=pointCmp);

	array_reserve(*originalVertices, cloud_of_points.count);
	originalVertices.count = cloud_of_points.count;

	for cloud_of_points
	{
		v := New(Vertex);
		v.edges = null;
		v.point = points[it_index];
		v.copy = -1;
		originalVertices[it_index] = v;
	}

	array_reset(*points);
	usedEdgePairs = 0;
	maxUsedEdgePairs = 0;
	mergeStamp = -3;

	temp_hull: IntermediateHull;
	computeInternal(hull, 0, cast(s32)cloud_of_points.count, *temp_hull);
	vertexList = temp_hull.minXy;
}

toBtVector :: (using this: *btConvexHullInternal, v: Point32) -> Vector3
{
	p: Vector3;
	p.component[medAxis] = cast,no_check(float)(v.x);
	p.component[maxAxis] = cast,no_check(float)(v.y);
	p.component[minAxis] = cast,no_check(float)(v.z);
	return p * scaling;
}

getBtNormal :: (using this: *btConvexHullInternal, face: *Face) -> Vector3{
	a := toBtVector(this, face.dir0);
	b := cross(a, toBtVector(this, face.dir1));
	return normalize(b);
}

getCoordinates :: (using h1: *btConvexHullInternal, v: *Vertex) -> Vector3{
	p: Vector3;
	p.component[medAxis] = xvalue(v);
	p.component[maxAxis] = yvalue(v);
	p.component[minAxis] = zvalue(v);
	return p * scaling + center;
}

shrink :: (using this: *btConvexHullInternal, amount: float, clampAmount: float) -> float
{
	if (!vertexList)
	{
		return 0;
	}
	mergeStamp -= 1;
	stamp := mergeStamp;
	stack: [..]*Vertex;
	vertexList.copy = stamp;
	array_add(*stack, vertexList);
	faces: [..]*Face;

	ref := vertexList.point;
	hullCenterX: Int128;
	hullCenterY: Int128;
	hullCenterZ: Int128;
	volume: Int128;

	while (stack.count > 0)
	{
		v := stack[stack.count - 1];
		pop(*stack);
		e := v.edges;
		if (e)
		{
			while true {
			     defer if (e == v.edges) break;
				if (e.target.copy != stamp)
				{
					e.target.copy = stamp;
					array_add(*stack, e.target);
				}
				if (e.copy != stamp)
				{
					face := New(Face);
					init(face,e.target, e.reverse.prev.target, v);
					array_add(*faces, face);
					f :*Edge = e;

					a :*Vertex = null;
					b :*Vertex = null;
					while true
					{
						defer if (f == e) break;
						if (a && b) {
                        // int64_t vol = (v.point - ref).dot((a.point - ref).cross(b.point - ref));


							vol : s64= dot((v.point - ref),cross((a.point - ref), b.point - ref, Point64),s64);
							btAssert(vol >= 0);
							c: Point32= v.point + a.point + b.point + ref;
							vol_128 := get_Int128(vol);
							volume += vol_128;
							hullCenterX += vol_128 * c.x;
							hullCenterY += vol_128 * c.y;
							hullCenterZ += vol_128 * c.z;
						}

						btAssert(f.copy != stamp);
						f.copy = stamp;
						f.face = face;

						a = b;
						b = f.target;

						f = f.reverse.prev;
					} ;
				}
				e = e.next;
			} ;
		}
	}

	if (getSign(volume) <= 0)
	{
		return 0;
	}

	hullCenter: Vector3;
	hullCenter.component[medAxis] = toScalar(hullCenterX);
	hullCenter.component[maxAxis] = toScalar(hullCenterY);
	hullCenter.component[minAxis] = toScalar(hullCenterZ);
	hullCenter /= 4 * toScalar(volume);
	hullCenter *= scaling;

	faceCount: s32 = cast(s32)faces.count;

	if (clampAmount > 0)
	{
		minDist := FLOAT32_INFINITY;
		for i: 0..faceCount - 1 {
			normal := getBtNormal(this, faces[i]);
			dist: float = dot(normal, toBtVector(this, faces[i].origin) - hullCenter);
			if (dist < minDist) {
				minDist = dist;
			}
		}

		if (minDist <= 0) {
			return 0;
		}

		amount = min(amount, minDist * clampAmount);
	}

	seed: s32 = 243703;
	for i: 0..faceCount-1 {
		defer seed = 1664525 * seed + 1013904223;
		slot_s := seed % faceCount;
		faces[i], faces[slot_s] =faces[slot_s], faces[i];
        // btSwap(faces[i], faces[seed % faceCount]);
	}

	for i: 0..faceCount-1 {
		if (!shiftFace(this, faces[i], amount, stack))
		{
			return -amount;
		}
	}

	return amount;
}

shiftFace :: (using this: *btConvexHullInternal, face: *Face, amount: float, stack: [..]*Vertex) -> bool
{
	origShift := getBtNormal(this, face) * -amount;
	if (scaling.component[0] != 0){
		origShift.component[0] /= scaling.component[0];
	}
	if (scaling.component[1] != 0){
		origShift.component[1] /= scaling.component[1];
	}
	if (scaling.component[2] != 0){
		origShift.component[2] /= scaling.component[2];
	}
	shift:=Point32.{cast(s32)origShift.component[medAxis], cast(s32)origShift.component[maxAxis], cast(s32)origShift.component[minAxis], -1};
	if (isZero(shift))	{
		return true;
	}
	normal: Point64 = getNormal(face);
	origDot: s64 = dot(face.origin, normal, s64);
	shiftedOrigin: Point32 = face.origin + shift;
	shiftedDot: s64 = dot(shiftedOrigin,normal, s64);
	btAssert(shiftedDot <= origDot);
	if (shiftedDot >= origDot)	{
		return false;
	}

	intersection: *Edge = null;

	startEdge: *Edge = face.nearbyVertex.edges;
	optDot: Rational128  = dot(face.nearbyVertex, normal);
	cmp: s32 = compare(optDot, shiftedDot);
	if (cmp >= 0) {
		e: *Edge = startEdge;
		while true
		{
			defer if (e == startEdge) break;
			comp: Rational128  = dot(e.target, normal);
			btAssert(compare(comp, origDot) <= 0);
			if (compare(comp, optDot) < 0)
			{
				c: s32 = compare(comp, shiftedDot);
				optDot = comp;
				e = e.reverse;
				startEdge = e;
				if (c < 0)
				{
					intersection = e;
					break;
				}
				cmp = c;
			}
			e = e.prev;
		};

		if (!intersection)	{
			return false;
		}
	} else {
		e := startEdge;
		while true
		{
			defer if (e == startEdge) break;
			component := dot(e.target, normal);
			btAssert(compare(component, origDot) <= 0);
			if compare(component, optDot) > 0 {
				cmp = compare(component, shiftedDot);
				if cmp >= 0 {
					intersection = e;
					break;
				}
				optDot = component;
				e = e.reverse;
				startEdge = e;
			}
			e = e.prev;
		};

		if (!intersection) {
			return true;
		}
	}


	if (cmp == 0) {
		e := intersection.reverse.next;
		while compare(dot(e.target,normal), shiftedDot) <= 0
		{
			e = e.next;
			if (e == intersection.reverse) {
				return true;
			}
		}
	}

	firstIntersection: *Edge;
	faceEdge: *Edge;
	firstFaceEdge: *Edge;

	while (true)
	{
		if (cmp == 0)
		{
			e := intersection.reverse.next;
			startEdge = e;
			while (true)
			{
				if compare(dot(e.target, normal), shiftedDot) >= 0 {
					break;
				}
				intersection = e.reverse;
				e = e.next;
				if (e == startEdge) {
					return true;
				}
			}
		}


		if (!firstIntersection) {
			firstIntersection = intersection;
		} else if (intersection == firstIntersection) {
			break;
		}

		prevCmp: s32 = cmp;
		prevIntersection:= intersection;
		prevFaceEdge:= faceEdge;

		e:= intersection.reverse;
		while (true)
		{
			e = e.reverse.prev;
			btAssert(e != intersection.reverse);
			cmp = compare(dot(e.target, normal),shiftedDot);
			if (cmp >= 0) {
				intersection = e;
				break;
			}
		} if (cmp > 0) {
			removed :*Vertex = intersection.target;
			e = intersection.reverse;
			if (e.prev == e) {
				removed.edges = null;
			} else {
				removed.edges = e.prev;
				link(e.prev, e.next);
				link(e, e);
			}

			n0 :Point64= getNormal(intersection.face);
			n1 :Point64= getNormal(intersection.reverse.face);
			m00 :s64= dot(face.dir0, n0,s64);
			m01 :s64= dot(face.dir1, n0,s64);
			m10 :s64= dot(face.dir0, n1,s64);
			m11 :s64= dot(face.dir1, n1,s64);
			r0 :s64= dot(intersection.face.origin - shiftedOrigin, n0,s64);
			r1 :s64= dot(intersection.reverse.face.origin - shiftedOrigin, n1,s64);
            // Int128 det = Int128::mul(m00, m11) - Int128::mul(m01, m10);
			det: Int128 = mul(m00, m11) - mul(m01, m10);
			btAssert(getSign(det) != 0);
			v := New(Vertex);
			v.point.index = -1;
			v.copy = -1;
			v.point128 = PointR128.{mul(face.dir0.x * r0, m11) - mul(face.dir0.x * r1, m01) + mul(face.dir1.x * r1, m00) - mul(face.dir1.x * r0, m10) + det * shiftedOrigin.x,
									mul(face.dir0.y * r0, m11) - mul(face.dir0.y * r1, m01) + mul(face.dir1.y * r1, m00) - mul(face.dir1.y * r0, m10) + det * shiftedOrigin.y,
									mul(face.dir0.z * r0, m11) - mul(face.dir0.z * r1, m01) + mul(face.dir1.z * r1, m00) - mul(face.dir1.z * r0, m10) + det * shiftedOrigin.z,
									det};
			v.point.x = cast(s32)xvalue(v.point128);
			v.point.y = cast(s32)yvalue(v.point128);
			v.point.z = cast(s32)zvalue(v.point128);
			intersection.target = v;
			v.edges = e;

			array_add(*stack, v);
			array_add(*stack, removed);
			array_add(*stack, null);
		}

		if (cmp || prevCmp || (prevIntersection.reverse.next.target != intersection.target))
		{
			faceEdge = newEdgePair(this, prevIntersection.target, intersection.target);
			if (prevCmp == 0) {
				link(faceEdge, prevIntersection.reverse.next);
			}
			if ((prevCmp == 0) || prevFaceEdge) {
				link(prevIntersection.reverse, faceEdge);
			}
			if (cmp == 0) {
				link(intersection.reverse.prev, faceEdge.reverse);
			}
			link(faceEdge.reverse, intersection.reverse);
		} else {
			faceEdge = prevIntersection.reverse.next;
		}

		if prevFaceEdge {
			if prevCmp > 0 {
				link(faceEdge, prevFaceEdge.reverse);
			} else if faceEdge != prevFaceEdge.reverse {
				array_add(*stack, prevFaceEdge.target);
				while faceEdge.next != prevFaceEdge.reverse {
					removed := faceEdge.next.target;
					removeEdgePair(this, faceEdge.next);
					array_add(*stack, removed);
				}
				array_add(*stack, null);
			}
		}
		faceEdge.face = face;
		faceEdge.reverse.face = intersection.face;

		if (!firstFaceEdge) {
			firstFaceEdge = faceEdge;
		}
	}

	if (cmp > 0) {
		firstFaceEdge.reverse.target = faceEdge.target;
		link(firstIntersection.reverse,firstFaceEdge);
		link(firstFaceEdge,faceEdge.reverse);
	} else if (firstFaceEdge != faceEdge.reverse) {
		array_add(*stack, faceEdge.target);
		while (firstFaceEdge.next != faceEdge.reverse) {
			removed := firstFaceEdge.next.target;
			removeEdgePair(this, firstFaceEdge.next);
			array_add(*stack, removed);
		}
		array_add(*stack, null);
	}

	btAssert(stack.count > 0);
	vertexList = stack[0];

	pos:s32 = 0;
	while (pos < stack.count)
	{
		end :s32= cast(s32)stack.count;
		while (pos < end)
		{
			kept := stack[pos];
			pos += 1;
			deeper := false;
			removed: *Vertex = (stack[pos]);
			while ( removed != null) {
			    defer {
			      pos += 1;
			      (removed = stack[pos]);
			    }
				receiveNearbyFaces(kept, removed);
				while (removed.edges) {
					if (!deeper)
					{
						deeper = true;
						array_add(*stack, kept);
					}
					array_add(*stack, removed.edges.target);
					removeEdgePair(this, removed.edges);
				}
			}
			if (deeper)
			{
				array_add(*stack, null);
			}
		}
	}

	array_reset(*stack);
	face.origin = shiftedOrigin;

	return true;
}

getVertexCopy :: (vertex: *Vertex, vertices: *[..]*Vertex) -> s32 {
	index := vertex.copy;
	if (index < 0)
	{
		index = xx vertices.count;
		vertex.copy = index;
		array_add(vertices, vertex);
	}
	return index;
}


compute_convex_hull :: (using c: *btConvexHullComputer, cloud_of_points: []Vector3, shrink_Factor: float, shrinkClamp: float) -> float {
	count := cloud_of_points.count;
	if (count <= 0)
	{
		array_reset(*vertices_list);
		array_reset(*edges_list);
		array_reset(*faces_list);
		return 0;
	}
	hull: btConvexHullInternal;
	compute(*hull, cloud_of_points);

	shift: float = 0;
    if (shrink_Factor > 0)
    {
      shift = shrink(*hull, shrink_Factor, shrinkClamp);
      if shift < 0 {
       array_reset(*vertices_list);
       array_reset(*edges_list);
       array_reset(*faces_list);
       return shift;
     }
    }

	array_reset(*original_vertex_index);
	array_reset(*vertices_list);
	array_reset(*edges_list);
	array_reset(*faces_list);

	oldVertices: [..]*Vertex;
	getVertexCopy(hull.vertexList, *oldVertices);
	copied: s32 = 0;
	while (copied < oldVertices.count)
	{
		v := oldVertices[copied];
		array_add(*vertices_list, getCoordinates(*hull,v));
		array_add(*original_vertex_index, v.point.index);
		firstEdge := v.edges;
		if firstEdge {
			firstCopy: s32 = -1;
			prevCopy: s32 = -1;
			e := firstEdge;
			while true {
				defer if (e == firstEdge) break;
				if (e.copy < 0)
				{
					s: s32= xx edges_list.count;
					c := array_add(*edges_list);
					r := array_add(*edges_list);
                    // c := *edges[s];
                    // r := *edges[s + 1];
					e.copy = s;
					e.reverse.copy = s + 1;
					c.reverse = 1;
					r.reverse = -1;
					c.targetVertex = getVertexCopy(e.target, *oldVertices);
					r.targetVertex = copied;
				}
				if (prevCopy >= 0)
				{
					edges_list[e.copy].next = prevCopy - e.copy;
				}
				else
				{
					firstCopy = e.copy;
				}
				prevCopy = e.copy;
				e = e.next;
			} ;
			edges_list[firstCopy].next = prevCopy - firstCopy;
		}
		copied+=1;
	}

	for i: 0..copied - 1 {
		v := oldVertices[i];
		firstEdge := v.edges;
		if (firstEdge)
		{
			e := firstEdge;
			while true {
				defer if (e == firstEdge) break;
				if (e.copy >= 0) {
					array_add(*faces_list, e.copy);
					f := e;
					while true {
						defer if (f == e) break;
						f.copy = -1;
						f = f.reverse.prev;
					};
				}
				e = e.next;
			};
		}
	}

	return shift;
}