Discrete Orthonormal Moments
Originally posted by Dan on 15:59 Sun 27 April 2008, last modified 08:05 Mon 17 November 2008.
File under: image processing math moments objective-c phd programming shape description
Shape description by image moments is a popular topic in image processing. Standard geometric moments are implemented as:
which is a specific implementation of the general moment function:
in this case, the kernel of the basis function Ψpq = {xpyq} is non-orthogonal, which means that there is a high degree of redundant information in the set of moments. Additional problems with such an implementation arise from approximating continuous integrals with discrete summations, which lead to numerical inaccuracies, and the very large dynamic range between moments of different orders as a result from the powers of p and q.
Orthogonal basis functions overcome some of these problems. The information contained in each moment is unique, which means that reconstruction is possible from a finite set. Additionally the variation in dynamic range is reduced. However they often suffer from being defined in an unsuitable domain. For example, Legendre polynomials are defined in the range [-1,1] and Zernike polynomials only inside the unit circle. To overcome this problem requires a transformation which again introduces further inaccuracies.
In light of this, we require a basis kernel that is a discrete orthogonal polynomial defined directly in the image coordinate space. There are a number of polynomials that satisfies these requirements, however the simplest is Tchebichef, also called Chebychev or Tchebicheff.
Given a square image of size N, the complete set of Tchebichef moments is computed as:
where the t functions are Tchebicef polynomials calculated over the domain [0,N-1]. Given such a set of moments, the original image may be fully reconstructed. The value of an individual pixel is reconstructed by:
So what about these Tchebichef Polynomials then?
The polynomials are shown in the image below for a 64x64 image.
They are computed using a recurrence scheme over n and x, as discussed in detail in Mukundan’s paper Some Computational Aspects of Discrete Orthonormal Moments. The code extract written using the Cocoa frameworks in Objective-C shows the details:
Consider the following class definition:
#import <Cocoa/Cocoa.h>
@interface TchebichefPolynomial : NSObject <NSCoding> {
NSInteger N;
// We'll be using the symmetry property, so we want half N
NSInteger half;
NSMutableData *storage;
}
+ (NSString *)pathToArchiveForSize:(NSInteger)n;
+ (TchebichefPolynomial *)polynomialOfGivenSize:(NSUInteger)n;
- (id)initGivenSize:(NSInteger)n;
- (void)populate;
- (BOOL)errorCheck;
- (void)nRecurrenceUpto:(NSInteger)upto;
- (void)xRecurrenceUpto:(NSInteger)upto;
- (void)outputToFile:(NSString *)filename;
- (double)valueForOrder:(NSUInteger)n atPosition:(NSUInteger)x;
@property NSInteger N;
@property NSInteger half;
@property(assign,readonly) NSMutableData *storage;
@end
- (id)initGivenSize:(NSInteger)n {
self = [super init];
if(self != nil) {
N = n;
half = (NSInteger)(self.N/2.0);
storage = [[NSMutableData alloc] initWithCapacity:(half*N*sizeof(double))];
[self populate];
[self errorCheck];
}
return self;
}
The populate function simply calls the recurrence methods in the correct order.
- (void)populate {
// Use n recurrence relation to compute values for x=0
[self nRecurrenceUpto:1.0];
// Now use x recurrence relation for the rest
[self xRecurrenceUpto:1.0];
NSInteger x;
double *ptr = (double *)[storage mutableBytes];
double NSqR = sqrt(N);
for(x=0; x<half; x++) {
ptr[x] = 1.0 / NSqR;
}
}
The real guts of this class implementation are the recurrence schemes based on Mukundan's paper:
- (void)nRecurrenceUpto:(NSInteger)upto {
NSInteger index, prevNindex, prev2Nindex;
double *ptr = (double *)[storage mutableBytes];
double value, n, x, NSq, NSqR, nSq, alpha1, alpha2, alpha3;
NSq = pow(N, 2);
NSqR = sqrt(N);
for(n=0.0; n<N; n++) {
nSq = pow(n, 2);
for(x=0.0; x<upto; x++) {
index = (n*half) + x;
if(n == 0.0) {
value = 1.0 / sqrt(N);
} else if(n == 1.0) {
value = ((2*x) + 1 - N)*sqrt( (3.0 / (N*(NSq - 1))) );
} else if(n > 1.0) {
alpha1 = (2.0/n) * sqrt( ((4*nSq) - 1.0) / (NSq-nSq) );
alpha2 = ((1.0 - N)/n) * sqrt( ((4*nSq) - 1.0) / (NSq-nSq) );
alpha3 = ((n-1.0)/n) * sqrt( (((2*n) + 1.0) / ((2*n) - 3)) ) * sqrt( ((NSq - pow(n-1,2)) / (NSq-nSq)) );
prevNindex = ((n-1)*half) + x;
prev2Nindex = ((n-2)*half) + x;
value = alpha1*x*ptr[prevNindex] + alpha2*ptr[prevNindex] + alpha3*ptr[prev2Nindex];
}
ptr[index] = value;
}
}
}
- (void)xRecurrenceUpto:(NSInteger)upto {
NSInteger index, prevXindex, prev2Xindex;
double *ptr = (double *)[storage mutableBytes];
double value, n, x, NSq, nSq, gamma1, gamma2;
NSq = pow(N, 2);
for(n=1.0; n<N; n++) {
nSq = pow(n, 2);
for(x=0.0; x<half; x++) {
index = (n*half) + x;
if(x == 0.0) {
value = -sqrt( (N-n)/(N+n) ) * sqrt( ((2*n)+1)/((2*n)-1) ) * ptr[(NSInteger)((n-1)*half)];
} else if(x == 1.0) {
value = (1+((n*(1+n))/(1-N)))*ptr[(NSInteger)(n*half)];
} else if (x > 1.0) {
gamma1 = ((-n*(n+1)) - (((2*x) - 1)*(x-N-1)) - x) / (x*(N-x));
gamma2 = ((x-1)*(x-N-1)) / (x*(N-x));;
prevXindex = index - 1;
prev2Xindex = index - 2;
value = gamma1*ptr[prevXindex] + gamma2*ptr[prev2Xindex];
}
// Write the actual value
ptr[index] = value;
}
}
}
Finally, we can perform an error check to make sure that the computed polynomial values are actually orthogonal.
- (BOOL)errorCheck {
// Check values conform to Orthogonality Property
double tmp, tmp2, s = 0.0;
NSInteger n, x;
for(n=1.0; n<N; n++) {
tmp = 0;
for(x=0; x<N; x++) {
tmp2 = [self valueForOrder:n atPosition:x];
tmp += tmp2;
}
s += tmp;
}
if(abs(s) == 0.0) {
return YES;
} else {
return NO;
}
}
The code at the end of the populate function normalizes the polynomial values. This ensures that that for large values of N the errors are minimized during image reconstruction. We can see that the polynomials are normalized from the graph, as the zeroth order polynomial is not 1.0, but actually ~0.125.
Okay... So how do we use these polynomials?
Well, given a image NxN pixels in size, we first need to compute (or retrieve) a TchebichefPolynomial of size N. Then we iterate variables p and q from 0 through to N-1, and compute the Tchebichef moment Tpq as given by the equation above. This essentially is done by iterating through the image pixels, and for each pixel, multiplying the image intensity with the value at x and y for Tchebichef polynomials of orders p and q, and summing those together, to return a double. Once we've done this for all orders up to N, we effectively have an image of moments. We can get the original image back by doing the same thing, but instead of multiplying the image intensity, use the Tchebichef moment computed for the given x and y.