region.c

seismic software,very useful

int min, max, seed1, seed2, seed3, neighborhood;
	{
	int stride1, stride2, stride3, count=0;
	Shadow endold;
	register Shadow old, new, listend;
	register int cell1;
	register Buffer cell;
	int neighbor[26], *face1, *edge1, *corner1;
	register int *nbr, *face2, *edge2, *corner2;

	/* precondition data boundaries */
	region->size = 0;
	RegionClear ();
	region->size = 0;
	RegionMarkBorder (data,size1,size2,size3);
/*
	seed1 = seed1 > 0 ? seed1 : 1;
	seed2 = seed2 > 0 ? seed2 : 1;
	seed3 = seed3 > 0 ? seed3 : 1;
	seed1 = seed1 <= size1 ? seed1 : size1 - 2;
	seed2 = seed2 <= size2 ? seed2 : size2 - 2;
	seed3 = seed3 <= size3 ? seed3 : size3 - 2;
*/
	/* enummerate the neighbor offsets */
	nbr = neighbor;
	/* face neighbors (single stride combinations) */
	stride1 = 1;
	stride2 = size1;
	stride3 = size1 * size2;
	face1 = nbr;
	*nbr++ = stride1;
	*nbr++ = stride2;
	*nbr++ = stride3;
	*nbr++ = -stride1;
	*nbr++ = -stride2;
	*nbr++ = -stride3;
	face2 = nbr;
	/* edge neighbors (two stride combinations) */
	edge1 = nbr;
	*nbr++ = stride1 + stride2;
	*nbr++ = stride1 - stride2;
	*nbr++ = -stride1 + stride2;
	*nbr++ = -stride1 - stride2;
	*nbr++ = stride2 + stride3;
	*nbr++ = stride2 - stride3;
	*nbr++ = -stride2 + stride3;
	*nbr++ = -stride2 - stride3;
	*nbr++ = stride3 + stride1;
	*nbr++ = stride3 - stride1;
	*nbr++ = -stride3 + stride1;
	*nbr++ = -stride3 - stride1;
	edge2 = nbr;
	/* corner neighbors (three stride combinations) */
	corner1 = nbr;
	*nbr++ = stride3 + stride2 + stride1;
	*nbr++ = stride3 + stride2 - stride1;
	*nbr++ = stride3 - stride2 + stride1;
	*nbr++ = stride3 - stride2 - stride1;
	*nbr++ = -stride3 + stride2 + stride1;
	*nbr++ = -stride3 + stride2 - stride1;
	*nbr++ = -stride3 - stride2 + stride1;
	*nbr++ = -stride3 - stride2 - stride1;
	corner2 = nbr;

	cell = data;
	old = list;
	endold = list + 1;
	new = list + 1;
	listend = list + nlist;
	list[0] = seed1 * stride1 + seed2 * stride2 + seed3 * stride3;
	UNMARK(cell[list[0]]);
/*	printf ("list[0]=%d min=%d max=%d value=%d\n",list[0],min,max,cell[list[0]]);*/
	while (old<new && new<listend) {
/*		printf ("old=%d endold=%d new=%d\n",old-list,endold-list,new-list);*/
		/* make room in list by slide new entries over old */
		for (endold=new, new=old, old=list; new<endold;) *old++ = *new++;
		count += old-list;
		for (new=old, endold=old, old=list; old<endold && new<listend; old++) {
			/* examine face cells */
			if (neighborhood & MARK_FACE) {
				for (nbr=face1; nbr<face2 && new<listend; nbr++) {
					*new = *old + *nbr;
					cell1 = cell[*new];
					if (INRANGE(cell1)) {
						MARK(cell[*new++]);
						}
					}
				}
			/* examine edge cells */
			if (neighborhood & MARK_EDGE) {
				for (nbr=edge1; nbr<edge2 && new<listend; nbr++) {
					*new = *old + *nbr;
					cell1 = cell[*new];
					if (INRANGE(cell1)) {
						MARK(cell[*new++]);
						}
					}
				}
			/* examine corner cells */
			if (neighborhood & MARK_CORNER) {
				for (nbr=corner1; nbr<corner2 && new<listend; nbr++) {
					*new = *old + *nbr;
					cell1 = cell[*new];
					if (INRANGE(cell1)) {
						MARK(cell[*new++]);
						}
					}
				}
			}
		}
	/* unmark boundaries */
	RegionUnMarkBorder (data,size1,size2,size3,min,max);
	/* return list size */
	return (count);
	}

/* mark the borders of the data cube */
RegionMarkBorder (data,size1,size2,size3)
Buffer data;
int size1, size2, size3;
	{
	Buffer d3, d2, e2;
	register Buffer d1, e1;
	int stride1, stride2, stride3, sz3, sz2, sz1, strd3, strd2;
	register int strd1;

	stride1 = 1;
	stride2 = size1;
	stride3 = size1 * size2;
	/* n3 edges */
	sz3 = size3;
	sz2 = size2;
	sz1 = size1;
	strd3 = stride3;
	strd2 = stride2;
	strd1 = stride1;
	d3 = data;
	for (d2=d3, e2=d2+sz2*strd2; d2<e2; d2+=strd2)
	for (d1=d2, e1=d1+sz1*strd1; d1<e1; d1+=strd1)
		MARK (*d1);
	d3 = data + (sz3-1)*strd3;
	for (d2=d3, e2=d2+sz2*strd2; d2<e2; d2+=strd2)
	for (d1=d2, e1=d1+sz1*strd1; d1<e1; d1+=strd1)
		MARK (*d1);

	/* n2 edges */
	sz3 = size2;
	sz2 = size1;
	sz1 = size3;
	strd3 = stride2;
	strd2 = stride1;
	strd1 = stride3;
	d3 = data;
	for (d2=d3, e2=d2+sz2*strd2; d2<e2; d2+=strd2)
	for (d1=d2, e1=d1+sz1*strd1; d1<e1; d1+=strd1)
		MARK (*d1);
	d3 = data + (sz3-1)*strd3;
	for (d2=d3, e2=d2+sz2*strd2; d2<e2; d2+=strd2)
	for (d1=d2, e1=d1+sz1*strd1; d1<e1; d1+=strd1)
		MARK (*d1);

	/* n1 edges */
	sz3 = size1;
	sz2 = size3;
	sz1 = size2;
	strd3 = stride1;
	strd2 = stride3;
	strd1 = stride2;
	d3 = data;
	for (d2=d3, e2=d2+sz2*strd2; d2<e2; d2+=strd2)
	for (d1=d2, e1=d1+sz1*strd1; d1<e1; d1+=strd1)
		MARK (*d1);
	d3 = data + (sz3-1)*strd3;
	for (d2=d3, e2=d2+sz2*strd2; d2<e2; d2+=strd2)
	for (d1=d2, e1=d1+sz1*strd1; d1<e1; d1+=strd1)
		MARK (*d1);
	}

/* unmark border of voxel cube */
RegionUnMarkBorder (data,size1,size2,size3,Min,Max)
Buffer data;
int size1, size2, size3, Min, Max;
	{
	Buffer d3, d2, e2;
	register Buffer d1, e1;
	int stride1, stride2, stride3, sz3, sz2, sz1, strd2;
	register int strd1, strd3, min, max;

	stride1 = 1;
	stride2 = size1;
	stride3 = size1 * size2;
	min = Min;
	max = Max;
	/* n3 edges */
	sz3 = size3;
	sz2 = size2;
	sz1 = size1;
	strd3 = stride3;
	strd2 = stride2;
	strd1 = stride1;
	d3 = data;
	for (d2=d3, e2=d2+sz2*strd2; d2<e2; d2+=strd2)
	for (d1=d2, e1=d1+sz1*strd1; d1<e1; d1+=strd1)
		{if (NOTMARK(*(d1+strd3)) || OUTRANGE(*d1)) UNMARK(*d1);}
	d3 = data + (sz3-1)*strd3;
	for (d2=d3, e2=d2+sz2*strd2; d2<e2; d2+=strd2)
	for (d1=d2, e1=d1+sz1*strd1; d1<e1; d1+=strd1)
		{if (NOTMARK(*(d1-strd3)) || OUTRANGE(*d1)) UNMARK(*d1);}

	/* n2 edges */
	sz3 = size2;
	sz2 = size1;
	sz1 = size3;
	strd3 = stride2;
	strd2 = stride1;
	strd1 = stride3;
	d3 = data;
	for (d2=d3, e2=d2+sz2*strd2; d2<e2; d2+=strd2)
	for (d1=d2, e1=d1+sz1*strd1; d1<e1; d1+=strd1)
		{if (NOTMARK(*(d1+strd3)) || OUTRANGE(*d1)) UNMARK(*d1);}
	d3 = data + (sz3-1)*strd3;
	for (d2=d3, e2=d2+sz2*strd2; d2<e2; d2+=strd2)
	for (d1=d2, e1=d1+sz1*strd1; d1<e1; d1+=strd1)
		{if (NOTMARK(*(d1-strd3)) || OUTRANGE(*d1)) UNMARK(*d1);}

	/* n1 edges */
	sz3 = size1;
	sz2 = size3;
	sz1 = size2;
	strd3 = stride1;
	strd2 = stride3;
	strd1 = stride2;
	d3 = data;
	for (d2=d3, e2=d2+sz2*strd2; d2<e2; d2+=strd2)
	for (d1=d2, e1=d1+sz1*strd1; d1<e1; d1+=strd1)
		{if (NOTMARK(*(d1+strd3)) || OUTRANGE(*d1)) UNMARK(*d1);}
	d3 = data + (sz3-1)*strd3;
	for (d2=d3, e2=d2+sz2*strd2; d2<e2; d2+=strd2)
	for (d1=d2, e1=d1+sz1*strd1; d1<e1; d1+=strd1)
		{if (NOTMARK(*(d1-strd3)) || OUTRANGE(*d1)) UNMARK(*d1);}
	}

/* write region coordinates and values to a file */
RegionDump ()
	{
	extern REgion region;
	extern Data data;
	string filename;
	FILE *fd;
	int n1, n2, n3, i1, i2, i3, i, ii=1, value;
	float *values0, *values1, *values2, *values3;
	register Buffer save;
	register Buffer buf, datap, edata;

	if (region->size <= 0) return;
	n1 = AxisSize(DataAxis(data,DATA_AXIS1));
	n2 = AxisSize(DataAxis(data,DATA_AXIS2));
	n3 = AxisSize(DataAxis(data,DATA_AXIS3));
	values0 = AxisValues(DataAxis(data,DATA_VALUE));
	values1 = AxisValues(DataAxis(data,DATA_AXIS1));
	values2 = AxisValues(DataAxis(data,DATA_AXIS2));
	values3 = AxisValues(DataAxis(data,DATA_AXIS3));
	sprintf (filename,"region.list.%d.txt",region->size);
	fd = fopen (filename,"w");
	fprintf (fd,"%5s %5s %5s %5s %5s %10s %10s %10s %10s\n","i","i3","i2","i1","i0","v3","v2","v1","v0");
	fprintf (fd,"%5d %5d %5d %5d %5d %10g %10g %10g %10g\n",0,0,0,0,region->bound[0],values3[0],values2[0],values1[0],values0[region->bound[0]]);
	fprintf (fd,"%5d %5d %5d %5d %5d %10g %10g %10g %10g\n",region->size,n3,n2,n1,region->bound[1],values3[n3-1],values2[n2-1],values1[n1-1],values0[region->bound[1]]);
	buf = DataBuffer(data);
	for (datap=buf, edata=datap+DataSize(data); datap<edata; datap++) {
		if (ISMARK(*datap)) {
			i = datap - buf;
			i1 = i % n1;
			i = i / n1;
			i2 = i % n2;
			i3 = i / n2;
			fprintf (fd,"%5d %5d %5d %5d %5d %10g %10g %10g %10g\n",ii++,i3,i2,i1,(*datap) & 0x7f,values3[i3],values2[i2],values1[i1],values0[(*datap) & 0x7f]);
			}
		}
	fclose (fd);
	UIMessage ("region dumped to file");
	}

/* test
#define N 16
unsigned char x[16][16][16];
#define	NLIST	1024
int list[NLIST];
main () {
	int i, j, k;
	for (i=4; i<12; i++)
	for (j=4; j<12; j++)
	for (k=4; k<12; k++)
		x[i][j][k] = 8;
	i = RegionMark (x,N,N,N,6,9,2180,MARK_FACE|MARK_CORNER,list,NLIST);
	printf ("nlist=%d\n",i);
	write (creat("tmp",0644),x,N*N*N);
	}
*/