/****************************************************************************
This program is used to compute L_2 soothness with dimension s=1,2,3
and to compute the L_infinty smoothness with positive symbol and s=1
with the dilation matrix 2I_s for the dimmensions 1, 2, 3.

!!!
Note that: For this program, I haven't had a chance to add more
documentation to make it more readable or to clean the program a little
bit to make it better. I will do so in the future. So watch out to get a
better progam than this program from my web.
If you have difficult in using this program or find some problems, please
email me at : bhan@ualberta.ca
The structure of this program is essential the same as the one d2sm.c.
So you may read the readme file for d2sm.c to understand what we have
done here.
!!! 

Set Rs_Dim=1,2,3. If Rs_Dim=2, you need to choose symmetric pattern!
                  If Rs_Dim=1, you need to choose L_2 or L_infinity

Usage: l2sm  maskfile
   or: l2sm maskfile output file
In the first case, the output is assumed to be the screen.
 
When s=1, it will tell you the L_2 smoothness of any mask.
When s=1, if the symbol is positive, you can get its
             L_infinity smoothness
Your input file must be:
-------------
Sx  Sumrule
i  Mask[i]
.
.
.
-1
-------------
The mask is supported on [0, Sx-1]
============================================================================
When s=2:
Find L2 smoothness for any bivariate mask with dialtion 2I
The method using here can deal with any mask and any sum rules.
There are four cases about the symmetry of the original mask.

X_Y_XYAxes_SYM   the mask is symmetric about y=x, y=-x, x=0, and y=0.
XYAxes_SYM       the mask is symmetric about x=0, and y=0.
X_Y_SYM          the mask is symmetric about y=x, and y=-x.
NO_SYM           the mask has no symmetry at all.
Hex_SYM          The mask has hexagonal symmetric, i.e., rotation by 120C
Hex_FullSYM      The mask has hexagonal full symmetry,
                 that is, its symmetry group has 12 elements.

You need to define one of the above to find the L_2 smoothness of the
scaling function.

You  need to put your mask into a file as
------------------------
Sx   Sy   SumRule
i    j    Mask[i,j]
.
.
.
-1
--------------------
Use -1 to end the input. Only nonzero Mask[i,j] is need to be put into the
file. The Mask is supported on [0,Sx-1][0,Sy-1] with sum rules of order SumRule.
Here we shift the mask so that the program is easier.
The program automatically checks the symmetry of the resulting convolution
of the mask as you choose among the four knids of symmetries.
==========================================================================

When s=3:
Find L2 smoothness for any 3-D mask with dialtion 2I
The mathod using here can deal with any mask with symmetry and any sum rules.
Just look at the results and find the first non integer value
which will be the result in most case. And we present actually
8 more value here for conveneience

Method : We use the full symmetry of the mask. symmtric about x, y,z-axis,
x=y, x=-y, ....  ie. all possible symmetry.
and we reduce the size of the transition matrix to be approx 1/48 of
the original one ( a matrix (r by r) is said to has size r).

maskfile must take the following format:
Sx  Sy  Sz
i   j   k   Mask[i][j][k]
.
.
.
-1
----
Mas is supported on [0,Sx-1][0,Sy-1][0,Sz-1]

The result is put into resultfile.

Precondition:
	The 3-D mask must have full symmetry with support [-Sx,Sx]^3.
	We don't check the symmetry of the mask in this program!!!!
****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include "/usr/local/matlab5/extern/include/engine.h"
#define MAX(x,y)  ((x)>(y) ? (x):(y))
#define MIN(x,y)  ((x)<(y) ? (x):(y))
#define MAXELM  50000
#define TOL  0.0001

/** you need to define the following before compile this program ***/
#undef Rs_Dim1
#define Rs_Dim2
#undef Rs_Dim3

/** when Rs_Dim1 is defined, you need to define or undefine Dim1_L0 **/
#define Dim1_L2   /* undefine Dim1_L2 to compute L_infinity smoothness */

/** when Rs_Dim2 is define, you must define the following macro **/
#undef NO_SYM         /* the mask has no symmetry at all */
#undef XYAxes_SYM     /* the mask is symmetric about x=0, and y=0 */
#undef X_Y_SYM        /* the mask is symmetric about y=x, and y=-x  */
#undef X_Y_XYAxes_SYM /*the mask is symmetric about y=x, y=-x, x=0, and y=0 */
#define Hex_FullSYM    /*the mask is symmetric about rotated by 30C */
#undef Hex_SYM        /*the mask is symmetric about rotated by 120C */

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

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

double **m2_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 m2_alloc */

double ***m3_alloc( row, col, height)
int row, col, height;
{
int i, j;
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]=m2_alloc(col, height);
  if ( a[i]==NULL ) {
    fprintf( stderr, " Could not allocate memery!\n" );
    exit(1);
  } /* end of if */
} /* end of for(i=..) */
return a;
} /* end of m3_alloc */

void delete2(a, row)
double **a;
int row;
{
int i;

for(i=0; i<row; ++i)
  free(a[i]);
free(a);
} /* delete2 */


void main(argc, argv)
int argc;
char *argv[];
{

int i,i1,j1,k1, ii, iii, j, jj, jjj,k,kk,kkk, Edge, Edge1, Edge2, Edge3, Surface;
int Sx, Sy, Sz, SZ, Delta, Dim, index, Num, Bound;
                       /* Bound is introduced here to further reduce the size */
int Px[MAXELM], Py[MAXELM], Pz[MAXELM];
int Listx[MAXELM], Listy[MAXELM];

double *MatDTPr, **DT;
double sum, Sm, *jsp;
Engine *ep;
mxArray *MatDT;
FILE *outfp, *infp;
#ifdef Rs_Dim1
    double *Mask, *ConMask;
#endif
#ifdef Rs_Dim2
   double **Mask, **ConMask;
#endif
#ifdef Rs_Dim3
   double ***Mask, ***ConMask;
#endif

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: l2sm maskfile (outputfile)! exit!\n");
  exit(0);
}

if ( (infp==NULL) || (outfp==NULL) ) {
  fprintf(stderr, "Can not open I/O files! exit!\n");
  exit(0);
}

if ( !( ep=engOpen("\0")) ){
  printf( "\n Can't start MATLAB engine\n"); exit(-1);
}

/*************************************************************
  Compute the L_2 smoothness for 1 dimensional mask.
**************************************************************/
#ifdef Rs_Dim1

fscanf(infp, "%d\n", &Sx);  fscanf(infp, "%d\n", &Delta);
#ifdef Dim1_L2
    Mask=(double *)calloc( Sx, sizeof( double ));
    ConMask=(double *)calloc( 2*Sx+1, sizeof( double ));
    if ( (Mask==NULL) || (ConMask==NULL) ) {
       fprintf( stderr, " Could not allocate memery!\n" );
       exit(1);
    }
    iii=0;
    for(ii=0;ii<9999999;++ii){
        fscanf(infp,"%d\n", &i); if ( i < 0 ) break; /* use -1 to end input */
        fscanf(infp, "%lf\n", &Mask[i]);
    ++iii;
}
--Sx;
fprintf(outfp,"Mask is supported in [0..%d] with sum rule=%d\n", Sx,Delta);

for ( i=0; i<(Sx+1); ++i ) {
  ConMask[i]=0.0;
  for ( ii=MAX(i-Sx,0); ii<=MIN(i, Sx); ++ii)
    ConMask[i]+=Mask[ii]*Mask[Sx+ii-i];
  ConMask[2*Sx-i]=ConMask[i];
}
#else
    ConMask=(double *)calloc( Sx, sizeof( double ));
    if ( (ConMask==NULL) ) {
       fprintf( stderr, " Could not allocate memery!\n" );
       exit(1);
    }
    iii=0;
    for(ii=0;ii<9999999;++ii){
        fscanf(infp,"%d\n", &i); if ( i < 0 ) break; /* use -1 to end input */
        fscanf(infp, "%lf\n", &ConMask[i]);
    ++iii;
}
--Sx;
fprintf(outfp,"Mask is supported in [0..%d] with sum rule=%d\n", Sx,Delta);
Sx/=2;   Delta/=2;
#endif

SZ=Sx+1;
Delta=Delta*Delta;
DT=m2_alloc(SZ,SZ);
MatDT=mxCreateDoubleMatrix( (SZ), (SZ), mxREAL );
mxSetName(MatDT, "MatDT");
MatDTPr=mxGetPr(MatDT);
/*** End of allocate memory ***/
 
for ( i=0; i<=Sx; ++i)
for ( j=0; j<=Sx; ++j) {
  ii=2*j-i;
  if ( abs(ii)<=Sx )
    DT[i][j]=ConMask[ii+Sx];
  else
    DT[i][j]=0.0;
  ii=2*j+i;
  if ( abs(ii)<=Sx )
    DT[i][j] += ConMask[ii+Sx];
}

/* modify it to do the calcluation */
for ( i=0; i<=Sx; i++)
  DT[0][i]/=2.0;

#endif
/*============================================================
  End of Computing L_2 smoothness for 1 dimensional mask.
============================================================*/


/*************************************************************
  Compute the L_2 smoothness for one dimensional mask.
**************************************************************/
#ifdef Rs_Dim2

Px[0]=0;Py[0]=0;
fscanf(infp, "%d\n", &Sx);
fscanf(infp, "%d\n", &Sy);
fscanf(infp, "%d\n", &Delta);
Mask=m2_alloc(Sx, Sy);
iii=0;
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]);
  ++iii;
}
--Sx; 		--Sy;
fprintf(outfp,"Mask is supported in [0..%d][0..%d] with sum rule=%d\n", Sx,Sy,Delta);
fprintf(outfp,"Nonzero elements in Mask is %d\n\n", iii);

ConMask=m2_alloc(2*Sx+1,2*Sy+1);

/***   Do convolution here  ***/
Bound=0;
for ( i=0; i<(Sx+1); ++i )
for ( j=0; j<(2*Sy+1); ++j ){
  ConMask[i][j]=0.0;
  for ( ii=MAX(i-Sx,0); ii<=MIN(i, Sx); ++ii)
  for ( jj=MAX(j-Sy,0); jj<=MIN(j, Sy); ++jj)
    ConMask[i][j]+=Mask[ii][jj]*Mask[Sx+ii-i][Sy+jj-j];
  if ( fabs( ConMask[i][j] != 0 ) ) {
    kk=MAX( abs(i-Sx+j-Sy), abs(i-Sx-(j-Sy)) );
    if( Bound < kk)
      Bound=kk;
  }
  ConMask[2*Sx-i][2*Sy-j]=ConMask[i][j];
}

SZ=1;
sum=0.0;

#ifdef X_Y_SYM
/* add the edge here */
ii=Bound/2;
for( i=1; i<= ii; i++ ) {
  Px[SZ]=i;	Py[SZ]=i;   SZ++;
  Px[SZ]=-i;	Py[SZ]=i;   SZ++;
}
Edge=SZ-1;

/* add inner points in the invariant region */
for(j=1; j<=Sx;  j++)
for(i=-j+1; i<j; i++) {
  if( ( abs(i)+abs(j)) <= Bound ) {
    Px[SZ]=i;   Py[SZ]=j;   ++SZ;
  }
}

/* check the desired symmetric for the convolution of the mask her */
for(i=0; i<SZ; i++)
    sum += (fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx-Px[i]][Sy-Py[i]])+
	    fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx+Py[i]][Sy+Px[i]])+
	    fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx-Py[i]][Sy-Px[i]]));	
if( sum<TOL )
    fprintf(outfp,"The symmetric (y=x,y=-x) invariant space has %d points\n\n", SZ);
else {
    fprintf("The mask is NOT symmetric about (y=x,y=-x)!! exit!\n");
    exit(0);
}
#endif

#ifdef XYAxes_SYM
/* add the edge here */
for( i=1; i<= Sx; i++ ) {
  Px[SZ]=i;	Py[SZ]=0;   SZ++;
}
for( i=1; i<=Sy; i++ ) {
  Px[SZ]=0;	Py[SZ]=i;   SZ++;
}
Edge=SZ-1;

/* add inner points in the invariant region */
for(i=1; i<=Sx; i++)
for(j=1; j<=Sy; j++) {
  if( ( abs(i+j)) <= Bound ) {
    Px[SZ]=i;   Py[SZ]=j;   ++SZ;
  }
}
for(i=0; i<SZ; i++)
    sum += (fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx-Px[i]][Sy+Py[i]])+
	    fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx+Px[i]][Sy-Py[i]])+
	    fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx-Px[i]][Sy-Py[i]]));	
if( sum<TOL )
    fprintf(outfp,"The symmetric (x=0,y=0) invariant space has %d points\n\n", SZ);
else {
    fprintf(outfp,"The mask is NOT symmetric about (x=0,y=0)!! exit!\n");
    exit(0);
}
#endif

#ifdef Hex_FullSYM
/* add the edge here, here we consider the area in 4th quadrant*/
/* here we play a small trick here */
for( i=1; i<= (Sx-1); i++ ) {
  Px[SZ]=i;	Py[SZ]=0;   SZ++;
}
for( i=1; i<=((int)((Sx-1)/2)); i++ ) {
  Px[SZ]=i;	Py[SZ]=-i;   SZ++;
}
Edge=SZ-1;

/* add inner points in the invariant region */
for(i=1; i<Sx; i++)
for(j=1; j<i; j++) {
  if( (i+j)<=(Sx-1) ) {
    Px[SZ]=i;   Py[SZ]=-j;   ++SZ;
  }
}

for(i=0; i<SZ; i++)
sum +=(fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx-Px[i]][Sy-Py[i]])+
      fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx-Py[i]][Sy+Px[i]-Py[i]])+
      fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx+Py[i]][Sy-Px[i]+Py[i]])+
      fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx+Py[i]-Px[i]][Sy-Px[i]])+
      fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx-Py[i]+Px[i]][Sy+Px[i]])+
      fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx+Py[i]][Sy+Px[i]])+
      fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx-Py[i]][Sy-Px[i]])+
      fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx+Px[i]-Py[i]][Sy-Py[i]])+
      fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx-Px[i]+Py[i]][Sy+Py[i]])+
      fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx-Px[i]][Sy+Py[i]-Px[i]])+
      fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx+Px[i]][Sy-Py[i]+Px[i]]));
  if( sum<TOL )
    fprintf(outfp,"The symmetric (rotate about 30C) invariant space has %d points\n\n", SZ);
else {
    fprintf(outfp,"The mask is NOT symmetric with rotation 30C!! exit!\n");
    exit(0);
}
#endif

#ifdef NO_SYM
for(i=-Sx; i<=0; i++)
for(j=MAX(i-Bound, -Sy); j<=MIN(Bound-i, Sy); j++) {
  if( (i==0) && (j==0) )
    break;
  Px[SZ]=i;   Py[SZ]=j;   ++SZ;
}
fprintf(outfp,"The origin-symmetric invariant space has %d points\n\n", SZ);
#endif

#ifdef X_Y_XYAxes_SYM
/* add the edge here */
for( i=1; i<= Sx; i++ ) {
  Px[SZ]=i;	Py[SZ]=0;   SZ++;
}
ii=Bound/2;
for( i=1; i<= ii; i++ ) {
  Px[SZ]=i;	Py[SZ]=i;   SZ++;
}
Edge=SZ-1;

/* add inner points in the invariant region */
for(i=1; i<=Sx; i++)
for(j=1; j<=MIN(i-1, Bound-i); j++) {
  Px[SZ]=i;   Py[SZ]=j;   ++SZ;
}

for(i=0; i<SZ; i++)
    sum += (fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx-Px[i]][Sy+Py[i]])+
	    fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx+Px[i]][Sy-Py[i]])+
	    fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx-Px[i]][Sy-Py[i]])+
 	    fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx+Py[i]][Sy+Px[i]])+
	    fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx-Py[i]][Sy+Px[i]])+
 	    fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx+Py[i]][Sy-Px[i]])+
	    fabs(ConMask[Sx+Px[i]][Sy+Py[i]]-ConMask[Sx-Py[i]][Sy-Px[i]])); 
if( sum<TOL ) {
fprintf(outfp,"The symmetric (x=0,y=0,y=x,y=-x) invariant space has %d points\n\n",SZ);
} else {
    fprintf(outfp,"The mask is NOT symmetric about (x=0,y=0,y=x,y=-x)! exit!\n");
    exit(0);
}
#endif

/**** Allocate memory here ***/
Delta=Delta*Delta;
DT=m2_alloc(SZ,SZ);
MatDT=mxCreateDoubleMatrix( (SZ), (SZ), mxREAL );
mxSetName(MatDT, "MatDT");
MatDTPr=mxGetPr(MatDT);
/*** End of allocate memory ***/

/*** Generate the basic matrix here ***/
#ifdef X_Y_SYM
for ( i=0; i<SZ; ++i)
for ( j=0; j<SZ; ++j){
  ii=2*Px[j]-Px[i];	 jj=2*Py[j]-Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j]=ConMask[ii+Sx][jj+Sy];
  else
    DT[i][j]=0;
  ii=2*Px[j]+Px[i];	 jj=2*Py[j]+Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];
  ii=2*Px[j]+Py[i];	 jj=2*Py[j]+Px[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];
  ii=2*Px[j]-Py[i];	 jj=2*Py[j]-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 XYAxes_SYM
for ( i=0; i<SZ; ++i)
for ( j=0; j<SZ; ++j){
  ii=2*Px[j]-Px[i];	 jj=2*Py[j]-Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j]=ConMask[ii+Sx][jj+Sy];
  else
    DT[i][j]=0;
  ii=2*Px[j]+Px[i];	 jj=2*Py[j]+Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];
  ii=2*Px[j]+Px[i];	 jj=2*Py[j]-Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];
  ii=2*Px[j]-Px[i];	 jj=2*Py[j]+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 NO_SYM
for ( i=0; i<SZ; ++i)
for ( j=0; j<SZ; ++j){
  ii=2*Px[j]-Px[i];	 jj=2*Py[j]-Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j]=ConMask[ii+Sx][jj+Sy];
  else
    DT[i][j]=0.0;
  ii=2*Px[j]+Px[i];	 jj=2*Py[j]+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 X_Y_XYAxes_SYM
for ( i=0; i<SZ; ++i)
for ( j=0; j<SZ; ++j){
  ii=2*Px[j]-Px[i];	 jj=2*Py[j]-Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j]=ConMask[ii+Sx][jj+Sy];
  else
    DT[i][j]=0.0;
  ii=2*Px[j]+Px[i];	 jj=2*Py[j]+Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];
  ii=2*Px[j]+Px[i];	 jj=2*Py[j]-Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];
  ii=2*Px[j]-Px[i];	 jj=2*Py[j]+Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];
  ii=2*Px[j]-Py[i];	 jj=2*Py[j]-Px[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];	
  ii=2*Px[j]+Py[i];	 jj=2*Py[j]+Px[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];
  ii=2*Px[j]+Py[i];	 jj=2*Py[j]-Px[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];
  ii=2*Px[j]-Py[i];	 jj=2*Py[j]+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 Hex_FullSYM
for ( i=0; i<SZ; ++i)
for ( j=0; j<SZ; ++j){
  ii=2*Px[j]-Px[i];	 jj=2*Py[j]-Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j]=ConMask[ii+Sx][jj+Sy];
  else
    DT[i][j]=0.0;
  ii=2*Px[j]+Px[i];	 jj=2*Py[j]+Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];

  ii=2*Px[j]-Py[i];	 jj=2*Py[j]+Px[i]-Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];
  ii=2*Px[j]+Py[i];	 jj=2*Py[j]-Px[i]+Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];

  ii=2*Px[j]+Py[i]-Px[i];	 jj=2*Py[j]-Px[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];	
  ii=2*Px[j]-Py[i]+Px[i];	 jj=2*Py[j]+Px[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];

  ii=2*Px[j]+Py[i];	 jj=2*Py[j]+Px[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];
  ii=2*Px[j]-Py[i];	 jj=2*Py[j]-Px[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];


  ii=2*Px[j]+Px[i]-Py[i];	 jj=2*Py[j]-Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];
  ii=2*Px[j]-Px[i]+Py[i];	 jj=2*Py[j]+Py[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];

  ii=2*Px[j]-Px[i];	 jj=2*Py[j]+Py[i]-Px[i];
  if ( abs(ii)<=Sx && abs(jj)<=Sy )
    DT[i][j] += ConMask[ii+Sx][jj+Sy];
  ii=2*Px[j]+Px[i];	 jj=2*Py[j]-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


#endif
/*=================================================================
  End of computing L_2 smoothness of 2 dimensional mask
==================================================================*/


#ifdef Rs_Dim3

fscanf(infp,"%d\n", &Sx); fscanf(infp,"%d\n", &Sy); fscanf(infp,"%d\n", &Sz);
fscanf(infp, "%d\n", &Delta);
Mask=m3_alloc(Sx, Sy, Sz);
--Sx; 		--Sy;		--Sz;
i1=Sx/2;        j1=Sy/2;        k1=Sz/2;
Num=0; /* Before shift, Mask is supported on {(i,j,k) : |i|+|j|+|k|<=Num} */
iii=0;
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, "%d\n", &k);
  fscanf(infp, "%lf\n", &Mask[i][j][k]);
  if( (fabs(Mask[i][j][k]) != 0.0) && ((abs(i-i1)+abs(j-j1)+abs(k-k1))>Num) )
    Num=abs(i-i1)+abs(j-j1)+abs(k-k1);
  ++iii;
}

fprintf(outfp, "Mask support [0..%d][0..%d][0..%d], sumrule=%d\n", Sx,Sy,Sz,Delta);
fprintf(outfp, "Nonzero elements in Mask is %d\n", iii);
ConMask=m3_alloc(2*Sx+1, 2*Sy+1, 2*Sz+1);

Px[0]=0;Py[0]=0; Pz[0]=0;
Num =2*Num; /* ConMask is supported on |i|+|j|+|j|<= 2Num */

SZ=1;
for( i=1; i<= Sx; i++ ) {
  Px[SZ]=i;	Py[SZ]=0;	Pz[SZ]=0;
  SZ++;
}
Edge1=SZ;

for( i=1; i<= MIN(Sx, (Num/3)); i++ ) {
  Px[SZ]=i;	Py[SZ]=i;	Pz[SZ]=i;
  SZ++;
}
Edge2=SZ;

for( i=1; i<= MIN(Sx, (Num/2)); i++ ) {
  Px[SZ]=i;	Py[SZ]=i;	Pz[SZ]=0;
  SZ++;
}
Edge3=SZ;

/* surfaces */

for( i=1; i<= Sx; i++ )
for( j=1; j< i; j++){
  if( (i+2*j) <= Num ) {
    Px[SZ]=i;	Py[SZ]=j;	Pz[SZ]=j;
    SZ++;
  }
}
 
for( i=1; i<= Sx; i++ )
for( k=1; k< i; k++){
  if( (2*i+k) <= Num ) {
    Px[SZ]=i;	Py[SZ]=i;	Pz[SZ]=k;
    SZ++;
  }
}
for( i=1; i<= Sx; i++ )
for( j=1; j< i; j++){
  if( (i+j) <= Num ) {
    Px[SZ]=i;	Py[SZ]=j;	Pz[SZ]=0;
    SZ++;
  }
}
Surface=SZ;

/* solid volume */
 
for( i=1; i<= Sx; i++ )
for( j=1; j< i; j++)
for( k=1; k< j; k++){
  if( (i+j+k)<= Num ) {
    Px[SZ]=i;	Py[SZ]=j;	Pz[SZ]=k;
    SZ++;
  }
}

fprintf(outfp, "The symmetric invariant space has %d points\n\n", SZ);

/* A symmetric version of convolution, first calcuate the value
in (Px[], Py[], Pz[]), then extend them                 */

for ( index=0; index<SZ; ++index ){
  i=Sx-Px[index];	j=Sy-Py[index];	k=Sz-Pz[index];
  sum=0.0;
  for ( ii=MAX(i-Sx,0); ii<=MIN(i, Sx); ++ii)
  for ( jj=MAX(j-Sy,0); jj<=MIN(j, Sy); ++jj)
  for ( kk=MAX(k-Sz,0); kk<=MIN(k, Sz); ++kk)
    sum+=Mask[ii][jj][kk]*Mask[Sx+ii-i][Sy+jj-j][Sz+kk-k];
  
  i=Px[index];	j=Py[index];	k=Pz[index];
  for(ii=-1; ii<2; ii=ii+2)
  for(jj=-1; jj<2; jj=jj+2)
  for(kk=-1; kk<2; kk=kk+2)
    ConMask[Sx+ii*i][Sy+jj*j][Sz+kk*k]=sum;
  
  i=Py[index];	j=Pz[index];	k=Px[index];
  for(ii=-1; ii<2; ii=ii+2)
  for(jj=-1; jj<2; jj=jj+2)
  for(kk=-1; kk<2; kk=kk+2)
    ConMask[Sx+ii*i][Sy+jj*j][Sz+kk*k]=sum;

  i=Pz[index];	j=Px[index];	k=Py[index];
  for(ii=-1; ii<2; ii=ii+2)
  for(jj=-1; jj<2; jj=jj+2)
  for(kk=-1; kk<2; kk=kk+2)
    ConMask[Sx+ii*i][Sy+jj*j][Sz+kk*k]=sum;
	
  i=Py[index];	j=Px[index];	k=Pz[index];
  for(ii=-1; ii<2; ii=ii+2)
  for(jj=-1; jj<2; jj=jj+2)
  for(kk=-1; kk<2; kk=kk+2)
    ConMask[Sx+ii*i][Sy+jj*j][Sz+kk*k]=sum;
		
  i=Px[index];	j=Pz[index];	k=Py[index];
  for(ii=-1; ii<2; ii=ii+2)
  for(jj=-1; jj<2; jj=jj+2)
  for(kk=-1; kk<2; kk=kk+2) 
    ConMask[Sx+ii*i][Sy+jj*j][Sz+kk*k]=sum;

  i=Pz[index];	j=Py[index];	k=Px[index];
  for(ii=-1; ii<2; ii=ii+2)
  for(jj=-1; jj<2; jj=jj+2)
  for(kk=-1; kk<2; kk=kk+2) 
    ConMask[Sx+ii*i][Sy+jj*j][Sz+kk*k]=sum;	
}

/**** Allocate memory here ***/
Delta=Delta*Delta*Delta;
DT=m2_alloc(SZ,SZ);
MatDT=mxCreateDoubleMatrix( (SZ), (SZ), mxREAL );
mxSetName(MatDT, "MatDT");
MatDTPr=mxGetPr(MatDT);
/*** End of allocate memory ***/

/*** Generate the basic matrix here ***/
for ( i=0; i<SZ; ++i)
for ( j=0; j<SZ; ++j){
  DT[i][j]=0;
  
  iii=Px[i];	jjj=Py[i];	kkk=Pz[i];
  for(i1=-1; i1<2; i1 = i1+2)
    for(j1=-1; j1<2; j1 = j1+2)
    for(k1=-1; k1<2; k1 = k1+2) {
      ii=2*Px[j]+i1*iii; jj=2*Py[j]+j1*jjj;	kk=2*Pz[j]+k1*kkk;
      if ( abs(ii)<=Sx && abs(jj)<=Sy && abs(kk)<=Sz )
	DT[i][j]+=ConMask[ii+Sx][jj+Sy][kk+Sz];
    }
  
  iii=Py[i];	jjj=Pz[i];	kkk=Px[i];
  for(i1=-1; i1<2; i1 = i1+2)
  for(j1=-1; j1<2; j1 = j1+2)
  for(k1=-1; k1<2; k1 = k1+2) {
    ii=2*Px[j]+i1*iii; jj=2*Py[j]+j1*jjj;	kk=2*Pz[j]+k1*kkk;
    if ( abs(ii)<=Sx && abs(jj)<=Sy && abs(kk)<=Sz )
      DT[i][j]+=ConMask[ii+Sx][jj+Sy][kk+Sz];
  }
  
  iii=Pz[i];	jjj=Px[i];	kkk=Py[i];
  for(i1=-1; i1<2; i1 = i1+2)
  for(j1=-1; j1<2; j1 = j1+2)
  for(k1=-1; k1<2; k1 = k1+2) {
    ii=2*Px[j]+i1*iii; jj=2*Py[j]+j1*jjj;	kk=2*Pz[j]+k1*kkk;
    if ( abs(ii)<=Sx && abs(jj)<=Sy && abs(kk)<=Sz )
      DT[i][j]+=ConMask[ii+Sx][jj+Sy][kk+Sz];
  }

/**************************************************/
  iii=Py[i];	jjj=Px[i];	kkk=Pz[i];
  for(i1=-1; i1<2; i1 = i1+2)
  for(j1=-1; j1<2; j1 = j1+2)
  for(k1=-1; k1<2; k1 = k1+2) {
    ii=2*Px[j]+i1*iii; jj=2*Py[j]+j1*jjj;	kk=2*Pz[j]+k1*kkk;
    if ( abs(ii)<=Sx && abs(jj)<=Sy && abs(kk)<=Sz )
      DT[i][j]+=ConMask[ii+Sx][jj+Sy][kk+Sz];
  }
  
  iii=Px[i];	jjj=Pz[i];	kkk=Py[i];
  for(i1=-1; i1<2; i1 = i1+2)
  for(j1=-1; j1<2; j1 = j1+2)
  for(k1=-1; k1<2; k1 = k1+2) {
    ii=2*Px[j]+i1*iii; jj=2*Py[j]+j1*jjj;	kk=2*Pz[j]+k1*kkk;
    if ( abs(ii)<=Sx && abs(jj)<=Sy && abs(kk)<=Sz )
      DT[i][j]+=ConMask[ii+Sx][jj+Sy][kk+Sz];
  }
  
  iii=Pz[i];	jjj=Py[i];	kkk=Px[i];
  for(i1=-1; i1<2; i1 = i1+2)
  for(j1=-1; j1<2; j1 = j1+2)
  for(k1=-1; k1<2; k1 = k1+2) {
    ii=2*Px[j]+i1*iii; jj=2*Py[j]+j1*jjj;	kk=2*Pz[j]+k1*kkk;
    if ( abs(ii)<=Sx && abs(jj)<=Sy && abs(kk)<=Sz )
      DT[i][j]+=ConMask[ii+Sx][jj+Sy][kk+Sz];
  }
}

 
/* modify it to do the calcluation */
for ( i=0; i<SZ; i++)
  DT[0][i]/=48.0;
for ( j=1; j<Edge1; ++j)
for ( i=0; i<SZ; i++)
  DT[j][i]/=8.0;
for ( j=Edge1; j<Edge2; ++j)
for ( i=0; i<SZ; i++)
  DT[j][i]/=6.0;
for ( j=Edge2; j<Edge3; ++j)
for ( i=0; i<SZ; i++)
  DT[j][i]/=4.0;
for ( j=Edge3; j<Surface; ++j)
for ( i=0; i<SZ; i++)
  DT[j][i]/=2.0;
#endif
/*==========================================================
  End of computing L_2 smoothness of 3 dimensional mask.
==========================================================*/

/* push into matlab engine */
index=0;
for ( i=0; i<SZ; ++i )
for ( j=0; j<SZ; ++j )
  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 );

sum=0.0;  ii=1;
#ifdef Dim1_L2
   sum=0.5; ii=2;
#endif
#ifdef Rs_Dim2
    sum=1.0; ii=2;
#endif
#ifdef Rs_Dim3
    sum=1.5; ii=2;
#endif
Edge=0;
j=0;
Dim=1;
index=0;
Dim=0;
Delta=MIN(Delta+32, SZ);
fprintf(outfp, "eigenvalue      L_2 smoothness\n");
for( i=0; i< Delta; i++) {
  if( fabs(jsp[i])> 0.0000000001 )
    Sm=sum-log(jsp[i])/log(2.0)/ii;
fprintf(outfp, "%20.16f    %20.16f\n", jsp[i], Sm);
if( Edge>=32 ) break; 
  if((index==0)&&( fabs( Sm- ( (int)(Sm+0.01) ) ) > 0.01)){ 
    fprintf(outfp, "==================================================\n");
    index=1;
  }
  if( index==1) 
	Edge++;
  if( Sm >16.0) break;
}
engClose(ep);
} /* end of main */