head	1.8;
access
	monica
	thoth
	ljr;
symbols
	Eglin_alpha_1:1.1;
locks
	monica:1.8; strict;
comment	@// @;


1.8
date	94.12.23.00.28.57;	author ljr;	state Exp;
branches;
next	1.7;

1.7
date	94.12.10.20.32.04;	author monica;	state Exp;
branches;
next	1.6;

1.6
date	94.12.08.22.30.29;	author monica;	state Exp;
branches;
next	1.5;

1.5
date	94.12.08.04.59.46;	author monica;	state Exp;
branches;
next	1.4;

1.4
date	94.11.09.20.25.21;	author monica;	state Exp;
branches;
next	1.3;

1.3
date	94.11.02.22.51.27;	author monica;	state Exp;
branches;
next	1.2;

1.2
date	94.08.17.04.12.11;	author monica;	state Exp;
branches;
next	1.1;

1.1
date	93.01.12.14.48.25;	author thoth;	state Exp;
branches;
next	;


desc
@more complete fast-convolution method (experimental)
@


1.8
log
@Added compat.h include.
Added instantiation directive for IA_Array<frame> for BC++.
Added workaround for temporary in a conditional expression.
@
text
@// Emacs: -*- C++ -*-
static char rcsid[] = "$Id: test015.c,v 1.7 1994/12/10 20:32:04 monica Exp ljr $";

//
// $Log: test015.c,v $
// Revision 1.7  1994/12/10  20:32:04  monica
// reincluded some .h files
//
// Revision 1.6  1994/12/08  22:30:29  monica
// removed refs to ImageDisplay
//
// Revision 1.5  1994/12/08  04:59:46  monica
// removed references to timing
//
// Revision 1.4  1994/11/09  20:25:21  monica
// removed include devfrags & any debugging calls to malloc
//
// Revision 1.3  1994/11/02  22:51:27  monica
// removed display refs
//
// Revision 1.2  1994/08/17  04:12:11  monica
// added int main & return 0
//
// Revision 1.1  1993/01/12  14:48:25  thoth
// Initial revision
//
// Revision 1.1  92/10/07  22:46:07  thoth
// Initial revision
// 
//

 #include        <iostream.h>
 #include        <stdio.h>
 #include        <stdlib.h>

// #include <sys/time.h>
// #include <sys/resource.h>

#include	"IntPS.h"
#include	"PSIter.h"
// #include "ImageDisplayProtocol.h"
// #include "tutil.h"

#include "compat.h"

IA_Set<IA_Point<int> >
slurp_pgm_file(const char *fname, unsigned char **data, unsigned *maxval)
{
    FILE	*fp = fopen(fname, "r");
    if (!fp) {
	cerr << "error opening PGM file " << fname << "\n";
	return IA_Set<IA_Point<int> >();
    }

    unsigned width, height;
    if (3!=fscanf(fp, "P5 %u %u %u", &width, &height, maxval)) {
        cerr << "Error: was expecting raw PGM (P5) image in file "
             << fname << "\n";
	fclose(fp);
        return IA_Set<IA_Point<int> >();
    }

    int ch;
    if ((ch=getc(fp))!='\n')
        ungetc(ch, fp);

    const unsigned size = width*height;

    *data = new unsigned char [size];

    if (!fread(*data, size, 1, fp)) {
	cerr << "short data read on PGM file " << fname << "\n";
	fclose(fp);
        exit(1);
    }

    fclose(fp);

    return IA_Set<IA_Point<int> >(IA_Point<int>(0,0), IA_Point<int>(height-1, width-1));
}

struct frame {
  IA_SetStructure	ssp;	// what SetStructure are we scanning through?
  unsigned	idx;	// our index into this SetStructure
  unsigned	howmany; // how many times we should repeat this SS
  unsigned	count;	 // how many times we already have
  frame() { idx=howmany=count=0; }
};

#ifdef __BCPLUSPLUS__
#include "Array.h"
#include "Array.c"
typedef IA_Array<frame> IA_array_frame_bcpp_pstest015;
#endif

void almost_fast_boxy_conv
(IA_Set<IA_Point<int> > src_ps, const unsigned char *src_data,
 IA_Set<IA_Point<int> > templ_ps, const double *templ_data,
 IA_Set<IA_Point<int> > dest_ps, unsigned char **dest_data)
{
    *dest_data=0;
    if (!src_ps.boxy()) {
	cerr << "source pointset not boxy\n";
	return;
    }

    if (dest_ps.dim() != src_ps.dim()) {
	cerr << "dimension mismatch\n";
	return;
    }

    // calculate the pointset that we will iterate over the fastest
    IA_Set<IA_Point<int> >	quickly;
    IA_SetStructure	ss;
    IA_Point<int>	inf_(src_ps.inf()-templ_ps.inf());
    {
	IA_Point<int>	sup_(src_ps.sup()-templ_ps.sup());

#ifndef IA_NO_TEMPS_IN_CONDITIONAL_EXPRS
	quickly = intersect_with_structure
	    (dest_ps,
	     (sup_<inf_) ? IA_BlackHole(src_ps.dim()) :IA_Set<IA_Point<int> >(inf_,sup_),
	     &ss);
#else
    if (sup_ < inf_)
        quickly = intersect_with_structure(dest_ps, IA_BlackHole(src_ps.dim()), &ss);
    else
        quickly = intersect_with_structure(dest_ps, IA_Set<IA_Point<int> >(inf_, sup_), &ss);
#endif
	if (ss.nintervals()<=0) {
	    cout << "too few intervals" << endl;
	}
    }

    // allocate space for the result image data
    *dest_data = new unsigned char[dest_ps.card()];

    // precalculate offsets for the invariant template pixels
    const int n=templ_ps.card();
    int *offsets = new int [n];
    IA_Point<int>	*points = new IA_Point<int> [n];
    {
	unsigned	base = src_ps.index(inf_);

	IA_PSIter<IA_Point<int> >	iter(templ_ps);
	IA_Point<int>	ip;

	unsigned i=0;
	while (iter(ip)) {
	    points[i] = ip;
	    offsets[i] = (int)src_ps.index(inf_+ip) - (int)base;
	    i++;
	}
    }

    IA_PSIter<IA_Point<int> > iter(dest_ps);
    IA_Point<int>	ip;
    unsigned char	*valp=*dest_data;

    // I wanted to declare frame here, but NOOOO.  CFront pukes.
    IA_Array<frame>	stack;
    stack[0].ssp=ss;
    stack[0].idx=0;
    stack[0].howmany=1;
    stack[0].count=0;
    int	stackp=0;
    while (stackp>=0) {
	frame	&curr=stack[stackp];
	if (curr.idx >= curr.ssp.nintervals()) {
	    // cout << curr.idx << " " << curr.ssp.nintervals() << endl;
	    if (++curr.count>=curr.howmany)
		stackp--;
	    else
		curr.idx=0;
	    continue;
	}
	IA_ss_interval	ssi(curr.ssp.retrieve_interval(curr.idx));
	if (ssi.substructure == IA_SetStructure::FIRST_ONLY) {
	    for (unsigned i=0; i<ssi.count; i++) {
		iter(ip);
		// next line is erroneous (index)
		const unsigned char	*base = src_data + src_ps.index(ip);
		double	a=0;
		for (unsigned i=0; i<n; i++) {
		    if (src_ps.contains(ip+points[i]))
			a += templ_data[i]*base[offsets[i]];
		}
		*(valp++) = (unsigned char)a;
	    }
	    curr.idx++;
	} else if (ssi.substructure == IA_SetStructure::BOTH) {
	    for (unsigned i=0; i<ssi.count; i++) {
		iter(ip);
		const unsigned char	*base = src_data + src_ps.index(ip);
		double	a=0;
		for (unsigned i=0; i<n; i++) {
		    a += templ_data[i]*base[offsets[i]];
		}
		*(valp++) = (unsigned char)a;
	    }
	    curr.idx++;
	} else {
	    curr.idx++;
	    stackp++;
	    frame	&newf = stack[stackp];
	    newf.ssp = ssi.substructure;
	    newf.idx = 0;
	    newf.howmany = ssi.count;
	    newf.count = 0;
	    continue;
	}
    }

    delete[] offsets;
    delete[] points;
}

int main(int argc, char**argv)
{
//    char	*disp=getenv("DISPLAY");
//    if (disp==0) {
//	cerr << "No DISPLAY environment variable, hasta.\n";
//	exit(1);
//    }

//    if (argc!=2) {
//	cerr << "Usage: "<<argv[0]<<" pgm-file\n";
//	exit(1);
//    }

//    const char	*fname=argv[1];
    const char	*fname="sr71_noisy.pgm";

    unsigned char	*data;
    unsigned	maxval;
    IA_Set<IA_Point<int> >	domain;
    domain = slurp_pgm_file(fname, &data, &maxval);

    const IA_Point<int>	gauss1[] = {
	IA_Point<int>(-3,0),
	    IA_Point<int>(-2,0),
	    IA_Point<int>(-1,0),
	    IA_Point<int>(0, 0),
	    IA_Point<int>( 1,0),
	    IA_Point<int>( 2,0),
	    IA_Point<int>( 3,0),
	};
    const IA_Point<int>	gauss2[] = {
	IA_Point<int>(0,-3),
	    IA_Point<int>(0,-2),
	    IA_Point<int>(0,-1),
	    IA_Point<int>(0, 0),
	    IA_Point<int>(0, 1),
	    IA_Point<int>(0, 2),
	    IA_Point<int>(0, 3),
	};
    static const double	gvalues[] = {
	0.00443,
	    0.05399,
	    0.24197,
	    0.39894,
	    0.24197,
	    0.05399,
	    0.00443,
	};

    IA_Set<IA_Point<int> >	templ(2,gauss1,7);

//    IDP::ClientID	cid;
//    cid = IDP::NewClient("iac++ test013");
//    IDP::LinkID		lid;
//    struct rusage	now, later;

    unsigned char	*data1;
    unsigned char	*data2;
//    double dt;
    IA_Array<unsigned>	dims(2);

//    lid = IDP::CreateDisplay(cid, disp, IDP::READ_WRITE);

//    getrusage(RUSAGE_SELF ,&now);
    templ = IA_Set<IA_Point<int> >(2,gauss1,7);
    almost_fast_boxy_conv(domain, data, templ, gvalues, domain, &data1);
//    getrusage(RUSAGE_SELF ,&later);

//    dt = (Ttof(later.ru_utime)- Ttof(now.ru_utime));
//    cout << "first convolution took a total of " << dt << " user seconds\n";

//    getrusage(RUSAGE_SELF ,&now);
    templ = IA_Set<IA_Point<int> >(2,gauss2,7);
    almost_fast_boxy_conv(domain, data1, templ, gvalues, domain, &data2);
//    getrusage(RUSAGE_SELF ,&later);

//    dt = (Ttof(later.ru_utime)- Ttof(now.ru_utime));
//    cout << "second convolution took a total of " << dt << " user seconds\n";

    {
	IA_Point<int>	width = domain.sup() - domain.inf();
	dims[0] = width[0]+1;
	dims[1] = width[1]+1;
    }
//    sleep(5);
//    IDP::WriteImage(lid, "gaussed twice", IDP::U_BYTE, 1, dims, data2);
    delete[] data1;
    delete[] data2;

    return 0;

}
@


1.7
log
@reincluded some .h files
@
text
@d2 1
a2 1
static char rcsid[] = "$Id: test015.c,v 1.6 1994/12/08 22:30:29 monica Exp monica $";
d6 3
d32 3
a34 3
 #include	<iostream.h>
 #include	<stdio.h>
 #include	<stdlib.h>
d44 2
d90 6
d119 1
d124 6
d302 1
a302 1
    sleep(5);
@


1.6
log
@removed refs to ImageDisplay
@
text
@d2 1
a2 1
static char rcsid[] = "$Id: test015.c,v 1.5 1994/12/08 04:59:46 monica Exp monica $";
d6 3
d29 3
a31 3
// #include	<iostream.h>
// #include	<stdio.h>
// #include	<stdlib.h>
@


1.5
log
@removed references to timing
@
text
@d2 1
a2 1
static char rcsid[] = "$Id: test015.c,v 1.4 1994/11/09 20:25:21 monica Exp monica $";
d6 3
d35 1
a35 1
#include "ImageDisplayProtocol.h"
d199 5
a203 5
    char	*disp=getenv("DISPLAY");
    if (disp==0) {
	cerr << "No DISPLAY environment variable, hasta.\n";
	exit(1);
    }
@


1.4
log
@removed include devfrags & any debugging calls to malloc
@
text
@d2 1
a2 1
static char rcsid[] = "$Id: test015.c,v 1.3 1994/11/02 22:51:27 monica Exp monica $";
d6 3
d23 3
a25 3
#include	<iostream.h>
#include	<stdio.h>
#include	<stdlib.h>
d27 2
a28 2
#include <sys/time.h>
#include <sys/resource.h>
d33 1
a33 3
#include "tutil.h"

/* #include "with-ia-pointset.h" */
d248 1
a248 1
    struct rusage	now, later;
d252 1
a252 1
    double dt;
d257 1
a257 1
    getrusage(RUSAGE_SELF ,&now);
d260 1
a260 1
    getrusage(RUSAGE_SELF ,&later);
d262 2
a263 2
    dt = (Ttof(later.ru_utime)- Ttof(now.ru_utime));
    cout << "first convolution took a total of " << dt << " user seconds\n";
d265 1
a265 1
    getrusage(RUSAGE_SELF ,&now);
d268 1
a268 1
    getrusage(RUSAGE_SELF ,&later);
d270 2
a271 2
    dt = (Ttof(later.ru_utime)- Ttof(now.ru_utime));
    cout << "second convolution took a total of " << dt << " user seconds\n";
@


1.3
log
@removed display refs
@
text
@d2 1
a2 1
static char rcsid[] = "$Id: test015.c,v 1.2 1994/08/17 04:12:11 monica Exp monica $";
d6 3
a142 6
#if 0
	struct frob {
           ~frob() { malloc_verify(); }
        } goop;
#endif

a194 2
    // malloc_debug(2);

@


1.2
log
@added int main & return 0
@
text
@d2 1
a2 1
static char rcsid[] = "$Id: test015.c,v 1.1 1993/01/12 14:48:25 thoth Exp monica $";
d6 3
d249 3
a251 3
    IDP::ClientID	cid;
    cid = IDP::NewClient("iac++ test013");
    IDP::LinkID		lid;
d259 1
a259 1
    lid = IDP::CreateDisplay(cid, disp, IDP::READ_WRITE);
d283 1
a283 1
    IDP::WriteImage(lid, "gaussed twice", IDP::U_BYTE, 1, dims, data2);
@


1.1
log
@Initial revision
@
text
@d2 1
a2 1
static char rcsid[] = "$Id: test009.c,v 1.1 92/10/07 22:46:07 thoth Exp Locker: thoth $";
d5 4
a8 1
// $Log:	test009.c,v $
d21 2
a22 2
#include	"IntPointSet.h"
#include	"IPSIter.h"
d26 1
a26 1
#include "with-ia-pointset.h"
d28 1
a28 1
IntPointSet
d34 1
a34 1
	return IntPointSet();
d42 1
a42 1
        return IntPointSet();
d61 1
a61 1
    return IntPointSet(IntPoint(0,0), IntPoint(height-1, width-1));
d65 1
a65 1
  SetStructure	ssp;	// what SetStructure are we scanning through?
d73 3
a75 3
(IntPointSet src_ps, const unsigned char *src_data,
 IntPointSet templ_ps, const double *templ_data,
 IntPointSet dest_ps, unsigned char **dest_data)
d89 3
a91 3
    IntPointSet	quickly;
    SetStructure	ss;
    IntPoint	inf_(src_ps.inf()-templ_ps.inf());
d93 1
a93 1
	IntPoint	sup_(src_ps.sup()-templ_ps.sup());
d97 1
a97 1
	     (sup_<inf_) ? BlackHole(src_ps.dim()) :IntPointSet(inf_,sup_),
d110 1
a110 1
    IntPoint	*points = new IntPoint [n];
d114 2
a115 2
	IPSIter	iter(templ_ps);
	IntPoint	ip;
d125 2
a126 2
    IPSIter iter(dest_ps);
    IntPoint	ip;
d130 1
a130 1
    Array<frame>	stack;
d152 2
a153 2
	ss_interval	ssi(curr.ssp.retrieve_interval(curr.idx));
	if (ssi.substructure == SetStructure::FIRST_ONLY) {
d166 1
a166 1
	} else if (ssi.substructure == SetStructure::BOTH) {
d203 4
a206 4
    if (argc!=2) {
	cerr << "Usage: "<<argv[0]<<" pgm-file\n";
	exit(1);
    }
d208 2
a209 1
    const char	*fname=argv[1];
d213 1
a213 1
    IntPointSet	domain;
d216 8
a223 8
    const IntPoint	gauss1[] = {
	IntPoint(-3,0),
	    IntPoint(-2,0),
	    IntPoint(-1,0),
	    IntPoint(0, 0),
	    IntPoint( 1,0),
	    IntPoint( 2,0),
	    IntPoint( 3,0),
d225 8
a232 8
    const IntPoint	gauss2[] = {
	IntPoint(0,-3),
	    IntPoint(0,-2),
	    IntPoint(0,-1),
	    IntPoint(0, 0),
	    IntPoint(0, 1),
	    IntPoint(0, 2),
	    IntPoint(0, 3),
d244 1
a244 1
    IntPointSet	templ(2,gauss1,7);
d254 1
a254 1
    Array<unsigned>	dims(2);
d259 1
a259 1
    templ = IntPointSet(2,gauss1,7);
d267 1
a267 1
    templ = IntPointSet(2,gauss2,7);
d275 1
a275 1
	IntPoint	width = domain.sup() - domain.inf();
d283 3
@
