[go: up one dir, main page]

Menu

[r999]: / libmona / doc / userguide.tex  Maximize  Restore  History

Download this file

1169 lines (967 with data), 44.8 kB

   1
   2
   3
   4
   5
   6
   7
   8
   9
  10
  11
  12
  13
  14
  15
  16
  17
  18
  19
  20
  21
  22
  23
  24
  25
  26
  27
  28
  29
  30
  31
  32
  33
  34
  35
  36
  37
  38
  39
  40
  41
  42
  43
  44
  45
  46
  47
  48
  49
  50
  51
  52
  53
  54
  55
  56
  57
  58
  59
  60
  61
  62
  63
  64
  65
  66
  67
  68
  69
  70
  71
  72
  73
  74
  75
  76
  77
  78
  79
  80
  81
  82
  83
  84
  85
  86
  87
  88
  89
  90
  91
  92
  93
  94
  95
  96
  97
  98
  99
 100
 101
 102
 103
 104
 105
 106
 107
 108
 109
 110
 111
 112
 113
 114
 115
 116
 117
 118
 119
 120
 121
 122
 123
 124
 125
 126
 127
 128
 129
 130
 131
 132
 133
 134
 135
 136
 137
 138
 139
 140
 141
 142
 143
 144
 145
 146
 147
 148
 149
 150
 151
 152
 153
 154
 155
 156
 157
 158
 159
 160
 161
 162
 163
 164
 165
 166
 167
 168
 169
 170
 171
 172
 173
 174
 175
 176
 177
 178
 179
 180
 181
 182
 183
 184
 185
 186
 187
 188
 189
 190
 191
 192
 193
 194
 195
 196
 197
 198
 199
 200
 201
 202
 203
 204
 205
 206
 207
 208
 209
 210
 211
 212
 213
 214
 215
 216
 217
 218
 219
 220
 221
 222
 223
 224
 225
 226
 227
 228
 229
 230
 231
 232
 233
 234
 235
 236
 237
 238
 239
 240
 241
 242
 243
 244
 245
 246
 247
 248
 249
 250
 251
 252
 253
 254
 255
 256
 257
 258
 259
 260
 261
 262
 263
 264
 265
 266
 267
 268
 269
 270
 271
 272
 273
 274
 275
 276
 277
 278
 279
 280
 281
 282
 283
 284
 285
 286
 287
 288
 289
 290
 291
 292
 293
 294
 295
 296
 297
 298
 299
 300
 301
 302
 303
 304
 305
 306
 307
 308
 309
 310
 311
 312
 313
 314
 315
 316
 317
 318
 319
 320
 321
 322
 323
 324
 325
 326
 327
 328
 329
 330
 331
 332
 333
 334
 335
 336
 337
 338
 339
 340
 341
 342
 343
 344
 345
 346
 347
 348
 349
 350
 351
 352
 353
 354
 355
 356
 357
 358
 359
 360
 361
 362
 363
 364
 365
 366
 367
 368
 369
 370
 371
 372
 373
 374
 375
 376
 377
 378
 379
 380
 381
 382
 383
 384
 385
 386
 387
 388
 389
 390
 391
 392
 393
 394
 395
 396
 397
 398
 399
 400
 401
 402
 403
 404
 405
 406
 407
 408
 409
 410
 411
 412
 413
 414
 415
 416
 417
 418
 419
 420
 421
 422
 423
 424
 425
 426
 427
 428
 429
 430
 431
 432
 433
 434
 435
 436
 437
 438
 439
 440
 441
 442
 443
 444
 445
 446
 447
 448
 449
 450
 451
 452
 453
 454
 455
 456
 457
 458
 459
 460
 461
 462
 463
 464
 465
 466
 467
 468
 469
 470
 471
 472
 473
 474
 475
 476
 477
 478
 479
 480
 481
 482
 483
 484
 485
 486
 487
 488
 489
 490
 491
 492
 493
 494
 495
 496
 497
 498
 499
 500
 501
 502
 503
 504
 505
 506
 507
 508
 509
 510
 511
 512
 513
 514
 515
 516
 517
 518
 519
 520
 521
 522
 523
 524
 525
 526
 527
 528
 529
 530
 531
 532
 533
 534
 535
 536
 537
 538
 539
 540
 541
 542
 543
 544
 545
 546
 547
 548
 549
 550
 551
 552
 553
 554
 555
 556
 557
 558
 559
 560
 561
 562
 563
 564
 565
 566
 567
 568
 569
 570
 571
 572
 573
 574
 575
 576
 577
 578
 579
 580
 581
 582
 583
 584
 585
 586
 587
 588
 589
 590
 591
 592
 593
 594
 595
 596
 597
 598
 599
 600
 601
 602
 603
 604
 605
 606
 607
 608
 609
 610
 611
 612
 613
 614
 615
 616
 617
 618
 619
 620
 621
 622
 623
 624
 625
 626
 627
 628
 629
 630
 631
 632
 633
 634
 635
 636
 637
 638
 639
 640
 641
 642
 643
 644
 645
 646
 647
 648
 649
 650
 651
 652
 653
 654
 655
 656
 657
 658
 659
 660
 661
 662
 663
 664
 665
 666
 667
 668
 669
 670
 671
 672
 673
 674
 675
 676
 677
 678
 679
 680
 681
 682
 683
 684
 685
 686
 687
 688
 689
 690
 691
 692
 693
 694
 695
 696
 697
 698
 699
 700
 701
 702
 703
 704
 705
 706
 707
 708
 709
 710
 711
 712
 713
 714
 715
 716
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741
 742
 743
 744
 745
 746
 747
 748
 749
 750
 751
 752
 753
 754
 755
 756
 757
 758
 759
 760
 761
 762
 763
 764
 765
 766
 767
 768
 769
 770
 771
 772
 773
 774
 775
 776
 777
 778
 779
 780
 781
 782
 783
 784
 785
 786
 787
 788
 789
 790
 791
 792
 793
 794
 795
 796
 797
 798
 799
 800
 801
 802
 803
 804
 805
 806
 807
 808
 809
 810
 811
 812
 813
 814
 815
 816
 817
 818
 819
 820
 821
 822
 823
 824
 825
 826
 827
 828
 829
 830
 831
 832
 833
 834
 835
 836
 837
 838
 839
 840
 841
 842
 843
 844
 845
 846
 847
 848
 849
 850
 851
 852
 853
 854
 855
 856
 857
 858
 859
 860
 861
 862
 863
 864
 865
 866
 867
 868
 869
 870
 871
 872
 873
 874
 875
 876
 877
 878
 879
 880
 881
 882
 883
 884
 885
 886
 887
 888
 889
 890
 891
 892
 893
 894
 895
 896
 897
 898
 899
 900
 901
 902
 903
 904
 905
 906
 907
 908
 909
 910
 911
 912
 913
 914
 915
 916
 917
 918
 919
 920
 921
 922
 923
 924
 925
 926
 927
 928
 929
 930
 931
 932
 933
 934
 935
 936
 937
 938
 939
 940
 941
 942
 943
 944
 945
 946
 947
 948
 949
 950
 951
 952
 953
 954
 955
 956
 957
 958
 959
 960
 961
 962
 963
 964
 965
 966
 967
 968
 969
 970
 971
 972
 973
 974
 975
 976
 977
 978
 979
 980
 981
 982
 983
 984
 985
 986
 987
 988
 989
 990
 991
 992
 993
 994
 995
 996
 997
 998
 999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
1126
1127
1128
1129
1130
1131
1132
1133
1134
1135
1136
1137
1138
1139
1140
1141
1142
1143
1144
1145
1146
1147
1148
1149
1150
1151
1152
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164
1165
1166
1167
1168
\documentclass[english, 10pt, a4paper,headsepline]{scrbook}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
\usepackage[left=2cm,right=1cm,top=2cm,bottom=2cm,twoside]{geometry}
\usepackage{array}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage[numbers]{natbib}
\usepackage{listings}
\usepackage{color}
\usepackage{graphics}
\usepackage{nonfloat}
\usepackage{fancyheadings}
\include{version}
%\makeatletter
\definecolor{listinggray}{gray}{0.9}
\lstset{backgroundcolor=\color{listinggray}}
\usepackage{babel}
\makeatother
\begin{document}
\vfill{}
\title{Mona User and Programming Guide \\Software Version: \libmonaversion}
\vfill{}
\date{12. June 2008}
\author{Gert Wollny}
\maketitle
This is the Mona User and Programming Guide.
In the first part this document describes, how mona2d can be installed, and how the tools provided by this package can be used
to segment micro CT tooth data.
The second part is dedicated how to write software based on the infrastructure provided by mona2d.
\tableofcontents{}
\chapter*{Preface}
\section*{License}
Copyright (c) 2006 Gert Wollny, MPI for Evolutionary Anthropology
Permission is granted to copy, distribute and/or modify this document under the terms of the GNU Free Documentation License, Version 1.1
or any later version published by the Free Software Foundation; with no Invariant Sections,
with no Front-Cover Texts and with no Back-Cover Texts.
A copy of the license is available at http://www.gnu.org/copyleft/fdl.html
\section*{Changes}
\begin{center}
\begin{tabular}{|c|c|}
\hline
Date & Description of changes\\
\hline
\hline
07/29/2006 & Beta\\
\hline
04/03/2006 & Easter Egg Release\\
\hline
03/20/2006 & Initial Release\\
\hline
\end{tabular}
\end{center}
\pagestyle{fancyplain}
\chapter{Introduction}
This is the Mona User and Programming Guide. This document describes, how mona2d can be installed,
and how the tools provided by this package can be used to segment micro CT tooth data.
This document is maintained by Gert Wollny <wollny@eva.mpg.de>.
Additions, modifications or corrections may be mailed there for inclusion in the next release.
Mona2d provides the means to do various image processing tasks:
\begin{itemize}
\item 2D image filtering
\item pseudo 3D image processing based on 2D image stacks
\item Histogram analysis
\item histogram based classification
\item read and write different image formats
\end{itemize}
\section{Installation}
In order to use Mona, it has currently to be installed from source code.
To do so, your software environment has to meet the following pre-requisites:
\begin{enumerate}
\item You will need a POSIX-compatible environment - any flavor of UNIX, LINUX , BSD or MacOSX based systems should work out of the box.
For Microsoft Windows Computers you will need to install a compatibility layer like CygWin (http://www.cygwin.com/) or MingW (http://www.mingw.org/).
Be warned, however, since this has not yet been tested.
\item You need a ANSI-compatible C++ compiler - GNU g++ (>=3.3) (http://gcc.gnu.org) is known to work.
\item The BOOST (http://www.boost.org) library needs to be installed.
\item You will need the libmona library. It is provided by the MPI for Human Cognitive and Brain Sciences (http://www.cbs.mpg.de)
(Download-Page to be announced) - this software requires popt, a command line parser library.
\item Out of the box, only Windows bitmap image (BMP) are supported. To add additional image types, make sure that
libtiff (http://www.remotesensing.org/libtiff/), OpenEXR (http://www.openexr.com), and/or libpng (http://www.libpng.org) are installed.
\item To support frequency domain filtering, libfftw (http://www.fftw.org) needs also to be installed.
\end{enumerate}
If all of the above pre-requisites are met, the installation of mona2d boils down to:
\lstset{language=bash}
\begin{lstlisting}
tar -zxvf mona2d-0.1.0.tgz
cd mona2d-0.1.0
./configure --prefix=<where-to-install>
make
make install
\end{lstlisting}
In order to test some of the components, you may run {}``make check''.
However, currently only a subset of the functionality is captured
by this test.
%%\part{User Guide}
\chapter{Using the software }
Currently, all the tools provided are command line only.
The big advantage of this is, that you can put everything into scripts, run the program on large data sets, and check the results later.
Note, the notation of the examples is UNIX-centric, i.e. the directory delimiter is notated as ``/''.
\section{Filtering 2D image }
\subsection{A single image at a time}
2D Image filters take one gray scale or binary image at a time, apply a set of filters that has been defined on the command line,
and output the filtered image, i.e. all the commands run like:
\begin{lstlisting}
mona-2dimagefilter -i <input image> -o <output image> \
filter1 filter2 ...
\end{lstlisting}
The output image will be of the same type as the input image (this behavior can be changed by the -t option).
A filter is described by its name followed by a colon and a coma separated list of optional parameters, e.g.:
\begin{lstlisting}
filter:parm1=1.0,parm2=string
\end{lstlisting}
Given we want to binarize a gray-scale image ``input.png'' by setting all pixels with intensities
in the range {[}a,b{]} to 1 and all other pixels to 0, and we want to write the result to ``output.png'' then we run the filter
like this:
\begin{lstlisting}
mona-2dimagefilter -i input.png -o output.png \
binarize:min=a,max=b
\end{lstlisting}
Assume now, we want to output the image in the TIFF format, then we may run:
\begin{lstlisting}
mona-2dimagefilter -i input.png -t tif -o output.tif \
binarize:min=a,max=b
\end{lstlisting}
Note: Currently, the suffix of the output file should match the type provided by the {}``-t'' option,
neither an automatic deduction of the output type from the output file name is implemented nor an
automatic adding of the appropriate suffix to the output file name.
Assume now, we want to run two filters in a row, e.g. first a median filter of width 3 and then a binarize like above.
This goes like
\begin{lstlisting}
mona-2dimagefilter -i input.png -t tif -o output.tif \
median:w=3 binarize:min=a,max=b
\end{lstlisting}
If you want to know, which filters are available, and you want to know a short description of their parameters, you can call
\begin{lstlisting}
mona-2dimagefilter --help-plugins
\end{lstlisting}
or look it up in Section \ref{sec:plugins}.
\subsection{Filtering in the frequency domain}
The only specific thing about filtering in the frequency domain is, that it is a subset of 2D filtering where one has
to provide an additional filtering kernel like this:
\begin{lstlisting}
mona-2dimagefilter -i input.png -o output.png \
fft:kernel=[lowpass:f=100,m=gauss]
\end{lstlisting}
Here we filter input.png by using a lowpass filter with the cut-off frequency 100 and a Gaussian filter shape.
Note, currently only one level of brackets is supported (and needed).
So far lowpass, highpass bandpass and bandreject kernels are available.
\subsection{Morphological filters}
Morphological operations like erode, dilate, open and close are supported as well.
The structuring element is given similar to the kernel in the frequency domain filters:
\begin{lstlisting}
mona-2dimagefilter -i input.png -o output.png \
erode:shape=[circle:r=5]
\end{lstlisting}
Here, an erosion with a circular shape of radius 5 is applied to the input image.
Note, gray scale images can also be subject to morphological filtering.
\subsection{Filtering multiple images independently }
Now, with the tooth data, we are usually dealing with a large amount of images numbered consecutive.
To run a 2D filter on all of these, we could write a shell script that would call the above
command for all the input images.
However, loading the filter program and all its filters takes time, therefore, Mona provides a filter command that
will run the very same filter on all input files that match a certain numbering pattern and where the numbers are consecutive.
Given, e.g., you have a stack of images \texttt{input0012.png} through \texttt{input0912.png}, you can also run the binarize filter like this:
\begin{lstlisting}
mona-2dimagefilterstack -i input0000.png -o output -t tif \
binarize:min=a,max=b
\end{lstlisting}
The output will be stored in files named \texttt{output0012.tif} through \texttt{output0912.tif}.
If you don't give the type parameter, then the image file type of the input file will be used and the apropriate suffix will be added.
There are a few things to note about this:
\begin{enumerate}
\item The given input image name doesn't have to exist, the program will search for the range of consecutive numbered files on its own,
in our example, the filtering will start at \texttt{input0012.png} and end at \texttt{input0912.png}.
\item Since image \texttt{input0913.png} is missing additional images like \texttt{input0914.png} will \emph{not} be filtered.
Also note, if you have two ranges with the same pattern, like \texttt{input0012.png} to \texttt{input0912.png},
and \texttt{input0914.png} to \texttt{input1014.png}, then the current version of the software will \textit{only} filter the first range.
This behavior is in place to be compatible with \texttt{mona-2dstackfilter}, which we will discuss later.
\end{enumerate}
\section{Filtering stacks of images in a 3D manner}
Some of the filters are also provided for processing images in a 3D manner, that is, neighboring slices of an image are taken into account
in the filtering process.
Because we target for large data sets, only filters are provided for this kind of processing that require only
a limited support and can be run in {}``one sweep''.
This rules out frequency domain filters and anisotropic diffusion.
Running a stack based 3D filter is similar to running a 2D filter on the image stack:
\begin{lstlisting}
mona-2dstackfilter -i input0000.png -o output -t png \
median:w=3
\end{lstlisting}
All the limitations discussed above concerning the image ranges are still valid. \footnote{%
Especially, it becomes clear, why multiple ranges are not (yet) supported.
The gap would need to be handles somehow, and this can be handled a lot easier by applying an appropriate naming
scheme to the different image ranges.}
Just like above \texttt{mona-2dimagefilter -{}-help-plugins} provides you with a list of supported filters. See also Section \ref{sec:2dstackfilters}.
\section{Running an Analysis across an image stack}
Some of the tools provided by Mona are tailored for the analysis of the image data in stacks of images.
Amongst them are Histogram analysis and ``counting pixels''. These programs are run similar to above filters:
\begin{lstlisting}
mona-multihisto -i input0000.png -o histo.dat
\end{lstlisting}
will evaluate the histogram of the input images and write it to histo.dat.
Additionally, it will write a threshold value to the standard output that provides some hint of the intensity range that is only background.
\chapter{Running a Segmentation}
Now that we know about the basic filtering, we might take it to the next level, and run a full segmentation.
We assume, that the input images come as 16 bit gray scale valued TIFF images \texttt{recon0000.tif}.
Coincidently, that is one possible output of the SkyScan reconstruction software.
The first step is to remove most of the background by setting it to Zero.
It is more of a convenience step for the further processing that helps to reduce the file size considerably,
and also may help to speed up some of the later processing:
\begin{lstlisting}
thresh=`mona-multihisto -i recon0000.tif -o histo_start.dat`
mona-2dimagefilterstack -i recon0000.tif -o thresh -t png \
thresh:min=$thresh
\end{lstlisting}
%$ This comment just closes the math code above for the xemacs syntac highlighting
%%
Now, the image background should be zero and, since we store the output images in the PNG format,
their size should be considerably smaller then the size of the original TIF files.
This is a good thing, because the time needed to (de-)compress these images nowadays is considerably
smaller then the disk-IO time, especially, if the files reside somewhere on the net.
In the nest step we smooth the input images.
With some experimenting you may find the right filters and parameters.
So far it seems that a combination of a median filter followed by an extended Kuwahara filter provide the best results.
The filter size is, of course, depended on the input image size. The values given here give good results for input images of size 2000x2000:
\begin{lstlisting}
mona-2dstackfilter -i thresh0000.png -o filtered \
median:w=7 extkuwa:w=7
\end{lstlisting}
Depending on the number of images and their size and your computer, you may now go home and return tomorrow or
in a few days to have a look at the results.
Now, that we have smooth images, it's time to have another look at the histogram:
\begin{lstlisting}
mona-multihisto -i filtered0000.png -o histo_filtered.dat
\end{lstlisting}
This histogram is used to get a classification of the input image into the different tissues.
\begin{lstlisting}
mona-cmeans <histo_filtered.dat -k 256 -c 0,10000,30000 \
-V message >cmeans.txt
\end{lstlisting}
Be sure to use the \texttt{-V message} option here, and have a look at the output. This step is rather unstable,
and it does not always result in a classification into three classes.
If the classification fails, one can increase the -k parameter, have a look at the histogram
to see where the peaks really are and set a better class center start
vector with the -c option.
When we got three classes, we can now proceed to the next segmentation step. First we segment the enamel:
\begin{lstlisting}
mkdir enamel
mona-2dstackfilter -i filtered0000.png -o enamel \
regiongrow:map=cmeans.txt,class=2,seed=0.99,low=0.65
\end{lstlisting}
Then we segment the background:
\begin{lstlisting}
mkdir background
mona-2dstackfilter -i filtered0000.png -o background \
regiongrow:map=cmeans.txt,class=0,seed=0.99,low=0.65
\end{lstlisting}
Finally, we obtain that what is supposed to be dentine by taking all
the not yet labeled pixels:
\begin{lstlisting}
cd enamel
files=`ls *.png`
cd ..
mkdir dentin
for f in $files; do
mona-maskdentine -1 enamel/$f.png -2 background/$f -o dentine/$f
done
\end{lstlisting}
Now we have binary masks for all the three tissue classes, enamel, dentine and the background.
However, because usually matrix is attached to the tooth at scan time, and matrix most of the times
comes with similar intensities like dentine, we might try to get rid of some of the matrix parts that are
segmented as dentine, if they are not attached to the main dentine core.
To do so, we label connected components in the dentine mask, select the largest connected component only and redo the dentine mask:
\begin{lstlisting}
mona-2dstackfilter -i dentine/0000.png -o label \
label:map=joins.txt
mona-2dimagefilterstack -i label0000.png -o relabela \
labelmap:map=joins.txt
mona-multihisto -i relabel0000.png -o label.dat
mona-labelsort <label.dat >sort.dat
mona-2dimagefilterstack -i relabela0000.png -o relabelb \
labelmap:map=sort.txt
mkdir dentine2
mona-2dimagefilterstack -i relabelb0000.png -o dentine2 \
binarize:min=1,max=1
\end{lstlisting}
The first step labels connected components.
Since some labels may join during the sweep and these labels can not be corrected on the fly, we save the joining labels
in joins.txt and use this information in a second sweep to relabel joint areas.
Then we get the histogram of the labelled images and sort that histogram.
After a second relabelling step, the connected components will be sorted according to their size,
with background usually being at label number 0.
In the last step we extract label 1, as the second larges component.
If the tooth is split into multiple parts, one may also select more then just one label, by issuing, e.g,
\begin{lstlisting}
mona-2dimagefilterstack -i relabelb0000.png -o dentine2 \
binarize:min=1,max=3
\end{lstlisting}
\noindent
in order to select labels 1, 2, and 3.
\include{programs}
\include{pluginlist}
%%\part{Programming Guide}
\chapter{Adding new filters and image types}
\label{ch:addplugins}
\section{The basic implementation of image IO and filtering}
\label{sec:addioplug}
Mona makes heavy use of the plugin structure provided by the \emph{libmona}
library. The core programs mona-2dimagefilter, mona-2dimagefilterstack,
and mona-2dstackfilter only provide the functionality to collect parameters,
file names and plugin names, the actual loading, filtering and storing
of the image is delegated to the plugins. To give you an impression
of how this works, we should take a look at the main routine of mona-2dimagefilter.
At first plugin handlers for the 2D image filters and for the image
IO are needed:
\lstset{language=c++}
\begin{lstlisting}
C2DImageFilterHandler filter_plugins;
C2DImageIOPluginHandler imageio;
\end{lstlisting}
Then, the usual command line parameter definition is inserted, and
the command line is parsed to obtain the input and output file name, the output file type and the filter list:
\begin{lstlisting}
using namespace popt;
string inf;
string outf;
string outt;
vector<string> filter_chain;
COptions opt;
opt.push_back(option( inf, "in-file", 'i', "input image(s) ", NULL));
opt.push_back(option( outf, "out-base", 'o', "output file base", NULL));
opt.push_back(option( outt, imageio.get_set(), "type", 't', "out file type", NULL));
opt.push_back(option( help_plugins, "help-plugins", 0,
"give some help about the filter plugins"));
parse_options(argc, argv, options, filter_chain);
\end{lstlisting}
After checking the validity of the command line parameters obtained,
the input image is loaded by using the imageio plugin handler:
\begin{lstlisting}
auto_ptr<C2DImageList> in_image_list(imageio.load(in_filename));
\end{lstlisting}
And then filter chain is created from the filter parameters stored
in \texttt{filter\_chain}:
\begin{lstlisting}
for (std::vector<string>::const_iterator i = filter_chain.begin();
i != filter_chain.end(); ++i) {
C2DImageFilterFactory::ProductPtr f =
C2DImageFilterFactory::produce(i->c_str(), filter_plugins);
filters.push_back(f);
}
\end{lstlisting}
Now the filter chain is applied to the input images:
\begin{lstlisting}
vector<string>::const_iterator filter_name = filter_chain.begin();
for (list<C2DImageFilterFactory::ProductPtr>::const_iterator
f = filters.begin();
f != filters.end(); ++f) {
for (C2DImageList::iterator i = in_image_list->begin();
i != in_image_list->end(); ++i)
(*f)->filter_inplace(*i);
}
\end{lstlisting}
Finally, the filtered images are stored:
\begin{lstlisting}
imageio.save(out_type, *in_image_list, out_filename);
\end{lstlisting}
All the error checking has been removed in this display. At the core,
\texttt{mona-2dimagefilterstack} and \texttt{mona-2dstackfilter} provide
the same basic functionality. Therefore, adding filters and image
types to mona2d boils down to just writing new plugins. This can
be done without touching and/or recompiling the original software.
\section{Writing a new 2D Image IO Plugin}
\subsection*{Setting up the Plugin Frame}
The first thing to implement with a new 2D image IO plugin is, we need to include some headers:
\begin{lstlisting}
// include the compile time specifics
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
// some standard headers might be needed
#include <stdexcept>
#include <cassert>
// include all the basic definitions for 2D images
#include <libmona/2DImageWrap.hh>
// include the test routines
#include <mona2d/testImageIO.hh>
\end{lstlisting}
\noindent
Then we can define the namespace, in which the plugin code will reside (to avoid run-time name conflicts),
and define the name of the plugin, which should be equal to the file extension.
\begin{lstlisting}
namespace example_2dimage_io {
using namespace std;
using namespace mona;
static char const *format = "example";
//further code goes here
..
//
}
\end{lstlisting}
With this setup, the image load/save class can be defined as a derivative of the abstract class \texttt{C2DImageListIO}:
\begin{lstlisting}
class CExample2DImageIO: public C2DImageListIO {
public:
CExample2DImageIO();
virtual C2DImageList *load(string const& filename,bool,
CProgressCallback *cb)const;
virtual bool save(const C2DImageList& data,
string const& filename,
bool, CProgressCallback *cb) const;
const string short_descr()const;
private:
virtual int do_test()const;
};
\end{lstlisting}
Now lets take a closer look at the class declaration:
\texttt{CExample2DImageIO} is the constructor supposed to initialize the plugin, \texttt{load} and \texttt{save} are the
actual image IO routines.
In the base class \texttt{C2DImageListIO} these are abstract functions, and here we have to put life into them%
\footnote{In future versions these routines may be renamed to something like \texttt{do\_{*}} and moved to the private part of the class. %
The base class would then define \texttt{load} and \texttt{save} as non-virtual functions to separate the interface better %
from its implementation}
\texttt{short\_descr} is supposed to return a description of the plugin and \texttt{do\_test} is responsible for the plugin self test.
Now we will go for the implementation of these methods.
The constructor initializes the base class with the plugin name:
\begin{lstlisting}
CExample2DImageIO::CExample2DImageIO():
C2DImageListIO(format)
{
}
\end{lstlisting}
The file loader takes the filename, a flag indicating whether the input is actually read from stdin, an a call-back function that
can be used to show the loading progress.
For 2D images this is usually not useful, because the loading times are small.
\begin{lstlisting}
C2DImageList *CExample2DImageIO::load(string const& filename,
bool from_stdin,
CProgressCallback *cb) const
{
// open the file, and if it does not exist,
// then return without loading something.
CInputFile f(filename, from_stdin);
if (!f)
return NULL;
// create the result image list
C2DImageList *result = new C2DImageList;
// now the actual image loading takes place -
// this is very file type dependent
..
// return the loading result;
return result:
}
\end{lstlisting}
The standard saving procedure consists of two parts, a save functor to handle the different image types, and
the implementation of the C2DImageListIO::save routine:
\begin{lstlisting}
struct CImageSaver:public TUnaryImageFilter<bool> {
CImageSaver(CFile& f):
_M_f(f)
{
}
template <typename Data2D>
result_type operator ()(const Data2D& image);
private:
CFile& _M_f;
};
\end{lstlisting}
\noindent
The save functor is initialized with the output file, the actual savings done in the templated \texttt{operator ()}.
Its definition is omitted here, because it depends heavily on the target file format.
Finally, the \texttt{CExample2DImageIO::save} opens the output file, sets up the save functor and calls it:
\begin{lstlisting}
bool CExample2DImageIO::save(const C2DImageList& data,
string const& filename,
bool to_stdout,
CProgressCallback *cb) const
{
COutputFile f(filename, to_stdout);
if (!f)
return false;
CImageSaver saver(f);
for (C2DImageList::const_iterator iimg = data.begin();
iimg != data.end(); ++iimg)
wrap_filter(saver, *iimg);
return true;
}
\end{lstlisting}
\noindent
Not all output formats support storing multiple files, this should be handled here.
The core of the IO self test is implemented in the \texttt{test\_imageio\_plugin} function.
All is needed to do a self-test of the plugin, is to define which image types are supported and call the test function:
\begin{lstlisting}
int CExample2DImageIO::do_test()const
{
set<EPixelType> pixeltypes;
// list all pixel types supported by this output format
pixeltypes.insert(it_ushort);
pixeltypes.insert(it_ubyte);
pixeltypes.insert(it_bit);
return test_imageio_plugin(*this, pixeltypes) ? 1 : -1;
}
\end{lstlisting}
\noindent
Finally, in order to make the plugin loadable at run-time and callable from the main application, the interface function needs
to be implemented; it returns an instance of the image IO handler:
\begin{lstlisting}
extern "C" CPluginBase *get_plugin_interface()
{
return new CExample2DImageIO;
}
\end{lstlisting}
In order to implement the parts of a real image IO plugin that are not discussed here, a look at the implementation of the
other IO plugins may be helpful.
\section{Writing a new 2D Filter Plugin}
\label{sec:2dfilterplug}
A filter plugin consists of four basic parts:
\begin{enumerate}
\item \emph{The filter class} as the work horse of the filter as a specialized derivative of the templated \texttt{TUnaryImageFilter}.
It must be able to deal with different pixel types, therefore, the main filter routine is always templated with respect to the pixel type.
\item \emph{The filter dispatch class} derived from \texttt{C2DImageFilterBase}
dispatches the filter calls to the work horse by defining the inherited pure virtual interface.
Splitting this class from the above one is done,
because C++ does not provide the possibility to define virtual member templates%
\footnote{One can put the two together by using multiple inheritance. However, providing two classes results in a cleaner interface%
}.
\item \emph{The filter factory} that provides the callable interface to create filters as needed.
\item One standard entry point function to get the interface to the filter factory.
\end{enumerate}
Below, we will use the example of an intensity band pass filter to explain how these parts above are implemented.
\subsection*{The formal stuff }
The plugin implementation file will be called \texttt{bandpass\_2dimage\_filter.cc.}
Before we start with the implementation of above parts, the right environment has to be created.
We need to include the definition of the \texttt{C2DImageFilterBase}:
\begin{lstlisting}
#include <libmona/filter_plugin2d.hh>
\end{lstlisting}
and we will need some more of the standard definitions provided by
C++, in order to create the plugin:
\begin{lstlisting}
#include <limits>
\end{lstlisting}
Then it is best to put the whole code in an extra namespace, to avoid naming collisions with other plugins.
This could also achieved with proper export attributes.
However, they may not be portable across platforms.
\begin{lstlisting}
namespace bandpass_2dimage_filter {
//Further code goes here
...
//
}
\end{lstlisting}
All code that we will add later will be put within the boundaries of the namespace, unless otherwise noted.
Of course, if we may want to use some other namespaces we have to specify them:
\begin{lstlisting}
using namespace mona;
using namespace std;
\end{lstlisting}
\subsection*{Defining and Implementing the Workhorse Class}
Now we are ready to define the workhorse class:
\begin{lstlisting}
class CBandPass: public TUnaryImageFilter<C2DImageWrap> {
public:
CBandPass(float min, float max);
template <class Data2D>
typename CBandPass::result_type
operator () (const Data2D& data) const;
private:
float _M_min;
float _M_max;
};
\end{lstlisting}
First, \texttt{CBandPass} is derived from \texttt{TUnaryImageFilter<C2DImageWrap>} that provides
the typedef for \texttt{result\_type}.
Then a constructor is needed to initialize the the boundaries, that will later be used for
filtering the image.
Finally, the templated \texttt{operator ()} provides the interface that will be used to call the filter. Let's have a look at its implementation:
\begin{lstlisting}
template <class Data2D>
typename CBandPass::result_type
CBandPass::operator () (const Data2D& data) const
{
Data2D *result = new Data2D(data.get_size());
// do the filtering
...
return C2DImageWrap(result);
}
\end{lstlisting}
At first, a result image of the same type and size as the input image is created.
This might be different for other filters.
For example, a binarize filter would return a bit valued image.
Then the actual filter is applied (code omitted), and finally, the resulting image is wrapped in a way
that makes it possible to pass the image around without templating everything,
and the image is returned in that way.
The implementation of the constructor is straightforward, and needs
no further explanation:
\begin{lstlisting}
CBandPass::CBandPass(float min, float max):
_M_min(min),
_M_max(max)
{}
\end{lstlisting}
\subsection*{The dispatch class}
The dispatch class is derived from the abstract base class \texttt{C2DImageFilterBase} and implements its pure virtual
function that is used to call the actual filter.
\begin{lstlisting}
class C2DBandPassImageFilter: public C2DImageFilterBase {
CBandPass _M_filter;
public:
C2DBandPassImageFilter(float min, float max);
virtual C2DImageWrap do_filter(const C2DImageWrap& image) const;
};
\end{lstlisting}
The constructor of the class needs to initialize the work horse class:
\begin{lstlisting}
C2DBandPassImageFilter::C2DBandPassImageFilter(float min, float max):
_M_filter(min,max)
{}
\end{lstlisting}
and the \texttt{do\_filter} function applies some hidden magic, to call the work horse class operator () with the propper image type:
\begin{lstlisting}
C2DImageWrap C2DBandPassImageFilter::do_filter(const C2DImageWrap& image) const
{
return wrap_filter(_M_filter, image);
}
\end{lstlisting}
The magic is all hidden in \texttt{wrap\_filter} that checks of which pixel type the actual image is in order to
call the appropriate version of the filter.
\subsection*{The filter factory}
In order to define the filter factory some additional constants need
to be provided:
\begin{lstlisting}
static const char *plugin_name ="bandpass";
static const CFloatOption param_min("min", "lower bound of filter", 0,
-numeric_limits<float>::max(),
numeric_limits<float>::max());
static const CFloatOption param_max("max", "upper bound of filter",
numeric_limits<float>::max(),
-numeric_limits<float>::max(),
numeric_limits<float>::max());
\end{lstlisting}
First, we define the name of the plugin \texttt{``bandpass''} that will be used to specify the filter at the command line,
and then we define the two parameters of the filter, the lower and upper bound.
Now, we can have a look at the definition of the filter factory:
\begin{lstlisting}
class C2DBandPassImageFilterFactory: public C2DImageFilterFactory {
public:
C2DBandPassImageFilterFactory();
virtual C2DImageFilterFactory::ProductPtr
create(const CParsedOptions& options) const;
virtual const string short_descr()const;
private:
virtual int do_test() const;
};
\end{lstlisting}
This is not very special, more interesting is the implementation.
First the constructor:
\begin{lstlisting}
C2DBandPassImageFilterFactory::C2DBandPassImageFilterFactory():
C2DImageFilterFactory(plugin_name)
{
add_help(param_min);
add_help(param_max);
}
\end{lstlisting}
The base class is initialized with the plugin name, and in the body of the constructor, the filter parameter
are added to the plugin help.
Next we turn our focus at the filter creator function:
\begin{lstlisting}
C2DImageFilterFactory::ProductPtr
C2DBandPassImageFilterFactory::create(const CParsedOptions& options)
const
{
float min = param_min.get_value(options);
float max = param_max.get_value(options);
return C2DImageFilterFactory::ProductPtr(
new C2DBandPassImageFilter(min, max));
}
\end{lstlisting}
Here, the the values of the parameters are obtained from the options that are provided by the caller of this filter creator.
Then the Filter is created and passed back to the caller\texttt{}%
\footnote{\texttt{C2DImageFilterFactory::ProductPtr} is actually a \texttt{shared\_ptr}
from the boost library that helps to avoid memory leaks.%
}.
The \texttt{short\_descr} method just returns some string to provide a basic information about the plugin, and finally,
the do\_test method of the \texttt{C2DBandPassImageFilterFactory} provides some test routines,
to check, whether the filter does the right thing.
\subsection*{The Standard Entry Point}
The last thing to implement is the standard entry point function:
\begin{lstlisting}
extern "C" CPluginBase *get_plugin_interface()
{
return new C2DBandPassImageFilterFactory();
}
\end{lstlisting}
This function is declared as \texttt{extern {}``C''} since the plugin handler obtains its address by using its name,
and C++ linkage would result in name mangling, making the actual function name compiler dependant and therefore more or less unpredictable.
\subsection*{Compiling and linking}
In the last step, the dynamic library comprising the plugin needs to be build.
All plugins follow a specific naming schema. For 2d image filter, the final library must match
the pattern {}``\texttt{{*}\_2dimage\_filter.so}''.
If the plugin is build within the mona2d build tree, then the file bandpass\_2dimage\_filter.cc needs to be placed into the directory
\emph{plugins}, and the following additions are to be made to the
\emph{Makefile.am}:
\begin{itemize}
\item To the list \texttt{plugin\_2DIMAGE\_FILTERS} add \texttt{bandpass\_2dimage\_filter.cc}
\item Add the following lines somewhere in the \emph{Makefile.am}:
\end{itemize}
\begin{lstlisting}
bandpass_2dimage_filter_la_SOURCES = bandpass_2dimage_filter.cc
bandpass_2dimage_filter_la_LDFLAGS = -module -avoid-version
\end{lstlisting}
If the build was configured by using the configure command switch --enable-maintainer-mode then the
next execution of \emph{make} will update the Makefile accordingly, and the new filter will be compiled and linked.
After installation it can be used like all the other filters.
If maintainer mode is not enabled, change to the project root directory and run
\lstset{language=bash}
\begin{lstlisting}
automake --gnu
./config.status plugins/Makefile
\end{lstlisting}
In case the new filter is created separately, within its own project please refer to Chapter (FIXME),
in order to learn how to create a proper \emph{Makefile}.
\section{Writing a new stack filter plugin}
Stack based filtering is used to be able to apply 3D filters to very large data sets,
where only a few slices of the data fit into the working memory.
Therefore, implementing such a filter requires good planning on how the buffering should be done.
The best strategy to write a stack filter is first, to write a test function
that uses a simple 3D filter on a small data set, and then implement the real stack filter,
until its results match the result of the (easier to implement) real 3D filter.
In the following we will discuss the implementation of a 3D median filter.
We will skip some details, like test functions and description implementation.
\subsection*{The basic setup}
Just like the 2D filter, a 2D stack filter consists of 4 crucial parts, the \emph{work horse class},
the \emph{data dispatcher} , the \emph{filter factory}, and the \emph{plugin interface function}.
However, before these come into play, we need to setup the environment:
\begin{lstlisting}
// compile time specifics
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
// need some system headers
#include <iomanip>
#include <limits>
// need the definitions for the stack filtering
#include <mona2d/filter2dstack.hh>
// open a namespace
namespace median_2dstack_filter {
// use some provided namespaces
using namespace mona;
using namespace std;
// additional code goes here
..
}// end the namespace
\end{lstlisting}
\noindent
Explaining this filter is easier with a top-down approach.
Hence we will start with the interface function and wor down our way to the filter class:
\begin{lstlisting}
extern "C" CPluginBase *get_plugin_interface()
{
return new C2DMedianStackFilterFactory();
}
\end{lstlisting}
\noindent
First, we need to define the filter factory. A absolute minimum that needs to be declared and implement is given below:
\begin{lstlisting}
class C2DMedianStackFilterFactory:
public C2DStackFilterFactoryWithTest {
public:
C2DMedianStackFilterFactory();
virtual C2DStackFilterFactory::ProductPtr
create(const CParsedOptions& options) const;
};
static char const * plugin_name = "median";
static const CIntOption param_hw("w", "half filter width",
1, 1, numeric_limits<int>::max());
C2DMedianStackFilterFactory::C2DMedianStackFilterFactory():
C2DStackFilterFactoryWithTest(plugin_name)
{
add_help(param_hw);
}
C2DStackFilterFactory::ProductPtr
C2DMedianStackFilterFactory::create(const CParsedOptions& options)
const
{
int hw = param_hw.get_value(options);
return C2DStackFilterFactory::ProductPtr(
new C2DMedianStackFilter(hw));
}
\end{lstlisting}
\noindent
If you have read Section \ref{sec:2dfilterplug} above lines should need no further explanations.
Now we go for the filter dispatch class \texttt{C2DMedianStackFilter} that implements the
\texttt{C2DStackFilterBase} interface.
\begin{lstlisting}
class C2DMedianStackFilter: public C2DStackFilterBase {
C2DStackMedian *_M_filter;
size_t _M_hw;
public:
C2DMedianStackFilter(int hw);
virtual ~C2DMedianStackFilter();
private:
virtual void prepare(C2DImageWrap const& image);
virtual void do_push(C2DImageWrap const& image);
virtual bool do_filter(C2DImageWrap& dest,
C2DImageWrap const& src, size_t start, size_t end);
virtual void cleanup();
};
\end{lstlisting}
\noindent
We now look at the method implementation and discuss the details - first the constructor:
\begin{lstlisting}
C2DMedianStackFilter::C2DMedianStackFilter(int w):
C2DStackFilterBase(2 * w + 1, w + 1, w + 1),
_M_filter(NULL),
_M_w(w)
{
}
\end{lstlisting}
\noindent
The class is initialised with the filter width parameter $w$ - the actual filter covers a $(2w+1)^3$ cube.
The base class gets initalised with the size of the slice buffer needed $2w+1$, the number of slices needed, until
the first output slice can be generated $w+1$ and the number of (remaining) slices in the buffer $w+1$,
needed to generate the last output slice.
The filter \texttt{\_M\_filter} itself will not yet be initialized, since the input image size and the pixel type is not yet known.
Finally, the filter width parameter is stored.
The destructor just calls the cleanup function that might also be issued for other reasons.
Cleanto destroys the filter and frees its memory:
\begin{lstlisting}
C2DMedianStackFilter::~C2DMedianStackFilter()
{
cleanup();
}
void C2DMedianStackFilter::cleanup()
{
delete _M_filter;
_M_filter = NULL;
}
\end{lstlisting}
\noindent
\texttt{C2DMedianStackFilter::prepare} is called with a prototype image and it is used to create the actual filter:
\begin{lstlisting}
void C2DMedianStackFilter::prepare(C2DImageWrap const& image)
{
_M_filter = new C2DStackMedian(_M_hw, image);
}
\end{lstlisting}
\noindent
The inner workings of the filter class \texttt{C2DStackMedian} will be discussed later.
Filtering an image stack is then done in a tweo step procedure per input slice:
\begin{enumerate}
\item The input image is pushed into the buffer
\item The filter function is called to get a result slice
\end{enumerate}
The two functions are outlined below:
\begin{lstlisting}
void C2DMedianStackFilter::do_push(C2DImageWrap const& image)
{
wrap_filter(*_M_filter, image);
}
bool C2DMedianStackFilter::do_filter(C2DImageWrap& dest,
C2DImageWrap const& src, size_t start, size_t n)
{
if (n == 2 * _M_w + 1)
_M_filter->set_slice_range(0, n);
else {
if (start == 0)
_M_filter->set_slice_range(2*_M_w+1-n, 2*_M_w+1);
else
_M_filter->set_slice_range(start, start + n);
}
return wrap_filter(*_M_filter, dest);
}
\end{lstlisting}
\noindent
\texttt{C2DMedianStackFilter::do\_push} simply sends the image to the filter by unwrapping it based on pixel type.
The method \texttt{C2DMedianStackFilter::do\_filter}, on the other hand, first sets the valid slice range
depending on the number of slices $n$ that should be used for filtering.
Then the result is requested from the filter.
Note, both function use \texttt{wrap\_filter} to call the filter, however, depending on whether the \texttt{C2DImageWrap}
is passed as constant reference (\texttt{do\_push}) or as writable reference (\texttt{do\_filter}), the apropriate
operator of the filter will be called.
So let's now turn to the actual filter:
\begin{lstlisting}
class C2DStackMedian: public TUnaryImageFilter<bool> {
public:
C2DStackMedian(int hw, C2DImageWrap const& image);
template <typename Data2D>
typename C2DStackMedian::result_type
operator () (Data2D const& src);
template <typename Data2D>
typename C2DStackMedian::result_type
operator () (Data2D& result);
void set_slice_range(size_t start, size_t end);
private:
size_t _M_hw;
size_t _M_start;
size_t _M_end;
size_t _M_bufsize;
C3DImageWrap _M_buffer;
};
\end{lstlisting}
\noindent
The constructor initialises the filter width parameter and the slice range.
Then an 3D image is created that has the same width and height as the input image and the number of slices size needed to
hold everything that is needed for the buffer.
If you are interested in the details, you may want to look up the class definition of \texttt{C3DImgCreator} in
\texttt{filter2dstack.hh}.
\begin{lstlisting}
C2DStackMedian::C2DStackMedian(int w, C2DImageWrap const& image):
_M_hw(w),
_M_start(0),
_M_end(0),
_M_bufsize((2*w+1)*(2*w+1)*(2*w+1))
{
C3DImgCreator c(C3DBounds(image.get_size().x,
image.get_size().y, 2 * w + 1));
_M_buffer = wrap_filter(c, image);
}
\end{lstlisting}
\noindent
With the filter creared, \texttt{C2DMedianStackFilter} will start to push slices through
its \texttt{do\_push} method.
Since in this method onyl holds a constant reference to the input image the images are constant,
the \texttt{C2DStackMedian::operator()} will be called that takes a constant input image:
\begin{lstlisting}
template <typename T>
typename C2DStackMedian::result_type
C2DStackMedian::operator () (T2DImage<T> const& src)
{
T3DImage<T> *buf = _M_buffer.template get<T>();
copy(buf->begin_at(0,0,1),
buf->begin_at(0,0,buf->get_size().z),
buf->begin());
copy(src.begin(),
src.end(),
buf->begin_at(0,0,buf->get_size().z - 1));
return true;
}
\end{lstlisting}
\noindent
First, we get a pointer to the actual buffer with the right pixel type.
Then all the data is shifted one slice forward, and finally, the new slice is copied into the last slice of the buffer.
After the minimum amount of slices has been pushed down this pipeline, the filtering mechanism implemented in the library will start to issue
calls to \linebreak[3]\texttt{C2DMedianStackFilter::do\_filter} to get the filtered data.
Now that version of \linebreak[2]\texttt{C2DStackMedian::operator()} is called that uses a writable input image reference.
Here, for each pixel of the result image, the median filtered value is obtained from
the slice buffer:
\begin{lstlisting}
template <typename T>
typename C2DStackMedian::result_type
C2DStackMedian::operator () (T2DImage<T>& result)
{
T3DImage<T> *src = _M_buffer.template get<T>();
typename T2DImage<T>::iterator i = result.begin();
for (size_t y = 0; y < result.get_size().y; ++y)
for (size_t x = 0; x < result.get_size().x; ++x, ++i)
*i = median(*src, x, y, _M_start, _M_end, _M_hw);
return true;
}
\end{lstlisting}
\noindent
The implementation of the \texttt{median} function is not dependend on the fact that we discuss stack based
filtering here, and therefore its implementation will not be discussed.
\bibliographystyle{plainnat}
\cleardoublepage\addcontentsline{toc}{chapter}{\bibname}\bibliography{segment}
\end{document}