/****************************************************************************
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
#include
#include
#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; iNum) )
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 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 */