//
//	Copyright 1993, Center for Computer Vision and Visualization,
//	University of Florida.  All rights reserved.
//
//
// $Log:	float-convolutions.desc,v $
// Revision 1.2  93/05/27  11:41:03  thoth
// Copyright Notices
// 
// Revision 1.1  93/05/26  17:16:30  thoth
// Initial revision
// 
// Revision 1.2  93/05/18  21:44:57  thoth
// New more descriptive names for standard convolutions.
// 
// Revision 1.1  93/03/18  11:23:20  thoth
// Initial revision
// 
// float linear_product float float zero=0
// Emacs: -*- C++ -*-

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


//
// $Log:	absorber,v $
// Revision 1.5  93/05/27  11:48:55  thoth
// Copyright Notices
// 
// Revision 1.4  93/04/29  11:20:56  thoth
// Faster(?) extension for VectorDIs.
// 
// Revision 1.3  93/04/17  18:56:44  jnw
// Fixed to match IA_IntPoint membername (dim)
// 
// Revision 1.2  93/04/08  13:21:59  thoth
// internal helper functions are now static.
// 
// Revision 1.1  93/03/18  11:39:42  thoth
// Initial revision
// 

#ifndef zero_extend_scan_float
#define zero_extend_scan_float

static void zero_extend_vec_scan(const IA_SetStructure &ss,
				 const float **src,
				 float **dest,
				 float zero)
{
    for (unsigned i=0; i<ss.nintervals(); i++) {
	IA_ss_interval	temp(ss.retrieve_interval(i));
	if (temp.substructure == IA_SetStructure::FIRST_ONLY) {
	    (*src) += temp.count;
	} else if (temp.substructure == IA_SetStructure::BOTH) {
	    for (unsigned j=0; j<temp.count; j++) {
		*((*dest)++) = *((*src)++);
	    }
	} else if (temp.substructure == IA_SetStructure::SECOND_ONLY) {
	    for (unsigned j=0; j<temp.count; j++)
		*((*dest)++) = zero;
	} else {
	    for (unsigned j=0; j<temp.count; j++)
		zero_extend_vec_scan(temp.substructure, src, dest, zero);
	}
    }
}

static void zero_extend_iter_scan(const IA_SetStructure &ss,
				  IA_DIVIter<float> *srciter,
				  float **dest,
				  float zero)
{
    for (unsigned i=0; i<ss.nintervals(); i++) {
	IA_ss_interval	temp(ss.retrieve_interval(i));
	if (temp.substructure == IA_SetStructure::FIRST_ONLY) {
	    float	blah;
	    for (unsigned j=0; j<temp.count; j++) 
		(*srciter)(blah);
	} else if (temp.substructure == IA_SetStructure::BOTH) {
	    float	blah;
	    for (unsigned j=0; j<temp.count; j++) {
		(*srciter)(blah);
		*((*dest)++) = blah;
	    }
	} else if (temp.substructure == IA_SetStructure::SECOND_ONLY) {
	    for (unsigned j=0; j<temp.count; j++)
		*((*dest)++) = zero;
	} else {
	    for (unsigned j=0; j<temp.count; j++)
		zero_extend_iter_scan(temp.substructure, srciter, dest, zero);
	}
    }
}

static void zero_extend(const IA_FloatDiscreteImage &srcimg,
			const IA_IntPointSet &dest_ps, float *dest_data,
			float zero)
{
    IA_SetStructure	ss;
    intersect_with_dualstruct(srcimg.domain(), dest_ps, &ss);

    if (srcimg.type() == IA_VectorDI<float>::s_type()) {
	float	*srcdata = ((IA_VectorDI<float>*)srcimg.bdip)->vec;
	zero_extend_vec_scan(ss, &srcdata, &dest_data, zero);
    } else {
	IA_DIVIter<float>	iter(srcimg);
	zero_extend_iter_scan(ss, &iter, &dest_data, zero);
    }
}

#endif

static IA_FloatDiscreteImage
linear_product_inv_core(IA_IntPoint src_infimum,
	      IA_IntPoint src_width,
	      const float *src_data, // length is prod(src_width)
	      const IA_FloatDiscreteImage &templ,
	      IA_IntPointSet dest_ps)
{
    const int	dimen = src_width.dim();

    IA_IntPointSet	templ_ps = templ.domain();
    int		templ_sz = templ_ps.card();
    float	*const templ_data = new float[templ_sz];
    int	*const templ_offsets = new int[templ_sz];

    {
	float	*d_scan = templ_data;
	int	*o_scan = templ_offsets;
	IA_DIPIter<float>	iter(templ);
	IA_IntPoint	ip;
	while (iter(ip, *d_scan)) {
	    *o_scan = ip[0];
	    for (unsigned i=1; i<dimen; i++) {
		*o_scan *= src_width[i];
		*o_scan += ip[i];
	    }
	    d_scan++;
	    o_scan++;
	}
    }

    IA_IPSIter	iter(dest_ps);
    IA_IntPoint	ip;
    float *const	dest_data = new float[dest_ps.card()];
    float *	valp = dest_data;
    while (iter(ip)) {
	int	offset= ip[0] - src_infimum[0];
	unsigned i;
	for (i=1; i<dimen; i++) {
	    offset *= src_width[i];
	    offset += ip[i] - src_infimum[i];
	}
	const float *const base = src_data + offset;

	// _IVAL_=(base[templ_offsets[i]])
	// _TVAL_=(templ_data[i])
	// _IRESULT_=(*valp)

	float	a = 0; ;

	for (i=0; i<templ_sz; i++) {
	    a += (base[templ_offsets[i]]) * (templ_data[i]); ;
	}
	{
	    (*valp) = a; ;
	}
	valp++;
    }
    delete[] templ_offsets;
    delete[] templ_data;

    return IA_FloatDiscreteImage(dest_ps, dest_data, dest_ps.card(), 1);
}

static IA_FloatDiscreteImage
linear_product_inv(const IA_FloatDiscreteImage &img,
	 const IA_FloatDiscreteImage &invtempl,
	 IA_IntPointSet dest_ps)
{
    IA_IntPoint	inf_ = dest_ps.inf()+invtempl.domain().inf();
    IA_IntPoint	sup_ = dest_ps.sup()+invtempl.domain().sup();
    IA_IntPointSet	src_ps = IA_IntPointSet(inf_,sup_);

    float	*const src_data = new float[src_ps.card()];

    zero_extend(img, src_ps, src_data, 0);
    IA_FloatDiscreteImage	rval =
	linear_product_inv_core(inf_, sup_-inf_ + 1, src_data, invtempl, dest_ps);

    delete[] src_data;

    return rval;
}

IA_FloatDiscreteImage linear_product(const IA_FloatDiscreteImage &img,
			     const IA_FloatDDTemplate &templ,
			     IA_IntPointSet dest_ps)
{
    if (templ.type() == IA_InvariantDT<IA_FloatDiscreteImage >::s_type()) {
	return linear_product_inv(img, ((IA_InvariantDT<IA_FloatDiscreteImage >*)templ.bdtp)->value, dest_ps);
    } else {
	float *const	dest_data = new float[dest_ps.card()];
	float	*valp = dest_data;
	IA_IPSIter	dest_iter(dest_ps);
	IA_IntPoint	base_ip;
	while (dest_iter(base_ip)) {
	    IA_FloatDiscreteImage	tv = templ(base_ip);

	    IA_DIPIter<float>	templ_iter(tv);
	    IA_IntPoint	templ_ip;
	    float	templ_val;

	    // _IVAL_=img(ip)
	    // _TVAL_=templ_val
	    // _IRESULT_=(*valp)

	    float	a = 0; ;

	    while ( templ_iter(templ_ip, templ_val) ) {
		const IA_IntPoint	ip = templ_ip; //+base_ip;
		if (! img.domain().contains(ip))
		    continue;
		a += img(ip) * templ_val; ;
	    }
	    {
		(*valp) = a; ;
	    }
	    valp++;
	}
	return IA_FloatDiscreteImage(dest_ps, dest_data, dest_ps.card(), 1);
    }
}
// float addmax_product float float  zero=-infinity()
// Emacs: -*- C++ -*-

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


//
// $Log:	absorber,v $
// Revision 1.5  93/05/27  11:48:55  thoth
// Copyright Notices
// 
// Revision 1.4  93/04/29  11:20:56  thoth
// Faster(?) extension for VectorDIs.
// 
// Revision 1.3  93/04/17  18:56:44  jnw
// Fixed to match IA_IntPoint membername (dim)
// 
// Revision 1.2  93/04/08  13:21:59  thoth
// internal helper functions are now static.
// 
// Revision 1.1  93/03/18  11:39:42  thoth
// Initial revision
// 

#ifndef zero_extend_scan_float
#define zero_extend_scan_float

static void zero_extend_vec_scan(const IA_SetStructure &ss,
				 const float **src,
				 float **dest,
				 float zero)
{
    for (unsigned i=0; i<ss.nintervals(); i++) {
	IA_ss_interval	temp(ss.retrieve_interval(i));
	if (temp.substructure == IA_SetStructure::FIRST_ONLY) {
	    (*src) += temp.count;
	} else if (temp.substructure == IA_SetStructure::BOTH) {
	    for (unsigned j=0; j<temp.count; j++) {
		*((*dest)++) = *((*src)++);
	    }
	} else if (temp.substructure == IA_SetStructure::SECOND_ONLY) {
	    for (unsigned j=0; j<temp.count; j++)
		*((*dest)++) = zero;
	} else {
	    for (unsigned j=0; j<temp.count; j++)
		zero_extend_vec_scan(temp.substructure, src, dest, zero);
	}
    }
}

static void zero_extend_iter_scan(const IA_SetStructure &ss,
				  IA_DIVIter<float> *srciter,
				  float **dest,
				  float zero)
{
    for (unsigned i=0; i<ss.nintervals(); i++) {
	IA_ss_interval	temp(ss.retrieve_interval(i));
	if (temp.substructure == IA_SetStructure::FIRST_ONLY) {
	    float	blah;
	    for (unsigned j=0; j<temp.count; j++) 
		(*srciter)(blah);
	} else if (temp.substructure == IA_SetStructure::BOTH) {
	    float	blah;
	    for (unsigned j=0; j<temp.count; j++) {
		(*srciter)(blah);
		*((*dest)++) = blah;
	    }
	} else if (temp.substructure == IA_SetStructure::SECOND_ONLY) {
	    for (unsigned j=0; j<temp.count; j++)
		*((*dest)++) = zero;
	} else {
	    for (unsigned j=0; j<temp.count; j++)
		zero_extend_iter_scan(temp.substructure, srciter, dest, zero);
	}
    }
}

static void zero_extend(const IA_FloatDiscreteImage &srcimg,
			const IA_IntPointSet &dest_ps, float *dest_data,
			float zero)
{
    IA_SetStructure	ss;
    intersect_with_dualstruct(srcimg.domain(), dest_ps, &ss);

    if (srcimg.type() == IA_VectorDI<float>::s_type()) {
	float	*srcdata = ((IA_VectorDI<float>*)srcimg.bdip)->vec;
	zero_extend_vec_scan(ss, &srcdata, &dest_data, zero);
    } else {
	IA_DIVIter<float>	iter(srcimg);
	zero_extend_iter_scan(ss, &iter, &dest_data, zero);
    }
}

#endif

static IA_FloatDiscreteImage
addmax_product_inv_core(IA_IntPoint src_infimum,
	      IA_IntPoint src_width,
	      const float *src_data, // length is prod(src_width)
	      const IA_FloatDiscreteImage &templ,
	      IA_IntPointSet dest_ps)
{
    const int	dimen = src_width.dim();

    IA_IntPointSet	templ_ps = templ.domain();
    int		templ_sz = templ_ps.card();
    float	*const templ_data = new float[templ_sz];
    int	*const templ_offsets = new int[templ_sz];

    {
	float	*d_scan = templ_data;
	int	*o_scan = templ_offsets;
	IA_DIPIter<float>	iter(templ);
	IA_IntPoint	ip;
	while (iter(ip, *d_scan)) {
	    *o_scan = ip[0];
	    for (unsigned i=1; i<dimen; i++) {
		*o_scan *= src_width[i];
		*o_scan += ip[i];
	    }
	    d_scan++;
	    o_scan++;
	}
    }

    IA_IPSIter	iter(dest_ps);
    IA_IntPoint	ip;
    float *const	dest_data = new float[dest_ps.card()];
    float *	valp = dest_data;
    while (iter(ip)) {
	int	offset= ip[0] - src_infimum[0];
	unsigned i;
	for (i=1; i<dimen; i++) {
	    offset *= src_width[i];
	    offset += ip[i] - src_infimum[i];
	}
	const float *const base = src_data + offset;

	// _IVAL_=(base[templ_offsets[i]])
	// _TVAL_=(templ_data[i])
	// _IRESULT_=(*valp)

	float	a=-infinity(); ;

	for (i=0; i<templ_sz; i++) {
	    register float ival = (base[templ_offsets[i]]); if (ival == -infinity()) continue; float temp = ival+(templ_data[i]); if (temp > a) a = temp; ;
	}
	{
	    (*valp) = a; ;
	}
	valp++;
    }
    delete[] templ_offsets;
    delete[] templ_data;

    return IA_FloatDiscreteImage(dest_ps, dest_data, dest_ps.card(), 1);
}

static IA_FloatDiscreteImage
addmax_product_inv(const IA_FloatDiscreteImage &img,
	 const IA_FloatDiscreteImage &invtempl,
	 IA_IntPointSet dest_ps)
{
    IA_IntPoint	inf_ = dest_ps.inf()+invtempl.domain().inf();
    IA_IntPoint	sup_ = dest_ps.sup()+invtempl.domain().sup();
    IA_IntPointSet	src_ps = IA_IntPointSet(inf_,sup_);

    float	*const src_data = new float[src_ps.card()];

    zero_extend(img, src_ps, src_data, -infinity());
    IA_FloatDiscreteImage	rval =
	addmax_product_inv_core(inf_, sup_-inf_ + 1, src_data, invtempl, dest_ps);

    delete[] src_data;

    return rval;
}

IA_FloatDiscreteImage addmax_product(const IA_FloatDiscreteImage &img,
			     const IA_FloatDDTemplate &templ,
			     IA_IntPointSet dest_ps)
{
    if (templ.type() == IA_InvariantDT<IA_FloatDiscreteImage >::s_type()) {
	return addmax_product_inv(img, ((IA_InvariantDT<IA_FloatDiscreteImage >*)templ.bdtp)->value, dest_ps);
    } else {
	float *const	dest_data = new float[dest_ps.card()];
	float	*valp = dest_data;
	IA_IPSIter	dest_iter(dest_ps);
	IA_IntPoint	base_ip;
	while (dest_iter(base_ip)) {
	    IA_FloatDiscreteImage	tv = templ(base_ip);

	    IA_DIPIter<float>	templ_iter(tv);
	    IA_IntPoint	templ_ip;
	    float	templ_val;

	    // _IVAL_=img(ip)
	    // _TVAL_=templ_val
	    // _IRESULT_=(*valp)

	    float	a=-infinity(); ;

	    while ( templ_iter(templ_ip, templ_val) ) {
		const IA_IntPoint	ip = templ_ip; //+base_ip;
		if (! img.domain().contains(ip))
		    continue;
		float temp = img(ip)+templ_val; if (temp > a) a = temp; ;
	    }
	    {
		(*valp) = a; ;
	    }
	    valp++;
	}
	return IA_FloatDiscreteImage(dest_ps, dest_data, dest_ps.card(), 1);
    }
}
// float addmin_product float float  zero=infinity()
// Emacs: -*- C++ -*-

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


//
// $Log:	absorber,v $
// Revision 1.5  93/05/27  11:48:55  thoth
// Copyright Notices
// 
// Revision 1.4  93/04/29  11:20:56  thoth
// Faster(?) extension for VectorDIs.
// 
// Revision 1.3  93/04/17  18:56:44  jnw
// Fixed to match IA_IntPoint membername (dim)
// 
// Revision 1.2  93/04/08  13:21:59  thoth
// internal helper functions are now static.
// 
// Revision 1.1  93/03/18  11:39:42  thoth
// Initial revision
// 

#ifndef zero_extend_scan_float
#define zero_extend_scan_float

static void zero_extend_vec_scan(const IA_SetStructure &ss,
				 const float **src,
				 float **dest,
				 float zero)
{
    for (unsigned i=0; i<ss.nintervals(); i++) {
	IA_ss_interval	temp(ss.retrieve_interval(i));
	if (temp.substructure == IA_SetStructure::FIRST_ONLY) {
	    (*src) += temp.count;
	} else if (temp.substructure == IA_SetStructure::BOTH) {
	    for (unsigned j=0; j<temp.count; j++) {
		*((*dest)++) = *((*src)++);
	    }
	} else if (temp.substructure == IA_SetStructure::SECOND_ONLY) {
	    for (unsigned j=0; j<temp.count; j++)
		*((*dest)++) = zero;
	} else {
	    for (unsigned j=0; j<temp.count; j++)
		zero_extend_vec_scan(temp.substructure, src, dest, zero);
	}
    }
}

static void zero_extend_iter_scan(const IA_SetStructure &ss,
				  IA_DIVIter<float> *srciter,
				  float **dest,
				  float zero)
{
    for (unsigned i=0; i<ss.nintervals(); i++) {
	IA_ss_interval	temp(ss.retrieve_interval(i));
	if (temp.substructure == IA_SetStructure::FIRST_ONLY) {
	    float	blah;
	    for (unsigned j=0; j<temp.count; j++) 
		(*srciter)(blah);
	} else if (temp.substructure == IA_SetStructure::BOTH) {
	    float	blah;
	    for (unsigned j=0; j<temp.count; j++) {
		(*srciter)(blah);
		*((*dest)++) = blah;
	    }
	} else if (temp.substructure == IA_SetStructure::SECOND_ONLY) {
	    for (unsigned j=0; j<temp.count; j++)
		*((*dest)++) = zero;
	} else {
	    for (unsigned j=0; j<temp.count; j++)
		zero_extend_iter_scan(temp.substructure, srciter, dest, zero);
	}
    }
}

static void zero_extend(const IA_FloatDiscreteImage &srcimg,
			const IA_IntPointSet &dest_ps, float *dest_data,
			float zero)
{
    IA_SetStructure	ss;
    intersect_with_dualstruct(srcimg.domain(), dest_ps, &ss);

    if (srcimg.type() == IA_VectorDI<float>::s_type()) {
	float	*srcdata = ((IA_VectorDI<float>*)srcimg.bdip)->vec;
	zero_extend_vec_scan(ss, &srcdata, &dest_data, zero);
    } else {
	IA_DIVIter<float>	iter(srcimg);
	zero_extend_iter_scan(ss, &iter, &dest_data, zero);
    }
}

#endif

static IA_FloatDiscreteImage
addmin_product_inv_core(IA_IntPoint src_infimum,
	      IA_IntPoint src_width,
	      const float *src_data, // length is prod(src_width)
	      const IA_FloatDiscreteImage &templ,
	      IA_IntPointSet dest_ps)
{
    const int	dimen = src_width.dim();

    IA_IntPointSet	templ_ps = templ.domain();
    int		templ_sz = templ_ps.card();
    float	*const templ_data = new float[templ_sz];
    int	*const templ_offsets = new int[templ_sz];

    {
	float	*d_scan = templ_data;
	int	*o_scan = templ_offsets;
	IA_DIPIter<float>	iter(templ);
	IA_IntPoint	ip;
	while (iter(ip, *d_scan)) {
	    *o_scan = ip[0];
	    for (unsigned i=1; i<dimen; i++) {
		*o_scan *= src_width[i];
		*o_scan += ip[i];
	    }
	    d_scan++;
	    o_scan++;
	}
    }

    IA_IPSIter	iter(dest_ps);
    IA_IntPoint	ip;
    float *const	dest_data = new float[dest_ps.card()];
    float *	valp = dest_data;
    while (iter(ip)) {
	int	offset= ip[0] - src_infimum[0];
	unsigned i;
	for (i=1; i<dimen; i++) {
	    offset *= src_width[i];
	    offset += ip[i] - src_infimum[i];
	}
	const float *const base = src_data + offset;

	// _IVAL_=(base[templ_offsets[i]])
	// _TVAL_=(templ_data[i])
	// _IRESULT_=(*valp)

	float	a=infinity(); ;

	for (i=0; i<templ_sz; i++) {
	    register float ival = (base[templ_offsets[i]]); if (ival == infinity()) continue; float temp = ival+(templ_data[i]); if (temp < a) a = temp; ;
	}
	{
	    (*valp) = a; ;
	}
	valp++;
    }
    delete[] templ_offsets;
    delete[] templ_data;

    return IA_FloatDiscreteImage(dest_ps, dest_data, dest_ps.card(), 1);
}

static IA_FloatDiscreteImage
addmin_product_inv(const IA_FloatDiscreteImage &img,
	 const IA_FloatDiscreteImage &invtempl,
	 IA_IntPointSet dest_ps)
{
    IA_IntPoint	inf_ = dest_ps.inf()+invtempl.domain().inf();
    IA_IntPoint	sup_ = dest_ps.sup()+invtempl.domain().sup();
    IA_IntPointSet	src_ps = IA_IntPointSet(inf_,sup_);

    float	*const src_data = new float[src_ps.card()];

    zero_extend(img, src_ps, src_data, infinity());
    IA_FloatDiscreteImage	rval =
	addmin_product_inv_core(inf_, sup_-inf_ + 1, src_data, invtempl, dest_ps);

    delete[] src_data;

    return rval;
}

IA_FloatDiscreteImage addmin_product(const IA_FloatDiscreteImage &img,
			     const IA_FloatDDTemplate &templ,
			     IA_IntPointSet dest_ps)
{
    if (templ.type() == IA_InvariantDT<IA_FloatDiscreteImage >::s_type()) {
	return addmin_product_inv(img, ((IA_InvariantDT<IA_FloatDiscreteImage >*)templ.bdtp)->value, dest_ps);
    } else {
	float *const	dest_data = new float[dest_ps.card()];
	float	*valp = dest_data;
	IA_IPSIter	dest_iter(dest_ps);
	IA_IntPoint	base_ip;
	while (dest_iter(base_ip)) {
	    IA_FloatDiscreteImage	tv = templ(base_ip);

	    IA_DIPIter<float>	templ_iter(tv);
	    IA_IntPoint	templ_ip;
	    float	templ_val;

	    // _IVAL_=img(ip)
	    // _TVAL_=templ_val
	    // _IRESULT_=(*valp)

	    float	a=infinity(); ;

	    while ( templ_iter(templ_ip, templ_val) ) {
		const IA_IntPoint	ip = templ_ip; //+base_ip;
		if (! img.domain().contains(ip))
		    continue;
		float temp = img(ip)+templ_val; if (temp < a) a = temp; ;
	    }
	    {
		(*valp) = a; ;
	    }
	    valp++;
	}
	return IA_FloatDiscreteImage(dest_ps, dest_data, dest_ps.card(), 1);
    }
}
