zernike.java

求图形Zernike矩的java程序

package ImageTool;

/**
 * Title:   zernike.java
 * Description:  My collection of java image processing tools.
 * Copyright:    Copyright (c) 2001
 * Company:      Columbia Univ.
 * @author Lijun Tang
 * @version 1.0
 */
import java.net.*;
import java.io.*;
import java.awt.*;
import java.util.*;
import java.awt.image.*;


import java.awt.event.*;
import javax.swing.*;

import java.awt.geom.AffineTransform;
import java.awt.font.TextLayout;

import java.awt.image.*;
import com.sun.image.codec.jpeg.*;

public class zernike implements ImgToolInterface{

  public static final int  NULL=0;
  public static final double M_PI=3.14159265358979323846;
  public static final double ZER_INFINITY=1.223e10;
  public static final double ZER_TINY=1e-6;

  public static final int ZER_NO_MEMORY=-10;
  public static final int ZER_ERROR=-1;
  public static final int ZER_OK=0;
//  char[] IBAND;

  public zernike() {
  }
/**
 * zer_pol_R() compute the Rnm(p) in polarnomial definition of V(n,m,x,y)
 * Definition of Rnm(p) refer to &[1]
 * @return res[0] as the value of Rnm(p)
 */
  public static int zer_pol_R(int n, int m_in, double x, double y, double[] res)
  {
  int i;
  int m;
  int s;
  int a; /* (n-s)! */
  int b; /*   s!   */
  int c; /* [(n+|m|)/2-s]! */
  int d; /* [(n-|m|)/2-s]! */
  int sign;

  m = Math.abs(m_in);

  if ((n-m)%2!=0)
    System.out.print(
	    "zer_pol_R: Warning. R(%i,%i,%5.2f,%5.2f). n-|m| = %i is odd!\n"+n+m+x+y+m);

  /* The code is optimized with respect to the faculty operations */

  res[0] = 0.0;
  if ((x*x + y*y) <= 1.0)
    {
      sign = 1;
      a = 1;
      for (i=2; i<=n; i++)
	a*=i;
      b=1;
      c = 1;
      for (i=2; i <= (n+m)/2; i++)
	c*=i;
      d = 1;
      for (i=2; i <= (n-m)/2; i++)
	d*=i;

      /* Before the loop is entered, all the integer variables (sign, a, */
      /* b, c, d) have their correct values for the s=0 case. */
      for (s=0; s<= (n-m)/2; s++)
	{
	  /*printf("zer_pol_R: s=%i, n=%i, m=%i, x=%6.3f, y=%6.3f, a=%i, */
	  /*b=%i, c=%i, d=%i, sign=%i\n", s,n,m,x,y,a,b,c,d,sign); */
	  res[0] += sign * (a*1.0/(b*c*d)) * Math.pow((x*x + y*y),(n/2.0)-s);

	  /* Now update the integer variables before the next iteration of */
	  /* the loop. */

	  if (s < (n-m)/2)
	    {
	      sign = -sign;
	      a = a/(n-s);
	      b = b*(s+1);
	      c = c / ((n+m)/2 - s);
	      d = d / ((n-m)/2 - s);
	    }
	}
    }
  return ZER_OK;
}

/**
 *  zer_pol() computes the zernike basis function
 *    V(n,m,x,y).
 * @return res[1] is the dcomplex for V(n,m,x,y)
*/
public static int zer_pol(int n, int m, double x, double y, DCOMPLEX[] res)
{
  double[] R=new double[1];
  double arg;

  if ((x*x + y*y) > 1.0)
    {
      res[0].re = 0.0;
      res[0].im = 0.0;
    }
  else
    {
      zer_pol_R(n,m,x,y, R);
      arg = m*Math.atan2(y,x);
      res[0].re = R[0]*Math.cos(arg);
      res[0].im = R[0]*Math.sin(arg);
    }

  return ZER_OK;
}



/*F:zer_con=zer_mom*/
/*F:zer_rec=zer_mom*/
/*F:zer_pol=zer_mom*/
/*F:zer_mom*

________________________________________________________________

		zer_mom
________________________________________________________________

Name:		zer_mom, zer_pol, zer_rec, zer_con - Zernike moments

Syntax:		| #include <xite/zernike.h>
                |
                | int zer_mom( BufferedImage inband, int n, int m,
		|    DCOMPLEX *res );
		|
		| int zer_pol( int n, int m, double x, double y,
		|    DCOMPLEX *res );
		|
		| int zer_con( int n, BufferedImage inband, BufferedImage outband );
		|
		| int zer_rec( int order, BufferedImage inband, BufferedImage outband );
		|

Description:    The Zernike moment of order 'n' and repetition 'm' of
                an image f(x,y) is defined as follows:
		|          n+1                             *
                | A(n,m) = ---- Sum Sum f(x,y)[V(n,m, x,y)]
		|           pi   x   y
		|
		| where x^2+y^2 <= 1
		|
		The image V(n,m, x,y) is the Zernike basis images of
		order 'n' and repetition 'm'. These basis images are
		complex and orthogonal. The Zernike moments are
		essentially the projections of the input image onto
		these basis images.

		The original image can be reconstructed from the
		Zernike moments. The N-th order approximation is given
		by
		|  ^         N
		|  f(x,y) = Sum Sum A(n,m) V(n,m, x,y)
		|           n=0  m
		|
		The contribution or information content of the n-th
		order moments is
		|
		|  I(x,y, n) = Sum A(n,m) V(n,m, x,y)
		|              m

		'zer_mom' computes A(n,m).

		'zer_pol' computes the zernike basis function
		V(n,m,x,y).

		'zer_con' computes the absolute	value of the contribution
		of the n-th order moments, i.e. the absolute value of
		I(x,y, n).

		'zer_rec' demonstrates how the input image can be
		reconstructed from the Zernike basis functions and the
		Zernike moments.

Return value:   | 0 : OK.

Author:		豬vind Due Trier, Dept. Informatics, Univ. Oslo
________________________________________________________________

*/
/**
 * 	'zer_mom' computes A(n,m).
 */
public static int zer_mom(BufferedImage inband, int n, int m, DCOMPLEX[] res)
{
  int i,j;
  int i_0, j_0;
  double i_scale, j_scale;
  double x,y;
  DCOMPLEX[] v=new DCOMPLEX[1];
  v[0]=new DCOMPLEX();
  int isize, jsize;
//  double[] m2=ImageOps.mmt(inband);

  if (inband.getType()!= BufferedImage.TYPE_BYTE_GRAY)
    System.out.print("zer_mom: Input pixel type must be unsigned byte.\n");

  if ((n<0) || (Math.abs(m) > n) || ((n-Math.abs(m))%2!=0))
    System.out.print("zer_mom: n=%i, m=%i, n-|m|=%i\n"+ n+ m+(n-Math.abs(m)));

  isize = inband.getWidth();
  jsize = inband.getHeight();

  i_0 = (isize+1)/2;
  j_0 = (jsize+1)/2;
  i_scale = isize/2.0;
  j_scale = jsize/2.0;

//  System.out.print("("+i_0+","+j_0+")");

  res[0].re = 0.0;
  res[0].im = 0.0;
  int[] src_1d=new int[isize*jsize];

  Raster rt=inband.getData();
  rt.getSamples(0,0,isize,jsize,0,src_1d);
//  getRGB(0,0,isize,jsize,src_1d,0,1);

  for (i=1; i<=isize; i++){
    for (j=1; j<= jsize; j++)
      {
	x = (i-i_0)/i_scale;
	y = (j-j_0)/j_scale;

/*	printf ("zer_mom: x=%6.3f y=%6.3f i=%i j=%i\n", x,y,i,j);
*/
//        int value=inband.getData().getSample(j ,i,0);
//        int value=src_1d[(i-1)*jsize+j-1]==BKGRND?0:1;
        if(src_1d[(i-1)*jsize+j-1]==BKGRND) continue;
	if (((x*x + y*y) <= 1.0))
	  {
	    zer_pol(n,m,x,y, v);
	    res[0].re += v[0].re;
	    res[0].im += (-v[0].im);
	  }
      }
//      System.out.print("i="+i+"\n");
      }
  res[0].re = res[0].re*(n+1)/M_PI;
  res[0].im = res[0].im*(n+1)/M_PI;

  return ZER_OK;

} /* END zer_mom */

/**
 * 	'zer_mom' computes A(n,m).
 */
public static int zer_mom(int[] src_1d, int ww,int hh,Point centroid, int n, int m, DCOMPLEX[] res)
{
  int i,j;
  int i_0, j_0;
  double i_scale, j_scale;
  double x,y;
  DCOMPLEX[] v=new DCOMPLEX[1];
  v[0]=new DCOMPLEX();
  int isize, jsize;
//  double[] m2=ImageOps.mmt(inband);

//  if (inband.getType()!= BufferedImage.TYPE_BYTE_GRAY)
//    System.out.print("zer_mom: Input pixel type must be unsigned byte.\n");

  if ((n<0) || (Math.abs(m) > n) || ((n-Math.abs(m))%2!=0))
    System.out.print("zer_mom: n=%i, m=%i, n-|m|=%i\n"+ n+ m+(n-Math.abs(m)));

  isize = ww;//inband.getWidth();
  jsize = hh;//inband.getHeight();

  i_0 = centroid.x;//(isize+1)/2;
  j_0 = centroid.y;//(jsize+1)/2;
//  i_scale = isize/2.0;
//  j_scale = jsize/2.0;

  int radius=ww-i_0;
  if(radius<i_0) radius=i_0;
  i_scale=Math.sqrt(2)*radius;
  radius=hh-j_0;
  if(radius<j_0) radius=j_0;
  j_scale=Math.sqrt(2)*radius; //note we want to construct a circle to contain the rectangle


  res[0].re = 0.0;
  res[0].im = 0.0;
//  int[] src_1d=new int[isize*jsize];
//
//  Raster rt=inband.getData();
//  rt.getSamples(0,0,isize,jsize,0,src_1d);
//  getRGB(0,0,isize,jsize,src_1d,0,1);

  for (i=0; i<isize; i++){
    for (j=0; j<jsize; j++)
      {
	x = (i-i_0)/i_scale;
	y = (j-j_0)/j_scale;

/*	printf ("zer_mom: x=%6.3f y=%6.3f i=%i j=%i\n", x,y,i,j);
*/
//        int value=inband.getData().getSample(j ,i,0);
        int value=src_1d[i*jsize+j];
        if(value<=0) continue;

        value=1;  //we regard all images as binary

	if (((x*x + y*y) <= 1.0))
	  {
	    zer_pol(n,m,x,y, v);
	    res[0].re += value*v[0].re;
	    res[0].im += value*(-v[0].im);
	  }
      }
//      System.out.print("i="+i+"\n");
      }
  res[0].re = res[0].re*(n+1)/M_PI;
  res[0].im = res[0].im*(n+1)/M_PI;

  return ZER_OK;

} /* END zer_mom */
/**
 * 		'zer_con' computes the absolute	value of the contribution
 *		of the n-th order moments, i.e. the absolute value of
 *		I(x,y, n).
 */
public static int zer_con(int n, BufferedImage inband, BufferedImage outband)
{
  DCOMPLEX[][] transformed;
  int m;
  int i, j;
  double i_0, j_0;
  double i_scale, j_scale;
  int isize, jsize;
  double x, y;
  double min, max, scale;
  DCOMPLEX[] v=new DCOMPLEX[1];
  v[0]=new DCOMPLEX();
  DCOMPLEX[] a=new DCOMPLEX[1];
  a[0]=new DCOMPLEX();

  if (inband.getType()!= BufferedImage.TYPE_BYTE_GRAY)// Iu_byte_typ)
    System.out.print( "zer_con: Input pixel type must be unsigned byte.\n");
  if (outband.getType() != BufferedImage.TYPE_BYTE_GRAY)
    System.out.print("zer_con: Output pixel type must be unsigned byte.\n");

  isize = inband.getWidth();
  jsize = inband.getHeight();

  i_0 = (isize+1)/2.0;
  j_0 = (jsize+1)/2.0;
  i_scale = isize/2.0;
  j_scale = jsize/2.0;

  transformed = new DCOMPLEX[isize+1][jsize+1];//(IDC_BAND) Imake_band(Id_complex_typ, isize, jsize);


  for (i=1; i<=isize; i++)
    for (j=1; j<=jsize; j++)
      {
        transformed[i][j]=new DCOMPLEX();
	transformed[j][i].re = 0.0;
	transformed[j][i].im = 0.0;
      }
  for (m=-n; m<=n; m++)
    {
      if ((n-Math.abs(m)) %2 == 0)
	{
	  zer_mom(inband, n,m, a);
	  for (i=1; i<=isize; i++)
	    for (j=1; j<=jsize; j++)
	      {
		x = (i-i_0)/i_scale;
		y = (j-j_0)/j_scale;
/*		if (x*x + y*y <= 1.0) */
/* 		printf ("zer_con: x=%6.3f y=%6.3f i=%i j=%i\n", x,y,i,j);
*/
		zer_pol(n,m,x,y, v);
		transformed[j][i].re += a[0].re*v[0].re - a[0].im*v[0].im;
		transformed[j][i].im += a[0].re*v[0].im + a[0].im*v[0].re;
/*		printf ("zer_con: x,y=(%6.3f,%6.3f), contr=(%9.6f,%9.6f),",
			x,y,transformed[j][i].re,
			transformed[j][i].im);*/
/*		printf (" a=(%9.6f,%9.6f), v=(%9.6f,%9.6f)\n",
			a.re, a.im, v.re, v.im);*/
	      }
	}
    }
  for (i=1; i<=isize; i++)
    for (j=1; j<=jsize; j++)
      {
	transformed[j][i].re = Math.sqrt
	  (  transformed[j][i].re*transformed[j][i].re
	   + transformed[j][i].im*transformed[j][i].im);
	transformed[j][i].im = 0.0;
      }
  min = max = transformed[1][1].re;
  for (i=1; i<=isize; i++)
    for (j=1; j<=jsize; j++)
      {
	if (transformed[j][i].re > max)
	  max = transformed[j][i].re;
	if  (transformed[j][i].re < min)
	  min = transformed[j][i].re;
      }
  if (max == min)
    scale = 0.0;
  else
    scale = 255.0/(max-min);

  System.out.print("zer_con: min=%9.6f, max=%9.6f, sacle=%f\n"+ min+ max+ scale);

  int[] dest=new int[isize*jsize];

  for (i=1; i<=isize; i++)
    for (j=1; j<=jsize; j++)
      {
        dest[(i-1)*jsize+j-1]=(int)(transformed[j][i].re * scale);
//	outband.getRaster().setSample(j,i,0,(byte) (transformed[j][i].re * scale));
/*	printf ("zer_con: outband[%3i][%3i]=%3i\n", j, i, outband[j][i]);*/
      }
    outband.getRaster().setSamples(0,0,isize,jsize,0,dest);

  return ZER_OK;
} /* END zer_con */
/**
 * 		'zer_rec' demonstrates how the input image can be
 *		reconstructed from the Zernike basis functions and the
 *		Zernike moments.
 */
public static BufferedImage zer_rec(int order, BufferedImage inband)
{
  DCOMPLEX[][] transformed;
  int n, m;
  int i, j;
  double i_0, j_0;
  double i_scale, j_scale;
  int isize, jsize;
  double x, y;
  DCOMPLEX[] a=new DCOMPLEX[1];
  a[0]=new DCOMPLEX();
  DCOMPLEX[] v=new DCOMPLEX[1];
  v[0]=new DCOMPLEX();

  BufferedImage outband=new BufferedImage(inband.getWidth(),inband.getHeight(),BufferedImage.TYPE_BYTE_GRAY);