package yi.iu.processor.canny;

import java.awt.*;
import yi.iu.ui.DisplayFrame;
import yi.iu.processor.channel.*;

public class CannyEdgeDetector extends yi.iu.processor.ImageProcessor
{
    public double sigma= 1.6;        // Standard deviation of the gaussian kernel.
    public double tlow= 0.4;       // Fraction of the high threshold in hysteresis.
    public double thigh= 0.8;      // High hysteresis threshold control. The actual
                            // threshold is the (100 * thigh) percentage point
                            // in the histogram of the magnitude of the
                            // gradient image that passes non-maximal
                            // suppression.
                            // if thigh>=1, [1..255], it is the absolute value
    public static int NOEDGE= 0;
    public static int POSSIBLE_EDGE= 128;
    public static int EDGE= 255;
    static double BOOSTBLURFACTOR= 90.0;

    protected boolean VERBOSE=false;

    public CannyEdgeDetector() {}

    public CannyEdgeDetector(double s, double l, double h) {
    	sigma=s; tlow=l; thigh=h;
    }
    public int getArguNum() { return 3;}

    public void setArgument(int i, String st)
    {
        if (st.length()==0) return;
        if (i==0) sigma= Double.valueOf(st).doubleValue();
        else if (i==1) tlow= Double.valueOf(st).doubleValue();
        else if (i==2) thigh= Double.valueOf(st).doubleValue();
    }

    public String getArgument(int i)
    {
        if (i==0) return (new Double(sigma)).toString();
        else if (i==1) return (new Double(tlow)).toString();
        else if (i==2) return (new Double(thigh)).toString();
        return null;
    }

    public String getArguName(int i)
    {
        if (i==0) return "Sigma";
        else if (i==1) return "Low";
        else if (i==2) return "High";
        return null;
    }

    public int[] getProcessedPseudoColor(int source[], Dimension d)
    {
        int[] rst= getProcessed(source,d);
        yi.util.ImageUtil.makeImageFromBinary(rst);
        DisplayFrame f= new DisplayFrame(this, rst, d.width, d.height);
        f.show();
        //yi.util.io.OutputResult.outputGrayImage(rst, d.width, d.height);

        return rst;
    }

    public int[] getProcessed(int source[], Dimension d)
    {
        int rows=d.height, cols=d.width;        // The dimensions of the image
        IntensityMaker im= new IntensityMaker();
        int[] image= im.getProcessed(source, d);

        // Perform the edge detection. All of the work takes place here.
        return canny(image, d);
    }

    // PURPOSE: To perform canny edge detection.
    public int[] canny(int[] image, Dimension d)
    {
        // Perform gaussian smoothing on the image using the input standard deviation.
        if(VERBOSE) System.err.println("Smoothing the image using a gaussian kernel.");
        int[] smoothedim= (sigma>0)?gaussian_smooth(image, d):image;

        // Compute the first derivative in the x and y directions.
        if(VERBOSE) System.err.println("Computing the X and Y first derivatives.");
        int[] delta_x= new int[d.width*d.height];
        int[] delta_y= new int[d.width*d.height];
        derrivative_x_y(smoothedim, d, delta_x, delta_y);

        // This option to write out the direction of the edge gradient was added
        // to make the information available for computing an edge quality figure
        // of merit.
        // Compute the direction up the gradient, in radians that are
        // specified counteclockwise from the positive x-axis.
        // double[] dir_radians= radian_direction(delta_x, delta_y, d, -1, -1);

        // Compute the magnitude of the gradient.
        if(VERBOSE) System.err.println("Computing the magnitude of the gradient.");
        int[] magnitude= magnitude_x_y(delta_x, delta_y, d);

        // Perform non-maximal suppression.
        if(VERBOSE) System.err.println("Doing the non-maximal suppression.");
        int[] nms= non_max_supp(magnitude, delta_x, delta_y, d);

        // Use hysteresis to mark the edge pixels.
        if(VERBOSE) System.err.println("Doing hysteresis thresholding.");
        int[] edge= apply_hysteresis(magnitude, nms, d);

        //for output the direction and magnitude of pixels on edges
        // this part is used to make the pic in paper
        /*
        double[] dir_radians= radian_direction(delta_x, delta_y, d, -1, -1);
        for (int i=0;i<d.width*d.height;i++) {
            if (edge[i]!=EDGE) {
                magnitude[i]=0;
                dir_radians[i]=0;
            }
            else {
                dir_radians[i]=dir_radians[i]*180/Math.PI;
                if (dir_radians[i]>=180) dir_radians[i]=360-dir_radians[i];
                dir_radians[i]=dir_radians[i]-90;
            }
        }
        double dir_189[]= new double[d.width];
        double dir_190[]= new double[d.width];
        for (int i=0;i<d.width;i++) {
            dir_189[i]=dir_radians[189*d.width+i];
            dir_190[i]=dir_radians[190*d.width+i];
        }
        yi.util.io.OutputResult.outputResult("dir",dir_189,1,d.width);
        yi.util.io.OutputResult.outputResult("dir",dir_190,1,d.width);
        //yi.util.io.OutputResult.outputResult("mag",magnitude, d.width, d.height);
        */

        delta_x= null; delta_y= null;
        magnitude= null; nms=null;
        System.gc();
        return edge;
    }

    // PURPOSE: To perform canny edge detection.
    public void canny(int[] image, Dimension d,
      int[] edges, double[] directions, int[] magnitudes)
    {
        // Perform gaussian smoothing on the image using the input standard deviation.
        if(VERBOSE) System.err.println("Smoothing the image using a gaussian kernel.");
        int[] smoothedim= (sigma>0)?gaussian_smooth(image, d):image;

        // Compute the first derivative in the x and y directions.
        if(VERBOSE) System.err.println("Computing the X and Y first derivatives.");
        int[] delta_x= new int[d.width*d.height];
        int[] delta_y= new int[d.width*d.height];
        derrivative_x_y(smoothedim, d, delta_x, delta_y);

        // This option to write out the direction of the edge gradient was added
        // to make the information available for computing an edge quality figure
        // of merit.
        // Compute the direction up the gradient, in radians that are
        // specified counteclockwise from the positive x-axis.
        // double[] dir_radians= radian_direction(delta_x, delta_y, d, -1, -1);

        // Compute the magnitude of the gradient.
        if(VERBOSE) System.err.println("Computing the magnitude of the gradient.");
        int[] magnitude= magnitude_x_y(delta_x, delta_y, d);

        // Perform non-maximal suppression.
        if(VERBOSE) System.err.println("Doing the non-maximal suppression.");
        int[] nms= non_max_supp(magnitude, delta_x, delta_y, d);

        // Use hysteresis to mark the edge pixels.
        if(VERBOSE) System.err.println("Doing hysteresis thresholding.");
        int[] edge= apply_hysteresis(magnitude, nms, d);

        //for output the direction and magnitude of pixels on edges
        // this part is used to make the pic in paper
        double[] dir_radians= radian_direction(delta_x, delta_y, d, -1, -1);
        for (int i=0;i<d.width*d.height;i++) {
            dir_radians[i]=dir_radians[i]*180/Math.PI;
        }

        for (int i=0;i<d.width*d.height;i++) {
            if (edge[i]!=EDGE) {
                magnitude[i]=0;
                dir_radians[i]=0;
            }
            else {
                //dir_radians[i]=dir_radians[i]*180/Math.PI;
                // convert to Qie1-Vector
                if (dir_radians[i] <270) dir_radians[i] = dir_radians[i] + 90;
                else dir_radians[i] = dir_radians[i] - 270;
            }
        }

        /*
        // output the direction information of horizontal line
        double dir_189[]= new double[d.width];
        double dir_190[]= new double[d.width];
        for (int i=0;i<d.width;i++) {
            dir_189[i]=dir_radians[189*d.width+i];
            dir_190[i]=dir_radians[190*d.width+i];
        }
        yi.util.io.OutputResult.outputResult("dir",dir_189,1,d.width);
        yi.util.io.OutputResult.outputResult("dir",dir_190,1,d.width);
        */

        for (int i=0; i<d.width*d.height; i++) {
          edges[i] = edge[i];
          directions[i] = dir_radians[i];
          magnitudes[i] = magnitude[i];
        }

        delta_x= null; delta_y= null;
        magnitude= null; nms=null;
        System.gc();
    }

    // Procedure: radian_direction
    // Purpose: To compute a direction of the gradient image from component dx and
    // dy images. Because not all derriviatives are computed in the same way, this
    // code allows for dx or dy to have been calculated in different ways.
    //
    // FOR X:  xdirtag = -1  for  [-1 0  1]
    //         xdirtag =  1  for  [ 1 0 -1]
    //
    // FOR Y:  ydirtag = -1  for  [-1 0  1]'
    //         ydirtag =  1  for  [ 1 0 -1]'
    //
    // The resulting angle is in radians measured counterclockwise from the
    // xdirection. The angle points "up the gradient".
    protected double[] radian_direction(int[] delta_x, int[] delta_y, Dimension d,
                int xdirtag, int ydirtag)
    {
        int rows= d.height, cols= d.width;
        int r,c,pos;
        double[] dirim= new double[rows*cols];

        for(r=0,pos=0;r<rows;r++){
            for(c=0;c<cols;c++,pos++){
                double dx = (double)delta_x[pos];
                double dy = (double)delta_y[pos];
                if(xdirtag == 1) dx = -dx;
                if(ydirtag == -1) dy = -dy;
                dirim[pos] = (double)angle_radians(dx, dy);
            }
        }
        return dirim;
    }

    // FUNCTION: angle_radians
    // PURPOSE: This procedure computes the angle of a vector with components x and
    // y. It returns this angle in radians with the answer being in the range
    // 0 <= angle <2*PI.
    public double angle_radians(double x, double y)
    {
        double xu = Math.abs(x);
        double yu = Math.abs(y);

        if((xu == 0) && (yu == 0)) return(0);
        double ang = Math.atan(yu/xu);
        if(x >= 0){
            if(y >= 0) return(ang);
            else return(2*Math.PI - ang);
        }
        else{
            if(y >= 0) return(Math.PI - ang);
            else return(Math.PI + ang);
        }
    }

    // PROCEDURE: magnitude_x_y
    // PURPOSE: Compute the magnitude of the gradient. This is the square root of
    // the sum of the squared derivative values.
    public int[] magnitude_x_y(int[] delta_x, int[] delta_y, Dimension d)
    {
        int rows= d.height, cols= d.width;
        int r, c, pos, sq1, sq2;
        int[] magnitude= new int[rows*cols];
        for(r=0,pos=0;r<rows;r++){
            for(c=0;c<cols;c++,pos++){
                sq1 = (int)delta_x[pos] * (int)delta_x[pos];
                sq2 = (int)delta_y[pos] * (int)delta_y[pos];
                magnitude[pos] = (int)Math.round(Math.sqrt((double)sq1 + (double)sq2));
            }
        }
        return magnitude;
    }

    // PROCEDURE: derrivative_x_y
    // PURPOSE: Compute the first derivative of the image in both the x any y
    // directions. The differential filters that are used are:
    //
    //                                   -1
    // dx =  -1 0 +1     and       dy =  0
    //                                   +1
    public void derrivative_x_y(int[] smoothedim, Dimension d, int[] delta_x, int[] delta_y)
    {
        int rows= d.height, cols= d.width;
        int r, c, pos;
        // Compute the x-derivative. Adjust the derivative at the borders to avoid losing pixels.
        if(VERBOSE) System.err.println("   Computing the X-direction derivative.");
        for(r=0;r<rows;r++){
            pos = r * cols;
            delta_x[pos] = smoothedim[pos+1] - smoothedim[pos];
            pos++;
            for(c=1;c<(cols-1);c++,pos++){
                delta_x[pos] = smoothedim[pos+1] - smoothedim[pos-1];
            }
            delta_x[pos] = smoothedim[pos] - smoothedim[pos-1];
        }
        // Compute the y-derivative. Adjust the derivative at the borders to avoid
        // losing pixels.
        if(VERBOSE) System.err.println("   Computing the Y-direction derivative.");
        for(c=0;c<cols;c++){
            pos = c;
            delta_y[pos] = smoothedim[pos+cols] - smoothedim[pos];
            pos += cols;
            for(r=1;r<(rows-1);r++,pos+=cols){
                delta_y[pos] = smoothedim[pos+cols] - smoothedim[pos-cols];
            }
            delta_y[pos] = smoothedim[pos] - smoothedim[pos-cols];
        }
    }

    // PROCEDURE: gaussian_smooth
    // PURPOSE: Blur an image with a gaussian filter.
    public int[] gaussian_smooth(int[] image, Dimension d)
    {
        int rows= d.height, cols= d.width;
        int r, c, rr, cc;   // Counter variables.
        double dot,         // Dot product summing variable.
            sum;            // Sum of the kernel weights variable.

        // Create a 1-dimensional gaussian smoothing kernel.
        if(VERBOSE) System.err.println("   Computing the gaussian smoothing kernel.");
        double[] kernel=make_gaussian_kernel();  // A one dimensional gaussian kernel.
        int windowsize= kernel.length;
        int center = windowsize / 2;

        // Allocate a temporary buffer image and the smoothed image.
        double[] tempim = new double[rows*cols];
        int[] smoothedim = new int[rows*cols];

        // Blur in the x - direction.
        if(VERBOSE) System.err.println("   Bluring the image in the X-direction.");
        for(r=0;r<rows;r++){
            for(c=0;c<cols;c++){
                dot = 0.0;
                sum = 0.0;
                for(cc=(-center);cc<=center;cc++){
                    if(((c+cc) >= 0) && ((c+cc) < cols)){
                        dot += (double)image[r*cols+(c+cc)] * kernel[center+cc];
                        sum += kernel[center+cc];
                    }
                }
                tempim[r*cols+c] = dot/sum;
            }
        }

        // Blur in the y - direction.
        if(VERBOSE) System.err.println("   Bluring the image in the Y-direction.");
        for(c=0;c<cols;c++){
            for(r=0;r<rows;r++){
                sum = 0.0;
                dot = 0.0;
                for(rr=(-center);rr<=center;rr++){
                    if(((r+rr) >= 0) && ((r+rr) < rows)){
                        dot += tempim[(r+rr)*cols+c] * kernel[center+rr];
                        sum += kernel[center+rr];
                    }
                }
                smoothedim[r*cols+c] = (int)(dot*BOOSTBLURFACTOR/sum + 0.5);
            }
        }
        tempim= null;
        System.gc();
        return smoothedim;
    }

    // PROCEDURE: make_gaussian_kernel
    // PURPOSE: Create a one dimensional gaussian kernel.
    double[] make_gaussian_kernel()
    {
        int windowsize = (int)(1 + 2 * Math.ceil(2.5 * sigma)+0.5);
        int center = windowsize / 2;
        int i;
        double x, fx, sum=0.0;

        if(VERBOSE) System.err.println("      The kernel has %d elements."+windowsize);
        double[] kernel = new double[windowsize];
        for(i=0;i<windowsize;i++){
            x = (double)(i - center);
            fx = Math.exp(-0.5*x*x/(sigma*sigma)) / (sigma * Math.sqrt(2*Math.PI));
            kernel[i] = fx;
            sum += fx;
        }

        for(i=0;i<windowsize;i++) kernel[i] /= sum;
        if(VERBOSE){
            System.err.println("The filter coefficients are:");
            for(i=0;i<windowsize;i++)
                System.err.println("kernel["+i+"] = "+kernel[i]);
        }
        return kernel;
    }

    // PROCEDURE: follow_edges
    // PURPOSE: This procedure edges is a recursive routine that traces edgs along
    // all paths whose magnitude values remain above some specifyable lower
    // threshhold.
    int stack_depth = 0;
    final static int MAX_STACK_DEPTH = 2048;
    public void follow_edges(int[] edge, int[] mag, int edgemapptr, int lowval, int cols)
    {
        int tempmagptr, tempmapptr;

        stack_depth++;
        if (stack_depth>MAX_STACK_DEPTH) return;

        int i;
        double thethresh;
        int[] x= new int[8];
        int[] y= new int[8];
        x[0]=x[1]=x[7]=1; x[2]=x[6]=0; x[3]=x[4]=x[5]=-1;
        y[0]=y[4]=0; y[1]=y[2]=y[3]=1; y[5]=y[6]=y[7]=-1;

        for(i=0;i<8;i++){
            tempmapptr = edgemapptr - y[i]*cols + x[i];
            if((edge[tempmapptr] == POSSIBLE_EDGE) && (mag[tempmapptr] > lowval)){
                edge[tempmapptr] = EDGE;
                follow_edges(edge, mag, tempmapptr, lowval, cols);
            }
        }
    }

    // PROCEDURE: apply_hysteresis
    // PURPOSE: This routine finds edges that are above some high threshhold or
    // are connected to a high pixel by a path of pixels greater than a low
    // threshold.
    public int[] apply_hysteresis(int[] mag, int[] nms, Dimension d)
    {
        int rows= d.height, cols= d.width;
        int r, c, pos, numedges, lowcount, highcount, lowthreshold, highthreshold,
            i,  rr, cc;
        int maximum_mag=0, sumpix;

        // Initialize the edge map to possible edges everywhere the non-maximal
        // suppression suggested there could be an edge except for the border. At
        // the border we say there can not be an edge because it makes the
        // follow_edges algorithm more efficient to not worry about tracking an
        // edge off the side of the image.
        int[] edge= new int[rows*cols];
        for(r=0,pos=0;r<rows;r++){
            for(c=0;c<cols;c++,pos++){
                if(nms[pos] == POSSIBLE_EDGE) edge[pos] = POSSIBLE_EDGE;
                else edge[pos] = NOEDGE;
            }
        }

        for(r=0,pos=0;r<rows;r++,pos+=cols){
            edge[pos] = NOEDGE;
            edge[pos+cols-1] = NOEDGE;
        }
        pos = (rows-1) * cols;
        for(c=0;c<cols;c++,pos++){
            edge[c] = NOEDGE;
            edge[pos] = NOEDGE;
        }

        if (thigh<1) {
            // Compute the histogram of the magnitude image. Then use the histogram to
            // compute hysteresis thresholds.
            int[] hist= new int[32768];
            for(r=0,pos=0;r<rows;r++){
                for(c=0;c<cols;c++,pos++){
                    if(edge[pos] == POSSIBLE_EDGE) {
                      if (mag[pos]>=32768)
                         System.err.println("mag out of range: row="+r+", col="+c+", mag="+mag[pos]);

                      hist[mag[pos]]++;
                    }
                }
            }

            // Compute the number of pixels that passed the nonmaximal suppression.
            for(r=1,numedges=0;r<32768;r++){
                if(hist[r] != 0) maximum_mag = r;
                numedges += hist[r];
            }
            highcount = (int)(numedges * thigh + 0.5);

            // Compute the high threshold value as the (100 * thigh) percentage point
            // in the magnitude of the gradient histogram of all the pixels that passes
            // non-maximal suppression. Then calculate the low threshold as a fraction
            // of the computed high threshold value. John Canny said in his paper
            // "A Computational Approach to Edge Detection" that "The ratio of the
            // high to low threshold in the implementation is in the range two or three
            // to one." That means that in terms of this implementation, we should
            // choose tlow ~= 0.5 or 0.33333.
            r = 1;
            numedges = hist[1];
            while((r<(maximum_mag-1)) && (numedges < highcount)){
                r++;
                numedges += hist[r];
            }
            highthreshold = r;
        }
        else highthreshold = (int)Math.round(thigh*thigh);

        lowthreshold = (int)(highthreshold * tlow + 0.5);
        if(VERBOSE){
            System.err.println("The input low and high fractions of "
                +tlow+" and "+thigh+" computed to");
            System.err.println("magnitude of the gradient threshold values of: "
                +lowthreshold+", "+ highthreshold);
        }

        // This loop looks for pixels above the highthreshold to locate edges and
        // then calls follow_edges to continue the edge.
        for(r=0,pos=0;r<rows;r++){
            for(c=0;c<cols;c++,pos++){
                if((edge[pos] == POSSIBLE_EDGE) && (mag[pos] >= highthreshold)){
                    edge[pos] = EDGE;
                    stack_depth = 0;
                    follow_edges(edge, mag, pos, lowthreshold, cols);
                }
            }
        }

        // Set all the remaining possible edges to non-edges.
        for(r=0,pos=0;r<rows;r++){
            for(c=0;c<cols;c++,pos++) if(edge[pos] != EDGE) edge[pos] = NOEDGE;
        }
        return edge;
    }

    // PROCEDURE: non_max_supp
    // PURPOSE: This routine applies non-maximal suppression to the magnitude of
    // the gradient image.
    public int[] non_max_supp(int[] mag, int[] gradx, int[] grady, Dimension d)
    {
        int nrows= d.height, ncols= d.width;
        int rowcount, colcount,count;
        int magrowptr, magptr;
        int gxrowptr, gxptr;
        int gyrowptr, gyptr;
        int z1,z2;
        int m00,gx,gy;
        double mag1,mag2,xperp,yperp;
        int resultrowptr, resultptr;

        int[] result= new int[nrows*ncols];

        // Suppress non-maximum points.
        for(rowcount=1, magrowptr=ncols+1, gxrowptr=ncols+1,
            gyrowptr=ncols+1, resultrowptr=ncols+1;
            rowcount<nrows-2;
            rowcount++,magrowptr+=ncols,gyrowptr+=ncols,gxrowptr+=ncols,
            resultrowptr+=ncols){
            for(colcount=1,magptr=magrowptr,gxptr=gxrowptr,gyptr=gyrowptr,
                resultptr=resultrowptr;colcount<ncols-2;
                colcount++,magptr++,gxptr++,gyptr++,resultptr++){
                m00 = mag[magptr];
                if(m00 == 0){ result[resultptr] = NOEDGE; continue; }
                xperp = -(gx = gradx[gxptr])/((double)m00);
                yperp = (gy = grady[gyptr])/((double)m00);
                if(gx >= 0){
                    if(gy >= 0){
                        if (gx >= gy)
                        {
                            // 111
                            // Left point
                            z1 = mag[magptr-1];
                            z2 = mag[magptr - ncols - 1];
                            mag1 = (m00 - z1)*xperp + (z2 - z1)*yperp;
                            // Right point
                            z1 = mag[magptr + 1];
                            z2 = mag[magptr + ncols + 1];
                            mag2 = (m00 - z1)*xperp + (z2 - z1)*yperp;
                        }
                        else
                        {
                            // 110
                            // Left point
                            z1 = mag[magptr - ncols];
                            z2 = mag[magptr - ncols - 1];
                            mag1 = (z1 - z2)*xperp + (z1 - m00)*yperp;
                            // Right point
                            z1 = mag[magptr + ncols];
                            z2 = mag[magptr + ncols + 1];
                            mag2 = (z1 - z2)*xperp + (z1 - m00)*yperp;
                        }
                    }
                    else
                    {
                        if (gx >= -gy)
                        {
                            // 101
                            // Left point
                            z1 = mag[magptr-1];
                            z2 = mag[magptr + ncols - 1];
                            mag1 = (m00 - z1)*xperp + (z1 - z2)*yperp;
                            // Right point
                            z1 = mag[magptr + 1];
                            z2 = mag[magptr - ncols + 1];
                            mag2 = (m00 - z1)*xperp + (z1 - z2)*yperp;
                        }
                        else
                        {
                            // 100
                            // Left point
                            z1 = mag[magptr + ncols];
                            z2 = mag[magptr + ncols - 1];
                            mag1 = (z1 - z2)*xperp + (m00 - z1)*yperp;
                            // Right point
                            z1 = mag[magptr - ncols];
                            z2 = mag[magptr - ncols + 1];
                            mag2 = (z1 - z2)*xperp  + (m00 - z1)*yperp;
                        }
                    }
                }
                else
                {
                    if ((gy = grady[gyptr]) >= 0)
                    {
                        if (-gx >= gy)
                        {
                            // 011
                            // Left point
                            z1 = mag[magptr + 1];
                            z2 = mag[magptr - ncols + 1];
                            mag1 = (z1 - m00)*xperp + (z2 - z1)*yperp;
                            // Right point
                            z1 = mag[magptr-1];
                            z2 = mag[magptr + ncols - 1];
                            mag2 = (z1 - m00)*xperp + (z2 - z1)*yperp;
                        }
                        else
                        {
                            // 010
                            // Left point
                            z1 = mag[magptr - ncols];
                            z2 = mag[magptr - ncols + 1];
                            mag1 = (z2 - z1)*xperp + (z1 - m00)*yperp;
                            // Right point
                            z1 = mag[magptr + ncols];
                            z2 = mag[magptr + ncols - 1];
                            mag2 = (z2 - z1)*xperp + (z1 - m00)*yperp;
                        }
                    }
                    else
                    {
                        if (-gx > -gy)
                        {
                            // 001
                            // Left point
                            z1 = mag[magptr + 1];
                            z2 = mag[magptr + ncols + 1];
                            mag1 = (z1 - m00)*xperp + (z1 - z2)*yperp;
                            // Right point
                            z1 = mag[magptr-1];
                            z2 = mag[magptr - ncols - 1];
                            mag2 = (z1 - m00)*xperp + (z1 - z2)*yperp;
                        }
                        else
                        {
                            // 000
                            // Left point
                            z1 = mag[magptr + ncols];
                            z2 = mag[magptr + ncols + 1];
                            mag1 = (z2 - z1)*xperp + (m00 - z1)*yperp;
                            // Right point
                            z1 = mag[magptr - ncols];
                            z2 = mag[magptr - ncols - 1];
                            mag2 = (z2 - z1)*xperp + (m00 - z1)*yperp;
                        }
                    }
                }
                // Now determine if the current point is a maximum point
                if ((mag1 > 0.0) || (mag2 > 0.0)) {
                    result[resultptr] = NOEDGE;
                }
                else {
                    if (mag2 == 0.0)
                        result[resultptr] = NOEDGE;
                    else
                        result[resultptr] = POSSIBLE_EDGE;
                }
                }
            }
            return result;
    }

    /*******************************************************************************
    * Procedure: radian_direction
    * Purpose: To compute a direction of the gradient image from component dx and
    * dy images. Because not all derriviatives are computed in the same way, this
    * code allows for dx or dy to have been calculated in different ways.
    *
    * FOR X:  xdirtag = -1  for  [-1 0  1]
    *         xdirtag =  1  for  [ 1 0 -1]
    *
    * FOR Y:  ydirtag = -1  for  [-1 0  1]'
    *         ydirtag =  1  for  [ 1 0 -1]'
    *
    * The resulting angle is in radians measured counterclockwise from the
    * xdirection. The angle points "up the gradient".
    *******************************************************************************/
    public double[] radian_direction(int[] delta_x, int[] delta_y)
    {
        double[] dir_radians = new double[delta_x.length];

        for(int i=0;i<dir_radians.length;i++){
            double dx = (double)delta_x[i];
            double dy = (double)delta_y[i];
            dir_radians[i] = (double)angle_radians(dx, dy);
        }
        return dir_radians;
    }


}