APPLICATION OF GENETIC ALGORITHMS
IN
VECTOR QUANTIZATION
by
Gina Leyba Cone
B.S. Colorado State University, 1991
A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado at Denver
in partial fulfillment
of the requirement for the degree of
Master of Science
Electrical Engineering
1995
A
This thesis for the Master of Science
degree by
Gina Leyba Cone
has been approved for the
Department of
Electrical Engineering
by
l. John Clark
Shelly Goggin
Date
Cone, Gina Leyba (M.S., Electrical Engineering)
Application of Genetic Algorithms in Vector Quantization
Thesis directed by Professor Tamal Bose
ABSTRACT
There has been much research recently into the use of Genetic Algorithms
(GAs) for complex, largescale problems. GAs are adaptive search techniques
based on evolutionary rules. This thesis investigates GAs as a method for the
Vector Quantization (VQ) problem of image compression.
The VQ design problem and the wellknown LBG algorithm for VQ are
reviewed. A brief history and the basic components of GAs are described with a
few variants of the algorithm explained.
Implementation of a VQGA is illustrated and performance is tracked for
two different test cases. An approximate fitness function is utilized in the
experimental design, which reduces the computational effort. Performance of the
VQGA design is analyzed and compared with LBG algorithm performance.
Conclusions are presented about the feasibility of VQGA design for practical use
and items for future research are suggested.
in
This abstract accurately represents the content of the candidate's thesis,
recommend its publication.
Signed
Tamal Bose
iv
CONTENTS
Chapter
1.Introduction ............................................................1
2. Vector Quantization....................................................5
2.1 LBG Algorithm.......................................................11
3. Genetic Algorithms....................................................13
3.2 Other Methods........................................................14
3.3 GA Basics..........................................................15
3.3.1 A GA Example......................................................17
3.4 GA Theory..........................................................22
3.5 GA Variations......................................................30
4. Vector Quantization Using Genetic Algorithms .........................33
4.1 Experimental Results................................................35
5. Conclusion............................................................64
Appendix
A. VQGA Program.........................................................66
B. VQLBG Program .......................................................88
References ...............................................................99
v
FIGURES
Figure
1.1 Basic Information Coding Model....................................... .2
1.2 Shannon's Cascaded Coding Model........................................3
2.1 Partitioning of Vectors................................................7
2.2 VQ Encoder and Decoder.................................................8
2.3 LBG VQ Design Algorithm...............................................11
3.1 Roulette Wheel........................................................20
3.2 Crossover Operation ..................................................21
3.3 Mutation Operation....................................................21
3.4 Schema Example .......................................................23
3.5 Hyperplanes for a 3bit String ...................................... 24
3.6 Population of Strings.................................................25
3.7 Crossover Disruption of Schema........................................28
4.1 Example ft\ Training Sequence.........................................38
4.2 Distortion Performance, Example #1,1/16 Sampling Rate ................39
4.3 Distortion Performance, Example #1,1/8 Sampling Rate..........40
4.4 Distortion Performance, Example #1,1/4 Sampling Rate..........41
4.5 Distortion Performance, Example #\,100% Sampling Rate ................42
4.6 Normalized Distortions, Example #1,1/4 Sampling Rate..........44
vi
4.7 Normalized Distortions, Example ^1,1/16 Sampling Rate ...........45
4.8 LBG Distortion Performance, Example #1............................46
4.9 Example Â§2 Training Sequence......................................47
4.10 Distortion Performance, Example #2,1/16 Sampling Rate ...........48
4.11 Distortion Performance, Example 2,1/8 Sampling Rate..............49
4.12 Performance Distortion, Example #2,1/4 Sampling Rate..............50
4.13 Distortion Performance, Example #2,100% Sampling Rate ............51
4.14 LBG Distortion Performance, Example ffl...........................52
4.15 Example #1,a) Original Images b) Images Compressed by LBG codebook
54
4.16 Example #2, a) Original Images b) Images Compressed by LBG Codebook
55
4.17 Example 1Images Compressed by VQGA Codebook................56
4.18 Example it\, Images Compressed by VQGA Codebook................57
4.19 Example #2, Images Compressed by VQGA Codebook................58
4.20 Example #2, Images Compressed by VQGA Codebook................59
4.21 Fixed Computation Sampling Comparison, Example it\ ...............61
4.22 Fixed Computation Sampling Comparison, Example #2.................62
vn
1.Introduction
The need for digital images has exploded recently due to the evolution of
digital processing, storage, and transmission techniques. Digital images have
many advantages such as accurate transmission quality, reduced complexity and
flexibility in processing, and the capability of storage on computers as opposed to
some type of analog medium, such as video tape. Digital images can be easily
obtained by taking horizontal and vertical samples (raster scan) across the analog
image. Unfortunately, the number of samples required to achieve a quality image
can be overwhelming for storage or transmission. Here lies the need for image
compression.
Image compression is the mapping of the original image into a coded image
such that the reconstructed image has the best possible fidelity for the available
storage or transmission capability [5]. The coded image is represented by fewer
bits than the original image representation. Compression can be a lossless process
where the original digital image can be reconstructed exactly, or a lossy process
where the reconstructed image will be have a loss of fidelity.
The basic information coding model, as shown in Figure 1.1, converts
information to a form that is suitable for transmission over a channel. The goal of
the encoder is to code the information in an efficient manner for accurate
1
Reconstructed
Image
Codes
Figure 1.1 Basic Information Coding Model
transmission. For digital images, the information could be a rectangular block of
pixels that represent the color of the image at each discrete point in space. The
information could also be in a form where each image pixel has been transformed
to another domain, such as a Discrete Fourier Transform (DFT) of the image.
The choice of image representation depends on the type of processing that will be
done since each form has certain advantages. Whether the image is of a pixel
based or frequencybased form, the image is just a collection of numbers to be
transmitted through a channel. Shannon introduced the coding system that
separates the coding into two different functions, as shown in Figure 1.2. The
encoding function compresses the original image with an irreversible process. The
entropy coder compacts the compressed image and transforms the information into
codes appropriate for transmission (or storage).
2
Image Symbols 
Encoder 1
Original
Image
Reconstructed
Image
Symbols Entropy ^ Codes 1
Decoder
Figure 1.2 Shannon's Cascaded Coding Model
There are several different image compression methods. These methods
fall into two basic categories: lossless and lossy. Lossless methods work with a
digitizedimageandtreattheimageasastreamofbits(sameasatextfile).
Variable length codes are used for the different pixel values. Any of the well
known entropy coding methods, such as Huffman, LizZempel, and Arithmetic
coding, take advantage of redundancies that occur. The simplest form of this type
of compression is Run Length Encoding (RLE). This method simply looks for
strings of pixels with the same values and codes them as a block that consists of a
count with the pixel value. For example, a string of 8 pixels with the value 188
would be coded as (8 188). This method works very well with computer
generated graphics since these type of images contains many strings of identical
pixel values. The other entropy compression schemes assign symbols with the
xauceq3
3
least number of bits to the pixels with the highest probabilities. These entropy
coding methods are easily implemented and are commonly used. Graphic file
formats such PCX, TIFF, TARGA use these lossless methods as a standard.
Typical compressions ratios for an 8bit grayscale phototype image can be as
high as 4:1[5].
The other image compression technique, lossless, introduces a irreversible
transform of information. The goal of a lossy compression technique is to
compress an image while keeping a certain fidelity or quality. This method can
only be used if the original image is not needed. If one is able to accept a
decrease in fidelity, very high compression ratios can be achieved.
Commonly used lossy techniques include [5][7]:
1) Scalar quantizers with entropy coding (PCM)
2) Scalar quantizers with predictive coding (DPCM)
3) Transform and subband coding (Discrete Transform
Coding, DCT)
4) Wavelets
5) Fractals
6) Vector Quantizers (VQs)
Typical compression ratios, using some of these methods, for an 8bit grayscale
phototype image can be from 4:1 to as high as 32:1.
4
2. Vector Quantization
The VQ compression technique is derived straight from the second part of
Shannon's Rate Distortion Theory which states that better performance can always
be achieved by coding vectors instead of scalars [7]. Most of the current image
compression methods used are based on quantizing individual pixels of images.
These methods produce samples that are highly correlated or dependent and are
suboptimal methods [14]. VQ has the potential to exploit the dependence of the
scalars within a vector [11].This allows for a bit rate of less than one bit per
pixel,a compression rate that cannot be achieved by scalar quantization. The
improvement in performance does come at a computational price. Some of the
other potential advantages of VQ are fast decompression times and easy
incorporation of image enhancement techniques and classification methods for
pattern recognition.
A VQ, or multidimensional quantizer, is a system that maps a sequence of
vectors into symbols. This pap>er will focus on memoryless VQ systems for the
purpose of image compression. A VQ, 0 is a mapping of input vectors, x, from
the kdimensional Euclidian space, /f4, to a finite set, C of vectors called a
codebook [4], that is, for
5
X Rk ,
(2.1)
Q Rk .
(2.2)
The codebook consists of code vectors y that can be represented by symbols
of indexes. For example, if the codebook consists of 4096 vectors, then the /rt
vector can be represented by a 12bit symbol i instead of the actual vector. This
coding technique assumes that the receiving end of the transmission will have the
codebook.
The quantizer is completely described by the codebook and the selection of
the appropriate code vector for an input vector, called partitioning [4]. This
process is shown in Figure 2.1, where the input vector space is partitioned into
groups /P with each group corresponding to a code vector. The process of
determining whether a given x vector belongs in a particular region /P depends on
the distortion measurement used. This is a measurement of dissimilarity between
the input vector and the code vector [11].There are many different distortion
measurements that can be used, but the most commonly used one is the mean
squared error (MSE)
6
(=E,_2 (23)
The MSE measurement is commonly used because it is mathematically simple to
compute. The overall performance of a VO is measured by the average distortion
i M
D = M^d{X^}
The VQ encoding and decoding model is shown in Figure 2.2. The
encoder receives a vector and produces an index /. The decoder receives the index
and looks up the Ith output vector in the codebook.
7
Encoder
x
codebook vectors y T,
! yL
Decoder
Figure 2.2 VQ Encoder and Decoder
Shannon, however, did not disclose any techniques for the actual design of
these VQ systems. Until Lloyd extended an algorithm used for scalar quantization
design to the design of VQ in the 1970s, not much practical work had been
presented [12]. Since then, there have been variations of Lloyd's simple algorithm
presented such the LBG or kmeans algorithm. Application of VQs to actual data
sources has been limited to low bitrate systems due to the complexity of the
mathematics and the large computational effort required to design a VQ and
encode the images [1].Within the past few years, however, increased computer
speeds have made VQ applications more feasible.
The design of an optimal VQ is a problem of satisfying two conditions:
8
1) An optimal quantizer must choose the ccxle vector yt which yields
the smallest distortion measurement for a given vector v(. This is the
nearest neighbor condition where
yQi'v) iff d(v,yt)
2) Each code vector must be chosen such that the average distortion D
is minimized. For a given cluster, y, is the centroid of that cluster.
Each of these conditions can easily be solved independently. If the partitions, or
clusters, are known, the code vectors will be the centroids. If the code vectors are
known, the clusters are found by using the distortion measurement.
The design of a VQ becomes a nonlinear problem when neither the
partitions nor the codebook exists and the two conditions must be solved jointly.
An iterative approach must be used that starts with all possible vectors of k
dimension v and an initial codebook. First, the clusters are determined from
initial ccxle vectors using condition #\ and then the new code vectors are
determined using condition #2. This results in a new codebook and an average
distortion that, by definition, is less than or equal to the previous value [4]. This
process is repeated until an optimal codebook is produced.
Unfortunately, this process has two elements that are impractical to
implement. First, the set of all possible vectors v is an extremely large set of
9
vectors, especially for higher dimensions. Secondly, the computation of the
centroids requires that the statistical distribution p(v) be known. A modification of
these requirements is to use a training sequence. The training sequence is a set of
training images several times larger than the size of the codebook to be designed.
Typically the size of the training sequence is at least 10 to 50 times larger than the
codebook size and should be a suitable representation of the class of images to be
compressed with the VQ [11].
One design issue is the selection of an initial codebook. This becomes a
problem itself since the initial codebook determines how well the final codebook
will be. There are several different techniques for obtaining an initial codebook.
The simplest technique is to randomly pick vectors, but this does not usually yield
a useful codebook [11].Another technique employs a splitting algorithm that
starts with one vector that is the centroid of all training vectors and then splits into
two vectors by adding a small perturbation vector [1].This process is repeated
until the desired number of vectors are produced. Another method used is a
"pruning" method which travels through the training sequence and picks training
vectors that are dissimilar.
10
2.1 LBG Algorithm
The LBG algorithm for VQ design is shown in the flow diagram in Figure
2.3. The clustering algorithm classifies each of the M training vectors into one of
the code vectors yt by finding the vector that yields the smallest distortion
(condition #1). The new code vectors are determined by condition #2, where the
average distortion is minimized. If the distortion measurement is the MSE in
Equation 2.4, then the new code vectors are the average of all the vectors in each
cluster. This simple calculation is the reason why MSE is widely used. This
process is repeated until the average distortion is small enough, usually tested by
D D
e
D
(2.6)
where e is a small value and D., is the average distortion of the last iteration. It
has been shown by Linde, Buzo, and Gray that this algorithm will converge to a
local minimum in D, given a particular initial codebook [12]. For a global
minimum, the approach is to use several different initial codebooks and choose the
VQ that produces the smallest average distortion.
11
1 ^ i $ L
>
Quantize M training vectors
to L clusters
Detennioe new y by
computing centriod
of each cluster
Signifii
Changei
inD?
No
t
Done
Figure 2.3 LBG VQ Design Algorithm
The quantization of the training vectors involves the majority of the
computations. The distortion measurement is computed L times for each of the M
training vectors, or ML number of times. The MSE calculations take k arithmetic
operations (k additions and k multiplications). The total number of arithmetic
operations for each iteration is MLK. For a codebook size of 1024, M = 20L,
and vector dimension of 16 (4x4 blocks), this grows to over 300 million
operations. The quantization of a given vector involves NL arithmetic operations.
This determines the speed of the encoder. The decoder, however, is simply a
table lookup operation and does not involve any computations.
12
3. Genetic Algorithms
A genetic algorithm is a search method based on the science of natural
selection where the fittest creatures tend to survive. In nature, the ability of a
creature to adapt to its changing environment is crucial for survival. Individual
characteristics are defined by an individual's genes, which are grouped into
chromosomes [17]. Characteristics have been observed in nature that have an
optimal design, simply produced by the evolution process. If an optimization
problem can be converted to a gene (or string) and subjected to the rules of
evolution where the fittest creatures are selected for reproduction, pairs of genes
recombine to create new offspring, and parts of genes are sometimes mutated, then
the continuous evolution of these genes should result in an optimal gene after
several generations.
Genetic Algorithms (GAs) were first proposed by John Holland in 1975
who designed computer programs to mimic the adaptive behavior of natural
selection [10]. There had been other evolutionary strategies which used mutation
as the primary operator for the introduction of new material. Holland determined
that the crossover function was the primary operator for introducing new material
and mutation should only be used to preserve diversity [10].
13
3.1 Other Methods
There are three basic optimization schemes: calculusbased, enumerative,
and random [15]. Calculusbased techniques typically are descent methods based
on the gradient of the objective function. These methods rely on the computation
of derivatives which may be very complex and computationally intensive to find.
Derivatives also require continuity which may not exist in many problems. These
methods also lack robustness in some cases. For multioptima problems, hill
climbing will only find the local optima, or peak. The LBG algorithm for vector
quantization is this type of method. Enumerative techniques are basically
exhaustive searches which are very inefficient for problems with any moderately
sized search space. Random search methods are variations of the enumerative
techniques. In summary, the traditional methods do not work well for problems
with multioptima, large search spaces, and complex functions.
GAs are different from traditional optimization methods in many ways.
GAs use codings of the problem variables, called strings, instead of the variables
themselves. GAs work with a group of these strings at one time which produces
many points in the multimodal search space. This is a contrast to the hillclimbing
techniques which use a single point. Other methods that transition from one point
to another point within a local region will only locate the local peak. Since GAs
work in parallel in several regions of the search space, they have a greater
possibility of finding the optimal peak.
14
GAs use very little problemspecific information so they are very simple to
implement and easy to integrate into existing models. Calculusbased techniques
require many calculations. The only calculation GAs need to compute is a fitness
evaluation for each string, which is typically the objective function. In other
words, GAs perform blindly.
GAs' search direction is guided by probabilistic rules as opposed to
deterministic rules used in other methods. These type of rules choose solutions that
improve the performance but may not be in the same local region as the previous
solution. These characteristics make GAs very robust and attractive for complex
problems with many optima.
3.2 GA Basics
GAs have four basic components:
1) Evaluation The performance for each string
2) Reproduction The selection of individuals for the next
generation
3) Crossover The pairing of individual strings and exchange of
substrings
4) Mutation The alteration of individual bits in a string
First, the optimization problem must be defined. A nonlinear problem may
15
contain several variables, x,, x2, x3,...., that are not independent of each other.
Each variable has a defined domain of values, or search space. The goal of the
optimization problem is to find the combination of variable values that maximizes
(or minimizes) the performance. GA uses strings which are encodings of the
variables in the problem. Typically this is a string of bits, although there has been
some use of nonbinary alphabets [18]. This paper will assume a binary alphabet.
Each variable must be confined to a finite number of values, or a finite search
space. For example, a sample encoding of a variable that can take on integers
from 0 to 31 would be a 5bit string where the variable value would be equal to
the binary value. For a continuous variable, the search space must be discretized
and the number of discrete values must be a power of 2.
A fitness function/W is needed to evaluate the performance of a particular
string. The fitness evaluation may be the exact function to be optimized, or it may
be an approximation. The computation time of the fitness evaluation must be very
short since the GA requires evaluations for each string in the population.
The process of GA starts with the creation of an initial population. This
may be a random generation of binary strings, or could be seeded fully or partially
with initial material[10]. The population of strings is subjected to three
operations: reproduction, crossover, and mutation.
Reproduction is the process of selecting the best strings in a particular
population for survival into the next generation. The Roulette Wheel Selection
16
Scheme [5is the simplest stochastic method. Each string is allocated a wedge on
a wheel, as shown in Fig 3.1, with a size proportional to the value of f(x). A
string is selected by a spin of the roulette wheel and added to the new population
of strings.
The next step is the crossover process. This process takes pairs of strings
and exchanges substrings. This is analogous to gene recombination. The strings
are randomly mated and made available for crossover. A crossover point along
the strings is randomly generated. The crossover is only performed for each pair
only if the crossover probability pc is met. Then the substrings are exchanged to
produce a new pair of strings for the population.
The third step is mutation. The simplest scheme is to flip each bit with a
certain mutation probability pmy usually a very small fixed value. If a good
substring is destroyed through either reproduction or crossover, mutation may be
the only opportunity of that string appearing again in the population.
3.2.1 A GA Example
To demonstrate each step in the GA, here is a simple example for the
problem of maximizing the function f(x,y) = x2 + y2 where the domain of x,y =
(0,31). The variables can be encoded as a combination of two substrings, each 5
bits long. A population size of N=8 will be used.
First, an initial population is obtained by random generation. The initial
17
population for this example is shown in Table 3.1.A fitness evaluation for this
problem is obtained by inserting the integer value from each substring into the
f(x,y) equation.
In the reproduction process, the probability of selection is
Pselect
fsum
(3.1)
whereis the fitness value for the Ith string, and/s/w is the sum of all the
fitness values. The values for pselea become the proportional sizes for each wedge in
the Roulette Wheel scheme, as shown in Figure 3.1. The spin of the wheel is
simulated by generating a random fractional number between 0 and 1,and
comparing the cumulative totals on the outside of the wheel to determine where the
random fraction lies. The corresponding string is added to the new population. The
larger wedges have a greater probability of being selected, while the smaller
wedges have a lesser probability. These "spins" are continued until the population
is full. This is a stochastic selection method where the actual number of copies of
each string may differ than the expected number. These numbers are shown in
Table 3.1.
18
61
3iduiBX3 v re aiq^l
Sad
906Â£ wns
9Â£t 0010101100 duon 0010101100 9 9 0011101100
L6 0010010010 6 0110010010 t z 0010110010
LS9 0001110010 8 001111010 9 8 OOUJIIOOIO
06Z 1000010001 3U0D 1000010001 E t 1101011001
Z8t noiOTiooi z 1101011011 Â£ S 1000010011
869 1011011101 L 1010011101 8 l oil nioi
6Z6 001011HOT soon 0010111101 V L 0TI00TIT0I
LK 0111011010 V 0111010010 8 Â£ 1011010010
WJ uoijBindoj 8mpU3 SUOllISOJ J3AOSSOJ3 uoijBndd J3AOSSOJ3 3 aoipnpojdsy jv uo ndoj
0*1 SZV 8.6 S(w
0_8 00*1 86Sf wns
I Z8*1 IT 0Z8 St9 0011101100
Z L0 1 EI_ m 06 0010110010
T A0* I cr zst n6i TIOIOIIOOI
0 A0. m OI'ZI 0101000110
I 6e*i LY 9Z9 V9Z 1000010011
z 9T\ 9I_ S9S 9Kil onoomoi
0 6Z W on ITÂ£ 1101011000
I 95' L0' 09Z Â£16 1011010010
junojjjy jurtojdx^ w/ uoijB[ndd IHiJiui
0.51
Figure 3.1 Roulette Wheel
For the crossover step, pairs of strings are mated randomly. A random
fractional number is generated and must be less than the crossover probability pc in
order for crossover to occur for each pair. If so, then the crossover point is
chosen by generating a random number between 1 and Ll, where L is the length
of the strings. This crossover point divides the strings into two substrings. The
substrings are exchanged and a new pair of strings replace the "parent" strings.
This crossover process for one pair of strings is shown in Fig 3.2. The crossover
process for the example in Table 3.1, where pc = l and crossover is always
performed. Multipoint crossover points can also be used [18].
20
Old Pair
New Pair
10011 010 110 1 00 11 0000 1
001011 i 001 0010110110
Figure 3.2 Crossover Operation
The mutation step is a very simple process. For each bit, a random
fractional number is generated and compared against the probability of mutation
pm. If the number is less than pm, then that particular bit is flipped. The value of pm
must be very small, otherwise the disruption to the population is too great and the
positive effects of reproduction and crossover processes are lost. For example, this
process is shown in Figure 3.3 for. a /?= 0.3, where the arrows point to bits that
have been selected for mutation.
Original____ ______Mutated
1011010110 1111000100
t t t
Figure 3.3 Mutation Operation
After an entire GA cycle, a new population is obtained. Note that in Table
3.1 the average fitness value has increased. Although this is only one iteration of
the GA, this is a typical result of GA. It has been proven that GAs search in the
direction of improvement, and are not just random searches [5].
21
One generation of a simple GA has been demonstrated. These generations,
or iterations, would be continued until a stopping criteria is met. The criteria may
be a fixed number of iterations, when a certain threshold fitness value is reached,
or when the strings have converged such that all the strings in the population are
almost identical.
3.3 GA Theory
In the previous section, one generation cycle was illustrated with favorable
results, However, a more rigorous theoretical analysis is needed. Most of the
basic work that has been presented is based on Holland's work [6]. This section
will cover the sampling of the search space and Holland's Schema Theorem. This
analysis should answer the fundamental question of how the use of GA operators
lead to a robust search.
The sampling of the search space can be demonstrated using similarity
templates called schemata. Schemata are special strings that define the similarity
between a set of strings. For example, Figure 3.4 shows a set of four different
strings of length L=5. The template for these strings is f* 7 0), where the
represents a "don't care" symbol and all other positions are fixed. The number
of fixed positions in a schema is called the order o(H). The schema in Figure 3.4
is an order 2 schema.
22
Strings: 0 0 110
0 1110
11110
0 1100
Schema: **1*0
Figure 3.4 Schema Example
Schemata can also be visualized as hyperplanes in a Ldimensional search
space, where L is the string, or schema length. For string length L=3, the search
space is threedimensional.A cube shown in Fig 3.5, represents all the possible
strings and schemata. The corners of the cube are the actual bit strings where
adjacent corners contain strings that differ only in one position. Each point is also
a schema of order 3. Each line in the cube is an order 2 schema and each face of
the cube is an order 1 schema.
23
x2
x
Figure 3.5 Hyperplanes for a 3bit String
Here are some other definitions needed for further investigation of
schemata:
A: Population of individual strings
A(t): Population at time t
H: Schema
a,: String Allele (Particular position in a string or schema)
Defining length of Schema (Distance between the outermost fixed position)
24
o(H): Order of Schema (Number of fixed positions)
m(H,t): Number of appearances of schema H with population Aft)
Schemata can also be used to analyze how GAs guide a search for better
strings. If a group of high fitness strings has some similarity, perhaps that
similarity is a characteristic of the optimal string. For example, the population of
strings shown in Fig 3.6 is ordered by the fitness values. Notice that there are
similarities between the three strings with the highest fitness values. The schema
for these three strings is (1 * 11 * 1 0 *). If this schema contains fixed
positions that are part of the optimal string, then this schema should be preserved
throughout the GA operations. The opposite argument can be made for strings
with low fitness values. It would not be desirable to preserve these schemata
throughout the GA operations. This concept of maintaining schemata from high
fitness strings and discarding schemata from low fitness strings is the basis for
GAs.
Figure 3.6 Population of Strings
J o 3 8 p 5 o
^7 7 9 2^
9 8 6 4 2 c
<1 tl
 o o o 1 o 1
25
Schemata also add another layer of data to the GA search. The only
information GAs have to work with are the actual strings and their fitness values.
If similarities between strings are also considered, then GAs have more
information in the search for an optimal string. For strings of length L, the total
number of possible schemata is (2+l)L = 3L. Consider a Nsized population. For
each individual string, there are 2L possible schemata. Therefore, the population
has from T to NT schemata, where the lower boundry is the case where all strings
are identical and the upper boundry is the case where all strings are dissimilar.
Here is a stepbystep analysis of the schema theorem first proposed by
Holland [5]. The schema theorem considers the effects of reproduction, crossover,
and mutation on schemata from generation to generation and quantifies the growth
(or decay) of the schemata.
At the beginning of a generation, the population A(t) contains m instances
of a particular schema H, which is written as m(H,t). After reproduction, there
are m(H,t+l) instances of a schema H in the population. A string is selected
during the reproduction process, based on its fitness, with the probability of
The probability that schema H is selected during the reproduction process is
(3.2)
26
Psr
N,
fsum
(3.3)
where/(7/J is the average fitness of string in population A(t) that contains schema
H. Substituting the average string fitness/avg = fsum/N, the number of schemata
to survive the reproduction process is
From this equation, it is observed that schemata that exist in highly fit strings are
more likely to be selected for reproduction.
The effects of crossover on schemata are not as easy to analyze as the
reproduction process. Since mates are chosen randomly, crossover points are
chosen randomly, and the decision to perform crossover is a random probability,
an exact survival probability cannot be derived but a lower boundary for schema
growth can be determined.
For illustration, consider two different schemata H,, and H2 and string A
shown in Figure 3.7. Suppose the onepoint crossover point chosen is 6, chosen
from points 1 to Li.fhe schema H2 will survive the crossover process only if the
mate for string A is identical, which is very unlikely. The probability of schema
H being disrupted in the crossover process is higher for H2 than H, since H2 has a
larger defining length. The probability of the destruction of schema H, is
AH)
(3.4)
favg
27
Pile
8(H{)
(1)
(3.5)
For H,, pdc = 0.25, and for H2, pdc = 0.725.
Figure 3.7 Crossover Disruption of Schema
The probability of crossover survival for schema H would be
Psc
Pc
8(H)
LI
(3.6)
where pc is the probability of crossover occurring. The inequality accounts for
the possibility of the creation of schema H through the crossover process.
The mutation process alters each allele of a string with probability pm. The
probability that a single position of a schema survives is The number of
positions in a schema is o(H). Since each mutation is statistically independent, the
probability of survival for schema H is
Psm  (! Pj
pm
(3.7)
where the total probability is equal to the multiplication of each position
//,'
28
probability.
The schema growth rate is obtained by the multiplication of the survival
schemata from reproduction processthe crossover survival probability, and the
mutation survival probability. This total generation growth is
m(H,t^\) > ~ Pc^](l ~ Pm){H) (3.8)
Schemata with high fitness values and short defining lengths should have a
high growth rate, observing the effects of reproduction and crossover. If mutation
is also considered, then schemata should also have a low order. These types of
schemata are the building blocks of Holland's Building Block Hypothesis. This
hypothesis proposes that the combination of good building blocks through
reproduction and crossover should result in strings with high fitness values.
Mutation is simply seen as a device for adding new building blocks into the
population. There have been numerous experimental trials that support this
hypothesis, but a rigorous proof has not yet been derived [6].
Recall that there are T to N21 schemata present in a population of size N.
A number of these schemata are destroyed by the crossover process, due to long
defining lengths. Holland estimated that N3 schemata are usefully processed by
GAs [5]. This can also be viewed as high rate hyperplane sampling where several
hyperplanes are evaluated simultaneously for each string. Holland called this
implicit parallelism, where several schemata are simultaneously evaluated each
29
time a single string is evaluated [18].
3.4 GA Variations
Since Holland's introduction of the simple GA (demonstrated in a previous
section), there have been many variations of the algorithm and its operators.
These modifications are attempts to improve the performance and robustness of
GAs.
The selection scheme for the reproduction process is based on proportionate
value of the fitness value to the average fitness value. Strings with high fitness
values relative to other strings, tend to dominate the selection scheme. This is
usually the case in the initial generation, as shown in the previous example.
Figure 3.1 illustrated the ratio of fitness values where four of the strings have high
fitness values compared to the other four strings. Strings that have extremely high
fitness value could take over the population. The opposite problem arises when
the variance in fitness values becomes small. The selection scheme now allocates
one copy of each string and the difference in fitness values is essentially ignored.
Scaling of the fitness values attempts to overcome both of these problems. There
are several types of fitness scaling. Linear scaling uses a firstorder equation
f'{x) af(x) + b, (3_9)
where f'(x) is the adjusted fitness. The coefficients a and b must be chosen so
30
that the maximum fitness is 1.52.0 times the average fitness [17]. Unfortunately,
this could result in negative adjusted fitness values if the fitness values are far
below the average fitness value. This problem is overcome in the Sigma
Truncation Scheme, where strings with fitness values below the standard deviation,
a, are discarded [17]. The adjusted fitness value for f(x) > a are calculated by
f\x) = f(x) (favg ca), (310)
where c is a small constant usually in the range from 1.0 to 3.0.
The selection of the crossover and mutation probability values is an
optimization problem itself. The crossover probability determines the rate of
recombination of building blocks. If the value is too high, then there is too much
disruption of high fitness strings. If the value is too low, then the belowaverage
fitness strings will not be modified. The mutation probability determines the rate
that material is reintroduced into the population. If the value is too high, then
high fitness strings are destroyed and the GA turns into a random search. If the
value is too low, then there is not enough diversity to find the optimal string.
Dynamically varying these parameters is an attempt to achieve optimal
performance for each generation and prevent premature convergence. One
adaptive scheme by Srinivas and Patnaik disrupts strings that have a belowaverage
fitness [16]. Crossover and mutation probabilities for each string are based on the
difference between the string's fitness and the maximum fitness in the population.
31
The adaptive crossover probability is
, (frnax fL) "f
Pc ^ kl {fmax favg) {/~^8) (3.U)
:3 if
where fL is the larger of the pairs' fitness values and where k k} < 1.0.. If
kl =kJ, a string with belowaverage fitness is given a higher probability that those
with above average fitness and high fitness strings are given the lowest
probabilities. A similar equation for the mutation probabilities is
_ k iPnax /)
m 2 {fmax favg)'
=4
(f>favg)
(/)
where k2, k4 < 1.0.
(3.12)
32
4. Vector Quantization Using Genetic Algorithms
This chapter will present the optimization of a vector quantization problem
using a genetic algorithm.
Two GA components must be defined. First, the problem variables must be
encoded to a binary string. Recall from Chapter 2 that a codebook is defined by
the collection of the vectors. Each vector could be represented by a binary string
of pixel intensities, where each pixel is an 8bit binary number. For a 4pixel
vector, the vector string would be 32 bits in length. For a 32vector codebook,
the GA string would be 1024 bits.
Next, a fitness function must be defined to evaluate each string in the
population. The performance of.a codebook can be measured by the mean squared
error, also called the average distortion. This measurement would have to be
computed for all the codebooks in the population, at each generation. Another
option would be to use an approximate fitness evaluation which would take a
fraction of the computational time and still be effective in determining
performance [8]. This approach takes advantage of the capability of GAs to
perform in noisy environments [3].
Approximate functions have been successfully utilized in GAs in the
problem of image registration in the field of medical diagnostics [2][13]. The
33
image registration performance evaluation is similar to the MSE measurement in
the VQ problem. The performance of a particular image registration solution is
measured by the calculation of the average of the differences between the two
different images. An approximate measurement can be obtained by using a sample
of the differences between the two images and calculating the average. Two
different techniques have been employed to determine the sampling. One method
is a simple random selection of a fixed number of pixels for sampling [2]. The
other selection method uses a stratification technique where the image is divided
into uniform sub regions, or blocks, and a maximum of one pixel per block is
selected for sampling [13]. Pixels with a maximum variance within a block are
chosen for sampling if the variance meets a fixed threshold.
The same stratification sampling technique can be used in the VQGA
problem. The training sequence is divided into equal blocks of vectors. The
variance of each vector is computed within each block and the vector with the
maximum variance is added to the sampling set if the variance is greater than a
fixed threshold. The approximate fitness function is the sum of the sampled vector
differences divided by the number of samples, or
f(x) = Z5 = j2d(xry,). (41)
n U
The same set of vectors would be used for each evaluation of a string in the
34
genetic algorithm operations.
4.1 Experimental Results
VQ design problems using a GA were implemented for two different
training sequences. Two performance experiments were tracked for each of the
VQ design problems. A vector size of 4 pixels (2x2 block) and a codebook size
of 32 vectors were chosen for both problems. This size of codebook creates a VQ
problem with a search space size of 21024, which is a fairly small size for VQ, but a
rather large GA search space. The fitness function for each string, or codebook, is
the approximate distortion measurement described in the previous section. GA
performance was tracked for a fixed number of generations, while varying the
number of samples used in the approximate evaluations. The other set of
experiments fixes the number of total evaluation computations by varying both the
number of generations and the number of samples.
Since VQ optimization is a minimization problem instead of a maximization
problem, the fitness function is actually the inverse of the approximate distortion.
The linear scaling scheme in Eq 3.9, with C^ = 2.0, was employed to help
prevent premature convergence. The constants a and b have to be computed for
each generation. If finin > (C^ favg fmax)I{C^ 1.0) then
35
(4.2)
(Cw g
fmax fmin
b favg*(frnax ~ favg*CmuIt)
fmax fmin
Otherwise,
a
favg
favg fmin
fmin favg
fmin favg
(4.3)
This method will always produce positive fitness values [5]. The mutation
probability was determined by the adaptive scheme given in Eq 3.12, with
parameters k2 = .001 and k4 = .05. The crossover probability was set to a fixed
value of pc = 0.7 after initial runs using an adaptive probability did not
demonstrate an improvement over using a fixed value.
A population seeded with vectors sampled from the actual training
sequence was used as the initial population. Since the search space is extremely
large, a randomly generated initial population starts the performance of the
algorithm at a very high average distortion and requires a larger number of
generations for convergence. Use of a seeded initial population gives the VQGA
a similar starting point as in the LBG algorithm. An initial population was
36
obtained by selecting vectors that had a distortion difference greater than a fixed
threshold. The same initial population was utilized for each run.
The first experiment used the training sequence, shown in Fig 4.1,
consisting of three different images with a total size of 45132 vectors. The
population size was fixed at a value of N = 50. The C program code is listed in
Appendix A. The performance progress of this algorithm was tracked for four
different approximate distortion sampling rates, with each run iterated for 100
generations. The actual distortion measurement and the approximate distortion
measurement were computed after each generation and the best values (minimal
distortion) from the population are plotted in Figures 4.2 4.5. Note that even
though the approximate distortion is greater than the actual distortion, there
appears to be some correlation between the behavior of the two measurements.
This behavior could be a result of the use of fitness values in GAs where only the
relative value of a fitness function is important. The Roulette Scheme, shown in
Fig. 3.1, is one GA operation example of this characteristic. Each GA run
displays a distortion which has a very gradual decrease. This is an illustration that
the GA is working, although very slowly. One observation is that the experimental
runs that sample the distortion measurement at a lower rate display a noisier
distortion track than the runs that sample at a higher rate. This was an expected
result since less accurate measurements could cause good strings, or codebooks, to
be thrown out and bad strings to be kept in the population.
37
Figuxje 4.1 Example 1 Training Sequence
38
enl^A swell11
39
Figure 4.2 Distortion Performance, Example #1,1/16 Sampling Rate
SMaLml
40
Figure 4.3 Distortion Performance, Example #1,1/8 Sampling Rate
41
Figure 4.4 Distortion Performance, Example #1,1/4 Sampling Rate
1600
6
1500
1400
1300
CD 3 1200
15
> O) 1100
v
1000
900
800
700
600
 . Exact A g Distor tion
0 10 20 30 40 50 60 70 80 90 100
Generations
Figure 4.5 Distortion Performance, Example #1,100% Sampling Rate
The normalized distortions for 1/4 sampling rate and the 1/16 sampling rate
are plotted is Figure 4.6 and Figure 4.7, where each distortion series is normalized
by the maximum value in the series. The normalized approximate and actual
distortion tracks have a very similar pattern for the 1/4 sampling rate, while the
normalized tracks for the 1/16 sampling rate are not as similar. The interesting
observation is that a noisier distortion track did not necessarily result in worse
performance. Notice that the sampling rate that achieved the minimal distortion
(best performance) was not the 100% sampling rate, but the 1/8 sampling rate.
The noisy environment may be preventing a premature convergence. In addition,
the processing time was 1/8 the processing time of the run that used 100%
sampling rate. This demonstrates the power that GAs have to perform well in
noisy environments.
The LBG algorithm was used against the same training sequence for a
comparison of performance. The same procedure for obtaining an initial
population for the GA was used for obtaining an initial codebook for the LBG
algorithm. The C program code is listed in Appendix B. The algorithm was run
for 30 iterations and the performance results are plotted in Figure 4.8. Notice that
the LBG has a much smoother descent and outperforms the GA for this test case.
43
eeA sssu'^pezllEEJON
44
Figure 4.6 Normalized Distortions, Example #1,1/4 Sampling Rate
0.75H
0 10 20 30 40 50 60 70 80 90 100
Generations
Approx Distortion Actual Distortion
Figure 4.7 Normalized Distortions, Example #1,1/16 Sampling Rate
1200
1100
1000
900
800
700
600
500
400
0
5
10
15
Iterations
20
25
30
uqtolsIaCTAV
46
Figure 4.8 LBG Distortion Performance, Example #1
This series of experiments was repeated for the training sequence, shown in
Figure 4.9, with a size of 31912 vectors. The performance tracks for four
different distortion sampling rates are plotted in Figures 4.10 4.13 using a
population size of N = 20. Again, there is a gradual decrease in the average
distortion with noisier tracks for the higher sampling runs. The LBG algorithms
results for the same training sequence are plotted in Figure 4.14 and again
outperforms any of the GA runs.
Figure 4.9 Example #2 Training.Sequence
47
1100
3
Id
>
B
1000
900
800
700
600
500
400
300
A
^ Approx \yg Dist >rtion
r~^T Exact A\ g Distorl ion
0 10 20 30 40 50 60 70 80 90 100
Generations
Figure 4.10 Distortion Performance, Example #2,1/16 Sampling Rate
1100
to
(/)
a)
1000
900
800
700
600
500
400
300
0 10 20 30 40 50 60 70 80 90 100
Generations
_VV"
Approx
\vg Dist >rtion
Exact A
to Distor tion
Figure 4.11 Distortion Performance, Example #2,1/8 Sampling Rate
1100
1000
900
800
700
600
500 =
400
300I
0 10 20 30 40 50 60 70 80 90 100
Generations
Figure 4.12 Performance Distortion, Example #2,1/4 Sampling Rate
1100
1000
900
800
700
600
500
400
300
Exact Avg Disl A ortion '
0 10 20 30 40 50 60 70 80 90 100
GGn^reitio s
Figure 4.13 Distortion Performance, Example #2,100% Sampling Rate
1100
1000
900
800
700
600
500
400
300
0
2
4
6 8 10
Iterations
12 14
16
uo'olsjQCTAV
52
Figure 4.14 LBG Distortion Performance, Example #2
Two different test images were used for a comparison of compression
quality. Figures 4.15 and 4.16 display the original images along with compressed
images using the two codebooks designed from the LBG algorithm. Figures 4.17
4.20 display compressed images using the two codebooks designed from the VQ
GA for four different sampling rates. The LBG codebooks and the VQGA
codebooks have the same compression ratio of 32:5 for both images. The actual
average distortion for each compressed image is also displayed. The LBG
codebooks almost always compressed images with a better quality than the VQGA
codebooks. This is a direct effect of compressing images with a codebook that
had a lower average distortion in the design process.
53
Original
LBG, AvgD = 343
a) b)
Figure 4.15 Example #1,a) Original Images b) Images Compressed by LBG
Codebook
a) b)
Figure 4.16 Example #2, a) Original Images b) Images Compressed by LBG
Codebook
55
6. 2S% Sampling
AvgD = 522
12.5 % Sampling
AvgD = 545
25% Sampling
AvgD =517
Figure 4.17 Example #1,Images Compressed by VQGA Ccxlebook
56
6.25% Sampling
AvgD = bJb
25% Sampling
AvgD = 849
100% Sampling
AvgD = 809
Figure 4.18 Example #1,Images Compressed by VQGA Codebook
57
6.25% Sampling
AvgD 582
12.5% Sampling
AvgD = 549
100% Sampling
AvgD =02
Figure 4.19 Example #2,Images Compressed by VQGA Codebook
58
6.25% Sampling
AvgD =762
12.5% Sampling
AvgD =5D
25% Sampling
AvgD = 727
100% Sampling
AvgD =73
Figure 4.20 Example #2, Images Compressed by VQGA Codebook
59
The next experiment was a comparison of performance using different
sampling rates while varying the number of generations. This was an investigation
to determine the best sampling rate, given a fixed computation time. The
population size was left as a fixed size, but the sampling rate and the number of
generations were restricted by fixing the total number of distortion calculations in
each run. For example, for a 500,000 total number of distortion calculations
using a population size of N =10, would it better to sample 1000 vectors each
fitness evaluation for 50 generations or to sample 50 vectors each fitness
evaluation for 1000 generations? Figure 4.21 displays the final actual distortion
for eight different sampling rates using the first training sequence and a fixed
population size of N = 50. The distortion values are an average of six runs. The
total number of distortion computations is held fixed at a value of ten times the
number of vectors in the training sequence, with a fixed population size of
N = 50. The best performance comes from the experimental runs that utilized a
lower sampling rate, while the worst performance comes from the run which used
100% sampling. The same experiment is carried out for the second training
sequence. The results are shown in Figure 4.22, with less dramatic results. The
lower sampling runs have a better performance than the higher sampling runs, but
several of the sampling rates have similar performances. In any case, the choice
would be to use a lower sampling rate and not sacrifice any performance. For
these two test cases, GAs seem to actually perform better in noisier environments.
60
(Generations)
(10)
6.25 9.26 12.5 18.5 25 38.5 50 100
Avg Distortion Sampling Percent
o 5
9 8
0(5(0(
8 7 7
UOIJJ01nBAVIBUI
/V /l
5 o
6 6
61
Figure 4.21 Fixed Computation Sampling Comparison, Example #1
(Generations)
6.25 9.26 12.5 18.5 25 38.5 50 100
Avg Distortion Sampling Percent
o
60
550
o
50
450
8
350
o
30
uoltolsKlsAV BUI
62
Figure 4.22 Fixed Computation Sampling Comparison, Example #2
Processing times for the GA are much longer than for the LBG algorithm.
The VQGA codebook designed for the first training sequence at a 1/8 sampling
rate took on the order of 50 times longer to compute as the LBG codebook. The
processing time for the second training sequence was 25 times longer than the
LBG processing time.
63
5. Conclusions
The test cases of a VQGA design illustrated some of the advantages and
disadvantages of GAs. One advantage is that the GA required very little problem
specific information for implementation. Only the code vector variables and the
distortion measurement were used in the GA. In contrast, the LBG algorithm is a
calculusbased method which could require complex calculations in the
determination of the centroids, depending on the performance measurement. Since
the MSE was used in the test cases, the LBG centroid computation was simple.
However, if a different measurement is used such as the Holder norm distortion,
the centroid computation may become mathematically very complex in the LBG
lgorithm [5.A different performance measurement would only change the
evaluation process in the GA. The simple GA operations would not have to be
modified.
One big disadvantage of the GA was the long processing times. The use
of an approximate fitness function did help in reducing the processing time,
without sacrificing performance. For the first test case, the run using a 1/8
sampling rate took 1/8 of the processing time of the run using a 100% sampling
rate and had a lower distortion. Unfortunately, the processing time was still much
longer than the LBG processing times. The GA processing time could also be
64
reduced further by using parallel processing in the fitness evaluation since each
string evaluation is independent.
The use of an approximate fitness function in both series of experiments
demonstrated the advantage GAs have in noisy environments. In fact, the sampled
fitness function seemed to perform better than the actual fitness function. The
noisy environment may actually prevent the algorithm from getting stuck at local
minima.
The VQGA, however, produced suboptimal codebooks, compared with
the LBG algorithm codebooks. GAs are supposed to have the ability to find global
minima, as opposed to the LBG algorithm which will only find the local minima.
VQ is a largescale problem which may exhibit a high level of epistasis, or
nonlinearity. Simple GAs may not perform optimally for these type of problems.
More adaptive GA components or other variants may have to be implemented.
In summary, these experimental results show some promise for the use of
GA in VQ design. However, the computational effort has to be reduced before a
practical implementation would be feasible. Further study is needed in the area of
premature convergence and additional GA components. Hybrid algorithms that
combine a GA and a conventional VQ method are also another possibility for
investigation.
65
Appendix A VQGA Program
/* This version uses a GA to design a VQ codebook. */
/* This distortion samples vectors by the stratifcation method */
/* by using the vector with the top pv variance for each subblock */
/* Linear scaling is used. */
/* Adaptive mutation propability is used */
include
^include
include
include
include
^include
^include
^include
< string. h>
< malloc.h
^define sbsize 2
^define bsize 2
^define VSIZE 4
#defme CSIZE 32
^define S 50
^define ga_iter 26
^define pc 0.9
^define k2 0.001
^define k4 0.05
^define pv 0.28
^define Cmult 2.0
^define FI 2
#define DI 2
Global Constants*************/
/* Sampling Block size */
/* Block size */
/* Vector size */
/* Codebook size */
/* Population Size */
I* ft oi generations */
/* probability of crossover threshold */
/* probability constant */
/* probability constant */
/* variance threshold */
/* Linear scaling variable */
/* Output file interval */
/* Display interval */
unsigned long xsizel, ysizel, nxsizel, nysizel, bsize 1;
unsigned long xsize2, ysize2, nxsize2, nysize2, bsize2;
unsigned long xsize3, ysize3, nxsize3, nysize3, bsize3;
/* image sizes */
66
unsigned long m; /* Number of training blocks*/
unsigned long R; /* Number of Sampled blocks*/
unsigned int bits; /* Bits in a pixel */
unsigned long BSIZE=VSIZE*CSIZE; /* Pixels in a codebook */
typedef struct / /* GA Array*/
i float fx[S]; float fxs[S]; float pselect[S]; int actct [S]; int rflag[S]; float a; float b; float fsum; float favg; float fmin; float fmax; unsigned long fmin_index; float fmaxs; float fsums; float favgs; float fminm; float fsumm; float favgm; float fmaxm; } GASTRUCT; /* Actual fitness */ /* Scaled fitness */ /* Prob of repro selection */ /* Copies for repro */ /* Repro process flag */ /* Linear scale constant */ 1* Linear scale constant */ /* Total fitness sum */ /* Fitness avg */ /* Best string fitness */ /* Worst string fitness */ /* Best string index */ /* Max scalef fitness */ /* Scaled fitness sum */ /* Scaled fitness avg */ /* Prescale min fitness */ /* Prescale fitness sum */ /* Prescale fitness avg */ /* Prescale max fitness */
typedef struct { char manufacturer; char version; char encoding; char bitsperpixel; int xmin; int ymin; int xmax; int ymax; int hres; int vres; char colormap[48]; char reserved; /* PCX File Format */
67
char nplanes;
int bytesperline;
int palettetype;
char filler[58];
char palette[768];
} PCXHEADER;
PCXHEADER pcxhdr;
int readpcxblock( int *, int *, FILE *);
unsigned long distortion (unsigned char huge *, unsigned char huge *);
float prob_xover (GA_STRUCT *, int int);
int readhdr( FILE *pcxfile)
Reads an entire file and stores it in a (large) buffer,
pointed to by the variable bufr
Function from Tricks of the Graphics Gurus
)!(/
pcxhdr. version = getc(pcxfile);
pcxndr_ encoding = getc(pcxfiie);
pcxndr.bitsperpixel=fgetc(pcxfile);
pcxhdr. xmin = getw(pcxfile);
pcxndr.ymin = getw(pcxfile);
pcxhdr. xmax = getw(pcxfile);
pcxhdr. ymax = getw(pcxfile);
pcxhdr. hres = getw(pcxfile);
pcxhdr.vres = getw(pcxfile);
fread(pcxhdr.colormap, 48, sizeof(char), pcxfile);
pcxhdr. reserved = getc(pcxfile);
pcxhdr. nplanes = getc(pcxfile);
pcxhdr. bytesperline = getw(pcxfile);
pcxhdr. palettetype = getw(pcxfile);
fread(pcxhdr.filler, 58, sizeof(char), pcxfile);
pcxhdr. manufacturer = getc(pcxnle);
if (pcxhdr. manufacturer !=10)
{ fclose(pcxfile);
retum(O);
68
retum(l);
unsigned long readpcx(unsigned char huge *buf, FILE *pcxfile, unsigned long
*linesize)
/})()((*(]^}e})(})(}](*
Read image pixel info rrom pcx file into buffer
Function from Tncks of the Graphics Gurus
{ unsigned int 1;
int a;
unsigned long 1=0;
unsigned char huge *bufr;
int blockl=0 countl=0;
unsigned char block =1,count =1;
bufr = buf;
for (1=0;1 < *linesize; ) { /* increment by count below */
a = readpcxblock(&blockl, &countl, pcxfile);
block = block l&0x00ff;
count = countl;
if (a = = EOF) break;
for (i = 0; i < count; i + +) {
(*bufr++) = block;
}
1 += count;
}
a = readpcxblock(&blockl, Sccountl, pcxfile);
block = block l&OxOOff;
count = countl;
fclose(pcxnle);
retum(l);
int readpcxblock( int *pblock, int *pcount, FILE *pcxfile)
/})Â£)e))<
Read one encoded block from the pcx
file and store a count and data byte.
If the code is greater than OxcO (192),
then read count. Otherwise, read a byte
Function from Tncks of the Graphics Gurus
})(**/
69
{ intj;
*pcount =1;
f (0 = getc(pcxfile)) == EOF) return(EOF);
if (OxcO = = (OxcO & j)) /* is it greater than 192? */
{ *pcount = Ox3f&j; /* subtract 192 to get count */
if((j = getc(pcxfile)) == EOF) return(EOF);
}
*pblock = j;
return(l);
int imagetovector( unsigned long nxsize, unsigned long nysize, unsigned long
xsize, unsigned char huge *buff, unsigned char huge *mbuff)
Convert pixel buffer to vector buffer
{ unsigned long r,n,d,i,j,k;
unsigned char huge *mbufr;
mbufr = mbuff;
for (k=0; k<(nxsize)*(nysize); k++) {
r = fmod(k,(nxsize));
d = kr;
n = d*VSIZE +r*bsize;
for (i=0; icbsize; i + +)
for (j=0; j
*mbufr++ = *(buff + (unsigned long)(n + i +j*xsize));
return (1);
int imagetoblock( unsigned long nxsize, unsigned long nysize, unsigned long xsize,
unsigned char huge *buff, unsigned char huge *mbuff)
Convert vector buffer to block buffer
)}f(/
{ unsigned long r,h,d,i,j,k,l,n;
unsigned char huge *mbutr;
mbufr = mbuff;
70
for (k=0; k<(nysize/sbsize); k++) {
for (h=0; h<(nxsize/sbsize); h++) {
d = k*sbsize*xsize*bsize + h*sbsize*bsize;
for (j=0; j
for (i=0; i
for (n=0; ncbsize; n++) {
for (1=0; lcbsize; 1 + +) {
*mbufr++ = *(buff + (unsigned long)(d + i*bsize +j*xsize*bsize
+ 1 + n*xsize));
return (1);
unsigned long findsamples(unsigned char huge *mbuff, unsigned long huge
*sample, unsigned long nxsize, unsigned long b, unsigned long offset)
Find vectors to be used in sampling distortion measurement
{ unsigned long row, column,block,var_index=0;
unsigned long k=0;
unsigned long i,j,l;
double mean[VS IZE], mean_mag, temp;
double var, var_max=0,var_total,var_per;
block = sbsize*sbsize;
for (i=0; i
for=0;j
mean[j] = 0;
for (j =0;j
for (I=0;1
mean[l]+= *(mbuff + (unsigned long)(VSIZE*blcx:k*i + j*VSIZE
+ D)
mean_mag = 0;
for (j=0;j
mean[j] /= block;
71
mean_mag += mean(j]*mean[j];
var_max = 0;
var_total=0;
for (j=0;j< block; j + +) {
var = 0;
for (1=0;1
temp = *(mbuff + (unsigned long)(VSIZE*block*i + j*VSIZE + 1));
if (temp = 0)
temp = fabs((pow(temp,2))mean[l]*mean[l]);
else
temp = mean ]*mean[l];
var += temp;
}
if(var > var_max) {
var_max = var;
var_index = j;
}
var_total += var;
}
if (var_total==0.0 )
var_per = 0;
else
varper = var_max/var_total;
if (var_per > = pv) {
row = sbsize*(i/(nxsize/sbsize)) + (var_index/sbsize);
column = sbsize*(i%(nxsize/sbsize)) + var_index%sbsize;
*(sample + k) = row*nxsize + column + offset;
k+ + ;
retum(k);
int copyvector(unsigned char huge *vl, unsigned char huge *v2)
{ int i;
for (i = 0; i< VSIZE; i + +)
*(v2 + i) = *(vl + i);
retum(l);
72
int copycodebook(unsigned char huge *codel, unsigned char huge *code2)
{ unsigned long i,j,k;
for (j = j< CSIZE; j + +)
for (k = 0; k< VSIZE; k++)
*(code2 + j*VSIZE + k) = *(codel + j*VSIZE + k);
retum(l);
int writebyte(unsigned char byt, FILE *pcxfile)
/* Write byte to pcx file */
{
if(EOF = = putc((int)byt, pcxfile))
retum(O); /* disk write error (probably full)*/
else
return
int popvalues (GA_STRUCT *ga, unsigned char huge *pop, unsigned char huge
*vl, unsigned long huge *sample)
/)c))e)le*
Compute GA pertormance values (distortion: avg, min, max, sum
}K)(}c/
{ unsigned long j,k,l;
unsigned long l;
float avgd,d,dt;
float fx;
ga>fsum=0;
ga>fmax=0;
ga>fmin=99999;
for(i=0; i
avgd = 0;
for(j=0; j
1=*(sample + j);
dt = distortion(pop + i^BSIZE, vl + PVSIZE);
for (k=l; k
if ((d=distortion(pop + i*BSIZE + k*VSIZE,vl + 1*VSIZE))<= dt){
dt = d;
73
avgd += dt;
avgd /= R* 1.00;
ga>fx[i] = avgd;
ga>fsum += ga>fx[i];
if (ga > fx[i] < ga > fmin) {
ga>fmin = ga>fx[i];
ga> fmin_index = i;
}
if (ga > fx[i] > ga > fmax) {
ga>fmax = ga>fx[i];
Compute fitness inverse and fitness avg, sum, max, min for scaling process
{ int i;
ga>fsumm=0;
ga>fmaxm=0;
ga > fminm=9999999;
for(i=0; i
ga>fx[i]=1.0/ga>fx[i];
if (ga > fx[i] < ga > fminm) {
ga> fminm = ga>fx[i];
ga>favg = ga>fsum/(S*1.00);
return(l);
int prescale (GA_STRUCT *ga)
if (ga>fx[i] > ga>fmaxm) {
ga>fmaxm = ga>fx[i];
ga>fsumm += ga>fx[i];
ga>favgm = ga>fsumm/(S*1.00);
retum(l);
74
int findscalecoeff (GA_STRUCT *ga)
/})Â£!)e%})(})cj)c
Find a and b scale coefficients for linear scaling
)c))(})e}e})C}(/
ir (ga > fmaxm = = ga > fminm   ga > favgm = = ga > fminm) {
pnntf("Error in findscalecoeff tunction");
return(O);
}
else {
if (ga> fminm > ((Cmult*ga> favgm*1.00 ga>fmaxm)/(Cmult1.0))) {
ga>a = (Cmult l)*ga>favgm/(ga> fmaxm ga> fminm);
ga>b = ga> favgm *(ga> fmaxm
Cmult*1.00*ga>favgm)/(ga> fmaxm ga> fminm);
}
else {
ga>a = ga> favgm/(ga> favgm ga> fminm);
ga>b = (ga>fminm*ga>favgm)/(ga> favgm ga>rminm);
}
return(l);
int select (GA.STRUCT *ga)
Select string index using Roulette Wheel Scheme
)(){)(()(9)()(()(()(()(>(}((()C9C)()((/
{ int k;
float sum,rand_num;
k=0; sum=0;
rand_num = random(101)*l.00/100;
do {
sum += ga>pselect[k];
k+ + ;
}
while (sum < rand_num && k < S);
retum(kl);
75
int crossover3 (unsigned char huge *pop)
/
Penorm twopt xover at a vector joint
}je)(e}(()((})()c/
{ unsigned char buft;
int mate[S], xoverl[S],xover2[S];
int i j;
for (i=0; i
for (i=0; i
if (mate[i] < 0) {
do
mate[i] = random(S);
while(mate[mate[i]] > =0  mate[i] = =i);
mate[mate[i]] = i;
xoverl[i] = random(CSIZEl)+l;
xover2[i] = random(CSIZEl) + l;
if (xoverl[i] > xover2[i]) {
j = xover l[i];
xover l[i] = xover2[iJ;
xover2[i] = j;
}
xoverl[mate[i]] = xover l[i];
xover2[mate[i]] = xover2[i];
if ((random(101)*l.00/100) < pc) {
for (j=xoverl[i]*VSIZE; j
buff = *(pop + BSIZE*i +j);
*(pop + BSIZE*i +j) = *(pop + BSIZE*mate[i] +j);
*(pop + BSIZE*mate[i] +j) = buff;
retum(l);
float prob_mut (GA_STRUCT *ga, unsigned long int i)
/)!(
Compute adaptive mutation probability
)K/
76
{ float prob;
if ((ga>fx[i] < ga>favgm)   (ga>fmaxm == ga>favgm))
prob = k4;
else
prob = k2*(ga>fmaxm ga>fx[i]*1.00)/(ga>fmaxm ga>favgm);
retum(prob);
^u*_***_**_*_
Print out performance progress to output file
{ int i;
printf("sum= %6.If avg = %5.If min= %6.lf'.ga>fsum,ga>favg,ga > fmin);
retum(O);
int printfile (FILE *fileptr, GA_STRUCT *ga, unsigned long g)
Print out performance progress to output file
*}{)^}()(}!()e}({)K/
{ int i;
fprintf(fileptr,"\n%lu %5Af %6.If ",g,ga>favg,ga>fmin);
retum(O);
unsigned long distortion (unsigned char huge *vl, unsigned char huge vz)
/))c)e
ComDute distortion between two vectors
}tnK)((#/
{ unsigned long sum = 0;
unsigned long diff;
long d;
int i;
for (i=0; i
d = ((*(vl+i)&0x00ff) (*(v2+i)&0x00ff));
77
diff = d*d;
sum += diff;
}
return (sum);
/* The main program will read in 3 pcx files that will be the training
/* sequence and an initial population.
/* The genetic algorithm will iterate for the max it of iterations
/* The final best codebook will be written to the codebook output file. */
main()
GA_STRUCT newga;
unsigned char mask;
char filenamel[30] = "d:\/tcwin\/images\/imagel .pcx"
d: \/tc wi n\/i mages\/ image2. pcx"
d:\/tcwin\/images\/image3.pcx"
char filename2[30]
char filename3[30]
*1
char filename4[40]
char filename5 [40]
d: \/tcwin\/output\outg.";
d: \/tcwi n\/codes\/";
char filename6[40] = "d:\/tcwin\/ccxles\/codeg."
char filename7[40] = "d:\/tcwin\/ccxies\/popo.
char filename[20];
unsigned char huge *buff;
unsigned char huge *start;
unsigned char huge *mbuff;
unsigned char huge *mstart;
unsigned char huge *newpop;
unsigned char huge *cbuff;
unsigned long huge *sample;
unsigned long i, j, k;
unsigned long p;
unsigned long g;
unsigned long ming;
unsigned long d, dt;
int c;
float minfx;
float pm;
int index;
/* 1st training image */
/* 2nd training image */
/* 3rd training image
/* results file */
I* initial pop */
/* ending codebook */
/* ending population */
/* Input image buffer */
/* Input image buffer pointer */
/* Vector image buffer*/
/* Vector image buffer pointer*/
/* Population buffer */
/* Codebook buffer pointer */
/* Sampling index buffer */
/* Index variables */
/* Total # of vectors */
/* Iteration index */
/* Min avgd Iteration */
/* Distortion variables */
/* File input char */
/* Min fitness */
/* Mutation prob */
/* Repro selection index *!
78
int flag; /* Flag *1
char ch; /* Keyboard input char */
float avgd, last_avgd; /* Avg Distortion variables */
unsigned long picsizel, picsize2, picsize3,offset; /* Image pixel sizes */
unsigned long npicsizel, npicsize2, npicsize3; /* Image block sizes */
FILE *fileptr 1,*fileptr2,*fileptr4; /* File pointers */
FILE *fileptr3,*fileptr5,*fileptr6, *fileptr7;
time_t timerl,timer2;
struct tm *tblockl;
/* File pointers */
/* Timer variables */
printf("Enter initial population file:");
gets( filename);
strcat(filename5,filename);
printf("Enter output version #:");
gets(filename);
strcat(filename7,filename);
strcat(fi lename4, filename);
strcat(filename6, filename);
if (((((fileptrl=
((fileptr2 =
((fileptr3 =
((((fileptr5 =
((fileptr6=
((fileptr4 =
fopen(filenamel, "rb"))
fopen(filename2, "rb"))=
fopen(filename3, "rb")):
fopen(filename5, "rb")):
fopen(filename6, "wb"))
fopen(filename4, "wt"))
NULL)  
NULL))  
NULL))  
NULL) j I
=NULL)))  
NULL))
printf(nError in file opening");
retum(O);
if ((fileptr7=fopen(filename7, "wb"))
printf("Error in file opening");
return (0);
NULL) {
/* Read image files and compute sizes */
if ((readhdr(fileptrl))!= 0) {
xsizel=pcxhdr.bytesperline;
ysizel=pcxhdr.ymax pcxhdr.ymin +1;
if ((fmod(xsizel,bsize)!=0)  (fmod(ysizel,bsize)! =0))
printf("\nPicture size 1 not divisible into blocks");
retum(O);
79
nxsizel=xsizel/bsize;
nysizel=ysizel/bsize;
picsizel=(long) (xsizel)* (pcxhdr.nplanes*l)* ysizel;
npicsizel=(long) (nxsizel)* (pcxhdr.nplanes*l)* nysizel;
bsizel=(long)(nxsizel/sbsize) (nysizel/sbsize);
}
else {
printf("Error in reading image header 1");
fclose(fileptrl); fclose(fileptr2); fclose(fileptr3);
return (0);
if ((readhdr(fileptr2))!= 0) {
xsize2 = pcxhdr.bytesperline;
ysize2 = pcxhdr.ymax pcxhdr.ymin +1;
if ((fmod(xsize2bsize)! =0)   (fmod(ysize2bsize)! =0)) {
printf(M\nPicture size2 not divisible into blocks");
retum(O);
}
nxsize2 = xsize2/bsize;
nysize2 = ysize2/bsize;
picsize2 =(long) (xsize2) (pcxhdr.nplanes*l)* ysize2;
npicsize2 =(long) (nxsize2) (pcxhdr.nplanes*l)* nysize2;
bsize2 =(long)(nxsize2/sbsize) (nysize2/sbsize);
}
else {
printf("Error in reading image header 2\n");
fclose(fileptrl); fclose(fileptr2); fclose(fileptr3);
return (0);
if ((readhdr(fileptr3))! = 0) {
xsize3 = pcxhdr.bytesperline;
ysize3 = pcxhdr.ymax pcxhdr.ymin +1;
if ((fmod(xsize3bsize)!=0)   (fmod(ysize3bsize)!=0)) {
printf(M\nPicture size3 not divisible into blocks");
return(O);
}
nxsize3 = xsize3/bsize;
nysize3 = ysize3/bsize;
picsize3 =(long) (xsize3) (pcxhdr.nplanes*l)* ysize3;
npicsize3 =(long) (nxsize3) (pcxhdr.nplanes*!) nysize3;
80
bsize3 =(long)(nxsize3/sbsize) (nysize3/sbsize);
}
else {
printf("Error in reading image header 3\n");
fclose(fileptrl); fclose(fileptr2); fclose(fileptr3);
return (0);
p = (unsigned long)(picsizel + picsize2 + picsize3 );
m = (unsigned long)(icsizel + icsize2 + npicsize3);
/*********** Mallbuffers *************/
(char huge *)farmalloc((p+l)*sizeof(char))) == (NULL))
if (((buff
  ((mbuff = (char huge *)farmalloc((p + l)*sizeof(char)))
printf("Error in buff malloc=\n");
fclose(fileptrl); fclose(fileptr2); fclose(fileptr3);
return(O);
(NULL)))
))()c/
/************* Mallocate sampling buffers
if ((sample = (unsigned long huge *)farmalloc((m+l)sizeof(long
= =(NULL)){
printf("Error in sample malloc=\n");
fclose(fileptrl); tclose(fileptr2); fclose(fileptri);
retum(O);
start = buff;
mstart = mbuff;
/*********** Read in training sequence *************/
if ((offset = readpcx(bufffileptrl&picsizel))==0)
retum(O);
buff = buff + (unsigned long)offset;
if ((offset =
return(O);
buff = buff
if ((offset =
retum(O);
buff = buff
readpcx(buff,fileptr2,&picsize2)) == 0)
+ (unsigned long)offset;
readpcx(buff,fileptr3,&picsize3))
+ (unsigned long)offset;
)
81
'***** Convert training sequence buffer to a sub block and find samples *****/
buff = start;
imagetoblock(nxsizelnysizel, xsizelbuff, mbuff);
buff = buff + (unsigned long)picsizel;
mbuff = mbuff + (unsigned long)bsizel*sbsize*sbsize*VSIZE;
imagetoblock(nxsize2, nysize2, xsize2, buff, mbuff);
buff = buff + (unsigned long)picsize2;
mbuff = mbuff + (unsigned long)bsize2*sbsize*sbsize*VSIZE;
imagetoblock(nxsize3, nysize3, xsize3, buff, mbuf;
/* Find vector samples for approx distortion */
mbuff = mstart;
if ((R = findsamples(mbuff,sample,nxsizel,bsizel,0)) < = 0) {
printf("\n No samples taken ");
retum(O);
mbuff = mbuff + bsizel*sbsize*sbsize*VSIZE;
if ((R + = findsamples(mbuff,sample + R,nxsize2,bsize2,npicsizel))< =0) {
printf("\11 No samples taken ");
retum(O);
mbuff = mbuff + bsize2*sbsize*sbsize*VSIZE;
if ((R += findsamples(mbuff,sample + R,nxsize2,bsize2,npicsizel+npicsize2))
< =0) {
printf("\n No samples taken ");
return(O);
buff = start;
mbuff = mstart;
/*******_*** Convert pixel sequence to vectors ****************/
imagetovector(nxsizel, nysizel, xsizel, buff, mbuff);
buff = buff + picsizel;
mbuff = mbuff + picsizel;
imagetovector(nxsize2, nysize2, xsize2, buff, mbuff);
buff = buff + picsize2;
mbuff = mbuff + picsize3;
82
imagetovector(nxsize3, nysize3, xsize3, buff, mbuff);
buff = start;
mbuff = mstart;
timerl=time(NULL);
farfree(start);
/****************** Mallow ga population array buffer *****************/
if ((newpop = (char huge *)farmalloc((CSIZE*VSIZE*S + l)*sizeof(char)))
==(NULL)){
printf("Error in maJloc\n");
fclose(fileptr 1);fclose(fileptr2); fclose(fileptr3);
return(O);
/* Use GA to find optimal codebook */
randomizeO;
g=
avgd = 0;
bits =(log(256))/log(2);
minfx = 999999;
ming=0;
while (g < ga_iter) {
if (g= =) { /* Starting generations */
/* Read in initial population file */
printf("reading pop file");
for (i=0; i
if ((c = getc(fileptr5)) = = EOF) {
printf("Error in reading codebook for i = %lu\n",i);
return(O);
}
else {
*(newpop+i) = c;
printf("\nS = %d pv=%3.3f C/V=%d/%d m=%lu R=%lu k2=%2.2f
k4=%2.2fv=%s\n",S,pv,CSIZE,VSIZE,m,R,k2,k4,filename);
popvalues(&newga, newpop, mbuff, sample);
printf("g=%lu:",g);
83
printga(&newga);
printfile(fileptr4&newgag);
/* Print out real avgd for min avg distortion codebook */
avgd = 0;
for (k=0; k
dt = distortion(newpop+newga.fmin_index*BSIZE, mbuff+k*VSIZE);
for (i = l; i
if ((d=distortion(newpop + newga.fmin_index*BSIZE + i*VSIZE,
mbuff + k*VSIZE)) < = dt) {
dt = d;
avgd += dt;
}
avgd /= m;
printf(" Avgd=%f',avgd);
fprintf(fileptr4," % 5.1 f" ,avgd);
/*********** Reproduction **********/
prescale(&newga);
findscalecoeff(&newga);
newga.fsums = 0;
newga.fmaxs = 9999999;
/******* Perform linear scaling process ********/
for(i=0; i
newga.fxs[i]=newga.a*newga.fx[i] + newga.b;
newga.fsums += newga.fxs[i];
if (newga.fxs[i] > newga.fmaxs) newga.fmaxs = newga.fxs[i];
}
/***** Compute probability of repro selection *******/
newga.favgs = newga.fsums/(S*1.00);
for(i=0; i
newga.actct[i] = 0;
newga.pselect[i] = newga.fxs[i]* 1.00/newga.fsums;
/********** Find actual # of copies for repro ************/
for(i=0; i
index = select(&newga);
newga.actct[index] + +;
84
newga.rflag[i]=0;
}
/******* Repro strings to intermediate population ********/
for(i=0; i
for(j =0;j +1 < (newga.actct[i]);j + +) {
flag=0;
for(k=0;k
if(newga.actct[k] = =0 && newga.rflag[k] = =0) {
newga.rflag[k = l;
copycodebook(newpop + i*BSIZE,newpop + k*BSIZE);
flag =1;
}
if (flag= = l)break;
/* Crossover Process*/
crossover3(newpop);
/* Mutation */
popvalues(&newga, newpop, mbuff, sample);
prescale(&newga);
for (i=0; iS; i + +) {
pm = prob_mut(&newga,i);
for (k=0; k
for (j=0; j< bits; j + +) {
if ((random(30001)*1.00 + random(30001)*1.00)/60000.00 < pm ) {
mask = *(newpop + i*BSIZE + k) < < (bitsj1);
mask > > = 7;
if (mask = = 0) {
mask + + ;
mask < < = j;
*(newpop + i*BSIZE + k) += mask;
}
else {
mask < < = j;
*(newpop + i*BSIZE + k) = mask;
85
g+ + ;
popvalues(&newga, newpop, mbuff, sample);
printfile(fileptr4,&newga,g);
/********** Compute actual distortion for output file interval **********/
avgd = 0;
if (g%(FI)==0) {
for (k=0; k
dt = distortion(newpop+newga.fmin_index*BSIZE, mbuff+k*VSIZE);
for (i = l; i
if ((d=distortion(newpop + newga.fmin_index*BSIZE + i*VSIZE,
mbuff + k*VSIZE)) < = dt) {
dt = d;
avgd += dt;
}
avgd /= m;
/****** Print results to results file ******
fpri ntf(fileptr4," % 5.1 f", avgd);
}
/******** Print results to display for display interval********/
if (g%(DI) = =0) {
printf("\ng=%lu:",g);
printga(&newga);
printf(" Avgd=%r,avgd);
timer2 = time(NULL);
printf("\nTotal Pr
/******* Find best codebook index ********/
for (i=0; i
if (newga.fx[i] == newga.fmin){
j=i
/* Write final codebook to file */
86
for (i=0; i
if(EOF == putc((int)*(newpop + j*BSIZE + i), fileptr6))
return(O); /* disk write error (probably full)*/
/* Write final population to file */
for (i=0; i CSIZE*VSIZE*S; i + +) {
if(EOF == putc((int)*(newpop + i), fileptr7))
retum(O); /* disk write error (probably full)*/
/***** **/
fclose(fileptr 1);fclose(fileptr2);fclose(fileptr3);close(fileptr4);
fclose(fileptr5); fclose (fileptr6); fclose(fi leptr7);
return(l);
87
Appendix B VQLBG Program
/* vq.c Turbo C program that finds a vq codebook given an initial */
/* codebook and a training sequence of images using the LBG */
/*.clustering algorithm and the MSE distortion measurement. */
^include
^include
include < string.h>
include
include < math.h>
include < time.h>
^include
/******************** Giobsl GonstBnsts
define bsize 2
^define vsize 4
define csize 32
^define epsilon 0.001
/* block size */
/* vector size */
/* codebook size */
/* Theshold value */
/*******V3.n3.blCS *******************/
unsigned long xsizel, vsize 1,nxsizel, nvsizei; /* image sizes */
unsigned long xsize2, ysize2, nxsize2, nysize2;
unsigned long xsize3, ysize3, nxsize3, nysize3;
unsigned long m; /* size of training blocks*/
typedef struct /* pcx file format */
{ char manufacturer;
char version;
char encoding;
char bitsperpixel;
int xmin;
int ymin;
int xmax;
int ymax;
int hres;
int vres;
char colormap[48];
char reserved;
88
char nplanes;
int bytesperline;
int palettetype;
char filler[58J;
char palette[768];
} PCXHEADER;
PCXHEADER pcxhdr;
/******* Program Functions ******/
int readpcxblock( int *int' *, FILE *);
/}e*K
Read an entire rue and store it in a (large) buffer,
pointed to by the variable butr
Function rrom Tricks of the Graphics Gurus
int readhdr( FILE *pcxfile)
pcxhdr. manufacturer = getc(pcxfile);
if (pcxhdr. manufacturer !=10)
{ fclose(pcxfile);
retum(O);
pcxhdr. version = getc(pcxfile);
pcxhdr. encoding = getc(pcxfile);
pcxhdr.bitsperpixel=fgetc(pcxfile);
pcxhdr.xmin = getw(pcxfile);
pcxhdr. ymin = getw(pcxfile);
pcxhdr. xmax = getw(pcxfile);
pcxhdr. ymax = getw(pcxfile);
pcxhdr. hres = getw(pcxfile);
pcxhdr. vres = getw(pcxfile);
fread(pcxhdr.colormap, 48, sizeof(char), pcxfile);
pcxhdr. reserved = getc(pcxnle);
pcxhdr. nplanes = getc(pcxfile);
pcxhdr.bytesperline = getw(pcxfile);
pcxhdr.palettetype = getw(pcxfile);
fread(pcxhdr.filler, 58, sizeof(char), pcxfile);
retum(l);
89
unsigned long readpcx(unsigned char huge *buf, FILE *pcxfile, unsigned long
*linesize)
/*
Reads lmaee mro from pcx rile format into bufier
Function from l'ncks of the Graphics Gurus
)(e>fcsfcsfeae)(c/
{ unsigned int 1;
int a;
unsigned long 1=0;
unsigned cnar huee *bufr;
int blockl=0 countl=0;
unsigned char block =1,count =1;
bufr = buf;
for (1=0;1 < *linesize; ) h increment by count below */
{ a = readpcxblock(&blockl, &countl, pcxfile);
block = block l&OxOOff;
count = countl;
if (a = = EOF) break;
for (i = 0; i < count; i + +) {
*bufr++ = block;
}
1 + = count;
}
a = readpcxblock(&blockl, &countl, pcxfile);
block = block li&OxOOff;
count = countl;
fclose(pcxfile);
retum(l);
int readpcxblock( int *pblock, int *pcount, FILE *pcxfile)
/^3((3{
Read one encoded block trom the pcx
file and store a count and data byte.
If the code is greater than OxcO (192), read count.
Otherwise, read a byte
Function from Tncks of the Graphics Gurus
_/
90
{ intj;
*pcount =1;
if (0 = getc(pcxfile)) = = EOF) return(EOF);
if (OxcO = = (OxcO & j)) /* is it greater than 192? */
{ *pcount = 0x3f&j; /* subtract 192 to get count */
if((j = getc(pcxfile)) == EOF) return(EOF);
}
*pblock = j;
retum(l);
int imagetovector( unsigned long nxsize, unsigned long nysize, unsigned long
xsize, unsigned char huge *buff, unsigned char huge *mbuff)
Transforms pixel buffer to vector buffer
{ unsigned long r,n,d,i,j,k;
unsigned char huge *mbufr;
mbufr = mbuff;
for (k=0; k<(nxsize)*(nysize); k++) {
r = fmod(k(nxsize));
d = kr;
n = d*vsize +r*bsize;
for (i=0; i
for (j=0; jcbsize; j + +) {
*mbufr++ = *(buff + n + i +j*xsize);
return (1);
int copyvector(unsigned char huge *vl, unsigned char huge *v2)
/)(
Copies vector vl to vector v2
/
for (1=0;i< vsize; i + +)
91
*(v2 + i) = *(vl + i);
return
unsigned long distortion (unsigned char huge *vl, unsigned char huge *v2)
/)(4{)c)e}(_}((})c*
Compute distortion between 2 vectors
}n)(})(}(>(e}K}je)(/
{ unsigned long sum = 0;
unsigned long diff;
long d;
for (i=0; iCvsize; i + +) {
d = ((*(vl+i)&0x00ff) (*(v2 + i)&OxOOf);
diff = d*d;
sum += diff;
}
return (sum);
int writebyte(unsigned char byt, FILE *pcxfile)
/}fcste})c})c)e})e}t(}e*
Writes bvte to pcx file
if(EOF == putc((int)byt, pcxfile))
retum(O); /* disk write error (probably full)*/
else
retum(l);
I* The main program will read in 3 pcx files that will be the training
/* sequence and an intitial codebook.
/* The LBG algorithm will iterate until threshold is met
/* or the max iterations is hit. The new codebook will be written to
/* the codebook output file. */
main()
92
char filename 1[30] = "c:\/tcwin\/images\/image 1.pcx";
char filename2[30] = "c:\/tcwin\/images\/image2.pcx";
char fiIename3[30] = "c:\/tcwin\/images\/image3.pcx";
char filename4[40] = "c:\/tcwin\/codes\/new.lbg\/";
char filename5[40] = "c:\/tcwin\/codes\/new.lbg\/";
char filename[20];
unsigned char huge *buff;
unsigned char huge *start;
unsigned char huge *mbuff;
unsigned char huge *mstart;
unsigned char huge *cstart;
unsigned char huge *cbuff;
unsigned long huge *mgroup;
unsigned long i, j, k;
unsigned long p;
unsigned long count;
unsigned long d, dt;
float test = 0.0;
int c;
int n,nt;
float avg[vsize], avgd, last_avgd;
unsigned long picsizel, picsize2, picsize3,offset;
/* 1st training image */
/* 2nd training image */
/* 3rd training image */
/* Initial codebook */
/* Output codebook */
/* filename buffer */
/* Input image buffer */
/* Input image buffer ptr */
/* Vector image buffer*/
/* Vector image buffer ptr */
/* Codebook buffer */
/* Codebook buffer ptr */
/* Clustering index buffer */
/* Index variables */
/* Total it of vectors */
/* Count for each cluster */
/* Distortion variables */
/* Threshold test */
/* File input char */
/* Iteration variables */
/* Avg Distortion values */
/*Image pixel size */
unsigned long npicsizel, npicsize2, npicsize3,npicsize4; /*Image vector size */
FILE *fileptrl,*fileptr2,*fileptr3,*fileptr4,*fileptr4,*fileptr6; /*File pts */
time_t timerl, timer2;
struct tm *tblockl;
/* Timer variables */
/* Timer variables */
printf("Enter initial codebook file:");
gets(fi lename);
strcat(filename4,filename);
printf("Enter codebook output file:");
gets(filename);
strcat(filename5, filename);
printf("Enter # of iterations:");
scanf("%d",&nt);
/* Open files */
if (((((fileptrl=fopen(filenamel, "rb")) == NULL)  
((fileptr2= fopen(filename2, "rb")) == NULL))  
((fileptr3= fopen(filename3, "rb")) == NULL))  
(((fileptr4= fopen(filename4, "rbM)) == NULL)  
93
