region.c

seismic software,very useful

    int tmp;

    int data_size;

/*--------------------------------------------------------------------*\
   Initialize the list using a single seed point. Any points in range
   will be added to the list for the next search pass.

   Note that stride1 == 1 so the list[0] is the index into the data
   array of the seed point.

   In an array of constant values, this leads to an expanding shell.
\*--------------------------------------------------------------------*/

    data_size = size1*size2*size3;

    stride1 = 1;
    stride2 = size1;
    stride3 = size1 * size2;

    cell = data;
    old = list;
    endold = list + 1;
    new = list + 1;
    listend = list + nlist;
    list[0] = seed1 * stride1 + seed2 * stride2 + seed3 * stride3;

    debug=0;

    if( debug ){
       printf( "Seed point values:\n"          );
       printf( "   list[0]=%d " ,list[0]       );
       printf( "   min=%d "     ,min           );
       printf( "   max=%d "     ,max           );
       printf( "   value=%d"    ,cell[list[0]] );
       printf( "\n"                            );
    }

    debug = 0;

    region->size = 0;
    RegionClear();
    region->size = 0;
    mark_count = 0;
    RegionMarkBorder(data, size1, size2, size3);
    UNMARK(cell[list[0]]);
   
/*--------------------------------------------------------------------*\
   Calculate the offsets in the data volume relative to the center 
   of a 3 x 3 x 3 cube assuming a 3D data volume w/ dimensions of
   size1 * size2 * size3.  Note that there are 26 neighbors in the
   volume surrounding the initial point.

   The basic logic is to start with a single point and expand the
   search until there are no more points to search, i.e. the points
   are either not in the data range selected or have already been
   marked.

   The preceding call to RegionMarkBorder, marks the edge cells in
   the data to prevent walking out of bounds.
\*--------------------------------------------------------------------*/

    memset( neighbor ,0 ,sizeof( neighbor ) );
    nbr = neighbor;

    /*---------------------------------------------*/
    /* face neighbors (single stride combinations) */
    /*---------------------------------------------*/

    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;


/*--------------------------------------------------------------------*\
   Continue the search so long as there were new points found during
   the last search pass.
\*--------------------------------------------------------------------*/

    while (old < new && new < listend ){

/*--------------------------------------------------------------------*\
   Shift seed values for this search to the start of the list to
   provide room for new values found during this search pass.
\*--------------------------------------------------------------------*/

        if( debug ){
           printf( "%d " ,endold - old );
           printf( "%d " ,new - endold );
        }

        for( endold = new, new = old, old = list; new < endold; ){
            *old++ = *new++;
        }
        count += old - list;

/*--------------------------------------------------------------------*\
   Search all of the cells adjacent to the current list of seed values
   and add any found within range that have not already been marked.

   Since the border cells were marked before we started, the search
   will terminate at the borders of the volume.
\*--------------------------------------------------------------------*/

        for( new = old, endold = old, old = list; 
            old < endold && new < listend && old < listend;
             old++ ){

            if( debug ){
               if( *old > size1*size2*size3 ){
                  printf( "stop" );
               }
            }

/*--------------------------------------------------------------------*\
   Examine the 26 neighbors of the current seed value denoted by the
   old pointer. Any unmarked voxels which are within the data range
   will be marked and placed in the list of voxels whose neighbors
   will be searched in the next pass.
\*--------------------------------------------------------------------*/

            /*--------------------*/
            /* examine face cells */
            /*--------------------*/

            if( neighborhood & MARK_FACE ){
                for( nbr = face1; nbr < face2 && new < listend; nbr++ ){
                    if( *old - *nbr >=0 && *old + *nbr < data_size ){
                        *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++ ){
                    if( *old - *nbr >=0 && *old + *nbr < data_size ){
                        *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++ ){
                    if( *old - *nbr >=0 && *old + *nbr < data_size ){
                        *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 (mark_count);
}

/* mark the borders of the data cube */
void RegionMarkBorder(Buffer data, int size1, int size2, int 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 */
void RegionUnMarkBorder(Buffer data, int size1, int size2, int size3, int Min,
                   int 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 */
void RegionDump(void)
{
    extern REgion region;
    extern Data data;
    string   filename;
    FILE    *fd;
    int      n1, n2, n3, i1, i2, i3, i, ii = 1;
    float   *values0, *values1, *values2, *values3;
    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 program

#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);
        }
\*====================================================================*/