/********************************************************************
Copyright 2001--200?, Bin Han, University of Alberta, June 18, 2001
All Rights Reserved

The use of this software is permitted for non-commercial, educational,
and research use only. The software and/or related materials are
provided "as-is" without warranty of any kind including any warranties
of performance or merchantability or fitness for a particular use or
purpose or for any purpose whatsoever, for the licensed product,
however used. In no event shall University of Alberta and/or Bin Han
be liable for any damages and/or costs, including but not
limited to incidental or consequential damages of any kind, including
economic damage or injury to property and lost profits, regardless of
whether University of Alberta shall be advised, have reason to know,
or in fact shall know of the possibility. User bears all risk relating
to quality and performance of the software and/or related materials.

Any use other than non-commercial, educational, or research, or any
redistribution in original or modified form requires prior written
authorization from the copyright holder.

Report errors, mistakes, and bugs to Bin Han at bhan@ualberta.ca
Send comments/suggestions to Bin Han at bhan@ualberta.ca

=============Documentation of the program=======================
This program is written in C while calling matlab engine function eig
to compute eigenvalues of matrices. The program is based on the method
described in the paper:
Bin Han, Computing the smoothness exponent of a symmetric
multivariate refinable function, preprint, (2001)

You can download the above paper at
http://www.ualberta.ca/~bhan/publ.htm

Though this program is not optimized from point of view of computation,
in general, it works fairly efficiently.

This program can be easily translated into matlab routine.
However, currently the matlab version of this program is NOT provided
yet.

====== Instruction for using this C program =====

Purpose:
      To compute the Sobolev smoothness (and also Holder smoothness when
      the symbol of the mask is nonnegative) for a two-dimensional
      refinable function with a finitely supported mask and a dilation
      matrix.

Put your 2D mask in a file using the following input format
(lines between -------).
-----------input mask as follows:------------
number_of_row   number_of_column   order_of_sum_rules
i  j Mask[i][j]
.
.
.
-1
----------end of mask-------
where min(i,j) >= 0 and only nonzero coefficients in a mask
are needed in the above format. So, the mask is supported on
[0, number_of_row-1][0, number_of_column-1] and satisfies the
sum rules of order order_of_sum_rules

For example, the mask for the Haar refinable function is put into
2  2  1
0  0  1.0
0  1  1.0
1  0  1.0
1  1  1.0
-1

---Since the mask for the Haar refinable function is supported on
[0,1]^2, so number_of_row=2, number_of_column=2, and since its mask
satisfies the sum rules of order 1, so order_of_sum_rule=1


You need to specify the following variables.

1. set the dilation matrix by modifying the variable dilation[2][2]

2. set the type of symmetry:
   FullAxesSym, HexagonalSym, AxesSym, XYSYM, OriginSym.
   by defining ONLY one of them

3. To compute Holder smoothness, define Holder
   otherwise undefine Holder to compute the Sobolev smoothness

4. IsIdentityDLT, is the dilation matrix a multiple of the identity matrix?
   This setting is used to find the invariant space quickly.
   So its only purpose to increase the efficiency of the program.
   You may leave this as undefined.

!!!!!!!!!!!!!!! Places that you must change in the program !!!!!
To find the places where necessary changes are needed, check the
comments with ??? sign

Explanation of types of symmetry:
   FullAxesSym: x=0, y=0, y=x, y=-x
   AxesSym: x=0, y=0
   XYSym: y=x, y=-x
   HexagonalSym: y=x,y=-x, plus 60 degree rotation, total 12 elements
   OriginSym: a(-i,-j)=a(i,j), symmetric about the origin, when the mask
   is not symmetric, choose this one.

   When the mask is not symmetric, you can define OriginSym since
   the mask after convolution is always symmetric about the origin. 

Condition assumed (!!!This version of program will not check it!!!):
   1. All the eigenvalues fo dilation must have the same modulus.
      That is, the dilation matrix is isotropic.

   2. The symmetry group FullAxesSym or HexagonalSym should compatible
      with the dilation matrix. I.e., M E M^{-1} \in G, E\in G.
      Where G is the symmetry group used.

   3. The mask (for computeing Holder smoothness) or the convolution of
      the mask (for computing Sobolev smoothness) must be invariant
      under the symmetry group.
   
   4. Stability of the refinable function is needed to obtain
      the exact smoothness exponent. However, the stability in many
      case can be deduced from the computed result. See the following
      paper for more details.

   For more details about the requirement, see the paper:
   Bin Han, Computing the smoothness exponent of a symmetric
   multivariate refinable function, preprint, (2001)
   Dowload the paper at: http://www.ualberta.ca/~bhan/publ.htm

Remarks:
   If the refinable function is not stable, then the computed Sobolev
   smoothness is a lower bound of the actual Sobolev smoothness exponent
   of the refinable function.
   When the refinable function is stable, but the symbol of the mask is
   not nonnegative and you are using this program to compute the
   Holder smoothness exponent, then the computed Holder smoothness is
   an upper bound of the actual Holder smoothness exponent.

Name of this program: d2sm.c
    If you put your mask into a file called cmask,
    after you compile d2sm.c using the makefile, 
    Then you can use
         d2sm cmask
    or use
         d2sm cmask outputfile
    to put the output into the file outputfile

    to obtain the smoothness exponent of its refinable function.

The makefile of the program is as follows:
=======makefile for the program d2sm.c=====
CC = gcc
CFLAGS = -Wall -g

INCDIR =-I/usr/local/matlabr12/extern/include
LIBINC = -L/usr/local/matlabr12/extern/lib/glnx86

d2sm: d2sm.c
	$(CC) $(INCDIR) -o d2sm d2sm.c $(LIBINC) -leng -lmx -lm

=======end of makefile for the program d2sm.c====

NOTICE: !!! You may have to change the paths in the makefile. !!!
        This program has been tested under Mandrake Linux 8.0 and
        Matlab version 12.

Modified by Bin Han from several other exisiting programs (these
existing programs were developed by Bin Han during the period
1997--2001 for special cases) on May 24, 2001.
Tested by Bin Han on June 28, 2001 on Linux platform.
Detailed documentation was added by Bin Han on June 18, 2001.

Posted on Web June 28, 2001 at
     http://www.ualberta.ca/~bhan/publ.htm
Report errors, mistakes, and bugs to
     Bin Han at bhan@ualberta.ca
Send comments/suggestions to
     Bin Han at bhan@ualberta.ca
*********************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
/*??? Change the path for the right location of the file engine.h
      Also change the paths in the makefile correspondingly!!!   ???*/
#include "/usr/local/matlabr12/extern/include/engine.h"

#define MAX(x)  ((x)>(0) ? (x):(0))
#define MIN(x,y)  ((x)<(y) ? (x):(y))
#define MAXELM  9000
#define TOL 0.0001

/*??? The maximum number of eigvaules to be displayed,
  don't change it unless your output is too large ???*/
#define MaxNumEig 128

/*??? If the dilation matrix is mI_2, then use this can speed up the
  program since the invariant space can be easily determined.
  If you are not sure, then leave it as undefined here   ???*/
#undef IsIdentityDLT

/*??? when define Holder, compute the Holder smoothness, otherwise,
  compute the Sobolev smoothness ??? */
#undef Holder

/*The symmetry group which must be compatible with the dilation matrix*/
/*??? choose ONLY one to the following symmetry type!!! ??? */
#define HexagonalSym  /* 12 elements in the group, 30 degree rotation */
#undef  FullAxesSym   /* 8 elements in the group, x=0, y=0, y=x, y=-x */
#undef AxesSym        /* 4 elements in the group, x=0, y=0            */
#undef XYSym          /* 4 elements in the group, y=x, y=-x           */
#undef OriginSym      /* 2 elements in the group, a(-i,-j)=a(i,j)     */
                      /* If mask is not symmetric, define OriginSym   */

extern int abs();
extern double fabs();
extern double log();
extern double sqrt();

int compare(a, b)
     double *a, *b;
{
  return ( *a < *b ?  1 : -1 );
}

double **m_alloc(row, col)
int row, col;
{
int i;
double **a;

a=(double **)calloc( row, sizeof( double *));
if ( a==NULL ) {
	fprintf( stderr, " Could not allocate memery!\n" );
	exit(1);
}
for(i=0; i<row; ++i) {
	a[i]=(double *)calloc( col, sizeof( double ));
	if ( a[i]==NULL ) {
		fprintf( stderr, " Could not allocate memery!\n" );
		exit(1);
	} /* end of if */
} /* end of for */
return a;
} /* end of m_alloc */

int main(argc, argv)
int argc;
char *argv[];
{
  /*you need to set this dilation matrix which is isotropic and
    compatible with the symmetry you choose at the first part */
  
  /*??? YOU NEED TO SET YOUR DILATION MATRIX HERE!!!
    THIS IS THE LAST PLACE FOR MODIFICATION      ???*/
  int dilation[2][2]={2,0,0,2};
  int i, ii, iii,i1, j, jj, jjj,j1, k, kk, kkk, Edge;
  int Sx, Sy, SZ, Delta, Dim, index, Num, Listx[MAXELM], Listy[MAXELM];
  int Px[MAXELM], Py[MAXELM];
  int det=dilation[0][0]*dilation[1][1]-dilation[0][1]*dilation[1][0], absdet=abs(det);
  double sum, sm, *jsp, *MatDTPr, **ConMask, **Mask, **DT, special_eig[MaxNumEig];
  Engine *ep;
  mxArray *MatDT;
  FILE *outfp, *infp;
  
  if( argc==1){
    infp=fopen("cmask", "r");
    outfp=stdout;
  } else if( argc==2) {
    infp=fopen(argv[1],"r"); 
    outfp=stdout;
  } else if (argc==3){
    infp=fopen(argv[1],"r"); 
    outfp=fopen(argv[2], "w");
  } else{
    fprintf(stderr, "Wrong! Usage: d2sm maskfile (outputfile)! exit!\n");
    exit(0);
  }
  
  if ( !( ep=engOpen("\0")) ){
    fprintf(stderr, "\n Can't start MATLAB engine\n"); exit(-1);
  }
  
  fscanf(infp, "%d\n", &Sx);
  fscanf(infp, "%d\n", &Sy);
  fscanf(infp, "%d\n", &Delta);
  
  fprintf(outfp, "The dilation matrix is=[%d %d; %d %d]\n", dilation[0][0],
                            dilation[0][1],dilation[1][0], dilation[1][1]);
  Dim=0;
#ifdef Holder
  ConMask=m_alloc(Sx, Sy);
  Sx/=2;  Sy/=2; 
  for(ii=0;ii<9999999;++ii){
    fscanf(infp, "%d\n", &i); if ( i < 0 ) break; /* use -1 end input */
    fscanf(infp, "%d\n", &j);
    fscanf(infp, "%lf\n", &ConMask[i][j]);
    if( fabs(ConMask[i][j])>0.000000000001 ){
      Listx[Dim]=i-Sx;  Listy[Dim]=j-Sy;   ++Dim;
    }
  }
  fprintf(outfp, "Computing the Holder (L_infinity) smoothness\n");
  fprintf(outfp,"Support of Mask=[%d,%d][%d,%d], SumRule=%d\n",-Sx,Sx,-Sy,Sy,Delta);
  fprintf(outfp,"Nonzero elments in the Mask is: %d\n", Dim);
#else
  Mask=m_alloc(Sx, Sy);
  for(ii=0;ii<9999999;++ii){
    fscanf(infp, "%d\n", &i); if (i<0) break; /* use -1 to end input */
    fscanf(infp, "%d\n", &j); fscanf(infp, "%lf\n", &Mask[i][j]);
  }
  --Sx; 	--Sy;
  ConMask=m_alloc(2*Sx+1,2*Sy+1);
  /***   Do convolution here  ***/
  for ( i=0; i<(2*Sx+1); ++i )
    for ( j=0; j<(2*Sy+1); ++j ){
      ConMask[i][j]=0.0;
      for ( ii=MAX(i-Sx); ii<=MIN(i, Sx); ++ii)
	for ( jj=MAX(j-Sy); jj<=MIN(j, Sy); ++jj)
	  ConMask[i][j]+=Mask[ii][jj]*Mask[Sx+ii-i][Sy+jj-j];
      if( fabs(ConMask[i][j])>0.000000000001 ){
	Listx[Dim]=i-Sx;	Listy[Dim]=j-Sy; ++Dim;
      }
    }
  fprintf(outfp, "Computing the Sobolev (L_2) smoothness\n");
  fprintf(outfp,"Support of Mask=[%d,%d][%d,%d], SumRule=%d\n",0,Sx,0,Sy,Delta);
  fprintf(outfp, "Nonzero elments in the convoluted Mask is: %d\n", Dim);
#endif
  
  /* add the basic position here, otherwise the result will be misleading ***/
  /* generate the invariant space without assuming any symmetry. */
  Px[0]=0;   Py[0]=0;   SZ=1;
#ifdef IsIdentityDLT /*the dilation is a multiple of the identity */
  Num=0;
  for(k=0; k<Dim; k++) {
    i=abs(Listx[k])+abs(Listy[k]);
    if (i>Num ) Num=i;
  }
  ii=Sx/(dilation[0][0]-1);  jj=Sy/(dilation[0][0]-1);
  Num=Num/(dilation[0][0]-1);
  for(i=-ii; i<=ii; i++)
    for(j=-jj; j<=jj; j++) {
      if( ((i!=0) || (j!=0)) && ((abs(i)+abs(j))<=Num ) ){
	Px[SZ]=i;  Py[SZ]=j; SZ++;
      }
    }
#else
  for( i=1;i<=Delta; ++i){
    Px[SZ]=i;  Py[SZ]=0;  ++SZ;
    Px[SZ]=0;  Py[SZ]=i;  ++SZ;
  }
  /*** generate the minimal invariant full symmtric space here ***/
  do {
    Num=SZ;
    for(k=0; k<Num; ++k)
      for(iii=0; iii<Dim; ++iii){
	i1=Listx[iii]+Px[k];	j1=Listy[iii]+Py[k];
	i=dilation[1][1]*i1-dilation[0][1]*j1; /* compute M^{-1}[i1, j1] */
	j=dilation[0][0]*j1-dilation[1][0]*i1;
	if( ((i % absdet)==0) && ((j % absdet)==0)) { /* add newpoints */ 
	  index=0; i=i/det;   j=j/det;
	  for(kk=0; kk<SZ; ++kk)
	    if( (i==Px[kk]) && (j==Py[kk]) ) {
	      index=1; break;
	    }
	  if(index==0) { Px[SZ]=i; Py[SZ]=j; ++SZ; }
	}
      }
  } while( Num!=SZ );
#endif
  
  fprintf(outfp, "The invariant space before symmetry has %d points\n", SZ);
  
  Edge=0;
#ifdef FullAxesSym
  /*generate the special_eig array, m(j)=floor(j/2)+1*/
  index=0;
  for(j=0;j<MaxNumEig; j++)
    for(i=0;i<=((int)(j/2)); i++)
      if(index<MaxNumEig)
	special_eig[index++]=j;

  for(i=1; i<SZ; i++) {
    /*check either Px[i]=0, or Py[i]=0, or Px[i]=Py[i], in first quadrant.
      those are the boudnary points*/
    if( ((Py[i]==0) || (Px[i]==Py[i])) && (Px[i]>=0) ) {
      Edge++; ii=Px[Edge]; jj=Py[Edge];
      Px[Edge]=Px[i];  Py[Edge]=Py[i];
      Px[i]=ii;        Py[i]=jj;
    }
  }
  /* add the inner points here */
  Num=Edge;
  for(i=Edge+1; i<SZ; i++) {
    if( (Py[i]>0) && (Px[i]>Py[i]) ) {
      Num++; Px[Num]=Px[i];  Py[Num]=Py[i];
    }
  }
  SZ=Num+1;
  printf("The FullAxes-symmetric invariant space has %d points\n", SZ);
#endif
  
#ifdef HexagonalSym
  /*generate the special_eig array, m(j)=floor(j/3)+1*/
  index=0;
  for(j=0;j<MaxNumEig; j++)
    for(i=0;i<=((int)(j/3)); i++)
      if(index<MaxNumEig)
	special_eig[index++]=j;
  
  for(i=1; i<SZ; i++) {
    /*check either Px[i]=0, or Py[i]=0, or Px[i]=Py[i], in fourth quadrant.
      those are the boudnary points*/
    if( ((Py[i]==0) || ((Px[i]+Py[i])==0)) && (Px[i]>=0) ) {
      Edge++; ii=Px[Edge]; jj=Py[Edge];
      Px[Edge]=Px[i];  Py[Edge]=Py[i];
      Px[i]=ii;        Py[i]=jj;
    }
  }
  /* add the inner points here */
  Num=Edge;
  for(i=Edge+1; i<SZ; i++) {
    if( (Py[i]<0) && ((Px[i]+Py[i])>0) ) {
      Num++; Px[Num]=Px[i];  Py[Num]=Py[i];
    }
  }
  SZ=Num+1;
  printf("The HexagonalSym symmetric invariant space has %d points\n", SZ);
#endif
  
#ifdef AxesSym
  /*generate the special_eig array, m(j)=j+1*/
  index=0;
  for(j=0;j<MaxNumEig; j++)
    for(i=0;i<=j; i++)
      if(index<MaxNumEig) 
	special_eig[index++]=j;

  for(i=1; i<SZ; i++) {
    if( ((Px[i]==0) && (Py[i]>=0)) || ((Py[i]==0) && Px[i]>=0) ) {
      Edge++; ii=Px[Edge]; jj=Py[Edge];
      Px[Edge]=Px[i];  Py[Edge]=Py[i];
      Px[i]=ii;        Py[i]=jj;
    }
  }
  /* add the inner points here */
  Num=Edge;
  for(i=Edge+1; i<SZ; i++) {
    if( (Px[i]>=0) && (Py[i]>=0) ) {
      Num++; Px[Num]=Px[i];  Py[Num]=Py[i];
    }
  }
  SZ=Num+1;
  printf("The AxesSym symmetric invariant space has %d points\n", SZ);
#endif
  
#ifdef XYSym
  /*generate the special_eig array, m(j)=j+1*/
  index=0;
  for(j=0;j<MaxNumEig; j++)
    for(i=0;i<=j; i++)
      if(index<MaxNumEig)
	special_eig[index++]=j;

  for(i=1; i<SZ; i++) {
    /*check either Py[i]=Px[i], or Py[i]=-Px[x] on Py[i]>=0,
      those are the boudnary points*/
    if( (Py[i]>=0) && ( abs(Px[i])==Py[i]) ) {
      Edge++; ii=Px[Edge]; jj=Py[Edge];
      Px[Edge]=Px[i];  Py[Edge]=Py[i];
      Px[i]=ii;        Py[i]=jj;
    }
  }
  /* add the inner points here */
  Num=Edge;
  for(i=Edge+1; i<SZ; i++) {
    if((Py[i]>=0) && ( abs(Px[i])<Py[i] ) ) {
      Num++;  Px[Num]=Px[i];  Py[Num]=Py[i];
    }
  }
  SZ=Num+1;
  printf("The XYSym symmetric invariant space has %d points\n", SZ);
#endif
  
#ifdef OriginSym
  /*generate the special_eig array, m(j)=2*j+1*/
  index=0;
  for(j=0;j<MaxNumEig; j++)
    for(i=0;i<=(2*j); i++)
      if(index<MaxNumEig)
	special_eig[index++]=j;

  /* add the inner points here */
  Num=Edge;
  for(i=Edge+1; i<SZ; i++) {
    if( (Py[i]>0) || ((Py[i]==0) && (Px[i]>0)) ) {
      Num++;  Px[Num]=Px[i];  Py[Num]=Py[i];
    }
  }
  SZ=Num+1;
  printf("The OriginSym symmetric invariant space has %d points\n", SZ);
#endif
  
  
  /**** Allocate memory here ***/
  DT=m_alloc(SZ,SZ);
  MatDT=mxCreateDoubleMatrix( (SZ), (SZ), mxREAL );
  mxSetName(MatDT, "MatDT");
  MatDTPr=mxGetPr(MatDT);
  /*** End of allocate memory ***/
  
  /*** Generate the basic matrix here ***/
#ifdef FullAxesSym
  for ( i=0; i<SZ; ++i)
    for ( j=0; j<SZ; ++j){
      /* (iii,jjj)=dilation*(Px[j],Py[j]) */
      iii=dilation[0][0]*Px[j]+dilation[0][1]*Py[j]; 
      jjj=dilation[1][0]*Px[j]+dilation[1][1]*Py[j]; 
      
      ii=iii-Px[i];	 jj=jjj-Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j]=ConMask[ii+Sx][jj+Sy];
      else
	DT[i][j]=0.0;
      ii=iii+Px[i];	 jj=jjj+Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      ii=iii+Px[i];	 jj=jjj-Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      ii=iii-Px[i];	 jj=jjj+Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      ii=iii-Py[i];	 jj=jjj-Px[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];	
      ii=iii+Py[i];	 jj=jjj+Px[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      ii=iii+Py[i];	 jj=jjj-Px[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      ii=iii-Py[i];	 jj=jjj+Px[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
    }
  
  /* modify it to do the calcluation */
  for ( i=0; i<SZ; i++)
    DT[0][i]/=8.0;
  for ( j=1; j<=Edge; ++j)
    for ( i=0; i<SZ; i++)
      DT[j][i]/=2.0;
#endif
  
#ifdef XYSym
  for ( i=0; i<SZ; ++i)
    for ( j=0; j<SZ; ++j){
      iii=dilation[0][0]*Px[j]+dilation[0][1]*Py[j]; 
      jjj=dilation[1][0]*Px[j]+dilation[1][1]*Py[j];   
      
      ii=iii-Px[i];	 jj=jjj-Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j]=ConMask[ii+Sx][jj+Sy];
      else
	DT[i][j]=0;
      ii=iii+Px[i];	 jj=jjj+Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      ii=iii+Py[i];	 jj=jjj+Px[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      ii=iii-Py[i];	 jj=jjj-Px[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
    }
  
  /* modify it to do the calcluation */
  for ( i=0; i<SZ; i++)
    DT[0][i]/=4.0;
  for ( j=1; j<=Edge; ++j)
    for ( i=0; i<SZ; i++)
      DT[j][i]/=2.0;
#endif
  
#ifdef AxesSym
  for ( i=0; i<SZ; ++i)
    for ( j=0; j<SZ; ++j){
      iii=dilation[0][0]*Px[j]+dilation[0][1]*Py[j]; 
      jjj=dilation[1][0]*Px[j]+dilation[1][1]*Py[j];
      
      ii=iii-Px[i];      jj=jjj-Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j]=ConMask[ii+Sx][jj+Sy];
      else
	DT[i][j]=0;
      ii=iii+Px[i];      jj=jjj+Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      ii=iii+Px[i];      jj=jjj-Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      ii=iii-Px[i];      jj=jjj+Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
    }
  /* modify it to do the calcluation */
  for ( i=0; i<SZ; i++)
    DT[0][i]/=4.0;
  for ( j=1; j<=Edge; ++j)
    for ( i=0; i<SZ; i++)
      DT[j][i]/=2.0;
#endif
  
#ifdef OriginSym
  for ( i=0; i<SZ; ++i)
    for ( j=0; j<SZ; ++j){
      iii=dilation[0][0]*Px[j]+dilation[0][1]*Py[j]; 
      jjj=dilation[1][0]*Px[j]+dilation[1][1]*Py[j];
      
      ii=iii-Px[i];      jj=jjj-Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j]=ConMask[ii+Sx][jj+Sy];
      else
	DT[i][j]=0.0;
      ii=iii+Px[i];      jj=jjj+Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
    }
  
  /* modify it to do the calcluation */
  for ( i=0; i<SZ; i++)
    DT[0][i]/=2.0;
#endif
   
#ifdef HexagonalSym
  for ( i=0; i<SZ; ++i)
    for ( j=0; j<SZ; ++j){
      iii=dilation[0][0]*Px[j]+dilation[0][1]*Py[j]; 
      jjj=dilation[1][0]*Px[j]+dilation[1][1]*Py[j];
      
      ii=iii-Px[i];	 jj=jjj-Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j]=ConMask[ii+Sx][jj+Sy];
      else
	DT[i][j]=0.0;
      ii=iii+Px[i];	 jj=jjj+Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      
      ii=iii-Py[i];	 jj=jjj+Px[i]-Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      ii=iii+Py[i];	 jj=jjj-Px[i]+Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      
      ii=iii+Py[i]-Px[i];	 jj=jjj-Px[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];	
      ii=iii-Py[i]+Px[i];	 jj=jjj+Px[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      
      ii=iii+Py[i];	 jj=jjj+Px[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      ii=iii-Py[i];	 jj=jjj-Px[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      

      ii=iii+Px[i]-Py[i];	 jj=jjj-Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      ii=iii-Px[i]+Py[i];	 jj=jjj+Py[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      
      ii=iii-Px[i];	 jj=jjj+Py[i]-Px[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy];
      ii=iii+Px[i];	 jj=jjj-Py[i]+Px[i];
      if ( abs(ii)<=Sx && abs(jj)<=Sy )
	DT[i][j] += ConMask[ii+Sx][jj+Sy]; 
    }
  
  /* modify it to do the calcluation */
  for ( i=0; i<SZ; i++)
    DT[0][i]/=12.0;
  for ( j=1; j<=Edge; ++j)
    for ( i=0; i<SZ; i++)
      DT[j][i]/=2.0;
#endif
  
  /* push int matlab engine */
  index=0;
  for ( j=0; j<SZ; ++j )
    for ( i=0; i<SZ; ++i )
      MatDTPr[index++]=DT[i][j];
  
  engPutArray(ep, MatDT);
  engEvalString(ep, "EV=abs(eig(MatDT))");
  jsp=mxGetPr(engGetArray(ep, "EV") );
  qsort( jsp, SZ, sizeof(double), compare );
  
#ifdef Holder
  for( i=0; i<MIN(MaxNumEig,SZ); ++i)
    jsp[i]=-2.0*log(jsp[i])/log(1.0*absdet);
  fprintf(outfp, "\nCompting Holder (L_infinity) Smoothness exponent...\n");
  Delta=Delta*Delta;
  sum=2.0;
#else
  for( i=0; i<MIN(MaxNumEig,SZ); ++i)
    jsp[i]=1.0-log(jsp[i])/log(1.0*absdet);
  fprintf(outfp, "\nComputing Sobolev (L_2) Smoothness exponent...\n");
  Delta=Delta*Delta;
  sum=1.0;
#endif
  Delta=2*Delta;
  fprintf(outfp, "Computed values      |  Special values to be deleted!\n");
  for( i=0; i<MIN(MaxNumEig, Delta); i++){
    fprintf(outfp, "%20.16f   %20.16f\n", jsp[i], sum*special_eig[i]);
    if( fabs(jsp[i]-sum*special_eig[i])>TOL ) {
      sm=jsp[i];
      fprintf(outfp, "--------------------------------------------\n");
      for( j=i+1; j<MIN(i+6, SZ); j++)
	fprintf(outfp, "%20.16f   %20.16f\n", jsp[j], sum*special_eig[j]);
      break;
    }
  }
#ifdef Holder
  fprintf(outfp, "\nThe Holder (L_infinity) smoothness exponent =%20.16f!\n\n", sm);
#else
  fprintf(outfp, "\nThe Sobolev (L_2) smoothness exponent =%20.16f!\n\n", sm);
#endif
  engClose(ep);
  return 0;
} /* end of main */