/********************************************************************
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
#include
#include
/*??? 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; i0.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; kNum ) 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=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; i0) && (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=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; i0) ) {
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=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=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=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=0) && ( abs(Px[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; iTOL ) {
sm=jsp[i];
fprintf(outfp, "--------------------------------------------\n");
for( j=i+1; j