region.c

seismic software,very useful

/*--------------------------------------------------------------------*\
   Identifies and processes continginous region between bounds
   starting at a seed. There are three definitions of continginous:
   6 faces +- 12 edges +- 8 corners. The eighth bit of the data
   bytes stores the region mark.
\*--------------------------------------------------------------------*/

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

#include "par.h"

#include "ui_window.h"
#include "ui_menu.h"
#include "main.h"
#include "axis.h"
#include "map.h"
#include "data.h"
#include "render.h"
#include "region.h"
#include "view.h"

REgion   region;

int mark_count = 0;

/* initialize region object */
void RegionInit(void)
{
    extern Data data;
    extern REgion region;
    int      n1, n2, n3, iaxis;
    memwatch = 0;

    {   extern int _alloc;

        region = (REgion) calloc((1) ,sizeof(region[0]));
        _alloc += (1) * sizeof(region[0]);
        if( region == 0 ){
            err("cant allocate %d bytes for  region; %d already allocated",
                (1) * sizeof(region[0]), _alloc);
        }
        if( memwatch ){
            (void) printf("malloc %s=%d\n", " region", (1) * sizeof(region[0]));
        }
    };
    region->live = 1;
    region->neighborhood = MARK_FACE | MARK_EDGE | MARK_CORNER;
    region->bound[0] = 0;
    region->bound[1] = 0;
    for( iaxis = DATA_AXIS0; iaxis <= DATA_AXIS3; iaxis++ ){
        region->seed[iaxis] = NO_INDEX;
    }
    n1 = AxisSize(DataAxis(data, DATA_AXIS0));
    n2 = AxisSize(DataAxis(data, DATA_AXIS1));
    n3 = AxisSize(DataAxis(data, DATA_AXIS2));
    region->nlist = 4 * (n1 * n2 + n2 * n3 + n3 * n1);
    {
        extern int _alloc;

        region->list =
            (Shadow) calloc((region->nlist) ,sizeof(region->list[0]));
            _alloc += (region->nlist) * sizeof(region->list[0]);
        if( region->list == 0 ){
            err
                ("cant allocate %d bytes for  region->list; %d already allocated",
                 (region->nlist) * sizeof(region->list[0]), _alloc);
        }
        if( memwatch ){
            (void) printf("malloc %s=%d\n", " region->list",
                          (region->nlist) * sizeof(region->list[0]));
        }
    };
    region->size = 0;
    region->save = 0;
    region->index = 0;
    memwatch=0;
}

/* turn on region picking callback */
void RegionSetLive(int live)
{
    extern REgion region;

    if( !region ){
        return;
    }
    region->live = live;
}

/* set bound */
void RegionSetBound(int index, int ibound)
{
    extern REgion region;

    if( !region || ibound < 0 || ibound > 1 ){
        return;
    }
    region->bound[ibound] = index;
}

/* set region neighborhood mode callback */
void RegionSetNeighborhood(int mode)
{
    extern REgion region;

    if( !region ){
        return;
    }
    if( region->neighborhood != mode ){
        region->neighborhood = mode;
        RegionMark0();
    }
}

/* set region seed callback */
void RegionSetSeed(int *seed)
{
    int      iaxis;
    extern REgion region;

    if( !region || !region->live ){
        return;
    }

    for( iaxis = DATA_AXIS0; iaxis <= DATA_AXIS3; iaxis++ ){
        region->seed[iaxis] = seed[iaxis];
    }
    UNMARK(region->seed[DATA_VALUE]);
}

/* return bound */
int RegionBound(int ibound)
{
    if( !region || ibound < 0 || ibound > 1 ){
        return (NO_INDEX);
    }
    return (region->bound[ibound]);
}

/* return region size */
int RegionSize(void)
{
    if( !region ){
        return (NO_INDEX);
    }
    return (region->size);
}

/* mark a region callback */
void RegionMark0(void)
{
    extern REgion region;
    extern Data data;
    Message  message;

    if( !region || !region->live || region->seed[DATA_AXIS1] == NO_INDEX ){
        return;
    }

    region->size = RegionMark(DataBuffer(data),
                              AxisSize(DataAxis(data, DATA_AXIS1)),
                              AxisSize(DataAxis(data, DATA_AXIS2)),
                              AxisSize(DataAxis(data, DATA_AXIS3)),
                              region->seed[DATA_AXIS1],
                              region->seed[DATA_AXIS2],
                              region->seed[DATA_AXIS3],
                              region->bound[0] <
                              region->bound[1] ? region->
                              bound[0] : region->bound[1],
                              region->bound[0] >
                              region->bound[1] ? region->
                              bound[0] : region->bound[1], region->neighborhood,
                              region->list, region->nlist);
    if( region->size ){
        ViewDrawAll();
    }
    sprintf(message, "Region: min=%g max=%g size=%d",
            DataValue(data, region->bound[0]),
            DataValue(data, region->bound[1]), region->size);
    UIMessage(message);
}

/* fill a region with a value callback */
void RegionFill(int *seed, int bound1, int bound2, int value)
{
    int      iaxis;
    extern REgion region;
    extern Data data;
    Message  message;

    if( !region ){
        return;
    }
    for( iaxis = DATA_AXIS0; iaxis <= DATA_AXIS3; iaxis++ ){
        region->seed[iaxis] = seed[iaxis];
    }
    region->bound[0] = bound1;
    region->bound[1] = bound2;
    region->size = RegionMark(DataBuffer(data),
                              AxisSize(DataAxis(data, DATA_AXIS1)),
                              AxisSize(DataAxis(data, DATA_AXIS2)),
                              AxisSize(DataAxis(data, DATA_AXIS3)),
                              region->seed[DATA_AXIS1],
                              region->seed[DATA_AXIS2],
                              region->seed[DATA_AXIS3],
                              region->bound[0],
                              region->bound[1],
                              region->neighborhood,
                              region->list, region->nlist);
    if( region->size ){
        RegionSetValue(value, value);
        RegionClear();
    }
    sprintf(message, "Region: min=%g max=%g size=%d",
            DataValue(data, region->bound[0]),
            DataValue(data, region->bound[1]), region->size);
    UIMessage(message);
}

/* clear region marks in data */
void RegionClear(void)
{
    extern Data data;
    register Buffer datap, edata;

    if( !region || !data ){
        return;
    }

    for( datap = DataBuffer(data), edata = datap + DataSize(data);
         datap < edata; ){
        UNMARK(*datap++);
    }
    if( region->size > 0 ){
        ViewDrawAll();
    }
}

/* set all marked cells to value */
void RegionSetValue(int bound1, int bound2)
{
    extern Data data;
    unsigned char _table[256];
    register Buffer datap, edata, savep, buf, table;
    register Shadow indexp;
    int      ibound;
    memwatch = 0;

    if( !region || !data || bound1 == NO_INDEX ){
        return;
    }

    /*------------------*/
    /* initialize table */
    /*------------------*/

    if( region->bound[0] < region->bound[1] ){
        for( ibound = region->bound[0]; ibound <= region->bound[1]; ibound++ ){
            _table[ibound] =
                bound1 + (ibound * (bound2 - bound1)) / (region->bound[1] -
                                                         region->bound[0]);
        }
    } else if( region->bound[0] > region->bound[1] ){
        for( ibound = region->bound[1]; ibound <= region->bound[0]; ibound++ ){
            _table[ibound] =
                bound1 + (ibound * (bound2 - bound1)) / (region->bound[1] -
                                                         region->bound[0]);
        }
    } else
        _table[region->bound[0]] = bound2;
    table = _table;

    if( region->index ){
        free(region->index);
        region->index = 0;
        if( memwatch ){
            printf("free %s\n", "region->index");
        }
    };
    if( region->save ){
        free(region->save);
        region->save = 0;
        if( memwatch ){
            printf("free %s\n", "region->save");
        }
    };
    {
        extern int _alloc;

        region->index =
            (Shadow) calloc((region->size) ,sizeof(region->index[0]));
            _alloc += (region->size) * sizeof(region->index[0]);
        if( region->index == 0 ){
            err
                ("cant allocate %d bytes for  region->index; %d already allocated",
                 (region->size) * sizeof(region->index[0]), _alloc);
        }
        if( memwatch ){
            (void) printf("malloc %s=%d\n", " region->index",
                          (region->size) * sizeof(region->index[0]));
        }
    };
    {
        extern int _alloc;

        region->save = (byte) calloc((region->size) ,sizeof(region->save[0]));
        _alloc += (region->size) * sizeof(region->save[0]);
        if( region->save == 0 ){
            err
                ("cant allocate %d bytes for  region->save; %d already allocated",
                 (region->size) * sizeof(region->save[0]), _alloc);
        }
        if( memwatch ){
            (void) printf("malloc %s=%d\n", " region->save",
                          (region->size) * sizeof(region->save[0]));
        }
    };
    buf = DataBuffer(data);
    indexp = region->index;
    savep = region->save;
    mark_count = 0;
    for( datap = buf, edata = datap + DataSize(data); datap < edata; datap++ ){
        if( ISMARK(*datap) ){
            *indexp++ = datap - buf;
            *savep = (*datap) & 0x7F;
            *datap = table[*savep++];
        }
    }
    memwatch = 0;
}

/* undo last region smoothing */
void RegionRestoreValue(void)
{
    extern REgion region;
    extern Data data;
    register Buffer datap, savep, endp;
    register Shadow indexp;
    memwatch = 0;
    if( !region || !region->save ){
        return;
    }
    if( region->size == 0 ){
        UIMessage("no smoothing or can only unsmooth once");
        return;
    }
    indexp = region->index;
    datap = DataBuffer(data);
    savep = region->save;
    endp = savep + region->size;
    for( ; savep < endp; ){
        datap[*indexp++] = *savep++;
    }
    region->size = 0;

    if( region->index ){
        free(region->index);
        region->index = 0;
        if( memwatch ){
            printf("free %s\n", "region->index");
        }
    };
    if( region->save ){
        free(region->save);
        region->save = 0;
        if( memwatch ){
            printf("free %s\n", "region->save");
        }
    };
    memwatch = 0;
}

/* return message about geion callback */
 void RegionInfo(void)
{
    Message  message;

    sprintf(message,
            "Region: live=%d size=%d bound1=%d bound2=%d seed=(%d,%d,%d) value=%d neighborhood=%d nlist=%d",
            region->live, region->size, region->bound[0], region->bound[1],
            region->seed[DATA_AXIS1], region->seed[DATA_AXIS2],
            region->seed[DATA_AXIS3], region->seed[DATA_VALUE],
            RegionNeighbors(), region->nlist);
    UIMessage(message);
}

/* count neighbors */
int RegionNeighbors(void)
{
    return (
            (region->neighborhood & MARK_FACE) / MARK_FACE * 6 +
            (region->neighborhood & MARK_EDGE) / MARK_EDGE * 12 +
            (region->neighborhood & MARK_CORNER) / MARK_CORNER * 8);
}

/* save region parameters */
void RegionSavePar(void)
{
    Message  message;

    sprintf(message, "Region: seed=%d,%d,%d value=%d neighborhood=%d",
            region->seed[DATA_AXIS1],
            region->seed[DATA_AXIS2],
            region->seed[DATA_AXIS3],
            region->seed[DATA_VALUE], region->neighborhood);
    UISaveMessage(message);
}

/*--------------------------------------------------------------------*\
   mark contiguous voxels, given seed, min and max; continginous
   defined by mode: MARK_FACE=6, MARK_EDGE=12, MARK_CORNER=8 list
   returned; nlist is maximum size
\*--------------------------------------------------------------------*/

int RegionMark(
    Buffer data
   ,int size1
   ,int size2
   ,int size3
   ,int seed1
   ,int seed2
   ,int seed3
   ,int min
   ,int max
   ,int neighborhood
   ,Shadow list
   ,int nlist 
   ){


    Buffer cell;

    Shadow   endold;
    Shadow   listend;
    Shadow   new;
    Shadow   old;

    int      count = 0;
    int      stride1;
    int      stride2;
    int      stride3;
    int      neighbor[26];
    int      cell1;

    int*   corner1;
    int*   corner2;
    int*   edge1;
    int*   edge2;
    int*   face1;
    int*   face2;
    int*   nbr;

    int debug = 0;