% File: boxcount.w [CWEB source code]
% Created: May 8, 2006 [v.1.0]
% Last change: November 26, 2011 [v.1.6]
% Author: Fredrik Jonsson
% Description: The CWEB source code for the BOXCOUNT simulator of nonlinear
% magneto-optical Bragg gratings. For information on the CWEB
% programming language, see http://www.literateprogramming.com.
% Compilation: Compile this program by using the enclosed Makefile, or use
% the blocks of the Makefile as listed in the documentation
% (file boxcount.ps or boxcount.pdf). The C source code (as
% generated from this CWEB code) conforms to the ANSI standard
% for the C programming language (which is equivalent to the
% ISO C90 standard for C).
%
% Copyright (C) 2006-2011, Fredrik Jonsson. Non-commercial copying welcome.
%
\input epsf
\def\version{1.6}
\def\lastrevdate{November 26, 2011}
\font\eightcmr=cmr8
\font\tensc=cmcsc10
\font\eightcmssq=cmssq8
\font\eightcmssqi=cmssqi8
\font\twentycmcsc=cmcsc10 at 20 truept
\def\boxcount{{\eightcmr BOXCOUNT\spacefactor1000}}
\def\poincare{{\eightcmr POINCARE\spacefactor1000}}
\def\ANSI{{\eightcmr ANSI\spacefactor1000}} % The language standard we stick to
\def\SI{{\eightcmr SI\spacefactor1000}} % Another standard for physical units
\def\GNU{{\eightcmr GNU\spacefactor1000}} % GNU is Not Unix
\def\GCC{{\eightcmr GCC\spacefactor1000}} % The GNU C-compiler
\def\CEE{{\eightcmr C\spacefactor1000}} % The C programming language
\def\ANSICEE{{\eightcmr ANSI~C\spacefactor1000}}% The ANSI-C standard language
\def\CWEB{{\eightcmr CWEB\spacefactor1000}} % The CWEB programming language
\def\MATLAB{{\eightcmr MATLAB\spacefactor1000}} % The MATLAB ditto
\def\UNIX{{\eightcmr UNIX\spacefactor1000}}
\def\CYGWIN{{\eightcmr CYGWIN\spacefactor1000}}
\def\CTANGLE{{\eightcmr CTANGLE\spacefactor1000}}
\def\CWEAVE{{\eightcmr CWEAVE\spacefactor1000}}
\def\OSX{{OS\,X}}
\def\re{\mathop{\rm Re}\nolimits} % The real part of a complex number
\def\im{\mathop{\rm Im}\nolimits} % The imaginary part of a complex number
\def\dollar{\char'044} % The `$' character
\def\tothepower{\char'136} % The `^' character
\def\onehalf{{\textstyle{{1}\over{2}}}}
\def\threefourth{{\textstyle{{3}\over{4}}}}
\def\endalg{\vrule height 1.4ex width .6ex depth .4ex} % Rectangular bullet
\def\ie{i.\thinspace{e.}~\ignorespaces}
\def\eg{e.\thinspace{g.}~\ignorespaces}
\let\,\thinspace
%% \def\subsection[#1]{\goodbreak\noindent{\it #1}\par\nobreak\noindent}
%
% Define a handy macro for listing the steps of an algorithm.
%
\newdimen\aitemindent \aitemindent=26pt
\newdimen\aitemleftskip \aitemleftskip=36pt
\def\aitem[#1]{\smallbreak\noindent\hbox to 10pt{}%
\hbox to\aitemindent{\bf #1\hfill}%
\hangindent\aitemleftskip\ignorespaces}
%
% Define a handy macro for bibliographic references.
%
\newdimen\refitemindent \refitemindent=18pt
\def\refitem[#1]{\smallbreak\noindent%
\hbox to\refitemindent{[#1]\hfill}%
\hangindent\refitemindent\ignorespaces}
\def\references{\medskip\goodbreak\noindent{\bf References}\nobreak\smallskip}
%
% Define a handy macro for nicely typeset variable descriptions.
%
\newdimen\varitemindent \varitemindent=100pt
\def\varitem[#1]#2{\smallbreak\noindent\hbox to 20pt{}%
\hbox to\varitemindent{#1\hfill}%
\hangindent 120pt\ignorespaces#2\par}
%
% Define a handy macro for nicely typeset descriptions of command line options.
%
\newdimen\optitemindent \optitemindent=80pt
\def\optitem[#1]#2{\smallbreak\noindent\hbox to\hsize{\hskip40pt{#1}\hfill}
\hangindent\optitemindent\ignorespaces{\hskip\optitemindent{#2}}\smallskip}
%
% Define a handy macro for the list of program revisions.
%
\newdimen\citemindent \citemindent=80pt
\newdimen\citemleftskip \citemleftskip=90pt
\def\citem[#1]{\smallbreak\noindent\hbox to 10pt{}%
\hbox to\citemindent{\bf #1\hfill}%
\hangindent\citemleftskip\ignorespaces}
\def\citindent{\smallbreak\noindent\hbox to 10pt{}%
\hbox to\citemindent{\hfil}%
\hangindent\citemleftskip\ignorespaces}
%
% Define the \beginvrulealign and \endvrulealign macros as described in
% Donald Knuth's The TeXbook, Appendix D: Dirty Tricks. These macros are
% used in typesetting nicely looking tables.
%
\def\beginvrulealign{\setbox0=\vbox\bgroup}
\def\endvrulealign{\egroup % now \box0 holds the entire alignment
\setbox0=\vbox{\setbox2=\hbox{\vrule height\ht0 depth\dp0 width0pt}
\unvbox0 \setbox0=\lastbox % now \box0 is the bottom row
\nointerlineskip \copy0 % put it back
\global\setbox1=\hbox{} % initialize box that will contain rules
\setbox4=\hbox{\unhbox0 % now open up the bottom row
\loop \skip0=\lastskip \unskip % remove tabskip glue
\advance\skip0 by-.4pt % rules are .4 pt wide
\divide\skip0 by 2
\global\setbox1=\hbox{\hskip\skip0\vrule\hskip\skip0
\unhbox2\unhbox1}%
\setbox2=\lastbox % remove alignment entry
\ifhbox2 \setbox2=\hbox{\kern\wd2}\repeat}}%
\hbox{\rlap{\box0}\box1}} % superimpose the alignment on the rules
%
% Make sure that the METAPOST logo can be loaded even in plain TeX.
%
\ifx\MP\undefined
\ifx\manfnt\undefined
\font\manfnt=logo10
\fi
\ifx\manfntsl\undefined
\font\manfntsl=logosl10
\fi
\def\MP{{\ifdim\fontdimen1\font>0pt \let\manfnt = \manfntsl \fi
{\manfnt META}\-{\manfnt POST}}\spacefactor1000 }%
\fi
\ifx\METAPOST\undefined \let\METAPOST=\MP \fi
\datethis
@*Introduction.
\vskip 100pt
\centerline{\twentycmcsc Boxcount}
\vskip 20pt
\centerline{A program for calculating box-counting estimates to the fractal
dimension of curves in the plane.}
\vskip 2pt
\centerline{(Version \version\ of \lastrevdate)}
\vskip 10pt
\centerline{Written by Fredrik Jonsson}
\vskip 10pt
\centerline{\epsfxsize=88mm\epsfbox{figures/figure-05.1}}
\vskip 10pt
\noindent
This \CWEB\footnote{${}^\dagger$}{For information on the \CWEB\ programming
language by Donald E.~Knuth, as well as samples of \CWEB\ programs, see
{\tt http://www-cs-faculty.stanford.edu/\~\ \kern -5pt knuth/cweb.html}.
For general information on literate programming, see
{\tt http://www.literateprogramming.com}.}
computer program calculates box-counting estimates of the fractal dimension
of curves in the two-dimensional plane.
In the box-counting estimate to the fractal dimension of a graph in the domain
$\{x,y:x_{\rm min}\le x\le x_{\rm max},y_{\rm min}\le y\le y_{\rm max}\}$, a
grid of squares, each of horizontal dimension $(x_{\rm max}-x_{\rm min})/2^m$
and vertical dimension $(y_{\rm max}-y_{\rm min})/2^m$, is superimposed onto
the graph for integer numbers $m$.
By counting the total number of such squares $N_m$ needed to cover the entire
graph at a given $m$ (hence the term ``box counting''), an estimate $D_m$ to
the fractal dimension $D$ (or Hausdorff dimension) is obtained as
$D_m=\ln(N_m)/\ln(2^m)$.
This procedure may be repeated may times, with $D_m\to D$ as $m\to\infty$
for real fractal sets. However, for finite-depth fractals (as generated by
a computer), some limit on $m$ is necessary in order to prevent trivial
convergence towards $D_m\to1$.
In addition to mere numerical calculation, the program also generates graphs
of the box distributions, in form of \METAPOST\ code which can be
post-processed by other programs.
\bigskip
\noindent Copyright \copyright\ Fredrik Jonsson, 2006--2011.
All rights reserved.
\vfill
@*The CWEB programming language.
For the reader who might not be familiar with the concept of the
\CWEB\ programming language, the following citations hopefully will
be useful. For further information, as well as freeware compilers for
compiling \CWEB\ source code, see {\tt http://www.literateprogramming.com}.
\bigskip
{\narrower\narrower\narrower\narrower\eightcmssqi\noindent
I believe that the time is ripe for significantly better documentation of
programs, and that we can best achieve this by considering programs to be
works of literature. Hence, my title: `Literate Programming.'
Let us change our traditional attitude to the construction of programs:
Instead of imagining that our main task is to instruct a computer what to
do, let us concentrate rather on explaining to human beings what we want
a computer to do.
The practitioner of literate programming can be regarded as an essayist,
whose main concern is with exposition and excellence of style. Such an
author, with thesaurus in hand, chooses the names of variables carefully
and explains what each variable means. He or she strives for a program
that is comprehensible because its concepts have been introduced in an
order that is best for human understanding, using a mixture of formal and
informal methods that reinforce each other.
\smallskip
{\eightcmssq --\,Donald Knuth,}
{\eightcmssqi The CWEB System of Structured Documentation}
{\eightcmssq (Addison-Wesley, Massachusetts, 1994)}
}
\bigskip
{\narrower\narrower\narrower\narrower\eightcmssqi\noindent
The philosophy behind CWEB is that an experienced system programmer, who
wants to provide the best possible documentation of his or her software
products, needs two things simultaneously: a language like \TeX\ for
formatting, and a language like C for programming. Neither type of language
can provide the best documentation by itself; but when both are appropriately
combined, we obtain a system that is much more useful than either language
separately.
The structure of a software program may be thought of as a `WEB' that is
made up of many interconnected pieces. To document such a program we want to
explain each individual part of the web and how it relates to its neighbors.
The typographic tools provided by \TeX\ give us an opportunity to explain the
local structure of each part by making that structure visible, and the
programming tools provided by languages like C make it possible for us to
specify the algorithms formally and unambiguously. By combining the two,
we can develop a style of programming that maximizes our ability to perceive
the structure of a complex piece of software, and at the same time the
documented programs can be mechanically translated into a working software
system that matches the documentation.
Besides providing a documentation tool, CWEB enhances the C language by
providing the ability to permute pieces of the program text, so that a
large system can be understood entirely in terms of small sections and
their local interrelationships. The CTANGLE program is so named because
it takes a given web and moves the sections from their web structure into
the order required by C; the advantage of programming in CWEB is that the
algorithms can be expressed in ``untangled'' form, with each section
explained separately. The CWEAVE program is so named because it takes a
given web and intertwines the \TeX\ and C portions contained in each
section, then it knits the whole fabric into a structured document.
\smallskip
{\eightcmssq --\,Donald Knuth, ``Literate Programming'', in}
{\eightcmssqi Literate Programming}
{\eightcmssq (CSLI Lecture Notes, Stanford, 1992)}
}
@*Revision history of the program.
\medskip
\citem[2006-05-08]{[v.1.0]} {\tt }\hfill\break
First properly working version of the \boxcount\ program, written in \CWEB\
and (\ANSI-conformant) \CEE. I have so far compiled the code with \GCC\ using
the {\tt --pedantic} option and verified that the box-covering routine
|get_num_covering_boxes_with_boxmaps()| and its more low-level core engine
|box_intersection()| both work properly, by direct inspection of the compiled
\METAPOST\ graphs generated by this routine.
That the numbers of counted boxes are correct has also been verified and the
only remaining blocks to be added are related to the extraction of the fractal
dimension as such.
This first version of the \boxcount\ program consists of 52671 bytes
(1292 lines) of \CWEB\ code, under \CYGWIN\ generating an executable file
of 22561 bytes and 33 pages of program documentation in 10~pt typeface.
\citem[2006-05-09]{[v.1.1]} {\tt }\hfill\break
Changed the block for the estimate of the fractal dimension, so that the
estimate now is obtained as the average of $\ln(N_{\rm boxes})/\ln(2^n)$
for a set of $n$ such that $N_{\rm min}\le n\le N_{\rm max}$, rather than
performing a linear curve fit to the data.
In order to sample statistical information on the estimate, such as standard
deviation, average deviation and skewness, I incorporated a slightly modified
version of the routine |moment()| from {\sl Numerical Recipes in C}.
Also added a preliminary section describing the test application of \boxcount\
to the Koch fractal, being a simple test case which is easily implemented by
means of a recursive generation of the data set for the input trajectory.
As of today, the \boxcount\ program consists of 59788 bytes
(1444 lines) of \CWEB\ code, under \CYGWIN\ generating an executable file
of 24253 bytes and 38 pages of program documentation in 10~pt typeface.
\citem[2006-05-14]{[v.1.2]} {\tt }\hfill\break
Added a section on the command-line options for supplying the \boxcount\
program with input parameters. Polished the section on the example of
estimation of the box-counting dimension of the Koch fractal, and in
particular changed the example from the Koch curve to the Koch snowflake
instead, just for the sake of visual beauty.
As of today, the \boxcount\ program consists of 72560 bytes (1699 lines) of
\CWEB\ code, under \OSX\ generating an executable file of 24996 bytes and 41
pages of program documentation in 10~pt typeface.
\citem[2006-05-17]{[v.1.3]} {\tt }\hfill\break
Added documentation on the |get_num_covering_boxes_with_boxmaps()| routine
and generally cleaned up in the documentation of the program.
As of today, the \boxcount\ program consists of 82251 bytes
(1844 lines) of \CWEB\ code, under \CYGWIN\ generating an executable file
of 29152 bytes and 40 pages of program documentation in 10~pt typeface.
\citem[2006-06-17]{[v.1.4]} {\tt }\hfill\break
Added the {\tt --graphsize} option, in order to override the default graph
size. Also changed the inclusion of the input trajectory in the boxmaps, so
that the \boxcount\ program rather than using a \MP\ call of the form
{\tt gdraw "input.dat"} now chops the trajectory into proper pieces and
includes the entire trajectory explicitly in the generated graph.
This way the output \MP\ code naturally increases considerably in size, but
is now at least self-sustaining even is separated from the original data file
containing the input trajectory.
\citem[2006-10-25]{[v.1.5]} {\tt }\hfill\break
Added two pages of text on the boxcounting estimate of the fractal dimension
of the Koch snowflake fractal in the example section.
\citem[2011-11-26]{[v.1.6]} {\tt }\hfill\break
Updated \.{Makefile}:s for the generation of figures. Also corrected a rather
stupid way of removing preceeding paths of file names.
@*Compiling the source code. The program is written in \CWEB, generating
\ANSICEE\ (ISO C90) conforming source code and documentation as plain
\TeX-source, and is to be compiled using the sequences as outlined in the
{\tt Makefile} listed below.
\bigskip
{\obeylines\obeyspaces\tt
\#
\# Makefile designed for use with ctangle, cweave, gcc, and plain TeX.
\#
\# Copyright (C) 2002-2011, Fredrik Jonsson
\#
\# The CTANGLE program converts a CWEB source document into a C program which
\# may be compiled in the usual way. The output file includes \#line specifica-
\# tions so that debugging can be done in terms of the CWEB source file.
\#
\# The CWEAVE program converts the same CWEB file into a TeX file that may be
\# formatted and printed in the usual way. It takes appropriate care of typo-
\# graphic details like page layout and the use of indentation, italics,
\# boldface, etc., and it supplies extensive cross-index information that it
\# gathers automatically.
\#
\# CWEB allows you to prepare a single document containing all the information
\# that is needed both to produce a compilable C program and to produce a well-
\# formatted document describing the program in as much detail as the writer
\# may desire. The user of CWEB ought to be familiar with TeX as well as C.
\#
PROJECT = boxcount
CTANGLE = ctangle
CWEAVE = cweave
CC = gcc
CCOPTS = -O2 -Wall -ansi -std=iso9899:1990 -pedantic
LNOPTS = -lm
TEX = tex
DVIPS = dvips
DVIPSOPT = -ta4 -D1200
METAPOST = mp
PS2PDF = ps2pdf
~ ~
all: \$(PROJECT) \$(PROJECT).pdf
~ @@echo "==============================================================="
~ @@echo " To verify the execution performance of the BOXCOUNT program"
~ @@echo " on the Koch snowflake fractal, run 'make verification'."
~ @@echo "==============================================================="
~ ~
\$(PROJECT): \$(PROJECT).o
~ \$(CC) \$(CCOPTS) -o \$(PROJECT) \$(PROJECT).o \$(LNOPTS)
~ ~
\$(PROJECT).o: \$(PROJECT).w
~ \$(CTANGLE) \$(PROJECT)
~ \$(CC) \$(CCOPTS) -c \$(PROJECT).c
~ ~
\$(PROJECT).pdf: \$(PROJECT).ps
~ \$(PS2PDF) \$(PROJECT).ps \$(PROJECT).pdf
~ ~
\$(PROJECT).ps: \$(PROJECT).dvi
~ \$(DVIPS) \$(DVIPSOPT) \$(PROJECT).dvi -o \$(PROJECT).ps
~ ~
\$(PROJECT).dvi: \$(PROJECT).w
~ @@make -C figures/
~ @@make -C kochxmpl/
~ @@make verification
~ \$(CWEAVE) \$(PROJECT)
~ \$(TEX) \$(PROJECT).tex
~ ~
verification:
~ @@echo "==============================================================="
~ @@echo " Verifying the performance of the \$(PROJECT) program on the Koch"
~ @@echo " snowflake fractal of iteration order 11."
~ @@echo "==============================================================="
~ make koch -C koch/
~ make kochgraphs -C koch/
~ make fractaldimension -C koch/
~ ~
tidy:
~ make -ik clean -C figures/
~ make -ik clean -C koch/
~ make -ik clean -C kochxmpl/
~ -rm -Rf *~ *.o *.exe *.dat *.tgz *.trj *.aux *.log *.idx *.scn *.dvi
~ ~
clean:
~ make -ik tidy
~ -rm -Rf \$(PROJECT) *.c *.pdf *mp *.toc *.tex *.ps
~ ~
archive:
~ make -ik tidy
~ tar --gzip --directory=../ -cf ../\$(PROJECT).tgz \$(PROJECT)
}
\bigskip
This {\tt Makefile} essentially executes two major calls. First, the \CTANGLE\
program parses the \CWEB\ source document {\tt boxcount.w} to extract a \CEE\
source file {\tt boxcount.c} which may be compiled in the usual way using any
\ANSICEE\ conformant compiler. The output source file includes {\tt \#line}
specifications so that any debugging can be done conveniently in terms of
the original \CWEB\ source file.
Second, the \CWEAVE\ program parses the same \CWEB\ source file to extract a
plain \TeX\ source file {\tt boxcount.tex} which may be compiled in the usual
way.
It takes appropriate care of typographic details like page layout and the use
of indentation, italics, boldface, and so on, and it supplies extensive
cross-index information that it gathers automatically.
After having executed {\tt make} in the same catalogue where the files
{\tt boxcount.w} and {\tt Makefile} are located, one is left with an
executable file {\tt boxcount}, being the ready-to-use compiled program,
and a PostScript file {\tt boxcount.ps}
which contains the full documentation of the program, that is to say the
document you currently are reading.
Notice that on platforms running Windows NT, Windows 2000, Windows ME, or any
other operating system by Microsoft, the executable file will instead
automatically be called {\tt boxcount.exe}.
This convention also applies to programs compiled under the \UNIX-like
environment \CYGWIN.
@*Running the program. The program is entirely controlled by the command
line options supplied when invoking the program. The syntax for executing the
program is\par
\medskip
\line{\hskip 40pt{\tt boxcount [options]}\hfill}\par
\medskip
\noindent
where {\tt options} include the following, given in their long as well as
their short forms (with prefixes `{\tt --}' and `{\tt -}', respectively):
\medskip
\optitem[{\tt --trajectoryfile}, {\tt -i} $\langle${\it trajectory filename}%
$\rangle$]{Specifies the input trajectory of the graph whose fractal
dimension is to be estimated. The input file should describe the trajectory
as two columns of $x$- and $y$-coordinates of the nodes, between which
straight lines will interpolate the trajectory. Unless the boundary of the
computational window is explicitly stated using the {\tt -w} or
{\tt --computationwindow} options, the minimum and maximum $x$- and
$y$-values of the specified trajectory will be used for the automatic
internal computation of the proper computational domain boundaries.}
\optitem[{\tt --outputfile}, {\tt -o [append{\char'174}overwrite]}
$\langle${\it output filename}$\rangle$]{Specifies the base name of the file
to which the calculated output data will be written.
If the {\tt --outputfile} or {\tt -o} option is followed by ``{\tt append}''
the estimate for the fractal dimension will be appended to the file named
$\langle${\it output filename}$\rangle$.dat, which will be created if it
does not exist. If the following word instead is ``{\tt overwrite}'' the
file will, of course, instead be overwritten.}
\optitem[{\tt -w}, {\tt --computationwindow llx} $\langle x_{\rm LL}\rangle$
{\tt lly} $\langle y_{\rm LL}\rangle$ {\tt urx} $\langle x_{\rm UR}\rangle$
{\tt ury} $\langle y_{\rm UR}\rangle$]{This option explicitly specifies the
domain over which the box-counting algorithm will be applied.
By specifying this option, any automatic calculation of computational window
will be neglected.}
\optitem[{\tt --verbose}, {\tt -v}]{Use this option to toggle verbose mode of
the program execution, controlling the amount of information written to
standard terminal output. (Default is off, that is to say quiet mode.)}
\optitem[{\tt --boxmaps}, {\tt -m}]{If this option is present, the \boxcount\
program will generate \METAPOST\ code for maps of the distribution of boxes,
so-called ``boxmaps''. In doing so, also the input trajectory is included in
the graphs. The convention for the naming of the output map files is that
they are saved as $\langle${\it output filename}$\rangle$.$N$.{\tt dat},
where $\langle${\it output filename}$\rangle$ is the base filename as
specified using the {\tt -o} or {\tt --outputfile} option, $N$ is the
automatically appended current level of resolution refinement in the
box-counting (that is to say, indicating the calculation performed for a
$[2^{N}\times2^{N}]$-grid of boxes), and where {\tt dat} is a file suffix
as automatically appended by the program. This option is useful for tracking
the performance of the program, and for verification of the box counting
algorithm.}
\optitem[{\tt --graphsize} $\langle w\rangle$ $\langle h\rangle$]{If the
{\tt -m} or {\tt --boxmaps} option is present at the command line, then
the {\tt --graphsize} option specifies the physical width $w$ and height
$h$ in millimetres of the generated graphs, overriding the default sizes
$w=80\,{\rm mm}$ and $h=80\,{\rm mm}$.}
\optitem[{\tt --minlevel}, {\tt -Nmin} $\langle N_{\rm min}\rangle$]{Specifies
the minimum level $N_{\rm min}$ of grid refinement that will be used in the
evaluation of the estimate of the fractal dimension.
With this option specified, the coarsest level used in the box-counting will
be a $[2^{N_{\rm min}}\times2^{N_{\rm min}}]$-grid of boxes.}
\optitem[{\tt --maxlevel}, {\tt -Nmax} $\langle N_{\rm max}\rangle$]{This
option is similar to the {\tt --minlevel} or {\tt -Nmin} options, with the
difference that it instead specifies the maximum level $N_{\rm max}$ of grid
refinement that will be used in the evaluation of the estimate of the fractal
dimension.
With this option specified, the finest level used in the box-counting will
be a $2^{N_{\rm max}}\times2^{N_{\rm max}}$-grid of boxes.
If this option is omitted, a default value of $N_{\rm max}=10$ will be
used as default.}
@*Application example: The Koch fractal.
The Koch curve is one of the earliest described fractal curves, appearing
in the 1904 article {\sl Sur une courbe continue sans tangente, obtenue par
une construction g\'eom\'etrique \'el\'ementaire}, written by the Swedish
mathematician Helge von Koch (1870--1924).
In this application example, we employ the ``snowflake'' variant of the Koch
fractal (merely for the sake of its beauty).
The Koch snowflake fractal is constructed recursively using the following
algorithm.
\medskip
\noindent{\bf Algorithm A} ({\it Construction of the Koch snowflake fractal}).
\aitem[{\bf A1.}]{[Create initiator.] Draw three line segments of equal length
so that they form an equilateral triangle.}
\aitem[{\bf A2.}]{[Line division.] Divide each of the line segments into three
segments of equal length.}
\aitem[{\bf A3.}]{[Replace mid segment by triangle.] Draw an equilateral
triangle that has the middle segment from step
one as its base.}
\aitem[{\bf A4.}]{[Remove base of triangle.] Remove the line segment that is
the base of the triangle from step A3.}
\aitem[{\bf A5.}]{[Recursion step.] For each of the line segments remaining,
repeat steps A2 through A4.}
\quad\endalg
\medskip
The triangle resulting after step A1 is called the {\sl initiator} of the
fractal.
After the first iteration of steps A1--A4, the result should be a shape similar
to the Star of David; this is called the {generator} of the fractal.
The Koch fractal resulting of the infinite iteration of the algorithm of
construction has an infinite length, since each time the steps above are
performed on each line segment of the figure there are four times as many line
segments, the length of each being one-third the length of the segments in the
previous stage.
Hence, the total length of the perimeter of the fractal increases by $4/3$
at each step, and for an initiator of total length $L$ the total length
of the perimeter at the $n$th step of iteration will be $(4/3)^nL$.
The fractal dimension is hence $D=\ln4/\ln3\approx1.26$, being greater
than the dimension of a line ($D=1$) but less than, for example, Peano's
space-filling curve ($D=2$).
The Koch fractal is an example of a curve which is continuous, but not
differentiable anywhere.
The area of the Koch snowflake is $8/5$ that of the initial triangle,
so an infinite perimeter encloses a finite area.
The stepwise construction of the snowflake fractal is illustrated in Fig.~1.
\bigskip
\centerline{\epsfxsize=54mm\epsfbox{koch/kochgraph-1.1}\hskip2mm
\epsfxsize=54mm\epsfbox{koch/kochgraph-2.1}\hskip2mm
\epsfxsize=54mm\epsfbox{koch/kochgraph-3.1}}
\centerline{\epsfxsize=54mm\epsfbox{koch/kochgraph-4.1}\hskip2mm
\epsfxsize=54mm\epsfbox{koch/kochgraph-5.1}\hskip2mm
\epsfxsize=54mm\epsfbox{koch/kochgraph-6.1}}
\medskip\noindent
{\bf Figure 1.} Construction of the Koch fractal of the ``snowflake'' type,
in this case as inscribed in the unitary circle. For this case, the length
of the initiator is $L=3\sqrt{3}$ while its area is
$A=L^2/(6\sqrt{3})=3\sqrt{3}/2$. For each fractal recursion depth $n$, the
trajectory consists of a total of $3\times4^{(n-1)}$ connected linear segments.
\vfill\eject
The data set for the Koch fractal is straightforward to generate by means of
recursion, as for example by using the following compact program (which in
fact was used for the generation of the data sets for the Koch fractals on
the previous page):
\bigskip
{\obeylines\obeyspaces\tt
/*-------------------------------------------------------------------------
{\char'174} File: koch.c [ANSI-C conforming source code]
{\char'174} Created: May 8, 2006, Fredrik Jonsson
{\char'174} Last modified: May 8, 2006, Fredrik Jonsson
{\char'174} Compile with: gcc -O2 -g -Wall -pedantic -ansi koch.c -o koch -lm
{\char'174} Description: The KOCH program creates data sets corresponding to
{\char'174} the Koch fractal, for the purpose of acting as test objects for
{\char'174} the BOXCOUNT program. The KOCH program is simply executed by
{\char'174} 'koch ', where is an integer describing the maximum
{\char'174} depth of recursion in the generation of the fractal data set.
{\char'174} If invoked without any arguments, =6 is used as default.
{\char'174} The generated trajectory is written to standard output.
{\char'174} Copyright (C) 2006 Fredrik Jonsson
=========================================================================*/
\#include
\#include
extern char *optarg;
\hskip\lineskip
void kochsegment(double xa,double ya,double xb,double yb,
\ int depth,int maxdepth) {\char'173}
\ double xca,yca,xcb,ycb,xcc,ycc;
\ if (depth==maxdepth) {\char'173}
\ fprintf(stdout,"\%2.8f \%2.8f{\char'134}n",xb,yb);
\ {\char'175} else {\char'173}
\ xca=xa+(xb-xa)/3.0;
\ yca=ya+(yb-ya)/3.0;
\ xcb=xb-(xb-xa)/3.0;
\ ycb=yb-(yb-ya)/3.0;
\ xcc=(xa+xb)/2.0-(yb-ya)/(2.0*sqrt(3.0));
\ ycc=(ya+yb)/2.0+(xb-xa)/(2.0*sqrt(3.0));
\ kochsegment(xa,ya,xca,yca,depth+1,maxdepth);
\ kochsegment(xca,yca,xcc,ycc,depth+1,maxdepth);
\ kochsegment(xcc,ycc,xcb,ycb,depth+1,maxdepth);
\ kochsegment(xcb,ycb,xb,yb,depth+1,maxdepth);
\ {\char'175}
{\char'175}
\hskip\lineskip
int main(int argc, char *argv[]) {\char'173}
\ int maxdepth=6;
\ if (argc>1) sscanf(argv[1],"\%d",{\char'046}maxdepth);
\ if (maxdepth>0) {\char'173}
\ fprintf(stdout,"\%2.8f \%2.8f{\char'134}n",0.0,1.0);
\ kochsegment(0.0,1.0,sqrt(3.0)/2.0,-0.5,1,maxdepth);
\ kochsegment(sqrt(3.0)/2.0,-0.5,-sqrt(3.0)/2.0,-0.5,1,maxdepth);
\ kochsegment(-sqrt(3.0)/2.0,-0.5,0.0,1.0,1,maxdepth);
\ {\char'175}
\ return(0);
{\char'175}
}\bigskip
\vfill\eject
The boxcounting dimension of the Koch snowflake fractal can now be investigated
with assistance of the \boxcount\ program. In the analysis as here presented,
this is done using the following Makefile:
\bigskip
{\obeylines\obeyspaces\tt
\#
\# Makefile designed for use with gcc, MetaPost and plain TeX.
\#
\# Copyright (C) 2002-2006, Fredrik Jonsson
\#
CC = gcc
CCOPTS = -O2 -Wall -ansi -std=iso9899:1990 -pedantic
LNOPTS = -lm
TEX = tex
DVIPS = dvips
METAPOST = mpost
\#
\# Define path and executable file for the BOXCOUNT program, used for
\# calculating estimates of the box-counting fractal dimension of a
\# trajectory in the plane.
\#
BOXCOUNTPATH = ../
BOXCOUNT = \$(BOXCOUNTPATH)/boxcount
\#
\# Define path and executable file for the KOCH program, used for generating
\# the trajectory of the Koch snowflake fractal.
\#
KOCHPATH = ../koch/
KOCH = \$(KOCHPATH)/koch
all: koch kochgen kochtab
koch:
\ make -C ../koch/
kochgen:
\ \$(KOCH) 7 > koch.trj
\ \$(BOXCOUNT) --verbose --boxmaps --graphsize 42.0 42.0 {\char'134}
\ --computationwindow llx -1.1 lly -1.1 urx 1.1 ury 1.1 {\char'134}
\ --minlevel 3 --maxlevel 8 {\char'134}
\ --trajectoryfile koch.trj --outputfile overwrite koch
\ for k in 03 04 05 06 07 08; do {\char'134}
\ \$(METAPOST) koch.\$\$k.mp ;{\char'134}
\ \$(TEX) -jobname=koch.\$\$k '{\char'134}input epsf{\char'134}nopagenumbers{\char'134}
\ {\char'134}centerline{{\char'134}epsfxsize=120mm{\char'134}%
epsfbox{koch.'\$\$k'.1}}{\char'134}bye';{\char'134}
\ \$(DVIPS) -D1200 -E koch.\$\$k -o koch.\$\$k.eps;{\char'134}
\ done
\hskip\lineskip
kochtab:
\ \$(KOCH) 7 > koch.trj
\ \$(BOXCOUNT) --verbose --minlevel 3 --maxlevel 14 {\char'134}
\ --computationwindow llx -1.1 lly -1.1 urx 1.1 ury 1.1 {\char'134}
\ --trajectoryfile koch.trj --outputfile overwrite koch
\hskip\lineskip
clean:
\ -rm -Rf koch *~ *.o *.exe *.dat *.mp *.mpx *.trj
\ -rm -Rf *.tex *.aux *.log *.toc *.idx *.scn *.dvi *.1 *.eps
}\bigskip
\vfill\eject
Having executed the {\tt Makefile} as displayed in the previous page, where a
recursion depth of $n=7$ is used for the generation of the Koch fractal, we are
left with a set of images of the consecutively refined grids in the boxcounting
algorithm, and a table containing the estimates of the boxcounting dimension
of the Koch snowflake fractal.
In Fig.~2 the resulting maps of the boxes used in the boxcounting algorithm
for refinement levels $m=3,4,\ldots,8$ are shown, and as the grid refinement
is finer and finer, the boxes finally will be barely visible in the limited
resolution of computer graphics, on screen as well as in the generated
Encapsulated PostScript code.
\bigskip
\centerline{\epsfxsize=54mm\epsfbox{kochxmpl/koch-03.1}\hskip2mm
\epsfxsize=54mm\epsfbox{kochxmpl/koch-04.1}\hskip2mm
\epsfxsize=54mm\epsfbox{kochxmpl/koch-05.1}}
\centerline{\epsfxsize=54mm\epsfbox{kochxmpl/koch-06.1}\hskip2mm
\epsfxsize=54mm\epsfbox{kochxmpl/koch-07.1}\hskip2mm
\epsfxsize=54mm\epsfbox{kochxmpl/koch-08.1}}
\medskip\noindent
{\bf Figure 2.} Consecutive steps of refinement $m=3,4,\ldots,8$ of the grid
of boxes covering the input trajectory, in this case a Koch snowflake fractal
of recursion depth $n=7$ (see Fig.~1). At each step $m$ of refinement, a
virtual grid of $2^m\times2^m$ boxes is superimposed on the trajectory
(fractal) and the number $N_m$ of boxes covering the trajectory are counted.
For the Koch snowflake fractal of recursion depth $n=7$, the trajectory
consists of a total of $3\times4^{(7-1)}=12288$ connected linear segments.
\vfill\eject
The estimated boxcounting dimension $D_m=\ln(N_m)/\ln(2^m)$, with $N_m$ as
previously denoting the number of boxes in a $2^m\times2^m$-grid required to
cover the entire curve, is displayed in Table~A.1.
The values for the estimates could be compared with the analytically obtained
value of $D=\ln4/\ln3\approx1.26$ for the fractal dimension.
However, it should be emphasized that the box counting dimension here just
is an estimate of one definition of the fractal dimension, which in many cases
do not equal to other measures, and that we in the computation of the estimate
always will arrive at a truncated result due to a limited precision and a
limited amount of memory resources and computational time.
As can be seen in the table, the initial estimates at lower resolutions are
pretty crude, but in the final estimate we nevertheless end up with the box
counting estimate of the fractal dimension as 1.29, which is reasonably close
to the analytically obtained value of $D\approx1.26$.
Of course, further refinements such as Richardson extrapolation could be
applied to increase the accuracy, but this is outside of the scope of the
\boxcount\ program, which only serves to provide the raw, basic algorithm
of boxcounting.\footnote{$\dagger$}{For discussions on
different definitions of the fractal dimension, see
the English Wikipedia section on the Minkowski--Bouligand dimension,
{\tt http://en.wikipedia.org/wiki/Minkowski-Bouligand{\char'137}dimension}.}
\bigskip
\noindent{\bf Table A.1} Boxcounting estimates of the fractal dimension of
the Koch snowflake fractal of recursion order $n=7$. In the table, $m$ is the
refinement order as indicated in the graphs in Fig.~2, $N_m$ is the number of
covering boxes counted at refinement level $m$, and $D_m=\ln(N_m)/\ln(2^m)$
is the estimate of the boxcounting dimension.
\smallskip
\noindent{\hrule width 328pt}\vskip1pt
\beginvrulealign
% \tabskip=0pt
\halign{
\hbox to 72pt{\hfil#\hfil}% first column centered
&\hbox to 128pt{\hfil#\hfil}% second column is centered
&\hbox to 128pt{\hfil#\hfil}\cr% third column is centered
\noalign{{\hrule width 328pt}\vskip 2pt}
$m$ & $N_m$ & $D_m=\ln(N_m)/\ln(2^m)$\cr
\noalign{\vskip 1pt{\hrule width 328pt}\vskip 1.4pt}
3 & 44 & 1.8198\cr
4 & 96 & 1.6462\cr
5 & 196 & 1.5229\cr
6 & 504 & 1.4962\cr
7 & 1180 & 1.4578\cr
8 & 2856 & 1.4350\cr
9 & 6844 & 1.4156\cr
10 & 15620 & 1.3931\cr
11 & 32320 & 1.3618\cr
12 & 66200 & 1.3345\cr
13 & 133600 & 1.3098\cr
14 & 268804 & 1.2883\cr
}\endvrulealign
\vskip-9.55pt
{\hskip-.4pt\hrule width 328pt}\vskip 1pt
\noindent{\hskip-.2pt\hrule width 328pt}
@*The main program. Here follows the general outline of the main program.
@c
@@;
@@;
@@;
@@;
int main(int argc, char *argv[])
{
@@;
@@;
@@;
@@;
@@;
@@;
@@;
@@;
@@;
@@;
@@;
@@;
return(SUCCESS);
}
@ Library dependencies.
The standard \ANSICEE\ libraries included in this program are:\par
\varitem[{\tt math.h}]{For access to common mathematical functions.}
\varitem[{\tt stdio.h}]{For file access and any block involving |fprintf|.}
\varitem[{\tt stdlib.h}]{For access to the |exit| function.}
\varitem[{\tt string.h}]{For string manipulation, |strcpy|, |strcmp| etc.}
\varitem[{\tt ctype.h}]{For access to the |isalnum| function.}
@=
#include
#include
#include
#include
#include
#include
@ Global definitions.
There are just a few global definitions present in the \boxcount\ program:\par
\varitem[{\tt VERSION}]{The current program revision number.}
\varitem[{\tt COPYRIGHT}]{The copyright banner.}
\varitem[{\tt SUCCESS}]{The return code for successful program termination.}
\varitem[{\tt FAILURE}]{The return code for unsuccessful program termination.}
\varitem[{\tt NCHMAX}]{The maximum number of characters allowed in strings for
storing file names, including path.}
\varitem[{\tt APPEND}]{Code for the flag |output_mode|, used to determine if
output should append to an existing file.}
\varitem[{\tt OVERWRITE}]{Code for the flag |output_mode|, used to determine
if output should overwrite an existing file.}
@=
#define VERSION "1.6"
#define COPYRIGHT "Copyright (C) 2006-2011, Fredrik Jonsson"
#define SUCCESS (0)
#define FAILURE (1)
#define NCHMAX (256)
#define APPEND 1
#define OVERWRITE 2
@ Declaration of global variables. The only global variables allowed in
my programs are |optarg|, which is a pointer to the the string of characters
that specified the call from the command line, and |progname|, which simply
is a pointer to the string containing the name of the program, as it was
invoked from the command line.
@=
extern char *optarg;
char *progname;
@ Declaration of local variables of the |main| program.
In \CWEB\ one has the option of adding variables along the program, for example
by locally adding temporary variables related to a given sub-block of code.
However, the philosophy in the \boxcount\ program is to keep all variables of
the |main| section collected together, so as to simplify tasks as, for example,
tracking down a given variable type definition.
The local variables of the program are as follows:
\medskip
\varitem[|verbose|]{Boolean flag which, if being nonzero, tells the program
to display information at terminal output during program execution.}
\varitem[|user_set_compwin|]{Boolean flag which indicates whether the user has
explicitly stated a window of computation or not.}
\varitem[|output_mode|]{Variable which states whether the estimate of
fractal dimension should be appended to an existing file which is
created if it does not exist (for $|output_mode|=|APPEND|$), or if
the data should overwrite already existing data in the file (for
$|output_mode|=|OVERWRITE|$).}
\varitem[|make_boxmaps|]{If nonzero, then graphs showing the distribution
of covering boxes will be generated.}
\varitem[|*num_boxes|]{Pointer to array holding the number of boxes at
various depths of division.}
\varitem[|initime|]{The time at which the \boxcount\ program was initialized.}
\varitem[|now|]{Dummy variable for extraction of current time from the system
clock.}
\varitem[|mm|]{The number of points $M$ in the input trajectory. This is the
length of the vectors |x_traj| and |y_traj| as described below.}
\varitem[|nn|]{Gives the number of boxes in the $x$- or $y$-directions as
$2^N$.}
\varitem[|nnmax|]{The maximum refinement depth $N_{\rm max}$ of $N$.}
\varitem[|nnmin|]{The minimum refinement depth $N_{\rm min}$ of $N$.}
\varitem[|global_llx|, |global_lly|]{Lower-left coordinates of global window
of computation.}
\varitem[|global_urx|, |global_ury|]{Upper-right coordinates of global window
of computation.}
\varitem[|*x_traj|,|*y_traj|]{Vectors containing the $x$- and $y$-values of
the coordinates along the input trajectory.}
\varitem[|*x|,|*y|]{Variables for keeping the refinement and number of boxes.
(This needs to be changed as they are easily confused with $x$ and $y$
coordinates of the trajectory.)}
\varitem[|*trajectory_file|]{Input file pointer, for reading the trajectory
whose fractal dimension is to be estimated.}
\varitem[|*frac_estimate_file|]{Output file pointer, for saving the estimated
fractal dimension of the input trajectory.}
\varitem[|*boxmap_file|]{Output file pointer, for saving maps of the spatial
locations of the boxes covering the trajectory.}
\varitem[|boxgraph_width|]{The physical width in millimetres of the generated
\MP\ boxmap graphs.}
\varitem[|boxgraph_height|]{The physical height in millimetres of the generated
\MP\ boxmap graphs.}
\varitem[|trajectory_filename|]{String for keeping the name of the file
containing the input trajectory.}
\varitem[|output_filename|]{String for keeping the base name of the set of
output files.}
\varitem[|frac_estimate_filename|]{String keeping the name of the file
to which the estimate of the fractal dimension will be saved.}
\varitem[|boxmap_filename|]{String for keeping the name of the file to
which \METAPOST\ code for maps of the spatial distribution of boxes will be
written.}
\varitem[|no_arg|]{Variable for extracting the number of input arguments
present on the command line as the program is executed.}
\varitem[|i|]{Dummy variable used in loops when reading the input trajectory.}
\varitem[|*fracdimen_estimates|]{Vector keeping the values characterizing
estimated fractal dimension for various orders of $N$.}
\varitem[|ave,adev,sdev|]{The average, average deviation, and standard
deviation of the estimated fractal dimension for various orders of
refinement $N$, as stored in |fracdimen_estimates|.}
\varitem[|var,skew,curt|]{The variance, skewness, and kurtosis of the
estimated fractal dimension for various orders of refinement $N$, as
stored in |fracdimen_estimates|.}
\medskip
\noindent
Generally in this program, the maximum number of characters a file name
string can contain is |NCHMAX|, as defined in the definitions section
of the program.
@=
short verbose,user_set_compwin,output_mode,make_boxmaps;
long int *num_boxes;
time_t initime;
time_t now;
long mm,nn,nnmin,nnmax;
double global_llx,global_lly,global_urx,global_ury;
double *x_traj,*y_traj,*x,*y;
FILE *trajectory_file,*frac_estimate_file,*boxmap_file;
char trajectory_filename[NCHMAX],output_filename[NCHMAX],frac_estimate_filename[NCHMAX];
char boxmap_filename[NCHMAX];
double boxgraph_width,boxgraph_height;
int no_arg;
long i;
double *fracdimen_estimates,ave,adev,sdev,var,skew,curt;
@ Initialization of variables.
@=
verbose=0; /* Verbose mode is off by default */
user_set_compwin=0; /* Before the command-line is parsed, nothing is known of these settings */
output_mode=OVERWRITE; /* default mode is to overwrite existing files */
make_boxmaps=0; /* Default is to not create graphs of the box distributions */
nnmin=0; /* Default value for $N_{\rm min}$ */
nnmax=10; /* Default value for $N_{\rm max}$ */
strcpy(output_filename,"out"); /* Default output file basename. */
strcpy(trajectory_filename,""); /* To indicate that no filename has been set yet. */
boxgraph_width=80.0; /* Default graph width in millimetres */
boxgraph_height=56.0; /* Default graph height in millimetres */
trajectory_file=NULL;
frac_estimate_file=NULL;
boxmap_file=NULL;
initime=time(NULL); /* Time of initialization of the program. */
@ Parsing command line options. All input parameters are passed to the
program through command line options and arguments to the program.
The syntax of command line options is listed
whenever the program is invoked without any options, or whenever any of the
{\tt --help} or {\tt -h} options are specified at startup.
@=
{
progname=strip_away_path(argv[0]);
fprintf(stdout,"This is %s v.%s. %s\n",progname,VERSION,COPYRIGHT);
no_arg=argc;
while(--argc) {
if(!strcmp(argv[no_arg-argc],"-o") ||
!strcmp(argv[no_arg-argc],"--outputfile")) {
--argc;
if(!strcmp(argv[no_arg-argc],"append") ||
!strcmp(argv[no_arg-argc],"a")) {
output_mode=APPEND;
} else if(!strcmp(argv[no_arg-argc],"overwrite") ||
!strcmp(argv[no_arg-argc],"o")) {
output_mode=OVERWRITE;
} else {
fprintf(stderr,"%s: Error in '-o' or '--outputfile' option!",
progname);
exit(FAILURE);
}
--argc;
strcpy(output_filename,argv[no_arg-argc]);
} else if(!strcmp(argv[no_arg-argc],"-w") ||
!strcmp(argv[no_arg-argc],"--computationwindow")) {
user_set_compwin=1;
--argc;
if(!strcmp(argv[no_arg-argc],"llx")) {
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&global_llx)) {
fprintf(stderr,"%s: Error in parsing lower-left x-value.\n",
progname);
exit(FAILURE);
}
} else {
fprintf(stderr,"%s: Error in computation window option\n",
progname);
fprintf(stderr,"%s: Expecting 'llx'\n",progname);
exit(FAILURE);
}
--argc;
if(!strcmp(argv[no_arg-argc],"lly")) {
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&global_lly)) {
fprintf(stderr,"%s: Error in parsing lower-left y-value.\n",
progname);
exit(FAILURE);
}
} else {
fprintf(stderr,"%s: Error in computation window option\n",
progname);
fprintf(stderr,"%s: Expecting 'lly'\n",progname);
exit(FAILURE);
}
--argc;
if(!strcmp(argv[no_arg-argc],"urx")) {
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&global_urx)) {
fprintf(stderr,"%s: Error in parsing lower-left x-value.\n",
progname);
exit(FAILURE);
}
} else {
fprintf(stderr,"%s: Error in computation window option\n",
progname);
fprintf(stderr,"%s: Expecting 'urx'\n",progname);
exit(FAILURE);
}
--argc;
if(!strcmp(argv[no_arg-argc],"ury")) {
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&global_ury)) {
fprintf(stderr,"%s: Error in parsing lower-left y-value.\n",
progname);
exit(FAILURE);
}
} else {
fprintf(stderr,"%s: Error in computation window option\n",
progname);
fprintf(stderr,"%s: Expecting 'ury'\n",progname);
exit(FAILURE);
}
} else if(!strcmp(argv[no_arg-argc],"-i") ||
!strcmp(argv[no_arg-argc],"--trajectoryfile")) {
--argc;
strcpy(trajectory_filename,argv[no_arg-argc]);
} else if(!strcmp(argv[no_arg-argc],"-v") ||
!strcmp(argv[no_arg-argc],"--verbose")) {
verbose=(verbose?0:1);
} else if(!strcmp(argv[no_arg-argc],"-h") ||
!strcmp(argv[no_arg-argc],"--help")) {
showsomehelp();
exit(SUCCESS);
} else if(!strcmp(argv[no_arg-argc],"-m") ||
!strcmp(argv[no_arg-argc],"--boxmaps")) {
make_boxmaps=1;
} else if(!strcmp(argv[no_arg-argc],"--graphsize")) {
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&boxgraph_width)) {
fprintf(stderr,"%s: Error in width of '--graphsize' option.\n",
progname);
exit(FAILURE);
}
--argc;
if(!sscanf(argv[no_arg-argc],"%lf",&boxgraph_height)) {
fprintf(stderr,"%s: Error in height of '--graphsize' option.\n",
progname);
exit(FAILURE);
}
} else if(!strcmp(argv[no_arg-argc],"--minlevel") ||
!strcmp(argv[no_arg-argc],"-Nmin")) {
--argc;
if(!sscanf(argv[no_arg-argc],"%ld",&nnmin)) {
fprintf(stderr,"%s: Error in '--minlevel' or '-Nmin' option.\n",
progname);
exit(FAILURE);
}
} else if(!strcmp(argv[no_arg-argc],"--maxlevel") ||
!strcmp(argv[no_arg-argc],"-Nmax")) {
--argc;
if(!sscanf(argv[no_arg-argc],"%ld",&nnmax)) {
fprintf(stderr,"%s: Error in '--maxlevel' or '-Nmax' option.\n",
progname);
exit(FAILURE);
}
} else {
fprintf(stderr,"%s: Specified option '%s' invalid!\n",
progname,argv[no_arg-argc]);
showsomehelp();
exit(FAILURE);
}
}
}
@ Display starting time of program execution.
@=
{
fprintf(stdout,"%s: Program execution started %s",progname,ctime(&initime));
}
@ Load input trajectory from file. This is the section where the trajectory
to be analyzed is loaded into the memory, and where the number~$M$ of input
coordinates present in the input data is determined by a single call to
the subroutine |num_coordinate_pairs()|.
@={
if (!strcmp(trajectory_filename,"")) {
fprintf(stderr,"%s: No input trajectory specified!\n",progname);
fprintf(stderr,"%s: Please use the '--trajectoryfile' option.\n",
progname);
fprintf(stderr,"%s: Use '--help' option to display a help message.\n",
progname);
exit(FAILURE);
}
if ((trajectory_file=fopen(trajectory_filename,"r"))==NULL) {
fprintf(stderr,"%s: Could not open %s for loading trajectory!\n",
progname,trajectory_filename);
exit(FAILURE);
}
mm=num_coordinate_pairs(trajectory_file);
x_traj=dvector(1,mm);
y_traj=dvector(1,mm);
for (i=1;i<=mm;i++) {
fscanf(trajectory_file,"%lf",&x_traj[i]); /* scan $x$-coordinate */
fscanf(trajectory_file,"%lf",&y_traj[i]); /* scan $y$-coordinate */
}
fclose(trajectory_file);
if (verbose) {
fprintf(stdout,"%s: Loaded %ld coordinate pairs from file '%s'.\n",
progname,mm,trajectory_filename);
}
}
@ Opening files for output by the program.
@={
if (!strcmp(output_filename,"")) {
fprintf(stderr,"%s: No output base name specified!\n",progname);
fprintf(stderr,"%s: Please use the '--outputfile' option.\n",
progname);
exit(FAILURE);
}
sprintf(frac_estimate_filename,"%s.dat",output_filename);
if (output_mode==APPEND) {
if ((frac_estimate_file=fopen(frac_estimate_filename,"a"))==NULL) {
fprintf(stderr,"%s: Could not open %s for output!\n",
progname,frac_estimate_filename);
exit(FAILURE);
}
fseek(frac_estimate_file,0L,SEEK_END);
/* set file pointer to the end of the file */
} else if (output_mode==OVERWRITE) {
if ((frac_estimate_file=fopen(frac_estimate_filename,"w"))==NULL) {
fprintf(stderr,"%s: Could not open %s for loading trajectory!\n",
progname,frac_estimate_filename);
exit(FAILURE);
}
fseek(frac_estimate_file,0L,SEEK_SET);
/* set file pointer to the beginning of the file */
} else {
fprintf(stderr,"%s: Error: Output mode (output_mode) undefined!",progname);
exit(FAILURE);
}
}
@ Extract global window of computation. This block will only be executed if
there was no explicit computational window stated at startup of the program
(that is to say, with the {\tt -w} or {\tt --computationwindow} option).
In order to determine the minimum area covering the entire input trajectory,
this section scans sequentially through the trajectory and finds the minimum
and maximum of the $x$- and $y$-coordinates.
These values then form the lower-left and upper-right corner coordinates
$(|global_llx|,|global_lly|)$ and $(|global_urx|,|global_ury|)$.
@=
{
if (!user_set_compwin) {
global_llx=x_traj[1];
global_lly=y_traj[1];
global_urx=global_llx;
global_ury=global_lly;
for (i=1;i<=mm;i++) {
if (x_traj[i]>global_urx) global_urx=x_traj[i];
if (y_traj[i]>global_ury) global_ury=y_traj[i];
if (x_traj[i]={
num_boxes=livector(1,nnmax);
nn=1;
for (i=1;i<=nnmin-1;i++) nn=2*nn; /* This leaves |nn| as $2^{(N_{\rm min}-1)}$ */
for (i=nnmin;i<=nnmax;i++) { /* For all levels in refinement of grid resolution */
nn=2*nn;
if (make_boxmaps) {
/* do we wish to generate \METAPOST\ graphs of the box distribution? */
sprintf(boxmap_filename,"%s-%02ld.mp",output_filename,i);
num_boxes[i]=get_num_covering_boxes_with_boxmaps(x_traj,y_traj,mm,i,
global_llx,global_lly,global_urx,global_ury,boxmap_filename,
boxgraph_width,boxgraph_height,trajectory_filename);
} else {
/* if not, just use the regular number-crunching routine */
num_boxes[i]=get_num_covering_boxes(x_traj,y_traj,mm,i,
global_llx,global_lly,global_urx,global_ury);
}
if (verbose) {
fprintf(stdout,"%s: N=%ld (%ldx%ld-grid of boxes): ",progname,i,nn,nn);
fprintf(stdout,"Trajectory covered by %ld boxes\n",num_boxes[i]);
}
}
}
@ Compute the logarithmic estimate of the fractal dimension. Having completed
the task of actually calculating the number of boxes necessary to cover the
trajectory, for varying box dimensions of width
$(|global_urx|-|global_llx|)/2^n$ and height $(|global_ury|-|global_lly|)/2^n$,
for $n=1,2,\ldots,N$, the only remaining arithmetics is to actually calculate
the estimate for the fractal dimension.
This is done by performing the fit of the linear function $y=ax+b$ to the data
set obtained with $2\ln(n)$ as abscissa and $\ln(|num_boxes|(n))$ as ordinata.
Also in this block, we analyze whether the option |make_boxmaps| is set as true
(1) or not (0), determining whether a graph of the number of boxes as function
of division depth should be written to file as well.
In the fitting of the linear function $y=ax+b$ to the data set $(x_1,y_1)$,
$(x_2,y_2)$,$\ldots$,$(x_N,y_N)$, the parameters $a$ and $b$ can be explicitly
obtained from the summation formulas
$$
a={{1}\over{D}}\Big(N\sum^N_{k=1}x_ky_k-\sum^N_{k=1}x_k\sum^N_{k=1}y_k\Big),
\qquad
b={{1}\over{N}}\Big(\sum^N_{k=1}y_k-a\sum^N_{k=1}x_k\Big),
$$
where
$$
D\equiv N\sum^N_{k=1}x^2_k-\Big(\sum^N_{k=1}x_k\Big)^2.
$$
@={
x=dvector(nnmin,nnmax);
y=dvector(nnmin,nnmax);
fracdimen_estimates=dvector(nnmin,nnmax);
nn=1;
for (i=1;i<=nnmax;i++) {
nn=2*nn;
if (nnmin<=i) x[i]=log((double)nn);
}
for (i=nnmin;i<=nnmax;i++) y[i]=log(num_boxes[i]);
for (i=nnmin;i<=nnmax;i++) fracdimen_estimates[i]=y[i]/x[i];
if (verbose) {
for (i=nnmin;i<=nnmax;i++) {
fprintf(stdout,"%s: N=%ld, fractal dimension estimate=%f\n",progname,
i,fracdimen_estimates[i]);
}
}
moment(fracdimen_estimates,nnmin,nnmax,&ave,&adev,&sdev,&var,&skew,&curt);
if (verbose) {
fprintf(stdout,"%s: Estimate of fractal dimension: %f\n",progname,ave);
fprintf(stdout,"%s: Average deviation: %f\n",progname,adev);
fprintf(stdout,"%s: Standard deviation: %f\n",progname,sdev);
fprintf(stdout,"%s: Skewness: %f\n",progname,skew);
}
free_livector(num_boxes,1,nnmax); /* release the memory occupied by |num_boxes| */
free_dvector(fracdimen_estimates,nnmin,nnmax);
free_dvector(x,nnmin,nnmax); /* release the memory occupied by |x| */
free_dvector(y,nnmin,nnmax); /* release the memory occupied by |y| */
}
@ Save the fractal dimension to file.
@={
if (output_mode==APPEND) {
fseek(frac_estimate_file,0L,SEEK_END);
} else if (output_mode==OVERWRITE) {
fseek(frac_estimate_file,0L,SEEK_SET);
}
fprintf(frac_estimate_file,"%f %f %f\n",ave,sdev,skew);
}
@ Close any open output files.
@={
fclose(frac_estimate_file);
}
@ Display elapsed execution time.
@={
now=time(NULL);
if (verbose)
fprintf(stdout,"%s: Total execution time: %d s\n",progname,
((int)difftime(now,initime)));
fprintf(stdout,"%s: Program execution closed %s",progname,ctime(&now));
}
@*Subroutines. In this section, all subroutines as used by the main program
are listed.
@=
@@;
@@;
@@;
@@;
@@;
@@;
@@;
@@;
@ Routine for computation of average, average deviation and standard deviation.
This routine is adopted from {\sl Numerical Recipes in C}.
@=
void moment(double data[],int nnmin,int nnmax,double *ave,double *adev,
double *sdev,double *var,double *skew,double *curt) {
int j;
double ep=0.0,s,p;
if (nnmax-nnmin <= 1) {
fprintf(stderr,"%s: Error in routine moment()! ",progname);
fprintf(stderr,"(nnmax-nnmin>1 is a requirement)\n");
exit(FAILURE);
}
s=0.0;
for (j=nnmin;j<=nnmax;j++) s += data[j];
*ave=s/(nnmax-nnmin+1);
*adev=(*var)=(*skew)=(*curt)=0.0;
for (j=nnmin;j<=nnmax;j++) {
*adev += fabs(s=data[j]-(*ave));
ep += s;
ep += s;
*var += (p=s*s);
*skew += (p *= s);
*curt += (p *= s);
}
*adev /= (nnmax-nnmin+1);
*var=(*var-ep*ep/(nnmax-nnmin+1))/(nnmax-nnmin);
*sdev=sqrt(*var);
if (*var) {
*skew /= ((nnmax-nnmin+1)*(*var)*(*sdev));
*curt=(*curt)/((nnmax-nnmin+1)*(*var)*(*var))-3.0;
} else {
fprintf(stderr,"%s: Error in routine moment()! ",progname);
fprintf(stderr,"No skew/kurtosis for zero variance\n");
exit(FAILURE);
}
}
@ Routine for obtaining the number of coordinate pairs in a file. This
routine is called prior to loading the trajectory, in order to get the
size needed for allocating the memory for the trajectory vector. As the
number of coordinates~$M$ has been established, two vectors |x_traj[1..M]|
and |y_traj[1..M]|, containing the coordinates $(x_m,y_m)$ of the trajectory.
@=
long int num_coordinate_pairs(FILE *trajectory_file) {
double tmp;
int tmpch;
long int mm=0;
fseek(trajectory_file,0L,SEEK_SET); /* rewind file to beginning */
while((tmpch=getc(trajectory_file))!=EOF) {
ungetc(tmpch,trajectory_file);
fscanf(trajectory_file,"%lf",&tmp); /* Read away the $x$ coordinate */
fscanf(trajectory_file,"%lf",&tmp); /* Read away the $y$ coordinate */
mm++;
tmpch=getc(trajectory_file); /* Read away blanks and linefeeds */
while ((tmpch!=EOF)&&(!isdigit(tmpch))) tmpch=getc(trajectory_file);
if (tmpch!=EOF) ungetc(tmpch,trajectory_file);
}
fseek(trajectory_file,0L,SEEK_SET); /* rewind file to beginning */
return(mm);
}
@ Routines for removing preceding path of filenames.
In this block all routines related to removing preceding path strings go.
Not really fancy programming, and no contribution to any increase of numerical
efficiency or precision; just for the sake of keeping a tidy terminal output
of the program. The |strip_away_path()| routine is typically called when
initializing the program name string |progname| from the command line string
|argv[0]|, and is typically located in the blocks related to parsing of the
command line options.
@=
@@;
@@;
@ Checking for a valid path character.
The |pathcharacter| routine takes one character |ch| as argument, and returns
1 (``true'') if the character is valid character of a path string, otherwise 0
(``false'') is returned.
@=
short pathcharacter(int ch) {
return(isalnum(ch)||(ch=='.')||(ch=='/')||(ch=='\\')||(ch=='_')||
(ch=='-')||(ch=='+'));
}
@ Routine for stripping away path string of a file name.
The |strip_away_path()| routine takes a character string |filename| as
argument, and returns a pointer to the same string but without any
preceding path segments. This routine is, for example, useful for
removing paths from program names as parsed from the command line.
@=
char *strip_away_path(char filename[]) {
int j,k=0;
while (pathcharacter(filename[k])) k++;
j=(--k); /* this is the uppermost index of the full path+file string */
while (isalnum((int)(filename[j]))) j--;
j++; /* this is the lowermost index of the stripped file name */
return(&filename[j]);
}
@ Subroutines for memory allocation. Here follows the routines for memory
allocation and deallocation of double precision real and complex valued
vectors, as used for storing the optical field distribution in the grating,
the refractive index distribution of the grating, etc.
@=
@@;
@@;
@@;
@@;
@@;
@@;
@ The |dvector| routine allocates a real-valued vector of double precision,
with vector index ranging from |nl| to |nh|.
@=
double *dvector(long nl, long nh) {
double *v;
v=(double *)malloc((size_t) ((nh-nl+2)*sizeof(double)));
if (!v) {
fprintf(stderr,"Error: Allocation failure in dvector()\n");
exit(FAILURE);
}
return v-nl+1;
}
@ The |free_dvector| routine release the memory occupied by the
real-valued vector |v[nl..nh]|.
@=
void free_dvector(double *v, long nl, long nh) {
free((char*) (v+nl-1));
}
@ The |ivector| routine allocates a real-valued vector of integer precision,
with vector index ranging from |nl| to |nh|.
@=
int *ivector(long nl, long nh) {
int *v;
v=(int *)malloc((size_t) ((nh-nl+2)*sizeof(int)));
if (!v) {
fprintf(stderr,"Error: Allocation failure in ivector()\n");
exit(FAILURE);
}
return v-nl+1;
}
@ The |free_ivector| routine release the memory occupied by the
real-valued vector |v[nl..nh]|.
@=
void free_ivector(int *v, long nl, long nh) {
free((char*) (v+nl-1));
}
@ The |livector| routine allocates a real-valued vector of long integer
precision, with vector index ranging from |nl| to |nh|.
@=
long int *livector(long nl, long nh) {
long int *v;
v=(long int *)malloc((size_t) ((nh-nl+2)*sizeof(long int)));
if (!v) {
fprintf(stderr,"Error: Allocation failure in livector()\n");
exit(FAILURE);
}
return v-nl+1;
}
@ The |free_livector| routine release the memory occupied by the
real-valued vector |v[nl..nh]|.
@=
void free_livector(long int *v, long nl, long nh) {
free((char*) (v+nl-1));
}
@ The |simatrix| routine allocates an array of short integer precision,
with array row index ranging from |nrl| to |nrh| and column index ranging
from |ncl| to |nch|.
@=
short int **simatrix(long nrl, long nrh, long ncl, long nch)
{
long i, nrow=nrh-nrl+1,ncol=nch-ncl+1;
short int **m;
m=(short int **) malloc((size_t)((nrow+1)*sizeof(short int*)));
if (!m) {
fprintf(stderr,"%s: Allocation failure 1 in simatrix()\n",progname);
exit(FAILURE);
}
m += 1;
m -= nrl;
m[nrl]=(short int *) malloc((size_t)((nrow*ncol+1)*sizeof(short int)));
if (!m[nrl]) {
fprintf(stderr,"%s: Allocation failure 2 in simatrix()\n",progname);
exit(FAILURE);
}
m[nrl] += 1;
m[nrl] -= ncl;
for(i=nrl+1;i<=nrh;i++) m[i]=m[i-1]+ncol;
return m;
}
@ The |free_simatrix| routine releases the memory occupied by the
short integer matrix |v[nl..nh]|, as allocated by |simatrix()|.
@=
void free_simatrix(short int **m,long nrl,long nrh,long ncl,long nch) {
free((char*) (m[nrl]+ncl-1));
free((char*) (m+nrl-1));
}
@ Routine for displaying help message to standard terminal output.
@=
void showsomehelp(void) {
fprintf(stderr,"Usage: %s [options]\n",progname);
fprintf(stderr,"Options:\n");
fprintf(stderr," -h, --help\n");
fprintf(stderr," Display this help message and exit clean.\n");
fprintf(stderr," -v, --verbose\n");
fprintf(stderr," Toggle verbose mode on/off.\n");
fprintf(stderr," -o, --outputfile \n");
fprintf(stderr," Specifies the base name of the output files where the program\n");
fprintf(stderr," is to save the calculated data. If the --outputfile or -o\n");
fprintf(stderr," option is followed by 'append' the estimate for the fractal\n");
fprintf(stderr," dimension will be appended to the file named .dat, which\n");
fprintf(stderr," will be created if it does not exist. If the following word\n");
fprintf(stderr," instead is `overwrite' the file will instead be overwritten.\n");
fprintf(stderr," -i, --trajectoryfile\n");
fprintf(stderr," Specifies the input trajectory of the graph whose fractal\n");
fprintf(stderr," dimension is to be estimated.\n");
fprintf(stderr," -w, --computationwindow llx lly urx ury \n");
fprintf(stderr," This option explicitly specifies the domain over which the\n");
fprintf(stderr," box-counting algorithm will be applied, in terms of the\n");
fprintf(stderr," lower-left and upper-right corners (llx,lly) and (urx,ury),\n");
fprintf(stderr," respectively. By specifying this option, any automatic\n");
fprintf(stderr," calculation of computational window will be neglected.\n");
fprintf(stderr," -m, --boxmaps\n");
fprintf(stderr," If this option is present, the program will generate\n");
fprintf(stderr," MetaPost code for maps of the distribution of boxes.\n");
fprintf(stderr," In doing so, also the input trajectory is included in\n");
fprintf(stderr," the graphs. The convention for the naming of the output\n");
fprintf(stderr," map files is that they are saved as ..dat,\n");
fprintf(stderr," where is the base filename as specified using the\n");
fprintf(stderr," -o or --outputfile option, is the automatically appended\n");
fprintf(stderr," current level of resolution refinement in the box-counting,\n");
fprintf(stderr," and where '.dat' is the file suffix as automatically appended\n");
fprintf(stderr," by the program.\n");
fprintf(stderr," --graphsize \n");
fprintf(stderr," If the -m or --boxmaps option is present at the command line,\n");
fprintf(stderr," then the --graphsize option will override the default graph\n");
fprintf(stderr," size of the generated boxmaps. (Default graph size is 80 mm\n");
fprintf(stderr," width and 56 mm height.)\n");
fprintf(stderr," -Nmin, --minlevel \n");
fprintf(stderr," Specifies the minimum level Nmin of grid refinement that \n");
fprintf(stderr," will be used in the evaluation of the estimate of the fractal\n");
fprintf(stderr," dimension. With this option specified, the coarsest level\n");
fprintf(stderr," used in the box-counting will be a [(2^Nmin)x(2^Nmin)]-grid\n");
fprintf(stderr," of boxes.\n");
fprintf(stderr," -Nmax, --maxlevel \n");
fprintf(stderr," This option is similar to the --minlevel or -Nmin options,\n");
fprintf(stderr," with the difference that it instead specifies the maximum\n");
fprintf(stderr," level Nmax of grid refinement that will be used in the\n");
fprintf(stderr," evaluation of the estimate of the fractal dimension.\n");
}
@ The |lines_intersect(p1x,p1y,q1x,q1y,p2x,p2y,q2x,q2y)| routine takes the
start and end points of two line segments, the first between $(|p1x|,|p1y|)$
and $(|q1x|,|q1y|)$ and the second between $(|p2x|,|p2y|)$ and $(|q2x|,|q2y|)$,
and returns $1$ (`true') if they are found to intersect, and $0$ (`false')
otherwise.
For a brief sumnmary of the algorithm behind this routine, consider two line
segments in the plane, the first one between the points ${\bf p}_1$ and
${\bf q}_1$ and the second one between ${\bf p}_2$ and ${\bf q}_1$.
In general, these segments can be written in parametric forms as
${\bf r}_1={\bf p}_1+t_1({\bf q}_1-{\bf p}_1)$ and
${\bf r}_2={\bf p}_2+t_2({\bf q}_2-{\bf p}_2)$
for $0\le t_1\le1$ and $0\le t_2\le1$.
These line segments intersect each other if they for these intervals for the
parameters $t_1$ and $t_2$ share a common point, or equivalently if the
solutions to
$$
{\bf p}_1+t_1({\bf q}_1-{\bf p}_1)={\bf p}_2+t_2({\bf q}_2-{\bf p}_2)
\qquad\Leftrightarrow\qquad
({\bf q}_1-{\bf p}_1)t_1+({\bf p}_2-{\bf q}_2)t_2={\bf p}_2-{\bf p}_1
$$
both simultaneously satisfy $0\le t_1\le1$ and $0\le t_2\le1$.
This vectorial equation can equivalently be written in component form as
$$
\eqalign{
(q_{1x}-p_{1x})t_1+(p_{2x}-q_{2x})t_2&=p_{2x}-p_{1x},\cr
(q_{1y}-p_{1y})t_1+(p_{2y}-q_{2y})t_2&=p_{2y}-p_{1y},\cr
}
$$
which after some straightforward algebra gives the explicit solutions for the
parameters $t_1$ and $t_2$ as
$$
t_1={{ed-bf}\over{ad-bc}},\qquad t_2={{af-ec}\over{ad-bc}},
$$
where
$$
\eqalign{
a\equiv(q_{1x}-p_{1x}),\qquad
b\equiv(p_{2x}-q_{2x}),\qquad
c\equiv(q_{1y}-p_{1y}),\cr
d\equiv(p_{2y}-q_{2y}),\qquad
e\equiv(p_{2x}-p_{1x}),\qquad
f\equiv(p_{2x}-p_{1x}).\cr
}
$$
Hence, the two line segments intersect if and only if
$$
0\le{{ed-bf}\over{ad-bc}}\le1,\qquad
{\rm and}\qquad 0\le{{af-ec}\over{ad-bc}}\le1.
$$
By observing that their denominators are equal, the evaluation of these quotes
involves in total 6 floating-point multiplications, 2 divisions and 3
subtractions.
Notice that whenever $ad-bd=0$, the two lines are parallell and will never
intersect, regardless of the values of $t_1$ and $t_2$.
@=
short int lines_intersect(double p1x,double p1y,double q1x,double q1y,
double p2x,double p2y,double q2x,double q2y) {
double a,b,c,d,e,f,det,tmp1,tmp2;
short int intersect;
a=q1x-p1x;
b=p2x-q2x;
c=q1y-p1y;
d=p2y-q2y;
e=p2x-p1x;
f=p2y-p1y;
det=a*d-b*c;
tmp1=e*d-b*f;
tmp2=a*f-e*c;
intersect=0;
if (det>0) {
if (((0.0<=tmp1)&&(tmp1<=det))&&((0.0<=tmp2)&&(tmp2<=det))) intersect=1;
} else if (det<0) {
if (((det<=tmp1)&&(tmp1<=0.0))&&((det<=tmp2)&&(tmp2<=0.0))) intersect=1;
}
return(intersect);
}
@ Routine for determining whether a line and a box intersect or not.
The |box_intersection()| routine simply divides the input box, being
characterized by its lower-left and upper-right corners $(|llx|,|lly|)$
and $(|urx|,|ury|)$, into the four line segments corresponding to its four
edges, followed by calling the routine |lines_intersect()| to determine if
any of these edges intersect the line segment. If an intersection is found,
the |box_intersection()| routine returns 1 (true), otherwise 0 (false).
\medskip\noindent
Input variables:
\varitem[|px|, |py|]{The coordinates of the first end point
${\bf p}=(p_x,p_y)$ of the line segment.}
\varitem[|qx|, |qy|]{The coordinates of the second end point
${\bf q}=(q_x,q_y)$ of the line segment.}
\varitem[|llx|, |lly|]{Coordinates of the lower-left corner of the box.}
\varitem[|urx|, |ury|]{Coordinates of the upper-right corner of the box.}
\medskip\noindent
Output variables:
\varitem[]{On exit, the routine returns 1 if an intersection is found,
otherwise 0 is returned, in either case the value are returned as
integers of |short| precision.}
\medskip
@=
short int box_intersection(double px,double py,double qx,double qy,
double llx,double lly,double urx,double ury) {
if (lines_intersect(px,py,qx,qy,llx,lly,urx,lly))
return(1); /* intersection with bottom edge */
if (lines_intersect(px,py,qx,qy,urx,lly,urx,ury))
return(1); /* intersection with right edge */
if (lines_intersect(px,py,qx,qy,urx,ury,llx,ury))
return(1); /* intersection with top edge */
if (lines_intersect(px,py,qx,qy,llx,ury,llx,lly))
return(1); /* intersection with left edge */
return(0); /* this happens only if no edge is intersecting the line segment */
}
@ Routines for calculation of number of boxes covering the trajectory. There
are two almost identical routines for the calculation of the number of boxes
covering the input trajectory at a given level of subdivision of the box sizes.
The first routine, |get_num_covering_boxes()| simply performs this task
without caring of documenting the box distributions as graphs, or ``box maps'',
while the second one, |get_num_covering_boxes_with_boxmaps()| also includes
the generation of these maps, with output in terms of \METAPOST\ code.
@=
@@;
@@;
@ Routine for calculation of number of boxes covering the trajectory, also
including the generation of documenting graphs of the box distributions.
The |get_num_covering_boxes_with_boxmaps()| routine takes a trajectory as
described by a discrete set of coordinates as input, and for a given grid
refinement $N$ returns the number of boxes needed to entirely cover the
trajectory. The grid refinement factor $N$ indicates that the domain of
computation is divided into a $[2^N\times2^N]$-grid of boxes.
The computational domain in which the box counting is to be performed is
explicitly stated by the coordinates of its lower-left and upper-right
corners, $(global_llx,global_lly)$ and $(global_urx,global_ury)$, respectively.
Parts of the trajectory which are outside of this domain are not included in
the box-counting.
If the entire input trajectory is to be used in the box counting, simply use
$(|global_llx|,|global_lly|)=({\rm min}[|x_traj|],{\rm min}[|y_traj|])$ and
$(|global_urx|,|global_ury|)=({\rm max}[|x_traj|],{\rm max}[|y_traj|])$ for the
specification of the computational domain.
\medskip\noindent
Input variables:
\varitem[|mm|]{The total number of coordinates forming the input trajectory,
or equivalently the length of the vectors |*x_traj| and |*y_traj|.}
\varitem[|*x_traj|, |*y_traj|]{Vectors of length |mm| containing the $x$- and
$y$-coordinates of the input trajectory.}
\varitem[|resolution|]{The grid refinement factor $N$.}
\varitem[|global_llx|, |global_lly|]{Coordinates of the lower-left corner
of the computational domain in which the box-counting is to be performed.}
\varitem[|global_urx|, |global_ury|]{Coordinates of the upper-right corner
of the computational domain in which the box-counting is to be performed.}
\varitem[|trajectory_filename|]{String containing the name of the file
containing the input trajectory.}
\varitem[|boxgraph_filename|]{String which, if non-empty, states the file name
to which the map of the spatial distribution of the covering boxes is to be
written, as \METAPOST\ code.}
\medskip\noindent
Output variables:
\varitem[]{On exit, the routine returns the number of covering boxes
as an integer of |long unsigned| precision.}
\medskip\noindent
Internal variables:
\varitem[|px|,|py|,|qx|,|qy|]{Keeps track of the $x$- and $y$-coordinates of
the start and end points of line segments, between ${\bf p}=(p_x,p_y)$ and
${\bf q}=(q_x,q_y)$.}
\medskip
@=
long unsigned int get_num_covering_boxes_with_boxmaps(
double *x_traj,double *y_traj,
long int mm,int resolution,double global_llx,double global_lly,
double global_urx,double global_ury,char boxgraph_filename[],
double boxgraph_width,double boxgraph_height,
char trajectory_filename[]) {
short int **box;
long unsigned int i,j,m,nn,istart,jstart,istop,jstop,iincr,jincr,num_boxes;
double *x_box,*y_box; /* Keeps track of the lower-left coordinates of
the boxes. */
double px,py,qx,qy;
FILE *boxgraph_file;
@@;
@@;
for (m=1;m<=mm-1;m++) {
@@;
@@;
}
@@;
@@;
@@;
@@;
@@;
return(num_boxes);
}
@ Set up the grid of $2^N\times2^N$ boxes covering the entire global window
of computation.
In this block, the resolution of the grid of boxes is defined. Notice that in
many cases, only a certain subset of boxes will actally intersect the input
trajectory. However, we do not \'a priori know this number of boxes, and in
order to safeguard against the possible danger of running out of allocated
memory, with the need for a dynamic memory allocation depending on the state
of computation, a matrix of short integer precision is allocated covering the
entire computational window.
In order to keep track of the coordinates of the boxes, two vectors
$|x_box|[1\ldots 2^N]$ and $|y_box|[1\ldots 2^N]$ are allocated to contain
the $x$- and $y$-coordinates of the lower-left corners of the boxes, with the
last elements $|x_box|[2^N]$ and $|y_box|[2^N]$ containing the upper-right
corner coordinates of the upper-right-most box of the global window.
@=
{
nn=1;
for (i=1;i<=resolution;i++) nn=2*nn;
box=simatrix(1,nn,1,nn);
for (i=1;i<=nn;i++) for (j=1;j<=nn;j++) box[i][j]=0;
x_box=dvector(1,nn+1);
y_box=dvector(1,nn+1);
for (m=1;m<=nn+1;m++) {
x_box[m]=global_llx+((double)(m-1))*(global_urx-global_llx)/((double)(nn));
y_box[m]=global_lly+((double)(m-1))*(global_ury-global_lly)/((double)(nn));
}
}
@ Find indices $(i_a,j_a)$ of the box covering the first coordinate of the
trajectory.
@=
{
istart=0;
jstart=0;
for (m=1;m<=nn;m++) {
if ((x_box[m]<=x_traj[1])&&(x_traj[1]<=x_box[m+1])) istart=m;
if ((y_box[m]<=y_traj[1])&&(y_traj[1]<=y_box[m+1])) jstart=m;
}
if ((istart==0)||(jstart==0)) {
fprintf(stderr,"%s: Error! Cannot find box indices of 1st coordinate!\n",
progname);
fprintf(stderr,"%s: Please check data for input trajetory.\n",progname);
exit(FAILURE);
}
}
@ Find indices $(i_b,j_b)$ of the box covering the end point of the $m$th
trajectory segment.
@={
px=x_traj[m];
py=y_traj[m];
qx=x_traj[m+1];
qy=y_traj[m+1];
istop=istart;
jstop=jstart;
if (pxqx) istop--;
}
if (pyqy) jstop--;
}
if (0==1) {
fprintf(stdout,"%s: Endpoint box indices: (i=%ld,j=%ld)\n",
progname,istop,jstop);
}
}
@ Scan the $m$th trajectory segment for intersecting boxes.
As the indices of the boxes covering the start and end points of the $m$th
trajectory segment have been previously determined, one may now use this
information in order to reduce the search for intersecting boxes to the
domain as covered by box indices $i_{\rm start}\le i\le i_{\rm stop}$ and
$j_{\rm start}\le j\le j_{\rm stop}$.
This way, the computational time needed is greatly reduced as compared to if
the entire window would be scanned for every segment.
@=
{
for (i=istart;i!=(istop+iincr);i+=iincr) {
for (j=jstart;j!=(jstop+jincr);j+=jincr) {
if (box_intersection(px,py,qx,qy,
x_box[i],y_box[j],x_box[i+1],y_box[j+1])) {
box[i][j]=1;
}
}
}
istart=istop;
jstart=jstop;
}
@ Open file for output of box distribution graph. In this block the preamble
of the output \METAPOST\ code is written to file. This preamble contains a
macro {\tt box(i,j)} which simply is a parameter-specific \METAPOST\ subroutine
which is used in order to reduce the final size of the code to be compiled
into a graph over the distribution of covering boxes.
@={
if (!strcmp(boxgraph_filename,"")) { /* is |boxgraph_filename| an empty string? */
boxgraph_file=NULL;
} else {
if ((boxgraph_file=fopen(boxgraph_filename,"w"))==NULL) {
fprintf(stderr,"%s: Could not open %s for box graphs!\n",
progname,boxgraph_filename);
exit(FAILURE);
}
fseek(boxgraph_file,0L,SEEK_SET);
fprintf(boxgraph_file,"input graph;\n");
fprintf(boxgraph_file,"def box(expr i,j)=\n");
fprintf(boxgraph_file,"begingroup\n");
fprintf(boxgraph_file,"llx:=%2.8f+(i-1)*%2.8f;\n",
global_llx,(global_urx-global_llx)/(nn));
fprintf(boxgraph_file,"lly:=%2.8f+(j-1)*%2.8f;\n",
global_lly,(global_ury-global_lly)/(nn));
fprintf(boxgraph_file,"urx:=%2.8f+(i)*%2.8f;\n",
global_llx,(global_urx-global_llx)/(nn));
fprintf(boxgraph_file,"ury:=%2.8f+(j)*%2.8f;\n",
global_lly,(global_ury-global_lly)/(nn));
fprintf(boxgraph_file,"gdraw (llx,lly)--(urx,lly);\n");
fprintf(boxgraph_file,"gdraw (urx,lly)--(urx,ury);\n");
fprintf(boxgraph_file,"gdraw (urx,ury)--(llx,ury);\n");
fprintf(boxgraph_file,"gdraw (llx,ury)--(llx,lly);\n");
fprintf(boxgraph_file,"endgroup\n");
fprintf(boxgraph_file,"enddef;\n");
fprintf(boxgraph_file,"beginfig(1);\n");
fprintf(boxgraph_file,"w:=%2.4fmm; h:=%2.4fmm;\n",
boxgraph_width,boxgraph_height);
fprintf(boxgraph_file,"draw begingraph(w,h);\n");
fprintf(boxgraph_file,"pickup pencircle scaled .3pt;\n");
fprintf(boxgraph_file,"setrange(%2.6f,%2.6f,%2.6f,%2.6f);\n",
global_llx,global_lly,global_urx,global_ury);
}
}
@ Write the input trajectory to the box distribution graph. Here there are two
possible choices for how the input trajectory is to be included in the box
graph.
First, we may use \MP\ to automatically scan and draw the trajectory
for us, simply by using a statment like {\tt gdraw input.dat}, assuming that
the file {\tt input.dat} contains properly formatted columnwise data.
However, this choice would have two major drawbacks, namely that the generated
code would be dependent on that the original input file always is in place,
hence not allowing the \MP\ code to be exported as a standalone code as such,
and also that this would put a limitation on the number of nodes allowed in the
input trajectory, as the {\tt graph.mp} macro package of \MP\ only accepts
roughly 4000 points before it cuts the mapping.
The other choice is to include the input trajectory directly into the generated
code, preferrably by chopping the trajectory into pieces small enough to easily
be mapped by \MP\ without reaching a critical limit of exhaustion.
This choice of course significantly increases the file size of the generated
code, but this is a price we will have to accept in order to get stand-alone
output.
In the \boxcount\ program, the second alternative was naturtally chosen, simply
because the author is a fan of self-sustaining code which can be exported for
later use.
In this block, the status of the pointer |boxgraph_file| is used to determine
whether to write the trajectory to file or not. If |boxgraph_file| equals to
|NULL| (NULL), then the \boxcount\ program will not attempt to write the input
trajectory to file.
@={
if (boxgraph_file!=NULL) {
fprintf(boxgraph_file,"pickup pencircle scaled .5pt;\n");
i=0;
for (m=1;m<=mm;m++) {
if (i==0) {
if (m==1) {
fprintf(boxgraph_file,"gdraw (%2.4f,%2.4f)",
x_traj[m],y_traj[m]);
} else if (2=
{
num_boxes=0;
for (i=1;i<=nn;i++) {
for (j=1;j<=nn;j++) {
if(box[i][j]==1) {
num_boxes++;
if (boxgraph_file!=NULL) {
fprintf(boxgraph_file,"box(%ld,%ld);\n",i,j);
}
}
}
}
free_simatrix(box,1,nn,1,nn);
}
@ Write closing blocks the box distribution graph. Here follows the
specification of tick marking (set as automatic for the sake of simplicity)
and axis labels, just before the closing statements of the \METAPOST\ code
for the graphs of box distributions.
@={
if (boxgraph_file!=NULL) {
fprintf(boxgraph_file,"autogrid(itick bot,itick lft);\n");
fprintf(boxgraph_file,"glabel.bot(btex $x$ etex,OUT);\n");
fprintf(boxgraph_file,"glabel.lft(btex $y$ etex,OUT);\n");
fprintf(boxgraph_file,"endgraph;\n");
fprintf(boxgraph_file,"endfig;\n");
fprintf(boxgraph_file,"end\n");
}
}
@ Close any open file for output of box distribution graph.
@={
if (boxgraph_file!=NULL) {
fprintf(stdout,
"%s: MetaPost output box distribution graph written to %s.\n",
progname,boxgraph_filename);
fclose(boxgraph_file);
}
}
@ Routine for calculation of number of boxes covering the trajectory. This
routine provides a simplified interface to the general boxcountng routine
|get_num_covering_boxes_with_boxmaps|, with the only difference that no
graph of the box distribution over the trajectory is generated.
The |get_num_covering_boxes()| routine takes a trajectory as
described by a discrete set of coordinates as input, and for a given grid
refinement $N$ returns the number of boxes needed to entirely cover the
trajectory. The grid refinement factor $N$ indicates that the domain of
computation is divided into a $[2^N\times2^N]$-grid of boxes.
The computational domain in which the box counting is to be performed is
explicitly stated by the coordinates of its lower-left and upper-right
corners, $(global_llx,global_lly)$ and $(global_urx,global_ury)$, respectively.
Parts of the trajectory which are outside of this domain are not included in
the box-counting.
If the entire input trajectory is to be used in the box counting, simply use
$(|global_llx|,|global_lly|)=({\rm min}[|x_traj|],{\rm min}[|y_traj|])$ and
$(|global_urx|,|global_ury|)=({\rm max}[|x_traj|],{\rm max}[|y_traj|])$ for the
specification of the computational domain.
\medskip\noindent
Input variables:
\varitem[|mm|]{The total number of coordinates forming the input trajectory,
or equivalently the length of the vectors |*x_traj| and |*y_traj|.}
\varitem[|*x_traj|, |*y_traj|]{Vectors of length |mm| containing the $x$- and
$y$-coordinates of the input trajectory.}
\varitem[|resolution|]{The grid refinement factor $N$.}
\varitem[|global_llx|, |global_lly|]{Coordinates of the lower-left corner
of the computational domain in which the box-counting is to be performed.}
\varitem[|global_urx|, |global_ury|]{Coordinates of the upper-right corner
of the computational domain in which the box-counting is to be performed.}
\medskip\noindent
Output variables:
\varitem[]{On exit, the routine returns the number of covering boxes
as an integer of |long unsigned| precision.}
\medskip
@=
long unsigned int get_num_covering_boxes(double *x_traj,double *y_traj,long int mm,
int i,double global_llx,double global_lly,
double global_urx,double global_ury) {
return(get_num_covering_boxes_with_boxmaps(x_traj,y_traj,mm,i,
global_llx,global_lly,global_urx,global_ury,"",0.0,0.0,""));
}
@*Index.