package yi.iu.processor.cannyox;

import yi.iu.ui.DisplayFrame;
import java.util.Vector;

/*
#define bytePixel(buf,x,y)   (*((byte*)buf->GetElementAddr(x,y)))
#define shortPixel(buf,x,y)  (*((short*)buf->GetElementAddr(x,y)))
#define colorPixel(buf,x,y)  (0.299*((byte*)buf->GetElementAddr(x,y))[0]+0.587*((byte*)buf->GetElementAddr(x,y))[1]+0.114*((byte*)buf->GetElementAddr(x,y))[2])
#define PixelSize(buf)       (buf->GetBitsPixel())
*/

//-----------------------------------------------------------------------------
//
// -- Constructor from an image.
//

public class CannyOX{
    int[][] _image;    // The original intensity image   //BufferXY
    int _xstart,_ystart; // The origin of the buffer in the image
    int _xsize,_ysize;   // The width of the image buffer

    float[][] _smooth;     // Smoothed intensity image
    float[][] _dx;         // Derivatives in x, and sub-pixel x coordinates
    float[][] _dy;         // Derivatives in y, and sub-pixel y coordinates
    float[][] _grad;       // Gradient image

    // Quite a few of the following could be done using hash tables
    float[][] _thick;      // Gradient image after NMS
    float[][] _thin;       // Gradient image after NMS and thinning
    float[][] _theta;      // Orientation image

    int[][] _junction;     // Binary image true only at junctions ends, and relevant lists
    int[][] _jx,_jy;     // Images of (x,y) coordinates of nearest cluster centre
    Vector _xjunc,_yjunc; //int
    Vector _vlist;   // The junction cluster centres //CoolListP<Vertex*>

    float _gauss_tail;  // The value of the kernel at its tail
    float _sigma;       // Smoothing sigma
    int _width;         // The smoothing kernel width - can change
    int _w0;            // Same as above, but does not change
    int _k_size;        // The kernel is 2*_width+1s
    float[] _kernel;     // 1-Dimensional convolution kernel of size k_size
    float _low;         // Low threshold for hysteresis
    float _high;        // High threshold for hysteresis

    float _jval;        // A dummy junction intensity step value
    int _chain_no;      // A dummy variable used in following
    boolean verbose;

    int _max_width_OX;       // The maximum smoothing kernel width
    float[] _sub_area_OX;     // Used in the smoothing the image near the borders
    int _edge_min_OX;        // Minimum edge pixel intensity
    int _min_length_OX;      // Minimum number of pixels in a curve
    boolean _join_flag_OX;      // True to enable pixel jumping
    int _border_size_OX;     // Border size around the image to be set
    float _border_value_OX;  //   to _border_value_OX (usually 0) to ensure
                            //   follow won't overrun.
    float _scale_OX;         // Value used in the follow part of canny to
                            //   scale image after the hysteresis part.
    int _follow_strategy_OX; // Flag used in the Final_followOX() to determined
                            //  the order of neighboring pixel checking
                            //  (see Final_followOX() function in CannyOX.C
                            //  When equal to 0, only NMS and Hysteresis
                            //  are performed; FollowerOX is not performed.
                            //  Locating Junctions in not performed either.
                            //  See Do_CannyOX(...)
    boolean _junction_option_OX; // True if we want to locate junctions

    int[] _x; // stores the pixel x-location of the edgels
            // (not the sub pixel accuracy)
    int[] _y; // stores the pixel y-location of the edgels
            // (not the sub pixel accuracy)


    //RectROI[] imageROI;

    public CannyOX(float sigma, float low, float high, boolean v)
    {
        verbose =v;
        _sigma = sigma;
        _low = low;
        _high = high;
    }
/*
    public CannyOX(int[][] image,float sigma, int max_width, float gauss_tail,
                    float low, float high, int follow_strategy, int edge_min,
                    int min_length, int border_size, float border_value,
                    float scale, boolean join_flag, boolean junction_option, boolean v)
    {
    this(sigma,low,high,v);

    // Set the maximum allowable convolution kernel
    _gauss_tail = gauss_tail;
    _max_width_OX = max_width;
    _kernel = Make_float_vector(_max_width_OX);
    _sub_area_OX = Make_float_vector(_max_width_OX);

    _edge_min_OX = edge_min;
    _min_length_OX = min_length;
    _join_flag_OX = join_flag;
    _border_size_OX = border_size;
    _border_value_OX = border_value;
    _scale_OX = scale;
    _follow_strategy_OX = follow_strategy;
    _junction_option_OX = junction_option;

    // Get the BufferXY
    imageROI = image.GetROI();
    _xsize = min(image.GetSizeX(),imageROI.GetSizeX());
    _ysize = min(image.GetSizeY(),imageROI.GetSizeY());
    _xstart = max(0,imageROI.GetOrigX());
    _ystart = max(0,imageROI.GetOrigY());
    _image = new BufferXY(_xsize,_ysize,image.GetBitsPixel());
    image.GetSection(_image.GetBuffer(),_xstart,_ystart,_xsize,_ysize);

    _smooth = Make_float_image(_image.GetSizeX(),_image.GetSizeY());
    _dx = Make_float_image(_image.GetSizeX(),_image.GetSizeY());
    _dy = Make_float_image(_image.GetSizeX(),_image.GetSizeY());
    _grad = Make_float_image(_image.GetSizeX(),_image.GetSizeY());
    _thick = Make_float_image(_image.GetSizeX(),_image.GetSizeY());
    _thin = Make_float_image(_image.GetSizeX(),_image.GetSizeY());
    _theta = Make_float_image(_image.GetSizeX(),_image.GetSizeY());
    initFloatImage(_theta, _image.GetSizeX(),_image.GetSizeY());

    _junction = Make_int_image(_image.GetSizeX(),_image.GetSizeY());
    Zero_int_image(_junction);
    _jx = Make_int_image(_image.GetSizeX(),_image.GetSizeY());
    Zero_int_image(_jx);
    _jy = Make_int_image(_image.GetSizeX(),_image.GetSizeY());
    Zero_int_image(_jy);

    _xjunc = new Vector(); //<int>;
    _yjunc = new Vector(); //<int>;
    _vlist = new Vector(); //P<Vertex*>;

    _jval = 2000;
    }

*/

    //-----------------------------------------------------------------------------
    //
    // -- Constructor from a BufferXY.
    //
    public CannyOX(int[][] image,float sigma, int max_width, float gauss_tail,
                    float low, float high, int follow_strategy, int edge_min,
                    int min_length, int border_size, float border_value,
                    float scale, boolean join_flag, boolean junction_option, boolean v)

    {
        this(sigma,low,high,v);
        // Set the maximum allowable convolution kernel
        _gauss_tail = gauss_tail;
        _max_width_OX = max_width;
        _kernel = Make_float_vector(_max_width_OX);
        _sub_area_OX = Make_float_vector(_max_width_OX);

        _edge_min_OX = edge_min;
        _min_length_OX = min_length;
        _join_flag_OX = join_flag;
        _border_size_OX = border_size;
        _border_value_OX = border_value;
        _scale_OX = scale;
        _follow_strategy_OX = follow_strategy;
        _junction_option_OX = junction_option;

        // Get the BufferXY
        //imageROI = 0;
        //_xsize  = image.GetSizeX();
        //_ysize = image.GetSizeY();
        _xsize  = image[0].length;
        _ysize = image.length;
        _xstart = 0;
        _ystart = 0;
        /*
        _image = new  BufferXY(image.GetSizeX(),
                                image.GetSizeY(),
                                image.GetBitsPixel());
        _image.WriteBytes(image.GetBuffer(), image.GetSize());
        */
            _image= (int[][]) image.clone();

        _smooth = Make_float_image(_xsize,_ysize);
        _dx = Make_float_image(_xsize,_ysize);
        _dy = Make_float_image(_xsize,_ysize);
        _grad = Make_float_image(_xsize,_ysize);
        _thick = Make_float_image(_xsize,_ysize);
        _thin = Make_float_image(_xsize,_ysize);
        _theta = Make_float_image(_xsize,_ysize);
        initFloatImage(_theta, _xsize,_ysize);

        _junction = Make_int_image(_xsize,_ysize);
        Zero_int_image(_junction);
        _jx = Make_int_image(_xsize,_ysize);
        Zero_int_image(_jx);
        _jy = Make_int_image(_xsize,_ysize);
        Zero_int_image(_jy);

        _xjunc = new Vector(); //<int>;
        _yjunc = new Vector(); //<int>;
        _vlist = new Vector(); //P<Vertex*>;

        _jval = 2000;
    }




    //-----------------------------------------------------------------------------
    //
    // -- Method to do canny on a CannyOX class that has all the members defined.
    //
    public Vector Do_cannyOX()
    {
        Vector edges=null;

        if (verbose) {
            System.err.println("Doing Canny on image region ");
            System.err.println( _xsize + " by " + _ysize);
            System.err.println( "Sigma               = " + _sigma);
            System.err.println( "Gaussian tail       = " + _gauss_tail);
            System.err.println( "Max kernel size     = " + _max_width_OX);
                System.err.println( "Upper threshold     = " + _high);
                System.err.println( "Lower threshold     = " + _low);
                System.err.println( "Min edgel intensity = " + _edge_min_OX);
                System.err.println( "Min edge length     = " + _min_length_OX);
                System.err.println( "Image border size   = " + _border_size_OX);
        }

        // Do the traditional Canny parts
        if (verbose)  System.err.println( "setting convolution kernel and zeroing images\n");
        Set_kernelOX();

        if (verbose)  System.err.println( "Kernel size     = " + _k_size );
        if (verbose)  System.err.println( "smoothing the image\n");
        Smooth_imageOX();
        //display_float_image("smooth",_smooth);

        if (verbose)  System.err.println( "computing x derivatives, y derivatives and norm of gradient\n");
        Compute_gradientOX();
        //display_float_image("dx",_dx);
        //display_float_image("dy",_dy);
        //display_float_image("grad",_grad);

        if (verbose)  System.err.println( "doing non-maximal supression\n");
        int n_edgels_NMS = Non_maximal_supressionOX();
        //if(verbose) System.err.println( "Number of edgels after NMS = " << n_edgels_NMS << endl);

        // (_x,_y) holds the pixel location and not the sub pixel accuracy
        // needed in Get_hysteresis_edgelsOX() to fill the _thin image
        // at the edgels locations
        // (_x,_y)'s are initialised in Get_NMS_edgelsOX()
        _x = new int[n_edgels_NMS];
        _y = new int[n_edgels_NMS];

        // Copy edgels from _thick image into and EdgelChain class.
        // Edgels after NMS are not connected, but still are stored in EdgelChain.
        // edgels_NMS is the input to Hysteresis
        EdgelChain edgels_NMS = new EdgelChain(n_edgels_NMS);
        Get_NMS_edgelsOX(edgels_NMS,n_edgels_NMS);

        //if(verbose) cerr << "doing hysteresis\n";
        int[] status = new int[n_edgels_NMS];
        int n_edgels_Hysteresis = HysteresisOX(edgels_NMS, status);
        //if(verbose) cerr << "Number of edgels after Hysteresis = " << n_edgels_Hysteresis << endl;
        EdgelChain edgels_Hysteresis = new EdgelChain(n_edgels_Hysteresis);
        Get_hysteresis_edgelsOX(edgels_NMS,status,edgels_Hysteresis);

        // delete the edgels that are output of Non_maximal_supression
        //delete [] status;
        status=null;
        //delete edgels_NMS;
        edgels_NMS=null;
        // delete *_x and *_y; they are not needed anymore
        //delete [] _x;
        _x=null;
        //delete [] _y;
        _y=null;

        if(_follow_strategy_OX == 0) // Don't do the follow stage of canny
            {
            if(edges==null){
                System.err.println("Do_cannyOX: First create an empty edge list!\n");
                return null;
            }
            NO_FollowerOX(edges, edgels_Hysteresis);
            // delete the edgels that are output of Hysteresis
            //delete edgels_Hysteresis;
            edgels_Hysteresis=null;
            return null;
            }

        // delete the edgels that are output of Hysteresis
        //delete edgels_Hysteresis;
        edgels_Hysteresis=null;

        // Multiply _thin image by _scaleOX. Values that are above 255 are set to 255
        Scale_imageOX(_thin, _scale_OX);
        // Set image border to _border_valueOX (default = 0) so follow can't overrun
        Set_image_borderOX(_thin, _border_size_OX, _border_value_OX);

        if (_junction_option_OX)
            {
            // Locate junctions in the edge image
            //if(verbose) cerr << "locating junctions in the edge image - ";
            Find_junctionsOX();
            //if(verbose) cerr << _xjunc->length() << " junctions found\n";
            Find_junction_clustersOX();
            //if(verbose)  cerr << _vlist->length() << " junction clusters found\n";
            //int[] la= yi.util.ImageUtil.twoD2oneD(_junction);
            //int[] img= yi.util.ImageUtil.makeNewImageFromLabel(la);
            //(new DisplayFrame("junction", img, _junction.length, _junction[0].length)).show();

            }

        // Finally do edge following to extract the edge data from the _thin image
        //if(verbose) cerr << "doing final edge following\n";
        if(edges!=null) {
            System.err.println("Do_cannyOX: First create an empty edge list!\n");
            return null;
        }
        edges= new Vector();
        FollowerOX(edges);

        //garbage collection
        _smooth=_dx=_dy=_grad=_thick=_thin=_theta=null;
        _junction=_jx=_jy=null;
        _xjunc=_yjunc=_vlist=null;
        _x=_y=null;
        System.gc();

        return edges;

        // Finally: check for mask in ROI

        //AdjustForMask(edges);


        //if(verbose) cerr << "finished CannyOX\n";
    }



    //-----------------------------------------------------------------------------
    //
    // -- Sets up the convolution kernel based on the first derivative of a
    // Gaussian.  Protected.
    //
    void Set_kernelOX()
    {
        int i;
        float area = 0;
        float s2 = (float)2.0*_sigma*_sigma;

        for (i=0;i<_max_width_OX;i++)
        {
            _kernel[i] = (float)0.0;
            for (int j=-5; j<=5; j++)
            {
                    float x = (float)i + (float)0.1*j;
                    _kernel[i] += Math.exp(-x*x/s2);
                }   // end for j
            _kernel[i] /= 11.0;                    // Averaging over six points
            if (_kernel[i] < _gauss_tail && i>0)
                break;
            area += _kernel[i];                   // Area under half profile curve
        }   // end for i

        float total_area = (float)2.0*area - _kernel[0]; // Total area under profile curve
        _width = i;                               // half Kernel width
        _k_size = 2*_width - 1;                   // kernel size

        for (i=0;i<_width;i++)
        {
            _sub_area_OX[i] = (total_area - area)/total_area;
            area -= _kernel[i];
            _kernel[i] /= total_area;
        }   // end for i
    }



    //-----------------------------------------------------------------------------
    //
    // -- Convolves the image with the smoothing kernel.
    //  _sub_area_OX is used to smooth pixels lying in a border of size _width_OX.
    //  Protected.
    //
    //  Meaning of _sub_area_OX[x] :
    //    It is the probability that the variable will assume a value <= x
    //    in the distribution function. It is the hashed area under the
    //    the distribution profile.
    // @{
    // \begin{verbatim}
    //                     . .
    //                   . ### .
    //                 .  #####| .
    //               . ########|   .
    //             . ##########|     .
    //         __ .############|      .__
    //       __________________|___________
    //                         x
    // \end{verbatim}
    // @}
    //
    void Smooth_imageOX()
    {
        int x,y,k;

        float[][] tmp = Make_float_image(_xsize,_ysize);
        //Zero_float_image(tmp);

        if (PixelSize(_image) == 8)
        {
            // x direction
            for (y=0;y < _ysize;y++)
                {
                    // left border of size _width
                    for (x=0;x < _width-1;x++)
                    {
                        tmp[x][y] = bytePixel(_image,x,y)*_kernel[0];
                        for (k=1; k < _width; k++)
                            {
                            if (x-k < 0)
                                tmp[x][y] += bytePixel(_image,x+k,y)*_kernel[k];
                            else
                                tmp[x][y] += bytePixel(_image,x-k,y)*_kernel[k] + bytePixel(_image,x+k,y)*_kernel[k];
                            }   // end for k
                        tmp[x][y] /= _sub_area_OX[x+1];
                    }  // end for x

                // Middle pixels along x direction
                for (x=_width-1;x < _xsize-_width+1;x++)
                    {
                    tmp[x][y] = bytePixel(_image,x,y)*_kernel[0];
                    for (k=1; k < _width; k++)
                        {
                            tmp[x][y] += bytePixel(_image,x-k,y)*_kernel[k] + bytePixel(_image,x+k,y)*_kernel[k];
                        }  // end for k
                    }   // end for x
                // Right border of size _width
                for (x=_xsize-_width+1;x < _xsize;x++)
                    {
                    tmp[x][y] = bytePixel(_image,x,y)*_kernel[0];
                    for (k=1; k < _width; k++)
                        {
                        if (x+k >= _xsize)
                            tmp[x][y] += bytePixel(_image,x-k,y)*_kernel[k];
                        else
                            tmp[x][y] += bytePixel(_image,x-k,y)*_kernel[k] + bytePixel(_image,x+k,y)*_kernel[k];
                        }   // end for k
                    tmp[x][y] /= _sub_area_OX[_xsize-x];
                    }   // end for x
                }   // end for y
            }   // end if

        else if (PixelSize(_image) == 16)
            {
            // x direction
            for (y=0;y < _ysize;y++)
                {
                // left border of size _width
                for (x=0;x < _width-1;x++)
                    {
                    tmp[x][y] = shortPixel(_image,x,y)*_kernel[0];
                    for (k=1; k < _width; k++)
                        {
                        if (x-k < 0)
                            tmp[x][y] += shortPixel(_image,x+k,y)*_kernel[k];
                        else
                            tmp[x][y] += shortPixel(_image,x-k,y)*_kernel[k] + shortPixel(_image,x+k,y)*_kernel[k];
                        }   // end for k
                    tmp[x][y] /= _sub_area_OX[x+1];
                    }  // end for x
                // Middle pixels along x direction
                for (x=_width-1;x < _xsize-_width+1;x++)
                    {
                    tmp[x][y] = shortPixel(_image,x,y)*_kernel[0];
                    for (k=1; k < _width; k++)
                        {
                            tmp[x][y] += shortPixel(_image,x-k,y)*_kernel[k] + shortPixel(_image,x+k,y)*_kernel[k];
                        }  // end for k
                    }   // end for x
                // Right border of size _width
                for (x=_xsize-_width+1;x < _xsize;x++)
                    {
                    tmp[x][y] = shortPixel(_image,x,y)*_kernel[0];
                    for (k=1; k < _width; k++)
                        {
                        if (x+k >= _xsize)
                            tmp[x][y] += shortPixel(_image,x-k,y)*_kernel[k];
                        else
                            tmp[x][y] += shortPixel(_image,x-k,y)*_kernel[k] + shortPixel(_image,x+k,y)*_kernel[k];
                        }   // end for k
                    tmp[x][y] /= _sub_area_OX[_xsize-x];
                    }   // end for x
                }   // end for y
            }   // end else if
        else if (PixelSize(_image) == 24) // assuming this is a color image
            {
            // x direction
            for (y=0;y < _ysize;y++)
                {
                // left border of size _width
                for (x=0;x < _width-1;x++)
                    {
                    tmp[x][y] = colorPixel(_image,x,y)*_kernel[0]; // average color value
                    for (k=1; k < _width; k++)
                        {
                        if (x-k < 0)
                            tmp[x][y] += colorPixel(_image,x+k,y)*_kernel[k];
                        else
                            tmp[x][y] += colorPixel(_image,x-k,y)*_kernel[k] + colorPixel(_image,x+k,y)*_kernel[k];
                        }   // end for k
                    tmp[x][y] /= _sub_area_OX[x+1];
                    }  // end for x
                // Middle pixels along x direction
                for (x=_width-1;x < _xsize-_width+1;x++)
                    {
                    tmp[x][y] = colorPixel(_image,x,y)*_kernel[0];
                    for (k=1; k < _width; k++)
                        {
                            tmp[x][y] += colorPixel(_image,x-k,y)*_kernel[k] + colorPixel(_image,x+k,y)*_kernel[k];
                        }  // end for k
                    }   // end for x
                // Right border of size _width
                for (x=_xsize-_width+1;x < _xsize;x++)
                    {
                    tmp[x][y] = colorPixel(_image,x,y)*_kernel[0];
                    for (k=1; k < _width; k++)
                        {
                        if (x+k >= _xsize)
                            tmp[x][y] += colorPixel(_image,x-k,y)*_kernel[k];
                        else
                            tmp[x][y] += colorPixel(_image,x-k,y)*_kernel[k] + colorPixel(_image,x+k,y)*_kernel[k];
                        }   // end for k
                    tmp[x][y] /= _sub_area_OX[_xsize-x];
                    }   // end for x
                }   // end for y
            }   // end if

        else
            System.err.println("Image Depth must be either 8, 16 or 24 bits. ");

        // y direction
        for (x=0;x < _xsize;x++)
            {
            // Top border of size _width
            for (y=0;y < _width-1;y++)
                {
                _smooth[x][y] = tmp[x][y]*_kernel[0];
                for (k=1; k < _width; k++)
                    {
                    if (y-k < 0)
                        _smooth[x][y] += tmp[x][y+k]*_kernel[k];
                    else
                        _smooth[x][y] += tmp[x][y-k]*_kernel[k] + tmp[x][y+k]*_kernel[k];
                    }   // end for k
                _smooth[x][y] /= _sub_area_OX[y+1];
                }  // end for y
            // Middle pixels along y direction
            for (y=_width-1;y < _ysize-_width+1;y++)
                {
                _smooth[x][y] = tmp[x][y]*_kernel[0];
                for (k=1; k < _width; k++)
                    {
                    _smooth[x][y] += tmp[x][y-k]*_kernel[k] + tmp[x][y+k]*_kernel[k];
                    }  // end for k
                }   // end for y
            // Bottom border of size _width
            for (y=_ysize-_width+1;y < _ysize;y++)
                {
                _smooth[x][y] = tmp[x][y]*_kernel[0];
                for (k=1; k < _width; k++)
                    {
                    if (y+k >= _ysize)
                        _smooth[x][y] += tmp[x][y-k]*_kernel[k];
                    else
                        _smooth[x][y] += tmp[x][y-k]*_kernel[k] + tmp[x][y+k]*_kernel[k];
                    }   // end for k
                _smooth[x][y] /= _sub_area_OX[_ysize-y];
                }   // end for y
            }   // end for x

        tmp=null;
    }



    //-----------------------------------------------------------------------------
    //
    //  Computes the absolute intensity surface gradient, _grad[][].  Protected.
    //
    void Compute_gradientOX()
    {
        int x,y;

        //  Topline
        _dy[0][0]=(_smooth[0][0]-_smooth[0][1])*2;   // left edge
        _dx[0][0]=(_smooth[1][0]-_smooth[0][0])*2;
        for (x=1;x<_xsize-1;x++)
            {
            _dy[x][0]=(_smooth[x][0]-_smooth[x][1])*2;
            _dx[x][0]=_smooth[x+1][0]-_smooth[x-1][0];
            }   // end for x
        _dy[_xsize-1][0]=(_smooth[_xsize-1][0]-_smooth[_xsize-1][1])*2;   // right edge
        _dx[_xsize-1][0]=(_smooth[_xsize-1][0]-_smooth[_xsize-2][0])*2;

        // Middle lines
        for (y=1;y<_ysize-1;y++)
            {
            _dy[0][y]=_smooth[0][y-1]-_smooth[0][y+1];   // left edge
            _dx[0][y]=(_smooth[1][y]-_smooth[0][y])*2;
            for (x=1;x<_xsize-1;x++)
                {
                _dy[x][y]=_smooth[x][y-1]-_smooth[x][y+1];
                _dx[x][y]=_smooth[x+1][y]-_smooth[x-1][y];
                }   // end for x
            _dy[_xsize-1][y]=_smooth[_xsize-1][y-1]-_smooth[_xsize-1][y+1];   // right edge
            _dx[_xsize-1][y]=(_smooth[_xsize-1][y]-_smooth[_xsize-2][y])*2;
            }   // end for y


        // Bottom line
        _dy[0][_ysize-1]=(_smooth[0][_ysize-2]-_smooth[0][_ysize-1])*2;   // left edge
        _dx[0][_ysize-1]=(_smooth[1][_ysize-1]-_smooth[0][_ysize-1])*2;
        for (x=1;x<_xsize-1;x++)
            {
            _dy[x][_ysize-1]=(_smooth[x][_ysize-2]-_smooth[x][_ysize-1])*2;
            _dx[x][_ysize-1]=_smooth[x+1][_ysize-1]-_smooth[x-1][_ysize-1];
            }   // end for x
        _dy[_xsize-1][_ysize-1]=(_smooth[_xsize-1][_ysize-2]-_smooth[_xsize-1][_ysize-1])*2;   // right edge
        _dx[_xsize-1][_ysize-1]=(_smooth[_xsize-1][_ysize-1]-_smooth[_xsize-2][_ysize-1])*2;

        //  Magnitude for entire image
        for (y=0;y<_ysize;y++)
            for (x=0;x<_xsize;x++)
                _grad[x][y] = (float)Math.sqrt(_dx[x][y]*_dx[x][y] + _dy[x][y]*_dy[x][y]);
    }



//-----------------------------------------------------------------------------
//
// -- Non-maximal supression.  Protected.
//
int Non_maximal_supressionOX()
{
  int x,y;
  float k = (float)(180.0/Math.PI);
  float newx, newy, theta, thick;
  int n_edgels_NMS = 0;

  for (y=_ysize-2;y>0;y--)
    {
      for (x=_xsize-2;x>0;x--)
        {
          float del;
          if (Math.abs(_dx[x][y])>Math.abs(_dy[x][y]))
            {
              if (_grad[x][y]<=_grad[x+1][y] || _grad[x][y]<_grad[x-1][y])
                continue;
	    }  	    }	// end if
          else if (_grad[x][y]<=_grad[x][y-1] || _grad[x][y]<_grad[x][y+1])
            continue;

          // we have an edge
          thick = _grad[x][y];
          theta = (float)(k*Math.atan2(_dx[x][y],_dy[x][y])); // Not to be used to define
                                                //  _theta[x][y].
                                                //  Only to define orient.
          int orient = ( (int) (theta+202.5) ) / 45;

          newx = 0; newy = 0; // Initialise

          // Identify quadrant:
          //                     3   2   1
          //
          //                     4   *   0(8)
          //
          //                     5   6   7
          switch(orient)
            {
            case 0:   // sort of horizontal
            case 4:
            case 8:
              newx=(x+(float)0.5);   // pixel centre
              del=_grad[x][y-1]-_grad[x][y+1];
              del/=(_grad[x][y+1]+_grad[x][y-1]-2*_grad[x][y])*2;
              if (del>0.5) continue;
              newy=(y+del+(float)0.5);
              break;
            case 2:   // sort of vertical
            case 6:
              newy=y+(float)0.5;
              del=_grad[x-1][y]-_grad[x+1][y];
              del/=(_grad[x-1][y]+_grad[x+1][y]-2*_grad[x][y])*2;
              if (del>0.5) continue;
              newx=x+del+(float)0.5;
              break;
            case 1:   // sort of left diagonal
            case 5:
              if (_grad[x][y]<=_grad[x+1][y-1] || _grad[x][y]<_grad[x-1][y+1])
                continue;
              del=_grad[x-1][y+1]-_grad[x+1][y-1];
              del/=(_grad[x-1][y+1]+_grad[x+1][y-1]-2*_grad[x][y])*2;
              if (del>0.5) continue;
              newy=y-del+(float)0.5;
              newx=x+del+(float)0.5;
              break;
            case 3:   // sort of right diagonal
            case 7:
              if (_grad[x][y]<=_grad[x-1][y-1] || _grad[x][y]<_grad[x+1][y+1])
                continue;
              del=_grad[x+1][y+1]-_grad[x-1][y-1];
              del/=(_grad[x+1][y+1]+_grad[x-1][y-1]-2*_grad[x][y])*2;
              if (del>0.5) continue;
              newy=y-del+(float)0.5;
              newx=x-del+(float)0.5;
              break;
            }   // end switch

          // _theta[x][y] as defined in the next line is compatible with
          //  the convention in TargetJr.
          //  The minus sign in front of _dy[x][y] is to change the way
          //  _dy is defined in this CannyOX (i.e, we want it to be
          //  [y(i+1) - y(i-1)] rather than [y(i-11) - y(i+1)]).
          //   See ComputeGradient above.
          //  _theta[x][y] now stores the normal to the edge tangent.
          //  Before it stored the tangent to the edge.
          //  _theta[x][y] = theta;  // This how it was defined previously
          _theta[x][y] = (float)(k*Math.atan2(-_dy[x][y],_dx[x][y]));

          _thick[x][y] = thick;
          _dx[x][y] = newx;
          _dy[x][y] = newy;

          n_edgels_NMS++;
        }   // end for x
    }   // end for y
  return n_edgels_NMS;
}



//-----------------------------------------------------------------------------
//
// -- Copy edgels from _thick image. These edgels are the result of NMS.
// Although the edgels are stored in EdgelChain edgels_NMS, they are
// not actually connected.
// Also initialises _x and _y: the pixel locations of the edgels. These
// are needed later in filling the image _thin with the edgels after
// hysteresis.  Protected.
//
void Get_NMS_edgelsOX(EdgelChain edgels_NMS,
                               int n_edgels_NMS)
{
  int x, y, i = 0;

  edgels_NMS.SetLength(n_edgels_NMS);

  for (y=1;y<=_ysize-2;y++)
    for (x=1;x<=_xsize-2;x++)
          if ( _thick[x][y] != 0 && i < n_edgels_NMS)
            {
              // fill edgels_NMS
              edgels_NMS.SetX(_dx[x][y],i);
              edgels_NMS.SetY(_dy[x][y],i);
              edgels_NMS.SetGrad(_thick[x][y],i);
              edgels_NMS.SetTheta(_theta[x][y],i);
              // remember the pixel locations (x,y) of the edgels
              _x[i] = x; _y[i] = y;
              i++;
            }
}



//-----------------------------------------------------------------------------
//
// -- Hysteresis follows edgels that lie above the low threshold and have at
// least one edgel above the high threshold.  Protected.
//
int HysteresisOX(EdgelChain edgels_NMS,
                          int[] status)
{
  float low, high;
  int[] rows;
  char[] row;
  char[] col;
  LINK[] links;

  int n_edgels_NMS = edgels_NMS.length();

  // Compute lower and upper threshold
  low  =(float) ((32.0/Math.log(2.0)) * Math.log(_low/100+1.0));     // lower threshold
  high = (float)((32.0/Math.log(2.0)) * Math.log(_high/100+1.0));    // upper threshold

  rows = Make_u_int_vectorOX(_ysize+1);
  row = Make_u_short_vectorOX(n_edgels_NMS);
  col = Make_u_short_vectorOX(n_edgels_NMS);

  links = new LINK[n_edgels_NMS];

    //prepare all the data structures
    Init_hyster_parametersOX(rows,status,links,row,col,edgels_NMS);

  Link_edgelsOX(col,rows,links, n_edgels_NMS); // for BOUNDS (ABR purify error)
  Do_hysteresisOX(edgels_NMS,links,status,low,high);
  int n_edgels_Hysteresis = Get_n_edgels_hysteresisOX(edgels_NMS,status);

  // Freeing memory that was allocated in Init_hyster_parametersOX()
  rows=null;
  row=null;
  col=null;
  links=null;

  return n_edgels_Hysteresis;
}



//-------------------------------------------------------------------------
//
// -- Perform Hysteresis part of canny.  Protected.
//
void Do_hysteresisOX(EdgelChain edgels_NMS,
                              LINK[] links,
                              int[] status,
                              float low,
                              float high)
{
  int i;
  LINK lptr;

  int n_edgels_NMS = edgels_NMS.length();

  for (i=0;i<n_edgels_NMS;++i)
    if (!(status[i]!=0) && edgels_NMS.GetGrad(i)>high)
      {
        status[i]=1;
        for (lptr=links[i];lptr!=null;lptr=lptr.nextl)
          Initial_followOX(lptr.to,i,edgels_NMS,links,status,low);
      } // end if
}



//-------------------------------------------------------------------------
//
// -- Initial follow. Used for the hysteresis part of canny.  Protected.
//
void Initial_followOX(int to,
                       int from,
                       EdgelChain edgels_NMS,
                       LINK[] links,
                       int[] status,
                       float low)
{
  LINK lptr;

  if (!(status[to]!=0) && edgels_NMS.GetGrad(to) > low)
    status[to]=1;
  else
    return;

  for (lptr=links[to];lptr!=null;lptr=lptr.nextl)
    if (lptr.to!=from)
      Initial_followOX(lptr.to,to,edgels_NMS,links,status,low);
}



    //-------------------------------------------------------------------------
    //
    // -- Initialise the parameter needed to do the hysteresis part of canny.
    // Memory allocated in this function is freed in the function
    // HysteresisOX(), from where this function is called.  Protected.
    //
    void Init_hyster_parametersOX(int[] rows,
                                        int[] status,
                                        LINK[] links,
                                        char[] row,
                                        char[] col,
                                        EdgelChain edgels_NMS)
    {
        int i,j;

        int n_edgels_NMS = edgels_NMS.length();
        /*
        // allocate memory
        rows = Make_u_int_vectorOX(_ysize+1);
        row = Make_u_short_vectorOX(n_edgels_NMS);
        col = Make_u_short_vectorOX(n_edgels_NMS);

        (*links) = new LINK*[n_edgels_NMS];
        */
        for (i=0;i<n_edgels_NMS;++i)
        {
            links[i]  = null;
            row[i]    = (char) edgels_NMS.GetY(i);
            col[i]    = (char) edgels_NMS.GetX(i);
            status[i] = 0;
        } // end for i

        for (i=0,j=0;i<=_ysize;++i)
        {
            // (j<n_edgels_NMS) && (row[j]<i)
            while (j<n_edgels_NMS ? row[j]<i : false)
            {
                    ++j;
                } // end while
                rows[i]=j;   // the first edgel after start of row i
            } // end for
    }



    //-------------------------------------------------------------------------
    //
    // -- Add link.  Protected.
    //
    void Add_linkOX(int edgel,
                            int to,
                            LINK[] links)
    {
        LINK lptr1=new LINK();
        LINK lptr2=new LINK();

        lptr1.to=to;
        lptr1.nextl=links[edgel];
        links[edgel]=lptr1;

        lptr2.to=edgel;
        lptr2.nextl=links[to];
        links[to]=lptr2;
    }



    //-------------------------------------------------------------------------
    //
    // -- Link edgels.  Protected.
    //
    void Link_edgelsOX(char c[],
                                int[] rows,
                                LINK[] links,
                                int csize) // for BOUNDS (ABR purify error)
    {
        char col;
        char e,w,s;

        for (int i=0;i<_ysize-1;++i)
        {
            for (int j=rows[i];j<(int)(rows[i+1]);++j)
            {
                col=c[j];
                e=0;
                w=0;
                s=0;

                if (j+1<(int)(rows[i+1]) && col+1==c[j+1])
                {
                    Add_linkOX(j,j+1,links);
                    e=1;
                }	// end if

                    if (j>(int)(rows[i]) && j>0 && col==c[j-1]+1) // BOUNDS (ABR purify error)
                        w=1;

                    int k;
                    for (k=rows[i+1]; k+1<(int)(rows[i+2]); ++k)
                        if (c[k]>col)
                        break;

                    if( csize > 0 && k >= csize) {
                        //cerr << "Warning: attempt to read beyond c["<<csize-1<<"] (k="<<k<<", c["<<csize-1<<"]="<<c[csize-1]<<", col="<<col<<")\n";
                        k = csize-1; // avoid ABR purify error below
                    }

                    if (col==c[k])       // was last on the row // BOUNDS (ABR purify error)
                        {
                        Add_linkOX(j,k,links);
                        s=1;
	                }           }	// end if

                    if (k>(int)(rows[i+1]) && k>0 && col==c[k-1])
                        {
                        Add_linkOX(j,k-1,links);
                        s=1;
	                }           }	// end if

                    if (e==0 && s==0 && col+1==c[k]) // BOUNDS
                        Add_linkOX(j,k,links);

                    if (w==0 && s==0 && col==c[k]+1) // BOUNDS
                        Add_linkOX(j,k,links);

                    if (w==0 && s==0 && k>(int)(rows[i+1]) && k>0 && col==c[k-1]+1)
                        Add_linkOX(j,k-1,links);
            } // end for j
        } // end for i
    }



//-------------------------------------------------------------------------
//
// -- Returns the number of edgels after hysteresis.  Protected.
//
int Get_n_edgels_hysteresisOX(EdgelChain edgels_NMS,
                                       int[] status)
{
  int i;			// loop control variables

  int n_edgels_Hysteresis = 0;
  int n_edgels_NMS = edgels_NMS.length();

  for (i=0;i<n_edgels_NMS;i++)
    {
      if (status[i]!=0)
        {
          n_edgels_Hysteresis++;
        } // end if
    } // end for i
  return n_edgels_Hysteresis;
}



//-------------------------------------------------------------------------
//
// -- Returns the edgels after hysteresis
// Also fill the image _thin with the edgels after hysteresis for
// further processing by the FollowerOX() part.  Protected.
//
void Get_hysteresis_edgelsOX(EdgelChain edgels_NMS,
                                      int[] status,
                                      EdgelChain edgels_Hysteresis)
{
  int i,j;			// loop control variables

  int n_edgels_NMS = edgels_NMS.length();

  // Initialise _thin to zero's
  Zero_float_image(_thin);

  for (i=0,j=0;i<n_edgels_NMS;i++)
    {
      if (status[i]!=0)
        {
          // Fill edgels_Hysteresis
          edgels_Hysteresis.SetX(edgels_NMS.GetX(i),j);
          edgels_Hysteresis.SetY(edgels_NMS.GetY(i),j);
          edgels_Hysteresis.SetGrad(edgels_NMS.GetGrad(i),j);
          edgels_Hysteresis.SetTheta(edgels_NMS.GetTheta(i),j);
          // Re-fill image _thin for further processing
          _thin[_x[i]][_y[i]] = edgels_NMS.GetGrad(i);
          j++;
        } // end if
    } // end for i
}



//-----------------------------------------------------------------------------
//
// -- Returns a vector of char length m.  Protected.
//
char[] Make_u_short_vectorOX(int m)
{
    return new char[m];
}

int[] Make_u_int_vectorOX(int m)
{
    return new int[m];
}



//-----------------------------------------------------------------------------
//
// -- Delete **links
// Delete an array of LINK *.  Protected.
//
void Delete_linksOX(LINK[] links, int n)
{
  int i;
  LINK tmp_link1 = null, tmp_link2 = null;

  for (i=0; i<n; i++) {
    tmp_link1 = links[i];
    if (tmp_link1 != null)
      tmp_link2 = tmp_link1.nextl;
    while (tmp_link1 != null) {
      tmp_link1 = tmp_link2;
      if (tmp_link1 != null)
        tmp_link2 = tmp_link1.nextl;
    }
  }
  //delete[] links; // memory leak removed ;-) PVR-24sept97
}



//----------------------------------------------------------------------------
//
// -- NO_FollowerOX : In the case of _follow_strategy_OX = 0, this function
// fills the (CoolListP<Edge*> *edges from EdgelChain *edgels_Hysteresis.
// i.e., the result of the Hysteresis part of Canny is written to
// (CoolListP<Edge*> *edges. The Follow part (FollowerOX) is not executed.
//  Protected.
//
// NOTE: The edgels from Hysteresis are stored in a CoolListP<Edge*> *edges.
//       Actually this List contains one Edge only that consists of a
//       DigitalCurve and two vertices.
//
void NO_FollowerOX(Vector edges, EdgelChain edgels_Hysteresis)
{
  Edge edge;
  DigitalCurve dc;
  int i;
  Vertex V1, V2;
  float[] px,py,pg,pt;
  int n_edgels_Hysteresis = edgels_Hysteresis.length();

  edges.removeAllElements();

  // V1 and V2 do not have any meaning in this case (dummies).
  // Just initialise them to zero's
  V1 = new Vertex((float)0.0,(float)0.0);
  V2 = V1;
  edge = new Edge(V1,V2);

  // Define a digital curve to store the edgels data
  //  dc only stores the list of edgels after the Hysteresis stage of canny.
  //  dc is used later in this function to define an Edge *edge
  dc = new DigitalCurve(n_edgels_Hysteresis);
  px = dc.GetX();    py = dc.GetY();
  pg = dc.GetGrad(); pt = dc.GetTheta();
int _px=0,_py=0,_pg=0,_pt=0;

  dc.SetStartX(edgels_Hysteresis.GetX(0)+_xstart);
  dc.SetStartY(edgels_Hysteresis.GetY(0)+_ystart);
  dc.SetEndX(edgels_Hysteresis.GetX(n_edgels_Hysteresis-1)+_xstart);
  dc.SetEndY(edgels_Hysteresis.GetY(n_edgels_Hysteresis-1)+_ystart);

  for(i=0;i<n_edgels_Hysteresis;i++)
    {
      px[_px++] = edgels_Hysteresis.GetX(i) + _xstart;
      py[_py++] = edgels_Hysteresis.GetY(i) + _ystart;
      pg[_pg++] = edgels_Hysteresis.GetGrad(i);
      pt[_pt++] = edgels_Hysteresis.GetTheta(i);
    }

  // Store the edgels data in an edge
  edge.SetCurve(dc);
  // Add edge to edges. Actually edges will have only one Edge (a list of all
  //  edgels after the Hyteresis stage of Canny.
  //  edges will be returned to the function where Do_CannyOX is called from.
  edges.insertElementAt(edge,0);
}



//----------------------------------------------------------------------------
//
// -- Go through every point in the image and for every one which is above a
// threshold:   follow from the point, in one direction and then the other.
// Protected.
//
void FollowerOX(Vector edges)
{
  Vector xcoords= new Vector(),ycoords= new Vector();   //CoolList<int>
  Vector grad= new Vector();              //CoolList<float>
  int x,y,tmpx=0,tmpy=0;
  float[] px,py,pg,pt;
  float val;
  DigitalCurve dc;
  Vertex v1,v2,V1,V2;
  Edge edge;

  if(_junction_option_OX)
    _chain_no = 10;    // Must be set to a number >= 1

  edges.removeAllElements();

  for(x=_border_size_OX;x<_xsize-_border_size_OX;x++)
    {
      for (y = _border_size_OX; y < _ysize-_border_size_OX; y++)
        {
          // Due to Initial_hysteresis we can follow everything > _edge_min_OX
          if( (_thin[x][y] >= _edge_min_OX) && (_junction[x][y]==0) )
            {
              if(_junction_option_OX)
                _chain_no++;

              Final_followOX(x,y,xcoords,ycoords,grad,0);

              // We may have picked up the edgel chain somewhere
              // away from its ends. Therefore, reverse the list
              // and try to follow again.
              //xcoords.reverse();  ycoords.reverse();  grad.reverse();
              reverse(xcoords);  reverse(ycoords);  reverse(grad);
              Final_followOX(x,y,xcoords,ycoords,grad,1);

              if (xcoords.size() < _min_length_OX)
                {
                  // cerr << "short list found in Final_followOX\n";
                  xcoords.removeAllElements();  ycoords.removeAllElements();  grad.removeAllElements();
                  continue;
                }

              // Create an EdgeChain and add to the list
              int count=xcoords.size();
              dc = new DigitalCurve(count);
              px = dc.GetX();     py = dc.GetY();
              pg = dc.GetGrad();  pt = dc.GetTheta();
              int _px=0,_py=0,_pg=0,_pt=0;


              // Write the points to the EdgelChain and the end points to the Curve
              dc.SetStartX(((Integer)xcoords.elementAt(0)).intValue()+_xstart);
              dc.SetStartY(((Integer)ycoords.elementAt(0)).intValue()+_ystart);
              while (count!=0)
                {
                  tmpx = ((Integer)xcoords.firstElement()).intValue();
                  tmpy = ((Integer)ycoords.firstElement()).intValue();
                  val = ((Float)grad.firstElement()).floatValue(); --count;
                  xcoords.removeElementAt(0);ycoords.removeElementAt(0);grad.removeElementAt(0);
                  // If we are not at a junction use sub-pixel value
                  if( val != _jval )
                    {
                      px[_px++] = _dx[tmpx][tmpy] + _xstart;
                      py[_py++] = _dy[tmpx][tmpy] + _ystart;
                      pg[_pg++] = val;
                    }
                  else
                    {
                      px[_px++] = tmpx + _xstart;
                      py[_py++] = tmpy + _ystart;
                      pg[_pg++] = 0; // Mark the gradient as zero at a junction
                    }
                  pt[_pt++] = _theta[tmpx][tmpy];
                }
              dc.SetEndX(tmpx+_xstart); dc.SetEndY(tmpy+_ystart);

              // Just check whether we have created a trivial edgechain
              // can happen if _min_length_OX = 2
              if( (dc.length()==2) && (dc.GetX(0)==dc.GetX(1)) &&
                 (dc.GetY(0)==dc.GetY(1)) )
                    dc=null;
              else if( dc.length() > 1 )
                {
                  // Create a edge for the image topology
                  v1 = new Vertex(dc.GetX(0),dc.GetY(0));
                  v2 = new Vertex(dc.GetX(dc.length()-1),dc.GetY(dc.length()-1));

                  if(_junction_option_OX)
                    {
                      // Check whether each vertex is a junction
                      V1=Find_vertex(v1);  V2=Find_vertex(v2);

                      // If neither are junctions we may have formed a single
                      // isolated chain that should have common vertex endpoints.
                      int single_chain = 0;
                      if( (V1==null) && (V2==null) )
                        {
                          float dx = dc.GetX(0) - dc.GetX(dc.length()-1);
                          float dy = dc.GetY(0) - dc.GetY(dc.length()-1);
                          if( (dx*dx+dy*dy)<4.0 ) // ie. dist < 2 pixels it is closed
                            {
                              V1 = v1;  V2 = v1;
                              v2=null;
                              single_chain = 1;
                            }
                        }
                      if( single_chain==0 )
                        {
                          if( V1==null )
                            V1 = v1;
                          else
                            v1=null;
                          if( V2==null )
                            V2 = v2;
                          else
                            v2=null;
                        }
                    }
                  else
                    {
                      // We may have formed a single isolated
                      // chain that should have common vertex endpoints.
                      float dx = dc.GetX(0) - dc.GetX(dc.length()-1);
                      float dy = dc.GetY(0) - dc.GetY(dc.length()-1);
                      if((dx*dx+dy*dy) < 4.0) // ie. dist < 2 pixels it is closed
                        {
                          V1 = v1;  V2 = v1;
                          v2=null;
                        }
                      else
                        {
                          V1 = v1; V2 = v2;
                        }
                    }
                  edge = new Edge(V1,V2);

                  // Note that the edge can start and end in the same place.
                  // However, if this is so the DigitalCurve has positive length
                  dc.SetStartX(dc.GetX(0));
                  dc.SetStartY(dc.GetY(0));
                  dc.SetEndX(dc.GetX(dc.length()-1));
                  dc.SetEndY(dc.GetY(dc.length()-1));

                  edge.SetCurve(dc);
                  edges.insertElementAt(edge,0);
                }
              // Clear-up in case there are any trailing dummy points
              xcoords.removeAllElements();  ycoords.removeAllElements();  grad.removeAllElements();
            }
        }
    }
}



//-----------------------------------------------------------------------------
// -- Adds point (x, y) to the current curve, sets its value in the image to
// zero (any point may be included at most once in at most one curve)
// then searches for adjacent pixels (in an order intended to make closed
// curves ordered in a clockwise direction).
// If none are found, and we are not at a junction Join_dotsOX() is called to
// search for non-adjacent pixels.
// If a new point is found, a recursive call is made to Final_followOX() with
// that point.  Protected.
//
// Two strategies can be followed according to _follow_strategy_OX
// if _follow_strategy_OX = 1 then neighbours are checked in the following
//    order (Charlie Rothwell's way)
// @{
// \begin{verbatim}
//         8  7  6
//         1  x  5
//         2  3  4
// \end{verbatim}
// @}
// but if _follow_strategy_OX = 2 then neighbours are checked in the following
//    order (Nic Pillow's way)
// @{
// \begin{verbatim}
//         8  6  7
//         1  x  4
//         3  2  5
// \end{verbatim}
// @}
//
void Final_followOX(int x,
                             int y,
                             Vector xc, //<int> *
                             Vector yc, //CoolList<int>
                             Vector grad,  //CoolList<float> *
                             int reverse)
{
  boolean junction_or_jump = false;

  // Make sure that we do not overun the border of the image
  if ( (x<=0) || (x>=_xsize-1) || (y<=0) || (y>=_ysize-1) )
    return;

  //if (!reverse)
  if (reverse==0)
    {
      xc.addElement(new Integer(x));
      yc.addElement(new Integer(y));
      grad.addElement(new Float( _thin[x][y]));
    }
  _thin[x][y] = 0;

  switch (_follow_strategy_OX)
    {
    case 1:
      // Find one adjacent pixel;  for a closed curve, the method `guarantees'
      // a clockwise ordering of the points in the image sense
      if ( (_thin[x-1][y] >= _edge_min_OX) && (_junction[x-1][y]==0) )
        Final_followOX(x-1,y,xc,yc,grad,0);
      else if ( (_thin[x-1][y+1] >= _edge_min_OX) && (_junction[x-1][y+1]==0) )
        Final_followOX(x-1,y+1,xc,yc,grad,0);
      else if ( (_thin[x][y+1] >= _edge_min_OX) && (_junction[x][y+1]==0) )
        Final_followOX(x,y+1,xc,yc,grad,0);
      else if ( (_thin[x+1][y+1] >= _edge_min_OX) && (_junction[x+1][y+1]==0) )
        Final_followOX(x+1,y+1,xc,yc,grad,0);
      else if ( (_thin[x+1][y] >= _edge_min_OX) && (_junction[x+1][y]==0) )
        Final_followOX(x+1,y,xc,yc,grad,0);
      else if ( (_thin[x+1][y-1] >= _edge_min_OX) && (_junction[x+1][y-1]==0) )
        Final_followOX(x+1,y-1,xc,yc,grad,0);
      else if ( (_thin[x][y-1] >= _edge_min_OX) && (_junction[x][y-1]==0) )
        Final_followOX(x,y-1,xc,yc,grad,0);
      else if ( (_thin[x-1][y-1] >= _edge_min_OX) && (_junction[x-1][y-1]==0) )
        Final_followOX(x-1,y-1,xc,yc,grad,0);
      else
        junction_or_jump = true;
      break;

    case 2: default:
      if ( (_thin[x-1][y] >= _edge_min_OX) && (_junction[x-1][y]==0) )
        Final_followOX(x-1,y,xc,yc,grad,0);
      else if ( (_thin[x][y+1] >= _edge_min_OX) && (_junction[x][y+1]==0) )
        Final_followOX(x,y+1,xc,yc,grad,0);
      else if ( (_thin[x-1][y+1] >= _edge_min_OX) && (_junction[x-1][y+1]==0) )
        Final_followOX(x-1,y+1,xc,yc,grad,0);
      else if ( (_thin[x+1][y] >= _edge_min_OX) && (_junction[x+1][y]==0) )
        Final_followOX(x+1,y,xc,yc,grad,0);
      else if ( (_thin[x+1][y+1] >= _edge_min_OX) && (_junction[x+1][y+1]==0) )
        Final_followOX(x+1,y+1,xc,yc,grad,0);
      else if ( (_thin[x][y-1] >= _edge_min_OX) && (_junction[x][y-1]==0) )
        Final_followOX(x,y-1,xc,yc,grad,0);
      else if ( (_thin[x+1][y-1] >= _edge_min_OX) && (_junction[x+1][y-1]==0) )
        Final_followOX(x+1,y-1,xc,yc,grad,0);
      else if ( (_thin[x-1][y-1] >= _edge_min_OX) && (_junction[x-1][y-1]==0) )
        Final_followOX(x-1,y-1,xc,yc,grad,0);
      else
        junction_or_jump = true;
      break;
    } // end switch

  if (!junction_or_jump)
    return;

  // Else see if there is a junction nearby, and record it. The _chain_no
  // variable is used to prevent the same junction being inserted at both
  // ends of the EdgelChains when reversel occurs next to the junction
  // (in that case there will only be two stored points: the edge and the
  // junction)
  if( _junction[x][y-1]!=0 && ((xc.size()>2)||(_junction[x][y-1]!=_chain_no)) )
    {
      xc.insertElementAt(new Integer(_jx[x][y-1]),0);
      yc.insertElementAt(new Integer(_jy[x][y-1]),0);
      grad.insertElementAt(new Float(_jval),0);
      _junction[x][y-1] = _chain_no;
    }
  else if( _junction[x-1][y]!=0 && ((xc.size()>2)||(_junction[x-1][y]!=_chain_no)) )
    {
      xc.insertElementAt(new Integer(_jx[x-1][y]),0);
      yc.insertElementAt(new Integer(_jy[x-1][y]),0);
      grad.insertElementAt(new Float(_jval),0);
      _junction[x-1][y] = _chain_no;
    }
  else if( _junction[x][y+1]!=0 && ((xc.size()>2)||(_junction[x][y+1]!=_chain_no)) )
    {
      xc.insertElementAt(new Integer(_jx[x][y+1]),0);
      yc.insertElementAt(new Integer(_jy[x][y+1]),0);
      grad.insertElementAt(new Float(_jval),0);
      _junction[x][y+1] = _chain_no;
    }
  else if( _junction[x+1][y]!=0 && ((xc.size()>2)||(_junction[x+1][y]!=_chain_no)) )
    {
      xc.insertElementAt(new Integer(_jx[x+1][y]),0);
      yc.insertElementAt(new Integer(_jy[x+1][y]),0);
      grad.insertElementAt(new Float(_jval),0);
      _junction[x+1][y] = _chain_no;
    }
  else if( _junction[x+1][y-1]!=0 && ((xc.size()>2)||(_junction[x+1][y-1]!=_chain_no)) )
    {
      xc.insertElementAt(new Integer(_jx[x+1][y-1]),0);
      yc.insertElementAt(new Integer(_jy[x+1][y-1]),0);
      grad.insertElementAt(new Float(_jval),0);
      _junction[x+1][y-1] = _chain_no;
    }
  else if( _junction[x-1][y-1]!=0 && ((xc.size()>2)||(_junction[x-1][y-1]!=_chain_no)) )
    {
      xc.insertElementAt(new Integer(_jx[x-1][y-1]),0);
      yc.insertElementAt(new Integer(_jy[x-1][y-1]),0);
      grad.insertElementAt(new Float(_jval),0);
      _junction[x-1][y-1] = _chain_no;
    }
  else if( _junction[x-1][y+1] !=0&& ((xc.size()>2)||(_junction[x-1][y+1]!=_chain_no)) )
    {
      xc.insertElementAt(new Integer(_jx[x-1][y+1]),0);
      yc.insertElementAt(new Integer(_jy[x-1][y+1]),0);
      grad.insertElementAt(new Float(_jval),0);
      _junction[x-1][y+1] = _chain_no;
    }
  else if( _junction[x+1][y+1]!=0 && ((xc.size()>2)||(_junction[x+1][y+1]!=_chain_no)) )
    {
      xc.insertElementAt(new Integer(_jx[x+1][y+1]),0);
      yc.insertElementAt(new Integer(_jy[x+1][y+1]),0);
      grad.insertElementAt(new Float(_jval),0);
      _junction[x+1][y+1] = _chain_no;
    }
  else if ( _join_flag_OX && (xc.size() > 1) )
    {
      // Try to find a pixel to jump to,
      //  and if successful, follow from it
      //int x_c = xc.pop();  int x_p = xc.pop();

      int x_c = ((Integer)xc.elementAt(0)).intValue();
      int x_p = ((Integer)xc.elementAt(0)).intValue();
      xc.removeElementAt(0);xc.removeElementAt(0);
      int y_c = ((Integer)yc.elementAt(0)).intValue();  int y_p = ((Integer)yc.elementAt(0)).intValue();
      yc.removeElementAt(0);yc.removeElementAt(0);
      int dx = x_c - x_p;
      int dy = y_c - y_p;
      xc.insertElementAt(new Integer(x_p),0);  xc.insertElementAt(new Integer(x_c),0);
      yc.insertElementAt(new Integer(y_p),0);  yc.insertElementAt(new Integer(y_c),0);
      int xNew, yNew;
      int[] _xNew= new int[1], _yNew= new int[1];
      boolean rst= (Join_dotsOX(x, y, dx, dy, _xNew, _yNew)!=0);
      xNew=_xNew[0]; yNew=_yNew[0];
      if ((Join_dotsOX(x, y, dx, dy, _xNew, _yNew)!=0) == true)
        Final_followOX(xNew,yNew,xc,yc,grad,0);
    }
}



//-----------------------------------------------------------------------------
//
// -- The point (x, y) is the current curve point;
// (dx, dy) is the increment in position to (x, y) from the preceding point.
// Based on dx and dy, it looks ahead in the same direction (or failing
// that,  similar directions) for another point.
// On success, xNew and yNew are set to the coordinates of the point found
// and true is returned;  otherwise false is returned.  Protected.
//
int Join_dotsOX(int x, int y, int dx, int dy, int[] xNew, int[] yNew)
{
  if ((Math.abs(dx) > 1) || (Math.abs(dy) > 1))      // If dx or dy is too large (> 1), meaning the last
    {                                        //   point was found by jumping, then jumping again
      return 0;                          //   will be too unreliable
    }

 // Make sure that we do not overun the border of the image
 if ( (x<=1) || (x>=_xsize-2) || (y<=1) || (y>=_ysize-2) )
   return 0;

  if (dx==0 && (Math.abs(dy) == 1))
    {
      if (_thin[x][y+2*dy] >= _edge_min_OX) Set_intsOX(xNew, yNew, x, y+2*dy);
      else if (_thin[x+1][y+2*dy] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+1, y+2*dy);
      else if (_thin[x-1][y+2*dy] >= _edge_min_OX) Set_intsOX(xNew, yNew, x-1, y+2*dy);
      else if (_thin[x+2][y+2*dy] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+2, y+2*dy);
      else if (_thin[x-2][y+2*dy] >= _edge_min_OX) Set_intsOX(xNew, yNew, x-2, y+2*dy);
      else if (_thin[x+2][y+dy] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+2, y+dy);
      else if (_thin[x-2][y+dy] >= _edge_min_OX) Set_intsOX(xNew, yNew, x-2, y+dy);
      else if (_thin[x+2][y] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+2, y);
      else if (_thin[x-2][y] >= _edge_min_OX) Set_intsOX(xNew, yNew, x-2, y);
      else return 0;
    }
  else if ((Math.abs(dx) == 1) && dy==0)
    {
      if (_thin[x+2*dx][y] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+2*dx, y);
      else if (_thin[x+2*dx][y+1] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+2*dx, y+1);
      else if (_thin[x+2*dx][y-1] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+2*dx, y-1);
      else if (_thin[x+2*dx][y+2] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+2*dx, y+2);
      else if (_thin[x+2*dx][y-2] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+2*dx, y-2);
      else if (_thin[x+dx][y+2] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+dx, y+2);
      else if (_thin[x+dx][y-2] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+dx, y-2);
      else if (_thin[x][y+2] >= _edge_min_OX) Set_intsOX(xNew, yNew, x, y+2);
      else if (_thin[x][y-2] >= _edge_min_OX) Set_intsOX(xNew, yNew, x, y-2);
      else return 0;
    }
  else if (Math.abs(dx*dy) == 1)
    {
      if (_thin[x+2*dx][y+2*dy] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+2*dx, y+2*dy);
      else if (_thin[x+2*dx][y+dy] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+2*dx, y+dy);
      else if (_thin[x+dx][y+2*dy] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+dx, y+2*dy);
      else if (_thin[x+2*dx][y] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+2*dx, y);
      else if (_thin[x][y+2*dy] >= _edge_min_OX) Set_intsOX(xNew, yNew, x, y+2*dy);
      else if (_thin[x+2*dx][y-dy] >= _edge_min_OX) Set_intsOX(xNew, yNew, x+2*dx, y-dy);
      else if (_thin[x-dx][y+2*dy] >= _edge_min_OX) Set_intsOX(xNew, yNew, x-dx, y+2*dy);
      else return 0;
    }
  else
    return 0;             // Should never be reached, but just in case...

  return 1;
}



//-----------------------------------------------------------------------------
//
// -- Multiply image by scale.  Protected.
//
void Scale_imageOX(float[][] image, float scale)
{
  int x, y;
  for (x=0; x<_xsize; x++)
    for (y=0; y<_ysize; y++)
      image[x][y] = Math.min( (int) (image[x][y]*scale), 255 );
}



//-----------------------------------------------------------------------------
//
// -- Set size of pixels around image (border) to value, so follow can't
// overrun.  Protected.
//
void Set_image_borderOX(float[][] image, int border_size, float value)
{
  int x, y, i;

  for (i=0; i<border_size; i++)
    {
      for (x=0; x<_xsize; x++)
        image[x][i] = image[x][_ysize-1-i] = value;
      for (y=0; y<_ysize; y++)
        image[i][y] = image[_xsize-1-i][y] = value;
    }
}



//-----------------------------------------------------------------------------
//
//  -- Searches for the junctions in the image.  Protected.
//
void Find_junctionsOX()
{
  int x,y,a,b,c,d,e,f,g,h;

  // Reset the junction variables
  _xjunc.removeAllElements();     _yjunc.removeAllElements();
  Zero_int_image(_junction);

  for(x=_border_size_OX;x<_xsize -_border_size_OX;x++)
    for(y=_border_size_OX;y<_ysize-_border_size_OX;y++)
      {
        if( _thin[x][y] < _edge_min_OX )
          continue;

        if( _thin[x-1][y-1] >= _edge_min_OX )  a = 1; else a = 0;
        if( _thin[x][y-1]   >= _edge_min_OX )  b = 1; else b = 0;
        if( _thin[x+1][y-1] >= _edge_min_OX )  c = 1; else c = 0;
        if( _thin[x+1][y]   >= _edge_min_OX )  d = 1; else d = 0;
        if( _thin[x+1][y+1] >= _edge_min_OX )  e = 1; else e = 0;
        if( _thin[x][y+1]   >= _edge_min_OX )  f = 1; else f = 0;
        if( _thin[x-1][y+1] >= _edge_min_OX )  g = 1; else g = 0;
        if( _thin[x-1][y]   >= _edge_min_OX )  h = 1; else h = 0;

        if( (a+b+c+d+e+f+g+h) > 2 )
          {
            _xjunc.insertElementAt(new Integer(x),0);  _yjunc.insertElementAt(new Integer(y),0);
            _junction[x][y] = 1;
          }
      }
}



//-----------------------------------------------------------------------------
//
// -- Locate junction clusters using the following method of hysteresis.
// Protected.
//
void Find_junction_clustersOX()
{
  Vector xcoords=new Vector(),ycoords=new Vector(),xvertices=new Vector(),yvertices=new Vector(),
        xjunc=new Vector(),yjunc=new Vector();   //int
  int x,y,x0=0,y0=0;

  // Find a junction and follow
  //xvertices.clear();  yvertices.clear();
  //xjunc.clear();      yjunc.clear();
  for(x=_border_size_OX;x<_xsize-_border_size_OX;x++)
    for(y=_border_size_OX;y<_ysize-_border_size_OX;y++)
      if( _junction[x][y]!=0 )
        {
          // Each cluster is written to (xcoords,ycooords)
          //xcoords.clear();  ycoords.clear();
          xcoords.removeAllElements();  ycoords.removeAllElements();
          Follow_junctions(x,y,xcoords,ycoords);

          // Find the `centre' of the cluster. This is defined as the
          // junction closest to the centre of gravity of the cluser
          int[] _x0={x0},_y0={y0};
          Cluster_centre_of_gravity(xcoords,ycoords,_x0,_y0);
          x0=_x0[0];y0=_y0[0];

          // Add both the junctions and the new cluster centre to
          // the main lists
          xvertices.insertElementAt(new Integer(x0),0);       yvertices.insertElementAt(new Integer(y0),0);
          for (int xi=0;xi<xcoords.size();xi++) {
            xjunc.insertElementAt(xcoords.elementAt(xi),xi);
            yjunc.insertElementAt(ycoords.elementAt(xi),xi);
          }
        }

  // Reset the junction image - this is order dependent because
  // the cluster centres appear in both lists
  // xjunc.reset();  yjunc.reset();
  while( xjunc.size()!=0 ) {
    int xjunc_v=((Integer)xjunc.elementAt(0)).intValue();
    int yjunc_v=((Integer)yjunc.elementAt(0)).intValue();
    xjunc.removeElementAt(0);yjunc.removeElementAt(0);
    _junction[xjunc_v][yjunc_v] = 1;
  }

  // Construct the list of junction cluster centres
  _vlist.removeAllElements();
  //for(xvertices.reset(),yvertices.reset();xvertices.next(),yvertices.next();)
  for(int vi=0;vi<xvertices.size();vi++)
    {
      Vertex v = new Vertex(((Integer)xvertices.elementAt(vi)).intValue()+_xstart,
            ((Integer)yvertices.elementAt(vi)).intValue()+_ystart);
      _vlist.insertElementAt(v,0);
      _junction[((Integer)xvertices.elementAt(vi)).intValue()][((Integer)yvertices.elementAt(vi)).intValue()] = 2;
    }

  xvertices.removeAllElements();  yvertices.removeAllElements();
}





//-- If the ROI has a mask, all edges are checked if they lie in the mask.
// If an edge does not lie completely in the mask, new edges are created: the parts
// of the original edge that do lie in the mask.


/*
void AdjustForMask(Vector edges)
  {
  if(imageROI && imageROI.GetMask())
    {

    Vector edges_old = edges;
    Vector edges_new = new Vector(); //P<Edge*>;

    RectROI* ROI = imageROI;
    int xo = ROI.GetOrigX();
    int yo = ROI.GetOrigY();

    for (int i=0;i<edges_old.length();++i)
      {
      Edge edge = edges_old.get(i);
      DigitalCurve dc = (DigitalCurve) edge.GetCurve();

      Vector curvelist = CurveOps.IntersectMask(dc,ROI.GetMask(),xo,yo);//CoolListP<Curve*>

      // curvelist is a list of digital curves that are in the mask.
      // See CurveOps for details.

      Vector nedgelist;//CoolListP<Edge*>

      for (curvelist.reset(); curvelist.next(); )
        {
        Edge* edge1 = new Edge(curvelist.value());
        edge1.Protect();
        nedgelist.push(edge1);
        }

      // nedgelist is a list of new edges, created from the original edge.


      // Check for topology:
      // All new edges have new vertices.
      // It's possible that the old vertices of the original edge can be reused.
      // If V1 of the original edge is in the mask, the first point of the first
      // edge in nedgelist is this vertex.
      // If V2 of the original edge is in the mask, the last point of the last
      // edge in nedgelist is this vertex.

      Vertex V1 = edge.GetV1();
      Vertex V2 = edge.GetV2();


      if (nedgelist.length()>0 && V1)
        if (bytePixel(ROI.GetMask(), ( (int) V1.GetX() - xo), ( (int) V1.GetY() - yo)))
          nedgelist.get(0).SetV1(V1);

      if (nedgelist.length()>0 && V2)
        if (bytePixel(ROI.GetMask(), ( (int) V2.GetX() - xo), ( (int) V2.GetY() - yo)))
          nedgelist.get(nedgelist.length()-1).SetV2(V2);

      for (nedgelist.reset(); nedgelist.next(); )
        edges_new.push(nedgelist.value());

      edge.UnProtect();

      }
    edges.clear();
    for (int ei=0;ei<edges_new.size();ei++) edges.addElement(edges_new.elementAt(ei));
    //*edges += *edges_new;
    //delete edges_new;
    edges_new=null;
    return;
    }
  else return;
  }
  */

void Set_intsOX(int[] int1, int[] int2, int val1, int val2)
{
  int1[0] = val1;
  int2[0] = val2;
}
void Zero_int_image(int[][] image)
{
    int x,y;
    int w=image[0].length, h=image.length;


    for(x=0;x<w;x++)
       {
        for(y=0;y<h;y++)
            image[y][x] = 0;
       }
}

static float DUMMYTHETA = 10000;
void initFloatImage(float[][] theta, int sizex, int sizey)
{
  for(int x =0; x < sizex; x++)
    for(int y =0; y < sizey; y++)
      theta[x][y] = DUMMYTHETA;
}


    //first m, then n.  means number of column go first
    float[][] Make_float_image(int m, int n)
    {
        return new float[m][n];


    }

    int[][] Make_int_image(int m, int n)
{
  return new int[m][n];
}

float[] Make_float_vector(int m)
{
    return new float[m];
}
//-----------------------------------------------------------------------------
//
// -- Following routine looking for searching out junction clusters.  Protected.
//
void Follow_junctions(int x, int y, Vector xc, Vector yc)    //xc,yc=int
{

  // Add the current junction to the coordinate lists, and delete from
  // the junction image
  xc.insertElementAt(new Integer(x),0);  yc.insertElementAt(new Integer(y),0);
  _junction[x][y] = 0;

  // Now recursively look for connected eight-neighbours
  if( _junction[x][y-1]!=0 )
    Follow_junctions(x,y-1,xc,yc);
  if( _junction[x-1][y]!=0 )
    Follow_junctions(x-1,y,xc,yc);
  if( _junction[x][y+1]!=0 )
    Follow_junctions(x,y+1,xc,yc);
  if( _junction[x+1][y]!=0 )
    Follow_junctions(x+1,y,xc,yc);
  if( _junction[x+1][y-1]!=0 )
    Follow_junctions(x+1,y-1,xc,yc);
  if( _junction[x-1][y-1]!=0 )
    Follow_junctions(x-1,y-1,xc,yc);
  if( _junction[x-1][y+1]!=0 )
    Follow_junctions(x-1,y+1,xc,yc);
  if( _junction[x+1][y+1] !=0)
    Follow_junctions(x+1,y+1,xc,yc);
}

//-----------------------------------------------------------------------------
//
// -- Finds which member of the lists lies closest to the centre
// of gravity of the lists.  Protected.
//
void Cluster_centre_of_gravity(Vector xc, Vector yc,
                             int[] _x0, int[] _y0)
{
  int x0=_x0[0],y0=_y0[0];
  if( xc==null||xc.size()==0 )
    return;

  // First find the CofG
  double x=0.0,y=0.0;

  for(int i=0;i<xc.size();i++) //xc.reset(),yc.reset();xc.next(),yc.next();)
    {
      x += ((Integer)xc.elementAt(i)).intValue();  y += ((Integer)yc.elementAt(i)).intValue();
    }
  x /= xc.size();  y /= yc.size();

  // Now find the point closest to the CofG
  float dist,newdist;
  dist = _image.length*_image[0].length; // A number larger than the image size
  for(int i=0;i<xc.size();i++){ //for(xc.reset(),yc.reset();xc.next(),yc.next();)
      int xv = ((Integer)xc.elementAt(i)).intValue(),  yv = ((Integer)yc.elementAt(i)).intValue();
    if( (newdist=(float)hypot((float)(x-xv),(float)(y-yv))) < dist )
      {
                x0 = xv;  y0 = yv;
                dist = newdist;
      }
}

  // Set up the (_jx,_jy) arrays to point to the cluster centre
  for(int i=0;i<xc.size();i++) //for(xc.reset(),yc.reset();xc.next(),yc.next();) {
    {
      int xv = ((Integer)xc.elementAt(i)).intValue(),  yv = ((Integer)yc.elementAt(i)).intValue();
      _jx[xv][yv] = x0;
      _jy[xv][yv] = y0;
    }
    _x0[0]=x0;_y0[0]=y0;
}


    class LINK {
        public int to;
        public LINK nextl;
    }

    float hypot(float x, float y)
    {
        return (float)Math.sqrt(x*x+y*y);
    }

    void reverse(Vector v)
    {
        int len= v.size();
        for (int i=0;i<len-1;i++) {
            Object o= v.lastElement();
            v.removeElementAt(len-1);
            v.insertElementAt(o, i);
        }
    }

    Vertex Find_vertex(Vertex vertex)
{
    if (_vlist.contains(vertex)) return vertex;
    return null;
  /*
  //for(_vlist->reset();_vlist->next();)
  for(int i=0 ;i<_vlist.size();i++)
    if( vertex == *(_vlist->value()) )
      return(_vlist->value());

  // If we get this far there is no match
  return(NULL);
  */

}

//-----------------------------------------------------------------------------
//
// -- Zeros a floating point image.  Protected.
//
void Zero_float_image(float[][] image)
{
    int h=image.length, w=image[0].length;
    int x,y;

    for(x=0;x<w;x++)
       {
        for(y=0;y<h;y++)
            image[y][x] = 0;
       }
}
int PixelSize(int[][] image)
{
    return 8; //8 means gray-level image
    //return 24;    //24 means color image
}
int bytePixel(int[][] image,int x,int y)
{
    return image[y][x];
}
int shortPixel(int[][] image,int x,int y)
{
    return image[y][x];
}

int colorPixel(int[][] image,int x,int y)
{
    int v=image[y][x];
    int r=(v>>16)&0xff;
    int g=(v>>8)&0xff;
    int b=(v)&0xff;
    return (int)(Math.round(0.299*r+0.587*g+0.114*b));
}

    void display_float_image(String name, float[][] data) {
        int h=data.length, w=data[0].length;
        int[] img= new int[w*h];
        int index=0;
        for (int y=0;y<h;y++) {
            for (int x=0;x<w;x++) {
                img[index++]=(int)Math.round(data[y][x]);
            }
        }
        yi.util.ImageUtil.makeImageFromGray(img);
        (new DisplayFrame(name,img,w,h)).show();
    }
}