Initial release of the baltrad rack software.
authorAnders Henja <anders@baltrad.eu>
Wed, 24 Aug 2011 12:09:12 +0000 (14:09 +0200)
committerAnders Henja <anders@baltrad.eu>
Wed, 24 Aug 2011 12:09:12 +0000 (14:09 +0200)
163 files changed:
.gitignore [new file with mode: 0644]
COPYING [new file with mode: 0644]
COPYING.LESSER [new file with mode: 0644]
def.mk [new file with mode: 0644]
drain/Makefile [new file with mode: 0644]
drain/image/.dep/CoordinateHandler.P [new file with mode: 0644]
drain/image/.dep/FileBinary.P [new file with mode: 0644]
drain/image/.dep/FileFQD.P [new file with mode: 0644]
drain/image/.dep/Geometry.P [new file with mode: 0644]
drain/image/CatenatorOp.h [new file with mode: 0644]
drain/image/CoordinateHandler.cpp [new file with mode: 0644]
drain/image/CoordinateHandler.h [new file with mode: 0644]
drain/image/CopyOp.h [new file with mode: 0644]
drain/image/Cumulator.h [new file with mode: 0644]
drain/image/DistanceTransformFillOp.h [new file with mode: 0644]
drain/image/DistanceTransformOp.h [new file with mode: 0644]
drain/image/DoubleSmootherOp.h [new file with mode: 0644]
drain/image/Drain.h [new file with mode: 0644]
drain/image/FastAverageOp.h [new file with mode: 0644]
drain/image/FastOpticalFlowOp.h [new file with mode: 0644]
drain/image/File.h [new file with mode: 0644]
drain/image/FileBinary.cpp [new file with mode: 0644]
drain/image/FileBinary.h [new file with mode: 0644]
drain/image/FileFQD.cpp [new file with mode: 0644]
drain/image/FileFQD.h [new file with mode: 0644]
drain/image/FilePng.h [new file with mode: 0644]
drain/image/FloodFillOp.h [new file with mode: 0644]
drain/image/FuzzyPeakOp.h [new file with mode: 0644]
drain/image/FuzzyThresholdOp.h [new file with mode: 0644]
drain/image/GammaOp.h [new file with mode: 0644]
drain/image/Geometry.cpp [new file with mode: 0644]
drain/image/Geometry.h [new file with mode: 0644]
drain/image/GradientOp.h [new file with mode: 0644]
drain/image/HighBoostOp.h [new file with mode: 0644]
drain/image/HighPassOp.h [new file with mode: 0644]
drain/image/Image.h [new file with mode: 0644]
drain/image/ImageOp.h [new file with mode: 0644]
drain/image/ImageView.h [new file with mode: 0644]
drain/image/Intensity.h [new file with mode: 0644]
drain/image/MagickDrain.h [new file with mode: 0644]
drain/image/MarginalStatisticOp.h [new file with mode: 0644]
drain/image/MathOpPack.h [new file with mode: 0644]
drain/image/PaletteOp.h [new file with mode: 0644]
drain/image/PixelVectorOp.h [new file with mode: 0644]
drain/image/Point.h [new file with mode: 0644]
drain/image/QuadraticSmootherOp.h [new file with mode: 0644]
drain/image/QualityOverrideOp.h [new file with mode: 0644]
drain/image/QuantizatorOp.h [new file with mode: 0644]
drain/image/RecursiveRepairerOp.h [new file with mode: 0644]
drain/image/RunLengthOp.h [new file with mode: 0644]
drain/image/SegmentAreaOp.h [new file with mode: 0644]
drain/image/SequentialImageOp.h [new file with mode: 0644]
drain/image/SlidingStripeOp.h [new file with mode: 0644]
drain/image/SlidingWindow.h [new file with mode: 0644]
drain/image/SlidingWindowHistogram.h [new file with mode: 0644]
drain/image/SlidingWindowMedianOp.h [new file with mode: 0644]
drain/image/SlidingWindowOp.h [new file with mode: 0644]
drain/image/ThresholdOp.h [new file with mode: 0644]
drain/image/Window.h [new file with mode: 0644]
drain/image/WindowOp.h [new file with mode: 0644]
drain/radar/.dep/Composite.P [new file with mode: 0644]
drain/radar/.dep/Coordinates.P [new file with mode: 0644]
drain/radar/.dep/Geometry.P [new file with mode: 0644]
drain/radar/.dep/PolarCoordinateHandler.P [new file with mode: 0644]
drain/radar/.dep/ProductXReader.P [new file with mode: 0644]
drain/radar/.dep/SubComposite.P [new file with mode: 0644]
drain/radar/Andre.h [new file with mode: 0644]
drain/radar/Composite.cpp [new file with mode: 0644]
drain/radar/Composite.h [new file with mode: 0644]
drain/radar/Constants.h [new file with mode: 0644]
drain/radar/Coordinates.cpp [new file with mode: 0644]
drain/radar/Coordinates.h [new file with mode: 0644]
drain/radar/Geometry.cpp [new file with mode: 0644]
drain/radar/Geometry.h [new file with mode: 0644]
drain/radar/PolarAttenuation.h [new file with mode: 0644]
drain/radar/PolarCappi.h [new file with mode: 0644]
drain/radar/PolarCoordinateHandler.cpp [new file with mode: 0644]
drain/radar/PolarCoordinateHandler.h [new file with mode: 0644]
drain/radar/PolarEchoTop.h [new file with mode: 0644]
drain/radar/PolarMaxEcho.h [new file with mode: 0644]
drain/radar/PolarProduct.h [new file with mode: 0644]
drain/radar/PolarToCartesian.h [new file with mode: 0644]
drain/radar/PolarToGeographical.h [new file with mode: 0644]
drain/radar/ProductXReader.cpp [new file with mode: 0644]
drain/radar/ProductXReader.h [new file with mode: 0644]
drain/radar/SubComposite.cpp [new file with mode: 0644]
drain/radar/SubComposite.h [new file with mode: 0644]
drain/util/.dep/Data.P [new file with mode: 0644]
drain/util/.dep/Debug.P [new file with mode: 0644]
drain/util/.dep/Options.P [new file with mode: 0644]
drain/util/.dep/Proj4.P [new file with mode: 0644]
drain/util/.dep/ProjectionFrame.P [new file with mode: 0644]
drain/util/.dep/RegExp.P [new file with mode: 0644]
drain/util/.dep/String.P [new file with mode: 0644]
drain/util/Data.cpp [new file with mode: 0644]
drain/util/Data.h [new file with mode: 0644]
drain/util/Debug.cpp [new file with mode: 0644]
drain/util/Debug.h [new file with mode: 0644]
drain/util/Histogram.h [new file with mode: 0644]
drain/util/MapReader.h [new file with mode: 0644]
drain/util/MapWrapper.h [new file with mode: 0644]
drain/util/MapWriter.h [new file with mode: 0644]
drain/util/Options.cpp [new file with mode: 0644]
drain/util/Options.h [new file with mode: 0644]
drain/util/Proj4.cpp [new file with mode: 0644]
drain/util/Proj4.h [new file with mode: 0644]
drain/util/ProjectionFrame.cpp [new file with mode: 0644]
drain/util/ProjectionFrame.h [new file with mode: 0644]
drain/util/Rectangle.h [new file with mode: 0644]
drain/util/RegExp.cpp [new file with mode: 0644]
drain/util/RegExp.h [new file with mode: 0644]
drain/util/String.cpp [new file with mode: 0644]
drain/util/String.h [new file with mode: 0644]
drain/util/StringMapper.h [new file with mode: 0644]
drain/util/Time.h [new file with mode: 0644]
drain/util/TreeNode.h [new file with mode: 0644]
pyrack/pyrack.c [new file with mode: 0644]
rack/Makefile [new file with mode: 0644]
rack/bin/rack.cpp [new file with mode: 0644]
rack/bin/ropo.c [new file with mode: 0644]
rack/bin/ropo_help.c [new file with mode: 0644]
rack/hi5/DataScaleOp.h [new file with mode: 0644]
rack/hi5/RaveConvert.cpp [new file with mode: 0644]
rack/hi5/RaveConvert.h [new file with mode: 0644]
rack/rack/Rack.cpp [new file with mode: 0644]
rack/rack/Rack.h [new file with mode: 0644]
rack/rack/RackIf.cpp [new file with mode: 0644]
rack/rack/RackIf.h [new file with mode: 0644]
rack/rack/RackLocalIf.h [new file with mode: 0644]
rack/ropo/Ropo.cpp [new file with mode: 0644]
rack/ropo/Ropo.h [new file with mode: 0644]
rack/ropo/convertImage.cpp [new file with mode: 0644]
rack/ropo/convertImage.h [new file with mode: 0644]
rack/ropo/fmi_image.c [new file with mode: 0644]
rack/ropo/fmi_image.h [new file with mode: 0644]
rack/ropo/fmi_image_arith.c [new file with mode: 0644]
rack/ropo/fmi_image_arith.h [new file with mode: 0644]
rack/ropo/fmi_image_filter.c [new file with mode: 0644]
rack/ropo/fmi_image_filter.h [new file with mode: 0644]
rack/ropo/fmi_image_filter_line.c [new file with mode: 0644]
rack/ropo/fmi_image_filter_line.h [new file with mode: 0644]
rack/ropo/fmi_image_filter_morpho.c [new file with mode: 0644]
rack/ropo/fmi_image_filter_morpho.h [new file with mode: 0644]
rack/ropo/fmi_image_filter_speck.c [new file with mode: 0644]
rack/ropo/fmi_image_filter_speck.h [new file with mode: 0644]
rack/ropo/fmi_image_filter_texture.c [new file with mode: 0644]
rack/ropo/fmi_image_filter_texture.h [new file with mode: 0644]
rack/ropo/fmi_image_histogram.c [new file with mode: 0644]
rack/ropo/fmi_image_histogram.h [new file with mode: 0644]
rack/ropo/fmi_image_restore.c [new file with mode: 0644]
rack/ropo/fmi_image_restore.h [new file with mode: 0644]
rack/ropo/fmi_meteosat.c [new file with mode: 0644]
rack/ropo/fmi_meteosat.h [new file with mode: 0644]
rack/ropo/fmi_radar_codes.h [new file with mode: 0644]
rack/ropo/fmi_radar_codes.inc [new file with mode: 0755]
rack/ropo/fmi_radar_image.c [new file with mode: 0644]
rack/ropo/fmi_radar_image.h [new file with mode: 0644]
rack/ropo/fmi_sunpos.c [new file with mode: 0644]
rack/ropo/fmi_sunpos.h [new file with mode: 0644]
rack/ropo/fmi_util.c [new file with mode: 0644]
rack/ropo/fmi_util.h [new file with mode: 0644]
rack/ropo/ropo_hdf.c [new file with mode: 0644]
rack/ropo/ropo_hdf.h [new file with mode: 0644]

diff --git a/.gitignore b/.gitignore
new file mode 100644 (file)
index 0000000..ddbe1f4
--- /dev/null
@@ -0,0 +1,8 @@
+*.o
+*.so
+*.a
+.cproject
+.project
+Debug/
+*~
+
diff --git a/COPYING b/COPYING
new file mode 100644 (file)
index 0000000..94a9ed0
--- /dev/null
+++ b/COPYING
@@ -0,0 +1,674 @@
+                    GNU GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+                            Preamble
+
+  The GNU General Public License is a free, copyleft license for
+software and other kinds of works.
+
+  The licenses for most software and other practical works are designed
+to take away your freedom to share and change the works.  By contrast,
+the GNU General Public License is intended to guarantee your freedom to
+share and change all versions of a program--to make sure it remains free
+software for all its users.  We, the Free Software Foundation, use the
+GNU General Public License for most of our software; it applies also to
+any other work released this way by its authors.  You can apply it to
+your programs, too.
+
+  When we speak of free software, we are referring to freedom, not
+price.  Our General Public Licenses are designed to make sure that you
+have the freedom to distribute copies of free software (and charge for
+them if you wish), that you receive source code or can get it if you
+want it, that you can change the software or use pieces of it in new
+free programs, and that you know you can do these things.
+
+  To protect your rights, we need to prevent others from denying you
+these rights or asking you to surrender the rights.  Therefore, you have
+certain responsibilities if you distribute copies of the software, or if
+you modify it: responsibilities to respect the freedom of others.
+
+  For example, if you distribute copies of such a program, whether
+gratis or for a fee, you must pass on to the recipients the same
+freedoms that you received.  You must make sure that they, too, receive
+or can get the source code.  And you must show them these terms so they
+know their rights.
+
+  Developers that use the GNU GPL protect your rights with two steps:
+(1) assert copyright on the software, and (2) offer you this License
+giving you legal permission to copy, distribute and/or modify it.
+
+  For the developers' and authors' protection, the GPL clearly explains
+that there is no warranty for this free software.  For both users' and
+authors' sake, the GPL requires that modified versions be marked as
+changed, so that their problems will not be attributed erroneously to
+authors of previous versions.
+
+  Some devices are designed to deny users access to install or run
+modified versions of the software inside them, although the manufacturer
+can do so.  This is fundamentally incompatible with the aim of
+protecting users' freedom to change the software.  The systematic
+pattern of such abuse occurs in the area of products for individuals to
+use, which is precisely where it is most unacceptable.  Therefore, we
+have designed this version of the GPL to prohibit the practice for those
+products.  If such problems arise substantially in other domains, we
+stand ready to extend this provision to those domains in future versions
+of the GPL, as needed to protect the freedom of users.
+
+  Finally, every program is threatened constantly by software patents.
+States should not allow patents to restrict development and use of
+software on general-purpose computers, but in those that do, we wish to
+avoid the special danger that patents applied to a free program could
+make it effectively proprietary.  To prevent this, the GPL assures that
+patents cannot be used to render the program non-free.
+
+  The precise terms and conditions for copying, distribution and
+modification follow.
+
+                       TERMS AND CONDITIONS
+
+  0. Definitions.
+
+  "This License" refers to version 3 of the GNU General Public License.
+
+  "Copyright" also means copyright-like laws that apply to other kinds of
+works, such as semiconductor masks.
+
+  "The Program" refers to any copyrightable work licensed under this
+License.  Each licensee is addressed as "you".  "Licensees" and
+"recipients" may be individuals or organizations.
+
+  To "modify" a work means to copy from or adapt all or part of the work
+in a fashion requiring copyright permission, other than the making of an
+exact copy.  The resulting work is called a "modified version" of the
+earlier work or a work "based on" the earlier work.
+
+  A "covered work" means either the unmodified Program or a work based
+on the Program.
+
+  To "propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy.  Propagation includes copying,
+distribution (with or without modification), making available to the
+public, and in some countries other activities as well.
+
+  To "convey" a work means any kind of propagation that enables other
+parties to make or receive copies.  Mere interaction with a user through
+a computer network, with no transfer of a copy, is not conveying.
+
+  An interactive user interface displays "Appropriate Legal Notices"
+to the extent that it includes a convenient and prominently visible
+feature that (1) displays an appropriate copyright notice, and (2)
+tells the user that there is no warranty for the work (except to the
+extent that warranties are provided), that licensees may convey the
+work under this License, and how to view a copy of this License.  If
+the interface presents a list of user commands or options, such as a
+menu, a prominent item in the list meets this criterion.
+
+  1. Source Code.
+
+  The "source code" for a work means the preferred form of the work
+for making modifications to it.  "Object code" means any non-source
+form of a work.
+
+  A "Standard Interface" means an interface that either is an official
+standard defined by a recognized standards body, or, in the case of
+interfaces specified for a particular programming language, one that
+is widely used among developers working in that language.
+
+  The "System Libraries" of an executable work include anything, other
+than the work as a whole, that (a) is included in the normal form of
+packaging a Major Component, but which is not part of that Major
+Component, and (b) serves only to enable use of the work with that
+Major Component, or to implement a Standard Interface for which an
+implementation is available to the public in source code form.  A
+"Major Component", in this context, means a major essential component
+(kernel, window system, and so on) of the specific operating system
+(if any) on which the executable work runs, or a compiler used to
+produce the work, or an object code interpreter used to run it.
+
+  The "Corresponding Source" for a work in object code form means all
+the source code needed to generate, install, and (for an executable
+work) run the object code and to modify the work, including scripts to
+control those activities.  However, it does not include the work's
+System Libraries, or general-purpose tools or generally available free
+programs which are used unmodified in performing those activities but
+which are not part of the work.  For example, Corresponding Source
+includes interface definition files associated with source files for
+the work, and the source code for shared libraries and dynamically
+linked subprograms that the work is specifically designed to require,
+such as by intimate data communication or control flow between those
+subprograms and other parts of the work.
+
+  The Corresponding Source need not include anything that users
+can regenerate automatically from other parts of the Corresponding
+Source.
+
+  The Corresponding Source for a work in source code form is that
+same work.
+
+  2. Basic Permissions.
+
+  All rights granted under this License are granted for the term of
+copyright on the Program, and are irrevocable provided the stated
+conditions are met.  This License explicitly affirms your unlimited
+permission to run the unmodified Program.  The output from running a
+covered work is covered by this License only if the output, given its
+content, constitutes a covered work.  This License acknowledges your
+rights of fair use or other equivalent, as provided by copyright law.
+
+  You may make, run and propagate covered works that you do not
+convey, without conditions so long as your license otherwise remains
+in force.  You may convey covered works to others for the sole purpose
+of having them make modifications exclusively for you, or provide you
+with facilities for running those works, provided that you comply with
+the terms of this License in conveying all material for which you do
+not control copyright.  Those thus making or running the covered works
+for you must do so exclusively on your behalf, under your direction
+and control, on terms that prohibit them from making any copies of
+your copyrighted material outside their relationship with you.
+
+  Conveying under any other circumstances is permitted solely under
+the conditions stated below.  Sublicensing is not allowed; section 10
+makes it unnecessary.
+
+  3. Protecting Users' Legal Rights From Anti-Circumvention Law.
+
+  No covered work shall be deemed part of an effective technological
+measure under any applicable law fulfilling obligations under article
+11 of the WIPO copyright treaty adopted on 20 December 1996, or
+similar laws prohibiting or restricting circumvention of such
+measures.
+
+  When you convey a covered work, you waive any legal power to forbid
+circumvention of technological measures to the extent such circumvention
+is effected by exercising rights under this License with respect to
+the covered work, and you disclaim any intention to limit operation or
+modification of the work as a means of enforcing, against the work's
+users, your or third parties' legal rights to forbid circumvention of
+technological measures.
+
+  4. Conveying Verbatim Copies.
+
+  You may convey verbatim copies of the Program's source code as you
+receive it, in any medium, provided that you conspicuously and
+appropriately publish on each copy an appropriate copyright notice;
+keep intact all notices stating that this License and any
+non-permissive terms added in accord with section 7 apply to the code;
+keep intact all notices of the absence of any warranty; and give all
+recipients a copy of this License along with the Program.
+
+  You may charge any price or no price for each copy that you convey,
+and you may offer support or warranty protection for a fee.
+
+  5. Conveying Modified Source Versions.
+
+  You may convey a work based on the Program, or the modifications to
+produce it from the Program, in the form of source code under the
+terms of section 4, provided that you also meet all of these conditions:
+
+    a) The work must carry prominent notices stating that you modified
+    it, and giving a relevant date.
+
+    b) The work must carry prominent notices stating that it is
+    released under this License and any conditions added under section
+    7.  This requirement modifies the requirement in section 4 to
+    "keep intact all notices".
+
+    c) You must license the entire work, as a whole, under this
+    License to anyone who comes into possession of a copy.  This
+    License will therefore apply, along with any applicable section 7
+    additional terms, to the whole of the work, and all its parts,
+    regardless of how they are packaged.  This License gives no
+    permission to license the work in any other way, but it does not
+    invalidate such permission if you have separately received it.
+
+    d) If the work has interactive user interfaces, each must display
+    Appropriate Legal Notices; however, if the Program has interactive
+    interfaces that do not display Appropriate Legal Notices, your
+    work need not make them do so.
+
+  A compilation of a covered work with other separate and independent
+works, which are not by their nature extensions of the covered work,
+and which are not combined with it such as to form a larger program,
+in or on a volume of a storage or distribution medium, is called an
+"aggregate" if the compilation and its resulting copyright are not
+used to limit the access or legal rights of the compilation's users
+beyond what the individual works permit.  Inclusion of a covered work
+in an aggregate does not cause this License to apply to the other
+parts of the aggregate.
+
+  6. Conveying Non-Source Forms.
+
+  You may convey a covered work in object code form under the terms
+of sections 4 and 5, provided that you also convey the
+machine-readable Corresponding Source under the terms of this License,
+in one of these ways:
+
+    a) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by the
+    Corresponding Source fixed on a durable physical medium
+    customarily used for software interchange.
+
+    b) Convey the object code in, or embodied in, a physical product
+    (including a physical distribution medium), accompanied by a
+    written offer, valid for at least three years and valid for as
+    long as you offer spare parts or customer support for that product
+    model, to give anyone who possesses the object code either (1) a
+    copy of the Corresponding Source for all the software in the
+    product that is covered by this License, on a durable physical
+    medium customarily used for software interchange, for a price no
+    more than your reasonable cost of physically performing this
+    conveying of source, or (2) access to copy the
+    Corresponding Source from a network server at no charge.
+
+    c) Convey individual copies of the object code with a copy of the
+    written offer to provide the Corresponding Source.  This
+    alternative is allowed only occasionally and noncommercially, and
+    only if you received the object code with such an offer, in accord
+    with subsection 6b.
+
+    d) Convey the object code by offering access from a designated
+    place (gratis or for a charge), and offer equivalent access to the
+    Corresponding Source in the same way through the same place at no
+    further charge.  You need not require recipients to copy the
+    Corresponding Source along with the object code.  If the place to
+    copy the object code is a network server, the Corresponding Source
+    may be on a different server (operated by you or a third party)
+    that supports equivalent copying facilities, provided you maintain
+    clear directions next to the object code saying where to find the
+    Corresponding Source.  Regardless of what server hosts the
+    Corresponding Source, you remain obligated to ensure that it is
+    available for as long as needed to satisfy these requirements.
+
+    e) Convey the object code using peer-to-peer transmission, provided
+    you inform other peers where the object code and Corresponding
+    Source of the work are being offered to the general public at no
+    charge under subsection 6d.
+
+  A separable portion of the object code, whose source code is excluded
+from the Corresponding Source as a System Library, need not be
+included in conveying the object code work.
+
+  A "User Product" is either (1) a "consumer product", which means any
+tangible personal property which is normally used for personal, family,
+or household purposes, or (2) anything designed or sold for incorporation
+into a dwelling.  In determining whether a product is a consumer product,
+doubtful cases shall be resolved in favor of coverage.  For a particular
+product received by a particular user, "normally used" refers to a
+typical or common use of that class of product, regardless of the status
+of the particular user or of the way in which the particular user
+actually uses, or expects or is expected to use, the product.  A product
+is a consumer product regardless of whether the product has substantial
+commercial, industrial or non-consumer uses, unless such uses represent
+the only significant mode of use of the product.
+
+  "Installation Information" for a User Product means any methods,
+procedures, authorization keys, or other information required to install
+and execute modified versions of a covered work in that User Product from
+a modified version of its Corresponding Source.  The information must
+suffice to ensure that the continued functioning of the modified object
+code is in no case prevented or interfered with solely because
+modification has been made.
+
+  If you convey an object code work under this section in, or with, or
+specifically for use in, a User Product, and the conveying occurs as
+part of a transaction in which the right of possession and use of the
+User Product is transferred to the recipient in perpetuity or for a
+fixed term (regardless of how the transaction is characterized), the
+Corresponding Source conveyed under this section must be accompanied
+by the Installation Information.  But this requirement does not apply
+if neither you nor any third party retains the ability to install
+modified object code on the User Product (for example, the work has
+been installed in ROM).
+
+  The requirement to provide Installation Information does not include a
+requirement to continue to provide support service, warranty, or updates
+for a work that has been modified or installed by the recipient, or for
+the User Product in which it has been modified or installed.  Access to a
+network may be denied when the modification itself materially and
+adversely affects the operation of the network or violates the rules and
+protocols for communication across the network.
+
+  Corresponding Source conveyed, and Installation Information provided,
+in accord with this section must be in a format that is publicly
+documented (and with an implementation available to the public in
+source code form), and must require no special password or key for
+unpacking, reading or copying.
+
+  7. Additional Terms.
+
+  "Additional permissions" are terms that supplement the terms of this
+License by making exceptions from one or more of its conditions.
+Additional permissions that are applicable to the entire Program shall
+be treated as though they were included in this License, to the extent
+that they are valid under applicable law.  If additional permissions
+apply only to part of the Program, that part may be used separately
+under those permissions, but the entire Program remains governed by
+this License without regard to the additional permissions.
+
+  When you convey a copy of a covered work, you may at your option
+remove any additional permissions from that copy, or from any part of
+it.  (Additional permissions may be written to require their own
+removal in certain cases when you modify the work.)  You may place
+additional permissions on material, added by you to a covered work,
+for which you have or can give appropriate copyright permission.
+
+  Notwithstanding any other provision of this License, for material you
+add to a covered work, you may (if authorized by the copyright holders of
+that material) supplement the terms of this License with terms:
+
+    a) Disclaiming warranty or limiting liability differently from the
+    terms of sections 15 and 16 of this License; or
+
+    b) Requiring preservation of specified reasonable legal notices or
+    author attributions in that material or in the Appropriate Legal
+    Notices displayed by works containing it; or
+
+    c) Prohibiting misrepresentation of the origin of that material, or
+    requiring that modified versions of such material be marked in
+    reasonable ways as different from the original version; or
+
+    d) Limiting the use for publicity purposes of names of licensors or
+    authors of the material; or
+
+    e) Declining to grant rights under trademark law for use of some
+    trade names, trademarks, or service marks; or
+
+    f) Requiring indemnification of licensors and authors of that
+    material by anyone who conveys the material (or modified versions of
+    it) with contractual assumptions of liability to the recipient, for
+    any liability that these contractual assumptions directly impose on
+    those licensors and authors.
+
+  All other non-permissive additional terms are considered "further
+restrictions" within the meaning of section 10.  If the Program as you
+received it, or any part of it, contains a notice stating that it is
+governed by this License along with a term that is a further
+restriction, you may remove that term.  If a license document contains
+a further restriction but permits relicensing or conveying under this
+License, you may add to a covered work material governed by the terms
+of that license document, provided that the further restriction does
+not survive such relicensing or conveying.
+
+  If you add terms to a covered work in accord with this section, you
+must place, in the relevant source files, a statement of the
+additional terms that apply to those files, or a notice indicating
+where to find the applicable terms.
+
+  Additional terms, permissive or non-permissive, may be stated in the
+form of a separately written license, or stated as exceptions;
+the above requirements apply either way.
+
+  8. Termination.
+
+  You may not propagate or modify a covered work except as expressly
+provided under this License.  Any attempt otherwise to propagate or
+modify it is void, and will automatically terminate your rights under
+this License (including any patent licenses granted under the third
+paragraph of section 11).
+
+  However, if you cease all violation of this License, then your
+license from a particular copyright holder is reinstated (a)
+provisionally, unless and until the copyright holder explicitly and
+finally terminates your license, and (b) permanently, if the copyright
+holder fails to notify you of the violation by some reasonable means
+prior to 60 days after the cessation.
+
+  Moreover, your license from a particular copyright holder is
+reinstated permanently if the copyright holder notifies you of the
+violation by some reasonable means, this is the first time you have
+received notice of violation of this License (for any work) from that
+copyright holder, and you cure the violation prior to 30 days after
+your receipt of the notice.
+
+  Termination of your rights under this section does not terminate the
+licenses of parties who have received copies or rights from you under
+this License.  If your rights have been terminated and not permanently
+reinstated, you do not qualify to receive new licenses for the same
+material under section 10.
+
+  9. Acceptance Not Required for Having Copies.
+
+  You are not required to accept this License in order to receive or
+run a copy of the Program.  Ancillary propagation of a covered work
+occurring solely as a consequence of using peer-to-peer transmission
+to receive a copy likewise does not require acceptance.  However,
+nothing other than this License grants you permission to propagate or
+modify any covered work.  These actions infringe copyright if you do
+not accept this License.  Therefore, by modifying or propagating a
+covered work, you indicate your acceptance of this License to do so.
+
+  10. Automatic Licensing of Downstream Recipients.
+
+  Each time you convey a covered work, the recipient automatically
+receives a license from the original licensors, to run, modify and
+propagate that work, subject to this License.  You are not responsible
+for enforcing compliance by third parties with this License.
+
+  An "entity transaction" is a transaction transferring control of an
+organization, or substantially all assets of one, or subdividing an
+organization, or merging organizations.  If propagation of a covered
+work results from an entity transaction, each party to that
+transaction who receives a copy of the work also receives whatever
+licenses to the work the party's predecessor in interest had or could
+give under the previous paragraph, plus a right to possession of the
+Corresponding Source of the work from the predecessor in interest, if
+the predecessor has it or can get it with reasonable efforts.
+
+  You may not impose any further restrictions on the exercise of the
+rights granted or affirmed under this License.  For example, you may
+not impose a license fee, royalty, or other charge for exercise of
+rights granted under this License, and you may not initiate litigation
+(including a cross-claim or counterclaim in a lawsuit) alleging that
+any patent claim is infringed by making, using, selling, offering for
+sale, or importing the Program or any portion of it.
+
+  11. Patents.
+
+  A "contributor" is a copyright holder who authorizes use under this
+License of the Program or a work on which the Program is based.  The
+work thus licensed is called the contributor's "contributor version".
+
+  A contributor's "essential patent claims" are all patent claims
+owned or controlled by the contributor, whether already acquired or
+hereafter acquired, that would be infringed by some manner, permitted
+by this License, of making, using, or selling its contributor version,
+but do not include claims that would be infringed only as a
+consequence of further modification of the contributor version.  For
+purposes of this definition, "control" includes the right to grant
+patent sublicenses in a manner consistent with the requirements of
+this License.
+
+  Each contributor grants you a non-exclusive, worldwide, royalty-free
+patent license under the contributor's essential patent claims, to
+make, use, sell, offer for sale, import and otherwise run, modify and
+propagate the contents of its contributor version.
+
+  In the following three paragraphs, a "patent license" is any express
+agreement or commitment, however denominated, not to enforce a patent
+(such as an express permission to practice a patent or covenant not to
+sue for patent infringement).  To "grant" such a patent license to a
+party means to make such an agreement or commitment not to enforce a
+patent against the party.
+
+  If you convey a covered work, knowingly relying on a patent license,
+and the Corresponding Source of the work is not available for anyone
+to copy, free of charge and under the terms of this License, through a
+publicly available network server or other readily accessible means,
+then you must either (1) cause the Corresponding Source to be so
+available, or (2) arrange to deprive yourself of the benefit of the
+patent license for this particular work, or (3) arrange, in a manner
+consistent with the requirements of this License, to extend the patent
+license to downstream recipients.  "Knowingly relying" means you have
+actual knowledge that, but for the patent license, your conveying the
+covered work in a country, or your recipient's use of the covered work
+in a country, would infringe one or more identifiable patents in that
+country that you have reason to believe are valid.
+
+  If, pursuant to or in connection with a single transaction or
+arrangement, you convey, or propagate by procuring conveyance of, a
+covered work, and grant a patent license to some of the parties
+receiving the covered work authorizing them to use, propagate, modify
+or convey a specific copy of the covered work, then the patent license
+you grant is automatically extended to all recipients of the covered
+work and works based on it.
+
+  A patent license is "discriminatory" if it does not include within
+the scope of its coverage, prohibits the exercise of, or is
+conditioned on the non-exercise of one or more of the rights that are
+specifically granted under this License.  You may not convey a covered
+work if you are a party to an arrangement with a third party that is
+in the business of distributing software, under which you make payment
+to the third party based on the extent of your activity of conveying
+the work, and under which the third party grants, to any of the
+parties who would receive the covered work from you, a discriminatory
+patent license (a) in connection with copies of the covered work
+conveyed by you (or copies made from those copies), or (b) primarily
+for and in connection with specific products or compilations that
+contain the covered work, unless you entered into that arrangement,
+or that patent license was granted, prior to 28 March 2007.
+
+  Nothing in this License shall be construed as excluding or limiting
+any implied license or other defenses to infringement that may
+otherwise be available to you under applicable patent law.
+
+  12. No Surrender of Others' Freedom.
+
+  If conditions are imposed on you (whether by court order, agreement or
+otherwise) that contradict the conditions of this License, they do not
+excuse you from the conditions of this License.  If you cannot convey a
+covered work so as to satisfy simultaneously your obligations under this
+License and any other pertinent obligations, then as a consequence you may
+not convey it at all.  For example, if you agree to terms that obligate you
+to collect a royalty for further conveying from those to whom you convey
+the Program, the only way you could satisfy both those terms and this
+License would be to refrain entirely from conveying the Program.
+
+  13. Use with the GNU Affero General Public License.
+
+  Notwithstanding any other provision of this License, you have
+permission to link or combine any covered work with a work licensed
+under version 3 of the GNU Affero General Public License into a single
+combined work, and to convey the resulting work.  The terms of this
+License will continue to apply to the part which is the covered work,
+but the special requirements of the GNU Affero General Public License,
+section 13, concerning interaction through a network will apply to the
+combination as such.
+
+  14. Revised Versions of this License.
+
+  The Free Software Foundation may publish revised and/or new versions of
+the GNU General Public License from time to time.  Such new versions will
+be similar in spirit to the present version, but may differ in detail to
+address new problems or concerns.
+
+  Each version is given a distinguishing version number.  If the
+Program specifies that a certain numbered version of the GNU General
+Public License "or any later version" applies to it, you have the
+option of following the terms and conditions either of that numbered
+version or of any later version published by the Free Software
+Foundation.  If the Program does not specify a version number of the
+GNU General Public License, you may choose any version ever published
+by the Free Software Foundation.
+
+  If the Program specifies that a proxy can decide which future
+versions of the GNU General Public License can be used, that proxy's
+public statement of acceptance of a version permanently authorizes you
+to choose that version for the Program.
+
+  Later license versions may give you additional or different
+permissions.  However, no additional obligations are imposed on any
+author or copyright holder as a result of your choosing to follow a
+later version.
+
+  15. Disclaimer of Warranty.
+
+  THERE IS NO WARRANTY FOR THE PROGRAM, TO THE EXTENT PERMITTED BY
+APPLICABLE LAW.  EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT
+HOLDERS AND/OR OTHER PARTIES PROVIDE THE PROGRAM "AS IS" WITHOUT WARRANTY
+OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO,
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
+PURPOSE.  THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE PROGRAM
+IS WITH YOU.  SHOULD THE PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF
+ALL NECESSARY SERVICING, REPAIR OR CORRECTION.
+
+  16. Limitation of Liability.
+
+  IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN WRITING
+WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MODIFIES AND/OR CONVEYS
+THE PROGRAM AS PERMITTED ABOVE, BE LIABLE TO YOU FOR DAMAGES, INCLUDING ANY
+GENERAL, SPECIAL, INCIDENTAL OR CONSEQUENTIAL DAMAGES ARISING OUT OF THE
+USE OR INABILITY TO USE THE PROGRAM (INCLUDING BUT NOT LIMITED TO LOSS OF
+DATA OR DATA BEING RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD
+PARTIES OR A FAILURE OF THE PROGRAM TO OPERATE WITH ANY OTHER PROGRAMS),
+EVEN IF SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF
+SUCH DAMAGES.
+
+  17. Interpretation of Sections 15 and 16.
+
+  If the disclaimer of warranty and limitation of liability provided
+above cannot be given local legal effect according to their terms,
+reviewing courts shall apply local law that most closely approximates
+an absolute waiver of all civil liability in connection with the
+Program, unless a warranty or assumption of liability accompanies a
+copy of the Program in return for a fee.
+
+                     END OF TERMS AND CONDITIONS
+
+            How to Apply These Terms to Your New Programs
+
+  If you develop a new program, and you want it to be of the greatest
+possible use to the public, the best way to achieve this is to make it
+free software which everyone can redistribute and change under these terms.
+
+  To do so, attach the following notices to the program.  It is safest
+to attach them to the start of each source file to most effectively
+state the exclusion of warranty; and each file should have at least
+the "copyright" line and a pointer to where the full notice is found.
+
+    <one line to give the program's name and a brief idea of what it does.>
+    Copyright (C) <year>  <name of author>
+
+    This program is free software: you can redistribute it and/or modify
+    it under the terms of the GNU General Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    (at your option) any later version.
+
+    This program is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU General Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+
+Also add information on how to contact you by electronic and paper mail.
+
+  If the program does terminal interaction, make it output a short
+notice like this when it starts in an interactive mode:
+
+    <program>  Copyright (C) <year>  <name of author>
+    This program comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
+    This is free software, and you are welcome to redistribute it
+    under certain conditions; type `show c' for details.
+
+The hypothetical commands `show w' and `show c' should show the appropriate
+parts of the General Public License.  Of course, your program's commands
+might be different; for a GUI interface, you would use an "about box".
+
+  You should also get your employer (if you work as a programmer) or school,
+if any, to sign a "copyright disclaimer" for the program, if necessary.
+For more information on this, and how to apply and follow the GNU GPL, see
+<http://www.gnu.org/licenses/>.
+
+  The GNU General Public License does not permit incorporating your program
+into proprietary programs.  If your program is a subroutine library, you
+may consider it more useful to permit linking proprietary applications with
+the library.  If this is what you want to do, use the GNU Lesser General
+Public License instead of this License.  But first, please read
+<http://www.gnu.org/philosophy/why-not-lgpl.html>.
diff --git a/COPYING.LESSER b/COPYING.LESSER
new file mode 100644 (file)
index 0000000..65c5ca8
--- /dev/null
@@ -0,0 +1,165 @@
+                   GNU LESSER GENERAL PUBLIC LICENSE
+                       Version 3, 29 June 2007
+
+ Copyright (C) 2007 Free Software Foundation, Inc. <http://fsf.org/>
+ Everyone is permitted to copy and distribute verbatim copies
+ of this license document, but changing it is not allowed.
+
+
+  This version of the GNU Lesser General Public License incorporates
+the terms and conditions of version 3 of the GNU General Public
+License, supplemented by the additional permissions listed below.
+
+  0. Additional Definitions.
+
+  As used herein, "this License" refers to version 3 of the GNU Lesser
+General Public License, and the "GNU GPL" refers to version 3 of the GNU
+General Public License.
+
+  "The Library" refers to a covered work governed by this License,
+other than an Application or a Combined Work as defined below.
+
+  An "Application" is any work that makes use of an interface provided
+by the Library, but which is not otherwise based on the Library.
+Defining a subclass of a class defined by the Library is deemed a mode
+of using an interface provided by the Library.
+
+  A "Combined Work" is a work produced by combining or linking an
+Application with the Library.  The particular version of the Library
+with which the Combined Work was made is also called the "Linked
+Version".
+
+  The "Minimal Corresponding Source" for a Combined Work means the
+Corresponding Source for the Combined Work, excluding any source code
+for portions of the Combined Work that, considered in isolation, are
+based on the Application, and not on the Linked Version.
+
+  The "Corresponding Application Code" for a Combined Work means the
+object code and/or source code for the Application, including any data
+and utility programs needed for reproducing the Combined Work from the
+Application, but excluding the System Libraries of the Combined Work.
+
+  1. Exception to Section 3 of the GNU GPL.
+
+  You may convey a covered work under sections 3 and 4 of this License
+without being bound by section 3 of the GNU GPL.
+
+  2. Conveying Modified Versions.
+
+  If you modify a copy of the Library, and, in your modifications, a
+facility refers to a function or data to be supplied by an Application
+that uses the facility (other than as an argument passed when the
+facility is invoked), then you may convey a copy of the modified
+version:
+
+   a) under this License, provided that you make a good faith effort to
+   ensure that, in the event an Application does not supply the
+   function or data, the facility still operates, and performs
+   whatever part of its purpose remains meaningful, or
+
+   b) under the GNU GPL, with none of the additional permissions of
+   this License applicable to that copy.
+
+  3. Object Code Incorporating Material from Library Header Files.
+
+  The object code form of an Application may incorporate material from
+a header file that is part of the Library.  You may convey such object
+code under terms of your choice, provided that, if the incorporated
+material is not limited to numerical parameters, data structure
+layouts and accessors, or small macros, inline functions and templates
+(ten or fewer lines in length), you do both of the following:
+
+   a) Give prominent notice with each copy of the object code that the
+   Library is used in it and that the Library and its use are
+   covered by this License.
+
+   b) Accompany the object code with a copy of the GNU GPL and this license
+   document.
+
+  4. Combined Works.
+
+  You may convey a Combined Work under terms of your choice that,
+taken together, effectively do not restrict modification of the
+portions of the Library contained in the Combined Work and reverse
+engineering for debugging such modifications, if you also do each of
+the following:
+
+   a) Give prominent notice with each copy of the Combined Work that
+   the Library is used in it and that the Library and its use are
+   covered by this License.
+
+   b) Accompany the Combined Work with a copy of the GNU GPL and this license
+   document.
+
+   c) For a Combined Work that displays copyright notices during
+   execution, include the copyright notice for the Library among
+   these notices, as well as a reference directing the user to the
+   copies of the GNU GPL and this license document.
+
+   d) Do one of the following:
+
+       0) Convey the Minimal Corresponding Source under the terms of this
+       License, and the Corresponding Application Code in a form
+       suitable for, and under terms that permit, the user to
+       recombine or relink the Application with a modified version of
+       the Linked Version to produce a modified Combined Work, in the
+       manner specified by section 6 of the GNU GPL for conveying
+       Corresponding Source.
+
+       1) Use a suitable shared library mechanism for linking with the
+       Library.  A suitable mechanism is one that (a) uses at run time
+       a copy of the Library already present on the user's computer
+       system, and (b) will operate properly with a modified version
+       of the Library that is interface-compatible with the Linked
+       Version.
+
+   e) Provide Installation Information, but only if you would otherwise
+   be required to provide such information under section 6 of the
+   GNU GPL, and only to the extent that such information is
+   necessary to install and execute a modified version of the
+   Combined Work produced by recombining or relinking the
+   Application with a modified version of the Linked Version. (If
+   you use option 4d0, the Installation Information must accompany
+   the Minimal Corresponding Source and Corresponding Application
+   Code. If you use option 4d1, you must provide the Installation
+   Information in the manner specified by section 6 of the GNU GPL
+   for conveying Corresponding Source.)
+
+  5. Combined Libraries.
+
+  You may place library facilities that are a work based on the
+Library side by side in a single library together with other library
+facilities that are not Applications and are not covered by this
+License, and convey such a combined library under terms of your
+choice, if you do both of the following:
+
+   a) Accompany the combined library with a copy of the same work based
+   on the Library, uncombined with any other library facilities,
+   conveyed under the terms of this License.
+
+   b) Give prominent notice with the combined library that part of it
+   is a work based on the Library, and explaining where to find the
+   accompanying uncombined form of the same work.
+
+  6. Revised Versions of the GNU Lesser General Public License.
+
+  The Free Software Foundation may publish revised and/or new versions
+of the GNU Lesser General Public License from time to time. Such new
+versions will be similar in spirit to the present version, but may
+differ in detail to address new problems or concerns.
+
+  Each version is given a distinguishing version number. If the
+Library as you received it specifies that a certain numbered version
+of the GNU Lesser General Public License "or any later version"
+applies to it, you have the option of following the terms and
+conditions either of that published version or of any later version
+published by the Free Software Foundation. If the Library as you
+received it does not specify a version number of the GNU Lesser
+General Public License, you may choose any version of the GNU Lesser
+General Public License ever published by the Free Software Foundation.
+
+  If the Library as you received it specifies that a proxy can decide
+whether future versions of the GNU Lesser General Public License shall
+apply, that proxy's public statement of acceptance of any version is
+permanent authorization for you to choose that version for the
+Library.
diff --git a/def.mk b/def.mk
new file mode 100644 (file)
index 0000000..cd01c32
--- /dev/null
+++ b/def.mk
@@ -0,0 +1,40 @@
+# Make sure I get access to all variables set when
+# installing the hlhdf package.
+include /opt/baltrad/hlhdf/mkf/hldef.mk
+
+srcdir =       .
+prefix =        /opt/rave
+
+
+CC=            gcc
+RANLIB=                ranlib
+AR=            ar
+
+OPTS=          -DNDEBUG -g -fwrapv -O3 -Wall -Wstrict-prototypes
+LDFLAGS=       
+LDSHARED=      $(CC) -shared
+CCSHARED=      -fPIC
+DEFS=          
+INCLUDE_PYTHON=        -I/opt/baltrad/third_party/include/python2.6
+SITEPACK_PYTHON= /opt/baltrad/third_party/lib/python2.6/site-packages
+
+CREATE_ITRUNC= -DCREATE_ITRUNC
+
+PROJ_INCLUDE_DIR=-I/opt/baltrad/third_party/include
+PROJ_LIB_DIR=-L/opt/baltrad/third_party/lib
+
+EXPAT_INCLUDE_DIR=
+EXPAT_LIB_DIR=
+EXPAT_SUPPRESSED=no
+
+NUMPY_INCLUDE_DIR=-I/opt/baltrad/third_party/lib/python2.6/site-packages/numpy/core/include/numpy
+
+HLHDF_ROOTDIR= /opt/baltrad/hlhdf
+HLHDF_INCLUDE_DIR= /opt/baltrad/hlhdf/include
+HLHDF_LIB_DIR= /opt/baltrad/hlhdf/lib
+HLHDF_INSTALL_BIN= /opt/baltrad/hlhdf/bin/hlinstall.sh
+HLHDF_HLDEF_MK_FILE= /opt/baltrad/hlhdf/mkf/hldef.mk
+
+RAVE_VERSION= 0.65a
+RAVE_PATCH_LEVEL=0
+
diff --git a/drain/Makefile b/drain/Makefile
new file mode 100644 (file)
index 0000000..4863d3e
--- /dev/null
@@ -0,0 +1,94 @@
+include ../def.mk
+
+CFLAGS=        $(OPTS) $(CCSHARED) $(DEFS) $(CREATE_ITRUNC) \
+       -I. $(NUMPY_INCLUDE_DIR) $(INCLUDE_PYTHON) -I$(HLHDF_INCLUDE_DIR) \
+       $(ZLIB_INCDIR) $(HDF5_INCDIR) $(PROJ_INCLUDE_DIR)
+
+INSTALL_HEADERS= 
+
+# --------------------------------------------------------------------
+# Fixed definitions
+
+IMAGE_SOURCES= image/CoordinateHandler.cpp image/FileBinary.cpp image/FileFQD.cpp image/Geometry.cpp
+RADAR_SOURCES= radar/Composite.cpp radar/Coordinates.cpp radar/Geometry.cpp radar/PolarCoordinateHandler.cpp \
+                               radar/ProductXReader.cpp radar/SubComposite.cpp
+UTIL_SOURCES= util/Data.cpp util/Debug.cpp util/Options.cpp util/Proj4.cpp util/ProjectionFrame.cpp \
+                               util/RegExp.cpp util/String.cpp
+
+IMAGE_OBJECTS= $(IMAGE_SOURCES:.cpp=.o)
+RADAR_OBJECTS= $(RADAR_SOURCES:.cpp=.o)
+UTIL_OBJECTS= $(UTIL_SOURCES:.cpp=.o)
+
+TARGET= libbdrain.so
+
+CPP=g++
+CPPFLAGS= -fPIC $(PROJ_INCLUDE_DIR) 
+
+MAKEDEPEND=g++ -MM $(CPPFLAGS) -MT '$(@D)/$(@F)' -o $(DF).d $<
+DEPDIR=.dep
+DF=$(DEPDIR)/$(@D)/$(*F)
+
+# --------------------------------------------------------------------
+# Rules
+
+# Contains dependency generation as well, so if you are not using
+# gcc, comment out everything until the $(CC) statement.
+image/%.o : image/%.cpp
+       @mkdir -p $(DEPDIR)/image; \
+       $(MAKEDEPEND); \
+       cp $(DF).d $(DF).P; \
+       sed -e 's/#.*//' -e 's/^[^:]*: *//' -e 's/ *\\$$//' \
+               -e '/^$$/ d' -e 's/$$/ :/' < $(DF).d >> $(DF).P; \
+       \rm -f $(DF).d
+       $(CPP) -c $(CPPFLAGS) $< -o $@
+
+radar/%.o : radar/%.cpp
+       @mkdir -p $(DEPDIR)/radar; \
+       $(MAKEDEPEND); \
+       cp $(DF).d $(DF).P; \
+       sed -e 's/#.*//' -e 's/^[^:]*: *//' -e 's/ *\\$$//' \
+               -e '/^$$/ d' -e 's/$$/ :/' < $(DF).d >> $(DF).P; \
+       \rm -f $(DF).d
+       $(CPP) -c $(CPPFLAGS) $< -o $@
+
+util/%.o : util/%.cpp
+       @mkdir -p $(DEPDIR)/util; \
+       $(MAKEDEPEND); \
+       cp $(DF).d $(DF).P; \
+       sed -e 's/#.*//' -e 's/^[^:]*: *//' -e 's/ *\\$$//' \
+               -e '/^$$/ d' -e 's/$$/ :/' < $(DF).d >> $(DF).P; \
+       \rm -f $(DF).d
+       $(CPP) -c $(CPPFLAGS) $< -o $@
+
+# Ensures that the .dep directory exists
+.PHONY=$(DEPDIR)
+$(DEPDIR):
+       +@[ -d $@ ] || mkdir -p $@
+
+.PHONY=all
+all:           $(TARGET)
+
+$(TARGET): $(DEPDIR) $(IMAGE_OBJECTS) $(RADAR_OBJECTS) $(UTIL_OBJECTS)
+       $(LDSHARED) -o $@ $(IMAGE_OBJECTS) $(RADAR_OBJECTS) $(UTIL_OBJECTS)
+
+.PHONY=install
+install:
+       @"$(HLHDF_INSTALL_BIN)" -f -o -C $(LIBRAVETRANSFORM) "$(prefix)/lib/$(LIBRAVETRANSFORM)"
+       @for i in $(INSTALL_HEADERS) ; \
+       do \
+               "$(HLHDF_INSTALL_BIN)" -f -o -m644 -C $$i "$(prefix)/include/$$i"; \
+       done
+
+.PHONY=clean
+clean:
+               @\rm -f image/*.o radar/*.o util/*.o core *~
+               @\rm -fr $(DEPDIR)
+
+.PHONY=distclean                
+distclean:     clean
+               @\rm -f *.so config.log config.status config.cache def.mk
+
+# NOTE! This ensures that the dependencies are setup at the right time so this should not be moved
+-include $(IMAGE_SOURCES:%.cpp=$(DEPDIR)/%.P)
+-include $(RADAR_SOURCES:%.cpp=$(DEPDIR)/%.P)
+-include $(UTIL_SOURCES:%.cpp=$(DEPDIR)/%.P)
diff --git a/drain/image/.dep/CoordinateHandler.P b/drain/image/.dep/CoordinateHandler.P
new file mode 100644 (file)
index 0000000..8318200
--- /dev/null
@@ -0,0 +1,2 @@
+CoordinateHandler.o: CoordinateHandler.cpp CoordinateHandler.h Point.h
+CoordinateHandler.cpp CoordinateHandler.h Point.h :
diff --git a/drain/image/.dep/FileBinary.P b/drain/image/.dep/FileBinary.P
new file mode 100644 (file)
index 0000000..17f6b0d
--- /dev/null
@@ -0,0 +1,6 @@
+FileBinary.o: FileBinary.cpp FileBinary.h Image.h ../util/Debug.h \
+ ../util/Options.h ../util/Data.h ../util/MapReader.h ../util/RegExp.h \
+ Geometry.h Intensity.h CoordinateHandler.h Point.h ImageView.h
+FileBinary.cpp FileBinary.h Image.h ../util/Debug.h :
+ ../util/Options.h ../util/Data.h ../util/MapReader.h ../util/RegExp.h :
+ Geometry.h Intensity.h CoordinateHandler.h Point.h ImageView.h :
diff --git a/drain/image/.dep/FileFQD.P b/drain/image/.dep/FileFQD.P
new file mode 100644 (file)
index 0000000..6e4d39e
--- /dev/null
@@ -0,0 +1,6 @@
+FileFQD.o: FileFQD.cpp FileFQD.h Image.h ../util/Debug.h \
+ ../util/Options.h ../util/Data.h ../util/MapReader.h ../util/RegExp.h \
+ Geometry.h Intensity.h CoordinateHandler.h Point.h ImageView.h
+FileFQD.cpp FileFQD.h Image.h ../util/Debug.h :
+ ../util/Options.h ../util/Data.h ../util/MapReader.h ../util/RegExp.h :
+ Geometry.h Intensity.h CoordinateHandler.h Point.h ImageView.h :
diff --git a/drain/image/.dep/Geometry.P b/drain/image/.dep/Geometry.P
new file mode 100644 (file)
index 0000000..c1872f8
--- /dev/null
@@ -0,0 +1,2 @@
+Geometry.o: Geometry.cpp Geometry.h
+Geometry.cpp Geometry.h :
diff --git a/drain/image/CatenatorOp.h b/drain/image/CatenatorOp.h
new file mode 100644 (file)
index 0000000..4a34a7e
--- /dev/null
@@ -0,0 +1,126 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef CATENATOROP_H_
+#define CATENATOROP_H_
+
+#include "ImageOp.h"
+
+namespace drain
+{
+
+namespace image
+{
+
+template <class T=unsigned char,class T2=unsigned char>
+class CatenatorOp : public ImageOp<T,T2>
+{
+public:
+       
+
+       CatenatorOp(const string & p="vert") : ImageOp<T,T2>("CatenatorOp","Catenates images, mode=vert|depth (vertically or in-depth, by adding channels). Horz not yet implemented.",
+               "mode",p)
+       {};
+
+       void filter(const Image<T> &src,Image<T2> &dst) const {
+               const string mode = this->parameters.get("mode","depth");
+               if (mode == "depth")
+                       catenateInDepth(src,dst);
+               else if (mode == "vert")
+                       catenateVertically(src,dst);
+               else
+                       throw runtime_error(mode + ": invalid option for Catenator.");
+       };
+       
+       /// The width of the image is the maximum of dst width, if src initialized, else src width.
+       void catenateInDepth(const Image<T> &src,Image<T2> &dst) const;
+       void catenateVertically(const Image<T> &src,Image<T2> &dst) const;
+};
+
+
+template <class T,class T2>
+void CatenatorOp<T,T2>::catenateInDepth(const Image<T> &src,Image<T2> &dst) const {
+       
+       const int imageChannelsSrc = src.getImageChannelCount();
+       
+       const int width  = dst.getWidth()>0 ? dst.getWidth() : src.getWidth();
+       const int height = dst.getHeight()>0 ? dst.getHeight() : src.getHeight();
+       const int imageChannelsDst = dst.getImageChannelCount();
+       const int alphaChannels    = max(src.getAlphaChannelCount(),dst.getAlphaChannelCount());
+       
+       dst.setGeometry(width,height,imageChannelsSrc + imageChannelsDst,alphaChannels);
+       
+       Point2D<int> p;
+       for (int k = 0; k < imageChannelsSrc; ++k) {
+               ImageView<T>  channelSrc(src,k);
+               ImageView<T2> channelDst(dst,imageChannelsDst + k);
+               const int w = min(width,(int)src.getWidth());
+               for (p.y = 0; p.y < height; ++p.y) {
+                       for (p.x = 0; p.x < w; ++p.x) {
+                               channelDst.at(p) = channelSrc.at(p);
+                       }
+               }
+       }
+       //cout << "GEOM NOW " << dst.getGeometry() << '\n';
+       
+};
+
+template <class T,class T2>
+void CatenatorOp<T,T2>::catenateVertically(const Image<T> &src,Image<T2> &dst) const {
+       
+       
+       const int width  = dst.getWidth()>0 ? dst.getWidth() : src.getWidth();
+       const int heightDst = dst.getHeight();
+       const int channelsDst = dst.getChannelCount();
+       
+       const int heightSrc = src.getHeight();
+       const int channelsSrc = src.getChannelCount();
+       //const int alphaChannels    = max(src.getAlphaChannelCount(),dst.getAlphaChannelCount());
+       
+       const int start = width*heightDst*channelsDst;
+       
+       dst.setGeometry(width,heightDst*channelsDst+heightSrc*channelsSrc,1,0);
+       
+       
+       //Point2D<int> p;
+       for (int k = 0; k < channelsSrc; ++k) {
+               //BufferedView<T>  channelSrc(src,k);
+               //BufferedView<T2> channelDst(dst,channelsDst + k);
+               const int w = min(width,(int)src.getWidth());
+               //for (p.y = 0; p.y < height; ++p.y) {
+               //      for (p.x = 0; p.x < w; ++p.x) {
+               for (int j=0; j<heightSrc; j++) {
+                       for (int i=0; i<w; i++) {
+                               //dst.at(start+src.address(i,j,k)) = src.at(i,j,k);
+                               dst.at(start+src.address(i,j,k)) = static_cast<T2>(src.at(i,j));
+                               //dst.at(i,j,k) = src.at(i,j,k);
+                       }
+               }
+       }
+       //cout << "GEOM NOW " << dst.getGeometry() << '\n';
+       
+};
+
+}
+
+}
+
+#endif /*CATENATOR_H_*/
diff --git a/drain/image/CoordinateHandler.cpp b/drain/image/CoordinateHandler.cpp
new file mode 100644 (file)
index 0000000..bfdbb6a
--- /dev/null
@@ -0,0 +1,185 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#include "CoordinateHandler.h"
+
+namespace drain
+{
+
+namespace image
+{
+
+CoordinateHandler::CoordinateHandler(const int &width, const int &height){
+       setBounds(width,height);
+       status = UNCHANGED;
+       name = "CoordinateHandler";
+}
+
+
+int CoordinateHandler::handle(int &x, int &y) const 
+{
+
+    status = UNCHANGED;
+
+    if (x < xMin)
+    {
+        x = xMin;
+        status = X_UNDERFLOW | IRREVERSIBLE;
+    }
+    else if (x > xMax)
+    {
+        x = xMax;
+        status = X_OVERFLOW | IRREVERSIBLE;
+    }
+
+    if (y < yMin)
+    {
+        y = yMin;
+        status |= Y_UNDERFLOW | IRREVERSIBLE;
+    }
+    else if (y > yMax)
+    {
+        y = yMax;
+        status |= Y_OVERFLOW | IRREVERSIBLE;
+    }
+
+    return status;
+}
+
+Mirror::Mirror(){
+        name = "Mirror";
+};
+
+int Mirror::handle(int &x, int &y) const
+{
+    status = UNCHANGED;
+
+    if (x < xMin)
+    {
+        status = X_UNDERFLOW;
+        x -= xMin;
+        x = (x%doubleWidth);
+
+        if (x < 0)
+            x += doubleWidth;
+        if (x >= width)
+            x = doubleWidth-x-1;
+        x += xMin;
+    }
+    else if (x > xMax)
+    {
+        status = X_OVERFLOW;
+        x -= xMin;
+        x = (x%doubleWidth);
+        if (x < 0)
+            x += doubleWidth;
+        if (x >= width)
+            x = doubleWidth-x-1;
+        x += xMin;
+    }
+
+    if (y < yMin)
+    {
+        status |= Y_UNDERFLOW;
+        y -= yMin;
+        y = (y%doubleHeight);
+        if (y < 0)
+            y += doubleHeight;
+        if (y >= height)
+            y = doubleHeight-y-1;
+        y += yMin;
+    }
+    else if (y > yMax)
+    {
+        status |= Y_OVERFLOW;
+        y -= yMin;
+        y = (y%doubleHeight);
+        if (y < 0)
+            y += doubleHeight;
+        if (y >= height)
+            y = doubleHeight-y-1;
+        y += yMin;
+    }
+
+    return status;
+}
+
+
+Wrapper::Wrapper(){
+       name = "Wrapper";
+};
+
+int Wrapper::handle(int &x, int &y) const
+{
+    status = UNCHANGED;
+    if (x < xMin)
+    {
+        status = X_UNDERFLOW;
+        x -= xMin;
+        x %= width;
+        if (x < 0)
+            x += width;
+        x += xMin;
+    }
+    else if (x > xMax)
+    {
+        status = X_OVERFLOW;
+        x -= xMin;
+        x %= width;
+        if (x < 0)
+            x += width;
+        x += xMin;
+    }
+
+    if (y < yMin)
+    {
+        status |= Y_UNDERFLOW;
+        y -= yMin;
+        y %= height;
+        if (y < 0)
+            y += height;
+        y += yMin;
+    }
+    else if (y > yMax)
+    {
+        status |= Y_OVERFLOW;
+        y -= yMin;
+        y %= height;
+        if (y < 0)
+            y += height;
+        y += yMin;
+    }
+    return status;
+}
+
+
+ostream & operator<<(ostream &ostr,const CoordinateHandler &handler){
+       ostr << handler.name << '[' << handler.xMin << ',' << handler.yMin << ',';
+       ostr << handler.xMax << ',' << handler.yMax << ']' << handler.getStatus();
+       return ostr;
+}
+
+}
+}
+
+
+
+
diff --git a/drain/image/CoordinateHandler.h b/drain/image/CoordinateHandler.h
new file mode 100644 (file)
index 0000000..d873e8f
--- /dev/null
@@ -0,0 +1,152 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef COORDINATEHANDLER_H_
+#define COORDINATEHANDLER_H_
+
+
+#include "Point.h"
+
+namespace drain
+{
+
+namespace image
+{
+
+// Facility for handling coordinate under- and overflows.
+/*!  
+ * This base implementation forces the point (i,j) to be inside the image.
+ * Any under- or overflow is tagged IRREVERSIBLE.
+ * 
+ */
+class CoordinateHandler
+{
+public:
+
+
+    CoordinateHandler(const int &width=0, const int &height=0);
+    virtual ~CoordinateHandler(){};
+    
+    /// No change in coordinates has taken place.
+    static const int UNCHANGED = 0;
+    static const int X_OVERFLOW = 1;
+    static const int X_UNDERFLOW = 2;
+    static const int Y_OVERFLOW = 4;
+    static const int Y_UNDERFLOW = 8;
+
+    /** Coordinates have been changed such that an inverse move would not return to original position.
+     */
+    static const int IRREVERSIBLE = 128;
+
+    //public Rectangle bounds;
+    int xMin;
+    /** Typically imageWidth-1.
+     */
+    int xMax;
+    int yMin;
+    /** Typically imageHeight-1.
+     */
+    int yMax;
+
+
+       
+
+    /** Sets the size of the image (or any rectangular area) to be accessed.
+     *  TODO: rename to setGeometry ?
+     * 
+     * @param width
+     * @param height
+     */
+    void setBounds(int width,int height)
+    {
+        this->width = width;
+        this->height = height;
+
+        doubleWidth  = 2*width;
+        doubleHeight = 2*height;
+
+        xMin = 0;
+        xMax = width-1;
+        yMin = 0;
+        yMax = height-1;
+
+        area = width * height;
+        // cerr << "CoordinateHandler: bounds set to: " << width << "," << height << endl;
+    }
+
+
+    /*!
+     * 
+     * 
+     */
+    virtual int handle(int &x, int &y) const;
+       
+       /// Returns status (zero if no edges crossed).
+       inline 
+       int handle(Point2D<> & p) const { return handle(p.x,p.y);       };
+       string name;
+       
+       int getStatus() const{ return status; };
+       
+ protected:
+       /** Width of the image(s) to be accessed.
+        */
+       int width;
+       int doubleWidth;
+
+       /** Height of the image(s) to be accessed.
+        */
+       int height;
+       int doubleHeight;
+
+       int area;
+
+       mutable int status; // = UNCHANGED;
+
+};
+
+
+/// Folds the coordinates back to the image.
+class Mirror :  public CoordinateHandler
+{
+public:
+       Mirror(); 
+       virtual int handle(int &x, int &y) const;
+    //int handle(Point2D<> &p);
+};
+
+/// Repeats the image in the infinities.
+class Wrapper :  public CoordinateHandler
+{
+public:
+       Wrapper();
+       virtual int handle(int &x, int &y) const;
+    //int handle(Point2D<> &p);
+};
+    
+ostream & operator<<(ostream &ostr,const CoordinateHandler &handler);
+
+}
+}
+
+#endif /*COORDINATEHANDLER_H_*/
diff --git a/drain/image/CopyOp.h b/drain/image/CopyOp.h
new file mode 100644 (file)
index 0000000..37a12d8
--- /dev/null
@@ -0,0 +1,134 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef COPYOP_H_
+#define COPYOP_H_
+
+#include "SequentialImageOp.h"
+
+namespace drain
+{
+
+namespace image
+{
+
+/** Copies intensities. Does not scale them. Uses
+ *  Internally, uses SequentialImageOp.
+ *  @see ScaleOp
+ */    
+template <class T=unsigned char,class T2=unsigned char>
+class CopyOp : public ImageOp<T,T2>
+{
+public:
+       CopyOp(const string & p = "f,f"){
+               this->setInfo("Copies channels. f=full image, i=image channels, a=alpha channel(s), 0=1st, 1=2nd,...",
+               "srcView,dstView",p);
+       };
+
+       virtual void filter(const Image<T> &src,Image<T2> &dst) const {
+               
+               ImageView<T>  srcView;
+               //const char s = this->parameters.get("srcView",'f');
+               const Data &s = this->getParameter("srcView","f");
+               switch (s[0]) {
+               case 'f': // full image
+                               srcView.viewImage(src);
+                               break;
+                       case 'i':  // image channels (ie. excluding alpha)
+                               srcView.viewChannels(src,src.getImageChannelCount());
+                               break;
+                       case 'a': // alpha channel(s)
+                               //srcView.viewImage(src.getAlphaChannel());
+                               srcView.viewChannels(src,src.getAlphaChannelCount(),src.getImageChannelCount());
+                               break;
+                       default:
+                               unsigned int k = s;
+                               //if ((k==0)&&(s[0]!=)
+                               //if (!s.isType<int>()){
+                               //cerr << "Unsupported channel option:" << this->parameters.get("srcView");
+                               //throw runtime_error("CopyOp: illegal source channel code.");
+                               //};
+                               srcView.viewChannel(src,k);
+               }
+               
+               //cerr << src.getGeometry() << '\n';
+               //cerr << srcView.getGeometry() << '\n';
+               
+               
+               makeCompatible(srcView,dst);
+               //cerr << dst.getGeometry() << " own="<< dst.hasOwnBuffer() <<'\n';
+               
+               
+               ImageView<T2> dstView;
+               const Data &d = this->getParameter("dstView","f");
+               
+               switch (d[0]) {
+               case 'f': // full image
+                       dstView.viewImage(dst);
+                               break;
+                       case 'i':  // image channels (ie. excluding alpha)//makeCompatible(srcView,dst);
+                               dstView.viewChannels(dst,dst.getImageChannelCount());
+                               break;
+                       case 'a': // alpha
+                               //cerr << "alpha\n";
+                               if (dst.getAlphaChannelCount() == 0)
+                                       dst.setAlphaChannelCount(srcView.getChannelCount());
+                               //cerr << "alpha\n";
+                               dstView.viewChannels(dst,dst.getAlphaChannelCount(),dst.getImageChannelCount());
+                               break;
+                       default:
+                               //int k = s;
+                               /*
+                               if (!d.isType<int>()){
+                                       cerr << "Unsupported channel option:" << d;
+                                       throw runtime_error("CopyOp: illegal source channel code.");    
+                               };
+                               */
+                               unsigned int kd = d;
+                               if (dst.getChannelCount() <= kd)
+                                       dst.setChannelCount(kd+1);
+                               dstView.viewChannel(dst,kd);
+                               break;
+               }
+
+               if (Debug > 0){
+                       cerr << srcView.getGeometry() << '\n';
+                       cerr << dstView.getGeometry() << '\n';
+               }
+
+               typename std::vector<T>::const_iterator si = srcView.begin();
+        typename std::vector<T2>::iterator di;
+        typename std::vector<T2>::const_iterator dEnd = dstView.end();
+               for (di=dstView.begin(); di!=dEnd; di++,si++)
+                       *di = static_cast<T2>(*si);
+
+
+
+       };
+
+       //inline        virtual void filterValue(const T &src, T2 &dst){ dst = static_cast<T2>(src);};
+};
+
+}
+
+}
+
+#endif /*COPYOP_H_*/
diff --git a/drain/image/Cumulator.h b/drain/image/Cumulator.h
new file mode 100644 (file)
index 0000000..3ffa9bd
--- /dev/null
@@ -0,0 +1,456 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef CUMULATOR_H_
+#define CUMULATOR_H_
+
+//#include "coordinates.h" // for site coords and bin coords.
+
+
+//#include <drain/util/proj.h>  // for geographical projection of radar data bins
+
+#include <math.h>
+
+#include <stdexcept>
+
+#include "Point.h"
+#include "Image.h"
+#include "CoordinateHandler.h" 
+
+
+
+// TODO: image/
+/** See also radar::Compositor
+ * 
+ */
+namespace drain
+{
+
+namespace image
+{
+
+
+/// General-purpose image compositing.
+/*!
+  Main features:
+
+  Injection. The main loop iterates radar (image) data points of type <Ti>, 
+  mapping each point to a single point in the target array of type <T>;
+  hence the name injector.
+
+  Cumulation. The compositing starts by initiating the target image,
+  see open().
+  Each input radar image is projected on the target image cumulatively,
+  see execute(). Finally, the cumulation must be completed by
+  normalizing the result with the number/weight of
+  cumulated inputs, see close().
+  Hence, the target which should be seen.
+
+  The weights are stored in the first alphaChannel (index 0) of the
+  target image. If the target has a second image channel, also
+  the (weighted) squared values are cumulated.  
+
+  In cumulation, ...
+  \f[
+  x = \sqrt[p]{\frac{\sum_i {q_i}^r x_i^p}{\sum_i {q_i}^r } }
+  \f]
+
+  The default type of target image, that is, cumulation array, is
+  double. If no (float valued) weighting  is applied, target image
+  could be long int as well.
+ */
+template <class T = double>
+class Cumulator 
+{
+public:
+
+       /// Default constructor. The channels are DATA, COUNT, WEIGHT, WEIGHT2
+       Cumulator(unsigned int width = 0, unsigned int height = 0, const string & method="WAVG", double p=1.0,double r=0.0){
+               setGeometry(width,height);
+               setMethod(method,p,r);
+       }
+
+       virtual ~Cumulator(){};
+
+       virtual
+       void setGeometry(unsigned int width,unsigned int height){
+               this->width = width;
+               this->height = height;
+
+               data.setGeometry(width,height);
+               weight.setGeometry(width,height);
+               count.setGeometry(width,height);
+               dataSquared.setGeometry(width,height);
+
+               coordinateHandler.setBounds(width,height);
+       }
+
+       void clear(){
+               data.clear();
+               weight.clear();
+               count.clear();
+               dataSquared.clear();
+       }
+
+       enum Method { AVG, WAVG, MAX, MAXW, OVERWRITE, UNDEFINED};
+
+       Method method;
+
+       void setMethod(const string & method, double p=1,double r=0){
+               if (method=="AVG")
+                       setMethod(AVG);
+               else if (method=="WAVG")
+                       setMethod(WAVG,p,r);
+               else if (method=="MAX")
+                       setMethod(MAX);
+               else if (method=="MAXW")
+                       setMethod(MAXW);
+               else if (method=="OVERWRITE")
+                       setMethod(OVERWRITE);
+               else
+                       setMethod(UNDEFINED);
+       }
+
+       void setMethod(Method m, double p=1,double r=0){
+                       this->method = m;
+                       this->p = p;
+                       this->r = r;    
+                       if (m == UNDEFINED)
+                               throw runtime_error(method + ": undefined method");
+       }
+
+       inline
+       const unsigned int & getWidth() const {return width;};
+
+       inline
+       const unsigned int & getHeight() const {return height;};
+
+       //void add(const Image &src,int x=0,int y=0);
+
+       //template <class T>
+       ///
+       inline
+       void add(const int &i, const int &j,double value, double weight);
+
+       //template <class T>    void putWeighted(int i, int j,int k,T value, T weight);
+       /// \deprecated Getting obsolete.ss
+       inline
+       //void addDirect(const T &s,const T &w,const int &i=0, const int &j=0);
+       void addDirect(const int &i, const int &j,const T &s,const T &w,const unsigned int &c,const T &s2);
+
+       //void add(const Image &src,int x=0,int y=0);
+
+       /// Copies the contents of the cumulation array to the target image.
+       // COnsider template based such that INTERNAL variables are of type T ? 
+       /*! d = data, scaled
+        *  w = weight, scaled
+        *  c = count
+        *  p = data, sum of squared ("power"), scaled
+        *  s = standard deviation, scaled
+        *  D = data, cumulated  (debugging)
+        *  W = weight, cumulated
+        *  C = count, scaled
+        *  S = standard deviation, unscaled
+        */
+       //void extract(Image &target,string channels = "dw") const;
+       template <class T2>
+       void extractTo(Image<T2> &dst,const string & channels="dw") const;
+
+       /// Drops an image
+       template <class T2>
+       void addImage(const Image<T2> &src,const Image<T2> &srcWeight,
+                       float globalWeight = 1.0, int i=0, int j=0);
+
+       // TODO: different projections
+
+
+
+       /// Overall weight for source image.
+       //double priorWeight;
+
+       /// Data weight (power)
+       double p;
+
+       /// Quality weight (power)
+       double r;
+
+       bool debug;
+
+
+       CoordinateHandler coordinateHandler;
+
+protected:
+       unsigned int width;
+       unsigned int height;
+
+       //private:
+       Image<T> data;
+       Image<T> weight;
+       Image<unsigned int> count;
+       Image<T> dataSquared;
+
+       Options properties;
+};
+
+
+template <class T>
+void Cumulator<T>::add(const int &i, const int &j, double s, double w){
+
+       //if (coordinateHandler.handle(int(i)));
+       if (i<0)
+               return;
+
+       if (j<0)
+               return;
+
+       if (i>=static_cast<int>(width))
+               return;
+
+       if (j>=static_cast<int>(height))
+               return;
+
+
+       const long int it = data.address(i,j);
+
+       static const T wd = 128;
+       //if (r != 1)
+       //  w = pow(w,r);
+       //cerr << (float)s << '\t';
+
+       switch (method) {
+       case AVG:
+               if (w > 0){ // Skip no-data areas
+                       data.at(it)   += wd*s;
+                       weight.at(it) += wd;
+                       count.at(it)  += 1;
+                       dataSquared.at(it) += wd*s*s;
+               }
+               break;
+       case WAVG:
+               if (p != 1)
+                       s = pow(s,p);
+               if ((r > 0.0)&&(r != 1.0))
+                       w = pow(w,r);
+               data.at(it)   += w * s;
+               weight.at(it) += w;
+               count.at(it)  += 1;
+               dataSquared.at(it) += w * s*s;
+               break;          
+       case MAX:
+               if ((weight.at(it)==0)||(s >= (data.at(it)/weight.at(it)))){
+                       data.at(it) = w * s;
+                       weight.at(it) = w;
+                       count.at(it)  = 1;
+               }
+               break;          
+       case MAXW:
+               if (w >= weight.at(it)){
+                       data.at(it) = w * s;
+                       weight.at(it) = w;
+                       count.at(it)  = 1;
+               }
+               break;
+       case OVERWRITE:
+               data.at(it) = w*s;
+               weight.at(it) = w;
+               count.at(it)  = 1;
+               dataSquared.at(it) = w*s*s;
+       //case UNDEFINED:
+       //      break;
+       default:
+               throw runtime_error(method + ": undefined method");
+       }
+
+       //dataSquared.at(it) = 57;
+
+}
+
+template <class T>
+void Cumulator<T>::addDirect(const int &i, const int &j,const T &s,const T &w,const unsigned int &c,const T &s2){
+
+       if ((i<0)||(j<0)||(i>=width)||(j>=height))
+               return;
+
+       const long int it = data.address(i,j);
+       data.at(it)   += s;
+       weight.at(it) += w;
+       count.at(it)  += c;
+       dataSquared.at(it) += s2;
+}
+
+template <class T>
+template <class T2>
+void Cumulator<T>::addImage(const Image<T2> &src,const Image<T2> &srcWeight,
+               float globalWeight, int iOffset, int jOffset){
+
+       const unsigned int w = src.getWidth();
+       const unsigned int h = src.getHeight();
+       //const unsigned int wCum = getWidth();
+       //const unsigned int hCum = getHeight();
+
+       int a;
+       Point2D<int> p;
+       //CoordinateHandler handler(getWidth(),getHeight());
+       //const int r = 64+(rand()&127);
+       for (unsigned int i = 0; i < w; ++i) {
+               for (unsigned int j = 0; j < h; ++j) {
+                       p.setLocation(iOffset+i,jOffset+j);
+                       if (coordinateHandler.handle(p) == CoordinateHandler::UNCHANGED){
+                               a = src.address(i,j);
+                               add(p.x,p.y,src.at(a),globalWeight * srcWeight.at(a));
+                       }
+                       // else overflow=true;
+
+                       //add(iOffset+i,jOffset+j,src.at(i,j),100);
+                       //add(i,j,src.at(i,j),globalWeight*10.0);
+                       //add(i,j,r,globalWeight*10.0);
+                       //debug
+                       //dataSquared.at(i,j) = (i+j)&255;
+                       //data.at(i,j) = 123;
+                       //weight.at(i,j) = 246;
+               }
+       }
+}
+
+
+template <class T>
+template <class T2>
+void Cumulator<T>::extractTo(Image<T2> &dst,const string &channels) const {
+
+
+       // const Geometry g = getGeometry();
+       // const int width  = g.getWidth();
+       // const int height = g.getHeight();
+       const int area = width * height;
+       const unsigned int channelCount = channels.length();
+       //const int area = g.getArea();
+       //const int alphaChannels =  g.getAlphaChannelCount();
+       unsigned int dstImageChannelCount = channelCount;
+
+       if ((dst.getWidth() != width) || (dst.getHeight() != height) ||
+                       (dst.getChannelCount() != channelCount)){
+               dst.setGeometry(width,height,channelCount,0);
+               // cerr warning...
+       }
+
+
+       if (drain::Debug > 2){
+               cerr << " Extracting cumulated image To dst:\n";
+               dst.debug(cerr);
+       }
+
+       double w = 0;
+       double c = 0;
+
+       for (unsigned int k=0; k<channelCount; ++k) {
+
+               char ch = channels.at(k);
+
+               if (drain::Debug > 3)
+                       cerr << "Extracting channel: " << ch << endl;
+
+               Image<T2> &channel = dst.getChannel(k);
+
+               if ((ch != 'D') && (ch != 'd') && (k<dstImageChannelCount))
+                       dstImageChannelCount = k;
+
+               switch (ch) {
+               case 'D':  // data, unscaled
+                       for (int i = 0; i < area; ++i) {
+                               channel.at(i) = static_cast<T2>(data.at(i));
+                       }
+                       break;
+               case 'd':  // data, scaled
+                       for (int i = 0; i < area; ++i) {
+                               w = weight.at(i);
+                               if (w > 0){  // or w!=0 ?
+                                       if ((p == 1)||(p == 0)) // unity??
+                                               channel.at(i) = static_cast<T2>(data.at(i)/w);
+                                       else //if (p > 0)
+                                               channel.at(i) = static_cast<T2>(pow(data.at(i)/w,1.0/p));
+                               }
+                               else
+                                       channel.at(i) = 0;
+                       }
+                       break;
+               case 'W':  // weight, scaled (average weight)
+                       for (int i = 0; i < area; ++i) {
+                               channel.at(i) = static_cast<T2>(weight.at(i));
+                       }
+                       break;
+               case 'w':  // weight, scaled (average weight)
+                       for (int i = 0; i < area; ++i) {
+                               c = count.at(i);
+                               if (c > 0){
+                                       if ((r == 1)||(r == 0)) // unity?
+                                               channel.at(i) = static_cast<T2>(weight.at(i)/c);
+                                       else //if (r>0)
+                                               channel.at(i) = static_cast<T2>(pow(weight.at(i)/c,1.0/r));
+                               }
+                               else
+                                       channel.at(i) = 0;
+                       }
+                       break;
+               case 'C':  // count
+                       for (int i = 0; i < area; ++i) {
+                               channel.at(i) = static_cast<T2>(count.at(i));
+                       }
+                       break;
+               case 'c':  // count, scaled
+                       for (int i = 0; i < area; ++i) {
+                               channel.at(i) = 255 - 255/(1+static_cast<T2>(count.at(i)));
+                       }
+                       break;
+               case 'S':  // st.deviation, unscaled (and quite good so) // TODO other is count-based?
+                       for (int i = 0; i < area; ++i) {
+                               double d  = data.at(i);
+                               double d2 = dataSquared.at(i);
+                               double w  = weight.at(i);
+                               //unsigned int c = count.at(i);
+                               if (w > 0){
+                                       //channel.at(i) = static_cast<T2>(d2);
+                                       channel.at(i) = static_cast<T2>(sqrt((double)(d2/w - (d*d)/(w*w))));
+                               }
+                               else
+                                       channel.at(i) = 0;
+                       }
+                       break;
+               case 's':  // count, scaled
+               default:
+                       throw runtime_error(string("Error: Cumulator: undefined channel extrator: ") +  ch);
+                       break;
+               }
+       }
+       dst.setChannelCount(dstImageChannelCount,channelCount-dstImageChannelCount);
+       //cerr << "setting dst geometry" << dst.getGeometry() << '\n';
+       dst.properties = properties;
+}
+
+
+
+
+
+
+}
+
+}
+
+#endif /* Cumulator_H_ */
diff --git a/drain/image/DistanceTransformFillOp.h b/drain/image/DistanceTransformFillOp.h
new file mode 100644 (file)
index 0000000..8bf5665
--- /dev/null
@@ -0,0 +1,314 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef DISTANCETRANSFORMFILLOP_H_
+#define DISTANCETRANSFORMFILLOP_H_
+       
+       //#include "BufferedImage<T>.h"
+       //#include "Operator.h"
+       //#include "Point.h"
+       
+#include "DistanceTransformOp.h"
+//#include "File.h"
+       
+       
+namespace drain
+{
+       
+namespace image
+{
+       
+       /**
+        *  Semantics. In filtering, if there are as many alpha channels as image channels, the channels
+        *  are treated separately. Else, the image channels are treated as colors, and the only
+        *  alpha channel (or the last of many) will serve as propagation criterion.
+        */
+       template <class T=unsigned char,class T2=unsigned char>
+       class DistanceTransformFillOp : public DistanceTransformOp<T,T2> 
+       {
+       public:
+               
+                       
+           DistanceTransformFillOp(const string & p = "5,-1,-1") : 
+               DistanceTransformOp<T,T2> (this->distanceLinear) {
+               this->setInfo("Spreads intensities to ranges defined by alpha intensities.",
+               "horz,vert,diag",p);
+               };
+               
+               /// For derived classes.
+               DistanceTransformFillOp(DistanceModel<T2> & distanceModel) : 
+               DistanceTransformOp<T,T2>(distanceModel) {
+           }; 
+           
+               ~DistanceTransformFillOp(){};
+       
+               
+               void filter(const Image<T> &src, Image<T2> &dest) const;
+               void filter(const Image<T> &src, const Image<T> &srcAlpha, 
+                       Image<T2> &dst, Image<T2> &dest) const;
+       
+               void filterDownRight(const Image<T> &src,  const Image<T> &srcAlpha,
+                       Image<T2> &dst, Image<T2> &destAlpha) const;
+               
+               void filterUpLeft(const Image<T> &src,  const Image<T> &srcAlpha,
+                       Image<T2> &dst, Image<T2> &destAlpha) const ;
+       
+       };
+               
+template <class T,class T2>
+void DistanceTransformFillOp<T,T2>::filter(
+       const Image<T> &src, Image<T2> &dst) const {
+               if (drain::Debug > 2){
+                       cerr << "DistanceTransformFillOp<T,T2>::filter\n" << this->getParameters() << '\n';
+               }
+               makeCompatible(src,dst);
+               const unsigned int iCh = src.getImageChannelCount();
+
+               const unsigned int aCh = src.getAlphaChannelCount();
+               if (aCh == iCh){
+                       for (unsigned int i = 0; i < iCh; ++i)
+                               filter(src.getChannel(i), src.getChannel(iCh+i), dst.getChannel(i), dst.getChannel(iCh+i));
+               }
+               else {
+                       const ImageView<T> srcView(src,0,iCh);
+                       ImageView<T2> dstView(dst,0,iCh);
+                       filter((const Image<T>&)srcView, src.getAlphaChannel(), (Image<T2> &)dstView, dst.getAlphaChannel());
+               }
+};
+               
+template <class T,class T2>
+void DistanceTransformFillOp<T,T2>::filter(
+               const Image<T> &src, const Image<T> &srcAlpha, 
+                       Image<T2> &dst, Image<T2> &dstAlpha) const {
+
+                       // TODO
+                       //if (filterWithTmp(src,srcAlpha,dst,dstAlpha)) return;
+                       
+                       dst.name = "DTF target";
+                       
+                       makeCompatible(src,dst);
+                       makeCompatible(srcAlpha,dstAlpha);
+                       
+                       this->initialize();
+                       
+                       filterDownRight(src,srcAlpha,dst,dstAlpha);
+                       //drain::image::File::write(dst,"dtd-i.jpg");
+                       //drain::image::File::write(dstAlpha,"dtd-i.png");
+                       
+                       filterUpLeft(dst,dstAlpha,dst,dstAlpha);
+                       
+                       //drain::image::File::write(dst,"dtu-i.jpg");
+                       //drain::image::File::write(dstAlpha,"dtu-i.png");
+                       
+               }
+               
+               
+template <class T,class T2>
+void DistanceTransformFillOp<T,T2>::filterDownRight(const Image<T> &src,const Image<T> &srcAlpha, 
+       Image<T2> &dst, Image<T2> &dstAlpha) const {
+       
+               
+                       const int width  = src.getWidth();
+                       const int height = src.getHeight();
+                       const CoordinateHandler & coordinateHandler = src.getCoordinateHandler();
+                       
+                       unsigned int iChannels = src.getImageChannelCount();
+       
+                       // proximity (inverted distance)
+                       float d;
+                   float dPrev;
+                               
+                       // Intensity (graylevel)
+                       //T g;
+               //cerr << "iChannels: " << iChannels << endl;
+               vector<T2> pixel(iChannels);
+       
+               //cerr << "koe" << endl;
+                       
+                       Point2D<> p; 
+                       
+                       Point2D<> t;
+                       int &tx = t.x;
+                       int &ty = t.y;
+        
+                       //coordinateOverflowHandler.setBounds(srcDist.getBounds());
+                       
+                       for (ty=0; ty<height; ty++){
+                               for (tx=0; tx<width; tx++){
+                                       
+                                       // Take source value as default
+                                       d = srcAlpha.at(t);
+                                       src.getPixel(t,pixel); 
+                                       
+                                       // Compare to previous value
+                                       dPrev = dstAlpha.at(t);
+                                       if (dPrev > d){
+                                               d = dPrev;
+                                               dst.getPixel(t,pixel);                                  
+                                       }
+                                       
+                                       // Compare to upper left neighbour       
+                                       p.setLocation(tx-1,ty-1);
+                                       coordinateHandler.handle(p);
+                                       dPrev = this->distanceModel.decreaseDiag(dstAlpha.at(p));
+                                       if (dPrev > d){
+                                               d = dPrev;
+                                               dst.getPixel(p,pixel);                          
+                                       } 
+       
+                                       // Compare to upper neighbour 
+                                       p.setLocation(tx,ty-1);
+                                       coordinateHandler.handle(p);
+                                       dPrev = this->distanceModel.decreaseVert(dstAlpha.at(p));
+                                       if (dPrev > d){
+                                               d = dPrev;
+                                               dst.getPixel(p,pixel);                  
+                                       } 
+       
+                                       // Compare to upper right neighbour 
+                                       p.setLocation(tx+1,ty-1);
+                                       coordinateHandler.handle(p);
+                                       dPrev = this->distanceModel.decreaseDiag(dstAlpha.at(p));
+                                       if (dPrev > d){
+                                               d = dPrev;
+                                               dst.getPixel(p,pixel);                                  
+                                       } 
+       
+                                       // Compare to left neighbour 
+                                       p.setLocation(tx-1,ty);
+                                       coordinateHandler.handle(p);
+                                       dPrev = this->distanceModel.decreaseHorz(dstAlpha.at(p));
+                                       if (dPrev > d){
+                                               d = dPrev;
+                                               dst.getPixel(p,pixel);                                  
+                                       } 
+       
+                                       if (d>0){
+                                               dstAlpha.at(t) = static_cast<T2>(d);
+                                               dst.putPixel(t, pixel);
+                                       }
+       
+                               }
+                       }
+                       // return dst;
+               };
+               
+template <class T,class T2>
+void DistanceTransformFillOp<T,T2>::filterUpLeft(
+       const Image<T> &src, const Image<T> &srcAlpha,
+                       Image<T2> &dst, Image<T2> &dstAlpha) const {
+       
+                       const int width  = src.getWidth();
+                       const int height = src.getHeight();
+                       const CoordinateHandler & coordinateHandler = src.getCoordinateHandler();
+                       
+                       // proximity (inverted distance)
+                       float d;
+                       float dPrev;
+                       
+                       vector<T> pixel(src.getImageChannelCount());
+                               
+                       Point2D<> p;
+                       
+                       Point2D<> t;
+                       int &tx = t.x;
+                       int &ty = t.y;
+        
+       
+                       //coordinateHandler.setBounds(width,height);
+                       for (ty=height-1; ty>=0; ty--){
+                               for (tx=width-1; tx>=0; tx--){
+                                       
+                                       // Source
+                                       d = srcAlpha.at(t); 
+                                       src.getPixel(t,pixel);
+                                       
+                                       // Compare to previous value
+                                       dPrev = dstAlpha.at(t);
+                                       if (dPrev > d){
+                                               d = dPrev;
+                                               dst.getPixel(t,pixel);
+                                       }
+                                       
+                                       // Compare to lower left neighbour       
+                                       p.setLocation(tx-1,ty+1);
+                                       coordinateHandler.handle(p);
+                                       dPrev = this->distanceModel.decreaseDiag(dstAlpha.at(p));
+                                       if (dPrev > d){
+                                               d = dPrev;
+                                               dst.getPixel(p,pixel);                                  
+                                       } 
+       
+                                       // Compare to lower neighbour 
+                                       p.setLocation(tx,ty+1);
+                                       coordinateHandler.handle(p);
+                                       dPrev = this->distanceModel.decreaseVert(dstAlpha.at(p));
+                                       if (dPrev > d){
+                                               d = dPrev;
+                                               dst.getPixel(p, pixel);                                 
+                                       } 
+       
+                                       // Compare to lower right neighbour 
+                                       p.setLocation(tx+1,ty+1);
+                                       coordinateHandler.handle(p);
+                                       dPrev = this->distanceModel.decreaseDiag(dstAlpha.at(p));
+                                       if (dPrev > d){
+                                               d = dPrev;
+                                               dst.getPixel(p,pixel);                                  
+                                       } 
+       
+                                       // Compare to right neighbour 
+                                       p.setLocation(tx+1,ty);
+                                       coordinateHandler.handle(p);
+                                       dPrev = this->distanceModel.decreaseHorz(dstAlpha.at(p));
+                                       if (dPrev > d){
+                                               d = dPrev;
+                                               dst.getPixel(p,pixel);                                  
+                                       } 
+       
+                                       if (d>0){
+                                               dstAlpha.at(t) = static_cast<T2>(d);
+                                               dst.putPixel(t, pixel);
+                                       }
+       
+       
+                       }
+               }
+       };
+
+       template <class T=unsigned char,class T2=unsigned char>
+       class DistanceTransformFillExponentialOp : public DistanceTransformFillOp<T,T2> 
+       {
+       public:
+                
+           DistanceTransformFillExponentialOp(const string & p = "5,-1,-1") : 
+               DistanceTransformFillOp<T,T2> (this->distanceExponential) {
+               this->setInfo("Spreads intensities to ranges defined by alpha intensities.",
+               "horz,vert,diag",p);
+               };
+       };
+
+       
+}
+}
+       
+       
+#endif /*DISTANCETRANSFORMFILL_H_*/
diff --git a/drain/image/DistanceTransformOp.h b/drain/image/DistanceTransformOp.h
new file mode 100644 (file)
index 0000000..7f03d39
--- /dev/null
@@ -0,0 +1,438 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef DISTANCETRANSFORMOP_H_
+#define DISTANCETRANSFORMOP_H_
+
+
+#include <math.h>
+
+
+#include "ImageOp.h"
+
+namespace drain
+{
+
+namespace image
+{
+
+template <class T=unsigned char> 
+class DistanceModel {
+       
+public:
+       /** Sets the horizontal and vertical radius for the distance function.
+        *  It is recommended to make \c vert and \c diag optional. ie. apply default values.
+        *  By default, the geometry is octagonal, applyin 8-distance, but setting diag 
+        *  a large value makes it square (4-distance, diamond). 
+       */
+       virtual ~DistanceModel(){};
+
+       /// The values are in float, as eg. half-width radius may be sharp, under 1.0.
+       virtual 
+       void setGeometry(float horz,float vert,float diag) = 0;
+       
+       virtual 
+       T decreaseHorz(const T &x) const = 0; 
+
+       virtual 
+       T decreaseVert(const T &x) const = 0;
+
+       virtual 
+       T decreaseDiag(const T &x) const = 0;
+
+
+};
+
+template <class T=unsigned char> 
+class DistanceModelLinear : public DistanceModel<T> {
+
+       
+
+       virtual 
+       void setGeometry(float horz,float vert = -1,float diag = -1){
+               
+               if (horz <= 1)
+                 horzDecrement = Intensity::max<T>();
+               else
+                 horzDecrement = static_cast<T>(Intensity::max<T>()/horz);
+               
+               if (vert < 0)
+                       vertDecrement = horzDecrement;
+               else if (vert == 0)
+                       vertDecrement = Intensity::max<T>();
+               else // (vert > 0)      
+                       vertDecrement = static_cast<T>(Intensity::max<T>()/vert);
+               
+               if (diag < 0)
+                       diagDecrement =
+                       static_cast<T>(sqrt(horzDecrement*horzDecrement + vertDecrement*vertDecrement)) ;
+               else if (diag == 0)  
+                       diagDecrement = Intensity::max<T>();
+               else
+                       diagDecrement = static_cast<T>(Intensity::max<T>()/diag);
+                               
+               //cerr << "Linear: " << (float)horzDecrement << ", " << (float)vertDecrement << ", " << (float)diagDecrement << "\n";
+                       
+       } 
+       
+       // TODO setDecrements alternative...
+       
+       virtual 
+       inline T decreaseHorz(const T &x) const { return Intensity::limit<T>(x - horzDecrement); }; 
+
+       virtual 
+       inline T decreaseVert(const T &x) const { return Intensity::limit<T>(x - vertDecrement); };
+
+       virtual 
+       inline T decreaseDiag(const T &x) const { return Intensity::limit<T>(x - diagDecrement); };
+       
+protected:
+       
+       T horzDecrement;
+       T vertDecrement;
+       T diagDecrement;
+       
+};
+
+
+template <class T=unsigned char> 
+class DistanceModelExponential : public DistanceModel<T> {
+       
+       virtual 
+       void setGeometry(float horz,float vert = -1,float diag = -1){
+               
+               // TODO: interpret handle 0 and -1 better
+               if (horz <= 0.0)
+                       horzDecay = 0;
+                 //horz = 1;
+               else
+                 horzDecay = pow(0.5,1.0/horz);
+               
+               if (vert < 0)
+                 vertDecay = horzDecay;
+               else if (vert == 0)  // label => spread to infinity ??
+                 vertDecay = 1.0;
+               else   
+                 vertDecay = pow(0.5,1.0/vert);
+               
+
+               // TODO Fix a math bug here
+               if (diag == -1){
+                       //if ((horzDecay > 0) && (horzDecay > 0))
+                       const double hLog = log(horzDecay);
+                       const double vLog = log(vertDecay);
+                       diagDecay = exp(-sqrt(hLog*hLog + vLog*vLog));
+               }
+               else if (diag == 0)
+                       diagDecay = 1.0;
+               else
+                       diagDecay = pow(0.5,1.0/diag);
+                       
+               cerr << "Exponential: " << horzDecay << ", " << vertDecay << ", " << diagDecay << "\n";
+       } 
+       
+       // TODO setDecays, as alternative...
+       
+       virtual 
+       inline T decreaseHorz(const T &x) const { return static_cast<T>(horzDecay * x); }; 
+
+       virtual 
+       inline T decreaseVert(const T &x) const { return static_cast<T>(vertDecay * x); };
+
+       virtual 
+       inline T decreaseDiag(const T &x) const { return static_cast<T>(diagDecay * x); };
+       
+protected:
+       
+       double horzDecay;
+    double vertDecay;
+    double diagDecay;
+};
+
+
+
+               //const double d = this->parameters.get("diag",sqrt(1.0/(h*h) + 1.0/(v*v)));
+               
+
+/** Fast distance transform using 8-directional distance function.
+ *  Class D is used for computing distances?
+ */
+template <class T=unsigned char,class T2=unsigned char> //,class D=int>
+class DistanceTransformOp : public ImageOp<T,T2>
+{
+
+public:
+    
+    DistanceTransformOp(DistanceModel<T2> & model) : distanceModel(model) {};
+    
+    /*
+    DistanceTransformOp(const string & p = "5,5"){
+       this->setInfo("Creates areas of decreasing intensities.","horz,vert,diag",p);
+               //this->parameters.setDefaultKeys("horz,vert,diag");
+               //this->parameters.set(p);
+       };
+       */
+           
+    virtual ~DistanceTransformOp(){};
+
+
+       virtual void initialize() const {
+               float h = this->parameters.get("horz",3.0);
+               float v = this->parameters.get("vert",-1.0);
+               float d = this->parameters.get("diag",-1.0); //sqrt(h*h + v*v));
+               distanceModel.setGeometry(h,v,d);
+       }
+
+       
+
+       void filter(const Image<T> &src, Image<T2> &dst) const ;
+
+       void filterDownRight(const Image<T> &src, Image<T2> &dst) const ;
+       void filterUpLeft(const Image<T> &src, Image<T2> &dst) const ;
+
+protected:
+       mutable DistanceModel<T2> & distanceModel;
+       mutable DistanceModelLinear<T2> distanceLinear;
+       mutable DistanceModelExponential<T2> distanceExponential;
+};
+  
+  
+  
+template <class T,class T2>
+void DistanceTransformOp<T,T2>::filter(const Image<T> &src, Image<T2> &dst) const
+{
+               if (filterWithTmp(src,dst))
+                       return;
+
+        //cerr << "DTrf" << endl; // << (float)rectDecrement << endl;
+        initialize();
+               makeCompatible(src,dst);
+                       
+        //this->coordinateHandler->setBounds(src.getWidth(), src.getHeight());
+
+        filterDownRight(src,dst);
+        filterUpLeft(dst,dst);
+        //return dst;
+}
+  
+template <class T,class T2>
+void DistanceTransformOp<T,T2>::filterDownRight(const Image<T> &src, Image<T2> &dst) const
+{
+
+
+    const int width  = src.getWidth();
+    const int height = src.getHeight();
+       const CoordinateHandler & coordinateHandler = src.getCoordinateHandler();
+                       
+        // proximity (inverted distance)
+               float d;
+               float dPrev;
+                       
+               //int k = 0;
+
+        cerr << "UR" << (int)src.at(0,0) << endl;
+        cerr << "UR dst" << (int)dst.at(0,0) << endl;
+
+
+               // Point in the source image
+        Point2D<> p; 
+
+               // Point in the target image
+        Point2D<int> t;
+               int &tx = t.x;
+        int &ty = t.y;
+        
+        for (ty=0; ty<height; ty++)
+        {
+            for (tx=0; tx<width; tx++)
+            {
+
+                // Take source value as default
+                d = src.at(t);
+
+                // Compare to previous value
+                dPrev = dst.at(t);
+                if (dPrev > d)
+                {
+                    d = dPrev;
+                }
+
+                // Compare to upper left neighbour
+                p.setLocation(tx-1,ty-1);
+                coordinateHandler.handle(p);
+                //dPrev = dst.at(p) - diagdistanceModel.decrement;
+                dPrev = distanceModel.decreaseDiag(dst.at(p));
+                if (dPrev > d)
+                {
+                    d = dPrev;
+                }
+
+                // Compare to upper neighbour
+                p.setLocation(tx,ty-1);
+                coordinateHandler.handle(p);
+                //dPrev = dst.at(p) - vertdistanceModel.decrement;
+                dPrev = distanceModel.decreaseVert(dst.at(p));
+                if (dPrev > d)
+                {
+                    d = dPrev;
+                }
+
+                // Compare to upper right neighbour
+                p.setLocation(tx+1,ty-1);
+                coordinateHandler.handle(p);
+                //dPrev = dst.at(p) - diagdistanceModel.decrement;
+                dPrev = distanceModel.decreaseDiag(dst.at(p));
+                if (dPrev > d)
+                {
+                    d = dPrev;
+                }
+
+                // Compare to left neighbour
+                p.setLocation(tx-1,ty);
+                coordinateHandler.handle(p);
+                //dPrev = dst.at(p) - horzdistanceModel.decrement;
+                dPrev = distanceModel.decreaseHorz(dst.at(p));
+                if (dPrev > d)
+                {
+                    d = dPrev;
+                }
+
+                if (d>0)
+                       dst.at(t) = static_cast<T2>(d);
+                    
+
+            };
+        };
+        //return dst;
+
+};
+    
+
+    //public WritableRaster filterUpLeft(Raster srcDist, WritableRaster dst){
+    //template <class T=unsigned char,class T2=unsigned char>
+template <class T,class T2>
+void DistanceTransformOp<T,T2>::filterUpLeft(const Image<T> &src, Image<T2> &dst) const {
+
+        const int width  = src.getWidth();
+        const int height = src.getHeight();
+               const CoordinateHandler & coordinateHandler = src.getCoordinateHandler();
+                       
+        // proximity (inverted distance)
+       float d;
+               float dPrev;
+               
+        Point2D<> p(0,0);
+        //coordinateHandler.setBounds(width,height);
+        
+        // TODO:  target:
+        Point2D<> t;
+        int &tx = t.x;
+        int &ty = t.y;
+        
+        
+        for (ty=height-1; ty>=0; ty--)
+        {
+            for (tx=width-1; tx>=0; tx--)
+            {
+                // Source
+                d = src.at(t);
+
+                // Compare to previous value
+                dPrev = dst.at(t);
+                if (dPrev > d)
+                {
+                    d = dPrev;
+                }
+
+                // Compare to lower left neighbour
+                p.setLocation(tx-1,ty+1);
+                coordinateHandler.handle(p);
+                //dPrev = dst.at(p) - diagdistanceModel.decrement;
+                dPrev = distanceModel.decreaseDiag(dst.at(p));
+                if (dPrev > d)
+                {
+                    d = dPrev;
+                }
+
+                // Compare to lower neighbour
+                p.setLocation(tx,ty+1);
+                coordinateHandler.handle(p);
+                //dPrev = dst.at(p) - vertdistanceModel.decrement;
+                dPrev = distanceModel.decreaseVert(dst.at(p));
+                if (dPrev > d)
+                {
+                    d = dPrev;
+                }
+
+                // Compare to lower right neighbour
+                p.setLocation(tx+1,ty+1);
+                coordinateHandler.handle(p);
+                //dPrev = dst.at(p) - diagdistanceModel.decrement;
+                dPrev = distanceModel.decreaseDiag(dst.at(p));
+                if (dPrev > d)
+                {
+                    d = dPrev;
+                }
+
+                // Compare to right neighbour
+                p.setLocation(tx+1,ty);
+                coordinateHandler.handle(p);
+                //dPrev = dst.at(p) - horzDecrement;
+                dPrev = distanceModel.decreaseHorz(dst.at(p));
+                if (dPrev > d)
+                {
+                    d = dPrev;
+                }
+
+                if (d>0)
+                    dst.at(t) = static_cast<T2>(d);
+
+
+            }
+        }
+        //return dst;
+};
+
+template <class T=unsigned char,class T2=unsigned char> //,class D=int>
+class DistanceTransformLinearOp : public DistanceTransformOp<T,T2>
+{
+       public:
+       DistanceTransformLinearOp(const string & p = "5,5,7") : DistanceTransformOp<T,T2>(this->distanceLinear) {
+       this->setInfo("Creates areas of linearly decreasing intensities.","horz,vert,diag",p);
+       };
+
+};
+
+template <class T=unsigned char,class T2=unsigned char> //,class D=int>
+class DistanceTransformExponentialOp : public DistanceTransformOp<T,T2>
+{
+public:
+   DistanceTransformExponentialOp(const string & p = "5,5,7") :
+       DistanceTransformOp<T,T2>(this->distanceExponential) {
+       this->setInfo("Creates areas of exponentially decreasing intensities. Set half-widths.",
+       "horz,vert,diag",p);
+   };
+};         
+
+}
+}
+       
+#endif /*DISTANCETRANSFORMOP_H_*/
diff --git a/drain/image/DoubleSmootherOp.h b/drain/image/DoubleSmootherOp.h
new file mode 100644 (file)
index 0000000..35021dc
--- /dev/null
@@ -0,0 +1,88 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef Double_SMOOTHEROP
+#define Double_SMOOTHEROP "DoubleSmootherOp Markus.Peura@iki.fi"
+
+#include "WindowOp.h"
+
+#include "FastAverageOp.h"
+#include "MathOpPack.h"
+//#include "SlidingStripeAverageOp.h"
+
+namespace drain
+{
+namespace image
+{
+
+
+//  TODO delicateSmoother (quadraticsmoother)
+// -------------------------------------------------
+/*! This operator produces
+ *  <CODE>
+ *  F2 = cF + (1-c)M{ cF + (1-c)M{F} }
+ *     = cF + (1-c)cM{F} + (1-c)^2M^2{F}
+ *  </CODE>
+ *  where 
+ *  <CODE>F</CODE> is an image,
+ *  <CODE>M{}</CODE> is a <CODE>W x W</CODE> WindowAverage operator, and
+ *  <CODE>c</CODE> is a coefficient.
+ */
+template <class T=unsigned char,class T2=unsigned char>
+class DoubleSmootherOp: public WindowOp<T,T2>
+{
+
+public:
+
+       DoubleSmootherOp(const string & p = "5,5,0.5") :
+               WindowOp<T,T2>("DoubleSmootherOp",
+                               "Smoothes image twice, mixing first result with coeff*100%.",
+                               "width,height,coeff", p){
+       };
+
+       virtual void initialize() const {
+       }
+
+       virtual void filter(const Image<T> &src,Image<T2> &dst) const {
+
+               const int width  = this->getParameter("width",3);
+               const int height = this->getParameter("height",width);
+               const double coeff  = this->getParameter("coeff",0.5);
+
+               Image<T2> tmp;
+               FastAverageOp<T,T2>(width,height).filter(src,tmp);
+               MixerOp<T,T2>(coeff).filter(src,tmp,dst);
+               //FastAverageOp<T,T2>(width,height).filter(tmp,dst);
+
+       }
+
+
+};
+
+
+
+
+} // namespace drain
+
+} //namespace image
+
+
+#endif
diff --git a/drain/image/Drain.h b/drain/image/Drain.h
new file mode 100644 (file)
index 0000000..1d75bfa
--- /dev/null
@@ -0,0 +1,295 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef DRAIN_H_
+#define DRAIN_H_
+
+#include <string>
+#include <list>
+#include <map>
+#include <set>
+
+#include <exception>
+
+
+#include "Image.h"
+#include "CatenatorOp.h"
+#include "CopyOp.h"
+#include "DoubleSmootherOp.h"
+#include "DistanceTransformOp.h"
+//#include "DistanceTransformLogOp.h"
+#include "DistanceTransformFillOp.h"
+#include "FastAverageOp.h"
+#include "FastOpticalFlowOp.h"
+#include "FloodFillOp.h"
+#include "FuzzyPeakOp.h"
+#include "FuzzyThresholdOp.h"
+#include "GammaOp.h"
+#include "GradientOp.h"
+#include "HighPassOp.h"
+#include "MarginalStatisticOp.h"
+#include "MathOpPack.h"
+#include "PaletteOp.h"
+#include "PixelVectorOp.h"
+#include "QuadraticSmootherOp.h"
+#include "QuantizatorOp.h"
+#include "RecursiveRepairerOp.h"
+#include "RunLengthOp.h"
+#include "SegmentAreaOp.h"
+#include "SlidingStripeOp.h"
+//#include "SlidingStripeAverageOp.h"
+//#include "SlidingStripeAverageWeightedOp.h"
+#include "SlidingWindowMedianOp.h"
+#include "ThresholdOp.h"
+
+//#include ""
+
+namespace drain
+{
+
+namespace image
+{
+
+///   A class that wraps various image processing operators into a service instance.
+/**
+ *   
+ * 
+ *  \section Examples
+ * 
+ *  \subsection Scaling intensities
+ * 
+ *  \code
+ *   drainage test-grey.png --threshold 128 -o threshold.png
+ *  \endcode
+ *  
+ *  \code
+ *   drainage test-grey.png --fuzzyThreshold 128,5 -o fuzzyThreshold.png
+ *  \endcode 
+ *  
+ *  Correspondingly, intensities around a given intensity can be highlighted 
+ *  with a \em fuzzy \em peak.
+ *    
+ *  \code
+ *   drainage test-grey.png --fuzzyPeak 128,32 -o fuzzyPeak.png
+ *  \endcode 
+ *
+ *  \subsection Window operations
+ * 
+ *   
+ *
+ */
+template <class T=unsigned char>  
+class Drain : private map<string,ImageOp<T,T> *> {
+       
+protected:
+       
+       AdditionOp<T,T> add;
+       CatenatorOp<T,T> cat;
+       CopyOp<T,T> copy;
+       DoubleSmootherOp<T,T> doubleSmoother;
+       FastAverageOp<T,T> fastAverage; //wavgWindow;   // window only
+       GradientHorizontalOp<T,T> gradientHorz;
+       GradientVerticalOp<T,T>   gradientVert;
+       DistanceOp<T,T> distance;
+       DistanceTransformLinearOp<> distanceTransform;
+       DistanceTransformExponentialOp<> distanceTransformLog;
+       DistanceTransformFillOp<> distanceTransformFill;
+       DistanceTransformFillExponentialOp<> distanceTransformFillExp;
+       DivisionOp<T,T> div;
+       GammaOp<T,T> gamma;
+       HighPassOp<T,T> highPass;
+       FloodFillOp<T,T> floodFill;
+       FuzzyPeakOp<T,T> fuzzyPeak;
+       FuzzyThresholdOp<T,T> fuzzyThreshold;
+       MaximumOp<T,T> max;
+       MagnitudeOp<T,T> magnitude;
+       MarginalStatisticOp<T,T> marginStat;
+       MinimumOp<T,T> min;
+       MixerOp<T,T> mix;
+       MultiplicationOp<T,T> mul;
+       NegateOp<T,T> negate;
+       // Palette?
+       ProductOp<T,T> product;
+       QuantizatorOp<T,T> quantize;
+       QuadraticSmootherOp<T,T> quadSmooth;
+       RecursiveRepairerOp<T,T> rec;
+       RemapOp<T,T> remap;
+       RunLengthOp<T,T> runLength;
+       ScaleOp<T,T> rescale;
+       SegmentAreaOp<T,T> segmentArea;
+       SlidingWindowMedianOp<T,T> median;
+       SubtractionOp<T,T> sub;
+       ThresholdOp<T,T> threshold;
+       
+public:
+
+       //: average(avgWindow) , averageWeighted(wavgWindow), median(medianWindow) {
+       Drain(){};
+       
+       virtual void addDefaultOps(){
+               // Simple math ops
+               // addOperator("",);
+               addOperator("add",add);
+               addOperator("average",fastAverage);
+               addOperator("cat",cat);
+               addOperator("copy",copy);
+               addOperator("gradientHorz",gradientHorz);
+               addOperator("gradientVert",gradientVert);
+               addOperator("distance",distance);
+               addOperator("distanceTransform",distanceTransform);
+               addOperator("distanceTransformLog",distanceTransformLog);
+               addOperator("distanceTransformFill",distanceTransformFill);
+               addOperator("distanceTransformFillExp",distanceTransformFillExp);
+               addOperator("div",div);
+               addOperator("doubleSmoother",doubleSmoother);
+               addOperator("gamma",gamma);
+               addOperator("highPass",highPass);
+               addOperator("floodFill",floodFill);
+               addOperator("fuzzyPeak",fuzzyPeak);
+               addOperator("fuzzyThreshold",fuzzyThreshold);
+               addOperator("max",max);
+               addOperator("magnitude",magnitude);
+               addOperator("marginalStatistic",marginStat);
+               addOperator("min",min);
+               addOperator("mix",mix);
+               addOperator("median",median);
+               addOperator("mul",mul);
+               addOperator("negate",negate);
+               addOperator("product",product);
+               addOperator("quantize",quantize);
+               addOperator("quadSmooth",quadSmooth);
+               addOperator("rescale",rescale);
+               addOperator("remap",remap);
+               addOperator("rec",rec);
+               addOperator("runLength",runLength);
+               addOperator("segmentArea",segmentArea);
+               addOperator("sub",sub);
+               addOperator("threshold",threshold);
+       };
+       
+       virtual ~Drain(){};
+       
+       void prefixAll(const string &prefix, bool upperCaseFirst=true){
+
+               set<string> keys;
+
+               for (typename Drain<T>::iterator it = this->begin(); it != this->end(); it++)
+                       keys.insert(it->first);
+
+               for (typename set<string>::iterator it = keys.begin();it != keys.end(); it++){
+                       string newKey = *it;
+                       if (upperCaseFirst)
+                               if (!newKey.empty())
+                                       newKey[0] = newKey[0]+'A'-'a';
+                       (*this)[prefix+newKey] = (*this)[*it];
+                       this->erase(*it);
+               }
+       }
+
+       /// Adds an image processing operator to the internal map. 
+       /// Thereafter the operator can be executed with executeOperator().
+       void addOperator(const string &name,ImageOp<T,T> &op){
+               (*this)[name] = & op;
+       }
+       
+       void removeOperator(const string &name){
+               this->erase(name);
+       }
+
+       /// 
+       bool hasOperator(const string &key) const {
+               return (this->find(key) != this->end());
+       }
+       
+       
+       /// An interface for single-image processing  
+       int executeOperator(const string & name, const string & parameters,
+               const Image<T> & src, Image<T> & dst){
+               if (hasOperator(name)){
+                       if (drain::Debug > 0)
+                               cerr << "executing: '" << name << "'\n";
+                       ImageOp<T,T> & op = *(*this)[name];
+                       op.parameters.set(parameters); 
+                       op.filter(src,dst);
+               }
+               else {
+                 throw runtime_error("drain::Drain: Could not find operator " + name);
+               }               
+               return 0;
+       };
+       
+
+       /// Lists the installed operator names, parameters, and descriptions. If key is given, returns only help for that command.
+       void help(ostream &ostr, const string & key = "", const string & prefix = "--", 
+                       const string & postfix = "\n\t") const {
+               if (key == ""){
+                       for (typename map<string,ImageOp<T,T> *>::const_iterator it = this->begin();
+                                       it != this->end(); it++){
+                               ostr << prefix << it->first;
+                               if ( !it->second->parameters.getParameterNames().empty() ) //it->second->parameters.isFlag )
+                                       //      ostr << it->second->parameters.info << '\n';
+                                       //else
+                                       ostr << " <" << it->second->parameters.getParameterNames() << '>';
+                               ostr << postfix;
+                               ostr << it->second->parameters.getDescription() << '\n';
+                       //it->second->help(ostr);
+                 }
+               }
+               else {
+                       if (hasOperator(key)){
+                               typename map<string,ImageOp<T,T> *>::const_iterator it = this->find(key);
+                               if (it != this->end()){
+                                       it->second->help(ostr);
+                               }
+                       }
+               }
+    } 
+       
+       /// Convenience
+    string help() const {
+       stringstream s;
+       help(s);
+       return s.str(); 
+    }
+       
+               
+       //int handle(const string & command, const string & parameters,
+       //      const list<BufferedImage<T> > & srcList, list<BufferedImage<T> > & dstList);
+};
+
+/*
+template <class T>  
+int Drain<T>::handleImage(const string & command, const string & parameters,
+               const BufferedImage<T> & src, BufferedImage<T> & dst){
+       list<BufferedView<T> > srcList(1);
+       list<BufferedView<T> > dstList(1);
+       srcList.front().viewImage(src);
+       dstList.front().viewImage(dst);
+       return handle(command,parameters,srcList,dstList);
+}
+*/
+               
+
+}
+
+}
+
+#endif /*DRAIN_H_*/
diff --git a/drain/image/FastAverageOp.h b/drain/image/FastAverageOp.h
new file mode 100644 (file)
index 0000000..6b654e9
--- /dev/null
@@ -0,0 +1,374 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef FASTAVERAGEOP_H_
+#define FASTAVERAGEOP_H_
+
+#include "SlidingStripeOp.h"
+#include "WindowOp.h"
+
+namespace drain
+{
+
+namespace image
+{
+
+
+template <class T=unsigned char,class T2=unsigned char>
+class SlidingStripeAverage : public SlidingStripe<T,T2>
+//public SlidingWindowAverage<T,T2>
+{
+public:
+
+       double sum;
+       float areaF;
+
+        SlidingStripeAverage(int width=1,int height=1) {
+          this->setSize(width,height);
+        };
+
+        //virtual ~SlidingWindowStripe(){};
+
+
+       /**
+        *   TODO: weight = > weightSUm
+        *   TODO: multiple channels (banks)?
+        */
+        virtual void initialize(){
+
+               this->location.setLocation(0,0);;
+               this->areaF = this->getArea();
+
+               cerr << "area" << this->areaF << "\n";
+               cerr << this->dst.getCoordinateHandler() << '\n';
+
+               this->sum = 0;
+               if (this->horizontal){
+                       for (int i = this->iMin; i <= this->iMax; i++) {
+                               this->p.setLocation(i,0);
+                               this->dst.getCoordinateHandler().handle(this->p);
+                               this->sum += this->src.at(this->p);
+                       }
+               }
+               else {
+                       for (int j = this->jMin; j <= this->jMax; j++) {
+                               this->p.setLocation(0,j);
+                               this->dst.getCoordinateHandler().handle(this->p);
+                               this->sum += this->src.at(this->p);
+                       }
+               }
+       }
+
+
+
+       /// Kesken
+       virtual void updateHorz(int dx){
+
+               const int xOld = this->getXOld(dx);
+               const int xNew = this->getXNew(dx);
+
+               if (this->horizontal){
+                       this->p.setLocation(xOld,this->location.y);
+                       this->dst.getCoordinateHandler().handle(this->p);
+                       this->sum -= this->src.at(this->p);
+
+                       this->p.setLocation(xNew,this->location.y);
+                       this->dst.getCoordinateHandler().handle(this->p);
+                       this->sum += this->src.at(this->p);
+
+               }
+               else
+               {
+                       this->sum = 0;
+                       for (int j = this->jMin; j <= this->jMax; j++) {
+                               this->p.setLocation(xNew,this->location.y+j);
+                               this->dst.getCoordinateHandler().handle(this->p);
+                               this->sum += this->src.at(this->p);
+                       }
+
+               }
+
+       }
+
+       /// Kesken
+       virtual void updateVert(int dy){
+
+               const int yOld = this->getYOld(dy);
+               const int yNew = this->getYNew(dy);
+
+               if (this->horizontal){
+                       this->sum = 0;
+                       for (int i = this->iMin; i <= this->iMax; i++) {
+                               this->p.setLocation(this->location.x+i,yNew);
+                               this->dst.getCoordinateHandler().handle(this->p);
+                               this->sum += this->src.at(this->p);
+                       }
+               }
+               else {
+
+                       this->p.setLocation(this->location.x,yOld);
+                       this->dst.getCoordinateHandler().handle(this->p);
+                       this->sum -= this->src.at(this->p);
+
+                       this->p.setLocation(this->location.x,yNew);
+                       this->dst.getCoordinateHandler().handle(this->p);
+                       this->sum += this->src.at(this->p);
+               }
+
+       }
+
+       virtual void write(){
+               this->dst.at(this->location) =  static_cast<T2>(sum/areaF);
+       }
+
+       //protected:
+       //      bool horizontal;
+
+};
+
+/// Supports using another image (channel) for weighting the pixels of the source image.
+/**
+ *
+ *
+ *
+ */
+template <class T=unsigned char,class T2=unsigned char>
+class SlidingStripeAverageWeighted : public SlidingStripeAverage<T,T2>
+{
+public:
+
+       double sum;
+       float sumW;
+       float areaF;
+
+       SlidingStripeAverageWeighted(int width=1,int height=1) : SlidingStripeAverage<T,T2>(width,height) {};
+
+
+       /**
+        *   TODO: weight = > weightSUm
+        *   TODO: multiple channels (banks)?
+        */
+        virtual void initialize(){
+
+               static T w;
+       
+               this->location.setLocation(0,0);
+               this->sum = 0;
+               this->sumW = 0;
+               this->areaF = this->getArea();
+               
+               if (this->horizontal){
+                       for (int i = this->iMin; i <= this->iMax; i++) {
+                               this->p.setLocation(i,0);
+                               this->dst.getCoordinateHandler().handle(this->p);
+                               w = this->srcWeight.at(this->p);
+                               this->sum  += w*this->src.at(this->p);
+                               this->sumW += w;
+                       }
+               }
+               else {
+                       for (int j = this->jMin; j <= this->jMax; j++) {
+                               this->p.setLocation(0,j);
+                               this->dst.getCoordinateHandler().handle(this->p);
+                               w = this->srcWeight.at(this->p);
+                               this->sum  += w*this->src.at(this->p);
+                               this->sumW += w;
+                       }
+               }
+       }
+
+
+
+       /// Kesken
+       virtual void updateHorz(int dx){
+               
+               static T w;
+               const int xOld = this->getXOld(dx);
+               const int xNew = this->getXNew(dx);
+
+               if (this->horizontal){
+                       this->p.setLocation(xOld,this->location.y);
+                       this->dst.getCoordinateHandler().handle(this->p);
+                       w = this->srcWeight.at(this->p);
+                       this->sum  -= w*this->src.at(this->p);
+                       this->sumW -= w;
+
+                       this->p.setLocation(xNew,this->location.y);
+                       this->dst.getCoordinateHandler().handle(this->p);
+                       w = this->srcWeight.at(this->p);
+                       this->sum  += w*this->src.at(this->p);
+                       this->sumW += w;
+               }
+               else
+               {
+                       this->sum  = 0;
+                       this->sumW = 0;
+                       for (int j = this->jMin; j <= this->jMax; j++) {
+                               this->p.setLocation(xNew,this->location.y+j);
+                               this->dst.getCoordinateHandler().handle(this->p);
+                               w = this->srcWeight.at(this->p);
+                               this->sum  += w*this->src.at(this->p);
+                               this->sumW += w;
+                       }
+
+               }
+
+       }
+
+       /// Kesken
+       virtual void updateVert(int dy){
+
+               const int yOld = this->getYOld(dy); // location.y + (dy>0 ? jMin-1 : jMax+1);
+               const int yNew = this->getYNew(dy); // location.y + (dy>0 ? jMax     : jMin);
+               static T w;
+
+               if (this->horizontal){
+                       this->sum  = 0;
+                       this->sumW = 0;
+                       for (int i = this->iMin; i <= this->iMax; i++) {
+                               this->p.setLocation(this->location.x+i,yNew);
+                               this->dst.getCoordinateHandler().handle(this->p);
+                               w = this->srcWeight.at(this->p);
+                               this->sum  += w*this->src.at(this->p);
+                               this->sumW += w;
+                       }
+               }
+               else {
+
+                       this->p.setLocation(this->location.x,yOld);
+                       this->dst.getCoordinateHandler().handle(this->p);
+                       w = this->srcWeight.at(this->p);
+                       this->sum  -= w*this->src.at(this->p);
+                       this->sumW -= w;
+
+                       this->p.setLocation(this->location.x,yNew);
+                       this->dst.getCoordinateHandler().handle(this->p);
+                       w = this->srcWeight.at(this->p);
+                       this->sum  += w*this->src.at(this->p);
+                       this->sumW += w;
+               }
+
+       }
+
+       virtual void write(){
+               if (sumW > 0)
+                       this->dst.at(this->location) =  static_cast<T2>(sum/sumW);
+               else
+                       this->dst.at(this->location) = 0;
+               this->dstWeight.at(this->location) =  static_cast<T2>(sumW/areaF);
+               //this->dst.at(this->location) =  128; 
+       }
+
+       /// Sets the source channel from which the pixel weights will be read.
+       /// Sets the target channel in which resulting weight will be written.
+
+};
+
+
+
+/// If source images contain alpha channels, adopts to weighted mode.
+
+template <class T=unsigned char,class T2=unsigned char>
+class FastAverageOp : public WindowOp<T,T2> //public SlidingStripeOp<T,T2>
+{
+public:
+       
+       FastAverageOp(const string & p = "")  { //: SlidingStripeOp<T,T2>(stripeAvg) {
+               this->setInfo("Weighted window averaging operator.","width,height",p);
+       };
+       
+       FastAverageOp(int width,int height) { //: SlidingStripeOp<T,T2>(stripeAvg) {
+               this->setInfo("Weighted window averaging operator.","width,height","0,0");
+               this->setParameter("height",height);
+               this->setParameter("width",width);
+       };
+
+       virtual void filter(const Image<T> &src,Image<T2> &dst) const {
+               
+               makeCompatible(src,dst);
+               
+               // UNWEIGHTED
+               if (src.getAlphaChannelCount() == 0){
+                       cerr << "ssaw:filter0\n";
+                       unsigned int w = this->getParameter("width",3);
+                       unsigned int h = this->getParameter("height",w);
+                       SlidingStripeAverage<T,T2> window1; //(w,h);
+                       SlidingStripeAverage<T2,T2> window2; //(w,h);
+                       SlidingStripeOp<T,T2> op(window1,window2,"width,height");
+                       op.setParameter("width",w);
+                       op.setParameter("height",h);
+                       //cerr << "SSAW calling SlidingStripeOp/f2 with " << window1 << "\n";
+                       //cerr << "SSAW calling SlidingStripeOp/f2 with " << window2 << "\n";
+                       for (unsigned int i=0; i<src.getChannelCount(); i++){
+                               op.filter(src.getChannel(i),dst.getChannel(i)); 
+                       }
+               }
+               // WEIGHTED
+               else {
+                       cerr << "SSAW calling f2->f4\n";
+                       filter(src.getImageChannels(),src.getAlphaChannels(),
+                               dst.getImageChannels(),dst.getAlphaChannels());
+               }
+               
+       };
+       
+       /// Calls SlidingStripeOp<T,T2>::filter() separately for each image channel. This is natural for many operations, such as averaging.
+
+       // Raise filterUnweighted       
+       virtual void filter(const Image<T> &src,const Image<T> &srcWeight,
+               Image<T2> &dst,Image<T2> &dstWeight) const {
+               
+               makeCompatible(src,dst);        
+               makeCompatible(srcWeight,dstWeight);    
+               
+               //src.setName("src");
+               //srcWeight.setName("srcWeight");
+               //dst.setName("dst");
+               //dstWeight.setName("dstWeight");
+               
+               unsigned int w = this->getParameter("width",3);
+               unsigned int h = this->getParameter("height",w);
+               SlidingStripeAverageWeighted<T,T2> window1; //(w,h);
+               SlidingStripeAverageWeighted<T2,T2> window2; //(w,h);
+               SlidingStripeOp<T,T2> op(window1,window2,"width,height");
+               op.setParameter("width",w);
+               op.setParameter("height",h);
+                       
+               const unsigned int imageChannels = src.getImageChannelCount();  
+               const unsigned int alphaChannels = srcWeight.getChannelCount();
+               for (unsigned int i=0; i<imageChannels; i++){   
+                       const unsigned int a = i % alphaChannels;
+                       cerr << "SSAW calling SlidingStripeOp/4 " << i << ':' << a << " \n";
+                       op.filter(src.getChannel(i),srcWeight.getChannel(a),
+                                               dst.getChannel(i),dstWeight.getChannel(a));
+               }       
+               
+       }
+       
+protected:
+       SlidingStripeAverageWeighted<T,T2> stripeAvg;   
+};
+
+}
+
+}
+
+#endif /*SLIDINGWINDOWAVERAGESTRIPEWEIGHTED_H_*/
diff --git a/drain/image/FastOpticalFlowOp.h b/drain/image/FastOpticalFlowOp.h
new file mode 100644 (file)
index 0000000..2d0239b
--- /dev/null
@@ -0,0 +1,500 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef SLIDINGWINDOWOPTICALFLOWOP_H_
+#define SLIDINGWINDOWOPTICALFLOWOP_H_
+
+#include <sstream>
+
+#include "SlidingWindow.h"
+#include "SlidingWindowOp.h"
+
+//#include "FastAverageOp.h"
+#include "DoubleSmootherOp.h"
+#include "FuzzyPeakOp.h"
+#include "GradientOp.h"
+#include "PixelVectorOp.h"
+#include "QuadraticSmootherOp.h"
+#include "RecursiveRepairerOp.h"
+
+namespace drain
+{
+
+namespace image
+{
+
+/// Consider using SlidingStripeOpticalFlow instead
+/**
+ *  T  = class of the input
+ *  T2 = class of the derivative  images, weight, and output (motion field and weight)
+ */
+ template <class T=double> // ,class T2=double> // NOTICE
+class SlidingOpticalFlow : public SlidingWindow<T,T> {
+
+  public:
+       
+       
+       double wXX;
+       double wXY;
+       double wYY;
+       double wXT;
+       double wYT;
+       //double area;
+
+       /// Precomputed horizontal diffential 
+       ImageView<T> Dx;
+       
+       /// Precomputed vertical diffential 
+       ImageView<T> Dy;
+       
+       /// Precomputed time diffential 
+       ImageView<T> Dt;
+       
+       /// Precomputed weight field (optional) 
+       ImageView<T> w;
+       
+       ImageView<T> u;
+       ImageView<T> v;
+       ImageView<T> q;
+       
+
+       SlidingOpticalFlow(int width = 0, int height = 0){
+               this->setSize(width,height);
+       }
+
+       ///
+       void setDx(Image<T> &image, unsigned int k = 0){
+               Dx.viewChannel(image,k);
+       }
+
+       void setDy(Image<T> &image, unsigned int k = 0){
+               Dy.viewChannel(image,k);
+       }
+
+
+       void setDt(Image<T> &image, unsigned int k = 0){
+               Dt.viewChannel(image,k);
+       }
+
+       /// Sets quality channel
+       void setW(Image<T> &image, unsigned int k = 0){
+               w.viewChannel(image,k);
+       }
+
+       /// Given 4-channel image, associates Dx, Dy, Dt and w to the channels.
+       void setDiff(const Image<T> &image){
+               Dx.viewChannel(image,0);
+               Dy.viewChannel(image,1);
+               Dt.viewChannel(image,2);
+               w.viewChannel(image,3);
+       }
+
+
+       // protected:
+       double areaD;
+       
+       
+       /** 
+        *   TODO: weight = > weightSUm
+        *   TODO: multiple channels (banks)?
+        */
+       virtual void initialize() {
+
+               u.viewChannel(this->dst,0);
+               v.viewChannel(this->dst,1);
+               q.viewChannel(this->dst,2);
+
+               areaD = this->getArea();
+               
+               //------------------------------
+               wXX = 0.0;
+       wXY = 0.0;
+       wYY = 0.0;
+       wXT = 0.0;
+       wYT = 0.0;
+               
+       static T DX, DY, DT, W;
+
+               this->location.setLocation(0,0);
+               // sum = 0;
+               // areaF = this->getArea();
+               // cerr << "area" << areaF << "\n";
+               cerr << this->dst.getCoordinateHandler();
+               
+               static unsigned int address;
+               
+               for (int i = this->iMin; i <= this->iMax; i++) {
+                       for (int j = this->jMin; j <= this->jMax; j++) {
+
+                               this->p.x = i;
+                               this->p.y = j;
+                               this->dst.getCoordinateHandler().handle(this->p);
+
+                               address = this->dst.address(this->p.x,this->p.y);
+
+                               DX = Dx.at(address);
+                               DY = Dy.at(address);
+                               DT = Dt.at(address);
+                               //W  = 1; // w.at(address);
+                               W = w.at(address);
+                               
+                               wXX += W*DX*DX;
+               wXY += W*DX*DX;
+               wYY += W*DY*DY;
+               wXT += W*DX*DT;
+               wYT += W*DY*DT;
+                               
+                       }
+               }
+               
+       };
+       
+       
+       
+       virtual void updateHorz(int dx){
+               
+               static unsigned int address;
+               static T DX, DY, DT, W;
+               static int xOld, xNew;
+
+               xOld = this->getXOld(dx); // location.x + (dx>0 ? iMin-1 : iMax+1);
+               xNew = this->getXNew(dx); // location.x + (dx>0 ? iMax     : iMin);
+
+               for (int j = this->jMin; j <= this->jMax; j++) {
+                       
+                       this->p.x = xOld;
+                       this->p.y = this->location.y+j;
+                       
+                       this->dst.getCoordinateHandler().handle(this->p);
+
+                       address = this->dst.address(this->p.x,this->p.y);
+
+                       DX = Dx.at(address);
+                       DY = Dy.at(address);
+                       DT = Dt.at(address);
+                       //W = 1;
+                       W  = w.at(address);
+                       
+                       wXX -= W*DX*DX;
+               wXY -= W*DX*DY;
+               wYY -= W*DY*DY;
+               wXT -= W*DX*DT;
+               wYT -= W*DY*DT;
+                       // sum -= this->src.at(this->p);
+                       
+                       this->p.x = xNew;
+                       this->p.y = this->location.y+j;
+                       
+                       this->dst.getCoordinateHandler().handle(this->p);
+
+                       address = this->dst.address(this->p.x,this->p.y);
+
+                       DX = Dx.at(address);
+                       DY = Dy.at(address);
+                       DT = Dt.at(address);
+                       W  =  w.at(address);
+
+                       wXX += W*DX*DX;
+               wXY += W*DX*DY;
+               wYY += W*DY*DY;
+               wXT += W*DX*DT;
+               wYT += W*DY*DT;
+                       
+               }
+       };
+       
+       // KESKEN
+       virtual void updateVert(int dy){
+               
+               static unsigned int address;
+               static T DX, DY, DT, W;
+               static int yOld, yNew;
+
+               yOld = this->getYOld(dy); // location.y + (dy>0 ? jMin-1 : jMax+1);
+               yNew = this->getYNew(dy); // location.y + (dy>0 ? jMax     : jMin);
+               
+               for (int i=this->iMin; i<=this->iMax; i++) {
+
+                       this->p.x = this->location.x+i;
+                       this->p.y = yOld;
+                       this->dst.getCoordinateHandler().handle(this->p);
+                       
+                       // tsekkaamatta!
+                       address = this->dst.address(this->p.x,this->p.y);
+
+                       DX = Dx.at(address);
+                       DY = Dy.at(address);
+                       DT = Dt.at(address);
+                       //W  =  1;
+                       W  =  w.at(address);
+
+                       wXX -= W*DX*DX;
+                       wXY -= W*DX*DY;
+                       wYY -= W*DY*DY;
+                       wXT -= W*DX*DT;
+                       wYT -= W*DY*DT;
+
+                       this->p.x = this->location.x+i;
+                       this->p.y = yNew;
+                       this->dst.getCoordinateHandler().handle(this->p);
+
+                       address = this->dst.address(this->p.x,this->p.y);
+
+                       DX = Dx.at(address);
+                       DY = Dy.at(address);
+                       DT = Dt.at(address);
+                       W  =  w.at(address);
+
+                       wXX += W*DX*DX;
+                       wXY += W*DX*DY;
+                       wYY += W*DY*DY;
+                       wXT += W*DX*DT;
+                       wYT += W*DY*DT;
+               }
+               
+       }
+
+       inline T nominator() const
+    {
+        return static_cast<T>(wXX*wYY - wXY*wXY + .01);  // TODO "epsilon", ensures positivity
+    };
+
+    /// Returns the horizontal component of motion. Must be scaled by nominator().
+    inline T uDenominator() const
+    {
+        return static_cast<T>(wXY*wYT - wYY*wXT);
+    };
+
+       /// Returns the vertical component of motion. Must be scaled by nominator().
+    inline T vDenominator() const
+    {
+        return static_cast<T>(wXY*wXT - wXX*wYT);
+    };
+
+
+       virtual void write(){
+               static T nom;
+               static double quality;
+
+               nom = nominator();
+               quality = sqrt(nom/(areaD*areaD*2048));  // hihasta kerroin
+
+               if (quality > 1.0){  // todo threshold
+                       u.at(this->location) = uDenominator()/nom;
+                       v.at(this->location) = vDenominator()/nom;
+                       q.at(this->location) = static_cast<T>( 255.0 - 255.0 / (1.0+quality) );
+               }
+               else {
+                       u.at(this->location) = 0;
+                       v.at(this->location) = 0;
+                       q.at(this->location) = 0;
+               }
+
+       }
+       
+       
+  };
+
+
+//------------------------------------------------------------------------------------------
+
+
+/// Detects motion between two subsequent images. Does not extrapolate images.
+/// Applies recursive correction for smoothing motion field.
+/**  Notice:
+ *   T  input data (images)
+ *   T2 all the other arrays (differentials, motion)
+ *
+ *   - gradPow - highlighting the role of gradient magnitude
+ *   - gradWidth -
+ */
+template <class T=unsigned char,class T2=double>
+class FastOpticalFlowOp : public SlidingWindowOp<T2,T2>
+{
+public:
+
+       FastOpticalFlowOp(const string & p = "5,5,0.5,2,16") :
+               SlidingWindowOp<T2,T2>(opticalFlowWindow,"OpticalFlow","A pipeline implementation of optical flow. ","width,height,smoothing,gradPow,gradWidth",p){
+
+       };
+
+       /// Creates a 2+1 channel target image for storing motion (u,v) and quality (q).
+       //  Notice T2,T2.
+       virtual void makeCompatible(const Image<T2> &src,Image<T2> &dst) const  {
+               unsigned int width = src.getWidth();
+               unsigned int height = src.getHeight();
+               dst.setGeometry(width,height,2,1);
+       dst.getCoordinateHandler().setBounds(width,height);  // ??
+       };
+
+       /// Computes an image with channels dx, dy, dt and w (quality of gradients). User may wish to redefine this.
+       ///
+       virtual void computeDerivativeImage(const Image<T> &src1,const Image<T> &src2,Image<T2> &dst) const {
+
+               // TODO: concentrate on the "middle image". Skip grad stability, let oflow use
+
+               // Window properties
+               const unsigned int width  = this->getParameter("width",3);
+               const unsigned int height = this->getParameter("height",width);
+               const double smoothingCoeff  = this->getParameter("smoothing",0.5);
+
+               const unsigned int imageWidth = src1.getWidth();
+               const unsigned int imageHeight = src1.getHeight();
+               dst.setGeometry(imageWidth,imageHeight,3,1);
+
+               const ImageView<T2> grad(dst,0,2);  // dx & dy
+               Image<T2> & dx = dst.getChannel(0);
+               Image<T2> & dy = dst.getChannel(1);
+               Image<T2> & dt = dst.getChannel(2);
+               Image<T2> &  w = dst.getChannel(3);
+
+               Image<T2> src1Smooth; //(width,height,2);
+               Image<T2> src2Smooth; //(width,height,2);
+
+               // Notice that weighted smoothing applies, if applicable
+               QuadraticSmootherOp<T,T2> smoother("5,5,0.5"); // bug
+               // FastAverageOp<T,T2> smoother("5,5"); // bug
+               //DoubleSmootherOp<T,T2> smoother("5,5,0.5");
+               smoother.setSize(width/2,height/2);  // half of oflow windows
+               smoother.setParameter("coeff",smoothingCoeff);
+               smoother.filter(src1,src1Smooth);
+               smoother.filter(src2,src2Smooth);
+
+               if (drain::Debug > 3){
+                       Image<T> tmp;
+                       CopyOp<T2,T>().filter(src1Smooth,tmp); File::write(tmp,"grad-s1.png");
+                       CopyOp<T2,T>().filter(src2Smooth,tmp); File::write(tmp,"grad-s2.png");
+               }
+
+               // Time derivative, dt
+               SubtractionOp<T2,T2>().filter(src2Smooth.getChannel(0),src1Smooth.getChannel(0),dt);
+
+
+               /// Gradient quality = stability * magnitude
+               // (part 1: gradient unchangedness)
+               // TEST ERAD
+
+               GradientHorizontalOp<T2,T2>().filter(src2Smooth.getChannel(0),dx);
+               GradientVerticalOp<T2,T2>().filter(src2Smooth.getChannel(0),dy);
+               Image<T2> gradTemp(imageWidth,imageHeight,2);
+               GradientHorizontalOp<T2,T2>().filter(src1Smooth.getChannel(0),gradTemp.getChannel(0));
+               GradientVerticalOp<T2,T2>().filter(src1Smooth.getChannel(0),gradTemp.getChannel(1));
+               //DistanceOp<T2,T2>("5.0,2").filter(grad,gradTemp,w);
+               DistanceOp<T2,T2> gradQ;
+               gradQ.setParameter("mapping","l"); // linear
+               const double gradPow = this->getParameter("gradPow",2.0); // TODO always 2 => skip
+               gradQ.setParameter("l",gradPow);
+               gradQ.setParameter("lInv",1.0/gradPow);
+               gradQ.filter(grad,gradTemp,w);
+
+               FuzzyPeakOp<T2,T2> peak;
+               peak.setParameter("location",0);
+               const double gradWidth = this->getParameter("gradWidth",16.0);
+               peak.setParameter("width",gradWidth);
+               peak.setParameter("scaleDst",255.0); // 255
+               peak.filter(w,w);
+
+
+               if (drain::Debug > 3){
+                       Image<T> tmp;
+                       ImageView<T2> view;
+
+                       view.viewImageFlat(grad); //
+                       ScaleOp<T2,T>("4.0,128").filter(view,tmp);
+                       File::write(tmp,"grad-raw1.png");
+                       view.viewImageFlat(gradTemp);
+                       ScaleOp<T2,T>("4.0,128").filter(view,tmp);
+                       File::write(tmp,"grad-raw2.png");
+                       //ScaleOp<T2,T>("1.0,0").filter(w,tmp);
+                       CopyOp<T2,T>().filter(w,tmp);
+                       File::write(tmp,"grad-stab.png");
+               }
+
+               // Intensity gradients should not be computed for the 1st or 2nd image, but "between" them.
+               // So images will be mixed here. To save memory, src2Smooth is recycled.
+               MixerOp<T2,T2>(0.5).filter(src2Smooth,src1Smooth,src2Smooth);
+               GradientHorizontalOp<T2,T2>().filter(src2Smooth.getChannel(0),dx);
+               GradientVerticalOp<T2,T2>().filter(src2Smooth.getChannel(0),dy);
+
+               // Update Gradient stability as well (emphasize strong gradients)
+               Image<T2> gradMagnitude;
+               MagnitudeOp<T2,T2>("l,2").filter(grad,gradMagnitude);
+               MultiplicationOp<T2,T2>("255").filter(w,gradMagnitude,w);
+
+               //MagnitudeOp<T2,T2>("l,2").filter(grad,w);
+               //FuzzyPeakOp<T2,T2>(0,width,8.0).filter(w,w);
+
+               if (drain::Debug > 3){
+                       Image<T> tmp;
+                       CopyOp<T2,T>().filter(src2Smooth,tmp); File::write(tmp,"grad-smix.png");
+                       ScaleOp<T2,T>(0.5,128).filter(dt,tmp); File::write(tmp,"grad-dt.png");
+                       ScaleOp<T2,T>(1.5,128).filter(dx,tmp); File::write(tmp,"grad-dx.png");
+                       ScaleOp<T2,T>(1.5,128).filter(dy,tmp); File::write(tmp,"grad-dy.png");
+                       ScaleOp<T2,T>(1.0,0).filter(gradMagnitude,tmp); File::write(tmp,"grad-magnitude.png");
+                       ScaleOp<T2,T>(1.0,0).filter(w,tmp); File::write(tmp,"grad-quality.png");
+               }
+
+
+       }
+
+       ///  Cre
+       /**
+        *   Notice T2, T2
+        *   @param src - difference image (dx,dy,dt,q)
+        *   @param dst - motion field (u,v,q);
+        */
+       virtual void filter(const Image<T2> &src,Image<T2> &dst) const {
+
+               const unsigned int width  = this->getParameter("width",3);
+               const unsigned int height = this->getParameter("height",width);
+               this->opticalFlowWindow.setSize(width,height);
+               this->opticalFlowWindow.setDiff(src);
+
+               SlidingWindowOp<T2,T2>::filter(src,dst);
+
+               if (drain::Debug > 3){
+                       Image<T> tmp;
+                       ScaleOp<T2,T>(2,128.0).filter(dst,tmp);
+                       File::write(tmp,"oflow-M.png");
+                       ImageView<T> tmpView;
+                       tmpView.viewImageFlat(tmp);
+                       File::write(tmpView,"oflow-MF.png");
+                       for (unsigned int i = 0; i < dst.getChannelCount(); ++i) {
+                               stringstream filename;
+                               filename << "oflow-M" << i << ".png";
+                               ScaleOp<T2,T>(0.5,128.0).filter(dst.getChannel(i),tmp);
+                               File::write(tmp,filename.str());
+                       };
+               }
+
+       };
+
+protected:
+
+       mutable SlidingOpticalFlow<T2> opticalFlowWindow;
+
+
+};
+
+}
+}
+
+#endif
diff --git a/drain/image/File.h b/drain/image/File.h
new file mode 100644 (file)
index 0000000..991da9c
--- /dev/null
@@ -0,0 +1,93 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef DRAIN_FILE_H_
+#define DRAIN_FILE_H_
+
+#include "MagickDrain.h"
+#include "FilePng.h"
+
+#include <string>
+
+#include "Image.h"
+
+
+namespace drain
+{
+namespace image
+{
+
+using namespace std;
+
+class File
+{
+public:
+       //((File(const string &filename);
+
+       //virtual ~File();
+
+       template <class T>
+       static void read(Image<T> &img,const string &path){
+#ifdef DRAIN_MAGICK_yes
+               Magick::Image magickImage;
+               magickImage.read(path);
+               MagickDrain::convert(magickImage,img);
+#else
+               // Todo PNM support
+               FilePng::read(img,path);
+#endif
+       };
+
+       //static void read(Image<unsigned char> &image,const string &path);
+
+       template <class T>
+       static void write(const Image<T> &img,const string &path){
+#ifdef DRAIN_MAGICK_yes
+
+               Magick::Image magickImage;
+               magickImage.size("1x1");
+               magickImage.magick("RGBA");
+
+               if (Debug > 1)
+                       cerr << "Converting image to Magick image." << endl;
+
+               MagickDrain::convert(img,magickImage);
+
+               if (Debug > 0)
+                       cerr << "Writing image" << img.getName() << " to "<< path << endl;
+
+               magickImage.write(path);
+#else
+               FilePng::write(img,path);
+#endif
+       };
+
+       //static void write(const Image<unsigned char> &image,const string &path);
+
+       
+};
+
+
+}  // image
+
+}  // drain
+
+#endif /*FILE_H_*/
diff --git a/drain/image/FileBinary.cpp b/drain/image/FileBinary.cpp
new file mode 100644 (file)
index 0000000..4cc1e5a
--- /dev/null
@@ -0,0 +1,56 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#include "FileBinary.h"
+
+
+#include <fstream>
+#include <iostream>
+#include <istream>
+
+#include <sstream>
+
+//#include <fstream>
+
+
+
+namespace drain
+{
+
+namespace image
+{
+
+/*
+File::File(const string &filename)
+{
+       this->filename = filename;
+}
+*/
+
+
+
+}
+
+
+
+}
+
+
diff --git a/drain/image/FileBinary.h b/drain/image/FileBinary.h
new file mode 100644 (file)
index 0000000..7baf348
--- /dev/null
@@ -0,0 +1,113 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef FILEBINARY_H_
+#define FILEBINARY_H_
+
+
+#include <string>
+#include <iostream>
+#include <fstream>
+
+
+
+#include "Image.h"
+
+
+namespace drain
+{
+namespace image
+{
+
+using namespace std;
+
+/** Format applied at FMI
+ * 
+ */
+class FileBinary
+{
+public:
+       //((FileBinary(const string &FileBinaryname);
+
+       virtual ~FileBinary();
+
+       /*
+       template <class T>
+       static void read(BufferedImage<T> &image,const string &path);
+
+       static void read(BufferedImage<unsigned char> &image,const string &path);
+       */
+
+       template <class T>
+       static void write(const Image<T> &image,const string &path, const string &header = "");
+
+       /*
+       static void write(const BufferedImage<unsigned char> &image,const string &path);
+       */
+       
+protected:
+       //string filename;
+       
+       
+};
+
+
+// Consider template <class T2> for on-the-fly cast T=>T2
+// Consider all permutations WHD dwH  
+template <class T>
+void FileBinary::write(const Image<T> &image,const string &path, const string &header){
+       ofstream ofstr;
+       ofstr.open(path.c_str(),ios::out);
+       
+       ofstr << header;
+       
+       const vector<T> & v = image.getVector();
+       const typename vector<T>::size_type s = v.size();
+       // for (int k=0; k<image.getChannelCount(); k++)
+       // for (int j=0; j<image.getChannelCount(); j++)
+       
+       for (typename vector<T>::size_type i=0; i<s;i++){
+               ofstr << v[i]; 
+       }
+       
+       
+       /*
+       string timestamp = image.parameters.get("TIMESTAMP","197001010000");
+       timestamp.append(14-timestamp.size(), '0');  // add seconds if needed
+       
+       string levels = image.parameters.get("LEVELS","197001010000");
+       string site = image.parameters.get("SITE","197001010000");
+
+105 300m 300
+dBz 123
+vantaantutka 1201
+567 623
+latlon:19.5,57.4,23.7,62.5
+4567456
+       */
+
+}
+
+}
+
+}
+
+#endif /*FILEBINARY_H_*/
diff --git a/drain/image/FileFQD.cpp b/drain/image/FileFQD.cpp
new file mode 100644 (file)
index 0000000..9afe084
--- /dev/null
@@ -0,0 +1,56 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#include "FileFQD.h"
+
+
+#include <fstream>
+#include <iostream>
+#include <istream>
+
+#include <sstream>
+
+//#include <fstream>
+
+
+
+namespace drain
+{
+
+namespace image
+{
+
+/*
+File::File(const string &filename)
+{
+       this->filename = filename;
+}
+*/
+
+
+
+}
+
+
+
+}
+
+
diff --git a/drain/image/FileFQD.h b/drain/image/FileFQD.h
new file mode 100644 (file)
index 0000000..5073db9
--- /dev/null
@@ -0,0 +1,95 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef FILEFQD_H_
+#define FILE_H_
+
+
+#include <string>
+
+#include "Image.h"
+
+
+namespace drain
+{
+namespace image
+{
+
+using namespace std;
+
+/** Format applied at FMI
+ * 
+ */
+class FileFQD
+{
+public:
+       //((FileFQD(const string &FileFQDname);
+
+       virtual ~FileFQD();
+
+       /*
+       template <class T>
+       static void read(BufferedImage<T> &image,const string &path);
+
+       static void read(BufferedImage<unsigned char> &image,const string &path);
+       */
+
+       template <class T>
+       static void write(const Image<T> &image,const string &path);
+
+       /*
+       static void write(const BufferedImage<unsigned char> &image,const string &path);
+       */
+       
+protected:
+       //string filename;
+       
+       
+};
+
+
+template <class T>
+void FileFQD::write(const Image<T> &image,const string &path){
+       ofstream ofstr;
+       //ofstr.open(path,"w");
+       
+       // 
+       string timestamp = image.parameters.get("TIMESTAMP","197001010000");
+       timestamp.append(14-timestamp.size(), '0');  // add seconds if needed
+       
+       string levels = image.parameters.get("LEVELS","197001010000");
+       string site = image.parameters.get("SITE","197001010000");
+       /*
+105 300m 300
+dBz 123
+vantaantutka 1201
+567 623
+latlon:19.5,57.4,23.7,62.5
+4567456
+       */
+
+}
+
+}
+
+}
+
+#endif /*FILEFQD_H_*/
diff --git a/drain/image/FilePng.h b/drain/image/FilePng.h
new file mode 100644 (file)
index 0000000..fa624d8
--- /dev/null
@@ -0,0 +1,348 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy  of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef FILEPNG_H_
+#define FILEPNG_H_
+
+
+#include <string>
+#include <iostream>
+#include <fstream>
+#include <exception>
+
+#include <png.h>
+
+#include "Image.h"
+#include "../util/Time.h"
+
+
+namespace drain
+{
+namespace image
+{
+
+using namespace std;
+
+/// For reading and writing images in PNG format
+/** A lighter alternative for Magick++ which supports tens of image formats.
+ *  Portable Network Graphics (PNG) format is a widely adapted, flexible image format
+ *  for gray level and color images. It supports lossless compression and alpha channels.
+ */
+class FilePng
+{
+public:
+
+       template <class T>
+       static void read(Image<T> &image,const string &path,
+                       int png_transforms = (PNG_TRANSFORM_PACKING || PNG_TRANSFORM_EXPAND));
+
+               template <class T>
+       static void write(const Image<T> &image,const string &path);
+       /// TODO: Write options (PNG_*)
+       
+};
+
+/** Reads a png file to drain::Image.
+ *  Converts indexed (palette) images to RGB or RGBA.
+ *  Scales data to 8 or 16 bits, according to template class.
+ *  Floating point images will be scaled as 16 bit integral (unsigned short int).
+ *
+ */
+template <class T>
+void FilePng::read(Image<T> &image,const string &path, int png_transforms ) {
+
+       // Try to open the file
+       FILE *fp = fopen(path.c_str(), "rb");
+       if (fp == NULL)
+               throw runtime_error(string("FilePng: could not open file: ") + path);
+
+       // For checking magic code (signature)
+       const unsigned int PNG_BYTES_TO_CHECK=4;
+       png_byte buf[PNG_BYTES_TO_CHECK];
+
+       /* Read in some of the signature bytes */
+       if (fread(buf, 1, PNG_BYTES_TO_CHECK, fp) != PNG_BYTES_TO_CHECK)
+               throw runtime_error(string("FilePng: suspicious size of file: ") + path);
+
+       /* Compare the first PNG_BYTES_TO_CHECK bytes of the signature.
+          Return nonzero (true) if they match */
+       if (png_sig_cmp(buf, (png_size_t)0, PNG_BYTES_TO_CHECK) != 0)
+               throw runtime_error(string("FilePng: not a png file: ")+path);
+
+       png_structp  png_ptr = png_create_read_struct(PNG_LIBPNG_VER_STRING,NULL,NULL,NULL);
+       if (!png_ptr){
+               throw runtime_error(string("FilePng: problem in allocating image memory for: ")+path);
+       }
+
+       png_infop info_ptr = png_create_info_struct(png_ptr);
+       if (!info_ptr){
+           png_destroy_read_struct(&png_ptr,(png_infopp)NULL, (png_infopp)NULL);
+               throw runtime_error(string("FilePng: problem in allocating info memory for: ")+path);
+       }
+
+       /*
+       png_infop end_info = png_create_info_struct(png_ptr);
+       if (!end_info){
+           png_destroy_read_struct(&png_ptr, &info_ptr,(png_infopp)NULL);
+           throw runtime_error(string("FilePng: problem in allocating end_info memory for: ")+path);
+       }
+       */
+
+       // This may be unstable. According to the documentation, if one uses the high-level interface png_read_png()
+       // one can only configure it with png_transforms flags (PNG_TRANSFORM_*)
+       png_set_palette_to_rgb(png_ptr);
+
+       png_init_io(png_ptr, fp);
+       png_set_sig_bytes(png_ptr, PNG_BYTES_TO_CHECK);
+
+       /// Main action
+       if (drain::Debug > 2)
+               cerr << "png read starting\n";
+       png_read_png(png_ptr, info_ptr, png_transforms, NULL);
+
+
+       /// Read comments
+       if (drain::Debug > 2)
+               cerr << "png image comments:\n";
+
+       int num_text;
+       png_textp text_ptr;
+       png_get_text(png_ptr, info_ptr,&text_ptr, &num_text);
+       for (int i = 0; i < num_text; ++i) {
+               if (drain::Debug > 2)
+                       cerr << text_ptr[i].key << '=' << text_ptr[i].text << '\n';
+               image.properties[text_ptr[i].key] = text_ptr[i].text;
+       }
+
+
+       /// Copy to drain::Image
+       const int width  = png_get_image_width(png_ptr, info_ptr);
+       const int height = png_get_image_height(png_ptr, info_ptr);
+       const int channels = png_get_channels(png_ptr, info_ptr);
+       const int bit_depth = png_get_bit_depth(png_ptr, info_ptr);
+
+       switch (channels){
+       case 4:
+               image.setGeometry(width,height,3,1);
+               break;
+       case 3:
+               image.setGeometry(width,height,3);
+               break;
+       case 2:
+               image.setGeometry(width,height,1,1);
+               break;
+       case 1:
+               image.setGeometry(width,height,1);
+               break;
+       default:
+               throw runtime_error(string("FilePng: invalid channel count in : ")+path);
+       }
+
+       if (drain::Debug > 2){
+               cerr << "png geometry ok\n";
+               cerr << "png channels =" << channels << "\n";
+               cerr << "png bit_depth=" << bit_depth << "\n";
+               cerr << image.getGeometry() << "\n";
+       }
+
+       if ((bit_depth!=8) && (bit_depth != 16)){
+               fclose(fp);
+               png_destroy_read_struct(&png_ptr,&info_ptr, (png_infopp)NULL);
+               //png_free_data(png_ptr,info_ptr,PNG_FREE_ALL,-1);  // ???
+               throw runtime_error(string("FilePng: unsupported bit depth in : ")+path);
+       }
+
+       png_bytep *row_pointers = png_get_rows(png_ptr, info_ptr);
+       png_bytep p;
+       int i0;
+       for (int j = 0; j < height; ++j) {
+               p = row_pointers[j];
+               for (int i = 0; i < width; ++i) {
+                       for (int k = 0; k < channels; ++k) {
+                               i0 = channels*i + k;
+                               if (bit_depth == 8) {
+                                       image.at(i,j,k) = p[i0];
+                               }
+                               else {
+                                       image.at(i,j,k) = p[i0*2] + (p[i0*2+1]<<8);
+                               }
+                       }
+               }
+       }
+
+       fclose(fp);
+       png_destroy_read_struct(&png_ptr,&info_ptr, (png_infopp)NULL);
+       //png_free_data(png_ptr,info_ptr,PNG_FREE_ALL,-1);  // ???
+
+       //png_destroy_read_struct(&png_ptr,(png_infopp)NULL, (png_infopp)NULL);
+       //png_destroy_info_struct(png_ptr,&info_ptr);
+
+
+}
+
+/** Writes drain::Image to a png image file applying G,GA, RGB or RGBA color model.
+ *  Writes in 8 or 16 bits, according to template class.
+ *  Floating point images will be scaled as 16 bit integral (unsigned short int).
+*/
+template <class T>
+void FilePng::write(const Image<T> &image,const string &path){
+
+
+       FILE *fp = fopen(path.c_str(), "wb");
+       if (!fp){
+               throw runtime_error(string("FilePng: could not open file : ")+path);
+       }
+
+       png_structp png_ptr = png_create_write_struct(PNG_LIBPNG_VER_STRING, NULL, NULL, NULL);
+       if (!png_ptr){
+               throw runtime_error(string("FilePng: could not allocate memory for: ")+path);
+       }
+
+       png_infop info_ptr = png_create_info_struct(png_ptr);
+       if (!info_ptr){
+          png_destroy_write_struct(&png_ptr,(png_infopp)NULL);
+          throw runtime_error(string("FilePng: could not allocate info memory for: ")+path);
+       }
+
+       // Required.
+       png_init_io(png_ptr, fp);
+
+       // Optional.
+       png_set_filter(png_ptr, 0,PNG_FILTER_HEURISTIC_DEFAULT);
+
+       // Optional
+       png_set_compression_level(png_ptr, Z_DEFAULT_COMPRESSION);
+
+       const int width  = image.getWidth();
+       const int height = image.getHeight();
+       const int channels = image.getChannelCount();
+
+       int color_type;
+       switch (channels) {
+               case 4:
+                       color_type = PNG_COLOR_TYPE_RGBA;
+                       break;
+               case 3:
+                       color_type = PNG_COLOR_TYPE_RGB;
+                       break;
+               case 2:
+                       color_type = PNG_COLOR_TYPE_GRAY_ALPHA;
+                       break;
+               case 1:
+                       color_type = PNG_COLOR_TYPE_GRAY;
+                       break;
+               default:
+                       throw runtime_error(string("FilePng: unsupported channel count: "));
+       }
+
+
+       const int byte_depth = sizeof(T);
+       const int bit_depth  = byte_depth <= 2 ? byte_depth*8 : 16;
+       if (drain::Debug > 1){
+               if (byte_depth > 2)
+                       cerr << "FilePng: source image " << byte_depth << " byte data, converting to 2 bytes (16bit uInt).\n";
+       }
+       // TODO: check integral vs float types, instead
+
+       if (drain::Debug > 2){
+               cerr << image.getGeometry() << "\n";
+               cerr << "png bit_depth=" << bit_depth << "\n";
+       }
+
+       // Set header information
+       png_set_IHDR(png_ptr, info_ptr, width, height,bit_depth, color_type,
+                       PNG_INTERLACE_NONE, PNG_COMPRESSION_TYPE_DEFAULT, PNG_FILTER_TYPE_DEFAULT);
+    // interlace_type, compression_type, filter_method
+
+       // Optional
+       // png_set_tIME(png_ptr, info_ptr, mod_time);
+       // png_set_text(png_ptr, info_ptr, text_ptr, num_text);
+
+       // Optional: text content
+       const int text_size = image.properties.size() + 2;
+       png_text text_ptr[text_size];
+       drain::Time time;
+       text_ptr[0].key = (char *)"Creation time";
+       text_ptr[0].text = (char *)time.str().c_str();
+       text_ptr[0].text_length = strlen(text_ptr[0].text); // ???
+       text_ptr[0].compression = PNG_TEXT_COMPRESSION_NONE;
+       //text_ptr[1].key = (char *)"Comment";
+       //text_ptr[1].text = (char *)"Test comment";
+       text_ptr[1].key = (char *)"Software";
+       text_ptr[1].text = (char *)"drain/image/FilePng Markus.Peura[c]fmi.fi";
+       text_ptr[1].text_length = strlen(text_ptr[1].text);
+       text_ptr[1].compression = PNG_TEXT_COMPRESSION_NONE;
+
+       int i=2;
+       for (map<string,Data>::const_iterator it = image.properties.begin(); it != image.properties.end(); it++){
+               text_ptr[i].key  = (char *)it->first.c_str();
+               text_ptr[i].text = (char *)(it->second.substr(0,79)).c_str();
+               text_ptr[i].text_length = it->second.length();
+               text_ptr[i].compression = PNG_TEXT_COMPRESSION_NONE;
+               i++;
+       }
+       png_set_text(png_ptr, info_ptr, text_ptr, 2);
+
+       // Create temporary image array.
+       png_byte data[height][width*channels*byte_depth];
+       for (int k = 0; k < channels; ++k) {
+               // 8 bit
+               if (byte_depth == 1){
+                       for (int i = 0; i < width; ++i) {
+                               //i0 = i*channels + k;
+                               for (int j = 0; j < height; ++j) {
+                                       data[j][i*channels + k] = static_cast<png_byte>(image.at(i,j,k));
+                               }
+                       }
+               }
+               // 16 bit
+               else {
+                       int i0;
+                       for (int i = 0; i < width; ++i) {
+                               i0 = (i*channels + k)*2;
+                               for (int j = 0; j < height; ++j) {
+                                       data[j][i0+1] = static_cast<png_byte>(image.at(i,j,k)&255);
+                                       data[j][i0  ] = static_cast<png_byte>((image.at(i,j,k)>>8)&255);
+                               }
+                       }
+               }
+       }
+       png_bytep row_pointers[height];
+       for (int j = 0; j < height; ++j) {
+               row_pointers[j] = data[j];
+       }
+       //cerr << "Preparing to write..." << endl;
+       png_set_rows(png_ptr, info_ptr, row_pointers);
+
+       // Main operation
+       int png_transforms = PNG_TRANSFORM_IDENTITY  || PNG_TRANSFORM_SHIFT;
+       png_write_png(png_ptr, info_ptr, png_transforms, NULL);
+
+       fclose(fp);
+       png_destroy_write_struct(&png_ptr,&info_ptr);
+       //png_free_data(png_ptr,info_ptr,PNG_FREE_ALL,-1);
+}
+
+}
+
+}
+
+#endif /*FILEPng_H_*/
diff --git a/drain/image/FloodFillOp.h b/drain/image/FloodFillOp.h
new file mode 100644 (file)
index 0000000..c3000f5
--- /dev/null
@@ -0,0 +1,148 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef FLOODFILL2_H
+#define FLOODFILL2_H
+
+#include "ImageOp.h"
+#include "CoordinateHandler.h"
+
+
+namespace drain
+{
+namespace image
+{
+       
+//template <class T = unsigned char, class T2 = unsigned char>
+//class FloodFillRecursion;  // see bottom of this file
+
+// The actual operator is defined at the bottom of this file.
+
+
+/// Helper class (not an operator itself) applied by FloodFillOp and SegmentAreaOp
+/**
+ *   \author Markus.Peura@fmi.fi
+ */
+template <class T = unsigned char, class T2 = unsigned char>
+//template <class T, class T2>
+class FloodFillRecursion {
+  public: 
+        FloodFillRecursion(const Image<T> &s,Image<T2> &d) : src(s), dst(d), h(s.getCoordinateHandler()) {};
+        
+        /// Fills the segment having constant intensity, that is, src.at(i0,j0).
+        unsigned long int fill(unsigned int i, unsigned int j, T2 fillValue){
+               dst.setGeometry(src.getWidth(),src.getHeight(),
+                       max(1u,dst.getImageChannelCount()),dst.getAlphaChannelCount()); 
+               return fill(i,j,fillValue,src.at(i,j),src.at(i,j));
+         };
+         
+         /// Fills the segment having intensity between min and max.
+         unsigned long int fill(unsigned int i, unsigned int j, T2 fillValue, T min, T max = Intensity::max<T>()){
+               value = fillValue;
+           anchorMin = min;
+               anchorMax = max;
+           size = 0;
+          // cout << "range:" << (int)anchorMin  << '-' << (int)anchorMax  << '\n';
+           fill8(i,j);
+           //cout << "size:" << size  << '\n';
+           return size;
+         };
+         
+  protected:
+        const Image<T> &src;
+        Image<T2> &dst;
+        CoordinateHandler &h;
+       mutable T  anchorMin;
+       mutable T  anchorMax;
+       mutable T2 value;
+       mutable Point2D<> p;
+       mutable unsigned long int size;
+       
+    void fill8(unsigned int i,unsigned int j){
+    
+       
+       // Is outside segment?
+       if (src.at(i,j) < anchorMin)
+               return;
+       if (src.at(i,j) > anchorMax)
+               return;
+       
+       // Is visited?
+       if (dst.at(i,j) == value)
+               return;
+       
+       
+       // MAIN ACTION:
+       dst.at(i,j) = value;
+       size++;
+       
+       p.setLocation(i-1,j);
+       if ((h.handle(p) & CoordinateHandler::IRREVERSIBLE) == 0)
+               fill8(p.x,p.y);
+       p.setLocation(i,j-1);   
+       if ((h.handle(p) & CoordinateHandler::IRREVERSIBLE) == 0)
+           fill8(p.x,p.y);
+       p.setLocation(i+1,j);   
+       if ((h.handle(p) & CoordinateHandler::IRREVERSIBLE) == 0)
+           fill8(p.x,p.y);
+       p.setLocation(i,j+1);   
+       if ((h.handle(p) & CoordinateHandler::IRREVERSIBLE) == 0)
+           fill8(p.x,p.y);
+    }
+
+};
+
+
+/**
+ *   \author Markus.Peura@fmi.fi
+ */
+template <class T = unsigned char, class T2 = unsigned char>
+class FloodFillOp : public ImageOp<T,T2> {
+       
+public:
+       FloodFillOp(const string & p = "1,1,128,255,255"){
+       this->setInfo("Fills an area starting at (i,j) having intensity in [min,max], with a value.",
+               "i,j,min,max,value",p);
+       };
+    
+    virtual void makeCompatible(const Image<T> &src,Image<T2> &dst) const  {
+               dst.setGeometry(src.getWidth(),src.getHeight(),
+                               max(1u,dst.getImageChannelCount()),dst.getAlphaChannelCount()); 
+       }
+       
+    virtual void filter(const Image<T> &src,Image<T2> &dst) const {
+                       makeCompatible(src,dst);
+                       const int i = this->parameters.get("i",0);
+                       const int j = this->parameters.get("j",0);
+                       const T min = this->parameters.get("min",128);
+                       const T max = this->parameters.get("max",(int)Intensity::max<T>());
+                       const T2 value = this->parameters.get("value",(int)Intensity::max<T2>());
+                       FloodFillRecursion<> f(src,dst);
+                       f.fill(i,j,value,min,max);
+       };
+       
+};
+
+
+}
+}
+#endif /* FLOODFILL_H_ */
+
diff --git a/drain/image/FuzzyPeakOp.h b/drain/image/FuzzyPeakOp.h
new file mode 100644 (file)
index 0000000..53a2892
--- /dev/null
@@ -0,0 +1,93 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef FUZZYPEAKOP_H_
+#define FUZZYPEAKOP_H_
+
+#include <math.h>
+
+
+#include "SequentialImageOp.h"
+
+using namespace std;
+
+namespace drain
+{
+
+namespace image
+{
+
+/// FuzzyPeak correction: intensitys is mapped as f' = f^(gamma)
+/*! FuzzyPeak correction
+ * 
+ *
+ */
+template <class T = unsigned char,class T2 = unsigned char>
+class FuzzyPeakOp : public SequentialImageOp<T,T2>
+{
+public:
+
+       //FuzzyPeakOp(const Data &option);
+       FuzzyPeakOp(const string & p = "0,128,255") :
+               SequentialImageOp<T,T2>("FuzzyPeakOp","Gaussian-like peak.","location,width,scaleDst",p){
+               //this->setInfo("Fuzzy Gaussian-like peak.","location,width,scaleDst",p);
+       }
+       
+       FuzzyPeakOp(double location,double width, double scaleDst=255.0) :
+               SequentialImageOp<T,T2>("FuzzyPeakOp","Gaussian-like peak.","location,width,scaleDst",""){
+               this->setParameter("location",location);
+               this->setParameter("width",width);
+               this->setParameter("scaleDst",scaleDst);
+               //this->setInfo("Fuzzy Gaussian-like peak.","location,width,scaleDst",p);
+       }
+
+       /// - scale of source image -
+       //FuzzyPeakOp(double width,double location){
+       void initialize() const {
+               double w = this->parameters.get("width",Intensity::max<T2>()/2.0);
+               width2 = w*w;
+               location = this->parameters.get("location",0);
+               scaleDst = this->parameters.get("scaleDst",static_cast<double>(Intensity::max<T2>()));
+               if (drain::Debug > 3)
+                       cerr << "FuzzyPeak: " << this->location << ',' << this->width2 << ',' << this->scaleDst << '\n';
+       };
+       
+       inline 
+       void filterValue(const T &src, T2 &dst) const {
+               static double x;
+               x = static_cast<double>(src) - this->location ;
+               
+               dst = static_cast<T2>(this->scaleDst*width2 / (width2 + x*x));
+                 
+       }; 
+       
+protected:
+
+       mutable double width2;
+       mutable double location;
+       mutable double scaleDst;
+       
+};
+
+}
+}
+
+#endif /*GAMMAOP_H_*/
diff --git a/drain/image/FuzzyThresholdOp.h b/drain/image/FuzzyThresholdOp.h
new file mode 100644 (file)
index 0000000..c2ab121
--- /dev/null
@@ -0,0 +1,93 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef FUZZYTHRESHOLDOP_H_
+#define FUZZYTHRESHOLDOP_H_
+
+#include <math.h>
+
+
+//#include "SequentialImageOp.h"
+#include "ThresholdOp.h"
+
+using namespace std;
+
+namespace drain
+{
+
+namespace image
+{
+
+/// FuzzyThreshold correction: intensity is mapped as f' = f^(gamma)
+/*! FuzzyThreshold correction
+ * 
+ *   @see
+ */
+template <class T = unsigned char,class T2 = unsigned char>
+class FuzzyThresholdOp : public SequentialImageOp<T,T2>
+       // public ThresholdOp<T,T2>
+{
+public:
+
+       /// Future versions may use dstMin, dstMax, srcMin, srcMax
+       FuzzyThresholdOp(const string & p = "0,1,0,255") :
+               SequentialImageOp<T,T2>("FuzzyThresholdOp","Smooth transition.","location,width,dstMin,dstMax",p){
+       }
+
+       /// - scale of source image -
+       void initialize() const {
+               //cout << this->parameters << '\n';
+               location = this->parameters.get("location",0);
+               width = this->parameters.get("width",(double)Intensity::max<T2>()/2.0);
+               float min = this->parameters.get("dstMin",(double)Intensity::min<T2>());
+               float max = this->parameters.get("dstMax",(double)Intensity::max<T2>());
+               scaleDstHalf  = (min+max)/2.0;
+               scaleDstCoeff = (max-min)/2.0;
+               if (drain::Debug > 3){
+                       cerr << "FuzzyTheshold: " << this->location << ',' << this->width << ',';
+                       cerr << " " << scaleDstHalf << ' ' << scaleDstCoeff << '\n';
+               }
+       };
+
+       inline 
+       void filterValue(const T &src, T2 &dst) const {
+               double x = static_cast<double>(src) - location;
+               if (x>0){
+                       dst = static_cast<T2>(scaleDstHalf + (scaleDstCoeff * x)/(width + x));  
+               }
+               else {
+                       dst = static_cast<T2>(scaleDstHalf + (scaleDstCoeff * x)/(width - x));  
+               }
+       }; 
+
+protected:
+
+       mutable double location;
+       mutable double width;
+       mutable double scaleDstHalf;
+       mutable double scaleDstCoeff;
+
+};
+
+}
+}
+
+#endif /*GAMMAOP_H_*/
diff --git a/drain/image/GammaOp.h b/drain/image/GammaOp.h
new file mode 100644 (file)
index 0000000..858fb39
--- /dev/null
@@ -0,0 +1,78 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef GAMMAOP_H_
+#define GAMMAOP_H_
+
+#include <math.h>
+
+//#include <limits>
+
+//#include "BufferedImage.h"
+#include "SequentialImageOp.h"
+
+using namespace std;
+
+namespace drain
+{
+
+namespace image
+{
+
+/// Gamma correction: intensitys is mapped as f' = f^(gamma)
+/*! Gamma correction
+ * 
+ *  NOTE. Design for paramters may vary, since multichannel could be handled by giving
+ *  a value for each: 1.2,1.4,0.7 for example. 
+ */
+template <class T = unsigned char,class T2 = unsigned char>
+class GammaOp : public SequentialImageOp<T,T2>
+{
+public:
+
+       GammaOp(const string & p = "1.0") : SequentialImageOp<T,T2>("GammaOp",
+                       "Gamma correction for brightness.","gamma,scaleSrc,scaleDst","1,255,255"){
+               this->setParameters(p);
+       };
+
+       virtual void initialize() const {
+               gammaInv = 1.0 / this->parameters.get("gamma",1.0);
+               scaleSrc = this->parameters.get("scaleSrc",Intensity::max<T>());
+               scaleDst = this->parameters.get("scaleDst",Intensity::max<T2>());
+       };
+
+       
+       inline 
+       void filterValue(const T &src, T2 &dst) const {
+               dst = static_cast<T2>(scaleDst * pow(static_cast<double>(src)/scaleSrc, gammaInv));
+       }; 
+       
+protected:
+       mutable double gammaInv;
+       mutable double scaleSrc;
+       mutable double scaleDst;
+       
+};
+
+}
+}
+
+#endif /*GAMMAOP_H_*/
diff --git a/drain/image/Geometry.cpp b/drain/image/Geometry.cpp
new file mode 100644 (file)
index 0000000..e38892d
--- /dev/null
@@ -0,0 +1,157 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+
+#include <algorithm>
+#include <sstream>
+
+#include "Geometry.h"
+//#include "DistanceTransformOp.h"
+
+namespace drain
+{
+
+namespace image
+{
+       
+using namespace std;
+       
+// ref ei onnistunut (width, height, chC)
+Geometry::Geometry() : 
+//     std::vector<int>(3), width(at(0)), height(at(1)), channelCount(at(2)), alphaChannelCount(0) 
+       std::vector<unsigned int>(3)
+{
+       setGeometry(0,0,0,0);
+}
+
+Geometry::Geometry(unsigned int width, unsigned int height, unsigned int imageChannelCount, unsigned int alphaChannelCount) :
+       vector<unsigned int>(max(3u,imageChannelCount+alphaChannelCount))
+       //vector<int>(max(3,channelCount)), width(at(0)), height(at(1)), channelCount(at(2)), alphaChannelCount(0)
+{
+       setGeometry(width,height,imageChannelCount,alphaChannelCount);
+       //DistanceTransformLinearOp<> foo;
+}
+
+Geometry::~Geometry()
+{
+}
+
+void Geometry::setGeometry(const Geometry &g){
+       setGeometry(g.getWidth(),g.getHeight(),g.getImageChannelCount(),g.getAlphaChannelCount());
+} 
+
+void Geometry::setGeometry(unsigned int width, unsigned int height, unsigned int imageChannelCount, unsigned int alphaChannelCount){
+  resize(std::max(static_cast<size_type>(3),size()));
+       
+       at(0) = width;
+       at(1) = height;
+       at(2) = imageChannelCount + alphaChannelCount;
+       this->alphaChannelCount = alphaChannelCount;
+       update();
+       //area = width * height;
+       //volume = area*(imageChannelCount + alphaChannelCount);
+} 
+       
+void Geometry::setWidth(unsigned int w){
+       //resize(std::max(3u,size()));
+       
+       at(0) = w;
+       update();
+};
+
+void Geometry::setHeight(unsigned int h){
+       //resize(std::max(3u,size()));
+       
+       at(1) = h;
+       update(); 
+};
+
+void Geometry::setChannelCount(unsigned int imageChannels, unsigned int alphaChannels){
+       
+       //resize(std::max(3u,size()));
+       
+       at(2) = imageChannels + alphaChannels;
+       alphaChannelCount = alphaChannels;
+       update();
+};
+
+
+void Geometry::setAlphaChannelCount(unsigned int alphaChannels){
+       setChannelCount(getImageChannelCount(),alphaChannels);
+};
+
+
+void Geometry::update(){
+       
+       //resize(std::max(3u,size()));
+       
+       width = at(0);
+       height = at(1);
+       channelCount = at(2); 
+       
+       imageChannelCount = channelCount - alphaChannelCount; 
+       area = width*height;
+       volume = 1;
+       for (unsigned int i = 0; i < this->size(); ++i) {
+               volume *= (*this)[i];
+       }
+}
+
+
+
+Geometry & Geometry::operator=(const Geometry &g)
+{
+       (*this) = (const std::vector<int> &)g;
+       return (*this); 
+};
+    
+
+ostream & operator<<(ostream &ostr,const Geometry &geometry) {
+
+       string separator("");
+       
+       ostr << "[";
+    for (unsigned int i=0; i<geometry.size(); i++)
+    {
+       ostr << separator;
+        ostr << geometry[i];
+        separator = "×";
+    }
+    ostr << "] ";
+    ostr << "iC=" << geometry.getImageChannelCount() << ", aC=" <<  geometry.getAlphaChannelCount()  << ", ";
+    ostr << "w=" << geometry.getWidth() << ", h=" <<  geometry.getHeight()  << ", ";
+    ostr << "A=" << geometry.getArea() << ", V=" <<  geometry.getVolume()  << " ";
+               
+       return ostr;
+}
+string &Geometry::toString(string & s) const
+{
+    stringstream sstr;
+       sstr << *this;
+       s = sstr.str();
+    return s;
+};
+
+
+}
+
+}
diff --git a/drain/image/Geometry.h b/drain/image/Geometry.h
new file mode 100644 (file)
index 0000000..c2b76be
--- /dev/null
@@ -0,0 +1,148 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef GEOMETRY_H_
+#define GEOMETRY_H_
+
+#include <iostream>
+#include <vector>
+
+namespace drain
+{
+
+namespace image
+{
+
+using namespace std;
+
+/*! The basic idea is to encode dimensions directly as a vector<int>; the number of elements is 
+ *  the number of dimensions. Each element states the discrete coordinate space.
+ * 
+ *  Three standard "copies" are defined for convenience; KLUDGE!
+ *  # geometry.width  = geometry[0]
+ *  # geometry.height = geometry[1]
+ *  # geometry.channelCount = geometry[2]  
+ *
+ *  To guarantee this, the minimum size of the vector will always be at least three (3).
+ *  
+ * 
+ *  For example, a three-channel image having width 640 and height 400 would have dimension vector 
+ *  [640,400,3]. The coordinates will be: x \in {0,1,...,639}, y \in {0,1,...,399} and z \in {0,1,2}.
+ *  The size of the vector of unlimited, allowing hypermatrices of any size. 
+ * 
+ */
+class Geometry : public std::vector<unsigned int> 
+// SP: itse en käyttäisi perintää vaan kompositiota, ts. vector-olio on luokan jäsenenä. 
+// MP: OK, muunnan kun ehdin. Muistaakseni ei ollut erityistä syytä tälle suoralle.
+{
+public:
+
+    Geometry();
+    
+       /*! Width, height, channels, imageChannels (alphaChannels)
+        */
+    Geometry(unsigned int width, unsigned int height, unsigned int channelCount=1, unsigned int alphaChannelCount=0);
+       
+       //Geometry(const std::vector<int> &vector, int a = 0); 
+
+    virtual ~Geometry();
+
+       void setGeometry(const Geometry &g);
+    
+    void setGeometry(unsigned int width,unsigned int height,
+               unsigned int imageChannelCount = 1,unsigned int alphaChannelCount = 0);
+    
+    void setWidth(unsigned int weight);
+    void setHeight(unsigned int height);
+    void setChannelCount(unsigned int imageChannelCount,unsigned  int alphaChannelCount = 0);
+    void setAlphaChannelCount(unsigned  int alphaChannelCount);
+    
+
+
+    // SP: Miksi palautetaan referenssit? Eikö alkeistyypeillä riitä arvon palautus?
+    // Hyvä kommentti tämäkin. En muista. Liittyneekö siihen, että jos kuva elää
+    // kesken suoritusta, muuttujakin pysyy ajan tasalla.
+       inline const unsigned int & getWidth() const { return width; };
+       inline const unsigned int & getHeight() const { return height; };
+       inline const unsigned int & getChannelCount() const { return channelCount; };
+       inline const unsigned int & getImageChannelCount() const { return imageChannelCount; };
+       inline const unsigned int & getAlphaChannelCount() const { return alphaChannelCount; };
+       
+       inline const unsigned int & getArea() const { return area; };
+       inline const unsigned int & getVolume() const { return volume; };
+
+    // to-be-protected? :
+
+
+    template <class T>
+    Geometry & operator=(const std::vector<T> &v);
+    
+    //inline
+    Geometry & operator=(const Geometry &g);
+    
+    //Geometry & operator=(const Geometry &g);
+    
+    
+    string & toString(string & s) const;
+    
+       
+    protected:
+       // alphaChannelCount?
+       // area and volume
+       void update();
+       
+    unsigned int width;
+    unsigned int height;
+    unsigned int channelCount;
+       
+       unsigned int imageChannelCount;
+       unsigned int alphaChannelCount;
+
+    unsigned int area;
+    unsigned int volume;
+    
+       
+};
+
+
+
+/// Not ready.
+template <class T>
+Geometry & Geometry::operator=(const std::vector<T> &v){
+       int d = v.size();
+       for (int i = 0; i < d; ++i) {
+               if (i<d)
+                       (*this)[i] = v[i];
+               else
+                       (*this)[i] = 0;
+       }
+       update();
+       return (*this);
+}
+    
+ostream & operator<<(ostream &ostr,const Geometry &geometry);
+    
+
+}
+
+}
+
+#endif /*GEOMETRY_H_*/
diff --git a/drain/image/GradientOp.h b/drain/image/GradientOp.h
new file mode 100644 (file)
index 0000000..0b7f6aa
--- /dev/null
@@ -0,0 +1,152 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef GradientOP_H_
+#define GradientOP_H_
+
+
+#include <math.h>
+
+#include "ImageOp.h"
+
+namespace drain
+{
+
+namespace image
+{
+
+
+/**! Computes spatial horizontal derivative, dx.
+ *    src(i,j)-src(i-1,j). 
+ * 
+ */
+template <class T=unsigned char,class T2=unsigned char>
+class GradientOp : public ImageOp<T,T2> 
+{
+public:
+
+       //GradientOp() : bias(0), scale(1), offsetLo(1), offsetHi(1) {};
+
+       
+       void filter(const Image<T> &src,Image<T2> &dst,int diLow,int djLow,int diHigh,int djHigh) const {
+
+                       makeCompatible(src,dst);
+
+                       const int width  = src.getWidth();
+                       const int height = src.getHeight();
+
+                       const float bias  = this->getParameter("bias",0.0);
+
+                       const int iSpan = diLow+diHigh;
+                       const int jSpan = djLow+djHigh;
+                       const float span = sqrt(iSpan*iSpan + jSpan*jSpan);
+
+                       if (span == 0.0)
+                               throw (runtime_error("GradientOp: zero span"));
+
+                       const float scale = this->getParameter("scale",1.0f) / span;
+
+                       if (drain::Debug > 0){
+                               cerr << this->name;
+                               cerr << " bias=" << bias;
+                               cerr << " scale=" << scale;
+                               cerr << '\n';
+                       }
+
+                       Point2D<> pLo;
+                       Point2D<> pHi;
+                       const CoordinateHandler & h = src.getCoordinateHandler();
+
+                       for (int j = 0; j < height; j++) {
+                               for (int i = 0; i < width; i++) {
+                                       pLo.setLocation(i-diLow,j-djLow);
+                                       h.handle(pLo);
+                                       pHi.setLocation(i+diHigh,j+djHigh);
+                                       h.handle(pHi);
+                                       dst.at(i,j) = static_cast<T2>(bias + scale*(src.at(pHi) - src.at(pLo)));
+                               }
+                       }
+               };
+
+
+/*
+               
+  protected:
+       mutable float bias;  // TODO full scale! Intensity::
+       mutable float scale;
+       mutable int offsetLo;
+       mutable int offsetHi;
+       */
+       //virtual void filter(const Image<T> &src,Image<T2> &dst) const = 0;
+};
+
+
+/**! Computes spatial horizontal derivative, dx.
+ *    src(i,j)-src(i-1,j). 
+ * 
+ */
+template <class T=unsigned char,class T2=unsigned char>
+class GradientHorizontalOp : public GradientOp<T,T2>
+{
+public:
+       
+       GradientHorizontalOp(){
+               this->setInfo("Horizontal intensity Gradient","bias,scale,span","0,1,2");
+       };
+
+       void filter(const Image<T> &src,Image<T2> &dst) const {
+               const int span = this->getParameter("span",2);
+               if (span>0)
+                       GradientOp<T,T2>::filter(src, dst, span/2, 0, (span+1)/2, 0);
+               else
+                       GradientOp<T,T2>::filter(src, dst, span/2, 0, (span-1)/2, 0);
+       }
+               
+};
+
+/**! Computes spatial vertical derivative, dy.
+ *    src(i,j)-src(i,j-1). 
+ * 
+ */
+template <class T=unsigned char,class T2=unsigned char>
+class GradientVerticalOp : public GradientOp<T,T2>
+{
+public:
+       
+       GradientVerticalOp(){
+               this->setInfo("Vertical intensity Gradient","bias,scale,span","0,1,2");
+       };
+
+       void filter(const Image<T> &src,Image<T2> &dst) const {
+               const int span = this->getParameter("span",2);
+               if (span>0)
+                       GradientOp<T,T2>::filter(src, dst, 0,span/2, 0, (span+1)/2);
+               else
+                       GradientOp<T,T2>::filter(src, dst, 0,span/2, 0, (span-1)/2);
+       };
+       
+};
+
+}
+}
+
+
+#endif /*GradientOP_H_*/
diff --git a/drain/image/HighBoostOp.h b/drain/image/HighBoostOp.h
new file mode 100644 (file)
index 0000000..6beae46
--- /dev/null
@@ -0,0 +1,87 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef HIGHBOOSTOP_H_
+#define HIGHBOOSTOP_H_
+
+//#include "SlidingWindowAverage.h"
+#include "ImageOp.h"
+
+#include "FastAverageOp.h"
+#include "MathOpPack.h"
+
+
+namespace drain
+{
+
+namespace image
+{
+
+
+// Returns src - scale*src2
+template <class T=unsigned char,class T2=unsigned char>
+class ScaledSubtractionOp: public MathOp<T,T2>
+{
+public:
+
+   ScaledSubtractionOp(const string & p = "1.0,0") : MathOp<T,T2>(p) {
+          this->setInfo("An arithmetic operator.","scale,bias",p);
+   };
+
+   void filterValue(const T &src, const T2 &src2, T2 &dst) const {
+               dst = MathOp<T,T2>::limit(static_cast<double>(src) - this->scale*src2 + this->bias);
+   };
+};
+
+/// Enhances details by mixing original image and result of HighPassOp op.
+/**
+ */
+template <class T=unsigned char,class T2=unsigned char>
+class HighBoostOp : public ImageOp<T,T2>
+{
+public:
+
+       HighBoostOp(const string & p = "3,3,0.5"){
+       this->setInfo("Mixture of original and high-pass filtered image.",
+               "width,height,mixCoeff",p);
+       };
+    
+       virtual void filter(const Image<T> &src,Image<T2> &dst) const {
+               makeCompatible(src,dst);
+               
+               const int  width = this->parameters.get("width",3);
+               const int height = this->parameters.get("height",width);
+               FastAverageOp<T,T2>(width,height).filter(src,dst);
+               
+               const string &coeff = this->parameters.get("mixCoeff","0.5"); // ????
+               //SubtractionOp<T,T2>(1,0).filter(src,dst,dst);
+               ScaledSubtractionOp<T,T2>(coeff).filter(src,dst,dst);
+       };
+
+
+};
+
+
+}
+
+}
+
+#endif /*HIGHPASSOP_H_*/
diff --git a/drain/image/HighPassOp.h b/drain/image/HighPassOp.h
new file mode 100644 (file)
index 0000000..4539a98
--- /dev/null
@@ -0,0 +1,82 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef HIGHPASSOP_H_
+#define HIGHPASSOP_H_
+
+//#include "SlidingWindowAverage.h"
+#include "ImageOp.h"
+
+#include "FastAverageOp.h"
+#include "MathOpPack.h"
+
+
+namespace drain
+{
+
+namespace image
+{
+
+/// Difference of image and its low-pass filtered copy. 
+/**
+ *   Finds details of one-pixel scale. 
+ *   Implements
+ *   The simplest square high-pass filter
+ *   -1  -1  -1
+ *   -1   8  -1
+ *   -1  -1  -1 
+ * 
+ *   The result can be scaled with coeff and transposed with bias.
+ *   Internally, applies the fast SlidingStripeAverage and SubtractOp .
+ */
+template <class T=unsigned char,class T2=unsigned char>
+class HighPassOp : public ImageOp<T,T2>
+{
+public:
+
+       HighPassOp(const string & p = "3,3,1.0,0"){
+       this->setInfo("High-pass filter for recognizing details.",
+               "width,height,scale,coeff",p);
+       };
+    
+       virtual void filter(const Image<T> &src,Image<T2> &dst) const {
+               makeCompatible(src,dst);
+               
+               const int  width = this->parameters.get("width",3);
+               const int height = this->parameters.get("height",width);
+               FastAverageOp<T,T2>(width,height).filter(src,dst);
+               
+               const float scale = this->parameters.get("scale",-1.0/(float)width*(float)height); // ????
+               const float bias  = this->parameters.get("bias",0);
+               //ScaledSubtractionOp<T,T2>(scale).filter(src,dst,dst);
+               SubtractionOp<T,T2>(scale,bias).filter(src,dst,dst);
+               
+       };
+
+
+};
+
+
+}
+
+}
+
+#endif /*HIGHPASSOP_H_*/
diff --git a/drain/image/Image.h b/drain/image/Image.h
new file mode 100644 (file)
index 0000000..61460fd
--- /dev/null
@@ -0,0 +1,650 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef IMAGE_H_
+#define IMAGE_H_ "Image 1.3,  2010.05.05 Markus.Peura@fmi.fi"
+
+
+// TODO: rename to WritableImage ?
+
+#include <stdlib.h>  // for exit(-1) TODO: obsolete, use exceptions.
+#include <stdexcept>
+#include <vector>
+#include <ostream>
+
+#include "../util/Debug.h"
+#include "../util/Options.h"
+#include "Geometry.h"
+
+#include "Intensity.h"
+#include "CoordinateHandler.h"
+#include "ImageView.h"
+
+#include "Point.h"
+
+
+namespace drain
+{
+
+/// General debugging // TODO: move to more general place
+//static unsigned int Debug;
+
+/// Namespace for images and image processing tools.
+namespace image
+{
+
+template <class T>
+class ImageView;
+
+/// A class for images having using an own memory resource.
+/*!
+ *  Essentially, BufferedImage is a std::vector, hence a 1D structure in which
+ *  further dimensions are implemented by Geometry.
+ *  
+ *  @see BufferedView
+ * 
+ */
+template <class T=unsigned char>
+class Image  //: public Image
+{
+public:
+
+
+    Image();
+
+       // Creates a new image having the geometry of \code image. 
+       Image(const Image<T> &image);
+
+    Image(int width,int height,int channelCount=1,int alphaChannelCount=0);
+
+       virtual ~Image(){buffer = NULL;};
+       inline
+    virtual const Geometry & getGeometry() const {
+               return geometry;
+       }
+
+       /// Collapses the image to zero size, releasing memory. @see clear().
+       void reset(){
+               setGeometry(0,0,0,0);
+       }
+
+       /// Defines the dimensions of the image: width x height x (imageChannelCoount + alphaChannelCount).
+       /// Convenience. Calls setGeometry(geom).
+       /// Final!
+       void setGeometry(unsigned int width,unsigned int height,unsigned int imageChannelCount = 1,
+               unsigned int alphaChannelCount = 0){
+               setGeometry( Geometry(width,height,imageChannelCount,alphaChannelCount) ); 
+       };
+       
+       //inline   
+       virtual void setGeometry(const Geometry &g)
+       {
+               geometry.setGeometry(g);        
+               initialize();
+       }
+
+               
+       /// Sets the number of image channels and alpha channels.
+       inline 
+       void setChannelCount(unsigned int k,unsigned int ak = 0){
+               //const Geometry &g = getGeometry();
+               setGeometry(getWidth(),getHeight(),k,ak);       
+       };
+       
+       inline 
+       void setAlphaChannelCount(unsigned int k){
+               setGeometry(getWidth(),getHeight(),getImageChannelCount(), k);  
+       };
+       
+       // Frequently needed convenience functions.
+       inline const unsigned int & getWidth() const { return geometry.getWidth();};
+       inline const unsigned int & getHeight() const { return geometry.getHeight();};
+       inline const unsigned int & getChannelCount() const { return geometry.getChannelCount();};
+       inline const unsigned int & getImageChannelCount() const { return geometry.getImageChannelCount();};
+       inline const unsigned int & getAlphaChannelCount() const { return geometry.getAlphaChannelCount();};
+       inline const unsigned int & getVolume() const { return geometry.getVolume();};
+       
+       inline
+       const vector< ImageView<T> > & getChannelVector() const {
+               updateChannelVector();
+               return channel;
+       };
+
+
+       virtual const Image<T> &getChannel(unsigned int index) const;
+    virtual       Image<T> &getChannel(unsigned int index);
+    virtual const Image<T> &getAlphaChannel(unsigned int index = 0) const;
+    virtual       Image<T> &getAlphaChannel(unsigned int index = 0);
+       
+       /// Returns a view to image channels of this image.
+       virtual const Image<T> &getImageChannels() const {
+               views.resize(2);
+               views[0].viewChannels(*this,geometry.getImageChannelCount());
+               return views[0];
+       }
+       
+       /// Returns a view to image channels of this image.
+    virtual       Image<T> &getImageChannels(){
+       views.resize(2);
+       views[0].viewChannels(*this,geometry.getImageChannelCount());
+       return views[0];        
+    }
+    
+    /// Returns a view to alpha channels of this image.
+    virtual const Image<T> &getAlphaChannels() const {
+       views.resize(2);
+               views[1].viewChannels(*this,geometry.getAlphaChannelCount(),geometry.getImageChannelCount());
+               return views[1]; //alphaChannels;
+       }
+       
+       /// Returns a view to alpha channels of this image.
+    virtual       Image<T> &getAlphaChannels(){
+       views.resize(2);
+       views[1].viewChannels(*this,geometry.getAlphaChannelCount(),geometry.getImageChannelCount());
+       return views[1]; //alphaChannels;       
+    }
+   
+    // Extensions (new methods):
+    inline
+    void setCoordinateHandler(CoordinateHandler &handler);
+
+    inline
+    CoordinateHandler &getCoordinateHandler() const;
+
+
+       inline
+    const T & at(int i) const {
+       return (*buffer)[ address(i) ];  // slow?
+    }
+    
+    inline 
+    T & at(int i) {
+       return (*buffer)[ address(i) ];  // slow?
+    }
+
+       inline
+       const T & at(const int &i, const int &j) const {
+       return (*buffer)[ address(i,j) ];
+       };
+
+       inline
+       T & at(const int &i, const int &j){
+       return (*buffer)[ address(i,j) ];
+       };
+
+       inline
+       const T & at(const int &i, const int &j, const int &k) const {
+       return (*buffer)[ address(i,j,k) ];
+       };
+
+       inline
+       T & at(const int &i, const int &j, const int &k) {
+       return (*buffer)[ address(i,j,k) ];
+       };
+
+       inline
+       const T & at(const Point2D<int> &p) const {
+       return (*buffer)[ address(p.x,p.y) ];
+       };
+
+       inline
+       T & at(const Point2D<int> &p) {
+       return (*buffer)[ address(p.x,p.y) ];
+       };
+       
+       void getPixel(const Point2D<> &p, std::vector<T> & pixel) const;
+       void putPixel(const Point2D<> &p, const vector<T> & pixel);
+
+    
+    
+       /// For sequential (in-place) computations
+    virtual typename vector<T>::const_iterator begin() const {
+               typename vector<T>::const_iterator it = buffer->begin();
+               advance(it, offset);
+       return it;
+       }
+       
+    virtual typename vector<T>::iterator begin() {
+       typename vector<T>::iterator it = buffer->begin();
+               advance(it, offset);
+       return it;
+       }
+       
+    virtual typename vector<T>::const_iterator end() const {
+       typename std::vector<T>::const_iterator it = buffer->begin();
+               advance(it, offset + geometry.getVolume());
+       return it;
+       }
+       
+    virtual typename vector<T>::iterator end() {
+       typename std::vector<T>::iterator it = buffer->begin();
+               advance(it, offset + geometry.getVolume());
+       return it;
+       }
+
+       /// The 1D address of a pixel in the vector 
+    virtual inline int address(const int &i) const {
+        return i;  
+    };
+
+       /// The 1D address of a pixel in the vector 
+    virtual inline int address(const int &i,const int &j) const {
+        return (i + geometry.getWidth()*j);  
+    };
+
+    /// The 1D address of a pixel in the vector 
+    virtual inline int address(const int &i, const int &j, const int &k) const {
+        return (i + geometry.getWidth()*j + geometry.getArea()*k);
+    };
+
+
+       const std::vector<T> & getBuffer() const;
+
+       // This is dangerous. Consider resize().
+       std::vector<T> & getBuffer();
+  
+       
+    ///  Fills all the channels with %value.
+       inline
+    void fill(T value){
+       for (typename vector<T>::iterator it=begin(); it!=end(); it++)
+                 *it = value;
+    }
+
+    /// Sets intensities to zero. Does not change geometry. @see reset();
+    inline
+    void clear(){fill(0);};
+
+       // Tells if this image has an own buffer or views the buffer of another image.
+       bool hasOwnBuffer() const;
+
+       // Tells if image have share memory image.
+
+       bool hasOverlap(const Image<T> & img) const {
+               typename vector<T>::const_iterator b = begin();
+               typename vector<T>::const_iterator e = end();
+               typename vector<T>::const_iterator bi = img.begin();
+               typename vector<T>::const_iterator ei = img.end();
+               return ((bi>=b) && (bi<e)) || ((ei>=b) && (ei<e));
+       }
+
+       // Test with T != T2 will always fail.
+       template <class T2>
+       bool hasOverlap(const Image<T2> & img) const {
+               return false;
+       }
+
+       // Tells if this image has a alpha channel. Single-channel image never has.
+       bool hasAlphaChannel() const { return (getAlphaChannelCount() > 0);};
+
+       inline
+       const unsigned int & getOffset() const {return offset;};
+
+       /// Prints image debugging data to the standard error buffer.
+       void debug(ostream &ostr = cout) const;
+       
+       void setName(const string &n){ name = n; };
+
+       const string & getName() const { return name; };
+
+       /// Probably kicked to properties, see below.
+       string name;
+
+    Options properties;
+    
+    /*! Expected maximum intensity. This value will be used when scaling intensities. User may change the value. 
+     *  By default, for \em char and \em short int this value will be 255 and 65535, respectively, else 1.0. 
+     */     
+    inline
+       T getMax() const {return max;};
+       
+       inline
+       void setMax(T max = Intensity::max<T>()){this->max = max;};
+       
+protected:
+
+       ///
+       virtual void updateChannelVector() const;
+
+       // needed?
+       T max;
+
+       Geometry geometry;
+
+       // Initializes the channel vector, if not yet initialized.
+       virtual void initialize();
+
+    CoordinateHandler * coordinateHandler;
+    CoordinateHandler defaultCoordinateHandler;
+
+    std::vector<T> *buffer;
+    std::vector<T> ownBuffer;
+
+       unsigned int offset;
+       
+private:
+
+       /// Container for storing referneces to each channel. Lazy update. 
+       /// Changed only when a channel is requested with getChannel().
+    mutable vector< ImageView<T> > channel;
+    
+    mutable vector< ImageView<T> > views;
+    
+    // experimental
+    //mutable BufferedView<T> imageChannels;
+    //mutable BufferedView<T> alphaChannels;
+};
+
+
+
+//unsigned int ImageDebug = 0;
+
+
+
+
+
+template <class T>
+Image<T>::Image() //:  imageArray(&defaultVector)
+{
+       setMax();
+       buffer = & ownBuffer;
+    setGeometry(0,0,0);
+    setCoordinateHandler(this->defaultCoordinateHandler);
+}
+
+template <class T>
+Image<T>::Image(const Image<T> &image) //:  imageArray(&defaultVector)
+{
+       max = image.max;
+       buffer = & ownBuffer;
+       setGeometry(image.getGeometry());
+       setCoordinateHandler(this->defaultCoordinateHandler);
+}
+
+
+template <class T>
+Image<T>::Image(int width,int height, int imageChannelCount, int alphaChannelCount) 
+       // : imageArray(&defaultVector)
+{
+       setMax();
+       buffer =& ownBuffer;
+    setGeometry(width,height,imageChannelCount,alphaChannelCount);
+       setCoordinateHandler(this->defaultCoordinateHandler);
+}
+
+
+template <class T>
+void Image<T>::setCoordinateHandler(CoordinateHandler &handler)
+{
+    coordinateHandler = &handler;
+    coordinateHandler->setBounds(getWidth(), getHeight());
+}
+
+
+template <class T>
+CoordinateHandler &Image<T>::getCoordinateHandler() const
+{
+    return *coordinateHandler;
+}
+
+
+
+
+template <class T>
+void Image<T>::getPixel(const Point2D<int> &p,std::vector<T> &pixel) const
+{
+       static unsigned int i;
+    for (i=0; i<pixel.size(); i++)
+       pixel[i] = at(p.x,p.y,i); 
+
+};
+
+template <class T>
+void Image<T>::putPixel(const Point2D<int> &p,const std::vector<T> &pixel)
+{
+       static unsigned int i;
+    for (i=0; i<pixel.size(); i++)
+       at(p.x,p.y,i) = pixel[i]; 
+};
+
+template <class T>
+void Image<T>::updateChannelVector() const {
+
+       const unsigned int channelCount = getChannelCount();
+
+       if (channel.size() != channelCount){
+               channel.resize(channelCount);
+               for (unsigned int i = 0; i < channelCount; ++i) {
+                       channel[i].viewChannel(*this,i);
+               }
+       }
+
+}
+
+template <class T>
+Image<T> &Image<T>::getChannel(unsigned int index)
+{
+       
+       const unsigned int channelCount = getChannelCount();
+       
+       if (index >= channelCount){
+               // TODO Exception
+               cerr << " getChannel() illegal index: " << index << endl;
+               exit(-1);
+       }
+       
+       if (channelCount == 1){
+               return *this;
+       }
+       
+       // ERROR? Risk of missing channels when re-allocating?
+       /*
+       if (channel.size() != channelCount){
+               channel.resize(channelCount);
+       }
+       channel[index].viewChannel(*this,index);
+    */
+       updateChannelVector();
+
+    return channel[index];
+};
+
+template <class T>
+const Image<T> &Image<T>::getChannel(unsigned int index) const
+{
+       const unsigned int channelCount = getGeometry().getChannelCount();
+       
+       if (index >= channelCount){
+               // TODO Exception
+               //cerr << " getChannel() illegal index: " << index << endl;
+               exit(-1);
+       }
+       
+       if (channelCount == 1){
+               return *this;
+       }
+
+       /* vanha
+       if (channel.size() != channelCount){
+               channel.resize(channelCount);
+       }
+    channel[index].viewChannel(*this,index);
+    */
+       
+    // uusi
+    updateChannelVector();
+    
+    return channel[index];
+    
+};
+
+template <class T>
+Image<T> &Image<T>::getAlphaChannel(unsigned int index)
+{
+       unsigned int alphaChannelCount = getGeometry().getAlphaChannelCount();
+
+       if (index < alphaChannelCount){
+               //cerr << "getAlphaChannel returns " << getGeometry().getChannelCount()-1 - index << endl;
+               return getChannel(getGeometry().getChannelCount()-1 - index);
+       }
+       else {
+               cerr << " getAlphaChannel() illegal alpha index: " << index << endl;
+               exit(-1);
+       }
+               
+};
+
+template <class T>
+const Image<T> &Image<T>::getAlphaChannel(unsigned int index) const
+{
+       unsigned int alphaChannelCount = getGeometry().getAlphaChannelCount();
+
+       if (index < alphaChannelCount){
+               return getChannel(getGeometry().getChannelCount()-1 - index);
+       }
+       else {
+               cerr << " getAlphaChannel() illegal index: " << index << endl;
+               exit(-1);
+       }
+};
+
+/**
+ *  Notice, this returns the whole vector, even when called
+ *
+ */
+template <class T>
+const vector<T> & Image<T>::getBuffer() const {
+       return *(this->buffer); 
+}
+
+template <class T>
+vector<T> & Image<T>::getBuffer() {
+       return *(this->buffer); 
+}
+
+
+template <class T>
+void Image<T>::initialize(){
+       
+       offset = 0;
+       channel.clear();
+       
+       //cerr << "BufferedImage<T>::initialize BEGIN " << this->geometry.toString(gStr) << endl;
+       //cerr << " name='" << this->name << "'" <<endl;
+       /*
+       if (!hasOwnBuffer()){
+               cerr << "WARNING: tried initializing (resizing) an extrenal buffer" << endl;
+               return; 
+       }
+       */
+
+       if ((*buffer).size() != getVolume()){
+
+               if (!hasOwnBuffer())
+                       throw runtime_error("WARNING: tried initializing (resizing) an extrenal buffer");
+
+               (*buffer).resize(getVolume());
+       }
+       
+       
+}
+
+template <class T>
+bool Image<T>::hasOwnBuffer() const
+{
+       return (buffer == &ownBuffer);
+}
+
+
+// DEBUGGING
+// TODO: change channels to getChVect
+template <class T>
+void Image<T>::debug(ostream &ostr) const
+{
+       int channelCount = getGeometry().getChannelCount();
+       vector<const Image<T> *> channels(channelCount+1);
+       
+       //cerr << "...debug: init THIS" << endl;
+       channels[0] = this;
+       for (int i = 0; i < channelCount; ++i) {
+               //cerr << "...debug: init channel #" << i << endl;
+               channels[i+1] = &getChannel(i);
+       }
+       
+       ostr << "\nImage debug info for "<< this->name << ", type=";
+       //cout << " Type: ";
+       if (hasOwnBuffer())
+               ostr << "BUFFER\n";
+       else
+           ostr<< "REFERENCE\n";
+       
+       if (this->coordinateHandler == NULL)
+               ostr << "WARNING: CoordHandler undefined!";
+         //ostream << " Coord Handler defined: " << (this->coordinateHandler != NULL) << "\n";
+       
+       for (int i = 0; i < channelCount+1; ++i) {
+               const Image<T> *img = channels[i];
+               string info;
+               long unsigned int a  = 0;
+               //long unsigned int a2 = 0;
+               long unsigned int m = 0;
+               if (img != NULL){
+                       img->getGeometry().toString(info);
+                       a  = img->address(0,0);
+                       //a2 = img->address(1,1);
+                       m = (long unsigned int)&(img->at(0,0));
+               }
+               else 
+                       info = "Null";
+                       
+               if (i == 0)
+                       ostr << "image:";
+               else    
+                       ostr << "ch #" << (i-1);
+                       
+               //ostr << ":\t" << info << "(" << a << ','<< a2 << ")" << "# " << m << endl;
+               ostr << ":\t" << info << " [" << a << "]" << " #" << m << endl;
+       }
+       ostr << " Default array: size= "<< ownBuffer.size( )
+                << " address= " << (long unsigned int)&ownBuffer[0] << endl;
+                
+               
+               
+       ostr << " Current array: size= " << buffer->size( )
+                << " address= " << (long unsigned int)&((*buffer)[0]) << endl
+                << "  = " << (long unsigned int)&(*begin()) << " -- " << (long unsigned int)&(*end()) << endl;
+       
+}
+
+/*
+template <class T>
+void BufferedImage<T>::setName(const string &n){
+       name = n;       
+}
+*/
+
+
+}
+
+}
+
+#endif /*BUFFEREDIMAGE_H_*/
diff --git a/drain/image/ImageOp.h b/drain/image/ImageOp.h
new file mode 100644 (file)
index 0000000..7839009
--- /dev/null
@@ -0,0 +1,250 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef IMAGE_OP_H_
+#define IMAGE_OP_H_
+
+//#include <exception>
+#include <stdexcept>
+#include <list>
+
+#include "Image.h"
+
+namespace drain
+{
+
+namespace image
+{
+
+/**  
+ *  Contract.
+ * 
+ *  If (dst == src), and computation involves writing to dst in a way that affects
+ *  succeeding computation, operator should apply an internal temporary image.   
+ * 
+ *  Experimental, optional: Using (Data) parameters only.
+ *  The allowed parameters are defined through the list parameterNames, which also
+ *  states the order of parameters, if given through setParameters(const Data & d). 
+ * 
+ */
+template <class T,class T2>
+class ImageOp  
+{
+public:
+
+       // \bug TODO: problem default value must be given in this call.
+   ImageOp(const string &name = "ImageOp",const string &description="Image processing operator.",
+                  const string & parameterNames="", const string & defaultValues = "") : name(name) {
+          setInfo(description,parameterNames,defaultValues);
+   };
+   
+   virtual ~ImageOp(){};
+
+   virtual void makeCompatible(const Image<T> &src,Image<T2> &dst) const  {
+       //cerr << "ImageOp\n";
+               dst.setGeometry(src.getGeometry());
+       //coordinateHandler->setBounds(src.getWidth(),src.getHeight());  // ?? both
+       // TODO: change handler to similar  
+       dst.getCoordinateHandler().setBounds(src.getWidth(),src.getHeight());  // ??
+       };
+   
+       /*
+       void setCoordinateHandler(CoordinateHandler &handler){
+       coordinateHandler = &handler;
+       };
+       */
+   
+
+   /// TODO: const
+   virtual void filter(const Image<T> &src,Image<T2> &dst) const = 0;
+   /*
+   {
+         throw runtime_error("(ImageOp unimplemented.)");
+   };
+    */
+
+
+   /// Container for storing operator parameters as KEY,VALUE pairs.  
+   /*! I am still hesitating with this---
+    *   Advantages:
+    *   The parameters can be more easily...
+    *   # read from a command file
+    *   # dumped, for debugging & reflection 
+    *   Disadvantages:
+    *   # using this hides the actual parameters from the interface
+    *   # type checks?  
+    */
+    //void setParameter(const string  
+    //Data info;
+    
+    //const string & getHelp(){ return parameters; };
+    
+    /// Throws exception if the number of parameters exceed the number.d 
+    void setParameters(const string & p){
+       parameters.set(p);
+    }
+    
+    bool hasParameter(const string & key) const {
+        return parameters.hasKey(key);
+    }
+       
+    void setParameter(const string & key,const Data & value){
+               parameters[key] = value;
+    }
+       
+
+    //template <>
+    const char * getParameter(const string & key,const char *defaultValue) const {
+               return parameters.get(key,defaultValue).c_str();
+       }
+
+
+       template <class T3>
+       T3 getParameter(const string & key,T3 defaultValue) const {
+               return parameters.get(key,defaultValue);
+    }
+
+       inline
+       const string & getParameterNames() const {
+               return parameters.getParameterNames();
+    }
+
+       inline
+       const string &getParameterValues() const {
+               return parameters.getParameterValues();
+       }
+       /*
+       inline
+       const string getParameterValues() const {
+               stringstream sstr;
+               string separator = "";
+               for (map<string,Data>::const_iterator it = opt.begin(); it != opt.end(); it++){
+                       sstr << separator << it->second;
+
+               }
+               return sstr.str();
+               //return parameters.getParameterNames();
+       }
+       */
+
+       inline
+       const Options & getParameters() const {
+               return parameters;
+       }
+
+       const string & getName() const { return this->name; };
+
+       // protect
+       Options parameters;
+       
+       /// Sets the applicable parameters, their order for setParameters() calls as well as usage information to be retrieved by help calls.
+       void setInfo(const string & description, const string & parameterNames, const string & defaultValues = ""){
+               parameters.clear();
+       parameters.setDescription(description);
+       parameters.setParameterNames(parameterNames);
+       parameters.set(defaultValues);
+    };
+   
+    void help(ostream &ostr) const {
+       ostr << name << ": ";
+       ostr << parameters.getDescription() << '\n';
+       ostr << "Parameters: " << parameters.getParameterNames() << '\n';
+       ostr << "Current values:\n" << parameters << '\n';      
+    } 
+    
+    string help() const {
+       stringstream s;
+       help(s);
+       return s.str(); 
+    }
+       // comma-separated list of parameters in the order specified by parameters.getAllowedParameters()
+       // void setParameters(const string & s);
+       
+
+
+protected:
+
+    /// Utility for re-calling filter() with a temporary dst image, when dst==src would corrupt the computation.
+    /**  Experimental. Place this in the first line of your filter() as simply as:
+     *  \code
+     *  if (filterWithTmp(src,dst))
+     *         return;
+     *  \endcode
+     */
+    bool filterWithTmp(const Image<T> &src,Image<T2> &dst) const {
+
+       if (src.hasOverlap(dst)){
+               Image<T2> tmp;
+               filter(src,tmp);
+               ///
+               dst.setGeometry(tmp.getGeometry());
+               /// Copy
+               typename std::vector<T2>::const_iterator si = tmp.begin();
+            typename std::vector<T2>::iterator di = dst.begin();
+            const typename std::vector<T2>::const_iterator dEnd = dst.end();
+               while (di!=dEnd){
+                       *di = *si;
+                       si++;
+                       di++;
+            }
+               return true;
+       }
+       else
+               return false;
+    }
+
+    /// Should the operation be carried out to each channel
+    bool filterEachChannel(const Image<T> &src,Image<T2> &dst) const {
+
+       const unsigned int c = src.getImageChannelCount();
+       dst.setGeometry(src.getWidth(), src.getHeight(), c, dst.getAlphaChannelCount());
+
+       if (c == 1)
+               return true;
+
+       for (unsigned int i = 0; i < c; ++i)
+               filter(src.getChannel(i),dst.getChannel(i));
+
+       return false;
+    }
+
+    // todo: make it const
+    string name;
+   
+    //drain::MapWrapper<string,Data> parameterWrapper;
+    /*
+       CoordinateHandler *coordinateHandler;
+       CoordinateHandler defaultCoordinateHandler;
+       Mirror mirror;
+       Wrapper wrapper;
+       */
+};
+
+//static unsigned int ImageOpDebug;
+//template <class T,class T2>
+//unsigned int ImageOp<T,T2>::debug = 0;
+   
+   
+   
+}
+}
+
+#endif /* IMAGE_OP_H_ */
diff --git a/drain/image/ImageView.h b/drain/image/ImageView.h
new file mode 100644 (file)
index 0000000..15dfcaf
--- /dev/null
@@ -0,0 +1,235 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef IMAGEVIEW_H_
+#define IMAGEVIEW_H_
+
+#include <sstream> 
+#include <stdexcept>
+#include "Image.h"
+
+namespace drain
+{
+namespace image
+{
+
+
+template <class T>
+class Image;
+
+/// A reference to one or several subsequent channels of a BufferedImage.
+/** It does not have memory itself. 
+ *  @see BufferedImage
+ *  Note: the name is not logical, as this view does NOT have a buffer...l
+ */
+template <class T=unsigned char>
+class ImageView : public Image<T>
+{
+public:
+       
+       ImageView(){ this->buffer=NULL; };
+       ImageView(const Image<T> &image, int channelStart=0, int channelCount=1);
+
+       /// Views all the channels of an image.
+       void viewImage(const Image<T> & src){
+               viewChannels(src,src.getChannelCount());
+               setGeometry(src.getGeometry());
+       }
+       
+       /// Views a single channel of an image.
+       void viewChannel(const Image<T> &image,unsigned int channel){
+               viewChannels(image,1,channel);
+       }
+       
+       /// Views \code channelCount channels of an image, starting from \code channelCount.
+       void viewChannels(const Image<T> &image,unsigned int channelCount, unsigned int channelStart = 0);
+       
+       
+       /// Concatenates all the channels to a single view.
+       void viewImageFlat(const Image<T> &src){
+               viewChannels(src,src.getChannelCount());
+               this->geometry.setGeometry(src.getWidth(),src.getHeight()*src.getChannelCount(),1);
+               cerr << "viewImageFlat" << this->geometry << endl;
+       };
+       
+       
+       //void setView(const Image<T> &image, unsigned int channelCount = 0, unsigned int channelStart = 0);
+
+       void setView(std::vector<T> & v,unsigned int width,unsigned int height, 
+               unsigned int imageChannelCount, unsigned int alphaChannelCount, unsigned int channelStart=0); 
+       
+       /// Overrides (Image). Does nothing.
+    //virtual void setGeometry(int width,int height,int imageChannelCount=1,int alphaChannelCount=0);
+
+       /// Overrides that of BufferedImage. Does nothing but checks if request would be illegal ie. would change geometry.
+    void setGeometry(const Geometry &geometry){
+               // todo: warning if tries to change?
+               if (this->getGeometry().getVolume() != geometry.getVolume()){
+                       stringstream sstr;
+                       sstr << "ImageView<T>::setGeometry: tried to change ";
+                       sstr << " '" << this->getGeometry() << "' to '" << geometry << "'\n";
+                       //throw (sstr.str());
+                       throw runtime_error(sstr.str());
+               }
+               // 
+       };
+       
+    
+    
+       /*
+       virtual const BufferedImage<T> &getChannel(int index) const;
+    virtual       BufferedImage<T> &getChannel(int index) ;
+       */
+
+       virtual inline int address(const int &i) const {
+       return (this->offset + i);  };
+       
+       virtual inline int address(const int &i,const int &j) const {
+       return (this->offset + i + this->geometry.getWidth()*j);  };
+
+    virtual inline int address(const int &i,const int &j,const int &k) const {
+               return (this->offset + i + this->geometry.getWidth()*j + this->geometry.getArea()*k);  };
+
+       
+       bool isSet(){ return !this->hasOwnBuffer(); }; 
+
+protected:
+
+       /// For ImageView.
+     //        int channelStart;
+       
+       
+       /// For ImageView.
+    //int offset;
+    
+       
+    /// Does nothing.
+       virtual void initialize(){
+               //this->channel.clear();  // is this needed?
+               cerr <<  "Bview init" << endl;
+       };
+       
+};
+
+template <class T>
+ImageView<T>::ImageView(const Image<T> &image,int channelStart, int channelCount){
+       //cerr << "ImageView() --> "  << image.name << endl;
+       viewChannels(image,channelCount,channelStart);
+}
+
+
+       
+/// Views \code channelCount channels of an image, starting from \code channelCount.
+template <class T>
+void ImageView<T>::viewChannels(const Image<T> &image,unsigned int channelCount, unsigned int channelStart)
+{
+       //void  ImageView<T>::setView(const BufferedImage<T> &image, unsigned int channelCount, unsigned int channelStart){
+       
+       const Geometry g = image.getGeometry();
+       
+       if (channelCount == 0)  // ???????????
+               channelCount = g.getChannelCount(); 
+       
+       const int i = g.getImageChannelCount();
+       const int a = g.getAlphaChannelCount();
+
+       int imageChannelCount = i - channelStart;
+       if (imageChannelCount < 0)
+               imageChannelCount = 0;
+       
+       int alphaChannelCount = channelCount - imageChannelCount;
+       if (alphaChannelCount > a)
+               alphaChannelCount = a;
+
+
+       /*
+       unsigned int imageChannelCount; 
+       unsigned int alphaChannelCount;
+       
+       if ((channelStart + channelCount) > i){
+               imageChannelCount = i - channelStart;
+       alphaChannelCount = channelCount - imageChannelCount;
+       }       
+    else {
+       imageChannelCount = channelCount;
+       alphaChannelCount = 0;
+    }
+    */
+       
+       
+       stringstream sstr;
+       sstr << image.name << "{" << channelStart;
+       if (channelCount > 1){
+               sstr <<  "..." << (channelStart+channelCount-1);
+       }
+       sstr <<  "}";
+       this->name = sstr.str();
+       
+       //cerr << "creating view " << this->name << endl;
+       
+       int currentChannelStart;
+       if (g.getArea() != 0)
+                currentChannelStart = image.getOffset()/g.getArea();
+       else
+               currentChannelStart = 0; // correct
+       
+       setView((vector<T> &)image.getBuffer(), g.getWidth(), g.getHeight(),
+               imageChannelCount,alphaChannelCount,channelStart + currentChannelStart);
+       
+       setCoordinateHandler(image.getCoordinateHandler()); 
+}
+
+
+template <class T>
+void ImageView<T>::setView(std::vector<T> & v,unsigned int width,unsigned int height,
+       unsigned int imageChannelCount, unsigned int alphaChannelCount,unsigned int channelStart) 
+{
+               Geometry &g = this->geometry;
+               
+               this->setMax();
+       
+               if (&(this->ownBuffer) != &v)
+                       this->ownBuffer.resize(0);
+               
+       this->buffer = &v;
+
+       g.setGeometry(width,height,imageChannelCount,alphaChannelCount);
+       setCoordinateHandler(this->defaultCoordinateHandler);
+
+               this->offset = g.getArea() * channelStart;
+       
+               if ((this->offset + g.getVolume()) > this->buffer->size()){
+                       stringstream sstr;
+                       sstr << "ImageView<T>::setView: requested allocation exceeds available buffer.\n";
+                       sstr << " buffer size :" << this->buffer->size() << endl;
+                       sstr << " offset :" << this->offset << endl;
+                       sstr << " volume :" << g.getVolume() << endl;
+                       throw runtime_error(sstr.str());
+               }
+       
+}
+
+
+}
+}
+
+#endif /*BUFFEREDVIEW_H_*/
diff --git a/drain/image/Intensity.h b/drain/image/Intensity.h
new file mode 100644 (file)
index 0000000..4dab2c3
--- /dev/null
@@ -0,0 +1,76 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef INTENSITY_H_
+#define INTENSITY_H_
+
+
+#include <limits>
+#include <cmath>
+
+namespace drain
+{
+
+namespace image
+{
+
+
+//template <class T>
+class Intensity
+{
+public:
+
+       /// For floating-point values, returns 0.0, otherwise the lower limit of the range.
+       template <class T>
+       inline
+       static T min(){ return std::numeric_limits<T>::is_integer ? std::numeric_limits<T>::min() : 0; };
+
+       /// For floating-point values, returns 1.0, otherwise the upper limit of the range.
+       template <class T>
+       inline
+       static T max(){ return std::numeric_limits<T>::is_integer ? std::numeric_limits<T>::max() : 1; };
+
+       /// Value to be used in inverting numbers; returns 0.0 for floating-point values.
+       template <class T>
+       inline
+       static T maxOrigin(){ return std::numeric_limits<T>::is_integer ? std::numeric_limits<T>::max() : 0; };
+
+
+       /// Ensures that the value stays within the numeric range of T2.
+       template <class T>
+       static T limit(double x){
+               static const double xMin = static_cast<double>(
+                               numeric_limits<T>::is_integer ?  std::numeric_limits<T>::min() : -std::numeric_limits<T>::max() );
+               static const double xMax = static_cast<double>( std::numeric_limits<T>::max() );
+               x = std::max(xMin,x);
+               x = std::min(x,xMax);
+               return static_cast<T>(x);
+       };
+
+};
+
+
+
+}
+
+}
+
+#endif /*INTENSITY_H_*/
diff --git a/drain/image/MagickDrain.h b/drain/image/MagickDrain.h
new file mode 100644 (file)
index 0000000..b4ac08a
--- /dev/null
@@ -0,0 +1,320 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+
+#ifndef MAGICKDRAIN_H_
+#define MAGICKDRAIN_H_
+
+
+
+
+
+#ifdef DRAIN_MAGICK_yes
+#include <Magick++.h>
+#endif
+
+/// In compilation, use "Magick*-config" to get libs and includes right.
+
+
+#include "Image.h"
+
+//#include "Point.h"
+
+namespace drain
+{
+
+namespace image
+{
+
+class MagickDrain
+{
+public:
+       //MagickDrain();
+       //virtual ~MagickDrain();
+       
+#ifdef DRAIN_MAGICK_yes
+       /// Converts drain::image to MagickImage
+       /// Does not support other types than <unsigned char>
+       template<class T>
+       static void convert(const Image<T> &drainImage, Magick::Image &magickImage);
+       
+       /// Converts MagickImage to drain::image
+       /// Does not support other types than <unsigned char>
+       template<class T>
+       static void convert(Magick::Image &magickImage, Image<T> &drainImage);
+       
+#endif
+       
+};
+
+
+
+#ifdef DRAIN_MAGICK_yes
+
+// Looks like Magick inverts alpha channel
+template <class T>
+void MagickDrain::convert(const Image<T> &drainImage, Magick::Image &magickImage) {
+
+       Magick::ImageType type = Magick::TrueColorMatteType;
+
+       int imageChannels = drainImage.getImageChannelCount();
+       int alphaChannels = drainImage.getAlphaChannelCount();
+
+       if (alphaChannels > 1){
+               std::cerr << "Warning: multiple alpha channel image, using 1st only \n";
+               std::cerr << "  Image geometry:" << drainImage.getGeometry() << '\n';
+               alphaChannels = 1;
+       }
+
+       const Image<T> *red = NULL, *green = NULL, *blue = NULL, *alpha = NULL;
+
+       //string info;
+       //drainImage.getGeometry().toString(info);
+       //cout << "MagickDrain: " << drainImage.getGeometry() << endl;
+       /*
+       cerr << "img.chs:" << imageChannels  << endl;
+       cerr << "alpha.chs:" << alphaChannels  << endl;
+       cerr << "Extra debug:"  << endl;
+       drainImage.debug();
+        */
+
+       switch (imageChannels){
+       case 0:
+               if (alphaChannels == 0){
+                       std::cerr << "Error: zero channel image!\n";
+               }
+               else {
+                       type = Magick::GrayscaleType;
+                       red   =& drainImage.getChannel(0);
+                       green =& drainImage.getChannel(0);
+                       blue  =& drainImage.getChannel(0);
+               }
+               break;
+       case 1:
+               type = Magick::GrayscaleType;
+               red   =& drainImage.getChannel(0);
+               green =& drainImage.getChannel(0);
+               blue  =& drainImage.getChannel(0);
+               if (alphaChannels > 0) {
+                       type = Magick::GrayscaleMatteType;
+                       alpha =& drainImage.getAlphaChannel();
+               }
+               break;
+       case 2:
+               if (!alphaChannels) {
+                       type  =  Magick::GrayscaleMatteType;
+                       std::cerr << "Notice: 2 channel image, storing 2nd as alpha channel \n";
+                       red   =& drainImage.getChannel(0);
+                       green =& drainImage.getChannel(0);
+                       blue  =& drainImage.getChannel(0);
+                       alpha =& drainImage.getChannel(1);
+               }
+               else {
+                       std::cerr << "Notice: (2+alpha ) channel image, creating 'RGAA' image instead of 'RGBA'.\n";
+                       type = Magick::TrueColorMatteType;
+                       red   =& drainImage.getChannel(0);
+                       green =& drainImage.getChannel(1);
+                       blue  =& drainImage.getAlphaChannel();
+                       alpha =& drainImage.getAlphaChannel();
+               }
+               break;
+       case 3:
+               type  = Magick::TrueColorType;
+               red   =& drainImage.getChannel(0);
+               green =& drainImage.getChannel(1);
+               blue  =& drainImage.getChannel(2);
+               if (alphaChannels){
+                       type  =  Magick::TrueColorMatteType;
+                       alpha =& drainImage.getAlphaChannel();
+               }
+               break;
+       default:
+               type  =  Magick::TrueColorMatteType;
+               red   =& drainImage.getChannel(0);
+               green =& drainImage.getChannel(1);
+               blue  =& drainImage.getChannel(2);
+
+               if (alphaChannels){
+                       cerr << "Warning: (" << imageChannels << "+alpha) channel image, using (3+alpha). \n";
+                       alpha =& drainImage.getAlphaChannel();
+               }
+               else if (imageChannels == 4){
+                       cerr << "Notice: 4 channel image, storing 4th as alpha channel \n";
+                       alpha =& drainImage.getChannel(3);
+               };
+               imageChannels = 3;      // WHY?
+       }
+
+
+       const int width = drainImage.getGeometry().getWidth();
+       const int height = drainImage.getGeometry().getHeight();
+
+
+       //Magick::Image magickImage(Magick::Geometry(width,height),Magick::Color("black"));
+
+
+       try {
+               magickImage.classType(Magick::DirectClass);  // here?
+               magickImage.size(Magick::Geometry(width,height));
+               //magickImage.read("cow.png");
+
+               magickImage.modifyImage(); // actually: _prepare_to_ modify
+               magickImage.type(type);   //  The order copied from Magick examples
+
+
+               Magick::PixelPacket *pixel_cache = magickImage.getPixels(0,0,width,height);
+
+
+               Point2D<> p;
+               int &i = p.x;
+               int &j = p.y;
+
+               // TODO: const unsigned drainBits = numeric_limits<unsigned char>::is_integer ?  numeric_limits<unsigned char>::max() :  //
+               const unsigned int drainMax = numeric_limits<unsigned char>::max();
+               const int shiftBits = magickImage.depth() - 8;  // ????
+
+
+               int rowAddress = 0;
+               int address = 0;
+
+               for (j=0; j<height; j++){
+                       rowAddress = j*width;
+
+                       for (i=0; i<width; i++){
+                               address = i + rowAddress;
+                               pixel_cache[address].red =   (red->at(i,j)   << shiftBits);
+                               pixel_cache[address].green = (green->at(i,j) << shiftBits);
+                               pixel_cache[address].blue =  (blue->at(i,j)  << shiftBits);
+                               if (alpha != NULL)
+                                       //pixel_cache[address].opacity = ((alpha->at(i,j))  << shiftBits);  //WARNING!!
+                                       // Looks like Magick does NOT invert alpha channel IN THIS DIRECTION
+                                       pixel_cache[address].opacity = ((drainMax-alpha->at(i,j))  << shiftBits);  //WARNING!!
+                       }
+               }
+               //cerr << "synching" << endl;
+               //pixel_cache[rowAddress + 10 + width] = Magick::Color("green");
+
+               magickImage.syncPixels();
+       }
+       catch (Magick::Error& e) {
+               // because 'Error' is derived from the standard C++ exception, it has a 'what()' method
+               cerr << "a Magick++ error occurred: " << e.what() << endl;
+       }
+       catch ( ... ) {
+               cerr << "MagickDrain: an unhandled error has occurred; exiting application." << endl;
+               exit(1);
+       }
+
+
+       // cerr << "magickImage.type = " << magickImage.type() << endl;
+
+       // Store comments as KEY=VALUE pairs
+       stringstream sstr;
+       map<string,Data>::const_iterator it;
+       for (it = drainImage.properties.begin(); it != drainImage.properties.end(); it++){
+               sstr << it->first << '=' << it->second << '\n';
+       }
+       magickImage.comment(sstr.str());
+
+
+}
+#endif
+
+
+
+#ifdef DRAIN_MAGICK_yes
+
+// Looks like Magick inverts alpha channel
+template <class T>
+void  MagickDrain::convert(Magick::Image &magickImage, Image<T> &drainImage) {
+
+       const int w = magickImage.columns();
+       const int h = magickImage.rows();
+
+       //drainImage.setGeometry(w,h)
+
+       // TODO: redChannel = &drainImage.at(0,0,0);
+
+       switch (magickImage.type()){
+       case Magick::GrayscaleType:
+               drainImage.setGeometry(w,h,1);
+               magickImage.write(0,0,w,h,"I",Magick::CharPixel,&drainImage.at(0,0));
+               break;
+       case Magick::GrayscaleMatteType:
+               drainImage.setGeometry(w,h,1,1);
+               magickImage.write(0,0,w,h,"I",Magick::CharPixel,&drainImage.at(0,0,0));
+               magickImage.write(0,0,w,h,"A",Magick::CharPixel,&drainImage.at(0,0,1));
+               break;
+               //    case Magick::RGBColorspace:
+       case Magick::PaletteType: // just test status
+       case Magick::TrueColorType:
+               drainImage.setGeometry(w,h,3);
+               magickImage.write(0,0,w,h,"R",Magick::CharPixel,&drainImage.at(0,0,0));
+               magickImage.write(0,0,w,h,"G",Magick::CharPixel,&drainImage.at(0,0,1));
+               magickImage.write(0,0,w,h,"B",Magick::CharPixel,&drainImage.at(0,0,2));
+               break;
+       case Magick::PaletteMatteType:
+       case Magick::TrueColorMatteType:
+               drainImage.setGeometry(w,h,3,1);
+               magickImage.write(0,0,w,h,"R",Magick::CharPixel,&drainImage.at(0,0,0));
+               magickImage.write(0,0,w,h,"G",Magick::CharPixel,&drainImage.at(0,0,1));
+               magickImage.write(0,0,w,h,"B",Magick::CharPixel,&drainImage.at(0,0,2));
+               magickImage.write(0,0,w,h,"A",Magick::CharPixel,&drainImage.at(0,0,3));
+               //      magickImage.write(0,0,w,h,"A",Magick::CharPixel,&drainImage.alphaChannel(0)[0]);
+               break;
+               //    default:
+       default:
+               stringstream sstr;
+               sstr << "operator<<(image,magickImage) : Magick type " << magickImage.type() << " not handled.";
+               throw runtime_error(sstr.str());
+       }
+
+       // TODO contradictory!
+       if (drainImage.getAlphaChannelCount()>0){
+               //Image<> & alpha = drainImage.getAlphaChannel();
+               //NegateOp<>().filter(alpha,alpha);  // Looks like Magick inverts alpha channel IN THIS DIRECTION.
+               //ScaleOp<>(-1.0,255).filter(alpha,alpha);  // Looks like Magick inverts alpha channel
+       }
+
+       std::stringstream sstr(magickImage.comment());  // dont touch sstr!!
+       drainImage.properties.reader.read(sstr);
+
+       if (drain::Debug > 5){ // TODO static (does not work)
+                       cerr << "read magickImage.type  = " << magickImage.type() << '\n';
+                       cerr << "comment='" << magickImage.comment() << "'\n";
+                       cerr << "::::::::::::::::::::\n";
+                       cerr << drainImage.properties;
+       }
+
+
+}
+
+#endif
+
+
+
+
+}  // image
+
+}  // drain
+
+#endif /*MAGICKDRAIN_H_*/
+//#endif // ImageMagick
diff --git a/drain/image/MarginalStatisticOp.h b/drain/image/MarginalStatisticOp.h
new file mode 100644 (file)
index 0000000..88d5972
--- /dev/null
@@ -0,0 +1,165 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef MARGINALSTATISTICOP_H_
+#define MARGINALSTATISTICOP_H_
+
+#include "math.h"
+
+#include "../util/Histogram.h"
+
+
+#include "ImageOp.h"
+
+
+namespace drain
+{
+
+namespace image
+{
+
+
+/** Computes horizontal or vertical intensity statistics: average, sum, ...
+ */
+template <class T=unsigned char,class T2=unsigned char> //,class D=int>
+class MarginalStatisticOp : public ImageOp<T,T2>
+{
+public:
+    
+    MarginalStatisticOp(const string & p = "horz,a"){
+       this->setInfo("Computes statistics on horizontal or vertical lines. ","mode,stat",p);
+       };
+
+
+    /*
+    virtual void makeCompatible(const Image<T> &src,Image<T2> &dst) const  {
+
+       const string &mode = this->parameters.get("mode","horz");
+       const string &stat = this->parameters.get("stat","avg");
+
+
+       if (mode == "horz"){
+                       dst.setGeometry(n,src.getHeight(),src.getImageChannelCount(),src.getAlphaChannelCount());
+               }
+               else if (mode == "vert"){
+                       dst.setGeometry(src.getWidth(),n,src.getImageChannelCount(),src.getAlphaChannelCount());
+               }
+               else {
+
+               }
+
+       /// dst.getCoordinateHandler(); TODO
+    };
+    */
+
+       void filter(const Image<T> &src, Image<T2> &dst) const {
+
+               Histogram<T,T2>  histogram(256);
+
+               const string &mode = this->parameters.get("mode","horz");
+               const string &stat = this->parameters.get("stat","a");
+       const unsigned int n = stat.size();
+
+               const unsigned int width  = src.getWidth();
+       const unsigned int height = src.getHeight();
+       const unsigned int iChannels = src.getImageChannelCount();
+       const unsigned int aChannels = src.getAlphaChannelCount();
+       const unsigned int channels = iChannels + aChannels;
+
+               if (mode == "horz"){
+                       dst.setGeometry(n,height,iChannels,aChannels);
+                       histogram.setSampleCount(height);
+                       for (unsigned int k=0; k<channels; k++){
+                               for (unsigned int j=0; j<height; j++){
+                                       histogram.clearBins();
+                                       //histogram.clear();
+                                       for (unsigned int i=0; i<width; i++)
+                                               histogram[src.at(i,j)]++;
+                                       for (unsigned int l=0; l<n; l++){
+                                               switch (stat[l]){
+                                               case 'M':
+                                                       dst.at(l,j,k) = histogram.getMax();
+                                                       break;
+                                               case 'm':
+                                                       dst.at(l,j,k) = histogram.getMin();
+                                                       break;
+                                               case 'd':
+                                                       dst.at(l,j,k) = histogram.getMedian();
+                                                       break;
+                                               case 'S':
+                                                       dst.at(l,j,k) = histogram.getStdDeviation();
+                                                       break;
+                                               case 's':
+                                                       dst.at(l,j,k) = histogram.getSum();
+                                                       break;
+                                               case 'a':
+                                                       dst.at(l,j,k) = histogram.getMean();
+                                                       break;
+                                               default:
+                                                       dst.at(l,j,k) = stat[l];
+                                                       //throw runtime_error(this->name + "undefined statistic: " + stat[l]);
+                                               }
+
+                                       }
+                               }
+                       }
+               }
+               else if (mode == "vert"){
+                       //dst.setGeometry(width,1);
+                       dst.setGeometry(width,n,iChannels,aChannels);
+                       histogram.setSampleCount(width);
+                       for (unsigned int k=0; k<channels; k++)
+                               for (unsigned int i=0; i<width; i++){
+                                       histogram.clear();
+                                       for (unsigned int j=0; j<height; j++)
+                                               histogram[src.at(i,j,k)]++;
+                                       for (unsigned int l=0; l<n; l++){
+                                               switch (stat[l]){
+                                               case 's':
+                                                       dst.at(i,l,k) = histogram.getSum();
+                                                       break;
+                                               case 'a':
+                                                       dst.at(i,l,k) = histogram.getMean();
+                                                       break;
+                                               default:
+                                                       throw runtime_error(this->name + "undefined statistic: " + stat[l]);
+                                               }
+
+                                       }
+                               }
+               }
+               else {
+                       throw runtime_error(this->name + ": geometry not horz or vert: " + mode);
+               }
+               //return dst;
+       };      
+       
+
+       
+protected:
+
+    
+};
+   
+}
+}
+
+#endif /*MARGINALSTATISTICOP_H_*/
diff --git a/drain/image/MathOpPack.h b/drain/image/MathOpPack.h
new file mode 100644 (file)
index 0000000..2c34b41
--- /dev/null
@@ -0,0 +1,399 @@
+/**
+
+    Copyright 2001 - 2010  Markus Peura, Finnish Meteorological Institute (First.Last@fmi.fi)
+
+
+    This file is part of Drain library for C++.
+
+    Drain is free software: you can redistribute it and/or modify
+    it under the terms of the GNU Lesser Public License as published by
+    the Free Software Foundation, either version 3 of the License, or
+    any later version.
+
+    Drain is distributed in the hope that it will be useful,
+    but WITHOUT ANY WARRANTY; without even the implied warranty of
+    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+    GNU Lesser Public License for more details.
+
+    You should have received a copy of the GNU General Public License
+    along with Drain.  If not, see <http://www.gnu.org/licenses/>.
+
+*/
+#ifndef MATHOP_H_
+#define MATHOP_H_
+
+#include <cmath>
+
+#include "SequentialImageOp.h"
+
+
+namespace drain
+{
+namespace image
+{
+
+/// A set of for mathematical operations.
+template <class T=unsigned char,class T2=unsigned char>
+class MathOp: public SequentialImageOp<T,T2>
+{
+public:
+
+       MathOp(const string &name = "MathOp",
+                       const string &description="Mathematical operator.",
+                       const string & parameterNames="scale,bias",
+                       const string & defaultValues = "1.0,0.0") :
+                               SequentialImageOp<T,T2>(name,description,parameterNames,defaultValues) {
+       }
+
+       MathOp(const string &name, const string &description, double scale, double bias = 0.0) :
+               SequentialImageOp<T,T2>(name,description,"scale,bias","1,0") {
+               this->setParameter("scale",scale);
+               this->setParameter("bias",bias);
+       }
+
+       virtual void setScale(float scale,float bias = 0.0){
+               this->setParameter("scale",scale);
+               this->setParameter("bias",bias);
+       }
+
+
+       void initialize() const {
+               scale = this->parameters.get("scale",scale);
+               bias  = this->parameters.get("bias",bias);
+       }
+
+       void filter(const Image<T> &src,Image<T2> &dst) const {
+