// Emacs: -*- C++ -*-

//
//	Copyright 1993, Center for Computer Vision and Visualization,
//	University of Florida.  All rights reserved.
//

static char copyright[] = "Copyright 1993, Center for Computer Vision and Visualization,\nUniversity of Florida.  All rights reserved.\n";

static char rcsid[] = "$Id: IntPointSet.c,v 1.30 1994/01/31 16:37:35 thoth Exp $";

// 
// $Log: IntPointSet.c,v $
// Revision 1.30  1994/01/31  16:37:35  thoth
// We now document internal errors when they occur.
// We can now transpose boxy 2D pointsets.
//
// Revision 1.29  1993/11/17  18:03:47  thoth
// IPSIter is now PSIter<IntPoint>
// Pointcmp is now pointcmp
// extensivep() is now extensive()
//
// Revision 1.28  1993/10/20  15:50:54  thoth
// New naming scheme.  Closures are named Closure.
//
// Revision 1.27  1993/09/15  12:36:18  thoth
// Pointset behavior controls have moved to IA::
//
// Revision 1.26  93/08/08  13:17:33  thoth
// Parametrization of the pointsets on the point type using
// template specialization.
// 
// Revision 1.25  93/07/15  11:45:17  thoth
// new lookup_and_set allows other classes (such as BasePSIter)
// to cooperate in refcounting the BasePSs.
// 
// Revision 1.24  93/05/28  15:05:29  thoth
// removed dependency on with-ia-pointset*
// 
// Revision 1.23  93/05/26  16:57:46  thoth
// Copyright Notices
// 
// Revision 1.22  93/05/19  14:44:37  thoth
// LatticeRelation is now in IA::
// 
// Revision 1.21  93/04/29  12:50:41  thoth
// Mondo PointSets are now templatized.
// 
// Revision 1.20  93/04/19  23:27:40  thoth
// Convert BaseIPS to IA_BasePS<IA_IntPoint>
// 
// Revision 1.19  93/04/17  15:19:51  thoth
// IA_IntPoint upgrades (dim, inf, sup)
// 
// Revision 1.18  93/04/16  15:50:18  thoth
// use IA_IntPoint instead of IntPoint.
// 
// Revision 1.17  93/02/12  16:23:41  thoth
// more IA_ prefixes.
// Removed #if 0 cruft.
// comp undefined bug is now gone from CFront
// 
// Revision 1.16  93/01/20  11:36:18  thoth
// new PredicatePS class to allow for families of comprehensive pointsets.
// 
// Revision 1.15  92/12/16  14:48:09  thoth
// conversion to IA_ mostly complete
// 
// Revision 1.14  92/12/15  21:16:21  thoth
// added boxy predicate method.
// 
// Revision 1.13  92/11/11  12:11:46  thoth
// added another restriction utility.
// 
// Revision 1.12  92/11/04  11:58:51  thoth
// New union_with_structure function for image extension.
// 
// Revision 1.11  92/10/25  20:53:24  thoth
// changed argument passing on set restriction routines.
// 
// Revision 1.10  92/10/22  15:11:13  thoth
// lazy comparison handling (untested)
// pointset restriction utilities.
// 
// Revision 1.9  92/10/07  23:15:31  thoth
// trim gratuitous includes from public headers.
// boxy boxy intersect should now be faster.
// 
// Revision 1.8  92/09/30  10:30:52  thoth
// Lazies are now templatized.
// White and Black holes are now semi-public classes, unfortunately.
// reference counts are protected from IntPointSet, unfortunately.
// type() system now based on address of static data member.
// 
// Revision 1.7  92/09/20  17:38:44  thoth
// Minkowski addition of point sets now supported.
// 
// Revision 1.6  92/09/14  21:56:27  thoth
// the IntPointSet is now responsible for deleting hashed pointers.
// dimension is now a parameter to the IntPoint array constructor.
// added special code for ext/comp difference.
// 
// Revision 1.5  92/09/07  13:21:31  thoth
// Black and White holes can be extensive
// |&^/ operators are now friends.
// intersection has untested code to handle comp&ext
// 
// Revision 1.4  92/09/03  13:34:27  thoth
// hash_self now deals with unsigneds
// I unnamed unused formal parameters
// lazy point sets now exist
// (but operators & and / need adjustment)
// 
// Revision 1.3  92/08/26  13:43:37  thoth
// meaningful exception types
// eliminated dummy parameter to constructor (was workaround for g++)
// aborts are clearly labeled and less common
// 
// Revision 1.2  92/08/25  15:28:48  thoth
// new offsetted virtual functions
// new converter from IntPoint to IntPointSet.
// CFront supports agregates in ternary, wheee!
// set difference is operator/, not -
// operator+ is set translation, using the offsetted() function of BaseIPS
// lattice operations on sets defined in terms of compare() friend
// 
// Revision 1.1  92/08/23  13:30:13  thoth
// Initial revision
// 
// 

#include <stdlib.h>
#include "IntPointSet.h"
#include "HashTable.h"
#include "PSIter.h"

#ifdef IA_PRII
#pragma implementation
#endif

#include "pointset_errors.h"
#include "BasePS.h"
#include "BoxyIPS.h"
#include "YoderIPS.h"
#include "LazyPS.h"
#include "MondoPS.h"
#include "ClosurePS.h"

IA_HashTable IA_Set<IA_IntPoint>::hash_table;

void IA_Set<IA_IntPoint>::lookup_and_set(IA_BasePS<IA_IntPoint> *ips)
{
    if (hash_table.reference(ips, &this->ps)) {
	ips->enref();
	if (ips->unref()<=0)
	    delete ips;
    } else {
	// ips->refcount=0;
    }

    ps->enref();
}

void IA_Set<IA_IntPoint>::disassociate_ps()
{
    if ( ps->unref() <= 0 ) {
	hash_table.remove(ps);
	delete ps;
    }
    ps=0;
}

IA_Set<IA_IntPoint>::IA_Set()
{
    this->lookup_and_set(new IA_BlackHole_<IA_IntPoint>(0));
}

IA_Set<IA_IntPoint>::IA_Set(IA_PointSetType t, unsigned dim)
{
    if (t == IA_WhiteHole_<IA_IntPoint>::s_type())
	this->lookup_and_set(new IA_WhiteHole_<IA_IntPoint>(dim));
    else if (t == IA_BlackHole_<IA_IntPoint>::s_type())
	this->lookup_and_set(new IA_BlackHole_<IA_IntPoint>(dim));
    else {
	ia_throw(IA::TYPE_MISMATCH, __FILE__, __LINE__);
	this->lookup_and_set(new IA_BlackHole_<IA_IntPoint>(0));
    }
}

IA_Set<IA_IntPoint>::IA_Set(const IA_IntPoint& p)
{
    this->lookup_and_set(new IA_BoxyIPS(p, p));
}

IA_Set<IA_IntPoint>::IA_Set(const IA_IntPoint& p0, const IA_IntPoint& p1)
{
    this->lookup_and_set(new IA_BoxyIPS(p0, p1));
}

IA_Set<IA_IntPoint>::IA_Set(unsigned dimen, int(*f)(IA_IntPoint)) {
    this->lookup_and_set(new IA_PredicatePS_V<IA_IntPoint>(dimen, f));
}

IA_Set<IA_IntPoint>::IA_Set(unsigned dimen, int(*f)(const IA_IntPoint&)) {
    this->lookup_and_set(new IA_PredicatePS_CR<IA_IntPoint>(dimen, f));
}

IA_Set<IA_IntPoint>::IA_Set(const IA_ClosurePS<IA_IntPoint>&ps)
{
    this->lookup_and_set(ps.clone_self());
}

IA_Set<IA_IntPoint>::IA_Set(const IA_BoxyIPS& bps)
{
    this->lookup_and_set(new IA_BoxyIPS(bps));
}

IA_Set<IA_IntPoint>::IA_Set(const IA_YoderIPS& yps)
{
    if (yps.isboxy()){
	if (yps.nslices<1) {
	    // We've got a black hole which is, of course, boxy
	    this->lookup_and_set(new  IA_BlackHole_<IA_IntPoint>(yps.dimen));
	} else {
	    this->lookup_and_set(new IA_BoxyIPS(yps.min(), yps.max()));
	}
    } else this->lookup_and_set(new IA_YoderIPS(yps));
}


extern "C"
/*static*/ int ipcmp(const void *a, const void *b)
{
  return pointcmp(*(IA_IntPoint*)a,*(IA_IntPoint*)b);
}

IA_Set<IA_IntPoint>::IA_Set(unsigned dimen, const IA_IntPoint *aip, unsigned len)
{
    if (len<1) {
	this->lookup_and_set(new IA_BlackHole_<IA_IntPoint>(dimen));
	return;
    }

    for (unsigned i=0; i<len; i++) {
	if (aip[i].dim() != dimen) {
	    ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	    this->lookup_and_set(new IA_BlackHole_<IA_IntPoint>(0));
	    return;
	}
    }

    if (dimen<1) {
	if (len>0)
	    this->lookup_and_set(new IA_WhiteHole_<IA_IntPoint>(0));
	else
	    this->lookup_and_set(new IA_BlackHole_<IA_IntPoint>(0));
	return;
    }

    for (i=0; i+1<len; i++) {
	if (pointcmp(aip[i],aip[i+1]) ==1)
	    break;
    }

    IA_IntPoint	*newaip;
    if (i+1<len) {
	newaip = new IA_IntPoint[len];
	for (i=0; i<len; i++)
	    newaip[i] = aip[i];
	SortIntPointArray(newaip, len);
#if 0
	cout << "sorting array\n" << flush;
	for (i=0; i<len; i++)
	    cout << newaip[i] << ' ';
	cout << '\n';
#endif
    } else
	newaip = (IA_IntPoint*)aip;

    IA_YoderIPS *temp = new IA_YoderIPS(newaip, len);
    if (temp->isboxy()) {
	if (temp->nslices<1)
	    this->lookup_and_set(new IA_BlackHole_<IA_IntPoint>(temp->dimen));
	else
	    this->lookup_and_set(new IA_BoxyIPS(temp->min(), temp->max()));
	delete temp;
    } else
	this->lookup_and_set(temp);

    if (newaip!=aip)
	delete [] newaip;
}

IA_Set<IA_IntPoint>::~IA_Set()
{
    if (ps==0) {
	cerr << "Internal Error: pointset destructed twice!\n";
	return;
    }
    this->disassociate_ps();
    ps = 0;
}


IA_Set<IA_IntPoint>& IA_Set<IA_IntPoint>::operator=(const IA_Set<IA_IntPoint> &ips)
{ 
    IA_BasePS<IA_IntPoint>	*temp;
    temp = ips.ps;
    temp->enref();
    this->disassociate_ps();
    ps = temp;

    return *this;
}

int IA_Set<IA_IntPoint>::boxy() const
{
    return ps->type() == IA_BoxyIPS::s_type();
}

IA_Set<IA_IntPoint> operator|(const IA_Set<IA_IntPoint>& lhs,const IA_Set<IA_IntPoint> &rhs)
{
    if (lhs.dim() != rhs.dim()) {
	// need better message
	// cout << lhs.dim() << "!=" << rhs.dim() << '\n' << flush;
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return IA_Set<IA_IntPoint>();
    }

    if (lhs.ps->type() == IA_BlackHole_<IA_IntPoint>::s_type())
	return rhs;
    if (lhs.ps->type() == IA_WhiteHole_<IA_IntPoint>::s_type())
	return lhs;

    if (rhs.ps->type() == IA_BlackHole_<IA_IntPoint>::s_type())
	return lhs;
    if (rhs.ps->type() == IA_WhiteHole_<IA_IntPoint>::s_type())
	return rhs;

    if (lhs.extensive() && rhs.extensive()) {
#if 0
	YoderIPS	*left, *right;
	switch (lhs.ps->type()) {
	case BOXY:
	    left = new YoderIPS(*(BoxyIPS*)lhs.ps);
	    break;
	case YODER:
	    left = (YoderIPS*)lhs.ps;
	    break;
	default:
	    // the other extensive pointsets were taken care of before.
	    IA::internal_error(__FILE__,__LINE__);
	    break;
	}

	switch (rhs.ps->type()) {
	case BOXY:
	    right = new YoderIPS(*(BoxyIPS*)rhs.ps);
	    break;
	case YODER:
	    right = (YoderIPS*)rhs.ps;
	    break;
	default:
	    // the other extensive pointsets were taken care of before.
	    IA::internal_error(__FILE__,__LINE__);
	    break;
	}

	YoderIPS result( *left | *right );

	if (left != lhs.ps)
	    delete left;
	if (right != rhs.ps)
	    delete right;
#else
	IA_YoderIPS result
	  ( ( (lhs.ps->type()==IA_BoxyIPS::s_type()) ?
	     IA_YoderIPS(*(IA_BoxyIPS*)lhs.ps) :
	     *(IA_YoderIPS*)lhs.ps )
	   |
	   ( (rhs.ps->type()==IA_BoxyIPS::s_type()) ?
	    IA_YoderIPS(*(IA_BoxyIPS*)rhs.ps) :
	    *(IA_YoderIPS*)rhs.ps )
	   );
#endif

	// this constructor checks for emptiness and boxiness
	return IA_Set<IA_IntPoint>(result);
    } else {
	return new IA_LazyUnionPS<IA_Set<IA_IntPoint>,IA_IntPoint>(lhs, rhs);
    }
}

//

// intersect an extensive point set with a comprehensive one
// precondition:
// ext.extensive() ==1
// ext.dim() == comp.dim()
IA_Set<IA_IntPoint> ext_comp_intersect(const IA_Set<IA_IntPoint> &ext, const IA_Set<IA_IntPoint> &comp)
{
    IA_PSIter<IA_IntPoint>	iter(ext);
    IA_IntPoint	ip;

    IA_IntPoint	*arr=new IA_IntPoint[ext.card()];
    unsigned	count=0;

    while (iter(ip)) {
	if (comp.contains(ip)) {
	    arr[count++] = ip;
	}
    }

    return IA_Set<IA_IntPoint>(ext.dim(), arr,count);
}

IA_Set<IA_IntPoint> operator&(const IA_Set<IA_IntPoint>& lhs, const IA_Set<IA_IntPoint> &rhs)
{
    if (lhs.dim() != rhs.dim()) {
	// need better message
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return IA_Set<IA_IntPoint>();
    }

    if (lhs.ps->type()==IA_BlackHole_<IA_IntPoint>::s_type())
	return lhs;
    if (lhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type())
	return rhs;

    if (rhs.ps->type()==IA_BlackHole_<IA_IntPoint>::s_type())
	return rhs;
    if (rhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type())
	return lhs;

    int lext, rext;
    lext = lhs.extensive();
    rext = rhs.extensive();
    if (lext && rext) {
	if (lhs.ps->type()==IA_BoxyIPS::s_type() &&
	    rhs.ps->type()==IA_BoxyIPS::s_type() ) {
	    IA_IntPoint	bottom(sup(lhs.ps->inf(), rhs.ps->inf()));
	    IA_IntPoint	top(inf (lhs.ps->sup(), rhs.ps->sup()));
	    if (bottom<=top)
		return IA_Set<IA_IntPoint>(bottom, top);
	    else
		return new IA_BlackHole_<IA_IntPoint>(lhs.dim());
	} else {
	    IA_YoderIPS result
		( ( (lhs.ps->type()==IA_BoxyIPS::s_type()) ?
		    IA_YoderIPS(*(IA_BoxyIPS*)lhs.ps) :
		    *(IA_YoderIPS*)lhs.ps )
		  &
		  ( (rhs.ps->type()==IA_BoxyIPS::s_type()) ?
		    IA_YoderIPS(*(IA_BoxyIPS*)rhs.ps) :
		    *(IA_YoderIPS*)rhs.ps )
		  );
	    return IA_Set<IA_IntPoint>(result);
	}
    } else if (lext) {
	return ext_comp_intersect(lhs, rhs);
    } else if (rext) {
	return ext_comp_intersect(rhs, lhs);
    } else {
	return new IA_LazyIntersectionPS<IA_Set<IA_IntPoint>,IA_IntPoint>(lhs,rhs);
    }
}


IA_Set<IA_IntPoint> intersect_with_structure
(const IA_Set<IA_IntPoint> &lhs, const IA_Set<IA_IntPoint> &rhs, IA_SetStructure *result)
{
    if (lhs.dim() != rhs.dim()) {
	// need better message
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return IA_Set<IA_IntPoint>();
    }

    if (lhs.ps->type()==IA_BlackHole_<IA_IntPoint>::s_type()) {
	return lhs;
    }
    if (lhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type()) {
	*result = (rhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type())
	    ? IA_SetStructure::BOTH
	    : IA_SetStructure::FIRST_ONLY;
	return rhs;
    }

    if (rhs.ps->type()==IA_BlackHole_<IA_IntPoint>::s_type()) {
	result->add_interval(lhs.card(), IA_SetStructure::FIRST_ONLY);
	return rhs;
    }
    if (rhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type()) {
	result->add_interval(lhs.card(), IA_SetStructure::BOTH);
	return lhs;
    }

    if (!lhs.extensive()) {
	// not rational
	return lhs & rhs;
    }

    if (!rhs.extensive()) {
	// deal with this later
	return lhs & rhs;
    }

    return ( (lhs.ps->type()==IA_BoxyIPS::s_type()) ?
	     IA_YoderIPS(*(IA_BoxyIPS*)lhs.ps) :
	     *(IA_YoderIPS*)lhs.ps ). intersect_with_structure(
		 ( (rhs.ps->type()==IA_BoxyIPS::s_type()) ?
		   IA_YoderIPS(*(IA_BoxyIPS*)rhs.ps) :
		   *(IA_YoderIPS*)rhs.ps ),
		 result);
}

IA_Set<IA_IntPoint> intersect_with_dualstruct
(const IA_Set<IA_IntPoint> &lhs, const IA_Set<IA_IntPoint> &rhs, IA_SetStructure *result)
{
    if (lhs.dim() != rhs.dim()) {
	// need better message
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return IA_Set<IA_IntPoint>();
    }

    if (lhs.ps->type()==IA_BlackHole_<IA_IntPoint>::s_type()) {
	result->add_interval(rhs.card(), IA_SetStructure::SECOND_ONLY);
	return lhs;
    }
    if (lhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type()) {
	*result = (rhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type())
	    ? IA_SetStructure::BOTH
	    : IA_SetStructure::FIRST_ONLY;
	return rhs;
    }

    if (rhs.ps->type()==IA_BlackHole_<IA_IntPoint>::s_type()) {
	result->add_interval(lhs.card(), IA_SetStructure::FIRST_ONLY);
	return rhs;
    }
    if (rhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type()) {
	result->add_interval(lhs.card(), IA_SetStructure::BOTH);
	return lhs;
    }

    if (!lhs.extensive()) {
	// not rational
	return lhs & rhs;
    }

    if (!rhs.extensive()) {
	// deal with this later
	return lhs & rhs;
    }

    return ( (lhs.ps->type()==IA_BoxyIPS::s_type()) ?
	     IA_YoderIPS(*(IA_BoxyIPS*)lhs.ps) :
	     *(IA_YoderIPS*)lhs.ps ). intersect_with_dualstruct(
		 ( (rhs.ps->type()==IA_BoxyIPS::s_type()) ?
		   IA_YoderIPS(*(IA_BoxyIPS*)rhs.ps) :
		   *(IA_YoderIPS*)rhs.ps ),
		 result);
}

IA_Set<IA_IntPoint> union_with_structure
(const IA_Set<IA_IntPoint> &lhs, const IA_Set<IA_IntPoint> &rhs, IA_SetStructure *result)
{
    if (lhs.dim() != rhs.dim()) {
	// need better message
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return IA_Set<IA_IntPoint>();
    }

    if (lhs.ps->type()==IA_BlackHole_<IA_IntPoint>::s_type()) {
	result->add_interval(rhs.card(), IA_SetStructure::SECOND_ONLY);
	return rhs;
    }
    if (lhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type()) {
	*result = (rhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type())
	    ? IA_SetStructure::BOTH
	    : IA_SetStructure::FIRST_ONLY;
	return lhs;
    }

    if (rhs.ps->type()==IA_BlackHole_<IA_IntPoint>::s_type()) {
	result->add_interval(lhs.card(), IA_SetStructure::FIRST_ONLY);
	return lhs;
    }
    if (rhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type()) {
	*result = (lhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type())
	    ? IA_SetStructure::BOTH
	    : IA_SetStructure::SECOND_ONLY;
	return rhs;
    }

    if (!lhs.extensive() || !rhs.extensive()) {
	return lhs | rhs;
    }

    return ( (lhs.ps->type()==IA_BoxyIPS::s_type()) ?
	     IA_YoderIPS(*(IA_BoxyIPS*)lhs.ps) :
	     *(IA_YoderIPS*)lhs.ps ). union_with_structure(
		 ( (rhs.ps->type()==IA_BoxyIPS::s_type()) ?
		   IA_YoderIPS(*(IA_BoxyIPS*)rhs.ps) :
		   *(IA_YoderIPS*)rhs.ps ),
		 result);
}


#if 1
IA_Set<IA_IntPoint> operator^(const IA_Set<IA_IntPoint>& lhs, const IA_Set<IA_IntPoint> &rhs)
{
    if (lhs.dim() != rhs.dim()) {
	// need better message
	// cout << lhs.dim() << "!=" << rhs.dim() << '\n' << flush;
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return IA_Set<IA_IntPoint>();
    }

    if (lhs.ps->type()==IA_BlackHole_<IA_IntPoint>::s_type())
	return rhs;
    if (lhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type())
	return ~rhs;

    if (rhs.ps->type()==IA_BlackHole_<IA_IntPoint>::s_type())
	return lhs;
    if (rhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type())
	return ~lhs;

    if (lhs.extensive() && rhs.extensive()) {
	IA_YoderIPS result
	  ( ( (lhs.ps->type()==IA_BoxyIPS::s_type()) ?
	     IA_YoderIPS(*(IA_BoxyIPS*)lhs.ps) :
	     *(IA_YoderIPS*)lhs.ps )
	   ^
	   ( (rhs.ps->type()==IA_BoxyIPS::s_type()) ?
	    IA_YoderIPS(*(IA_BoxyIPS*)rhs.ps) :
	    *(IA_YoderIPS*)rhs.ps )
	   );

	return IA_Set<IA_IntPoint>(result);
    } else {
#if 1
	return new IA_LazyXORPS<IA_Set<IA_IntPoint>,IA_IntPoint>(lhs, rhs);
    }
#else
    } // this is where the curly brace should NOT be :(
    /* The following code belongs in the above else clause.
       However, if it is placed there, CFront generates erroneous C code.
       Since the code MAY reside outside the scope of the else and still
       be semantically equivalent, and doing so suppresses the error, 
       we will do this until CFront is fixed
       */
	int	kludge;	// if this declaration is removed, CFront flakes out
	IA_BasePS<IA_IntPoint>	*temp =  new IA_LazyXORPS<IntPointSet,IA_IntPoint>(lhs, rhs);
	return temp;
//    } // this is where the curly brace should be
#endif
}
#endif


/*
  asymmetric set difference
  */

// subtract a comprehensive point set from an extensive one
// precondition:
// ext.extensive() ==1
// ext.dim() == comp.dim()
IA_Set<IA_IntPoint> ext_comp_diff(const IA_Set<IA_IntPoint> &ext, const IA_Set<IA_IntPoint> &comp)
{
    IA_PSIter<IA_IntPoint>	iter(ext);
    IA_IntPoint	ip;

    IA_IntPoint	*arr=new IA_IntPoint[ext.card()];
    unsigned	count=0;

    while (iter(ip)) {
	if (!comp.contains(ip)) {
	    arr[count++] = ip;
	}
    }

    return IA_Set<IA_IntPoint>(ext.dim(), arr,count);
}

IA_Set<IA_IntPoint> operator/(const IA_Set<IA_IntPoint>& lhs, const IA_Set<IA_IntPoint> &rhs)
{
    if (lhs.dim() != rhs.dim()) {
	// need better message
	// cout << lhs.dim() << "!=" << rhs.dim() << '\n' << flush;
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return IA_Set<IA_IntPoint>();
    }

    if (lhs.ps->type()==IA_BlackHole_<IA_IntPoint>::s_type())
	return lhs;
    if (lhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type())
	return ~rhs;

    if (rhs.ps->type()==IA_BlackHole_<IA_IntPoint>::s_type())
	return lhs;
    if (rhs.ps->type()==IA_WhiteHole_<IA_IntPoint>::s_type())
	return IA_BlackHole(lhs.dim());

    if (lhs.extensive() && rhs.extensive()) {
	IA_YoderIPS result
	  ( ( (lhs.ps->type()==IA_BoxyIPS::s_type()) ?
	     IA_YoderIPS(*(IA_BoxyIPS*)lhs.ps) :
	     *(IA_YoderIPS*)lhs.ps )
	   /
	   ( (rhs.ps->type()==IA_BoxyIPS::s_type()) ?
	    IA_YoderIPS(*(IA_BoxyIPS*)rhs.ps) :
	    *(IA_YoderIPS*)rhs.ps )
	   );

	return IA_Set<IA_IntPoint>(result);
    } else if (lhs.extensive()) {
	return ext_comp_diff(lhs, rhs);
    } else {
	return new IA_LazyMinusPS<IA_Set<IA_IntPoint>,IA_IntPoint>(lhs,rhs);
    }
}

IA_Set<IA_IntPoint> IA_Set<IA_IntPoint>::operator~() const
{
    if (ps->type() == IA_BlackHole_<IA_IntPoint>::s_type())
	return IA_WhiteHole(dim());
    else if (ps->type() == IA_WhiteHole_<IA_IntPoint>::s_type())
	return IA_BlackHole(dim());
    else 
	return new IA_LazyNotPS<IA_Set<IA_IntPoint>, IA_IntPoint>(*this);
}


IA_Set<IA_IntPoint> operator+(const IA_Set<IA_IntPoint>& ips, const IA_IntPoint &ip)
{
    if (ips.extensive())
	return ips.ps->offsetted(ip);
    else
	return new IA_LazyOffsettedPS<IA_Set<IA_IntPoint>, IA_IntPoint>(ips, ip);
}

static IA_Set<IA_IntPoint> Mink_aux(const IA_Set<IA_IntPoint>& small, const IA_Set<IA_IntPoint>& large)
{
    IA_PSIter<IA_IntPoint>	iter(small);
    IA_IntPoint	ip;
    IA_Set<IA_IntPoint>	rval(IA_BlackHole(small.dim()));

    while (iter(ip)) {
	rval |= large+ip;
    }

    return rval;
}

IA_Set<IA_IntPoint> operator+(const IA_Set<IA_IntPoint>& lhs, const IA_Set<IA_IntPoint>& rhs)
{
    if (lhs.dim() != rhs.dim()) {
	// cout << lhs.dim() << "!=" << rhs.dim() << '\n' << flush;
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return IA_Set<IA_IntPoint>();
    }

    if (lhs.extensive()) {
	if (rhs.extensive()) {
	    return (lhs.card() < rhs.card()) ?
		Mink_aux(lhs, rhs) : Mink_aux(rhs,lhs);

	} else {
	    return new IA_LazyMinkowskiPS<IA_PSIter<IA_IntPoint>,IA_Set<IA_IntPoint>, IA_IntPoint>(lhs,rhs);
	}
    } else {
	if (rhs.extensive()) {
	    return new IA_LazyMinkowskiPS<IA_PSIter<IA_IntPoint>,IA_Set<IA_IntPoint>, IA_IntPoint>(rhs,lhs);
	} else {
	    ia_throw(IA::PSET_REQUIRE_EXTENSIVE, __FILE__, __LINE__);
	}
    }
}


// this doesn't handle lazys
IA::LatticeRelation
compare(const IA_Set<IA_IntPoint> &lhs, const IA_Set<IA_IntPoint> &rhs)
{
    if (lhs.dim() != rhs.dim()) {
	// need better message
	// cout << lhs.dim() << "!=" << rhs.dim() << '\n' << flush;
	ia_throw(IA::PSET_DIMENSION_MISMATCH, __FILE__, __LINE__);
	return IA::UNKNOWN;
    }

    if (lhs.ps==rhs.ps)
	return IA::EQUAL;

    if (!lhs.extensive()) {
	if (!rhs.extensive())
	    return IA::UNKNOWN;
	// else, we can check if lhs>=rhs
	IA_PSIter<IA_IntPoint>	iter(rhs);
	IA_IntPoint	ip;
	while (iter(ip)) {
	    if (!lhs.contains(ip))
		return IA::UNKNOWN;
	}
	return IA::MP_SUPERSET;
    } else if (!rhs.extensive()) {
	// else, we can check if lhs<=rhs
	IA_PSIter<IA_IntPoint>	iter(lhs);
	IA_IntPoint	ip;
	while (iter(ip)) {
	    if (!rhs.contains(ip))
		return IA::UNKNOWN;
	}
	return IA::MP_SUBSET;
    }

    // invariant: lhs.extensive() && rhs.extensive()

    IA_PointSetType ltype=lhs.ps->type();
    IA_PointSetType rtype=rhs.ps->type();
    if (ltype==IA_WhiteHole_<IA_IntPoint>::s_type()) {
	if (rtype==IA_WhiteHole_<IA_IntPoint>::s_type())
	    return IA::EQUAL;
	else if (rtype==IA_BlackHole_<IA_IntPoint>::s_type() ||
		 rtype==IA_YoderIPS::s_type() ||
		 rtype==IA_BoxyIPS::s_type())
	    return IA::P_SUPERSET;
	else {
	    IA::internal_error(__FILE__,__LINE__);
	}
    } else if (ltype==IA_BlackHole_<IA_IntPoint>::s_type()) {
	if (rtype==IA_BlackHole_<IA_IntPoint>::s_type())
	    return IA::EQUAL;
	else if (rtype==IA_WhiteHole_<IA_IntPoint>::s_type() ||
		 rtype==IA_YoderIPS::s_type() ||
		 rtype==IA_BoxyIPS::s_type())
	    return IA::P_SUBSET;
	else {
	    IA::internal_error(__FILE__,__LINE__);
	}
    } else if (ltype==IA_BoxyIPS::s_type() ||
	       ltype==IA_YoderIPS::s_type()) {
	if (rtype==IA_WhiteHole_<IA_IntPoint>::s_type())
	    return IA::P_SUBSET;
	else if (rtype==IA_BlackHole_<IA_IntPoint>::s_type())
	    return IA::P_SUPERSET;
	else if (rtype==IA_BoxyIPS::s_type() ||
		 rtype==IA_YoderIPS::s_type())
	    return compare( (lhs.ps->type()==IA_BoxyIPS::s_type())?
			    IA_YoderIPS(*(IA_BoxyIPS*)lhs.ps) : *(IA_YoderIPS*)lhs.ps ,
			    (rhs.ps->type()==IA_BoxyIPS::s_type())?
			    IA_YoderIPS(*(IA_BoxyIPS*)rhs.ps) : *(IA_YoderIPS*)rhs.ps);
	else {
	    IA::internal_error(__FILE__,__LINE__);
	}
    }
    IA::internal_error(__FILE__,__LINE__);
}


int operator<(const IA_Set<IA_IntPoint> &lhs, const IA_Set<IA_IntPoint> &rhs)
{
    IA::LatticeRelation	result;
    result = compare(lhs,rhs);
    return result==IA::P_SUBSET;
}

int operator<=(const IA_Set<IA_IntPoint> &lhs, const IA_Set<IA_IntPoint> &rhs)
{
    IA::LatticeRelation	result;
    result = compare(lhs,rhs);
    return result==IA::P_SUBSET
	|| result==IA::MP_SUBSET
	    || result==IA::EQUAL;
}

int operator>=(const IA_Set<IA_IntPoint> &lhs, const IA_Set<IA_IntPoint> &rhs)
{
    IA::LatticeRelation	result;
    result = compare(lhs,rhs);
    return result==IA::P_SUPERSET
	|| result==IA::MP_SUPERSET
	    || result==IA::EQUAL;
}

int operator>(const IA_Set<IA_IntPoint> &lhs, const IA_Set<IA_IntPoint> &rhs)
{
    IA::LatticeRelation	result;
    result = compare(lhs,rhs);
    return result==IA::P_SUPERSET;
}

//
//
//


IA_Set<IA_IntPoint> IA_BlackHole(unsigned dim)
{
    return IA_Set<IA_IntPoint>(IA_BlackHole_<IA_IntPoint>::s_type(), dim);
}

IA_Set<IA_IntPoint> IA_WhiteHole(unsigned dim)
{
    return IA_Set<IA_IntPoint>(IA_WhiteHole_<IA_IntPoint>::s_type(), dim);
}

//
//

void SortIntPointArray(IA_IntPoint *aip, unsigned len)
{
  // We can use qsort because a swap of the bytes in
  // two IntPoints does not do any damage.
  qsort((char*)aip, len, sizeof(*aip), ipcmp);
}

//
//


IA_Set<IA_IntPoint> transpose(const IA_Set<IA_IntPoint> &ps)
{
    if (!ps.boxy() || ps.dim() != 2) {
	ia_throw(Pointset_BadTranspose_Exception(__FILE__, __LINE__));
	return IA_Set<IA_IntPoint>();
    }

    return IA_IntPointSet(ps.inf().derange(1,0), ps.sup().derange(1,0));
}
