import ij.*;
import java.io.*;
import ij.plugin.filter.PlugInFilter;
import ij.gui.*;
import ij.process.*;
import java.awt.*;
import java.awt.image.*;
import ij.text.TextWindow;
import java.util.*;


/** Implements ImageJ's Analyze Particles command.

*/
public class Descriptor_Moment implements PlugInFilter {

    private RuffleParticle[][] m_nParticleRec;
    private RuffleParticle[] m_lstParticle;
    private ImageStack StackParticle;
    private ImagePlus StackImp;

    private int m_nHeight, m_nWidth;
    private double m_dMinThr, m_dMaxThr;
    private TextWindow tw;


    private static int staticMinSize = 25;
    private static int staticMaxSize = 10000;
    private static double staticMinAxisLen = 1;
    private static int staticLowerThr = 1;

    private int minSize;
    private int maxSize;
    private double MinAxisLen;
    private int LowerThr;

	public int setup(String arg, ImagePlus imp)
        {
            String Para = arg;
            if(Para.length()>3)
            {
            }

	    if( !showDialog() )
                return DONE;

	    return DOES_8G+DOES_STACKS+SUPPORTS_MASKING;
	}

	/** Displays a modal options dialog. */
	public boolean showDialog()
        {
		GenericDialog gd = new GenericDialog("Analyze Particles");

                minSize = staticMinSize;
		maxSize = staticMaxSize;
                MinAxisLen = staticMinAxisLen;
                LowerThr = staticLowerThr;

                gd.addNumericField("Minimum Size (pixels):", minSize, 0);
		gd.addNumericField("Maximum Size (pixels):", maxSize, 0);
		gd.addNumericField("Min Axis Length:", MinAxisLen, 2);
		gd.addNumericField("Lower Threshold:", LowerThr, 0);


                gd.showDialog();
		if (gd.wasCanceled())
			return false;

		minSize = (int)gd.getNextNumber();
		maxSize = (int)gd.getNextNumber();
                MinAxisLen = (double)gd.getNextNumber();
                LowerThr = (int)gd.getNextNumber();

		if (gd.invalidNumber())
                {
			IJ.error("Minimum Size, Maximum Size or Bins invalid.");
			return false;
		}
		staticMinSize = minSize;
		staticMaxSize = maxSize;
                staticMinAxisLen = MinAxisLen;
                staticLowerThr = LowerThr;

		return true;
	}

        boolean Initialize(ImageProcessor ip)
        {
            //The height and width of the image
            m_nHeight = ip.getHeight();
            m_nWidth = ip.getWidth();

            if(m_nHeight<=0 || m_nWidth<=0)
            {
		IJ.showMessage("Particle Analyzer",
        			"Invalid Image size." );
                return false;
            }

            m_nParticleRec = new RuffleParticle[m_nWidth][m_nHeight];

            for(int i=0; i<m_nWidth; i++)
            {
                for(int j=0; j<m_nHeight; j++)
                {
                    m_nParticleRec[i][j] = null;
                }
            }

            //Allocate the space for the particles
            m_lstParticle = new RuffleParticle [m_nWidth*m_nHeight];
            for(int i=0; i<m_nWidth*m_nHeight; i++)
                m_lstParticle[i] = null;

            m_dMinThr = staticLowerThr;
            m_dMaxThr = 255;

            return true;
        }


	public void run(ImageProcessor ip)
        {

            if( !Initialize( ip ) )
                return;

            RunTag(ip);

            if( SetupStack() )
                JudgeRuffle(ip);
	}


        int RunTag(ImageProcessor ip)
        {
            int nLstPos = 0;
            for(int i=0; i<m_nWidth; i++)
            {
               //

                for(int j=0; j<m_nHeight; j++)
                {
                    int nBrightness = ip.getPixel(i, j);
                    if( nBrightness<m_dMinThr )
                        continue;

                    if( m_nParticleRec[i][j]!=null )
                        continue;

                    RuffleParticle particle = new RuffleParticle();

                    TagParticle(i, j, ip, particle);

                    m_lstParticle[nLstPos] = particle;
                    nLstPos++;
                }
            }

            return 0;
        }


        int TagParticle(int nPosX, int nPosY, ImageProcessor ip, 
RuffleParticle particle)
        {
            if(nPosX<0 || nPosX>=m_nWidth  ||
               nPosY<0 || nPosY>=m_nHeight)
            {
                return 1;
            }

            int nBrightness = ip.getPixel(nPosX, nPosY);
            if( nBrightness<m_dMinThr )
                return 2;

            boolean[][] nFoundFlag = new boolean [m_nWidth][m_nHeight];
            for(int u=0; u<m_nWidth; u++)
               for(int v=0; v<m_nHeight; v++)
                  nFoundFlag[u][v] = false;

            int nNewPosX, nNewPosY;
            int nStartX = nPosX;
            int nStartY = nPosY;
            nFoundFlag[nStartX][nStartY] = true;
            boolean bSearch = true;

            while(bSearch)
            {
                //Right
                nNewPosX = nStartX;
                nNewPosY = nStartY;
                for(;nNewPosX<m_nWidth; nNewPosX++)
                {
                    nBrightness = ip.getPixel(nNewPosX, nNewPosY);
                    if( nBrightness >= m_dMinThr )
                    {
                        m_nParticleRec[nNewPosX][nNewPosY] = particle;
                    }
                    else
                    {
                        break;
                    }
                }

                //Left
                nNewPosX = nStartX;
                nNewPosY = nStartY;
                for(;nNewPosX>=0; nNewPosX--)
                {
                    nBrightness = ip.getPixel(nNewPosX, nNewPosY);
                    if( nBrightness >= m_dMinThr )
                    {
                        m_nParticleRec[nNewPosX][nNewPosY] = particle;
                    }
                    else
                    {
                        break;
                    }
                }

                //Upper
                nNewPosX = nStartX;
                nNewPosY = nStartY;
                for(;nNewPosY<m_nHeight; nNewPosY++)
                {
                    nBrightness = ip.getPixel(nNewPosX, nNewPosY);
                    if( nBrightness >= m_dMinThr )
                    {
                        m_nParticleRec[nNewPosX][nNewPosY] = particle;
                    }
                    else
                    {
                        break;
                    }
                }

                //Lower
                nNewPosX = nStartX;
                nNewPosY = nStartY;
                for(;nNewPosY>=0; nNewPosY--)
                {
                    nBrightness = ip.getPixel(nNewPosX, nNewPosY);
                    if( nBrightness >= m_dMinThr )
                    {
                        m_nParticleRec[nNewPosX][nNewPosY] = particle;
                    }
                    else
                    {
                        break;
                    }
                }

                //Right-Upper
                nNewPosX = nStartX;
                nNewPosY = nStartY;
                for(;nNewPosX<m_nWidth; nNewPosX++)
                {
                    nBrightness = ip.getPixel(nNewPosX, nNewPosY);
                    if( nBrightness >= m_dMinThr )
                    {
                        m_nParticleRec[nNewPosX][nNewPosY] = particle;
                    }
                    else
                    {
                        break;
                    }

                    nNewPosY--;
                    if(nNewPosY<0)
                        break;
                }

                //Left-Upper
                nNewPosX = nStartX;
                nNewPosY = nStartY;
                for(;nNewPosX>=0; nNewPosX--)
                {
                    nBrightness = ip.getPixel(nNewPosX, nNewPosY);
                    if( nBrightness >= m_dMinThr )
                    {
                        m_nParticleRec[nNewPosX][nNewPosY] = particle;
                    }
                    else
                    {
                        break;
                    }

                    nNewPosY--;
                    if(nNewPosY<0)
                        break;
                }

                //Right-Lower
                nNewPosX = nStartX;
                nNewPosY = nStartY;
                for(;nNewPosX<m_nWidth; nNewPosX++)
                {
                    nBrightness = ip.getPixel(nNewPosX, nNewPosY);
                    if( nBrightness >= m_dMinThr )
                    {
                        m_nParticleRec[nNewPosX][nNewPosY] = particle;
                    }
                    else
                    {
                        break;
                    }

                    nNewPosY++;
                    if(nNewPosY>=m_nHeight)
                        break;
                }

                //Left-Lower
                nNewPosX = nStartX;
                nNewPosY = nStartY;
                for(;nNewPosX>=0; nNewPosX--)
                {
                    nBrightness = ip.getPixel(nNewPosX, nNewPosY);
                    if( nBrightness >= m_dMinThr )
                    {
                        m_nParticleRec[nNewPosX][nNewPosY] = particle;
                    }
                    else
                    {
                        break;
                    }

                    nNewPosY++;
                    if(nNewPosY>=m_nHeight)
                        break;
                }

                //Search for the point
                boolean bFound = false;
                for(int i=0; i<m_nWidth; i++)
                {
                    for(int j=0; j<m_nHeight; j++)
                    {
                        if( m_nParticleRec[i][j]== particle && 
nFoundFlag[i][j]==false )
                        {
                            nStartX = i;
                            nStartY = j;
                            nFoundFlag[i][j] = true;
                            bFound = true;
                            break;
                        }
                    }
                    if(bFound==true)
                        break;
                }

                if(bFound==false)
                    bSearch = false;
            }

            return 0;
        }

    boolean SetupStack()
    {
        byte[] rLUT = new byte[256];
        byte[] gLUT = new byte[256];
	byte[] bLUT = new byte[256];
	for(int i=0; i<256; i++)
        {
            rLUT[i]=(byte)i;
            gLUT[i]=(byte)i;
            bLUT[i]=(byte)i;
	}

	ColorModel cm = new IndexColorModel(8, 256, rLUT, gLUT, bLUT);
        StackParticle = new ImageStack(m_nWidth, m_nHeight, cm);

        boolean bSliceExist = false;
        int nPos = 0;
        int nParticleCount = 0;
        byte[] pixels = new byte[m_nWidth*m_nHeight];
        while(m_lstParticle[nPos]!=null)
        {
            int nPixelCount = 0;
            int nOffset=0;
            for(int i=0; i<m_nHeight; i++)
            {
                for(int j=0; j<m_nWidth; j++)
                {
                    if( m_nParticleRec[j][i]== m_lstParticle[nPos] )
                    {
                        nPixelCount++;
                        pixels[nOffset++] = (byte)255;
                    }
                    else
                    {
                        pixels[nOffset++] = 0;
                    }
                }
            }
            m_lstParticle[nPos].m_dActualArea = nPixelCount;
            m_lstParticle[nPos].m_dActualPerimeter = 
EvalPerimerter(m_lstParticle[nPos]);

            if( nPixelCount>minSize )
            {
                if( nPixelCount>maxSize )
                {
        		nPos++;
                        continue;
                }

                byte[] tmp = new byte[m_nWidth*m_nHeight];
                System.arraycopy(pixels, 0, tmp, 0, m_nWidth*m_nHeight);
                StackParticle.addSlice("", tmp);
                bSliceExist = true;
            }

            //Update the loop controller
	    nPos++;
        }

        if(bSliceExist)
        {
            StackImp = new ImagePlus("Particle", StackParticle);
            StackImp.show();
        }
        else
        {
            byte[] tmp = new byte[m_nWidth*m_nHeight];
            StackParticle.addSlice("", tmp);
            StackImp = new ImagePlus("Particle", StackParticle);
            StackImp.show();
        }

        return bSliceExist;
    }

    int EvalPerimerter(RuffleParticle Particle)
    {
        int nPerimeter = 0;
        for(int i=0; i<m_nWidth; i++)
        {
            for(int j=0; j<m_nHeight; j++)
            {
                int nRes = IsBoundary( i, j, Particle);
                if(nRes==1)
                    nPerimeter++;
            }
        }
        Particle.m_dActualPerimeter = nPerimeter;

        return nPerimeter;
    }

    int EvalPerimerter(boolean[][] Img)
    {
        int nPerimeter = 0;
        for(int i=0; i<m_nWidth; i++)
        {
            for(int j=0; j<m_nHeight; j++)
            {
                int nRes = IsBoundary( i, j, Img);
                if(nRes==1)
                    nPerimeter++;
            }
        }

        return nPerimeter;
    }

    int IsBoundary(int nPosX, int nPosY, RuffleParticle SelParticle)
    {
            if(nPosX<0 || nPosX>=m_nWidth  ||
               nPosY<0 || nPosY>=m_nHeight)
            {
                return -1;
            }

            if( m_nParticleRec[nPosX][nPosY]!=SelParticle )
                return -1;

            int nNbX, nNbY;
            int nBndSide = 0;
            //Upper
            nNbX = nPosX;
            nNbY = nPosY-1;
            if(nNbX<0 || nNbX>=m_nWidth  ||
               nNbY<0 || nNbY>=m_nHeight)
            {
                return 1;
            }
            else
            {
                if( m_nParticleRec[nNbX][nNbY]!=SelParticle )
                    return 1;
            }

            //Lower
            nNbX = nPosX;
            nNbY = nPosY+1;
            if(nNbX<0 || nNbX>=m_nWidth  ||
               nNbY<0 || nNbY>=m_nHeight)
            {
                return 1;
            }
            else
            {
                if( m_nParticleRec[nNbX][nNbY]!=SelParticle )
                    return 1;
            }

            //Left
            nNbX = nPosX-1;
            nNbY = nPosY;
            if(nNbX<0 || nNbX>=m_nWidth  ||
               nNbY<0 || nNbY>=m_nHeight)
            {
                return 1;
            }
            else
            {
                if( m_nParticleRec[nNbX][nNbY]!=SelParticle )
                    return 1;
            }

            //Right
            nNbX = nPosX+1;
            nNbY = nPosY;
            if(nNbX<0 || nNbX>=m_nWidth  ||
               nNbY<0 || nNbY>=m_nHeight)
            {
                return 1;
            }
            else
            {
                if( m_nParticleRec[nNbX][nNbY]!=SelParticle )
                    return 1;
            }

            return 0;
    }

    void JudgeRuffle(ImageProcessor ip)
    {
        int nSlice = StackParticle.getSize();
        String OutFile = "c:\\result.txt";
        String Space = "\t";

        String ColName ="No." + Space + "Area" + Space
                       + "M1" + Space + "M2";

        double[] dResultArray = new double[100];

        try
        {
            PrintWriter q = new PrintWriter( new FileWriter(OutFile), true 
);
            q.println(ColName);

            for(int i=1; i<=nSlice; i++)
            {
                for(int j=0; j<99; j++)
                    dResultArray[j] = 0.0;

                byte[] pixels = (byte[])StackParticle.getPixels(i);
                if(pixels==null)
                    break;

                ProcessOneParticle(pixels, dResultArray, ip);

                String RowData = i-1 + Space +
                                 dResultArray[0] + Space +
                                 dResultArray[5] + Space +
                                 dResultArray[6];

                q.println(RowData);
            }
        }
        catch (Exception e)
        {
            e.printStackTrace();    // this will print the error to standard 
output
        }

    }

    void ProcessOneParticle(byte[] pixels, double[] dResultArray, 
ImageProcessor ip)
    {
        boolean[][] ImgOrg = new boolean [m_nWidth][m_nHeight];
        int nPixelCount = 0;
        int nPos=0;
        for(int i=0; i<m_nHeight; i++)
        {
            for(int j=0; j<m_nWidth; j++)
            {
                byte nValue = pixels[nPos++];
                if(nValue==(byte)255)
                {
                    ImgOrg[j][i] = true;
                    nPixelCount++;
                }
                else
                    ImgOrg[j][i] = false;
            }
        }

        double[] dSaved = new double[20];
        calc_features_1( ImgOrg,
                         0,0,
                         m_nWidth, m_nHeight, dSaved, ip);


        //Ellipse axis ratio
        double[] dAxisLen = new double[2];
        EvalEllipseAxis( ImgOrg,
                         m_nWidth, m_nHeight,
                         dSaved[3],
                         dAxisLen );


        //
        dResultArray[0] = dSaved[0];        //area
        dResultArray[1] = dSaved[1];        //xmean
        dResultArray[2] = dSaved[2];        //ymean
        dResultArray[3] = dSaved[3];        //oreintation
        dResultArray[4] = dSaved[4];        //Nt11
        dResultArray[5] = dSaved[5];        //fai1
        dResultArray[6] = dSaved[6];        //fai2


        //Condition 1: size requirement
        if( nPixelCount>maxSize || nPixelCount<minSize )
        {
            //Clear image
            int nOffset=0;
            for(int i=0; i<m_nHeight; i++)
            {
                for(int j=0; j<m_nWidth; j++)
               {
                   pixels[nOffset++] = 0;
               }
            }
            return;
        }



        //Condition2: Long Axis requirement
        if(
            dAxisLen[0]*2 <= MinAxisLen)
        {
              //Clear image
              int nOffset=0;
              for(int i=0; i<m_nHeight; i++)
              {
                  for(int j=0; j<m_nWidth; j++)
                 {
                     pixels[nOffset++] = 0;
                 }
              }
              return;
        }


        return;
   }

    void GetConnectInWidth(int nW, int nStart, int[] nEnd, boolean[][] img)
    {
        nEnd[0] = nStart;
        boolean bFlag = img[nW][nStart];
        for(int i=nStart+1; i<m_nHeight; i++)
        {
            if( img[nW][i]==bFlag )
            {
                nEnd[0] = i;
            }
            else
                break;
        }
    }

    void GetConnectInHeight(int nH, int nStart, int[] nEnd, boolean[][] img)
    {
        nEnd[0] = nStart;
        boolean bFlag = img[nStart][nH];
        for(int i=nStart+1; i<m_nWidth; i++)
        {
            if( img[i][nH]==bFlag )
            {
                nEnd[0] = i;
            }
            else
                break;
        }
    }

    int IsBoundary(int nPosX, int nPosY, boolean[][] Img)
    {
            if(nPosX<0 || nPosX>=m_nWidth  ||
               nPosY<0 || nPosY>=m_nHeight)
            {
                return -1;
            }

            if( Img[nPosX][nPosY]==false )
                return -1;

            int nNbX, nNbY;
            int nBndSide = 0;
            //Upper
            nNbX = nPosX;
            nNbY = nPosY-1;
            if(nNbX<0 || nNbX>=m_nWidth  ||
               nNbY<0 || nNbY>=m_nHeight)
            {
                return 1;
            }
            else
            {
                if( Img[nNbX][nNbY]==false )
                    return 1;
            }

            //Lower
            nNbX = nPosX;
            nNbY = nPosY+1;
            if(nNbX<0 || nNbX>=m_nWidth  ||
               nNbY<0 || nNbY>=m_nHeight)
            {
                return 1;
            }
            else
            {
                if( Img[nNbX][nNbY]==false )
                    return 1;
            }

            //Left
            nNbX = nPosX-1;
            nNbY = nPosY;
            if(nNbX<0 || nNbX>=m_nWidth  ||
               nNbY<0 || nNbY>=m_nHeight)
            {
                return 1;
            }
            else
            {
                if( Img[nNbX][nNbY]==false )
                    return 1;
            }

            //Right
            nNbX = nPosX+1;
            nNbY = nPosY;
            if(nNbX<0 || nNbX>=m_nWidth  ||
               nNbY<0 || nNbY>=m_nHeight)
            {
                return 1;
            }
            else
            {
                if( Img[nNbX][nNbY]==false )
                    return 1;
            }

            return 0;
    }

    double find_moment( boolean[][] ImageX,
                        int p, int q,
                        int N1, int M1,
                        int N2, int M2,
                        double _x, double _y,
                        ImageProcessor ip, double[] dEvalArea)
    {
      int i, j, k;
      double a=0, h;

      dEvalArea[0] = 0.0;
      for( i=N1, a=0.; i<N2; i++ )
        for( j=M1; j<M2; j++ )
        {
          if(ImageX[i][j]==true)
          {
            //temp h = 1.0;
            dEvalArea[0] = dEvalArea[0] + 1.0;
            h = ip.getPixel(i,j);
            for( k=0; k<p; k++)
                h*=(i-_x);

            for(k=0; k<q; k++)
                h*=(-j-_y);

            a+=h;
          }
        }
      return (a);
    }


int calc_features_1( boolean[][] ImageX,
                     int N1, int M1,
                     int N2, int M2,
                     double[] dResult, ImageProcessor ip)
{
      double area, xmean, ymean, orientation, eccentricity1, eccentricity2, 
spread;
      double m00, m10, m01, mi11, mi20, mi02, th, nom, denom;
      double mi03, mi30, mi12, mi21;

      double[] dEvalArea = new double [1];

      m00 = find_moment(ImageX, 0, 0, N1, M1, N2, M2, 0., 0., ip, 
dEvalArea);
      m10 = find_moment(ImageX, 1, 0, N1, M1, N2, M2, 0., 0., ip, 
dEvalArea);
      m01 = find_moment(ImageX, 0, 1, N1, M1, N2, M2, 0., 0., ip, 
dEvalArea);

      area=m00;
      xmean=m10/m00;
      ymean=m01/m00;

      mi11 = find_moment(ImageX, 1, 1, N1, M1, N2, M2, xmean, ymean, ip, 
dEvalArea);
      mi20 = find_moment(ImageX, 2, 0, N1, M1, N2, M2, xmean, ymean, ip, 
dEvalArea);
      mi02 = find_moment(ImageX, 0, 2, N1, M1, N2, M2, xmean, ymean, ip, 
dEvalArea);
      mi03 = find_moment(ImageX, 0, 3, N1, M1, N2, M2, xmean, ymean, ip, 
dEvalArea);
      mi12 = find_moment(ImageX, 1, 2, N1, M1, N2, M2, xmean, ymean, ip, 
dEvalArea);
      mi21 = find_moment(ImageX, 2, 1, N1, M1, N2, M2, xmean, ymean, ip, 
dEvalArea);
      mi30 = find_moment(ImageX, 3, 0, N1, M1, N2, M2, xmean, ymean, ip, 
dEvalArea);

      th = (mi11!=0 || mi02-mi20!=0) ? (1./2)*Math.atan2(2*mi11, mi20-mi02) 
: 0;
      nom = mi02 * Math.cos(th) * Math.cos(th)
            + mi20 * Math.sin(th) *  Math.sin(th)
            - mi11 * Math.sin(2*th);
      denom = mi02 * Math.sin(th) * Math.sin(th)
              + mi20 * Math.cos(th) * Math.cos(th)
              + mi11 * Math.cos(2*th);
      orientation = th;

      //normalized central moment
      double nt11 = mi11/(m00*m00);
      double nt20 = mi20/(m00*m00);
      double nt02 = mi02/(m00*m00);
      double nt03 = mi03/( Math.pow(m00, 2.5) );
      double nt30 = mi30/( Math.pow(m00, 2.5) );
      double nt12 = mi12/( Math.pow(m00, 2.5) );
      double nt21 = mi21/( Math.pow(m00, 2.5) );;

      double fai1 = nt20 + nt02;
      double fai2 = (nt20-nt02)*(nt20-nt02)+4.0*nt11*nt11;

      dResult[0] = dEvalArea[0];
      dResult[1] = xmean;
      dResult[2] = ymean;
      dResult[3] = orientation;
      dResult[4] = nt11;

      dResult[5] = fai1;
      dResult[6] = fai2;


      return(0);
}

int EvalEllipseAxis( boolean[][] ImageX,
                     int Width, int Height,
                     double dAngle,
                     double[] dResult)
{
    boolean bFirstEnter = true;
    double dMinX=0.0, dMaxX=0.0;
    double dMinY=0.0, dMaxY=0.0;

    for(int i=0; i<Width; i++)
    {
        for(int j=0; j<Height; j++)
        {
            if(!ImageX[i][j])
                continue;
            double x0 = i;
            double y0 = -j;

            double dAngleMod=0.0;
            if(dAngle>0.0)
                dAngleMod = dAngle;
            else
                dAngleMod = dAngle + Math.PI*2.0;

            double r = Math.sqrt(x0*x0 + y0*y0);
            double xr =    x0 * Math.cos(dAngleMod)   + y0 * 
Math.sin(dAngleMod);
            double yr = -( x0 * Math.sin(dAngleMod) ) + y0 * 
Math.cos(dAngleMod);

            if(bFirstEnter)
            {
                bFirstEnter = false;
                dMinX = dMaxX = xr;
                dMinY = dMaxY = yr;
            }

            if(dMinX>xr)
                dMinX = xr;
            if(dMaxX<xr)
                dMaxX = xr;
            if(dMinY>yr)
                dMinY = yr;
            if(dMaxY<yr)
                dMaxY = yr;
        }
    }

    dResult[0] = (dMaxX-dMinX)*0.5;
    dResult[1] = (dMaxY-dMinY)*0.5;

    return(0);
}
}