From cbad1fc2ce6a7d9f53183a2f042a5aa80131d223 Mon Sep 17 00:00:00 2001 From: Matt McCormick Date: Tue, 28 May 2024 16:07:22 -0400 Subject: [PATCH] COMP: Add itkIRCommon.h includes --- include/IRException.h | 75 + include/IRLog.h | 169 ++ include/IRMutex.h | 76 + include/IRMutexInterface.h | 74 + include/IRPath.h | 113 ++ include/IRTerminator.h | 130 ++ include/IRThread.h | 130 ++ include/IRThreadInterface.h | 203 +++ include/IRThreadPool.h | 217 +++ include/IRThreadStorage.h | 91 + include/IRTransaction.h | 216 +++ include/itkIRCommon.h | 36 +- include/itkIRDynamicArray.h | 350 ++++ include/itkIRTerminator.h | 64 + include/itkIRText.h | 324 ++++ include/itkIRUtils.h | 936 ++++++++++ include/itkInverseTransform.h | 134 ++ include/itkLegendrePolynomialTransform.h | 410 +++++ include/itkLegendrePolynomialTransform.hxx | 556 ++++++ include/itkNormalizeImageFilterWithMask.h | 111 ++ include/itkNormalizeImageFilterWithMask.hxx | 98 + include/itkNumericInverse.h | 183 ++ include/itkStatisticsImageFilterWithMask.h | 192 ++ include/itkStatisticsImageFilterWithMask.hxx | 390 ++++ itk-module.cmake | 6 +- src/CMakeLists.txt | 13 +- src/IRLog.cxx | 200 +++ src/IRMutex.cxx | 113 ++ src/IRMutexInterface.cxx | 65 + src/IRTerminator.cxx | 232 +++ src/IRThread.cxx | 196 ++ src/IRThreadInterface.cxx | 510 ++++++ src/IRThreadPool.cxx | 735 ++++++++ src/IRThreadStorage.cxx | 102 ++ src/IRTransaction.cxx | 183 ++ src/itkIRCommon.cxx | 1673 ++++++++++++++++++ src/itkIRText.cxx | 580 ++++++ src/itkIRUtils.cxx | 314 ++++ 38 files changed, 10179 insertions(+), 21 deletions(-) create mode 100644 include/IRException.h create mode 100644 include/IRLog.h create mode 100644 include/IRMutex.h create mode 100644 include/IRMutexInterface.h create mode 100644 include/IRPath.h create mode 100644 include/IRTerminator.h create mode 100644 include/IRThread.h create mode 100644 include/IRThreadInterface.h create mode 100644 include/IRThreadPool.h create mode 100644 include/IRThreadStorage.h create mode 100644 include/IRTransaction.h create mode 100644 include/itkIRDynamicArray.h create mode 100644 include/itkIRTerminator.h create mode 100644 include/itkIRText.h create mode 100644 include/itkIRUtils.h create mode 100644 include/itkInverseTransform.h create mode 100644 include/itkLegendrePolynomialTransform.h create mode 100644 include/itkLegendrePolynomialTransform.hxx create mode 100644 include/itkNormalizeImageFilterWithMask.h create mode 100644 include/itkNormalizeImageFilterWithMask.hxx create mode 100644 include/itkNumericInverse.h create mode 100644 include/itkStatisticsImageFilterWithMask.h create mode 100644 include/itkStatisticsImageFilterWithMask.hxx create mode 100644 src/IRLog.cxx create mode 100644 src/IRMutex.cxx create mode 100644 src/IRMutexInterface.cxx create mode 100644 src/IRTerminator.cxx create mode 100644 src/IRThread.cxx create mode 100644 src/IRThreadInterface.cxx create mode 100644 src/IRThreadPool.cxx create mode 100644 src/IRThreadStorage.cxx create mode 100644 src/IRTransaction.cxx create mode 100644 src/itkIRCommon.cxx create mode 100644 src/itkIRText.cxx create mode 100644 src/itkIRUtils.cxx diff --git a/include/IRException.h b/include/IRException.h new file mode 100644 index 0000000..814b6a9 --- /dev/null +++ b/include/IRException.h @@ -0,0 +1,75 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_exception.hxx +// Author : Pavel A. Koshevoy +// Created : Sun Sep 24 18:06:00 MDT 2006 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : an exception convenience class + +#ifndef THE_EXCEPTION_HXX_ +#define THE_EXCEPTION_HXX_ + +// system includes: +#include +#include +#include + + +//---------------------------------------------------------------- +// the_exception_t +// +class the_exception_t : public std::exception +{ +public: + the_exception_t(const char * description = NULL, + const char * file = NULL, + const unsigned int & line = 0) + { + std::ostringstream os; + + if (file != NULL) + { + os << file << ':' << line << " -- "; + } + + if (description != NULL) + { + os << description; + } + + what_ = os.str(); + } + + virtual ~the_exception_t() throw () + {} + + // virtual: + const char * what() const throw() + { return what_.c_str(); } + + // data: + std::string what_; +}; + + +#endif // THE_EXCEPTION_HXX_ diff --git a/include/IRLog.h b/include/IRLog.h new file mode 100644 index 0000000..7036743 --- /dev/null +++ b/include/IRLog.h @@ -0,0 +1,169 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_log.hxx +// Author : Pavel A. Koshevoy +// Created : Fri Mar 23 10:34:12 MDT 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A text log object -- behaves almost like a std::ostream. + +#ifndef THE_LOG_HXX_ +#define THE_LOG_HXX_ + +// system includes: +#include +#include +#include +#include + +// local includes: +#include "IRMutexInterface.h" +#include "itkIRUtils.h" + + +//---------------------------------------------------------------- +// the_log_t +// +class the_log_t +{ +protected: + void log_no_lock(std::ostream & (*f)(std::ostream &)); + +public: + the_log_t(); + virtual ~the_log_t(); + + virtual the_log_t & + operator << (std::ostream & (*f)(std::ostream &)); + + template + the_log_t & + operator << (const data_t & data) + { + the_lock_t lock(mutex_); + line_ << data; + return *this; + } + + std::streamsize precision(); + std::streamsize precision(std::streamsize n); + + std::ios::fmtflags flags() const; + std::ios::fmtflags flags(std::ios::fmtflags fmt); + + void setf(std::ios::fmtflags fmt); + void setf(std::ios::fmtflags fmt, std::ios::fmtflags msk); + void unsetf(std::ios::fmtflags fmt); + + void copyfmt(std::ostream & ostm); + + std::ostringstream line_; + mutable the_mutex_interface_t * mutex_; +}; + + +//---------------------------------------------------------------- +// the_null_log_t +// +class the_null_log_t : public the_log_t +{ +public: + // virtual: + the_log_t & operator << (std::ostream & (*)(std::ostream &)) + { return *this; } + + template + the_log_t & operator << (const data_t &) + { return *this; } +}; + +//---------------------------------------------------------------- +// the_stream_log_t +// +class the_stream_log_t : public the_log_t +{ +public: + the_stream_log_t(std::ostream & ostm): + ostm_(ostm) + {} + + // virtual: + the_log_t & operator << (std::ostream & (*f)(std::ostream &)) + { + the_lock_t lock(the_log_t::mutex_); + the_log_t::log_no_lock(f); + ostm_ << the_log_t::line_.str(); + the_log_t::line_.str(""); + return *this; + } + + template + the_log_t & operator << (const data_t & data) + { return the_log_t::operator << (data); } + + std::ostream & ostm_; +}; + + +//---------------------------------------------------------------- +// the_text_log_t +// +class the_text_log_t : public the_log_t +{ +public: + // virtual: + the_log_t & operator << (std::ostream & (*f)(std::ostream &)) + { + the_lock_t lock(the_log_t::mutex_); + the_log_t::log_no_lock(f); + text_ += the_log_t::line_.str(); + the_log_t::line_.str(""); + return *this; + } + + template + the_log_t & operator << (const data_t & data) + { return the_log_t::operator << (data); } + + inline std::string text() + { return text_; } + + std::string text_; +}; + +//---------------------------------------------------------------- +// null_log +// +extern the_null_log_t * null_log(); + +//---------------------------------------------------------------- +// cerr_log +// +extern the_stream_log_t * cerr_log(); + +//---------------------------------------------------------------- +// cout_log +// +extern the_stream_log_t * cout_log(); + + +#endif // THE_LOG_HXX_ diff --git a/include/IRMutex.h b/include/IRMutex.h new file mode 100644 index 0000000..bad00f7 --- /dev/null +++ b/include/IRMutex.h @@ -0,0 +1,76 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_boost_mutex.hxx +// Author : Pavel A. Koshevoy +// Created : Sat Oct 25 12:33:43 MDT 2008 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A thin wrapper for Boost mutex class. + +#ifndef THE_BOOST_MUTEX_HXX_ +#define THE_BOOST_MUTEX_HXX_ + +// local includes: +#include "thread/the_mutex_interface.hxx" + +// Boost includes: +#include + + +//---------------------------------------------------------------- +// the_boost_mutex_t +// +class the_boost_mutex_t : public the_mutex_interface_t +{ +public: + the_boost_mutex_t(); + + // the destructor is protected on purpose, + // see delete_this for details: + virtual ~the_boost_mutex_t(); + + // In order to avoid memory management problems with shared libraries, + // whoever provides this interface instance (via it's creator), has to + // provide a way to delete the instance as well. This will avoid + // issues with multiple-instances of C runtime libraries being + // used by the app and whatever libraries it links against that + // either use or provide this interface: + virtual void delete_this(); + + // the creation method: + static the_mutex_interface_t * create(); + + // virtual: + void lock(); + + // virtual: + void unlock(); + + // virtual: + bool try_lock(); + +private: + boost::mutex mutex_; +}; + + +#endif // THE_BOOST_MUTEX_HXX_ diff --git a/include/IRMutexInterface.h b/include/IRMutexInterface.h new file mode 100644 index 0000000..0f4e7ad --- /dev/null +++ b/include/IRMutexInterface.h @@ -0,0 +1,74 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_mutex_interface.hxx +// Author : Pavel A. Koshevoy +// Created : Sun Feb 18 16:12:00 MST 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : An abstract mutex class interface. + +#ifndef THE_MUTEX_INTERFACE_HXX_ +#define THE_MUTEX_INTERFACE_HXX_ + + +//---------------------------------------------------------------- +// the_mutex_interface_t +// +class the_mutex_interface_t +{ +protected: + // the destructor is protected on purpose, + // see delete_this for details: + virtual ~the_mutex_interface_t(); + +public: + // In order to avoid memory management problems with shared libraries, + // whoever provides this interface instance (via it's creator), has to + // provide a way to delete the instance as well. This will avoid + // issues with multiple-instances of C runtime libraries being + // used by the app and whatever libraries it links against that + // either use or provide this interface: + virtual void delete_this() = 0; + + // mutex controls: + virtual void lock() = 0; + virtual void unlock() = 0; + virtual bool try_lock() = 0; + + //---------------------------------------------------------------- + // creator_t + // + typedef the_mutex_interface_t *(*creator_t)(); + + // specify a thread creation method: + static void set_creator(creator_t creator); + + // create a new instance of a thread: + static the_mutex_interface_t * create(); + +protected: + // an abstract mutex creator: + static creator_t creator_; +}; + + +#endif // THE_MUTEX_INTERFACE_HXX_ diff --git a/include/IRPath.h b/include/IRPath.h new file mode 100644 index 0000000..816101c --- /dev/null +++ b/include/IRPath.h @@ -0,0 +1,113 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: nil -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : IRPath.h +// Author : Bradley C. Grimm +// Created : 2009/09/15 11:24 AM +// Copyright : (C) 2009 University of Utah +// License : GPLv2 +// Description : Helpful path manipulation tools. + +#ifndef __IR_PATH_HXX__ +#define __IR_PATH_HXX__ + +#if defined(WIN32) + #pragma warning ( disable : 4996 ) +#endif +#include "dirent.h" +#if defined(WIN32) + #pragma warning ( default : 4996 ) +#endif + +#include + +class IRPath +{ +public: + static the_text_t CleanPath(const the_text_t &path) + { + the_text_t directory = path; + CleanSlashes(directory); + if ( !directory.is_empty() && !directory.match_tail("/") ) + directory += the_text_t("/"); + return directory; + } + + static the_text_t DirectoryFromPath(const the_text_t &path) + { + the_text_t directory = path; + CleanSlashes(directory); + + directory = directory.splitAt('/', std::numeric_limits::max())[0]; + if ( !directory.is_empty() && !directory.match_tail("/") ) + return directory + "/"; + else + return directory; + } + + static the_text_t FilenameFromPath(the_text_t &path) + { + CleanSlashes(path); + + the_text_t directory; + if ( path.contains('/') ) + directory = path.splitAt('/', std::numeric_limits::max())[1]; + else + directory = path; + return directory; + } + + static void CleanSlashes(the_text_t &path) + { + if ( path.contains('\\') ) + path.replace('\\','/'); + } + + static std::list DirectoryContents(const the_text_t &path) + { + std::list files; + + DIR *dp; + struct dirent *dirp; + if((dp = opendir(path.text())) == NULL) { + cout << "Error(" << errno << ") opening " << path << endl; + return files; + } + + while ((dirp = readdir(dp)) != NULL) { + files.push_back( dirp->d_name ); + } + closedir(dp); + + // Remove . + files.pop_front(); + + // Remove .. + files.pop_front(); + + return files; + } + +}; + +#endif + + diff --git a/include/IRTerminator.h b/include/IRTerminator.h new file mode 100644 index 0000000..5999f1f --- /dev/null +++ b/include/IRTerminator.h @@ -0,0 +1,130 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_terminator.hxx +// Author : Pavel A. Koshevoy +// Created : Sun Sep 24 18:05:00 MDT 2006 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : a thread terminator convenience class + +#ifndef THE_TERMINATOR_HXX_ +#define THE_TERMINATOR_HXX_ + +#if defined(USE_THE_TERMINATORS) || defined(USE_ITK_TERMINATORS) +#define WRAP(x) +#else +#define WRAP(x) +#endif + +// system includes: +#include +#include + +// local includes: +#include "IRTerminator.h" +#include "itkIRUtils.h" + + +//---------------------------------------------------------------- +// the_terminator_t +// +class the_terminator_t +{ +public: + the_terminator_t(const char * id); + virtual ~the_terminator_t(); + + // terminator id accessor: + inline const std::string & id() const + { return id_; } + + // read-only accessor to the termination request flag: + inline const bool & termination_requested() const + { return termination_requested_; } + + // this function may be called periodically from any time consuming + // function -- in case the user has decided to terminate its execution: + inline void terminate_on_request() const + { + // if (termination_requested_ || should_terminate_immediately()) + if (termination_requested_) + { + throw_exception(); + } + } + + // this is a no-op for simple terminators, itk terminators will + // override this to turn on the process AbortGenerateData flag: + virtual void terminate(); + + // this function will throw an exception: + void throw_exception() const; + + // make sure there are no terminators left: + static bool verify_termination(); + +protected: + // consult the thread regarding whether termination has been requested + // for all terminators in the current thread: + static bool should_terminate_immediately(); + + // a list of active terminators per current thread: + static std::list & terminators(); + + // add/remove a terminator to/from the list of active terminators + // in the current thread: + static void add(the_terminator_t * terminator); + static void del(the_terminator_t * terminator); + + // id of this terminator: + std::string id_; + + // flag indicating that termination was requested explicitly + // for this terminator: + bool termination_requested_; +}; + +//---------------------------------------------------------------- +// the_terminators_t +// +class the_terminators_t +{ +public: + virtual ~the_terminators_t(); + + virtual void terminate(); + virtual bool verify_termination(); + + virtual void add(the_terminator_t * terminator); + virtual void del(the_terminator_t * terminator); + + // concurrent access controls: + virtual void lock() = 0; + virtual void unlock() = 0; + +private: + // the list of terminators: + std::list terminators_; +}; + + +#endif // THE_TERMINATOR_HXX_ diff --git a/include/IRThread.h b/include/IRThread.h new file mode 100644 index 0000000..2ff2cba --- /dev/null +++ b/include/IRThread.h @@ -0,0 +1,130 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_boost_thread.hxx +// Author : Pavel A. Koshevoy +// Created : Sat Oct 25 12:35:09 MDT 2008 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A thin wrapper for Boost thread class. + +#ifndef THE_BOOST_THREAD_HXX_ +#define THE_BOOST_THREAD_HXX_ + +// std includes: +#include +#include + +// local includes: +#include "IRTerminator.h" +#include "IRThreadInterface.h" + +// forward declarations: +class the_mutex_interface_t; +class the_thread_storage_t; + + +//---------------------------------------------------------------- +// the_boost_terminators_t +// +class the_boost_terminators_t : public the_terminators_t +{ +public: + // virtual: concurrent access controls: + void lock() { mutex_.lock(); } + void unlock() { mutex_.unlock(); } + +private: + mutable std::mutex mutex_; +}; + + +//---------------------------------------------------------------- +// the_boost_thread_t +// +// 1. the thread will not take ownership of the transactions. +// 2. the thread will take ownership of the mutex. +// +class the_boost_thread_t : public the_thread_interface_t +{ +private: + struct callable_t + { + callable_t(the_boost_thread_t * thread): + thread_(thread) + {} + + void operator()() + { + thread_->run(); + } + + private: + the_boost_thread_t * thread_; + }; + +public: + the_boost_thread_t(); + + // the destructor is protected on purpose, + // see delete_this for details: + virtual ~the_boost_thread_t(); + + // In order to avoid memory management problems with shared libraries, + // whoever provides this interface instance (via it's creator), has to + // provide a way to delete the instance as well. This will avoid + // issues with multiple-instances of C runtime libraries being + // used by the app and whatever libraries it links against that + // either use or provide this interface: + virtual void delete_this(); + + // the creation method: + static the_thread_interface_t * create() + { return new the_boost_thread_t(); } + + // the thread storage accessor: + static the_thread_storage_t & thread_storage(); + + // virtual: start the thread: + void start(); + + // virtual: + void wait(); + + // virtual: put the thread to sleep: + void take_a_nap(const unsigned long & microseconds); + + // virtual: accessor to the transaction terminators: + the_terminators_t & terminators(); + +protected: + // virtual: + void run(); + + // the boost thread: + std::thread * boost_thread_; + + // a list of active terminators for this thread: + the_boost_terminators_t terminators_; +}; + + +#endif // THE_BOOST_THREAD_HXX_ diff --git a/include/IRThreadInterface.h b/include/IRThreadInterface.h new file mode 100644 index 0000000..8670839 --- /dev/null +++ b/include/IRThreadInterface.h @@ -0,0 +1,203 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_thread_interface.hxx +// Author : Pavel A. Koshevoy +// Created : Fri Feb 16 09:20:00 MST 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : An abstract thread class interface. + +#ifndef THE_THREAD_INTERFACE_HXX_ +#define THE_THREAD_INTERFACE_HXX_ + +// system includes: +#include + +// local includes: +#include "IRTransaction.h" + +// forward declarations: +class the_mutex_interface_t; +class the_terminators_t; +class the_thread_pool_t; +struct the_thread_pool_data_t; + + +//---------------------------------------------------------------- +// the_thread_interface_t +// +// 1. the thread will not take ownership of the transactions. +// 2. the thread will take ownership of the mutex. +// +class the_thread_interface_t : public the_transaction_handler_t +{ + friend class the_thread_pool_t; + +protected: + // the destructor is protected on purpose, + // see delete_this for details: + virtual ~the_thread_interface_t(); + +public: + // In order to avoid memory management problems with shared libraries, + // whoever provides this interface instance (via it's creator), has to + // provide a way to delete the instance as well. This will avoid + // issues with multiple-instances of C runtime libraries being + // used by the app and whatever libraries it links against that + // either use or provide this interface: + virtual void delete_this() = 0; + + // the thread will own the mutex passed down to it, + // it will delete the mutex when the thread is deleted: + the_thread_interface_t(the_mutex_interface_t * mutex = NULL); + + //---------------------------------------------------------------- + // creator_t + // + typedef the_thread_interface_t *(*creator_t)(); + + // specify a thread creation method: + static void set_creator(creator_t creator); + + // create a new instance of a thread: + static the_thread_interface_t * create(); + + inline const unsigned int & id() const + { return id_; } + + // start the thread: + virtual void start() = 0; + + // wait for the thread to finish: + virtual void wait() = 0; + + // put the thread to sleep: + virtual void take_a_nap(const unsigned long & microseconds) = 0; + + // accessor to the transaction terminators: + virtual the_terminators_t & terminators() = 0; + + // give this thread an execution synchronization control: + void set_mutex(the_mutex_interface_t * mutex); + + // this controls whether the thread will voluntarily terminate + // once it runs out of transactions: + void set_idle_sleep_duration(bool enable, unsigned int microseconds = 10000); + + // schedule a transaction: + inline void add_transaction(the_transaction_t * transaction) + { push_back(transaction); } + + void push_back(the_transaction_t * transaction); + void push_front(the_transaction_t * transaction); + + // add transactions to the list: + void push_back(std::list & schedule); + + // NOTE: it is the responsibility of the caller to secure a mutex lock + // on this thread prior to checking whether the thread has any work left: + bool has_work() const; + + // schedule a transaction and start the thread: + void start(the_transaction_t * transaction); + + // abort the current transaction and clear pending transactions; + // transactionFinished will be emitted for the aborted transaction + // and the discarded pending transactions: + void stop(); + + // accessors to the thread "stopped" flag: + inline const bool & stopped() const + { return stopped_; } + + inline void set_stopped(bool stopped) + { stopped_ = stopped; } + + // clear all pending transactions, do not abort the current transaction: + void flush(); + + // this will call terminate_all for the terminators in this thread, + // but it will not stop the thread, so that new transactions may + // be scheduled while the old transactions are being terminated: + void terminate_transactions(); + + // terminate the current transactions and schedule a new transaction: + void stop_and_go(the_transaction_t * transaction); + void stop_and_go(std::list & schedule); + + // flush the current transactions and schedule a new transaction: + void flush_and_go(the_transaction_t * transaction); + void flush_and_go(std::list & schedule); + + // execute the scheduled transactions, return true if all + // transactions had been executed successfully: + virtual bool work(); + + // virtual: default transaction communication handlers: + void handle(the_transaction_t * transaction, the_transaction_t::state_t s); + void blab(const char * message) const; + + // transaction callback accessor: + inline void + set_thread_pool_cb(the_thread_pool_t * pool, + the_thread_pool_data_t * cb_data) + { + thread_pool_ = pool; + thread_pool_cb_data_ = cb_data; + } + + // mutex accessor: + inline the_mutex_interface_t * mutex() + { return mutex_; } + +protected: + // an abstract creator of threads: + static creator_t creator_; + + // FIXME: this may have been a bad idea: + unsigned int id_; + + // thread synchronization control: + the_mutex_interface_t * mutex_; + + // execution control flag: + bool stopped_; + + // these attributes control the thread behavior once all transactions + // have been processed. If sleep_when_idle_ flag is set to true, the + // thread will put itself to sleep instead of terminating: + bool sleep_when_idle_; + unsigned int sleep_microsec_; + + // currently executing transaction: + the_transaction_t * active_transaction_; + + // scheduled transactions: + std::list transactions_; + + // the transaction callback: + the_thread_pool_t * thread_pool_; + the_thread_pool_data_t * thread_pool_cb_data_; +}; + + +#endif // THE_THREAD_INTERFACE_HXX_ diff --git a/include/IRThreadPool.h b/include/IRThreadPool.h new file mode 100644 index 0000000..0b616ef --- /dev/null +++ b/include/IRThreadPool.h @@ -0,0 +1,217 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_thread_pool.hxx +// Author : Pavel A. Koshevoy +// Created : Wed Feb 21 08:30:00 MST 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A thread pool class. + +#ifndef THE_THREAD_POOL_HXX_ +#define THE_THREAD_POOL_HXX_ + +// system includes: +#include +#include + +// local includes: +#include "IRThreadInterface.h" +#include "IRTransaction.h" + +// forward declarations: +class the_mutex_interface_t; +class the_thread_pool_t; + + +//---------------------------------------------------------------- +// transaction_wrapper_t +// +class the_transaction_wrapper_t +{ +public: + the_transaction_wrapper_t(const unsigned int & num_parts, + the_transaction_t * transaction); + ~the_transaction_wrapper_t(); + + static void notify_cb(void * data, + the_transaction_t * t, + the_transaction_t::state_t s); + + bool notify(the_transaction_t * t, + the_transaction_t::state_t s); + +private: + the_transaction_wrapper_t(); + the_transaction_wrapper_t & operator = (const the_transaction_wrapper_t &); + + the_mutex_interface_t * mutex_; + + the_transaction_t::notify_cb_t cb_; + void * cb_data_; + + const unsigned int num_parts_; + unsigned int notified_[the_transaction_t::DONE_E + 1]; +}; + + +//---------------------------------------------------------------- +// the_thread_pool_data_t +// +struct the_thread_pool_data_t +{ + the_thread_pool_data_t(): + parent_(NULL), + thread_(NULL), + id_(~0u) + {} + + ~the_thread_pool_data_t() + { thread_->delete_this(); } + + the_thread_pool_t * parent_; + the_thread_interface_t * thread_; + unsigned int id_; +}; + + +//---------------------------------------------------------------- +// the_thread_pool_t +// +class the_thread_pool_t : public the_transaction_handler_t +{ + friend class the_thread_interface_t; + +public: + the_thread_pool_t(unsigned int num_threads); + virtual ~the_thread_pool_t(); + + // start the threads: + void start(); + + // this controls whether the thread will voluntarily terminate + // once it runs out of transactions: + void set_idle_sleep_duration(bool enable, unsigned int microseconds = 10000); + + // schedule a transaction: + // NOTE: when multithreaded is true, the transaction will be scheduled + // N times, where N is the number of threads in the pool. + // This means the transaction will be executed by N threads, so + // the transaction has to support concurrent execution internally. + virtual void push_front(the_transaction_t * transaction, + bool multithreaded = false); + + virtual void push_back(the_transaction_t * transaction, + bool multithreaded = false); + + virtual void push_back(std::list & schedule, + bool multithreaded = false); + + // split the work among the threads: + void pre_distribute_work(); + + // check whether the thread pool has any work left: + bool has_work() const; + + // schedule a transaction and start the thread: + virtual void start(the_transaction_t * transaction, + bool multithreaded = false); + + // abort the current transaction and clear pending transactions; + // transactionFinished will be emitted for the aborted transaction + // and the discarded pending transactions: + void stop(); + + // wait for all threads to finish: + void wait(); + + // clear all pending transactions, do not abort the current transaction: + void flush(); + + // this will call terminate_all for the terminators in this thread, + // but it will not stop the thread, so that new transactions may + // be scheduled while the old transactions are being terminated: + void terminate_transactions(); + + // terminate the current transactions and schedule a new transaction: + void stop_and_go(the_transaction_t * transaction, + bool multithreaded = false); + void stop_and_go(std::list & schedule, + bool multithreaded = false); + + // flush the current transactions and schedule a new transaction: + void flush_and_go(the_transaction_t * transaction, + bool multithreaded = false); + void flush_and_go(std::list & schedule, + bool multithreaded = false); + + // virtual: default transaction communication handlers: + void handle(the_transaction_t * transaction, the_transaction_t::state_t s); + void blab(const char * message) const; + + // access control: + inline the_mutex_interface_t * mutex() + { return mutex_; } + + inline const unsigned int & pool_size() const + { return pool_size_; } + +private: + // intentionally disabled: + the_thread_pool_t(const the_thread_pool_t &); + the_thread_pool_t & operator = (const the_thread_pool_t &); + +protected: + // thread callback handler: + virtual void handle_thread(the_thread_pool_data_t * data); + + // helpers: + inline the_thread_interface_t * thread(unsigned int id) const + { + assert(id < pool_size_); + return pool_[id].thread_; + } + + void no_lock_flush(); + void no_lock_terminate_transactions(); + void no_lock_push_front(the_transaction_t * transaction, bool multithreaded); + void no_lock_push_back(the_transaction_t * transaction, bool multithreaded); + void no_lock_push_back(std::list & schedule, bool mt); + + // thread synchronization control: + the_mutex_interface_t * mutex_; + + // the thread pool: + the_thread_pool_data_t * pool_; + unsigned int pool_size_; + + // the working threads: + std::list busy_; + + // the waiting threads: + std::list idle_; + + // scheduled transactions: + std::list transactions_; +}; + + +#endif // THE_THREAD_POOL_HXX_ diff --git a/include/IRThreadStorage.h b/include/IRThreadStorage.h new file mode 100644 index 0000000..e2444d6 --- /dev/null +++ b/include/IRThreadStorage.h @@ -0,0 +1,91 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_thread_storage.hxx +// Author : Pavel A. Koshevoy +// Created : Fri Jan 2 09:30:00 MDT 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A thread storage abstract interface class. + +#ifndef THE_THREAD_STORAGE_HXX_ +#define THE_THREAD_STORAGE_HXX_ + +// forward declarations: +class the_thread_interface_t; +class the_terminators_t; + +//---------------------------------------------------------------- +// the_thread_observer_t +// +class the_thread_observer_t +{ +public: + the_thread_observer_t(the_thread_interface_t & thread): + thread_(thread) + {} + + the_thread_interface_t & thread_; +}; + +//---------------------------------------------------------------- +// the_thread_storage_t +// +class the_thread_storage_t +{ +public: + virtual ~the_thread_storage_t() {} + + // check whether the thread storage has been initialized: + virtual bool is_ready() const = 0; + + // check whether the thread has been stopped: + virtual bool thread_stopped() const = 0; + + // terminator access: + virtual the_terminators_t & terminators() = 0; + + // thread id: + virtual unsigned int thread_id() const = 0; +}; + +//---------------------------------------------------------------- +// the_thread_storage_provider_t +// +typedef the_thread_storage_t&(*the_thread_storage_provider_t)(); + +//---------------------------------------------------------------- +// set_the_thread_storage_provider +// +// Set the new thread storage provider, return the old provider. +// +extern the_thread_storage_provider_t +set_the_thread_storage_provider(the_thread_storage_provider_t p); + +//---------------------------------------------------------------- +// the_thread_storage +// +// Thread storage accessors. +// +extern the_thread_storage_t & the_thread_storage(); + + +#endif // THE_THREAD_STORAGE_HXX_ diff --git a/include/IRTransaction.h b/include/IRTransaction.h new file mode 100644 index 0000000..acbad3a --- /dev/null +++ b/include/IRTransaction.h @@ -0,0 +1,216 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_transaction.hxx +// Author : Pavel A. Koshevoy +// Created : Fri Feb 16 09:52:00 MST 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A thread transaction class. + +#ifndef THE_TRANSACTION_HXX_ +#define THE_TRANSACTION_HXX_ + +// system includes: +#include +#include + +// the includes: +#include +#include + +// forward declarations: +class the_thread_interface_t; +class the_transaction_handler_t; + + +//---------------------------------------------------------------- +// the_transaction_t +// +// NOTE: the transaction will not take ownership of the mutex. +// +class the_transaction_t +{ +public: + the_transaction_t(); + virtual ~the_transaction_t(); + + // execute the transaction: + virtual void execute(the_thread_interface_t * thread) = 0; + + // transaction execution state: + typedef enum { + PENDING_E, + SKIPPED_E, + STARTED_E, + ABORTED_E, + DONE_E + } state_t; + + inline state_t state() const + { return state_; } + + inline void set_state(const state_t & s) + { state_ = s; } + + inline bool done() const + { return state_ == DONE_E; } + + //---------------------------------------------------------------- + // notify_cb_t + // + typedef void(*notify_cb_t)(void *, the_transaction_t *, state_t s); + + inline notify_cb_t notify_cb() const + { return notify_cb_; } + + inline void set_notify_cb(notify_cb_t cb, void * cb_data) + { + notify_cb_ = cb; + notify_cb_data_ = cb_data; + } + + //---------------------------------------------------------------- + // status_cb_t + // + // NOTE: if status is NULL, it means the transaction is requesting + // a callback from the main thread. This could be used by the + // transaction for some GUI user interaction (such as finding replacement + // tiles for a mosaic, etc...) + // + typedef void(*status_cb_t)(void *, the_transaction_t *, const char * status); + + inline status_cb_t status_cb() const + { return status_cb_; } + + inline void set_status_cb(status_cb_t cb, void * cb_data) + { + status_cb_ = cb; + status_cb_data_ = cb_data; + } + + // notify the transaction about a change in it's state: + virtual void notify(the_transaction_handler_t * handler, + state_t s, + const char * message = NULL); + + // helper: + virtual void blab(the_transaction_handler_t * handler, + const char * message); + + // FIXME: this is a relic: + inline static const the_text_t tr(const char * text) + { return the_text_t(text); } + + inline static const the_text_t & tr(const the_text_t & text) + { return text; } + + // if a transaction needs to have something executed in the main thread + // context it should call callback_request. This requires that a valid + // status callback is set for this transaction, and that the status + // callback will acknowledge the main thread callback request. This + // call will block until the callback is executed in the main thread, + // and the request is removed: + bool callback_request(); + + // if a transaction needs to have something executed in the main thread + // context it should override this callback. This method will be executed + // in the main thread context. The default implementation will simply + // remove the request: + virtual void callback(); + +protected: + // when requesting a callback from the main thread + // the status of request will be set to WAITING_E: + typedef enum { + NOTHING_E, + WAITING_E + } request_t; + + // synchronization control: + the_mutex_interface_t * mutex_; + + // status of callback request: + request_t request_; + + // current state of the transaction: + state_t state_; + +public: + // the callbacks: + notify_cb_t notify_cb_; + void * notify_cb_data_; + + status_cb_t status_cb_; + void * status_cb_data_; +}; + +//---------------------------------------------------------------- +// operator << +// +extern std::ostream & +operator << (std::ostream & so, const the_transaction_t::state_t & state); + + +//---------------------------------------------------------------- +// the_transaction_handler_t +// +class the_transaction_handler_t +{ +public: + virtual ~the_transaction_handler_t() {} + + virtual void handle(the_transaction_t * transaction, + the_transaction_t::state_t s) = 0; + virtual void blab(const char * message) const = 0; +}; + + +//---------------------------------------------------------------- +// the_transaction_log_t +// +class the_transaction_log_t : public the_log_t +{ +public: + the_transaction_log_t(the_transaction_handler_t * handler): + handler_(handler) + {} + + // virtual: + the_log_t & operator << (std::ostream & (*f)(std::ostream &)) + { + the_log_t::operator << (f); + std::string text(the_log_t::line_.str()); + the_log_t::line_.str(""); + handler_->blab(text.c_str()); + return *this; + } + + template + the_log_t & operator << (const data_t & data) + { return the_log_t::operator << (data); } + +private: + the_transaction_handler_t * handler_; +}; + + +#endif // THE_TRANSACTION_HXX_ diff --git a/include/itkIRCommon.h b/include/itkIRCommon.h index c39230f..2231c76 100644 --- a/include/itkIRCommon.h +++ b/include/itkIRCommon.h @@ -28,18 +28,18 @@ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA // image preprocessing, convenience wrappers for // ITK file and ITK filters. -#ifndef COMMON_HXX_ -#define COMMON_HXX_ +#ifndef itkIRCommon_h +#define itkIRCommon_h // local includes: -#include "utils/the_text.hxx" -#include "utils/the_utils.hxx" -#include "itk/itk_terminator.hxx" -#include "itk/itkNormalizeImageFilterWithMask.h" -#include "itk/itkLegendrePolynomialTransform.h" -#include "thread/the_boost_thread.hxx" -#include "thread/the_transaction.hxx" -#include "thread/the_thread_pool.hxx" +#include "itkIRText.h" +#include "itkIRUtils.h" +#include "itkIRTerminator.h" +#include "itkNormalizeImageFilterWithMask.h" +#include "itkLegendrePolynomialTransform.h" +#include "IRThread.h" +#include "IRTransaction.h" +#include "IRThreadPool.h" //#include "utils/AsyncMosaicSave.h" // system includes: @@ -99,11 +99,9 @@ using std::flush; #include #include #include -#include +#include -#include - #include "IRPath.h" extern const char *g_Version; @@ -2190,7 +2188,7 @@ double double final_sb = 0.0; unsigned long int final_pixels = 0; - #pragma omp parallel for + // #pragma omp parallel for for(int iThread = 0; iThread < (int)list_fi_thread_itex.size(); iThread++) { //std::vector indexList = threadIndexList[iThread]; @@ -2698,7 +2696,7 @@ void // calculate the bounding boxes: std::vector vector_in(in.size()); vector_in.assign(in.begin(), in.end()); - #pragma omp parallel for + // #pragma omp parallel for for (int i = 0; i < (int)vector_in.size(); i++) { the_text_t *iter = &vector_in[i]; @@ -2811,7 +2809,7 @@ bool // calculate the bounding boxes bool ok = true; - #pragma omp parallel for + // #pragma omp parallel for for (int i = 0; i < (int)num_images; i++) { ok = ok & calc_tile_mosaic_bbox(xform[i], @@ -4978,7 +4976,7 @@ void typedef itk::RGBPixel composite_pixel_t; typedef itk::Image composite_image_t; - typedef itk::ComposeRGBImageFilter + typedef itk::ComposeImageFilter compose_filter_t; compose_filter_t::Pointer composer = compose_filter_t::New(); @@ -5278,7 +5276,7 @@ void typedef itk::RGBPixel composite_pixel_t; typedef itk::Image composite_image_t; - typedef itk::ComposeRGBImageFilter + typedef itk::ComposeImageFilter compose_filter_t; typename compose_filter_t::Pointer composer = compose_filter_t::New(); @@ -7119,4 +7117,4 @@ inline static void increment_minor_progress(double amount = 1) -#endif // COMMON_HXX_ +#endif // itkIRCommon_h diff --git a/include/itkIRDynamicArray.h b/include/itkIRDynamicArray.h new file mode 100644 index 0000000..62be2a0 --- /dev/null +++ b/include/itkIRDynamicArray.h @@ -0,0 +1,350 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_dynamic_array.hxx +// Author : Pavel A. Koshevoy +// Created : Fri Oct 31 17:16:25 MDT 2003 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : Implementation of a dynamically resizable array +// that grows automatically. + +#ifndef THE_DYNAMIC_ARRAY_HXX_ +#define THE_DYNAMIC_ARRAY_HXX_ + +// system includes: +#include +#include +#include + +// forward declarations: +template class the_dynamic_array_t; + +#undef min +#undef max + + +//---------------------------------------------------------------- +// the_dynamic_array_ref_t +// +template +class the_dynamic_array_ref_t +{ +public: + the_dynamic_array_ref_t(the_dynamic_array_t & array, + size_t index = 0): + array_(array), + index_(index) + {} + + inline the_dynamic_array_ref_t & operator << (const T & elem) + { + array_[index_++] = elem; + return *this; + } + +private: + // reference to the array: + the_dynamic_array_t & array_; + + // current index into the array: + size_t index_; +}; + + +//---------------------------------------------------------------- +// the_dynamic_array_t +// +template +class the_dynamic_array_t +{ +public: + the_dynamic_array_t(): + array_(NULL), + page_size_(16), + size_(0), + init_value_() + {} + + the_dynamic_array_t(const size_t & init_size): + array_(NULL), + page_size_(init_size), + size_(0), + init_value_() + {} + + the_dynamic_array_t(const size_t & init_size, + const size_t & page_size, + const T & init_value): + array_(NULL), + page_size_(page_size), + size_(0), + init_value_(init_value) + { + resize(init_size); + } + + // copy constructor: + the_dynamic_array_t(const the_dynamic_array_t & a): + array_(NULL), + page_size_(0), + size_(0), + init_value_(a.init_value_) + { + (*this) = a; + } + + // destructor: + ~the_dynamic_array_t() + { + clear(); + } + + // remove all contents of this array: + void clear() + { + size_t num = num_pages(); + for (size_t i = 0; i < num; i++) + { + delete (*array_)[i]; + } + + delete array_; + array_ = NULL; + + size_ = 0; + } + + // the assignment operator: + the_dynamic_array_t & operator = (const the_dynamic_array_t & array) + { + clear(); + + page_size_ = array.page_size_; + init_value_ = array.init_value_; + + resize(array.size_); + for (size_t i = 0; i < size_; i++) + { + (*this)[i] = array[i]; + } + + return *this; + } + + // resize the array, all contents will be preserved: + void resize(const size_t & new_size) + { + // bump the current size value: + size_ = new_size; + + // do nothing if resizing is unnecessary: + if (size_ <= max_size()) return; + + // we'll have to do something about the existing data: + size_t old_num_pages = num_pages(); + size_t new_num_pages = + std::max((size_t)(2 * old_num_pages), + (size_t)(1 + size_ / page_size_)); + + // create a new array: + std::vector< std::vector * > * new_array = + new std::vector< std::vector * >(new_num_pages); + + // shallow-copy the old content: + for (size_t i = 0; i < old_num_pages; i++) + { + (*new_array)[i] = (*array_)[i]; + } + + // initialize the new pages: + for (size_t i = old_num_pages; i < new_num_pages; i++) + { + (*new_array)[i] = new std::vector(page_size_); + for (size_t j = 0; j < page_size_; j++) + { + (*(*new_array)[i])[j] = init_value_; + } + } + + // get rid of the old array: + delete array_; + + // put the new array in place of the old array: + array_ = new_array; + } + + // the size of this array: + inline const size_t & size() const + { return size_; } + + inline const size_t & page_size() const + { return page_size_; } + + // maximum usable size of the array that does not require resizing the array: + inline size_t max_size() const + { return num_pages() * page_size_; } + + // number of pages currently allocated: + inline size_t num_pages() const + { return (array_ == NULL) ? 0 : array_->size(); } + + inline const T * page(const size_t & page_index) const + { return &((*(*array_)[page_index])[0]); } + + inline T * page(const size_t & page_index) + { return &((*(*array_)[page_index])[0]); } + + // return either first or last index into the array: + inline size_t end_index(bool last) const + { + if (last == false) return 0; + return size_ - 1; + } + + // return either first or last element in the array: + inline const T & end_elem(bool last) const + { return elem(end_index(last)); } + + inline T & end_elem(bool last) + { return elem(end_index(last)); } + + inline const T & front() const + { return end_elem(false); } + + inline T & front() + { return end_elem(false); } + + inline const T & back() const + { return end_elem(true); } + + inline T & back() + { return end_elem(true); } + + // non-const accessors: + inline T & elem(const size_t i) + { + if (i >= size_) resize(i + 1); + return (*(*array_)[i / page_size_])[i % page_size_]; + } + + inline T & operator [] (const size_t & i) + { return elem(i); } + + // const accessors: + inline const T & elem(const size_t & i) const + { return (*((*array_)[i / page_size_]))[i % page_size_]; } + + inline const T & operator [] (const size_t & i) const + { return elem(i); } + + // this is usefull for filling-in the array: + the_dynamic_array_ref_t operator << (const T & elem) + { + (*this)[0] = elem; + return the_dynamic_array_ref_t(*this, 1); + } + + // grow the array by one and insert a new element at the tail: + inline void push_back(const T & elem) + { (*this)[size_] = elem; } + + inline void append(const T & elem) + { push_back(elem); } + + // return the index of the first occurrence of a given element in the array: + size_t index_of(const T & element) const + { + for (size_t i = 0; i < size_; i++) + { + if (!(elem(i) == element)) continue; + return i; + } + + return ~0u; + } + + // check whether this array contains a given element: + inline bool has(const T & element) const + { return index_of(element) != ~0u; } + + // remove an element from the array: + bool remove(const T & element) + { + size_t idx = index_of(element); + if (idx == ~0u) return false; + + for (size_t i = idx + 1; i < size_; i++) + { + elem(i - 1) = elem(i); + } + + size_--; + return true; + } + + void assign(const size_t & size, const T & element) + { + resize(size); + for (size_t i = 0; i < size; i++) + { + elem(i) = element; + } + } + + // for debugging, dumps this list into a stream: + void dump(std::ostream & strm) const + { + strm << "the_dynamic_array_t(" << (void *)this << ") {\n"; + for (size_t i = 0; i < size_; i++) + { + strm << elem(i) << std::endl; + } + strm << '}'; + } + +protected: + // an array of pointers to arrays (pages) of data: + std::vector< std::vector *> * array_; + + // page size: + size_t page_size_; + + // current array size: + size_t size_; + + // init value used when resizing the array: + T init_value_; +}; + +//---------------------------------------------------------------- +// operator << +// +template +std::ostream & +operator << (std::ostream & s, const the_dynamic_array_t & a) +{ + a.dump(s); + return s; +} + + +#endif // THE_DYNAMIC_ARRAY_HXX_ diff --git a/include/itkIRTerminator.h b/include/itkIRTerminator.h new file mode 100644 index 0000000..a88a0d5 --- /dev/null +++ b/include/itkIRTerminator.h @@ -0,0 +1,64 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itk_terminator.hxx +// Author : Pavel A. Koshevoy +// Created : Sun Sep 24 18:16:00 MDT 2006 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : a ITK filter terminator convenience class + +#ifndef ITK_TERMINATOR_HXX_ +#define ITK_TERMINATOR_HXX_ + +// local includes: +#include "IRTerminator.h" + +#ifdef USE_ITK_TERMINATORS +#define itk_terminator_t the_terminator_t + + +//---------------------------------------------------------------- +// terminator_t +// +template +class terminator_t : public the_terminator_t +{ +public: + terminator_t(process_t * proc): + itk_terminator_t(proc->GetNameOfClass()), + process_(proc) + {} + + // virtual: + void terminate() + { + process_->AbortGenerateDataOn(); + itk_terminator_t::terminate(); + } + +private: + typename process_t::Pointer process_; +}; + + +#endif // USE_ITK_TERMINATORS +#endif // ITK_TERMINATOR_HXX_ diff --git a/include/itkIRText.h b/include/itkIRText.h new file mode 100644 index 0000000..3a03a01 --- /dev/null +++ b/include/itkIRText.h @@ -0,0 +1,324 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_text.hxx +// Author : Pavel A. Koshevoy +// Created : Sun Aug 29 14:53:00 MDT 2004 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : text convenience class. + +#ifndef itkIRText_h +#define itkIRText_h + +// system includes: +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +//---------------------------------------------------------------- +// the_text_t +// +class the_text_t +{ +public: + the_text_t(const char * text = ""); + the_text_t(const char * text, const size_t & size); + the_text_t(const the_text_t & text); + the_text_t(const std::list & text); + ~the_text_t(); + + // assignment operator: + inline the_text_t & operator = (const the_text_t & text) + { + if (this != &text) + { + assign(text.text_, text.size_); + } + + return *this; + } + + // clear the string: + inline void clear() + { + delete [] text_; + text_ = NULL; + size_ = 0; + } + + // shorthand: + inline bool is_empty() const + { return size_ == 0; } + + // assign a new string to this text: + inline void assign(const char * text) + { assign(text, strlen(text)); } + + void assign(const char * text, const size_t & text_size); + + // append a new string to this text: + inline void append(const char * text) + { append(text, strlen(text)); } + + void append(const char * text, const size_t & text_size); + + // replace + void replace(const char find, const char replace); + + // equality/inequality tests: + inline bool operator == (const the_text_t & text) const + { return ((size_ == text.size_) && (strcmp(text_, text.text_) == 0)); } + + inline bool operator != (const the_text_t & text) const + { return !(*this == text); } + + inline bool operator == (const char * text) const + { return (*this == the_text_t(text)); } + + inline bool operator != (const char * text) const + { return !(*this == text); } + + inline bool operator < (const the_text_t & text) const + { return (strcmp(text_, text.text_) < 0); } + + inline bool operator > (const the_text_t & text) const + { return (strcmp(text_, text.text_) > 0); } + + // arithmetic: + inline the_text_t & operator += (const the_text_t & text) + { + append(text.text_, text.size_); + return *this; + } + + inline the_text_t operator + (const the_text_t & text) const + { + the_text_t text_sum(*this); + text_sum += text; + return text_sum; + } + + // access operators: + template + inline const char & operator [] (const index_t & index) const + { return text_[index]; } + + template + inline char & operator [] (const index_t & index) + { return text_[index]; } + + // accessors: + inline const char * text() const + { return text_; } + + inline const size_t & size() const + { return size_; } + + // conversion operator: + inline operator const char * () const + { return text_; } + + inline static the_text_t pad(const char * str, + const size_t width = 0, + const char pad_char = ' ', + const bool pad_left = true) + { + the_text_t txt(str); + + if (width > txt.size()) + { + the_text_t padding; + padding.fill(pad_char, width - txt.size()); + txt = pad_left ? padding + txt : txt + padding; + } + + return txt; + } + + // helpers: + template + static the_text_t number(const number_t & number, + const size_t width = 0, + const char pad_char = ' ', + const bool pad_left = true) + { + std::ostringstream os; + os << number; + + std::string str = os.str(); + return pad(str.c_str(), width, pad_char, pad_left); + } + + inline static the_text_t number(const size_t & number, + const size_t width = 0, + const char pad_char = ' ', + const bool pad_left = true) + { +#ifdef _WIN32 +#ifndef snprintf +#define snprintf _snprintf_s +#endif +#endif + + static char buffer[256]; + snprintf(buffer, sizeof(buffer), "%llu", (long long unsigned int)(number)); + return pad(buffer, width, pad_char, pad_left); + } + + short int toShort(bool * ok = 0, int base = 10) const; + unsigned short int toUShort(bool * ok = NULL, int base = 10) const; + + int toInt(bool * ok = NULL, int base = 10) const; + unsigned int toUInt(bool * ok = NULL, int base = 10) const; + + long int toLong(bool * ok = NULL, int base = 10) const; + unsigned long int toULong(bool * ok = NULL, int base = 10) const; + + float toFloat(bool * ok = NULL) const; + double toDouble(bool * ok = NULL) const; + + void to_ascii(); + void to_lower(); + void to_upper(); + void fill(const char & c, const size_t size); + + void fill(const char & c) + { fill(c, size_); } + + // return true if the given text matches the head/tail of this text: + bool match_head(const the_text_t & t, bool ignore_case = false) const; + bool match_tail(const the_text_t & t, bool ignore_case = false) const; + + bool match_text(const the_text_t & t, + const size_t & index, + bool ignore_case = false) const; + + // remove leading/tailing white space, replace internal white space + // with a single space: + the_text_t simplify_ws() const; + + // split the text into a set of tokens, return the number of tokens: + size_t split(std::vector & tokens, + const char & separator, + const bool & empty_ok = false) const; + + // splitAt splits string into two parts upon a character. + std::vector splitAt(const char split_char, + unsigned int n) const; + + // count the number of occurrences of a given symbol in the text: + size_t contains(const char & symbol) const; + + // extract a portion of the string: + void extract(the_text_t & to, + const size_t & from, + const size_t & size) const + { + assert(from + size < size_ + 1); + to.assign(&text_[from], size); + } + + inline the_text_t extract(const size_t & from, + const size_t & size) const + { + the_text_t to; + extract(to, from, size); + return to; + } + + inline the_text_t reverse() const + { + the_text_t rev(*this); + for (size_t i = 0; i < size_; i++) + { + rev.text_[i] = text_[size_ - i - 1]; + } + + return rev; + } + + inline the_text_t cut(const char & separator, + size_t f0, + size_t f1 = 0) const + { + const char sep_str[2] = { separator, '\0' }; + + std::vector fields; + split(fields, separator, true); + size_t num_fields = fields.size(); + + if (f1 < f0) + { + f1 = f0; + } + else if (f1 >= num_fields) + { + f1 = num_fields - 1; + } + + the_text_t out; + for (size_t f = f0; f <= f1; f++) + { + out += fields[f]; + if (f + 1 <= f1) + { + out += sep_str; + } + } + + return out; + } + +private: + // the text itself: + char * text_; + + // the length of the text: + size_t size_; +}; + +extern std::ostream & +operator << (std::ostream & out, const the_text_t & text); + +extern std::istream & +operator >> (std::istream & in, the_text_t & text); + +extern std::istream & +getline(std::istream & in, the_text_t & text); + +//---------------------------------------------------------------- +// to_binary +// +// return a 0 and 1 string representation of a byte +// +extern the_text_t +to_binary(const unsigned char & byte, bool lsb_first = true); + + +#endif // itkIRText_h diff --git a/include/itkIRUtils.h b/include/itkIRUtils.h new file mode 100644 index 0000000..a72f678 --- /dev/null +++ b/include/itkIRUtils.h @@ -0,0 +1,936 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_utils.hxx +// Author : Pavel A. Koshevoy +// Created : 2006/04/15 11:25:00 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : utility functions for working with arrays, +// lists, numbers, angles, etc... + +#ifndef THE_UTILS_HXX_ +#define THE_UTILS_HXX_ + +// system includes: +#ifndef _USE_MATH_DEFINES +#define _USE_MATH_DEFINES +#endif + +#ifndef NOMINMAX +#define NOMINMAX +#endif + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// local includes: +#include "itkIRDynamicArray.h" + + +//---------------------------------------------------------------- +// array2d +// +#define array2d( T ) std::vector > + +//---------------------------------------------------------------- +// array3d +// +#define array3d( T ) std::vector > > + +//---------------------------------------------------------------- +// TWO_PI +// +static const double TWO_PI = 2.0 * M_PI; + + +//---------------------------------------------------------------- +// clamp_angle +// +inline double +clamp_angle(const double & absolute_angle) +{ + double a = fmod(absolute_angle + TWO_PI, TWO_PI); + assert(a >= 0.0 && a < TWO_PI); + return a; +} + +//---------------------------------------------------------------- +// calc_angle +// +inline double +calc_angle(const double & x, + const double & y, + const double & reference_angle = 0.0) +{ + return clamp_angle(fmod(atan2(y, x) + TWO_PI, TWO_PI) - + fmod(reference_angle, TWO_PI)); +} + +//---------------------------------------------------------------- +// sleep_msec +// +extern void sleep_msec(size_t msec); + +//---------------------------------------------------------------- +// drand +// +inline static double drand() +{ return double(rand()) / double(RAND_MAX); } + +//---------------------------------------------------------------- +// integer_power +// +template +inline scalar_t +integer_power(scalar_t x, size_t p) +{ + scalar_t result = scalar_t(1); + while (p != 0u) + { + if (p & 1) result *= x; + x *= x; + p >>= 1; + } + + return result; +} + +//---------------------------------------------------------------- +// closest_power_of_two_larger_than_given +// +template +inline scalar_t +closest_power_of_two_larger_than_given(const scalar_t & given) +{ + size_t n = sizeof(given) * 8; + scalar_t closest = scalar_t(1); + for (size_t i = 0; + (i < n) && (closest < given); + i++, closest *= scalar_t(2)) {} + + return closest; +} + + +//---------------------------------------------------------------- +// divide +// +template +data_t +divide(const data_t & numerator, const data_t & denominator) +{ + static data_t zero = data_t(0); + return (denominator != zero) ? (numerator / denominator) : zero; +} + +//---------------------------------------------------------------- +// clear_stack +// +template +void +clear_stack(std::stack & s) +{ + while (!s.empty()) + { + s.pop(); + } +} + +//---------------------------------------------------------------- +// resize +// +template +void +resize(array2d(data_t) & array, + const size_t & rows, + const size_t & cols) +{ + array.resize(rows); + for (size_t j = 0; j < rows; j++) + { + array[j].resize(cols); + } +} + +//---------------------------------------------------------------- +// assign +// +template +void +assign(array2d(data_t) & array, + const size_t & rows, + const size_t & cols, + const data_t & value) +{ + array.resize(rows); + for (size_t j = 0; j < rows; j++) + { + array[j].assign(cols, value); + } +} + +//---------------------------------------------------------------- +// resize +// +template +void +resize(array3d(data_t) & array, + const size_t & slices, + const size_t & rows, + const size_t & cols) +{ + array.resize(slices); + for (size_t i = 0; i < slices; i++) + { + array[i].resize(rows); + for (size_t j = 0; j < rows; j++) + { + array[i][j].resize(cols); + } + } +} + + +//---------------------------------------------------------------- +// push_back_unique +// +template +void +push_back_unique(container_t & container, + const data_t & data) +{ + typename container_t::const_iterator where = + std::find(container.begin(), container.end(), data); + if (where != container.end()) return; + + container.push_back(data); +} + +//---------------------------------------------------------------- +// push_front_unique +// +template +void +push_front_unique(container_t & container, + const data_t & data) +{ + typename container_t::const_iterator where = + std::find(container.begin(), container.end(), data); + if (where != container.end()) return; + + container.push_front(data); +} + +//---------------------------------------------------------------- +// remove_head +// +template +data_t +remove_head(std::list & container) +{ + data_t head = container.front(); + container.pop_front(); + return head; +} + +//---------------------------------------------------------------- +// remove_tail +// +template +data_t +remove_tail(std::list & container) +{ + data_t tail = container.back(); + container.pop_back(); + return tail; +} + +//---------------------------------------------------------------- +// remove_head +// +template +data_t +remove_head(std::vector & container) +{ + data_t head = container.front(); + container.erase(container.begin()); + return head; +} + +//---------------------------------------------------------------- +// remove_tail +// +template +data_t +remove_tail(std::vector & container) +{ + data_t tail = container.back(); + container.pop_back(); + return tail; +} + +//---------------------------------------------------------------- +// is_size_two_or_larger +// +template +inline bool +is_size_two_or_larger(const container_t & c) +{ + typename container_t::const_iterator i = c.begin(); + typename container_t::const_iterator e = c.end(); + return (i != e) && (++i != e); +} + +//---------------------------------------------------------------- +// is_size_three_or_larger +// +template +inline bool +is_size_three_or_larger(const container_t & c) +{ + typename container_t::const_iterator i = c.begin(); + typename container_t::const_iterator e = c.end(); + return (i != e) && (++i != e) && (++i != e); +} + +//---------------------------------------------------------------- +// is_size_one +// +template +inline bool +is_size_one(const container_t & c) +{ + typename container_t::const_iterator i = c.begin(); + typename container_t::const_iterator e = c.end(); + return (i != e) && (++i == e); +} + + +//---------------------------------------------------------------- +// replace +// +template +bool +replace(container_t & container, const data_t & a, const data_t & b) +{ + typename container_t::iterator end = container.end(); + typename container_t::iterator i = std::find(container.begin(), end, a); + if (i == end) return false; + + *i = b; + return true; +} + + +//---------------------------------------------------------------- +// iterator_at_index +// +template +typename std::list::const_iterator +iterator_at_index(const std::list & container, + const size_t & index) +{ + typename std::list::const_iterator iter = container.begin(); + for (size_t i = 0; i < index && iter != container.end(); i++, ++iter) ; + return iter; +} + +//---------------------------------------------------------------- +// iterator_at_index +// +template +typename std::list::iterator +iterator_at_index(std::list & container, + const size_t & index) +{ + typename std::list::iterator iter = container.begin(); + for (size_t i = 0; i < index && iter != container.end(); i++, ++iter) ; + return iter; +} + +//---------------------------------------------------------------- +// index_of +// +template +size_t +index_of(const std::list & container, const data_t & data) +{ + typename std::list::const_iterator iter = container.begin(); + for (size_t i = 0; iter != container.end(); i++, ++iter) + { + if (data == *iter) return i; + } + + return ~0; +} + +//---------------------------------------------------------------- +// has +// +template +bool +has(const std::list & container, const data_t & data) +{ + typename std::list::const_iterator iter = + std::find(container.begin(), container.end(), data); + + return iter != container.end(); +} + +//---------------------------------------------------------------- +// expand +// +template +container_t & +expand(container_t & a, + const container_t & b, + const bool unique = false) +{ + for (typename container_t::const_iterator i = b.begin(); i != b.end(); ++i) + { + if (unique) + { + push_back_unique(a, *i); + } + else + { + a.push_back(*i); + } + } + + return a; +} + +//---------------------------------------------------------------- +// next +// +template +iterator_t +next(const iterator_t & curr) +{ + iterator_t tmp(curr); + return ++tmp; +} + +//---------------------------------------------------------------- +// prev +// +template +iterator_t +prev(const iterator_t & curr) +{ + iterator_t tmp(curr); + return --tmp; +} + +//---------------------------------------------------------------- +// dump +// +template +stream_t & +dump(stream_t & so, const std::list & c) +{ + for (typename std::list::const_iterator + i = c.begin(); i != c.end(); ++i) + { + so << *i << ' '; + } + + return so; +} + + +//---------------------------------------------------------------- +// operator << +// +template +std::ostream & +operator << (stream_t & so, const std::list & c) +{ + return dump(so, c); +} + + +//---------------------------------------------------------------- +// operator + +// +// Construct an on-the-fly linked list containing two elements: +// +template +inline std::list +operator + (const T & a, const T & b) +{ + std::list ab; + ab.push_back(a); + ab.push_back(b); + return ab; +} + +//---------------------------------------------------------------- +// operator + +// +// Construct an on-the-fly linked list containing list a with item b appended: +template +inline std::list +operator + (const std::list & a, const T & b) +{ + std::list ab(a); + ab.push_back(b); + return ab; +} + +//---------------------------------------------------------------- +// inserter_t +// +template +class inserter_t +{ +public: + inserter_t(container_t & container, + const typename container_t::iterator & iter, + const bool & expand): + container_(container), + iter_(iter), + expand_(expand) + {} + + inline inserter_t & operator << (const data_t & data) + { + if (iter_ == container_.end()) + { + if (expand_) + { + iter_ = container_.insert(iter_, data); + } + else + { + assert(0); + } + } + else + { + *iter_ = data; + ++iter_; + } + + return *this; + } + +private: + // reference to the container: + container_t & container_; + + // current index into the container: + typename container_t::iterator iter_; + + // a flag indicating whether the container should be expanded to + // accomodate the insertions: + bool expand_; +}; + +//---------------------------------------------------------------- +// operator << +// +template +inserter_t, data_t> +operator << (std::vector & container, const data_t & data) +{ + inserter_t, data_t> + inserter(container, container.begin(), false); + return inserter << data; +} + +//---------------------------------------------------------------- +// operator << +// +template +inserter_t, data_t> +operator << (std::list & container, const data_t & data) +{ + inserter_t, data_t> + inserter(container, container.begin(), true); + return inserter << data; +} + + +//---------------------------------------------------------------- +// calc_euclidian_distance_sqrd +// +template +data_t +calc_euclidian_distance_sqrd(const std::vector & a, + const std::vector & b) +{ + data_t distance_sqrd = data_t(0); + for (size_t i = 0; i < dimensions; i++) + { + data_t d = a[i] - b[i]; + distance_sqrd += d * d; + } + + return distance_sqrd; +} + +//---------------------------------------------------------------- +// calc_euclidian_distance +// +template +data_t +calc_euclidian_distance(const std::vector & a, + const std::vector & b) +{ + data_t dist_sqrd = calc_euclidian_distance_sqrd(a, b); + double dist = sqrt(double(dist_sqrd)); + return data_t(dist); +} + +//---------------------------------------------------------------- +// calc_frobenius_norm_sqrd +// +template +data_t +calc_frobenius_norm_sqrd(const std::vector & vec) +{ + data_t L2_norm_sqrd = data_t(0); + + const size_t len = vec.size(); + for (size_t i = 0; i < len; i++) + { + L2_norm_sqrd += vec[i] * vec[i]; + } + + return L2_norm_sqrd; +} + +//---------------------------------------------------------------- +// calc_frobenius_norm +// +template +data_t +calc_frobenius_norm(const std::vector & vec) +{ + data_t norm_sqrd = calc_frobenius_norm_sqrd(vec); + double norm = sqrt(double(norm_sqrd)); + return data_t(norm); +} + +//---------------------------------------------------------------- +// normalize +// +template +void +normalize(std::vector & vec) +{ + data_t norm = calc_frobenius_norm(vec); + + const size_t len = vec.size(); + for (size_t i = 0; i < len; i++) + { + vec[i] /= norm; + } +} + + +//---------------------------------------------------------------- +// the_sign +// +template +inline T +the_sign(const T & a) +{ + if (a < 0) return -1; + if (a > 0) return 1; + return 0; +} + + +//---------------------------------------------------------------- +// copy_a_to_b +// +// list -> array: +// +template +void +copy_a_to_b(const std::list & container_a, + std::vector & container_b) +{ + container_b.assign(container_a.begin(), container_a.end()); +} + +//---------------------------------------------------------------- +// copy_a_to_b +// +// list -> array: +// +template +void +copy_a_to_b(const std::list & container_a, + the_dynamic_array_t & container_b) +{ + container_b.resize(container_a.size()); + + const size_t size = container_a.size(); + typename std::list::const_iterator iter = container_a.begin(); + for (size_t i = 0; i < size; i++, ++iter) + { + container_b[i] = *iter; + } +} + +//---------------------------------------------------------------- +// copy_a_to_b +// +// dynamic_array -> array +// +template +void +copy_a_to_b(const the_dynamic_array_t & container_a, + std::vector & container_b) +{ + container_b.resize(container_a.size()); + + const size_t & size = container_a.size(); + for (size_t i = 0; i < size; i++) + { + container_b[i] = container_a[i]; + } +} + + +//---------------------------------------------------------------- +// the_lock_t +// +template +class the_lock_t +{ +public: + the_lock_t(T * lock, bool lock_immediately = true): + lock_(lock), + armed_(false) + { if (lock_immediately) arm(); } + + the_lock_t(T & lock, bool lock_immediately = true): + lock_(&lock), + armed_(false) + { if (lock_immediately) arm(); } + + ~the_lock_t() + { disarm(); } + + inline void arm() + { + if (!armed_ && lock_ != NULL) + { + lock_->lock(); + armed_ = true; + } + } + + inline void disarm() + { + if (armed_ && lock_ != NULL) + { + lock_->unlock(); + armed_ = false; + } + } + +private: + the_lock_t(); + the_lock_t(const the_lock_t &); + the_lock_t & operator = (const the_lock_t &); + + T * lock_; + bool armed_; +}; + +//---------------------------------------------------------------- +// the_unlock_t +// +template +class the_unlock_t +{ +public: + the_unlock_t(T * lock): + lock_(lock) + { + if (lock_ != NULL) + { + assert(lock_->try_lock() == false); + } + } + + the_unlock_t(T & lock): + lock_(&lock) + { + if (lock_ != NULL) + { + assert(lock_->try_lock() == false); + } + } + + ~the_unlock_t() + { + if (lock_ != NULL) + { + lock_->unlock(); + } + } + +private: + the_unlock_t(); + the_unlock_t(const the_unlock_t &); + the_unlock_t & operator = (const the_unlock_t &); + + T * lock_; +}; + +//---------------------------------------------------------------- +// the_scoped_variable_t +// +template +class the_scoped_variable_t +{ +public: + the_scoped_variable_t(T & variable, + const T & in_scope_value, + const T & out_of_scope_value): + var_(variable), + in_(in_scope_value), + out_(out_of_scope_value) + { + var_ = in_; + } + + ~the_scoped_variable_t() + { + var_ = out_; + } + +private: + the_scoped_variable_t(const the_scoped_variable_t &); + the_scoped_variable_t & operator = (const the_scoped_variable_t &); + + T & var_; + const T in_; + const T out_; +}; + +//---------------------------------------------------------------- +// the_scoped_increment_t +// +template +class the_scoped_increment_t +{ +public: + the_scoped_increment_t(T & variable): + var_(variable) + { + var_++; + } + + ~the_scoped_increment_t() + { + var_--; + } + +private: + the_scoped_increment_t(const the_scoped_increment_t &); + the_scoped_increment_t & operator = (const the_scoped_increment_t &); + + T & var_; +}; + +//---------------------------------------------------------------- +// THROW_ARG2_IF_FALSE +// +#ifndef THROW_ARG2_IF_FALSE +#define THROW_ARG2_IF_FALSE(predicate, arg2) \ +if (predicate) {} else throw arg2 +#endif + +//---------------------------------------------------------------- +// restore_console_stdio +// +// Reopen stdin, stdout, stderr on windows, no-op everywhere else +// +extern bool +restore_console_stdio(); + + +//---------------------------------------------------------------- +// off_t +// +#ifdef _WIN32 +#define off_t __int64 +#endif + +namespace the +{ + extern int open_utf8(const char * filename_utf8, + int oflag, + int pmode); + + extern void open_utf8(std::fstream & fstream_to_open, + const char * filename_utf8, + std::ios_base::openmode mode); + + extern FILE * fopen_utf8(const char * filename_utf8, + const char * mode); + + extern int rename_utf8(const char * old_utf8, const char * new_utf8); + extern int remove_utf8(const char * filename_utf8); + + extern int rmdir_utf8(const char * path_utf8); + extern int mkdir_utf8(const char * path_utf8); + + extern int fseek64(FILE * file, off_t offset, int whence); + extern off_t ftell64(const FILE * file); + + inline static bool + close_enough(const float & ref, + const float & given, + const float tolerance = 1e-6f) + { + float err = fabsf(given - ref); + return err < tolerance; + } + + inline static bool + close_enough(const double & ref, + const double & given, + const double tolerance = 1e-6) + { + double err = fabs(given - ref); + return err < tolerance; + } +} + + +#endif // THE_UTILS_HXX_ diff --git a/include/itkInverseTransform.h b/include/itkInverseTransform.h new file mode 100644 index 0000000..276b815 --- /dev/null +++ b/include/itkInverseTransform.h @@ -0,0 +1,134 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkInverseTransform.h +// Author : Pavel A. Koshevoy +// Created : 2005/06/03 10:16 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A inverse transform class. + +#ifndef __itkInverseTransform_h +#define __itkInverseTransform_h + +// system includes: +#include + +// ITK includes: +#include +#include + + +//---------------------------------------------------------------- +// itk::InverseTransform +// +namespace itk +{ + //---------------------------------------------------------------- + // InverseTransform + // + template + class InverseTransform : + public Transform< typename ForwardTransform::ScalarType, + ForwardTransform::OutputSpaceDimension, + ForwardTransform::InputSpaceDimension > + { + public: + /** Standard class typedefs. */ + typedef InverseTransform Self; + + typedef Transform< typename ForwardTransform::ScalarType, + ForwardTransform::OutputSpaceDimension, + ForwardTransform::InputSpaceDimension > Superclass; + + typedef SmartPointer< Self > Pointer; + typedef SmartPointer< const Self > ConstPointer; + + /** Base inverse transform type. */ + typedef typename Superclass::InverseTransformType InverseTransformType; + typedef SmartPointer< InverseTransformType > InverseTransformPointer; + + /** New method for creating an object using a factory. */ + itkNewMacro(Self); + + /** Run-time type information (and related methods). */ + itkTypeMacro( InverseTransform, Transform ); + + /** Standard scalar type for this class. */ + typedef typename Superclass::ScalarType ScalarType; + + /** Type of the input parameters. */ + typedef typename Superclass::ParametersType ParametersType; + + /** Type of the Jacobian matrix. */ + typedef typename Superclass::JacobianType JacobianType; + + /** Standard coordinate point type for this class. */ + typedef typename Superclass::InputPointType InputPointType; + typedef typename Superclass::OutputPointType OutputPointType; + + /** Dimension of the domain space. */ + itkStaticConstMacro(InputSpaceDimension, + unsigned int, + ForwardTransform::OutputSpaceDimension); + itkStaticConstMacro(OutputSpaceDimension, + unsigned int, + ForwardTransform::InputSpaceDimension); + + /** Set the forward transform pointer. */ + void SetForwardTransform(const ForwardTransform * forward) + { forward_ = forward; } + + /** Method to transform a point. */ + virtual OutputPointType TransformPoint(const InputPointType & y) const + { + assert(forward_ != NULL); + return forward_->BackTransformPoint(y); + } + + virtual const JacobianType & GetJacobian(const InputPointType &) const + { + itkExceptionMacro(<< "GetJacobian is not implemented " + "for InverseTransform"); + return this->m_Jacobian; + }; + + virtual unsigned int GetNumberOfParameters() const + { return 0; } + + virtual InverseTransformPointer GetInverse() const + { return const_cast(forward_); } + + protected: + InverseTransform(): Superclass(0, 0) {} + + private: + // disable default copy constructor and assignment operator: + InverseTransform(const Self & other); + const Self & operator = (const Self & t); + + // the transform whose inverse we are trying to evaluate: + const ForwardTransform * forward_; + }; + +} // namespace itk + +#endif // __itkInverseTransform_h diff --git a/include/itkLegendrePolynomialTransform.h b/include/itkLegendrePolynomialTransform.h new file mode 100644 index 0000000..5dc98e0 --- /dev/null +++ b/include/itkLegendrePolynomialTransform.h @@ -0,0 +1,410 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkLegendrePolynomialTransform.h +// Author : Pavel A. Koshevoy +// Created : 2006/02/22 15:55 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A bivariate centered/normalized Legendre polynomial +// transform and helper functions. + +#ifndef __itkLegendrePolynomialTransform_h +#define __itkLegendrePolynomialTransform_h + +// system includes: +#include +#include + +// ITK includes: +#include +#include +#include + +// local includes: +#include "itkInverseTransform.h" + +// VXL includes: +#include + + +//---------------------------------------------------------------- +// itk::LegendrePolynomialTransform +// +// Let +// A = (u - uc) / Xmax +// B = (v - vc) / Ymax +// +// where uc, vc correspond to the center of the image expressed in +// the coordinate system of the mosaic. +// +// The transform is defined as +// x(u, v) = Xmax * Sa +// y(u, v) = Ymax * Sb +// +// where +// Sa = sum(i in [0, N], sum(j in [0, i], a_jk * Pj(A) * Qk(B))); +// Sb = sum(i in [0, N], sum(j in [0, i], b_jk * Pj(A) * Qk(B))); +// +// where k = i - j and (Pj, Qk) are Legendre polynomials +// of degree (j, k) respectively. +// +namespace itk +{ + template < class TScalar = double, unsigned int N = 2 > + class LegendrePolynomialTransform : + public Transform + { + public: + // standard typedefs: + typedef LegendrePolynomialTransform Self; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + typedef Transform Superclass; + + // Base inverse transform type: + typedef typename Superclass::InverseTransformType InverseTransformType; + typedef SmartPointer< InverseTransformType > InverseTransformPointer; + + // static constant for the degree of the polynomial: + itkStaticConstMacro(Degree, unsigned int, N); + + // static constant for the number of a_jk (or b_jk) coefficients: + itkStaticConstMacro(CoefficientsPerDimension, + unsigned int, + ((N + 1) * (N + 2)) / 2); + + // static constant for the length of the parameter vector: + itkStaticConstMacro(ParameterVectorLength, + unsigned int, + (N + 1) * (N + 2)); + + // RTTI: + itkTypeMacro(LegendrePolynomialTransform, Transform); + + // macro for instantiation through the object factory: + itkNewMacro(Self); + + /** Standard scalar type for this class. */ + typedef typename Superclass::ScalarType ScalarType; + + /** Dimension of the domain space. */ + itkStaticConstMacro(InputSpaceDimension, unsigned int, 2); + itkStaticConstMacro(OutputSpaceDimension, unsigned int, 2); + + // shortcuts: + typedef typename Superclass::ParametersType ParametersType; + typedef typename Superclass::JacobianType JacobianType; + + typedef typename Superclass::InputPointType InputPointType; + typedef typename Superclass::OutputPointType OutputPointType; + + // virtual: + OutputPointType TransformPoint(const InputPointType & x) const; + + // Inverse transformations: + // If y = Transform(x), then x = BackTransform(y); + // if no mapping from y to x exists, then an exception is thrown. + InputPointType BackTransformPoint(const OutputPointType & y) const; + + // virtual: + void SetFixedParameters(const ParametersType & params) + { this->m_FixedParameters = params; } + + // virtual: + const ParametersType & GetFixedParameters() const + { return this->m_FixedParameters; } + + // virtual: + void SetParameters(const ParametersType & params) + { this->m_Parameters = params; } + + // virtual: + const ParametersType & GetParameters() const + { return this->m_Parameters; } + + // virtual: + unsigned int GetNumberOfParameters() const + { return ParameterVectorLength; } + + // virtual: + const JacobianType & GetJacobian(const InputPointType & point) const; + + // virtual: return an inverse of this transform. + InverseTransformPointer GetInverse() const + { + typedef InverseTransform InvTransformType; + typename InvTransformType::Pointer inv = InvTransformType::New(); + inv->SetForwardTransform(this); + return inv.GetPointer(); + } + + // setup the fixed transform parameters: + void setup(// image bounding box expressed in the physical space: + const double x_min, + const double x_max, + const double y_min, + const double y_max, + + // normalization parameters: + const double Xmax = 0.0, + const double Ymax = 0.0) + { + double & uc_ = this->m_FixedParameters[0]; + double & vc_ = this->m_FixedParameters[1]; + double & xmax_ = this->m_FixedParameters[2]; + double & ymax_ = this->m_FixedParameters[3]; + double & a00_ = this->m_Parameters[index_a(0, 0)]; + double & b00_ = this->m_Parameters[index_b(0, 0)]; + + // center of the image: + double xc = (x_min + x_max) / 2.0; + double yc = (y_min + y_max) / 2.0; + uc_ = xc; + vc_ = yc; + + // setup the normalization parameters: + if (Xmax != 0.0 && Ymax != 0.0) + { + xmax_ = Xmax; + ymax_ = Ymax; + } + else + { + const double w = x_max - x_min; + const double h = y_max - y_min; + + // -1 : 1 + xmax_ = w / 2.0; + ymax_ = h / 2.0; + } + + // setup a00, b00 (local translation parameters): + a00_ = xc / xmax_; + b00_ = yc / ymax_; + } + + // setup the translation parameters: + void setup_translation(// translation is expressed in the physical space: + const double tx_Xmax = 0.0, + const double ty_Ymax = 0.0) + { + // incorporate translation into the (uc, vc) fixed parameters: + double & uc_ = this->m_FixedParameters[0]; + double & vc_ = this->m_FixedParameters[1]; + + // FIXME: the signs might be wrong here (20051101): + uc_ -= tx_Xmax; + vc_ -= ty_Ymax; + + // FIXME: untested: + /* + double & a00_ = this->m_Parameters[index_a(0, 0)]; + double & b00_ = this->m_Parameters[index_b(0, 0)]; + a00_ -= tx_Xmax / GetXmax(); + b00_ -= ty_Ymax / GetYmax(); + */ + } + + // helper required for numeric inverse transform calculation; + // evaluate F = T(x), J = dT/dx (another Jacobian): + void eval(const std::vector & x, + std::vector & F, + std::vector > & J) const; + + // setup a linear system to solve for the parameters of this + // transform such that it maps points uv to xy: + void setup_linear_system(const unsigned int start_with_degree, + const unsigned int degrees_covered, + const std::vector & uv, + const std::vector & xy, + vnl_matrix & M, + vnl_vector & bx, + vnl_vector & by) const; + + // find the polynomial coefficients such that this transform + // would map uv to xy: + void solve_for_parameters(const unsigned int start_with_degree, + const unsigned int degrees_covered, + const std::vector & uv, // mosaic + const std::vector & xy,// tile + ParametersType & params) const; + + inline void solve_for_parameters(const unsigned int start_with_degree, + const unsigned int degrees_covered, + const std::vector & uv, + const std::vector & xy) + { + ParametersType params = GetParameters(); + solve_for_parameters(start_with_degree, + degrees_covered, + uv, + xy, + params); + SetParameters(params); + } + + // calculate the number of coefficients of a given + // degree range (per dimension): + inline static unsigned int + count_coefficients(const unsigned int start_with_degree, + const unsigned int degrees_covered) + { + return + index_a(0, start_with_degree + degrees_covered) - + index_a(0, start_with_degree); + } + + // accessors to the warp origin expressed in the mosaic coordinate system: + inline const double & GetUc() const + { return this->m_FixedParameters[0]; } + + inline const double & GetVc() const + { return this->m_FixedParameters[1]; } + + // accessors to the normalization parameters Xmax, Ymax: + inline const double & GetXmax() const + { return this->m_FixedParameters[2]; } + + inline const double & GetYmax() const + { return this->m_FixedParameters[3]; } + + // generate a mask of shared parameters: + static void setup_shared_params_mask(bool shared, std::vector & mask) + { + mask.assign(ParameterVectorLength, shared); + mask[index_a(0, 0)] = false; + mask[index_b(0, 0)] = false; + } + + // Convert the j, k indeces associated with a(j, k) coefficient + // into an index that can be used with the parameters array: + inline static unsigned int index_a(const unsigned int & j, + const unsigned int & k) + { return j + ((j + k) * (j + k + 1)) / 2; } + + inline static unsigned int index_b(const unsigned int & j, + const unsigned int & k) + { return CoefficientsPerDimension + index_a(j, k); } + + // virtual: Generate a platform independent name: + std::string GetTransformTypeAsString() const + { + std::string base = Superclass::GetTransformTypeAsString(); + std::ostringstream name; + name << base << '_' << N; + return name.str(); + } + + protected: + LegendrePolynomialTransform(); + + // virtual: + void PrintSelf(std::ostream & s, Indent indent) const; + + private: + // disable default copy constructor and assignment operator: + LegendrePolynomialTransform(const Self & other); + const Self & operator = (const Self & t); + + // polynomial coefficient accessors: + inline const double & a(const unsigned int & j, + const unsigned int & k) const + { return this->m_Parameters[index_a(j, k)]; } + + inline const double & b(const unsigned int & j, + const unsigned int & k) const + { return this->m_Parameters[index_b(j, k)]; } + + }; // class LegendrePolynomialTransform + +} // namespace itk + + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkLegendrePolynomialTransform.hxx" +#endif + + +//---------------------------------------------------------------- +// setup_transform +// +template +typename transform_t::Pointer +setup_transform(const itk::Point & bbox_min, + const itk::Point & bbox_max) +{ + double w = bbox_max[0] - bbox_min[0]; + double h = bbox_max[1] - bbox_min[1]; + double Umax = w / 2.0; + double Vmax = h / 2.0; + + typename transform_t::Pointer t = transform_t::New(); + t->setup(bbox_min[0], + bbox_max[0], + bbox_min[1], + bbox_max[1], + Umax, + Vmax); + + return t; +} + + +//---------------------------------------------------------------- +// setup_transform +// +template +typename transform_t::Pointer +setup_transform(const image_t * image) +{ + typedef typename image_t::IndexType index_t; + + index_t i00; + i00[0] = 0; + i00[1] = 0; + + typename image_t::PointType origin; + image->TransformIndexToPhysicalPoint(i00, origin); + + index_t i11; + i11[0] = 1; + i11[1] = 1; + + typename image_t::PointType spacing; + image->TransformIndexToPhysicalPoint(i11, spacing); + spacing[0] -= origin[0]; + spacing[1] -= origin[1]; + + typename image_t::SizeType sz = + image->GetLargestPossibleRegion().GetSize(); + + typename image_t::PointType bbox_min = origin; + typename image_t::PointType bbox_max; + bbox_max[0] = bbox_min[0] + spacing[0] * double(sz[0]); + bbox_max[1] = bbox_min[1] + spacing[1] * double(sz[1]); + + return setup_transform(bbox_min, bbox_max); +} + + +#endif // __itkLegendrePolynomialTransform_h diff --git a/include/itkLegendrePolynomialTransform.hxx b/include/itkLegendrePolynomialTransform.hxx new file mode 100644 index 0000000..8fe8946 --- /dev/null +++ b/include/itkLegendrePolynomialTransform.hxx @@ -0,0 +1,556 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkLegendrePolynomialTransform.txx +// Author : Pavel A. Koshevoy +// Created : 2006/02/22 15:55 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A bivariate centered/normalized Legendre polynomial +// transform and helper functions. + +#ifndef _itkLegendrePolynomialTransform_txx +#define _itkLegendrePolynomialTransform_txx + +// local includes: +#include "itkNumericInverse.h" + +// system includes: +#include +#include + +// ITK includes: +#include + + +namespace itk +{ + namespace Legendre + { + // Legendre polynomials: + static double P0(const double & /* x */) + { return 1.0; } + + static double P1(const double & x) + { return x; } + + static double P2(const double & x) + { + const double xx = x * x; + return (3.0 * xx - 1.0) / 2.0; + } + + static double P3(const double & x) + { + const double xx = x * x; + return ((5.0 * xx - 3.0) * x) / 2.0; + } + + static double P4(const double & x) + { + const double xx = x * x; + return ((35.0 * xx - 30.0) * xx + 3.0) / 8.0; + } + + static double P5(const double & x) + { + const double xx = x * x; + return (((63.0 * xx - 70.0) * xx + 15.0) * x) / 8.0; + } + + static double P6(const double & x) + { + const double xx = x * x; + return (((231.0 * xx - 315.0) * xx + 105.0) * xx - 5.0) / 16.0; + } + + //---------------------------------------------------------------- + // polynomial_t + // + typedef double(*polynomial_t)(const double & x); + + //---------------------------------------------------------------- + // A partial table of the Legendre polynomials + // + static polynomial_t P[] = { P0, P1, P2, P3, P4, P5, P6 }; + + // first derivatives of the Legendre polynomials: + static double dP0(const double & /* x */) + { return 0.0; } + + static double dP1(const double & /* x */) + { return 1.0; } + + static double dP2(const double & x) + { return 3.0 * x; } + + static double dP3(const double & x) + { + const double xx = x * x; + return (15.0 * xx - 3.0) / 2.0; + } + + static double dP4(const double & x) + { + const double xx = x * x; + return ((35.0 * xx - 15.0) * x) / 2.0; + } + + static double dP5(const double & x) + { + const double xx = x * x; + return ((315.0 * xx - 210.0) * xx + 15.0) / 8.0; + } + + static double dP6(const double & x) + { + const double xx = x * x; + return (((693.0 * xx - 630.0) * xx + 105.0) * x) / 8.0; + } + + //---------------------------------------------------------------- + // A partial table of the derivatives of Legendre polynomials + // + static polynomial_t dP[] = { dP0, dP1, dP2, dP3, dP4, dP5, dP6 }; + } + + //---------------------------------------------------------------- + // LegendrePolynomialTransform + // + template + LegendrePolynomialTransform:: + LegendrePolynomialTransform(): + Superclass(2, ParameterVectorLength) + { + // by default, the parameters are initialized for an identity transform: + // initialize the parameters for an identity transform: + for (unsigned int i = 0; i < ParameterVectorLength; i++) + { + this->m_Parameters[i] = 0.0; + } + + this->m_Parameters[index_a(1, 0)] = 1.0; + this->m_Parameters[index_b(0, 1)] = 1.0; + + // allocate some space for the fixed parameters: + this->m_FixedParameters.SetSize(4); + setup(-1, 1, -1, 1); + + this->m_Parameters[index_a(0, 0)] = GetUc() / GetXmax(); + this->m_Parameters[index_b(0, 0)] = GetVc() / GetYmax(); + + // zero-out the Jacobian: + for (unsigned int i = 0; i < ParameterVectorLength; i++) + { + this->m_Jacobian(0, i) = 0.0; + this->m_Jacobian(1, i) = 0.0; + } + } + + //---------------------------------------------------------------- + // TransformPoint + // + template + typename LegendrePolynomialTransform::OutputPointType + LegendrePolynomialTransform:: + TransformPoint(const InputPointType & x) const + { + const double & uc = this->GetUc(); + const double & vc = this->GetVc(); + const double & Xmax = this->GetXmax(); + const double & Ymax = this->GetYmax(); + + const ScalarType & u = x[0]; + const ScalarType & v = x[1]; + + const double A = (u - uc) / Xmax; + const double B = (v - vc) / Ymax; + + double P[N + 1]; + double Q[N + 1]; + for (unsigned int i = 0; i <= N; i++) + { + P[i] = Legendre::P[i](A); + Q[i] = Legendre::P[i](B); + } + + double Sa = 0.0; + double Sb = 0.0; + + for (unsigned int i = 0; i <= N; i++) + { + for (unsigned int j = 0; j <= i; j++) + { + unsigned int k = i - j; + double PjQk = P[j] * Q[k]; + Sa += a(j, k) * PjQk; + Sb += b(j, k) * PjQk; + } + } + + OutputPointType y; + y[0] = Xmax * Sa; + y[1] = Ymax * Sb; + + return y; + } + + //---------------------------------------------------------------- + // BackTransformPoint + // + template + typename LegendrePolynomialTransform::InputPointType + LegendrePolynomialTransform:: + BackTransformPoint(const OutputPointType & y) const + { + NumericInverse > inverse(*this); + + std::vector vy(2); + std::vector vx(2); + vy[0] = y[0]; + vy[1] = y[1]; + + // initialize x: first guess - x is close to y: + vx = vy; + const double & uc = this->GetUc(); + const double & vc = this->GetVc(); + + vx[0] += uc; + vx[1] += vc; + bool ok = inverse.transform(vy, vx, true); + if (!ok) + { + itk::ExceptionObject e(__FILE__, __LINE__); + e.SetDescription("could not perform a numeric inverse transformation " + "for a Legendre polynomial transform"); + throw e; + } + + OutputPointType x; + x[0] = vx[0]; + x[1] = vx[1]; + + return x; + } + + //---------------------------------------------------------------- + // GetJacobian + // + template + const typename LegendrePolynomialTransform::JacobianType & + LegendrePolynomialTransform:: + GetJacobian(const InputPointType & x) const + { + const double & uc = this->GetUc(); + const double & vc = this->GetVc(); + const double & Xmax = this->GetXmax(); + const double & Ymax = this->GetYmax(); + + const ScalarType & u = x[0]; + const ScalarType & v = x[1]; + + const double A = (u - uc) / Xmax; + const double B = (v - vc) / Ymax; + + double P[N + 1]; + double Q[N + 1]; + for (unsigned int i = 0; i <= N; i++) + { + P[i] = Legendre::P[i](A); + Q[i] = Legendre::P[i](B); + } + + // derivatives with respect to a_jk, b_jk: + for (unsigned int i = 0; i <= N; i++) + { + for (unsigned int j = 0; j <= i; j++) + { + unsigned int k = i - j; + + double PjQk = P[j] * Q[k]; + this->m_Jacobian(0, index_a(j, k)) = Xmax * PjQk; + this->m_Jacobian(1, index_b(j, k)) = Ymax * PjQk; + } + } + + // done: + return this->m_Jacobian; + } + + //---------------------------------------------------------------- + // eval + // + template + void + LegendrePolynomialTransform:: + eval(const std::vector & x, + std::vector & F, + std::vector > & J) const + { + const double & uc = this->GetUc(); + const double & vc = this->GetVc(); + const double & Xmax = this->GetXmax(); + const double & Ymax = this->GetYmax(); + + const ScalarType & u = x[0]; + const ScalarType & v = x[1]; + + const double A = (u - uc) / Xmax; + const double B = (v - vc) / Ymax; + + double P[N + 1]; + double Q[N + 1]; + for (unsigned int i = 0; i <= N; i++) + { + P[i] = Legendre::P[i](A); + Q[i] = Legendre::P[i](B); + } + + double Sa = 0.0; + double Sb = 0.0; + + // derivatives with respect to a_jk, b_jk: + for (unsigned int i = 0; i <= N; i++) + { + for (unsigned int j = 0; j <= i; j++) + { + unsigned int k = i - j; + double PjQk = P[j] * Q[k]; + Sa += a(j, k) * PjQk; + Sb += b(j, k) * PjQk; + } + } + + F[0] = Xmax * Sa; + F[1] = Ymax * Sb; + + // derivatives with respect to u: + double dSa_du = 0.0; + double dSb_du = 0.0; + + for (unsigned int i = 1; i <= N; i++) + { + for (unsigned int j = 0; j <= i; j++) + { + unsigned int k = i - j; + double dPjQk = Legendre::dP[j](A) * Q[k]; + dSa_du += a(j, k) * dPjQk; + dSb_du += b(j, k) * dPjQk; + } + } + + // derivatives with respect to v: + double dSa_dv = 0.0; + double dSb_dv = 0.0; + + for (unsigned int i = 1; i <= N; i++) + { + for (unsigned int j = 0; j <= i; j++) + { + unsigned int k = i - j; + double dQkPj = Legendre::dP[k](B) * P[j]; + dSa_dv += a(j, k) * dQkPj; + dSb_dv += b(j, k) * dQkPj; + } + } + + // dx/du: + J[0][0] = dSa_du; + + // dx/dv: + J[0][1] = Xmax / Ymax * dSa_dv; + + // dy/du: + J[1][0] = Ymax / Xmax * dSb_du; + + // dy/dv: + J[1][1] = dSb_dv; + } + + //---------------------------------------------------------------- + // setup_linear_system + // + template + void + LegendrePolynomialTransform:: + setup_linear_system(const unsigned int start_with_degree, + const unsigned int degrees_covered, + const std::vector & uv, + const std::vector & xy, + vnl_matrix & M, + vnl_vector & bx, + vnl_vector & by) const + { + if (degrees_covered == 0) return; + + const double & uc = this->GetUc(); + const double & vc = this->GetVc(); + const double & Xmax = this->GetXmax(); + const double & Ymax = this->GetYmax(); + + static double P[N + 1]; + static double Q[N + 1]; + + const unsigned int & num_points = uv.size(); + assert(num_points == xy.size()); + + const unsigned int offset = + index_a(0, start_with_degree); + + const unsigned int extent = + index_a(0, start_with_degree + degrees_covered) - offset; + + M.set_size(num_points, extent); + bx.set_size(num_points); + by.set_size(num_points); + + for (unsigned int row = 0; row < num_points; row++) + { + const ScalarType & u = uv[row][0]; + const ScalarType & v = uv[row][1]; + const ScalarType & x = xy[row][0]; + const ScalarType & y = xy[row][1]; + + const double A = (u - uc) / Xmax; + const double B = (v - vc) / Ymax; + + bx[row] = x / Xmax; + by[row] = y / Ymax; + + for (unsigned int i = 0; i < start_with_degree + degrees_covered; i++) + { + P[i] = Legendre::P[i](A); + Q[i] = Legendre::P[i](B); + } + + for (unsigned int i = 0; i < start_with_degree; i++) + { + for (unsigned int j = 0; j <= i; j++) + { + const unsigned int k = i - j; + double PjQk = P[j] * Q[k]; + bx[row] -= a(j, k) * PjQk; + by[row] -= b(j, k) * PjQk; + } + } + + for (unsigned int i = start_with_degree; + i < start_with_degree + degrees_covered; + i++) + { + for (unsigned int j = 0; j <= i; j++) + { + const unsigned int k = i - j; + const unsigned int col = index_a(j, k) - offset; + + M[row][col] = P[j] * Q[k]; + } + } + } + } + + //---------------------------------------------------------------- + // solve_for_parameters + // + template + void + LegendrePolynomialTransform:: + solve_for_parameters(const unsigned int start_with_degree, + const unsigned int degrees_covered, + const std::vector & uv, + const std::vector & xy, + ParametersType & params) const + { + if (degrees_covered == 0) return; + + vnl_matrix M; + vnl_vector bx; + vnl_vector by; + + setup_linear_system(start_with_degree, degrees_covered, uv, xy, M, bx, by); + + // use SVD to find the inverse of M: + vnl_svd svd(M); + vnl_vector xa = svd.solve(bx); + vnl_vector xb = svd.solve(by); + + unsigned int offset = index_a(0, start_with_degree); + for (unsigned int i = start_with_degree; + i < start_with_degree + degrees_covered; + i++) + { + for (unsigned int j = 0; j <= i; j++) + { + const unsigned int k = i - j; + const unsigned int a_jk = index_a(j, k); + const unsigned int b_jk = index_b(j, k); + + params[a_jk] = xa[a_jk - offset]; + params[b_jk] = xb[a_jk - offset]; + } + } + } + + //---------------------------------------------------------------- + // PrintSelf + // + template + void + LegendrePolynomialTransform:: + PrintSelf(std::ostream & os, Indent indent) const + { + Superclass::PrintSelf(os,indent); + + for (unsigned int i = 0; i <= N; i++) + { + for (unsigned int j = 0; j <= i; j++) + { + unsigned int k = i - j; + os << indent << "a(" << j << ", " << k << ") = " << a(j, k) + << std::endl; + } + } + for (unsigned int i = 0; i <= N; i++) + { + for (unsigned int j = 0; j <= i; j++) + { + unsigned int k = i - j; + os << indent << "b(" << j << ", " << k << ") = " << b(j, k) + << std::endl; + } + } + os << indent << "uc = " << GetUc() << std::endl + << indent << "vc = " << GetVc() << std::endl + << indent << "Xmax = " << GetXmax() << std::endl + << indent << "Ymax = " << GetYmax() << std::endl; +#if 0 + for (unsigned int i = 0; i < ParameterVectorLength; i++) + { + os << indent << "params[" << i << "] = " << this->m_Parameters[i] << ';' + << std::endl; + } +#endif + } + +} // namespace itk + + +#endif // _itkLegendrePolynomialTransform_txx diff --git a/include/itkNormalizeImageFilterWithMask.h b/include/itkNormalizeImageFilterWithMask.h new file mode 100644 index 0000000..c0163e6 --- /dev/null +++ b/include/itkNormalizeImageFilterWithMask.h @@ -0,0 +1,111 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkNormalizeImageFilterWithMask.h +// Author : Pavel A. Koshevoy +// Created : 2006/06/15 17:09 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : An enhanced version of the itk::NormalizeImageFilter +// adding support for a mask image. + +#ifndef __itkNormalizeImageFilterWithMask_h +#define __itkNormalizeImageFilterWithMask_h + +// ITK includes: +#include +#include +#include +#include + +// local includes: +#include "itkStatisticsImageFilterWithMask.h" + +namespace itk { + +/** \class NormalizeImageFilterWithMask + * \brief Normalize an image by setting its mean to zero and variance to one. + * + * NormalizeImageFilterWithMask shifts and scales an image so that the pixels + * in the image have a zero mean and unit variance. This filter uses + * StatisticsImageFilterWithMask to compute the mean and variance of the input + * and then applies ShiftScaleImageFilter to shift and scale the pixels. + * + * NB: since this filter normalizes the data to lie within -1 to 1, + * integral types will produce an image that DOES NOT HAVE a unit variance. + * \ingroup MathematicalImageFilters + */ +template +class ITK_EXPORT NormalizeImageFilterWithMask : + public ImageToImageFilter +{ +public: + /** Standard Self typedef */ + typedef NormalizeImageFilterWithMask Self; + typedef ImageToImageFilter Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(NormalizeImageFilterWithMask, ImageToImageFilter); + + /** Image related typedefs. */ + typedef typename TInputImage::Pointer InputImagePointer; + typedef typename TOutputImage::Pointer OutputImagePointer; + + /** Type for the mask of the fixed image. Only pixels that are "inside" + this mask will be considered for the computation of the metric */ + typedef SpatialObject ImageMaskType; + typedef typename ImageMaskType::Pointer ImageMaskPointer; + + /** Set/Get the image mask. */ + itkSetObjectMacro( ImageMask, ImageMaskType ); + itkGetConstObjectMacro( ImageMask, ImageMaskType ); + +protected: + NormalizeImageFilterWithMask(); + + /** GenerateData. */ + void GenerateData (); + + // Override since the filter needs all the data for the algorithm + void GenerateInputRequestedRegion(); + +private: + NormalizeImageFilterWithMask(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + typename StatisticsImageFilterWithMask::Pointer m_StatisticsFilter; + typename ShiftScaleImageFilter::Pointer m_ShiftScaleFilter; + + mutable ImageMaskPointer m_ImageMask; +} ; // end of class + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkNormalizeImageFilterWithMask.hxx" +#endif + +#endif diff --git a/include/itkNormalizeImageFilterWithMask.hxx b/include/itkNormalizeImageFilterWithMask.hxx new file mode 100644 index 0000000..b2eef6a --- /dev/null +++ b/include/itkNormalizeImageFilterWithMask.hxx @@ -0,0 +1,98 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkNormalizeImageFilterWithMask.txx +// Author : Pavel A. Koshevoy +// Created : 2006/06/15 17:09 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : An enhanced version of the itk::NormalizeImageFilter +// adding support for a mask image. + +#ifndef _itkNormalizeImageFilterWithMask_txx +#define _itkNormalizeImageFilterWithMask_txx + +// ITK includes: +#include +#include +#include +#include + +namespace itk +{ + +template +NormalizeImageFilterWithMask +::NormalizeImageFilterWithMask() +{ + m_StatisticsFilter = 0; + m_StatisticsFilter = StatisticsImageFilterWithMask::New(); + m_ShiftScaleFilter = ShiftScaleImageFilter::New(); +} + +template +void +NormalizeImageFilterWithMask +::GenerateInputRequestedRegion() +{ + Superclass::GenerateInputRequestedRegion(); + if ( this->GetInput() ) + { + InputImagePointer image = + const_cast< typename Superclass::InputImageType * >( this->GetInput() ); + image->SetRequestedRegionToLargestPossibleRegion(); + } +} + +template +void +NormalizeImageFilterWithMask +::GenerateData() +{ + ProgressAccumulator::Pointer progress = ProgressAccumulator::New(); + progress->SetMiniPipelineFilter(this); + + progress->RegisterInternalFilter(m_StatisticsFilter,.5f); + progress->RegisterInternalFilter(m_ShiftScaleFilter,.5f); + + // Gather statistics + + m_StatisticsFilter->SetInput(this->GetInput()); + m_StatisticsFilter->GetOutput()->SetRequestedRegion(this->GetOutput()->GetRequestedRegion()); + m_StatisticsFilter->SetImageMask(this->m_ImageMask); + m_StatisticsFilter->Update(); + + // Set the parameters for Shift + m_ShiftScaleFilter->SetShift(-m_StatisticsFilter->GetMean()); + m_ShiftScaleFilter->SetScale(NumericTraits::RealType>::One + / m_StatisticsFilter->GetSigma()); + m_ShiftScaleFilter->SetInput(this->GetInput()); + + m_ShiftScaleFilter->GetOutput()->SetRequestedRegion(this->GetOutput()->GetRequestedRegion()); + m_ShiftScaleFilter->Update(); + + // Graft the mini pipeline output to this filters output + this->GraftOutput(m_ShiftScaleFilter->GetOutput()); +} + +} // end namespace itk + +#endif diff --git a/include/itkNumericInverse.h b/include/itkNumericInverse.h new file mode 100644 index 0000000..84bd026 --- /dev/null +++ b/include/itkNumericInverse.h @@ -0,0 +1,183 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkNumericTransform.h +// Author : Pavel A. Koshevoy +// Created : 2005/06/03 10:16 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A numerical inverse transform class, based on the +// Newton-Raphson method for nonlinear systems of equations. + +#ifndef __itkNumericInverse_h +#define __itkNumericInverse_h + +// system includes: +#include +#include +#include + +// ITK includes: +#include +#include + +// VXL includes: +#include + +namespace help +{ + //---------------------------------------------------------------- + // nonlinear_system_evaluator_t + // + template + class nonlinear_system_evaluator_t + { + public: + virtual ~nonlinear_system_evaluator_t() {} + + // evaluate the nonlinear system of equations + // and its Jacobian (dF/dx) at a given point: + virtual void + operator() (const std::vector & x, + std::vector & F, + std::vector > & J) const = 0; + }; + + //---------------------------------------------------------------- + // NewtonRaphson + // + template + bool + NewtonRaphson(const nonlinear_system_evaluator_t & usrfun, + std::vector & x, // estimated root point + const unsigned int & ntrial, // number of iterations + const ScalarType & tolx, // convergence tolerance in x + const ScalarType & tolf) // convergence tolerance in F + { + // shortcut: + const unsigned int n = x.size(); + + std::vector F(n); + std::vector > J(n); + for (unsigned int i = 0; i < n; i++) J[i].resize(n); + + vnl_matrix A(n, n); + vnl_vector b(n); + + for (unsigned int k = 0; k < ntrial; k++) + { + // evaluate the function at the current position: + usrfun(x, F, J); + + // check for function convergence: + ScalarType errf = ScalarType(0); + for (unsigned int i = 0; i < n; i++) + { + errf += fabs(F[i]); + } + if (errf <= tolf) break; + + // setup the left hand side of the linear system: + for (unsigned int i = 0; i < n; i++) + { + for (unsigned int j = 0; j < n; j++) + { + A[i][j] = J[i][j]; + } + } + + // setup the right hand side of the linear system: + for (unsigned int i = 0; i < n; i++) + { + b[i] = -F[i]; + } + + vnl_svd svd(A); + vnl_vector dx = svd.solve(b); + + // check for root convergence: + ScalarType errx = ScalarType(0); + for (unsigned int i = 0; i < n; i++) + { + errx += fabs(dx[i]); + x[i] += dx[i]; + } + if (errx <= tolx) break; + } + + return true; + } + +} // namespace help + +namespace itk +{ + //---------------------------------------------------------------- + // NumericInverse + // + template + class NumericInverse : + public help::nonlinear_system_evaluator_t + { + public: + typedef TTransform TransformType; + typedef typename TTransform::ScalarType ScalarType; + + NumericInverse(const TransformType & transform): + transform_(transform) + {} + + // virtual: + void operator() (const std::vector & x, + std::vector & F, + std::vector > & J) const + { + transform_.eval(x, F, J); + + const unsigned int n = x.size(); + for (unsigned int i = 0; i < n; i++) + { + F[i] -= y_[i]; + } + } + + // If y = Transform(x), then x = BackTransform(y). + // given y, find x: + bool transform(const std::vector & y, + std::vector & x, + bool x_is_initialized = false) const + { + y_ = y; + if (!x_is_initialized) x = y; + return NewtonRaphson(*this, x, 50, 1e-12, 1e-12); + } + + private: + // the transform whose inverse we are trying to evaluate: + const TransformType & transform_; + + // the point for which we are tryying to find the inverse mapping: + mutable std::vector y_; + }; + +} // namespace itk + +#endif // __itkNumericInverse_h diff --git a/include/itkStatisticsImageFilterWithMask.h b/include/itkStatisticsImageFilterWithMask.h new file mode 100644 index 0000000..1e111f8 --- /dev/null +++ b/include/itkStatisticsImageFilterWithMask.h @@ -0,0 +1,192 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkStatisticsImageFilterWithMask.h +// Author : Pavel A. Koshevoy +// Created : 2006/06/15 17:09 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : An enhanced version of the itk::StatisticsImageFilter +// adding support for a mask image. + +#ifndef __itkStatisticsImageFilterWithMask_h +#define __itkStatisticsImageFilterWithMask_h + +// ITK includes: +#include +#include +#include +#include +#include + + +namespace itk { + +/** \class StatisticsImageFilterWithMask + * \brief Compute min. max, variance and mean of an Image. + * + * StatisticsImageFilterWithMask computes the minimum, maximum, sum, mean, variance + * sigma of an image. The filter needs all of its input image. It + * behaves as a filter with an input and output. Thus it can be inserted + * in a pipline with other filters and the statistics will only be + * recomputed if a downstream filter changes. + * + * The filter passes its input through unmodified. The filter is + * threaded. It computes statistics in each thread then combines them in + * its AfterThreadedGenerate method. + * + * \ingroup MathematicalStatisticsImageFilters + */ +template +class ITK_EXPORT StatisticsImageFilterWithMask : + public ImageToImageFilter +{ +public: + /** Standard Self typedef */ + typedef StatisticsImageFilterWithMask Self; + typedef ImageToImageFilter Superclass; + typedef SmartPointer Pointer; + typedef SmartPointer ConstPointer; + + /** Method for creation through the object factory. */ + itkNewMacro(Self); + + /** Runtime information support. */ + itkTypeMacro(StatisticsImageFilterWithMask, ImageToImageFilter); + + /** Image related typedefs. */ + typedef typename TInputImage::Pointer InputImagePointer; + + typedef typename TInputImage::RegionType RegionType ; + typedef typename TInputImage::SizeType SizeType ; + typedef typename TInputImage::IndexType IndexType ; + typedef typename TInputImage::PixelType PixelType ; + typedef typename TInputImage::PointType PointType ; + + /** Image related typedefs. */ + itkStaticConstMacro(ImageDimension, unsigned int, + TInputImage::ImageDimension ) ; + + /** Type to use for computations. */ + typedef typename NumericTraits::RealType RealType; + + /** Smart Pointer type to a DataObject. */ + typedef typename DataObject::Pointer DataObjectPointer; + + /** Type of DataObjects used for scalar outputs */ + typedef SimpleDataObjectDecorator RealObjectType; + typedef SimpleDataObjectDecorator PixelObjectType; + + /** Type for the mask of the fixed image. Only pixels that are "inside" + this mask will be considered for the computation of the metric */ + typedef SpatialObject ImageMaskType; + typedef typename ImageMaskType::Pointer ImageMaskPointer; + + /** Set/Get the image mask. */ + itkSetObjectMacro( ImageMask, ImageMaskType ); + itkGetConstObjectMacro( ImageMask, ImageMaskType ); + + /** Return the computed Minimum. */ + PixelType GetMinimum() const + { return this->GetMinimumOutput()->Get(); } + PixelObjectType* GetMinimumOutput(); + const PixelObjectType* GetMinimumOutput() const; + + /** Return the computed Maximum. */ + PixelType GetMaximum() const + { return this->GetMaximumOutput()->Get(); } + PixelObjectType* GetMaximumOutput(); + const PixelObjectType* GetMaximumOutput() const; + + /** Return the computed Mean. */ + RealType GetMean() const + { return this->GetMeanOutput()->Get(); } + RealObjectType* GetMeanOutput(); + const RealObjectType* GetMeanOutput() const; + + /** Return the computed Standard Deviation. */ + RealType GetSigma() const + { return this->GetSigmaOutput()->Get(); } + RealObjectType* GetSigmaOutput(); + const RealObjectType* GetSigmaOutput() const; + + /** Return the computed Variance. */ + RealType GetVariance() const + { return this->GetVarianceOutput()->Get(); } + RealObjectType* GetVarianceOutput(); + const RealObjectType* GetVarianceOutput() const; + + /** Return the compute Sum. */ + RealType GetSum() const + { return this->GetSumOutput()->Get(); } + RealObjectType* GetSumOutput(); + const RealObjectType* GetSumOutput() const; + + /** Make a DataObject of the correct type to be used as the specified + * output. */ + virtual DataObjectPointer MakeOutput(unsigned int idx); + +protected: + StatisticsImageFilterWithMask(); + ~StatisticsImageFilterWithMask(){}; + void PrintSelf(std::ostream& os, Indent indent) const; + + /** Pass the input through unmodified. Do this by Grafting in the AllocateOutputs method. */ + void AllocateOutputs(); + + /** Initialize some accumulators before the threads run. */ + void BeforeThreadedGenerateData (); + + /** Do final mean and variance computation from data accumulated in threads. */ + void AfterThreadedGenerateData (); + + /** Multi-thread version GenerateData. */ + void ThreadedGenerateData (const RegionType& + outputRegionForThread, + int threadId) ; + + // Override since the filter needs all the data for the algorithm + void GenerateInputRequestedRegion(); + + // Override since the filter produces all of its output + void EnlargeOutputRequestedRegion(DataObject *data); + +private: + StatisticsImageFilterWithMask(const Self&); //purposely not implemented + void operator=(const Self&); //purposely not implemented + + mutable ImageMaskPointer m_ImageMask; + + Array m_ThreadSum; + Array m_SumOfSquares; + Array m_Count; + Array m_ThreadMin; + Array m_ThreadMax; + +} ; // end of class + +} // end namespace itk + +#ifndef ITK_MANUAL_INSTANTIATION +#include "itkStatisticsImageFilterWithMask.hxx" +#endif + +#endif diff --git a/include/itkStatisticsImageFilterWithMask.hxx b/include/itkStatisticsImageFilterWithMask.hxx new file mode 100644 index 0000000..284f359 --- /dev/null +++ b/include/itkStatisticsImageFilterWithMask.hxx @@ -0,0 +1,390 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : itkStatisticsImageFilterWithMask.txx +// Author : Pavel A. Koshevoy +// Created : 2006/06/15 17:09 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : An enhanced version of the itk::StatisticsImageFilter +// adding support for a mask image. + +#ifndef _itkStatisticsImageFilterWithMask_txx +#define _itkStatisticsImageFilterWithMask_txx + +// ITK includes: +#include +#include +#include +#include + +namespace itk { + +template +StatisticsImageFilterWithMask +::StatisticsImageFilterWithMask(): m_ThreadSum(1), m_SumOfSquares(1), m_Count(1), m_ThreadMin(1), m_ThreadMax(1) +{ + // first output is a copy of the image, DataObject created by + // superclass + // + // allocate the data objects for the outputs which are + // just decorators around pixel types + for (int i=1; i < 3; ++i) + { + typename PixelObjectType::Pointer output + = static_cast(this->MakeOutput(i).GetPointer()); + this->ProcessObject::SetNthOutput(i, output.GetPointer()); + } + // allocate the data objects for the outputs which are + // just decorators around real types + for (int i=3; i < 7; ++i) + { + typename RealObjectType::Pointer output + = static_cast(this->MakeOutput(i).GetPointer()); + this->ProcessObject::SetNthOutput(i, output.GetPointer()); + } + + this->GetMinimumOutput()->Set( NumericTraits::max() ); + this->GetMaximumOutput()->Set( NumericTraits::NonpositiveMin() ); + this->GetMeanOutput()->Set( NumericTraits::max() ); + this->GetSigmaOutput()->Set( NumericTraits::max() ); + this->GetVarianceOutput()->Set( NumericTraits::max() ); + this->GetSumOutput()->Set( NumericTraits::Zero ); +} + + +template +DataObject::Pointer +StatisticsImageFilterWithMask +::MakeOutput(unsigned int output) +{ + switch (output) + { + case 0: + return static_cast(TInputImage::New().GetPointer()); + break; + case 1: + return static_cast(PixelObjectType::New().GetPointer()); + break; + case 2: + return static_cast(PixelObjectType::New().GetPointer()); + break; + case 3: + case 4: + case 5: + case 6: + return static_cast(RealObjectType::New().GetPointer()); + break; + default: + // might as well make an image + return static_cast(TInputImage::New().GetPointer()); + break; + } +} + + +template +typename StatisticsImageFilterWithMask::PixelObjectType* +StatisticsImageFilterWithMask +::GetMinimumOutput() +{ + return static_cast(this->ProcessObject::GetOutput(1)); +} + +template +const typename StatisticsImageFilterWithMask::PixelObjectType* +StatisticsImageFilterWithMask +::GetMinimumOutput() const +{ + return static_cast(this->ProcessObject::GetOutput(1)); +} + + +template +typename StatisticsImageFilterWithMask::PixelObjectType* +StatisticsImageFilterWithMask +::GetMaximumOutput() +{ + return static_cast(this->ProcessObject::GetOutput(2)); +} + +template +const typename StatisticsImageFilterWithMask::PixelObjectType* +StatisticsImageFilterWithMask +::GetMaximumOutput() const +{ + return static_cast(this->ProcessObject::GetOutput(2)); +} + + +template +typename StatisticsImageFilterWithMask::RealObjectType* +StatisticsImageFilterWithMask +::GetMeanOutput() +{ + return static_cast(this->ProcessObject::GetOutput(3)); +} + +template +const typename StatisticsImageFilterWithMask::RealObjectType* +StatisticsImageFilterWithMask +::GetMeanOutput() const +{ + return static_cast(this->ProcessObject::GetOutput(3)); +} + + +template +typename StatisticsImageFilterWithMask::RealObjectType* +StatisticsImageFilterWithMask +::GetSigmaOutput() +{ + return static_cast(this->ProcessObject::GetOutput(4)); +} + +template +const typename StatisticsImageFilterWithMask::RealObjectType* +StatisticsImageFilterWithMask +::GetSigmaOutput() const +{ + return static_cast(this->ProcessObject::GetOutput(4)); +} + + +template +typename StatisticsImageFilterWithMask::RealObjectType* +StatisticsImageFilterWithMask +::GetVarianceOutput() +{ + return static_cast(this->ProcessObject::GetOutput(5)); +} + +template +const typename StatisticsImageFilterWithMask::RealObjectType* +StatisticsImageFilterWithMask +::GetVarianceOutput() const +{ + return static_cast(this->ProcessObject::GetOutput(5)); +} + + +template +typename StatisticsImageFilterWithMask::RealObjectType* +StatisticsImageFilterWithMask +::GetSumOutput() +{ + return static_cast(this->ProcessObject::GetOutput(6)); +} + +template +const typename StatisticsImageFilterWithMask::RealObjectType* +StatisticsImageFilterWithMask +::GetSumOutput() const +{ + return static_cast(this->ProcessObject::GetOutput(6)); +} + +template +void +StatisticsImageFilterWithMask +::GenerateInputRequestedRegion() +{ + Superclass::GenerateInputRequestedRegion(); + if ( this->GetInput() ) + { + InputImagePointer image = + const_cast< typename Superclass::InputImageType * >( this->GetInput() ); + image->SetRequestedRegionToLargestPossibleRegion(); + } +} + +template +void +StatisticsImageFilterWithMask +::EnlargeOutputRequestedRegion(DataObject *data) +{ + Superclass::EnlargeOutputRequestedRegion(data); + data->SetRequestedRegionToLargestPossibleRegion(); +} + + +template +void +StatisticsImageFilterWithMask +::AllocateOutputs() +{ + // Pass the input through as the output + InputImagePointer image = + const_cast< TInputImage * >( this->GetInput() ); + this->GraftOutput( image ); + + // Nothing that needs to be allocated for the remaining outputs +} + +template +void +StatisticsImageFilterWithMask +::BeforeThreadedGenerateData() +{ + int numberOfThreads = this->GetNumberOfThreads(); + + // Resize the thread temporaries + m_Count.SetSize(numberOfThreads); + m_SumOfSquares.SetSize(numberOfThreads); + m_ThreadSum.SetSize(numberOfThreads); + m_ThreadMin.SetSize(numberOfThreads); + m_ThreadMax.SetSize(numberOfThreads); + + // Initialize the temporaries + m_Count.Fill(NumericTraits::Zero); + m_ThreadSum.Fill(NumericTraits::Zero); + m_SumOfSquares.Fill(NumericTraits::Zero); + m_ThreadMin.Fill(NumericTraits::max()); + m_ThreadMax.Fill(NumericTraits::NonpositiveMin()); +} + +template +void +StatisticsImageFilterWithMask +::AfterThreadedGenerateData() +{ + int i; + long count; + RealType sumOfSquares; + + int numberOfThreads = this->GetNumberOfThreads(); + + PixelType minimum; + PixelType maximum; + RealType mean; + RealType sigma; + RealType variance; + RealType sum; + + sum = sumOfSquares = NumericTraits::Zero; + count = 0; + + // Find the min/max over all threads and accumulate count, sum and + // sum of squares + minimum = NumericTraits::max(); + maximum = NumericTraits::NonpositiveMin(); + for( i = 0; i < numberOfThreads; i++) + { + count += m_Count[i]; + sum += m_ThreadSum[i]; + sumOfSquares += m_SumOfSquares[i]; + + if (m_ThreadMin[i] < minimum) + { + minimum = m_ThreadMin[i]; + } + if (m_ThreadMax[i] > maximum) + { + maximum = m_ThreadMax[i]; + } + } + // compute statistics + mean = sum / static_cast(count); + + // unbiased estimate + variance = (sumOfSquares - (sum*sum / static_cast(count))) + / (static_cast(count) - 1); + sigma = sqrt(variance); + + // Set the outputs + this->GetMinimumOutput()->Set( minimum ); + this->GetMaximumOutput()->Set( maximum ); + this->GetMeanOutput()->Set( mean ); + this->GetSigmaOutput()->Set( sigma ); + this->GetVarianceOutput()->Set( variance ); + this->GetSumOutput()->Set( sum ); +} + +template +void +StatisticsImageFilterWithMask +::ThreadedGenerateData(const RegionType& outputRegionForThread, + int threadId) +{ + RealType realValue; + PixelType value; + ImageRegionConstIteratorWithIndex it (this->GetInput(), outputRegionForThread); + + // support progress methods/callbacks + ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels()); + + const TInputImage * image = this->GetInput(); + + // do the work + while (!it.IsAtEnd()) + { + bool skip = false; + if (m_ImageMask) + { + // test the mask: + PointType point; + image->TransformIndexToPhysicalPoint( it.GetIndex(), point ); + skip = !m_ImageMask->IsInside( point ); + } + + if (!skip) + { + value = it.Get(); + realValue = static_cast( value ); + if (value < m_ThreadMin[threadId]) + { + m_ThreadMin[threadId] = value; + } + if (value > m_ThreadMax[threadId]) + { + m_ThreadMax[threadId] = value; + } + + m_ThreadSum[threadId] += realValue; + m_SumOfSquares[threadId] += (realValue * realValue); + m_Count[threadId]++; + } + + ++it; + progress.CompletedPixel(); + } +} + +template +void +StatisticsImageFilterWithMask +::PrintSelf(std::ostream& os, Indent indent) const +{ + Superclass::PrintSelf(os,indent); + + os << indent << "Minimum: " + << static_cast::PrintType>(this->GetMinimum()) << std::endl; + os << indent << "Maximum: " + << static_cast::PrintType>(this->GetMaximum()) << std::endl; + os << indent << "Sum: " << this->GetSum() << std::endl; + os << indent << "Mean: " << this->GetMean() << std::endl; + os << indent << "Sigma: " << this->GetSigma() << std::endl; + os << indent << "Variance: " << this->GetVariance() << std::endl; +} + + +}// end namespace itk +#endif diff --git a/itk-module.cmake b/itk-module.cmake index 009cc43..f000248 100644 --- a/itk-module.cmake +++ b/itk-module.cmake @@ -15,11 +15,15 @@ itk_module(Nornir DEPENDS ITKCommon ITKStatistics + ITKIOImageBase COMPILE_DEPENDS ITKImageSources + ITKImageIntensity + ITKThresholding + ITKSmoothing + ITKRegistrationCommon TEST_DEPENDS ITKTestKernel - ITKMetaIO DESCRIPTION "${DOCUMENTATION}" EXCLUDE_FROM_DEFAULT diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 92c2241..adde20b 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,5 +1,16 @@ set(Nornir_SRCS - itkMinimalStandardRandomVariateGenerator.cxx + itkIRCommon.cxx + itkIRText.cxx + itkIRUtils.cxx + IRTerminator.cxx + IRThreadStorage.cxx + IRThread.cxx + IRThreadInterface.cxx + IRMutex.cxx + IRMutexInterface.cxx + IRTransaction.cxx + IRLog.cxx + IRThreadPool.cxx ) itk_module_add_library(Nornir ${Nornir_SRCS}) diff --git a/src/IRLog.cxx b/src/IRLog.cxx new file mode 100644 index 0000000..b8c81c5 --- /dev/null +++ b/src/IRLog.cxx @@ -0,0 +1,200 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_log.cxx +// Author : Pavel A. Koshevoy +// Created : Fri Mar 23 11:04:53 MDT 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A text log object -- behaves almost like a std::ostream. + +// local includes: +#include "the_log.hxx" + + +//---------------------------------------------------------------- +// the_log_t::the_log_t +// +the_log_t::the_log_t(): + mutex_(NULL) +{ + mutex_ = the_mutex_interface_t::create(); +} + +//----------------------------------------------------------------' +// the_log_t::~the_log_t +// +the_log_t::~the_log_t() +{ + mutex_->delete_this(); +} + +//---------------------------------------------------------------- +// the_log_t::log_no_lock +// +void +the_log_t::log_no_lock(std::ostream & (*f)(std::ostream &)) +{ + f(line_); +} + +//---------------------------------------------------------------- +// the_log_t::operator +// +the_log_t & +the_log_t::operator << (std::ostream & (*f)(std::ostream &)) +{ + the_lock_t lock(mutex_); + log_no_lock(f); + return *this; +} + +//---------------------------------------------------------------- +// the_log_t::precision +// +std::streamsize +the_log_t::precision() +{ + the_lock_t lock(mutex_); + std::streamsize p = line_.precision(); + return p; +} + +//---------------------------------------------------------------- +// the_log_t::precision +// +std::streamsize +the_log_t::precision(std::streamsize n) +{ + the_lock_t lock(mutex_); + std::streamsize p = line_.precision(n); + return p; +} + +//---------------------------------------------------------------- +// the_log_t::flags +// +std::ios::fmtflags +the_log_t::flags() const +{ + the_lock_t lock(mutex_); + std::ios::fmtflags f = line_.flags(); + return f; +} + +//---------------------------------------------------------------- +// the_log_t::flags +// +std::ios::fmtflags +the_log_t::flags(std::ios::fmtflags fmt) +{ + the_lock_t lock(mutex_); + std::ios::fmtflags f = line_.flags(fmt); + return f; +} + +//---------------------------------------------------------------- +// the_log_t::setf +// +void +the_log_t::setf(std::ios::fmtflags fmt) +{ + the_lock_t lock(mutex_); + line_.setf(fmt); +} + +//---------------------------------------------------------------- +// the_log_t::setf +// +void +the_log_t::setf(std::ios::fmtflags fmt, std::ios::fmtflags msk) +{ + the_lock_t lock(mutex_); + line_.setf(fmt, msk); +} + +//---------------------------------------------------------------- +// the_log_t::unsetf +// +void +the_log_t::unsetf(std::ios::fmtflags fmt) +{ + the_lock_t lock(mutex_); + line_.unsetf(fmt); +} + +//---------------------------------------------------------------- +// the_log_t::copyfmt +// +void +the_log_t::copyfmt(std::ostream & ostm) +{ + the_lock_t lock(mutex_); + line_.copyfmt(ostm); +} + + +//---------------------------------------------------------------- +// null_log +// +the_null_log_t * +null_log() +{ + static the_null_log_t * log = NULL; + if (log == NULL) + { + log = new the_null_log_t; + } + + return log; +} + + +//---------------------------------------------------------------- +// cerr_log +// +the_stream_log_t * +cerr_log() +{ + static the_stream_log_t * log = NULL; + if (log == NULL) + { + log = new the_stream_log_t(std::cerr); + } + + return log; +} + + +//---------------------------------------------------------------- +// cout_log +// +the_stream_log_t * +cout_log() +{ + static the_stream_log_t * log = NULL; + if (log == NULL) + { + log = new the_stream_log_t(std::cout); + } + + return log; +} diff --git a/src/IRMutex.cxx b/src/IRMutex.cxx new file mode 100644 index 0000000..c2fb8a2 --- /dev/null +++ b/src/IRMutex.cxx @@ -0,0 +1,113 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_boost_mutex.cxx +// Author : Pavel A. Koshevoy +// Created : Sat Oct 25 12:33:23 MDT 2008 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A thin wrapper Boost mutex class. + +// local includes: +#include "IRMutex.h" + +// system includes: +#include + +// namespace access: +using std::cerr; +using std::endl; + +//---------------------------------------------------------------- +// define +// +// #define DEBUG_MUTEX + + +//---------------------------------------------------------------- +// the_boost_mutex_t::the_boost_mutex_t +// +the_boost_mutex_t::the_boost_mutex_t(): + the_mutex_interface_t() +{} + +//---------------------------------------------------------------- +// the_boost_mutex_t::~the_boost_mutex_t +// +the_boost_mutex_t::~the_boost_mutex_t() +{} + +//---------------------------------------------------------------- +// the_boost_mutex_t::delete_this +// +void +the_boost_mutex_t::delete_this() +{ + delete this; +} + +//---------------------------------------------------------------- +// the_boost_mutex_t::create +// +the_mutex_interface_t * +the_boost_mutex_t::create() +{ + return new the_boost_mutex_t(); +} + +//---------------------------------------------------------------- +// the_boost_mutex_t::lock +// +void +the_boost_mutex_t::lock() +{ +#ifdef DEBUG_MUTEX + cerr << this << "\tlock" << endl; +#endif + + mutex_.lock(); +} + +//---------------------------------------------------------------- +// the_boost_mutex_t::unlock +// +void +the_boost_mutex_t::unlock() +{ +#ifdef DEBUG_MUTEX + cerr << this << "\tunlock" << endl; +#endif + + mutex_.unlock(); +} + +//---------------------------------------------------------------- +// the_boost_mutex_t::try_lock +// +bool +the_boost_mutex_t::try_lock() +{ +#ifdef DEBUG_MUTEX + cerr << this << "\ttry_lock" << endl; +#endif + + return mutex_.try_lock(); +} diff --git a/src/IRMutexInterface.cxx b/src/IRMutexInterface.cxx new file mode 100644 index 0000000..b0c4d4b --- /dev/null +++ b/src/IRMutexInterface.cxx @@ -0,0 +1,65 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_mutex_interface.cxx +// Author : Pavel A. Koshevoy +// Created : Wed Feb 21 08:22:00 MST 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : An abstract mutex class interface. + +// local includes: +#include "IRMutexInterface.h" + +// system includes: +#include + + +//---------------------------------------------------------------- +// the_mutex_interface_t::creator_ +// +the_mutex_interface_t::creator_t +the_mutex_interface_t::creator_ = NULL; + +//---------------------------------------------------------------- +// the_mutex_interface_t::~the_mutex_interface_t +// +the_mutex_interface_t::~the_mutex_interface_t() +{} + +//---------------------------------------------------------------- +// the_mutex_interface_t::set_creator +// +void +the_mutex_interface_t::set_creator(the_mutex_interface_t::creator_t creator) +{ + creator_ = creator; +} + +//---------------------------------------------------------------- +// the_mutex_interface_t::create +// +the_mutex_interface_t * +the_mutex_interface_t::create() +{ + if (creator_ == NULL) return NULL; + return creator_(); +} diff --git a/src/IRTerminator.cxx b/src/IRTerminator.cxx new file mode 100644 index 0000000..6f94671 --- /dev/null +++ b/src/IRTerminator.cxx @@ -0,0 +1,232 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_terminator.cxx +// Author : Pavel A. Koshevoy +// Created : Sun Sep 24 18:08:00 MDT 2006 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : a thread terminator convenience class + + +// local includes: +#include "IRTerminator.h" +#include "IRException.h" +#include "IRThreadStorage.h" + +// system includes: +#include +#include +#include +#include + +// namespace access: +using std::cerr; +using std::endl; + + +//---------------------------------------------------------------- +// DEBUG_TERMINATORS +// +// #define DEBUG_TERMINATORS + + +//---------------------------------------------------------------- +// the_terminator_t::the_terminator_t +// +the_terminator_t::the_terminator_t(const char * id): + id_(id), + termination_requested_(should_terminate_immediately()) +{ + terminate_on_request(); + add(this); +} + +//---------------------------------------------------------------- +// the_terminator_t:: +// +the_terminator_t::~the_terminator_t() +{ + del(this); +} + +//---------------------------------------------------------------- +// the_terminator_t::terminate +// +void +the_terminator_t::terminate() +{ +#ifdef DEBUG_TERMINATORS + cerr << "TERMINATE: " << id_ << ", this: " << this << endl; +#endif // DEBUG_TERMINATORS + termination_requested_ = true; +} + +//---------------------------------------------------------------- +// the_terminator_t::throw_exception +// +void +the_terminator_t::throw_exception() const +{ + // abort, termination requested: + std::ostringstream os; + os << id_ << " interrupted" << endl; + + the_exception_t e(os.str().c_str()); + throw e; +} + +//---------------------------------------------------------------- +// the_terminator_t::verify_termination +// +bool +the_terminator_t::verify_termination() +{ + return the_thread_storage().terminators().verify_termination(); +} + +//---------------------------------------------------------------- +// the_terminator_t::should_terminate_immediately +// +bool +the_terminator_t::should_terminate_immediately() +{ + return the_thread_storage().thread_stopped(); +} + +//---------------------------------------------------------------- +// the_terminator_t::add +// +void +the_terminator_t::add(the_terminator_t * terminator) +{ + the_thread_storage().terminators().add(terminator); +} + +//---------------------------------------------------------------- +// the_terminator_t::del +// +void +the_terminator_t::del(the_terminator_t * terminator) +{ + the_thread_storage().terminators().del(terminator); +} + + +//---------------------------------------------------------------- +// the_terminators_t::~the_terminators_t +// +the_terminators_t::~the_terminators_t() +{ + for (std::list::iterator i = terminators_.begin(); + i != terminators_.end(); ++i) + { + the_terminator_t * t = *i; + delete t; + } + + terminators_.clear(); +} + +//---------------------------------------------------------------- +// the_terminators_t::terminate +// +void +the_terminators_t::terminate() +{ + the_lock_t lock(*this); + +#ifdef DEBUG_TERMINATORS + cerr << endl << &terminators_ << ": terminate_all -- begin" << endl; +#endif // DEBUG_TERMINATORS + + for (std::list::iterator i = terminators_.begin(); + i != terminators_.end(); ++i) + { + the_terminator_t * t = *i; + t->terminate(); + } + +#ifdef DEBUG_TERMINATORS + cerr << &terminators_ << ": terminate_all -- end" << endl; +#endif // DEBUG_TERMINATORS +} + +//---------------------------------------------------------------- +// the_terminators_t::verify_termination +// +bool +the_terminators_t::verify_termination() +{ + the_lock_t lock(*this); + +#ifdef DEBUG_TERMINATORS + for (std::list::iterator iter = terminators_.begin(); + iter != terminators_.end(); ++iter) + { + the_terminator_t * t = *iter; + cerr << "ERROR: remaining terminators: " << t->id() << endl; + } +#endif + + return terminators_.empty(); +} + +//---------------------------------------------------------------- +// the_terminators_t::add +// +void +the_terminators_t::add(the_terminator_t * terminator) +{ + the_lock_t lock(*this); + terminators_.push_front(terminator); + +#ifdef DEBUG_TERMINATORS + cerr << &terminators_ << ": appended " << terminator->id() + << " terminator, addr " << terminator << endl; +#endif // DEBUG_TERMINATORS +} + +//---------------------------------------------------------------- +// the_terminators_t::del +// +void +the_terminators_t::del(the_terminator_t * terminator) +{ + the_lock_t lock(*this); + + std::list::iterator iter = + std::find(terminators_.begin(), terminators_.end(), terminator); + + if (iter == terminators_.end()) + { + cerr << "ERROR: no such terminator: " << terminator->id() << endl; + assert(0); + return; + } + + terminators_.erase(iter); + +#ifdef DEBUG_TERMINATORS + cerr << &terminators_ << ": removed " << terminator->id() + << " terminator, addr " << terminator << endl; +#endif // DEBUG_TERMINATORS +} diff --git a/src/IRThread.cxx b/src/IRThread.cxx new file mode 100644 index 0000000..cb6effa --- /dev/null +++ b/src/IRThread.cxx @@ -0,0 +1,196 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_boost_thread.cxx +// Author : Pavel A. Koshevoy +// Created : Sat Oct 25 12:33:09 MDT 2008 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A thin wrapper for Boost thread class. + +// local include: +#include "IRThread.h" +#include "IRThreadStorage.h" +#include "IRMutex.h" +#include "thread/the_mutex_interface.hxx" + +// system includes: +#include + +// namespace access: +using std::cout; +using std::cerr; +using std::endl; + +//---------------------------------------------------------------- +// DEBUG_THREAD +// +// #define DEBUG_THREAD + + +//---------------------------------------------------------------- +// THREAD_STORAGE +// +static the_boost_thread_storage_t THREAD_STORAGE; + +//---------------------------------------------------------------- +// the_boost_thread_t::the_boost_thread_t +// +the_boost_thread_t::the_boost_thread_t(): + the_thread_interface_t(the_boost_mutex_t::create()), + boost_thread_(NULL) +{ + if (THREAD_STORAGE.get() == NULL) + { + THREAD_STORAGE.reset(new the_thread_observer_t(*this)); + } +} + +//---------------------------------------------------------------- +// the_boost_thread_t::~the_boost_thread_t +// +the_boost_thread_t::~the_boost_thread_t() +{ + if (boost_thread_) + { + wait(); + } +} + +//---------------------------------------------------------------- +// the_boost_thread_t::delete_this +// +void +the_boost_thread_t::delete_this() +{ + delete this; +} + +//---------------------------------------------------------------- +// ImageProcessingThread::thread_storage +// +the_thread_storage_t & +the_boost_thread_t::thread_storage() +{ + return THREAD_STORAGE; +} + +//---------------------------------------------------------------- +// the_boost_thread_t::start +// +void +the_boost_thread_t::start() +{ + the_lock_t locker(mutex_); +#ifdef DEBUG_THREAD + cerr << "start of thread " << this << " requested" << endl; +#endif + + if (boost_thread_) + { + if (!stopped_) + { + // already running: +#ifdef DEBUG_THREAD + cerr << "thread " << this << " is already running" << endl; +#endif + return; + } + else + { + // wait for the shutdown to succeed, then restart the thread: +#ifdef DEBUG_THREAD + cerr << "waiting for thread " << this << " to shut down" << endl; +#endif + wait(); + } + } + +#ifdef DEBUG_THREAD + cerr << "starting thread " << this << endl; +#endif + + // we shouldn't have a Boost thread at this stage: + assert(!boost_thread_); + + // clear the termination flag: + stopped_ = false; + boost_thread_ = new boost::thread(callable_t(this)); +} + +//---------------------------------------------------------------- +// the_boost_thread_t::wait +// +void +the_boost_thread_t::wait() +{ + if (!boost_thread_) return; + + if (boost_thread_->get_id() == boost::this_thread::get_id()) + { + assert(false); + return; + } + + boost_thread_->join(); + delete boost_thread_; + boost_thread_ = NULL; +} + +//---------------------------------------------------------------- +// the_boost_thread_t::take_a_nap +// +void +the_boost_thread_t::take_a_nap(const unsigned long & microseconds) +{ + boost::xtime xt; + xt.sec = microseconds / 1000000; + xt.nsec = microseconds % 1000000; + boost::thread::sleep(xt); +} + +//---------------------------------------------------------------- +// the_boost_thread_t::terminators +// +the_terminators_t & +the_boost_thread_t::terminators() +{ + return terminators_; +} + +//---------------------------------------------------------------- +// the_boost_thread_t::run +// +void +the_boost_thread_t::run() +{ + // setup the thread storage: + { + the_lock_t locker(mutex_); + THREAD_STORAGE.reset(new the_thread_observer_t(*this)); + } + + // process the transactions: + work(); + + // clean up the thread storage: + THREAD_STORAGE.reset(NULL); +} diff --git a/src/IRThreadInterface.cxx b/src/IRThreadInterface.cxx new file mode 100644 index 0000000..182d4f9 --- /dev/null +++ b/src/IRThreadInterface.cxx @@ -0,0 +1,510 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_thread_interface.cxx +// Author : Pavel A. Koshevoy +// Created : Fri Feb 16 10:11:00 MST 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : An abstract thread class interface. + +// local includes: +#include "thread/the_thread_interface.hxx" +#include "thread/the_thread_pool.hxx" +#include "thread/the_mutex_interface.hxx" +#include "thread/the_terminator.hxx" +#include "thread/the_transaction.hxx" +#include "utils/the_utils.hxx" + +// system includes: +#include + +// namespace access: +using std::cout; +using std::cerr; +using std::endl; + +//---------------------------------------------------------------- +// DEBUG_THREAD +// +// #define DEBUG_THREAD + + +//---------------------------------------------------------------- +// MUTEX +// +static the_mutex_interface_t * MUTEX() +{ + static the_mutex_interface_t * mutex_ = NULL; + if (mutex_ == NULL) + { + mutex_ = the_mutex_interface_t::create(); + } + + return mutex_; +} + + +//---------------------------------------------------------------- +// the_thread_interface_t::creator_ +// +the_thread_interface_t::creator_t +the_thread_interface_t::creator_ = NULL; + +//---------------------------------------------------------------- +// the_thread_interface_t::the_thread_interface_t +// +the_thread_interface_t::the_thread_interface_t(the_mutex_interface_t * mutex): + mutex_(mutex), + stopped_(true), + sleep_when_idle_(false), + sleep_microsec_(10000), + active_transaction_(NULL), + thread_pool_(NULL), + thread_pool_cb_data_(NULL) +{ + the_lock_t locker(MUTEX()); + static unsigned int counter = 0; + id_ = counter++; +} + +//---------------------------------------------------------------- +// the_thread_interface_t::~the_thread_interface_t +// +the_thread_interface_t::~the_thread_interface_t() +{ + mutex_->delete_this(); +} + +//---------------------------------------------------------------- +// the_thread_interface_t::set_creator +// +void +the_thread_interface_t::set_creator(the_thread_interface_t::creator_t creator) +{ + creator_ = creator; +} + +//---------------------------------------------------------------- +// the_thread_interface_t::create +// +the_thread_interface_t * +the_thread_interface_t::create() +{ + if (creator_ == NULL) return NULL; + return creator_(); +} + +//---------------------------------------------------------------- +// the_thread_interface_t::set_mutex +// +void +the_thread_interface_t::set_mutex(the_mutex_interface_t * mutex) +{ + mutex_->delete_this(); + mutex_ = mutex; +} + +//---------------------------------------------------------------- +// the_thread_interface_t::set_idle_sleep_duration +// +void +the_thread_interface_t::set_idle_sleep_duration(bool enable, + unsigned int microseconds) +{ + sleep_when_idle_ = enable; + sleep_microsec_ = microseconds; +} + +//---------------------------------------------------------------- +// the_thread_interface_t::push_back +// +void +the_thread_interface_t::push_back(the_transaction_t * transaction) +{ + the_lock_t locker(mutex_); + transactions_.push_back(transaction); +} + +//---------------------------------------------------------------- +// the_thread_interface_t::push_front +// +void +the_thread_interface_t::push_front(the_transaction_t * transaction) +{ + the_lock_t locker(mutex_); + transactions_.push_front(transaction); +} + +//---------------------------------------------------------------- +// the_thread_interface_t::push_back +// +void +the_thread_interface_t::push_back(std::list & schedule) +{ + the_lock_t locker(mutex_); + transactions_.splice(transactions_.end(), schedule); +} + +//---------------------------------------------------------------- +// the_thread_interface_t::has_work +// +bool +the_thread_interface_t::has_work() const +{ + return (active_transaction_ != NULL) || !transactions_.empty(); +} + +//---------------------------------------------------------------- +// the_thread_interface_t::start +// +void +the_thread_interface_t::start(the_transaction_t * transaction) +{ + push_back(transaction); + start(); +} + +//---------------------------------------------------------------- +// the_thread_interface_t::stop +// +void +the_thread_interface_t::stop() +{ + the_lock_t locker(mutex_); + if (!stopped_) + { +#ifdef DEBUG_THREAD + cerr << "stopping thread: " << this << endl; +#endif + stopped_ = true; + terminators().terminate(); + } +} + +//---------------------------------------------------------------- +// the_thread_interface_t::flush +// +void +the_thread_interface_t::flush() +{ + the_lock_t locker(mutex_); + while (!transactions_.empty()) + { + the_transaction_t * t = remove_head(transactions_); +#ifdef DEBUG_THREAD + cerr << this << " flush " << t << endl; +#endif + t->notify(this, the_transaction_t::SKIPPED_E); + } +} + +//---------------------------------------------------------------- +// the_thread_interface_t::terminate_transactions +// +void +the_thread_interface_t::terminate_transactions() +{ + the_lock_t locker(mutex_); + + // remove any further pending transactions: + while (!transactions_.empty()) + { + the_transaction_t * t = remove_head(transactions_); +#ifdef DEBUG_THREAD + cerr << this << " terminate_transactions " << t << endl; +#endif + t->notify(this, the_transaction_t::SKIPPED_E); + } + + terminators().terminate(); +} + +//---------------------------------------------------------------- +// the_thread_interface_t::stop_and_go +// +void +the_thread_interface_t::stop_and_go(the_transaction_t * transaction) +{ + the_lock_t locker(mutex_); + + // remove any further pending transactions: + while (!transactions_.empty()) + { + the_transaction_t * t = remove_head(transactions_); +#ifdef DEBUG_THREAD + cerr << this << " stop_and_go (1) " << t << endl; +#endif + t->notify(this, the_transaction_t::SKIPPED_E); + } + + // terminate the currently executing transaction: + terminators().terminate(); + + // schedule the next transaction: + transactions_.push_back(transaction); + + // start the thread if it isn't already running: + start(); +} + +//---------------------------------------------------------------- +// the_thread_interface_t::stop_and_go +// +void +the_thread_interface_t::stop_and_go(std::list & schedule) +{ + the_lock_t locker(mutex_); + + // remove any further pending transactions: + while (!transactions_.empty()) + { + the_transaction_t * t = remove_head(transactions_); +#ifdef DEBUG_THREAD + cerr << this << " stop_and_go " << t << endl; +#endif + t->notify(this, the_transaction_t::SKIPPED_E); + } + + // terminate the currently executing transaction: + terminators().terminate(); + + // schedule the new transactions: + transactions_.splice(transactions_.end(), schedule); + + // start the thread if it isn't already running: + start(); +} + +//---------------------------------------------------------------- +// the_thread_interface_t::flush_and_go +// +void +the_thread_interface_t::flush_and_go(the_transaction_t * transaction) +{ + the_lock_t locker(mutex_); + + // remove any further pending transactions: + while (!transactions_.empty()) + { + the_transaction_t * t = remove_head(transactions_); +#ifdef DEBUG_THREAD + cerr << this << " flush_and_go (1) " << t << endl; +#endif + t->notify(this, the_transaction_t::SKIPPED_E); + } + + // schedule the next transaction: + transactions_.push_back(transaction); + + // start the thread if it isn't already running: + start(); +} + +//---------------------------------------------------------------- +// the_thread_interface_t::flush_and_go +// +void +the_thread_interface_t:: +flush_and_go(std::list & schedule) +{ + the_lock_t locker(mutex_); + + // remove any further pending transactions: + while (!transactions_.empty()) + { + the_transaction_t * t = remove_head(transactions_); +#ifdef DEBUG_THREAD + cerr << this << " flush_and_go " << t << endl; +#endif + t->notify(this, the_transaction_t::SKIPPED_E); + } + + // schedule the next transaction: + transactions_.splice(transactions_.end(), schedule); + + // start the thread if it isn't already running: + start(); +} + +//---------------------------------------------------------------- +// the_thread_interface_t::work +// +bool +the_thread_interface_t::work() +{ + the_lock_t lock_this(mutex_, false); + the_lock_t lock_pool(thread_pool_ ? + thread_pool_->mutex_ : NULL, + false); + bool all_transactions_completed = true; + + while (!stopped_) + { + // get the next transaction: + the_transaction_t * t = NULL; + { + lock_pool.arm(); + lock_this.arm(); + + if (thread_pool_ != NULL) + { + // call back the thread pool: +#ifdef DEBUG_THREAD + cerr << "thread " << this << " calling the pool" << endl; +#endif + thread_pool_->handle_thread(thread_pool_cb_data_); + } + + if (transactions_.empty()) + { + if (sleep_when_idle_ && !stopped_) + { + lock_this.disarm(); + lock_pool.disarm(); + take_a_nap(sleep_microsec_); +#ifdef DEBUG_THREAD + cerr << this << " sleeping" << endl; +#endif + continue; + } + else + { + // NOTE: the mutex will remain locked until the function returns: +#ifdef DEBUG_THREAD + cerr << "thread " << this << " is finishing up" << endl; +#endif + break; + } + } + else + { + t = remove_head(transactions_); +#ifdef DEBUG_THREAD + cerr << "thread " << this << " received " << t << endl; +#endif + lock_this.disarm(); + lock_pool.disarm(); + } + } + + try + { + active_transaction_ = t; + t->notify(this, the_transaction_t::STARTED_E); + t->execute(this); + all_transactions_completed = (all_transactions_completed && true); + active_transaction_ = NULL; + t->notify(this, the_transaction_t::DONE_E); + } + catch (std::exception & e) + { + active_transaction_ = NULL; +#ifdef DEBUG_THREAD + cerr << "FIXME: caught exception: " << e.what() << endl; +#endif + t->notify(this, + the_transaction_t::ABORTED_E, + e.what()); + } + catch (...) + { + active_transaction_ = NULL; +#ifdef DEBUG_THREAD + cerr << "FIXME: caught unknonwn exception" << endl; +#endif + t->notify(this, + the_transaction_t::ABORTED_E, + "unknown exception intercepted"); + } + } + + stopped_ = true; + + // make sure that all terminators have executed: +#ifndef NDEBUG + bool ok = the_terminator_t::verify_termination(); + assert(ok); +#endif + + all_transactions_completed = (all_transactions_completed && + transactions_.empty()); + + // abort pending transaction: + while (!transactions_.empty()) + { + the_transaction_t * t = remove_head(transactions_); +#ifdef DEBUG_THREAD + cerr << this << " work " << t << endl; +#endif + t->notify(this, the_transaction_t::SKIPPED_E); + } + +#ifdef DEBUG_THREAD + cerr << "thread " << this << " is finished" << endl; +#endif + + if (thread_pool_ != NULL) + { + lock_pool.arm(); + lock_this.arm(); + thread_pool_->handle_thread(thread_pool_cb_data_); + } + + return all_transactions_completed; +} + +//---------------------------------------------------------------- +// the_thread_interface_t::handle +// +void +the_thread_interface_t::handle(the_transaction_t * transaction, + the_transaction_t::state_t s) +{ + switch (s) + { + case the_transaction_t::SKIPPED_E: + case the_transaction_t::ABORTED_E: + case the_transaction_t::DONE_E: + delete transaction; + break; + + default: + break; + } +} + +//---------------------------------------------------------------- +// the_thread_interface_t::blab +// +void +the_thread_interface_t::blab(const char * message) const +{ + if (thread_pool_ == NULL) + { + cerr << message << endl; + } + else + { + thread_pool_->blab(message); + } +} + diff --git a/src/IRThreadPool.cxx b/src/IRThreadPool.cxx new file mode 100644 index 0000000..161b34e --- /dev/null +++ b/src/IRThreadPool.cxx @@ -0,0 +1,735 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_thread_pool.cxx +// Author : Pavel A. Koshevoy +// Created : Wed Feb 21 09:01:00 MST 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A thread pool class. + +// local includes: +#include "thread/the_thread_pool.hxx" +#include "thread/the_transaction.hxx" +#include "thread/the_mutex_interface.hxx" +#include "utils/the_utils.hxx" + +// system includes: +#include + +// namespace access: +using std::cout; +using std::cerr; +using std::endl; + +//---------------------------------------------------------------- +// DEBUG_THREAD +// +// #define DEBUG_THREAD + + +//---------------------------------------------------------------- +// the_transaction_wrapper_t::the_transaction_wrapper_t +// +the_transaction_wrapper_t:: +the_transaction_wrapper_t(const unsigned int & num_parts, + the_transaction_t * transaction): + mutex_(the_mutex_interface_t::create()), + cb_(transaction->notify_cb_), + cb_data_(transaction->notify_cb_data_), + num_parts_(num_parts) +{ + assert(mutex_ != NULL); + + for (unsigned int i = 0; i <= the_transaction_t::DONE_E; i++) + { + notified_[i] = 0; + } + transaction->set_notify_cb(notify_cb, this); +} + +//---------------------------------------------------------------- +// the_transaction_wrapper_t::~the_transaction_wrapper_t +// +the_transaction_wrapper_t::~the_transaction_wrapper_t() +{ + mutex_->delete_this(); + mutex_ = NULL; +} + +//---------------------------------------------------------------- +// the_transaction_wrapper_t::notify_cb +// +void +the_transaction_wrapper_t::notify_cb(void * data, + the_transaction_t * t, + the_transaction_t::state_t s) +{ + the_transaction_wrapper_t * wrapper = (the_transaction_wrapper_t *)(data); + bool done = wrapper->notify(t, s); + if (done) + { + delete wrapper; + } +} + +//---------------------------------------------------------------- +// the_transaction_wrapper_t::notify +// +bool +the_transaction_wrapper_t::notify(the_transaction_t * t, + the_transaction_t::state_t s) +{ + the_lock_t locker(mutex_); + notified_[s]++; + + unsigned int num_terminal_notifications = + notified_[the_transaction_t::SKIPPED_E] + + notified_[the_transaction_t::ABORTED_E] + + notified_[the_transaction_t::DONE_E]; + + switch (s) + { + case the_transaction_t::PENDING_E: + case the_transaction_t::STARTED_E: + if (notified_[s] == 1 && num_terminal_notifications == 0 && cb_ != NULL) + { + cb_(cb_data_, t, s); + } + return false; + + case the_transaction_t::SKIPPED_E: + case the_transaction_t::ABORTED_E: + case the_transaction_t::DONE_E: + { + if (num_terminal_notifications == num_parts_) + { + // restore the transaction callback: + t->set_notify_cb(cb_, cb_data_); + + // if necessary, consolidate the transaction state into one: + if (num_terminal_notifications != + notified_[the_transaction_t::DONE_E]) + { + the_transaction_t::state_t state = + notified_[the_transaction_t::ABORTED_E] > 0 ? + the_transaction_t::ABORTED_E : + the_transaction_t::SKIPPED_E; + t->set_state(state); + } + + if (cb_ != NULL) + { + cb_(cb_data_, t, t->state()); + } + else + { + delete t; + } + + return true; + } + return false; + } + + default: + assert(0); + } + + return false; +} + + +//---------------------------------------------------------------- +// the_thread_pool_t::the_thread_pool_t +// +the_thread_pool_t::the_thread_pool_t(unsigned int num_threads): + pool_size_(num_threads) +{ + mutex_ = the_mutex_interface_t::create(); + assert(mutex_ != NULL); + + pool_ = new the_thread_pool_data_t[num_threads]; + for (unsigned int i = 0; i < num_threads; i++) + { + pool_[i].id_ = i; + pool_[i].parent_ = this; + pool_[i].thread_ = the_thread_interface_t::create(); + assert(pool_[i].thread_ != NULL); + pool_[i].thread_->set_thread_pool_cb(this, &pool_[i]); + + // mark the thread as idle, initially: + idle_.push_back(i); + } +} + +//---------------------------------------------------------------- +// the_thread_pool_t::~the_thread_pool_t +// +the_thread_pool_t::~the_thread_pool_t() +{ + mutex_->delete_this(); + mutex_ = NULL; + + delete [] pool_; + pool_ = NULL; + pool_size_ = 0; +} + +//---------------------------------------------------------------- +// the_thread_pool_t::start +// +void +the_thread_pool_t::start() +{ + if (!transactions_.empty()) + { + while (true) + { + the_lock_t locker(mutex_); + if (transactions_.empty()) return; + + if (idle_.empty()) + { +#ifdef DEBUG_THREAD + // sanity test: + for (unsigned int i = 0; i < pool_size_; i++) + { + the_lock_t locker(pool_[i].thread_->mutex()); + if (!pool_[i].thread_->has_work()) + { + cerr << "WARNING: thread " << i << ", " << pool_[i].thread_ + << " is actually idle" << endl; + } + } +#endif + return; + } + + // get the next worker thread: + unsigned int id = remove_head(idle_); + busy_.push_back(id); + + // get the next transaction: + the_transaction_t * t = remove_head(transactions_); + + // start the thread: + thread(id)->start(t); + +#ifdef DEBUG_THREAD + cerr << "starting thread " << id << ", " << thread(id) << ", " << t + << endl; + + for (std::list::const_iterator i = idle_.begin(); + i != idle_.end(); ++i) + { + cerr << "idle: thread " << *i << ", " << thread(*i) << endl; + } + + for (std::list::const_iterator i = busy_.begin(); + i != busy_.end(); ++i) + { + cerr << "busy: thread " << *i << ", " << thread(*i) << endl; + } +#endif + } + } + else + { + // Transaction list is empty, perhaps they've already + // been distributed to the threads? Just start the threads + // and let them figure it out for themselves. + for (unsigned int i = 0; i < pool_size_; i++) + { + pool_[i].thread_->start(); + } + } +} + +//---------------------------------------------------------------- +// the_thread_pool_t::set_idle_sleep_duration +// +void +the_thread_pool_t::set_idle_sleep_duration(bool enable, unsigned int microsec) +{ + the_lock_t locker(mutex_); + for (unsigned int i = 0; i < pool_size_; i++) + { + pool_[i].thread_->set_idle_sleep_duration(enable, microsec); + } +} + +//---------------------------------------------------------------- +// notify_cb +// +static void +notify_cb(void * data, + the_transaction_t * transaction, + the_transaction_t::state_t s) +{ + the_thread_pool_t * pool = (the_thread_pool_t *)(data); + pool->handle(transaction, s); +} + +//---------------------------------------------------------------- +// status_cb +// +static void +status_cb(void * data, the_transaction_t *, const char * text) +{ + the_thread_pool_t * pool = (the_thread_pool_t *)(data); + pool->blab(text); +} + +//---------------------------------------------------------------- +// wrap +// +static void +wrap(the_transaction_t * transaction, + the_thread_pool_t * pool, + bool multithreaded) +{ + if (transaction->notify_cb() == NULL) + { + // override the callback: + transaction->set_notify_cb(::notify_cb, pool); + } + + if (transaction->status_cb() == NULL) + { + // override the callback: + transaction->set_status_cb(::status_cb, pool); + } + + if (multithreaded) + { + // silence the compiler warning: + // the_transaction_wrapper_t * wrapper = + new the_transaction_wrapper_t(pool->pool_size(), transaction); + } +} + +//---------------------------------------------------------------- +// the_thread_pool_t::push_front +// +void +the_thread_pool_t::push_front(the_transaction_t * transaction, + bool multithreaded) +{ + the_lock_t locker(mutex_); + no_lock_push_front(transaction, multithreaded); +} + +//---------------------------------------------------------------- +// the_thread_pool_t::push_back +// +void +the_thread_pool_t::push_back(the_transaction_t * transaction, + bool multithreaded) +{ + the_lock_t locker(mutex_); + no_lock_push_back(transaction, multithreaded); +} + +//---------------------------------------------------------------- +// the_thread_pool_t::push_back +// +void +the_thread_pool_t::push_back(std::list & schedule, + bool multithreaded) +{ + the_lock_t locker(mutex_); + no_lock_push_back(schedule, multithreaded); +} + +//---------------------------------------------------------------- +// the_thread_pool_t::pre_distribute_work +// +void +the_thread_pool_t::pre_distribute_work() +{ + the_lock_t locker(mutex_); + std::vector > split_schedule_(pool_size_); + + while (!transactions_.empty()) + { + for (unsigned int i = 0; i < pool_size_ && !transactions_.empty(); i++) + { + the_transaction_t * job = remove_head(transactions_); + split_schedule_[i].push_back(job); + } + } + + for (unsigned int i = 0; i < pool_size_; i++) + { + pool_[i].thread_->push_back(split_schedule_[i]); + } +} + +//---------------------------------------------------------------- +// the_thread_pool_t::has_work +// +bool +the_thread_pool_t::has_work() const +{ + the_lock_t locker(mutex_); + if (!transactions_.empty()) return true; + + // make sure the busy threads have work: + for (std::list::const_iterator + iter = busy_.begin(); iter != busy_.end(); ++iter) + { + unsigned int i = *iter; + if (pool_[i].thread_->has_work()) return true; + } + + return false; +} + +//---------------------------------------------------------------- +// the_thread_pool_t::start +// +void +the_thread_pool_t::start(the_transaction_t * transaction, bool multithreaded) +{ + push_back(transaction, multithreaded); + start(); +} + +//---------------------------------------------------------------- +// the_thread_pool_t::stop +// +void +the_thread_pool_t::stop() +{ + // remove any pending transactions: + flush(); + + the_lock_t locker(mutex_); + for (unsigned int i = 0; i < pool_size_; i++) + { + pool_[i].thread_->stop(); + } +} + +//---------------------------------------------------------------- +// the_thread_pool_t::wait +// +void +the_thread_pool_t::wait() +{ + // the_lock_t locker(mutex_); + for (unsigned int i = 0; i < pool_size_; i++) + { + pool_[i].thread_->wait(); + } +} + +//---------------------------------------------------------------- +// the_thread_pool_t::flush +// +void +the_thread_pool_t::flush() +{ + the_lock_t locker(mutex_); + no_lock_flush(); +} + +//---------------------------------------------------------------- +// the_thread_pool_t::terminate_transactions +// +void +the_thread_pool_t::terminate_transactions() +{ + the_lock_t locker(mutex_); + no_lock_terminate_transactions(); +} + +//---------------------------------------------------------------- +// the_thread_pool_t::stop_and_go +// +void +the_thread_pool_t::stop_and_go(the_transaction_t * transaction, + bool multithreaded) +{ + // safely manipulate the transactions: + { + the_lock_t locker(mutex_); + no_lock_terminate_transactions(); + no_lock_push_back(transaction, multithreaded); + } + + start(); +} + +//---------------------------------------------------------------- +// the_thread_pool_t::stop_and_go +// +void +the_thread_pool_t::stop_and_go(std::list & schedule, + bool multithreaded) +{ + // safely manipulate the transactions: + { + the_lock_t locker(mutex_); + no_lock_terminate_transactions(); + no_lock_push_back(schedule, multithreaded); + } + + start(); +} + +//---------------------------------------------------------------- +// the_thread_pool_t::flush_and_go +// +void +the_thread_pool_t::flush_and_go(the_transaction_t * transaction, + bool multithreaded) +{ + // safely manipulate the transactions: + { + the_lock_t locker(mutex_); + no_lock_flush(); + no_lock_push_back(transaction, multithreaded); + } + + start(); +} + +//---------------------------------------------------------------- +// the_thread_pool_t::flush_and_go +// +void +the_thread_pool_t::flush_and_go(std::list & schedule, + bool multithreaded) +{ + // safely manipulate the transactions: + { + the_lock_t locker(mutex_); + no_lock_flush(); + no_lock_push_back(schedule, multithreaded); + } + + start(); +} + +//---------------------------------------------------------------- +// the_thread_pool_t::blab +// +void +the_thread_pool_t::handle(the_transaction_t * transaction, + the_transaction_t::state_t s) +{ + switch (s) + { + case the_transaction_t::SKIPPED_E: + case the_transaction_t::ABORTED_E: + case the_transaction_t::DONE_E: + delete transaction; + break; + + default: + break; + } +} + +//---------------------------------------------------------------- +// the_thread_pool_t::handle +// +void +the_thread_pool_t::blab(const char * message) const +{ + cerr << message << endl; +} + +//---------------------------------------------------------------- +// the_thread_pool_t::handle_thread +// +// prior to calling this function the thread interface locks +// the pool mutex and it's own mutex -- don't try locking these +// again while inside this callback, and don't call any other +// pool or thread member functions which lock either of these +// because that would result in a deadlock +// +void +the_thread_pool_t::handle_thread(the_thread_pool_data_t * data) +{ + const unsigned int & id = data->id_; + the_thread_interface_t * t = thread(id); + + if (t->stopped_) + { +#ifdef DEBUG_THREAD + cerr << "thread has stopped: " << t << endl; +#endif + + if (!has(idle_, id)) + { + idle_.push_back(id); + } + + if (has(busy_, id)) + { + busy_.remove(id); + } + + return; + } + + if (t->has_work()) + { +#ifdef DEBUG_THREAD + cerr << "thread still has work: " << t << endl; +#endif + return; + } + + if (transactions_.empty()) + { + // tell the thread to stop: + t->stopped_ = true; + + if (!has(idle_, id)) + { + idle_.push_back(id); + } + + if (has(busy_, id)) + { + busy_.remove(id); + } + +#ifdef DEBUG_THREAD + cerr << "no more work for thread: " << t << endl; + + for (std::list::const_iterator i = idle_.begin(); + i != idle_.end(); ++i) + { + cerr << "idle: thread " << *i << ", " << thread(*i) << endl; + } + + for (std::list::const_iterator i = busy_.begin(); + i != busy_.end(); ++i) + { + cerr << "busy: thread " << *i << ", " << thread(*i) << endl; + } +#endif + + return; + } + + // execute another transaction: + the_transaction_t * transaction = remove_head(transactions_); +#ifdef DEBUG_THREAD + cerr << "giving " << transaction << " to thread " << t << endl; +#endif + + t->transactions_.push_back(transaction); +} + +//---------------------------------------------------------------- +// the_thread_pool_t::no_lock_flush +// +void +the_thread_pool_t::no_lock_flush() +{ + while (!transactions_.empty()) + { + the_transaction_t * t = remove_head(transactions_); + t->notify(this, the_transaction_t::SKIPPED_E); + } +} + +//---------------------------------------------------------------- +// the_thread_pool_t::no_lock_terminate_transactions +// +void +the_thread_pool_t::no_lock_terminate_transactions() +{ + // remove any pending transactions: + no_lock_flush(); + + for (unsigned int i = 0; i < pool_size_; i++) + { + pool_[i].thread_->terminate_transactions(); + } +} + +//---------------------------------------------------------------- +// the_thread_pool_t::no_lock_push_front +// +void +the_thread_pool_t::no_lock_push_front(the_transaction_t * transaction, + bool multithreaded) +{ + wrap(transaction, this, multithreaded); + transactions_.push_front(transaction); + for (unsigned int i = 1; multithreaded && i < pool_size_; i++) + { + transactions_.push_front(transaction); + } +} + +//---------------------------------------------------------------- +// the_thread_pool_t::no_lock_push_back +// +void +the_thread_pool_t::no_lock_push_back(the_transaction_t * transaction, + bool multithreaded) +{ + wrap(transaction, this, multithreaded); + transactions_.push_back(transaction); + for (unsigned int i = 1; multithreaded && i < pool_size_; i++) + { + transactions_.push_back(transaction); + } +} + +//---------------------------------------------------------------- +// the_thread_pool_t::no_lock_push_back +// +void +the_thread_pool_t::no_lock_push_back(std::list & schedule, + bool multithreaded) +{ + // preprocess the transactions: + for (std::list::iterator i = schedule.begin(); + i != schedule.end(); i++) + { + wrap(*i, this, multithreaded); + } + + if (multithreaded) + { + while (!schedule.empty()) + { + the_transaction_t * t = remove_head(schedule); + for (unsigned int i = 0; i < pool_size_; i++) + { + transactions_.push_back(t); + } + } + } + else + { + transactions_.splice(transactions_.end(), schedule); + } +} diff --git a/src/IRThreadStorage.cxx b/src/IRThreadStorage.cxx new file mode 100644 index 0000000..3142631 --- /dev/null +++ b/src/IRThreadStorage.cxx @@ -0,0 +1,102 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_thread_storage.cxx +// Author : Pavel A. Koshevoy +// Created : Fri Jan 9 12:27:00 MDT 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A thread storage abstract interface class + +// local includes: +#include +#include + + +//---------------------------------------------------------------- +// the_dummy_terminators_t +// +class the_dummy_terminators_t : public the_terminators_t +{ +public: + // virtual: + void lock() {} + void unlock() {} +}; + +//---------------------------------------------------------------- +// the_dummy_thread_storage_t +// +class the_dummy_thread_storage_t : public the_thread_storage_t +{ +public: + // virtual: + bool is_ready() const + { return true; } + + bool thread_stopped() const + { return false; } + + the_terminators_t & terminators() + { return terminators_; } + + unsigned int thread_id() const + { return ~0u; } + +private: + the_dummy_terminators_t terminators_; +}; + +//---------------------------------------------------------------- +// the_dummy_thread_storage +// +static the_thread_storage_t & +the_dummy_thread_storage() +{ + static the_dummy_thread_storage_t thread_storage; + return thread_storage; +} + +//---------------------------------------------------------------- +// thread_storage_provider_ +// +static the_thread_storage_provider_t +thread_storage_provider_ = the_dummy_thread_storage; + +//---------------------------------------------------------------- +// set_the_thread_storage_provider +// +the_thread_storage_provider_t +set_the_thread_storage_provider(the_thread_storage_provider_t p) +{ + the_thread_storage_provider_t old = thread_storage_provider_; + thread_storage_provider_ = p; + return old; +} + +//---------------------------------------------------------------- +// the_thread_storage +// +the_thread_storage_t & +the_thread_storage() +{ + return thread_storage_provider_(); +} diff --git a/src/IRTransaction.cxx b/src/IRTransaction.cxx new file mode 100644 index 0000000..4116fa9 --- /dev/null +++ b/src/IRTransaction.cxx @@ -0,0 +1,183 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_transaction.cxx +// Author : Pavel A. Koshevoy +// Created : Fri Feb 16 09:52:00 MST 2007 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : A thread transaction class. + +// system includes: +#include +#include + +// local includes: +#include "thread/the_transaction.hxx" +#include "thread/the_thread_interface.hxx" +#include "thread/the_mutex_interface.hxx" + + +//---------------------------------------------------------------- +// the_transaction_t::the_transaction_t +// +the_transaction_t::the_transaction_t(): + mutex_(the_mutex_interface_t::create()), + request_(NOTHING_E), + state_(PENDING_E), + notify_cb_(NULL), + notify_cb_data_(NULL), + status_cb_(NULL), + status_cb_data_(NULL) +{} + +//---------------------------------------------------------------- +// the_transaction_t::~the_transaction_t +// +the_transaction_t::~the_transaction_t() +{ + if (mutex_ != NULL) + { + mutex_->delete_this(); + mutex_ = NULL; + } +} + +//---------------------------------------------------------------- +// the_transaction_t::notify +// +void +the_transaction_t::notify(the_transaction_handler_t * handler, + state_t s, + const char * message) +{ + set_state(s); + blab(handler, message); + + if (notify_cb_ == NULL) + { + handler->handle(this, s); + } + else + { + notify_cb_(notify_cb_data_, this, s); + } +} + +//---------------------------------------------------------------- +// the_transaction_t::blab +// +void +the_transaction_t::blab(the_transaction_handler_t * handler, + const char * message) +{ + if (message == NULL) return; + + if (status_cb_ == NULL) + { + handler->blab(message); + } + else + { + status_cb_(status_cb_data_, this, message); + } +} + +//---------------------------------------------------------------- +// the_transaction_t::callback_request +// +bool +the_transaction_t::callback_request() +{ + if (status_cb_ == NULL) + { + return false; + } + + // change the request state: + { + the_lock_t locker(mutex_); + request_ = WAITING_E; + } + + // execute the status callback: + status_cb_(status_cb_data_, this, NULL); + + while (true) + { + sleep_msec(100); + + // check the request state: + the_lock_t locker(mutex_); + if (request_ == NOTHING_E) + { + break; + } + } + + return true; +} + +//---------------------------------------------------------------- +// the_transaction_t::callback +// +void +the_transaction_t::callback() +{ + // remove the callback request: + the_lock_t locker(mutex_); + request_ = NOTHING_E; +} + +//---------------------------------------------------------------- +// operator << +// +std::ostream & +operator << (std::ostream & so, const the_transaction_t::state_t & state) +{ + switch (state) + { + case the_transaction_t::PENDING_E: + so << "pending"; + return so; + + case the_transaction_t::SKIPPED_E: + so << "skipped"; + return so; + + case the_transaction_t::STARTED_E: + so << "started"; + return so; + + case the_transaction_t::ABORTED_E: + so << "aborted"; + return so; + + case the_transaction_t::DONE_E: + so << "done"; + return so; + + default: + so << int(state); + assert(0); + return so; + } +} diff --git a/src/itkIRCommon.cxx b/src/itkIRCommon.cxx new file mode 100644 index 0000000..9ac6b1a --- /dev/null +++ b/src/itkIRCommon.cxx @@ -0,0 +1,1673 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : common.cxx +// Author : Pavel A. Koshevoy +// Created : 2005/03/24 16:53 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : Helper functions for mosaicing, image warping, +// image preprocessing, convenience wrappers for +// ITK file and ITK filters. + +// local includes: +#include "itkIRCommon.h" +#include "itk/itkLegendrePolynomialTransform.h" +#include "itk/itkGridTransform.h" +#include "itk/itkMeshTransform.h" +#include "itk/itkRBFTransform.h" +#include "itk/itkRadialDistortionTransform.h" + +// ITK includes: +#include +#include +#include + +// system includes: +#include + +// namespace access: +using std::endl; +using std::flush; +using std::cerr; +using std::ios; +using std::setw; + +#if defined (__APPLE__) + #include // hack necessary to get omp to link even when there aren't any omp pragmas around. Osx-only. See http://stackoverflow.com/questions/8983038/linker-errors-after-enabling-openmp-on-mac. + pthread_attr_t gomp_thread_attr; +#endif + +//---------------------------------------------------------------- +// BREAK +// +int BREAK(unsigned int i) +{ + printf("\n\n\nBREAK BREAK BREAK BREAK BREAK BREAK BREAK: %i\n\n\n\n", i); + return 0; +} + +//---------------------------------------------------------------- +// register_transforms +// +static void +register_transforms() +{ + // try to avoid re-registering the same transforms over and over again: + static bool transforms_registered = false; + if (transforms_registered) return; + + // make sure the transforms I care about are known to the object factory: + itk::TransformFactoryBase::RegisterDefaultTransforms(); + + itk::TransformFactory >:: + RegisterTransform(); + + itk::TransformFactory >:: + RegisterTransform(); + + itk::TransformFactory >:: + RegisterTransform(); + + itk::TransformFactory >:: + RegisterTransform(); + + itk::TransformFactory >:: + RegisterTransform(); + + itk::TransformFactory >:: + RegisterTransform(); + + itk::TransformFactory >:: + RegisterTransform(); + + itk::TransformFactory >:: + RegisterTransform(); + + itk::TransformFactory:: + RegisterTransform(); + + itk::TransformFactory:: + RegisterTransform(); + + itk::TransformFactory< itk::RBFTransform>:: + RegisterTransform(); + + itk::TransformFactory< itk::RadialDistortionTransform >:: + RegisterTransform(); + + transforms_registered = true; +} + +//---------------------------------------------------------------- +// save_transform +// +void +save_transform(std::ostream & so, const itk::TransformBase * t) +{ + typedef itk::TransformBase::ParametersType params_t; + + // make sure the transforms I care about are known to the object factory: + register_transforms(); + + // save the transform type: + so << t->GetTransformTypeAsString(); + + // get the stream precision: + int prec = so.precision(); + + // save the variable parameters: + params_t vp = t->GetParameters(); + const unsigned int num_vp = vp.size(); + so << " vp " << num_vp; + for (unsigned int i = 0; i < num_vp; i++) + { + so << ' ' << setw(prec + 7) << vp[i]; + } + + // save the fixed parameters: + try + { + params_t fp = t->GetFixedParameters(); + + const unsigned int num_fp = fp.size(); + so << " fp " << num_fp; + for (unsigned int i = 0; i < num_fp; i++) + { + so << ' ' << setw(prec + 7) << fp[i]; + } + } + catch (itk::ExceptionObject &) + { + so << " no_fp"; + } +} + +//---------------------------------------------------------------- +// load_transform +// +// Load an ITK transform of specified type from a stream. +// +itk::TransformBase::Pointer +load_transform(std::istream & si, const std::string & transform_type) +{ + typedef itk::TransformBase::ParametersType params_t; + + // make sure the transforms I care about are known to the object factory: + register_transforms(); + + // FIXME: +#if 0 + std::list names = + itk::TransformFactoryBase::GetFactory()->GetClassOverrideWithNames(); + for (std::list::iterator it = names.begin(); + it != names.end(); + ++it) + { + cerr << "FIXME: " << *it << endl; + } +#endif + + itk::LightObject::Pointer tmp = + itk::ObjectFactoryBase::CreateInstance(transform_type.c_str()); + + itk::TransformBase::Pointer t = + dynamic_cast(tmp.GetPointer()); + + if (t.GetPointer() == NULL) + { + cerr << "could not instantiate " << transform_type + << ", giving up ..." << endl; + return t; + } + t->Register(); + + // load the variable parameters: + params_t vp; + std::string vp_token; + si >> vp_token; + if (vp_token != "vp") + { + cerr << "expected vp, received '" << vp_token + << "', aborting ..." << endl; + assert(false); + } + + unsigned int num_vp = 0; + si >> num_vp; + vp.SetSize(num_vp); + + for (unsigned int i = 0; i < num_vp; i++) + { +#ifdef __APPLE__ + // NOTE: OSX std::iostream is broken -- it can write out doubles + // so small that it can't read them back in. The workaround is + // to read them in as long doubles, and then cast them down to + // a double: + long double tmp = 0; + si >> tmp; + if (si.fail()) si.clear(); + vp[i] = double(tmp); +#else + si >> vp[i]; +#endif + } + + // load the fixed parameters: + params_t fp; + std::string fp_token; + si >> fp_token; + if (fp_token == "fp") + { + fp = t->GetFixedParameters(); + unsigned int num_fp = 0; + si >> num_fp; + + if (num_fp != fp.size()) + { + // some transforms (RBF) have variable number of fixed parameters: + fp.SetSize(num_fp); + } + + for (unsigned int i = 0; i < num_fp; i++) + { +#ifdef __APPLE__ + // NOTE: OSX std::iostream is broken -- it can write out doubles + // so small that it can't read them back in. The workaround is + // to read them in as long doubles, and then cast them down to + // a double: + long double tmp = 0; + si >> tmp; + if (si.fail()) si.clear(); + fp[i] = double(tmp); +#else + si >> fp[i]; +#endif + } + } + else if (fp_token != "no_fp") + { + cerr << "unexpected token: '" << fp_token + << "', aborting ..." << endl; + assert(false); + } + + // set the fixed parameters first: + try + { + t->SetFixedParameters(fp); + } + catch (itk::ExceptionObject &) + {} + + // set the variable parameters: + t->SetParametersByValue(vp); + + return t; +} + +//---------------------------------------------------------------- +// load_transform +// +itk::TransformBase::Pointer +load_transform(std::istream & si) +{ + // load the transform type: + std::string transform_type; + si >> transform_type; + + return load_transform(si, transform_type); +} + + +//---------------------------------------------------------------- +// save_mosaic +// +// Save image filenames and associated ITK transforms to a stream. +// +void +save_mosaic(std::ostream & so, + const unsigned int & num_images, + const double & pixel_spacing, + const bool & use_std_mask, + const the_text_t * image, + const itk::TransformBase ** transform) +{ + ios::fmtflags old_flags = so.setf(ios::scientific); + int old_precision = so.precision(); + so.precision(12); + + so << "format_version_number: " << 1 << endl + << "number_of_images: " << num_images << endl + << "pixel_spacing: " << pixel_spacing << endl + << "use_std_mask: " << use_std_mask << endl; + + for (unsigned int i = 0; i < num_images; i++) + { + // Find just the filename to save out to the mosaic. + so << "image:" << endl << image[i] << endl; + save_transform(so, transform[i]); + so << endl; + } + + so.setf(old_flags); + so.precision(old_precision); +} + +//---------------------------------------------------------------- +// load_mosaic +// +// Load image filenames and associated ITK transforms from a stream. +// +void +load_mosaic(std::istream & si, + double & pixel_spacing, + bool & use_std_mask, + std::vector & image, + std::vector & transform) +{ + // make sure the transforms I care about are known to the object factory: + register_transforms(); + + unsigned int num_images = ~0; + unsigned int i = 0; + unsigned int version = 0; + + // for backwards compatibility assume the standard mask is used: + use_std_mask = true; + + while (si.eof() == false && num_images != i) + { + std::string token; + si >> token; + + if (token == "format_version_number:") + { + si >> version; + } + else if (token == "number_of_images:") + { + si >> num_images; + image.resize(num_images); + transform.resize(num_images); + } + else if (token == "pixel_spacing:") + { + si >> pixel_spacing; + } + else if (token == "use_std_mask:") + { + si >> use_std_mask; + } + else if (token == "image:") + { + if (version == 0) + { + si >> image[i]; + } + else + { + si >> std::ws; + getline(si, image[i]); + } + + // the next token should be the transform type string: + std::string next_token; + si >> next_token; + + if (version == 0) + { + // Original .mosaic file format kept the filename + // and the transform type string on the same line, + // which made filenames with spaces difficult to handle: + + while (si.eof() == false) + { + itk::LightObject::Pointer tmp = + itk::ObjectFactoryBase::CreateInstance(next_token.c_str()); + + if (tmp.GetPointer()) + { + break; + } + + // NOTE: this is a dirty hack to work around filenames with + // spaces problem in the original file format. + + // assume the string that was just read in is a part of filename: + image[i] += the_text_t(" "); + image[i] += the_text_t(next_token.c_str()); + si >> next_token; + } + } + + transform[i] = load_transform(si, next_token); + i++; + } + else + { + cerr << "WARNING: unknown token: '" << token << "', ignoring ..." << endl; + } + } +} + +//---------------------------------------------------------------- +// is_empty_bbox +// +// Test whether a bounding box is empty (min > max) +// +bool +is_empty_bbox(const pnt2d_t & min, + const pnt2d_t & max) +{ + return min[0] > max[0] || min[1] > max[1]; +} + +//---------------------------------------------------------------- +// is_singular_bbox +// +// Test whether a bounding box is singular (min == max) +// +bool +is_singular_bbox(const pnt2d_t & min, + const pnt2d_t & max) +{ + return min == max; +} + +//---------------------------------------------------------------- +// clamp_bbox +// +// Restrict a bounding box to be within given limits. +// +void +clamp_bbox(const pnt2d_t & confines_min, + const pnt2d_t & confines_max, + pnt2d_t & min, + pnt2d_t & max) +{ + if (!is_empty_bbox(confines_min, confines_max)) + { + min[0] = std::min(confines_max[0], std::max(confines_min[0], min[0])); + min[1] = std::min(confines_max[1], std::max(confines_min[1], min[1])); + + max[0] = std::min(confines_max[0], std::max(confines_min[0], max[0])); + max[1] = std::min(confines_max[1], std::max(confines_min[1], max[1])); + } +} + + +//---------------------------------------------------------------- +// clip_histogram +// +// This is used by CLAHE to limit the contrast ratio slope. +// +void +clip_histogram(const double & clipping_height, + const unsigned int & pdf_size, + const unsigned int * pdf, + double * clipped_pdf) +{ + // count the number of pixels exceeding the clipping height: + double excess = double(0); + double deficit = double(0); + for (unsigned int i = 0; i < pdf_size; i++) + { + if (double(pdf[i]) <= clipping_height) + { + deficit += clipping_height - double(pdf[i]); + clipped_pdf[i] = double(pdf[i]); + } + else + { + excess += double(pdf[i]) - clipping_height; + clipped_pdf[i] = clipping_height; + } + } + + const double height_increment = excess / deficit; + for (unsigned int i = 0; i < pdf_size; i++) + { + if (clipped_pdf[i] >= clipping_height) continue; + + double d = clipping_height - clipped_pdf[i]; + clipped_pdf[i] += d * height_increment; + } +} + +//---------------------------------------------------------------- +// std_mask +// +// Given image dimensions, generate a mask image. +// use_std_mask -- the standard mask for the Robert Marc Lab data. +// +mask_t::Pointer +std_mask(const itk::Point & origin, + const itk::Vector & spacing, + const itk::Size<2> & sz, + bool use_std_mask) +{ + // setup the mask: + mask_t::Pointer mask = make_image(spacing, sz, 255); + mask->SetOrigin(origin); + + if (use_std_mask) + { + unsigned int roi_x = 0; + unsigned int roi_y = (unsigned int)(0.97 * sz[1]); + unsigned int roi_w = (unsigned int)(0.3025 * sz[0]); + unsigned int roi_h = sz[1] - roi_y; + ::fill(mask, roi_x, roi_y, roi_w, roi_h, 0); + + /* + roi_y = (unsigned int)(0.98 * sz[1]); + roi_w = sz[0]; + roi_h = sz[1] - roi_y; + fill(mask, roi_x, roi_y, roi_w, roi_h, 0); + */ + } + + return mask; +} + + +//---------------------------------------------------------------- +// to_wcs +// +inline static vec2d_t to_wcs(const vec2d_t & u_axis, + const vec2d_t & v_axis, + const double & u, + const double & v) +{ + return u_axis * u + v_axis * v; +} + + +//---------------------------------------------------------------- +// to_wcs +// +inline static pnt2d_t to_wcs(const pnt2d_t & origin, + const vec2d_t & u_axis, + const vec2d_t & v_axis, + const double & u, + const double & v) +{ + return origin + to_wcs(u_axis, v_axis, u, v); +} + +//---------------------------------------------------------------- +// calc_tile_mosaic_bbox +// +bool + calc_tile_mosaic_bbox(const base_transform_t * mosaic_to_tile, + + // image space bounding boxes of the tile: + const pnt2d_t & tile_min, + const pnt2d_t & tile_max, + + // mosaic space bounding boxes of the tile: + pnt2d_t & mosaic_min, + pnt2d_t & mosaic_max, + + // sample points along the image edges: + const unsigned int np) +{ + // initialize an empty bounding box: + mosaic_min[0] = std::numeric_limits::max(); + mosaic_min[1] = mosaic_min[0]; + mosaic_max[0] = -mosaic_min[0]; + mosaic_max[1] = -mosaic_min[0]; + + // it happens: + if (tile_min[0] == std::numeric_limits::max() || !mosaic_to_tile) + { + return true; + } + + base_transform_t::Pointer tile_to_mosaic = mosaic_to_tile->GetInverse(); + if (tile_to_mosaic.GetPointer() == NULL) + { + return false; + } + + double W = tile_max[0] - tile_min[0]; + double H = tile_max[1] - tile_min[1]; + + // a temporary vector to hold the sample points: + std::vector xy((np + 1) * 4); + + // corner points: + xy[0] = pnt2d(tile_min[0], tile_min[1]); + xy[1] = pnt2d(tile_min[0], tile_max[1]); + xy[2] = pnt2d(tile_max[0], tile_min[1]); + xy[3] = pnt2d(tile_max[0], tile_max[1]); + + // edge points: + for (unsigned int j = 0; j < np; j++) + { + const double t = double(j + 1) / double(np + 1); + double x = tile_min[0] + t * W; + double y = tile_min[1] + t * H; + + unsigned int offset = (j + 1) * 4; + xy[offset + 0] = pnt2d(x, tile_min[1]); + xy[offset + 1] = pnt2d(x, tile_max[1]); + xy[offset + 2] = pnt2d(tile_min[0], y); + xy[offset + 3] = pnt2d(tile_max[0], y); + } + + // find the inverse mapping for each point, if possible: + std::list uv_list; + for (unsigned int j = 0; j < xy.size(); j++) + { + pnt2d_t uv; + if (find_inverse(tile_min, + tile_max, + mosaic_to_tile, + tile_to_mosaic.GetPointer(), + xy[j], + uv)) + { + uv_list.push_back(uv); + } + } + + // calculate the bounding box of the inverse-mapped points: + for (std::list::const_iterator iter = uv_list.begin(); + iter != uv_list.end(); ++iter) + { + const pnt2d_t & uv = *iter; + mosaic_min[0] = std::min(mosaic_min[0], uv[0]); + mosaic_max[0] = std::max(mosaic_max[0], uv[0]); + mosaic_min[1] = std::min(mosaic_min[1], uv[1]); + mosaic_max[1] = std::max(mosaic_max[1], uv[1]); + } + + return !uv_list.empty(); +} + + +//---------------------------------------------------------------- +// make_colors +// +// Generate a scrambled rainbow colormap. +// +void +make_colors(const unsigned int & num_colors, std::vector & color) +{ + static const xyz_t EAST = xyz(1, 0, 0); + static const xyz_t NORTH = xyz(0, 1, 0); + static const xyz_t WEST = xyz(0, 0, 1); + static const xyz_t SOUTH = xyz(0, 0, 0); + + color.resize(num_colors); + + for (unsigned int i = 0; i < num_colors; i++) + { +#if 1 + double t = fmod(double(i % 2) / 2.0 + + double(i) / double(num_colors - 1), 1.0); +#else + double t = double(i) / double(num_colors); +#endif + + double s = 0.5 + 0.5 * fmod(double((i + 1) % 3) / 3.0 + + double(i) / double(num_colors - 1), 1.0); + color[i] = hsv_to_rgb(xyz(t, s, 1.0)) * 255.0; + } +} + +//---------------------------------------------------------------- +// find_inverse +// +// Given a forward transform and image space coordinates, +// find mosaic space coordinates. +// +bool +find_inverse(const pnt2d_t & tile_min, // tile space + const pnt2d_t & tile_max, // tile space + const base_transform_t * mosaic_to_tile, + const pnt2d_t & xy, // tile space + pnt2d_t & uv, // mosaic space + const unsigned int max_iterations, + const double min_step_scale, + const double min_error_sqrd, + const unsigned int pick_up_pace_steps) +{ + base_transform_t::Pointer tile_to_mosaic = mosaic_to_tile->GetInverse(); + return find_inverse(tile_min, + tile_max, + mosaic_to_tile, + tile_to_mosaic.GetPointer(), + xy, + uv, + max_iterations, + min_step_scale, + min_error_sqrd); +} + +//---------------------------------------------------------------- +// find_inverse +// +// Given a forward transform, an approximate inverse transform, +// and image space coordinates, find mosaic space coordinates. +// +bool +find_inverse(const pnt2d_t & tile_min, // tile space + const pnt2d_t & tile_max, // tile space + const base_transform_t * mosaic_to_tile, + const base_transform_t * tile_to_mosaic, + const pnt2d_t & xy, // tile space + pnt2d_t & uv, // mosaic space + const unsigned int max_iterations, + const double min_step_scale, + const double min_error_sqrd, + const unsigned int pick_up_pace_steps) +{ + // #define DEBUG_FIND_INVERSE + if (tile_to_mosaic == NULL) + { + return false; + } + + const itk::GridTransform *gridTransform = + dynamic_cast(tile_to_mosaic); + if (gridTransform != NULL) + { + // special case for the grid transform -- the inverse is either exact + // or it doesn't exist (maps to extreme coordinates): + pnt2d_t xy_plus; + + // Convert image space to mosaic space. + xy_plus[0] = xy[0] + gridTransform->transform_.tile_min_[0]; + xy_plus[1] = xy[1] + gridTransform->transform_.tile_min_[1]; + uv = tile_to_mosaic->TransformPoint(xy_plus); + return uv[0] != std::numeric_limits::max(); + } + + // the local coordinate system, at the center of the tile: + pnt2d_t p00 = tile_min + (tile_max - tile_min) * 0.5; + pnt2d_t p10 = p00; + pnt2d_t p01 = p00; + + double x_unit = 1.0; + double y_unit = 1.0; + p10[0] += x_unit; + p01[1] += y_unit; + + vec2d_t x_axis = p10 - p00; + vec2d_t y_axis = p01 - p00; + + // FIXME: this may fail, because inverse may be unstable: + // the same coordinate system in the mosaic space: + pnt2d_t q00 = tile_to_mosaic->TransformPoint(p00); + pnt2d_t q10 = tile_to_mosaic->TransformPoint(p10); + pnt2d_t q01 = tile_to_mosaic->TransformPoint(p01); + + vec2d_t u_axis = q10 - q00; + vec2d_t v_axis = q01 - q00; + +#ifdef DEBUG_FIND_INVERSE + cerr << endl + << "x_axis: " << x_axis << endl + << "u_axis: " << u_axis << endl << endl + << "y_axis: " << y_axis << endl + << "v_axis: " << v_axis << endl << endl + << "p00: " << p00 << endl + << "q00: " << q00 << endl << endl; +#endif + + // initial guess: + uv = to_wcs(q00, + u_axis, + v_axis, + (xy[0] - p00[0]) / x_unit, + (xy[1] - p00[1]) / y_unit); + + // estimate the error: + const pnt2d_t & p0 = xy; + const pnt2d_t p1 = mosaic_to_tile->TransformPoint(uv); + vec2d_t er = p1 - p0; + +#ifdef DEBUG_FIND_INVERSE + cerr << "initial error: " << er.GetSquaredNorm() + << ", dx: " << er[0] << ", dy: " << er[1] + << endl; +#endif + + const double max_step_scale = 1.0; + double step_scale = max_step_scale; + unsigned int iteration = 0; + unsigned int successful_steps = 0; + + double e0_sqrd = er.GetSquaredNorm(); + + // don't try to improve samples with poor initialization (map to NaN): + if (e0_sqrd != e0_sqrd) return false; + + while (true) + { + const vec2d_t uv_correction = to_wcs(u_axis, + v_axis, + -er[0] / x_unit, + -er[1] / y_unit); + + pnt2d_t q = uv; + q[0] += step_scale * uv_correction[0]; + q[1] += step_scale * uv_correction[1]; + + pnt2d_t p = mosaic_to_tile->TransformPoint(q); + vec2d_t e = p - p0; + double e1_sqrd = e.GetSquaredNorm(); + +#ifdef DEBUG_FIND_INVERSE + cerr << setw(2) << iteration << ": " + << e0_sqrd << " vs " << e1_sqrd << " -- "; +#endif + if (e1_sqrd < e0_sqrd) + { +#ifdef DEBUG_FIND_INVERSE + cerr << "ok" << endl; +#endif + uv = q; + er = e; + e0_sqrd = e1_sqrd; + successful_steps++; + + if (successful_steps % pick_up_pace_steps == 0) + { + step_scale = std::min(max_step_scale, step_scale * 2.0); + +#ifdef DEBUG_FIND_INVERSE + cerr << successful_steps + << " successful steps -- increasing pace, new step length: " + << step_scale << endl; +#endif + } + } + else + { + step_scale /= 2.0; + successful_steps = 0; +#ifdef DEBUG_FIND_INVERSE + cerr << "relaxing and backtracking, new step length: " + << step_scale << endl; +#endif + } + + iteration++; + if (iteration >= max_iterations) break; + if (e0_sqrd < min_error_sqrd) break; + if (step_scale < min_step_scale) break; + } +#ifdef DEBUG_FIND_INVERSE + cerr << endl; +#endif + + return true; +} + +//---------------------------------------------------------------- +// find_inverse +// +// Given a forward transform, an approximate inverse transform, +// and image space coordinates, find mosaic space coordinates. +// +bool +find_inverse(const base_transform_t * mosaic_to_tile, + const base_transform_t * tile_to_mosaic, + const pnt2d_t & xy, // tile space + pnt2d_t & uv, // mosaic space + const unsigned int max_iterations, + const double min_step_scale, + const double min_error_sqrd, + const unsigned int pick_up_pace_steps) +{ + // #define DEBUG_FIND_INVERSE + if (tile_to_mosaic == NULL) + { + return false; + } + + if (dynamic_cast(tile_to_mosaic) != NULL || dynamic_cast(tile_to_mosaic) != NULL) + { + // special case for the grid transform -- the inverse is either exact + // or it doesn't exist (maps to extreme coordinates): + uv = tile_to_mosaic->TransformPoint(xy); + return uv[0] != std::numeric_limits::max(); + } + + // initial guess: + uv = tile_to_mosaic->TransformPoint(xy); + + // calculate initial error: + const pnt2d_t & p0 = xy; + pnt2d_t p1 = mosaic_to_tile->TransformPoint(uv); + vec2d_t er = p1 - p0; + +#ifdef DEBUG_FIND_INVERSE + cerr << "initial error: " << er.GetSquaredNorm() + << ", dx: " << er[0] << ", dy: " << er[1] + << endl; +#endif + + const double max_step_scale = 1.0; + double step_scale = max_step_scale; + unsigned int iteration = 0; + unsigned int successful_steps = 0; + + double e0_sqrd = er.GetSquaredNorm(); + + // don't try to improve samples with poor initialization (map to NaN): + if (e0_sqrd != e0_sqrd) return false; + + while (true) + { + pnt2d_t uv1 = tile_to_mosaic->TransformPoint(xy - er); + const vec2d_t uv_correction = uv - uv1; + + pnt2d_t q = uv; + q[0] += step_scale * uv_correction[0]; + q[1] += step_scale * uv_correction[1]; + + pnt2d_t p = mosaic_to_tile->TransformPoint(q); + vec2d_t e = p - p0; + double e1_sqrd = e.GetSquaredNorm(); + +#ifdef DEBUG_FIND_INVERSE + cerr << setw(2) << iteration << ": " + << e0_sqrd << " vs " << e1_sqrd << " -- "; +#endif + if (e1_sqrd < e0_sqrd) + { +#ifdef DEBUG_FIND_INVERSE + cerr << "ok" << endl; +#endif + uv = q; + er = e; + e0_sqrd = e1_sqrd; + successful_steps++; + + if (successful_steps % pick_up_pace_steps == 0) + { + step_scale = std::min(max_step_scale, step_scale * 2.0); + +#ifdef DEBUG_FIND_INVERSE + cerr << successful_steps + << " successful steps -- increasing pace, new step length: " + << step_scale << endl; +#endif + } + } + else + { + step_scale /= 2.0; + successful_steps = 0; +#ifdef DEBUG_FIND_INVERSE + cerr << "relaxing and backtracking, new step length: " + << step_scale << endl; +#endif + } + + iteration++; + if (iteration >= max_iterations) break; + if (e0_sqrd < min_error_sqrd) break; + if (step_scale < min_step_scale) break; + } +#ifdef DEBUG_FIND_INVERSE + cerr << endl; +#endif + + return true; +} + +//---------------------------------------------------------------- +// generate_landmarks_v1 +// +// Given a forward transform, generate a set of image space +// coordinates and find their matching mosaic space coordinates. +// +// version 1 -- uniform jittered sampling over the tile +// +bool +generate_landmarks_v1(const pnt2d_t & tile_min, + const pnt2d_t & tile_max, + const mask_t * tile_mask, + const base_transform_t * mosaic_to_tile, + const unsigned int samples, + std::vector & xy, + std::vector & uv, + bool refine) +{ + // #define DEBUG_LANDMARKS + + base_transform_t::Pointer tile_to_mosaic = mosaic_to_tile->GetInverse(); + if (tile_to_mosaic.GetPointer() == NULL) + { + return false; + } + + // the local coordinate system, at the center of the tile: + pnt2d_t p00 = tile_min + (tile_max - tile_min) * 0.5; + pnt2d_t p10 = p00; + pnt2d_t p01 = p00; + + double W = tile_max[0] - tile_min[0]; + double H = tile_max[1] - tile_min[1]; + double x_unit = 1.0; + double y_unit = 1.0; + p10[0] += x_unit; + p01[1] += y_unit; + + vec2d_t x_axis = p10 - p00; + vec2d_t y_axis = p01 - p00; + + // FIXME: this may fail, because tile_to_mosaic may be unstable: + // the same coordinate system in the mosaic space: + pnt2d_t q00 = tile_to_mosaic->TransformPoint(p00); + pnt2d_t q10 = tile_to_mosaic->TransformPoint(p10); + pnt2d_t q01 = tile_to_mosaic->TransformPoint(p01); + + vec2d_t u_axis = q10 - q00; + vec2d_t v_axis = q01 - q00; + +#ifdef DEBUG_LANDMARKS + cerr << endl + << "x_axis: " << x_axis << endl + << "u_axis: " << u_axis << endl << endl + << "y_axis: " << y_axis << endl + << "v_axis: " << v_axis << endl << endl + << "p00: " << p00 << endl + << "q00: " << q00 << endl << endl; +#endif + + std::list xy_list; + std::list uv_list; + std::list er_list; + + // sample the transform space: + for (unsigned int i = 0; i < samples; i++) + { + for (unsigned int j = 0; j < samples; j++) + { + const double s = (double(i) + drand()) / double(samples); + const double t = (double(j) + drand()) / double(samples); + + double x = tile_min[0] + s * W; + double y = tile_min[1] + t * H; + pnt2d_t pt_xy = pnt2d(x, y); + + // check to make sure this point is actually inside the tile/mask: + if (!pixel_in_mask(tile_mask, pt_xy)) continue; + + // map (approximately) into the mosaic space: + pnt2d_t pt_uv = to_wcs(q00, + u_axis, + v_axis, + (pt_xy[0] - p00[0]) / x_unit, + (pt_xy[1] - p00[1]) / y_unit); + + // estimate the error: + const pnt2d_t & p0 = pt_xy; + const pnt2d_t p1 = mosaic_to_tile->TransformPoint(pt_uv); + vec2d_t pt_er = p1 - p0; +#ifdef DEBUG_LANDMARKS + cerr << "initial error: " << pt_er.GetSquaredNorm() + << ", dx: " << p1[0] - p0[0] << ", dy: " << p1[1] - p0[1] + << endl; +#endif + + xy_list.push_back(pt_xy); + uv_list.push_back(pt_uv); + er_list.push_back(pt_er); + } + } + + // try to refine the boundary: + xy.assign(xy_list.begin(), xy_list.end()); + uv.assign(uv_list.begin(), uv_list.end()); + std::vector er(er_list.begin(), er_list.end()); + + if (refine) + { + for (unsigned int j = 0; j < xy.size(); j++) + { + const unsigned int max_iterations = 100; + const double min_step_scale = 1e-12; + const double min_error_sqrd = 1e-16; + const unsigned int pick_up_pace_steps = 5; + const double max_step_scale = 1.0; + + double step_scale = max_step_scale; + unsigned int iteration = 0; + unsigned int successful_steps = 0; + + const pnt2d_t & p0 = xy[j]; + double e0_sqrd = er[j].GetSquaredNorm(); + + // don't try to improve samples that fail: + if (e0_sqrd != e0_sqrd) continue; + + while (true) + { + const vec2d_t uv_correction = to_wcs(u_axis, + v_axis, + -er[j][0] / x_unit, + -er[j][1] / y_unit); + + pnt2d_t q = uv[j]; + q[0] += step_scale * uv_correction[0]; + q[1] += step_scale * uv_correction[1]; + + pnt2d_t p = mosaic_to_tile->TransformPoint(q); + vec2d_t e = p - p0; + double e1_sqrd = e.GetSquaredNorm(); + +#ifdef DEBUG_LANDMARKS + cerr << setw(3) << j << " " << setw(2) << iteration << ": " + << e0_sqrd << " vs " << e1_sqrd << " -- "; +#endif + if (e1_sqrd < e0_sqrd) + { +#ifdef DEBUG_LANDMARKS + cerr << "ok" << endl; +#endif + uv[j] = q; + er[j] = e; + e0_sqrd = e1_sqrd; + successful_steps++; + + if (successful_steps % pick_up_pace_steps == 0) + { + step_scale = std::min(max_step_scale, step_scale * 2.0); + +#ifdef DEBUG_LANDMARKS + cerr << successful_steps + << " successful steps -- increasing pace, new step length: " + << step_scale << endl; +#endif + } + } + else + { + step_scale /= 2.0; + successful_steps = 0; +#ifdef DEBUG_LANDMARKS + cerr << "relaxing and backtracking, new step length: " + << step_scale << endl; +#endif + } + + iteration++; + if (iteration >= max_iterations) break; + if (e0_sqrd < min_error_sqrd) break; + if (step_scale < min_step_scale) break; + } +#ifdef DEBUG_LANDMARKS + cerr << endl; +#endif + } + } + + // clean up -- remove failed samples: + xy_list.clear(); + uv_list.clear(); + for (unsigned int j = 0; j < er.size(); j++) + { + if (er[j] != er[j]) continue; + + xy_list.push_back(xy[j]); + uv_list.push_back(uv[j]); + } + + xy.assign(xy_list.begin(), xy_list.end()); + uv.assign(uv_list.begin(), uv_list.end()); + + return true; +} + + +//---------------------------------------------------------------- +// generate_landmarks_v2 +// +// Given a forward transform, generate a set of image space +// coordinates and find their matching mosaic space coordinates. +// +// version 2 -- non-uniform sampling skewed towards the center +// of the tile radially. This may be useful when +// the transform is ill-behaved at the tile edges. +// +bool +generate_landmarks_v2(const pnt2d_t & tile_min, + const pnt2d_t & tile_max, + const mask_t * tile_mask, + const base_transform_t * mosaic_to_tile, + const unsigned int samples, + std::vector & xy, + std::vector & uv, + bool refine) +{ + // #define DEBUG_LANDMARKS + + base_transform_t::Pointer tile_to_mosaic = mosaic_to_tile->GetInverse(); + if (tile_to_mosaic.GetPointer() == NULL) + { + return false; + } + + // the local coordinate system, at the center of the tile: + pnt2d_t p00 = tile_min + (tile_max - tile_min) * 0.5; + pnt2d_t p10 = p00; + pnt2d_t p01 = p00; + + double W = tile_max[0] - tile_min[0]; + double H = tile_max[1] - tile_min[1]; + double x_unit = 1.0; + double y_unit = 1.0; + p10[0] += x_unit; + p01[1] += y_unit; + + vec2d_t x_axis = p10 - p00; + vec2d_t y_axis = p01 - p00; + + // FIXME: this may fail, because tile_to_mosaic may be unstable: + // the same coordinate system in the mosaic space: + pnt2d_t q00 = tile_to_mosaic->TransformPoint(p00); + pnt2d_t q10 = tile_to_mosaic->TransformPoint(p10); + pnt2d_t q01 = tile_to_mosaic->TransformPoint(p01); + + vec2d_t u_axis = q10 - q00; + vec2d_t v_axis = q01 - q00; + +#ifdef DEBUG_LANDMARKS + cerr << endl + << "x_axis: " << x_axis << endl + << "u_axis: " << u_axis << endl << endl + << "y_axis: " << y_axis << endl + << "v_axis: " << v_axis << endl << endl + << "p00: " << p00 << endl + << "q00: " << q00 << endl << endl; +#endif + + std::list xy_list; + std::list uv_list; + std::list er_list; + + // sample the transform space (stay away from the corners): + for (unsigned int i = 0; i < samples; i++) + { + for (unsigned int j = 0; j < samples; j++) + { + // non-uniform sampling of the radius: + const double r = sqrt(drand()); + + // uniform sampling of the theta angle: + const double t = TWO_PI * (double(i) + drand()) / double(samples); + + // convert the polar coordinates into cartesian coordinates: + const double x = p00[0] + 0.7 * r * cos(t) * W / 2.0; + const double y = p00[1] + 0.7 * r * sin(t) * H / 2.0; + + pnt2d_t pt_xy = pnt2d(x, y); + + // check to make sure this point is actually inside the tile/mask: + if (!pixel_in_mask(tile_mask, pt_xy)) continue; + + // map (approximately) into the mosaic space: + pnt2d_t pt_uv = to_wcs(q00, + u_axis, + v_axis, + (pt_xy[0] - p00[0]) / x_unit, + (pt_xy[1] - p00[1]) / y_unit); + + if (refine) + { + // estimate the error: + const pnt2d_t & p0 = pt_xy; + const pnt2d_t p1 = mosaic_to_tile->TransformPoint(pt_uv); + vec2d_t pt_er = p1 - p0; +#ifdef DEBUG_LANDMARKS + cerr << "initial error: " << pt_er.GetSquaredNorm() + << ", dx: " << p1[0] - p0[0] << ", dy: " << p1[1] - p0[1] + << endl; +#endif + + xy_list.push_back(pt_xy); + uv_list.push_back(pt_uv); + er_list.push_back(pt_er); + } + else + { + pt_xy = mosaic_to_tile->TransformPoint(pt_uv); + xy_list.push_back(pt_xy); + uv_list.push_back(pt_uv); + } + } + } + + // try to refine the boundary: + xy.assign(xy_list.begin(), xy_list.end()); + uv.assign(uv_list.begin(), uv_list.end()); + + if (refine) + { + std::vector er(er_list.begin(), er_list.end()); + + for (unsigned int j = 0; j < xy.size(); j++) + { + const unsigned int max_iterations = 100; + const double min_step_scale = 1e-12; + const double min_error_sqrd = 1e-16; + const unsigned int pick_up_pace_steps = 5; + const double max_step_scale = 1.0; + + double step_scale = max_step_scale; + unsigned int iteration = 0; + unsigned int successful_steps = 0; + + const pnt2d_t & p0 = xy[j]; + double e0_sqrd = er[j].GetSquaredNorm(); + + while (true) + { + const vec2d_t uv_correction = to_wcs(u_axis, + v_axis, + -er[j][0] / x_unit, + -er[j][1] / y_unit); + + pnt2d_t q = uv[j]; + q[0] += step_scale * uv_correction[0]; + q[1] += step_scale * uv_correction[1]; + + pnt2d_t p = mosaic_to_tile->TransformPoint(q); + vec2d_t e = p - p0; + double e1_sqrd = e.GetSquaredNorm(); + +#ifdef DEBUG_LANDMARKS + cerr << setw(3) << j << " " << setw(2) << iteration << ": " + << e0_sqrd << " vs " << e1_sqrd << " -- "; +#endif + if (e1_sqrd < e0_sqrd) + { +#ifdef DEBUG_LANDMARKS + cerr << "ok" << endl; +#endif + uv[j] = q; + er[j] = e; + e0_sqrd = e1_sqrd; + successful_steps++; + + if (successful_steps % pick_up_pace_steps == 0) + { + step_scale = std::min(max_step_scale, step_scale * 2.0); + +#ifdef DEBUG_LANDMARKS + cerr << successful_steps + << " successful steps -- increasing pace, new step length: " + << step_scale << endl; +#endif + } + } + else + { + step_scale /= 2.0; + successful_steps = 0; +#ifdef DEBUG_LANDMARKS + cerr << "relaxing and backtracking, new step length: " + << step_scale << endl; +#endif + } + + iteration++; + if (iteration >= max_iterations) break; + if (e0_sqrd < min_error_sqrd) break; + if (step_scale < min_step_scale) break; + } +#ifdef DEBUG_LANDMARKS + cerr << endl; +#endif + } + } + + return true; +} + +//---------------------------------------------------------------- +// generate_landmarks +// +// Given a forward transform, generate a set of image space +// coordinates and find their matching mosaic space coordinates. +// +// version 1 -- uniform jittered sampling over the tile +// version 2 -- non-uniform sampling skewed towards the center +// of the tile radially. This may be useful when +// the transform is ill-behaved at the tile edges. +// +bool +generate_landmarks(const pnt2d_t & tile_min, + const pnt2d_t & tile_max, + const mask_t * tile_mask, + const base_transform_t * mosaic_to_tile, + const unsigned int samples, + std::vector & xy, + std::vector & uv, + const unsigned int version, + const bool refine) +{ + switch (version) + { + case 1: + return generate_landmarks_v1(tile_min, + tile_max, + tile_mask, + mosaic_to_tile, + samples, + xy, + uv, + refine); + + case 2: + return generate_landmarks_v2(tile_min, + tile_max, + tile_mask, + mosaic_to_tile, + samples, + xy, + uv, + refine); + } + + return false; +} + +//---------------------------------------------------------------- +// generate_landmarks +// +// Given a forward transform, generate a set of image space +// coordinates and find their matching mosaic space coordinates. +// +// version 1 -- uniform jittered sampling over the tile +// version 2 -- non-uniform sampling skewed towards the center +// of the tile radially. This may be useful when +// the transform is ill-behaved at the tile edges. +// +bool +generate_landmarks(const image_t * tile, + const mask_t * mask, + const base_transform_t * mosaic_to_tile, + const unsigned int samples, + std::vector & xy, + std::vector & uv, + const unsigned int version, + const bool refine) +{ + image_t::SpacingType sp = tile->GetSpacing(); + image_t::SizeType sz = tile->GetLargestPossibleRegion().GetSize(); + const pnt2d_t min = tile->GetOrigin(); + const pnt2d_t max = min + vec2d(sp[0] * double(sz[0]), + sp[1] * double(sz[1])); + + switch (version) + { + case 1: + return generate_landmarks_v1(min, + max, + mask, + mosaic_to_tile, + samples, + xy, + uv, + refine); + + case 2: + return generate_landmarks_v2(min, + max, + mask, + mosaic_to_tile, + samples, + xy, + uv, + refine); + } + + return false; +} + +//---------------------------------------------------------------- +// calc_size +// +// Given a bounding box expressed in the image space and +// pixel spacing, calculate corresponding image size. +// +image_t::SizeType +calc_size(const pnt2d_t & min, + const pnt2d_t & max, + const double & spacing) +{ + image_t::SizeType sz; + sz[0] = (unsigned int)((ceilf(max[0]) - floor(min[0])) / spacing); + sz[1] = (unsigned int)((ceilf(max[1]) - floor(min[1])) / spacing); + return sz; +} + + +the_thread_pool_t * save_mosaic::_pthread_pool; + + +//---------------------------------------------------------------- +// hsv_to_rgb +// +// Convenience function for converting between HSV/RGB color +// spaces. This is used for colormapping. +// +xyz_t + hsv_to_rgb(const xyz_t & HSV) +{ + double H = HSV[0]; + double S = HSV[1]; + double V = HSV[2]; + + xyz_t RGB; + double & R = RGB[0]; + double & G = RGB[1]; + double & B = RGB[2]; + + if (S == 0.0) + { + // monochromatic: + R = V; + G = V; + B = V; + return RGB; + } + + H *= 6.0; + double i = floor(H); + double f = H - i; + + double p = V * (1.0 - S); + double q = V * (1.0 - S * f); + double t = V * (1.0 - S * (1 - f)); + + if (i == 0.0) + { + R = V; + G = t; + B = p; + } + else if (i == 1.0) + { + R = q; + G = V; + B = p; + } + else if (i == 2.0) + { + R = p; + G = V; + B = t; + } + else if (i == 3.0) + { + R = p; + G = q; + B = V; + } + else if (i == 4.0) + { + R = t; + G = p; + B = V; + } + else + { + // i == 5.0 + R = V; + G = p; + B = q; + } + + return RGB; +} + +//---------------------------------------------------------------- +// rgb_to_hsv +// +// Convenience function for converting between RGB/HSV color +// spaces. This is used for colormapping. +// +xyz_t + rgb_to_hsv(const xyz_t & RGB) +{ + double R = RGB[0]; + double G = RGB[1]; + double B = RGB[2]; + + xyz_t HSV; + double & H = HSV[0]; + double & S = HSV[1]; + double & V = HSV[2]; + + double min = std::min(R, std::min(G, B)); + double max = std::max(R, std::max(G, B)); + V = max; + + double delta = max - min; + if (max == 0) + { + S = 0; + H = -1; + } + else + { + S = delta / max; + + if (delta == 0) + { + delta = 1; + } + + if (R == max) + { + // between yellow & magenta + H = (G - B) / delta; + } + else if (G == max) + { + // between cyan & yellow + H = (B - R) / delta + 2; + } + else + { + // between magenta & cyan + H = (R - G) / delta + 4; + } + + H /= 6.0; + + if (H < 0.0) + { + H = H + 1.0; + } + } + + return HSV; +} diff --git a/src/itkIRText.cxx b/src/itkIRText.cxx new file mode 100644 index 0000000..3c9479a --- /dev/null +++ b/src/itkIRText.cxx @@ -0,0 +1,580 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_text.hxx +// Author : Pavel A. Koshevoy +// Created : Sun Aug 29 14:53:00 MDT 2004 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : ASCII text convenience class. + +// local includes: +#include "itkIRText.h" + +// system includes: +#include +#include +#include +#include +#include +#include +#include +#include + + +//---------------------------------------------------------------- +// the_text_t::the_text_t +// +the_text_t::the_text_t(const char * text): + text_(NULL), + size_(0) +{ + assign(text); +} + +//---------------------------------------------------------------- +// the_text_t::the_text_t +// +the_text_t::the_text_t(const char * text, const size_t & size): + text_(NULL), + size_(0) +{ + assign(text, size); +} + +//---------------------------------------------------------------- +// the_text_t::the_text_t +// +the_text_t::the_text_t(const the_text_t & text): + text_(NULL), + size_(0) +{ + assign(text.text_, text.size_); +} + +//---------------------------------------------------------------- +// the_text_t::the_text_t +// +the_text_t::the_text_t(const std::list & text): + text_(NULL), + size_(text.size()) +{ + text_ = new char [size_ + 1]; + size_t j = 0; + for (std::list::const_iterator i = text.begin(); + i != text.end(); ++i, ++j) + { + text_[j] = *i; + } +} + +//---------------------------------------------------------------- +// the_text_t::~the_text_t +// +the_text_t::~the_text_t() +{ + clear(); +} + +//---------------------------------------------------------------- +// the_text_t::assign +// +void +the_text_t::assign(const char * text, const size_t & text_size) +{ + char * text_copy = new char [text_size + 1]; + for (size_t i = 0; i < text_size; i++) + { + text_copy[i] = text[i]; + } + + text_copy[text_size] = '\0'; + + delete [] text_; + text_ = text_copy; + size_ = text_size; +} + +//---------------------------------------------------------------- +// the_text_t::append +// +void +the_text_t::append(const char * text, const size_t & text_size) +{ + char * text_copy = new char [size_ + text_size + 1]; + + for (size_t i = 0; i < size_; i++) + { + text_copy[i] = text_[i]; + } + + for (size_t i = 0; i < text_size; i++) + { + text_copy[i + size_] = text[i]; + } + + text_copy[size_ + text_size] = '\0'; + + delete [] text_; + text_ = text_copy; + size_ += text_size; +} + +//---------------------------------------------------------------- +// the_text_t::replace +// replaces all occurances of one character with another. +void +the_text_t::replace(const char find, const char replace) +{ + char * cursor = text_; + while (*cursor != '\0') + { + if ( *cursor == find ) + *cursor = replace; + cursor++; + } +} + +//---------------------------------------------------------------- +// the_text_t::toShort +// +short int +the_text_t::toShort(bool * ok, int base) const +{ + long int num = toULong(ok, base); + if (ok != NULL) *ok &= ((num >= SHRT_MIN) && (num <= SHRT_MAX)); + return (short int)(num); +} + +//---------------------------------------------------------------- +// the_text_t::toUShort +// +unsigned short int +the_text_t::toUShort(bool * ok, int base) const +{ + unsigned long int num = toULong(ok, base); + if (ok != NULL) *ok &= (num <= USHRT_MAX); + return (unsigned short int)(num); +} + +//---------------------------------------------------------------- +// the_text_t::toInt +// +int +the_text_t::toInt(bool * ok, int base) const +{ + long int num = toULong(ok, base); + if (ok != NULL) *ok &= ((num >= INT_MIN) && (num <= INT_MAX)); + return int(num); +} + +//---------------------------------------------------------------- +// the_text_t::toUInt +// +unsigned int +the_text_t::toUInt(bool * ok, int base) const +{ + unsigned long int num = toULong(ok, base); + if (ok != NULL) *ok &= (num <= UINT_MAX); + return (unsigned int)(num); +} + +//---------------------------------------------------------------- +// the_text_t::toLong +// +long int +the_text_t::toLong(bool * ok, int base) const +{ + char * endptr = NULL; + long int num = strtol(text_, &endptr, base); + if (ok != NULL) *ok = !(text_ == endptr || errno == ERANGE); + return num; +} + +//---------------------------------------------------------------- +// the_text_t::toULong +// +unsigned long int +the_text_t::toULong(bool * ok, int base) const +{ + char * endptr = NULL; + unsigned long int num = strtoul(text_, &endptr, base); + if (ok != NULL) *ok = !(text_ == endptr || errno == ERANGE); + return num; +} + +//---------------------------------------------------------------- +// the_text_t::toFloat +// +float +the_text_t::toFloat(bool * ok) const +{ + double num = toDouble(ok); + return float(num); +} + +//---------------------------------------------------------------- +// the_text_t::toDouble +// +double +the_text_t::toDouble(bool * ok) const +{ + char * endptr = NULL; + double num = strtod(text_, &endptr); + if (ok != NULL) *ok = !(text_ == endptr || errno == ERANGE); + return num; +} + +//---------------------------------------------------------------- +// the_text_t::to_ascii +// +void +the_text_t::to_ascii() +{ + for (size_t i = 0; i < size_; i++) + { + text_[i] = ((size_t)(text_[i]) < 128) ? text_[i] : '?'; + } +} + +//---------------------------------------------------------------- +// the_text_t::to_lower +// +void +the_text_t::to_lower() +{ + for (size_t i = 0; i < size_; i++) + { + text_[i] = tolower(text_[i]); + } +} + +//---------------------------------------------------------------- +// the_text_t::to_upper +// +void +the_text_t::to_upper() +{ + for (size_t i = 0; i < size_; i++) + { + text_[i] = toupper(text_[i]); + } +} + +//---------------------------------------------------------------- +// the_text_t::fill +// +void +the_text_t::fill(const char & c, const size_t size) +{ + char * text = new char [size + 1]; + for (size_t i = 0; i < size; i++) + { + text[i] = c; + } + text[size] = '\0'; + + delete [] text_; + text_ = text; + size_ = size; +} + +//---------------------------------------------------------------- +// the_text_t::match_head +// +bool +the_text_t::match_head(const the_text_t & t, bool ignore_case) const +{ + if (t.size_ > size_) return false; + return match_text(t, 0, ignore_case); +} + +//---------------------------------------------------------------- +// the_text_t::match_tail +// +bool +the_text_t::match_tail(const the_text_t & t, bool ignore_case) const +{ + if (t.size_ > size_) return false; + return match_text(t, size_ - t.size_, ignore_case); +} + +//---------------------------------------------------------------- +// the_text_t::match_text +// +bool +the_text_t::match_text(const the_text_t & t, + const size_t & start, + bool ignore_case) const +{ + size_t end = start + t.size_; + if (end > size_) return false; + + for (size_t i = start; i < end; i++) + { + char a = text_[i]; + char b = t.text_[i - start]; + + if (ignore_case) + { + a = tolower(a); + b = tolower(b); + } + + if (a != b) return false; + } + + return true; +} + +//---------------------------------------------------------------- +// the_text_t::simplify_ws +// +the_text_t +the_text_t::simplify_ws() const +{ + // find the first non-whitespace character: + int start = 0; + for (; start < int(size_) && isspace(text_[start]); start++) + { + // skipping whitespace... + } + + // NOTE: an all-whitespace string will simplify to an empty string: + if (start == int(size_)) + { + return the_text_t(""); + } + + // find the last non-whitespace character: + int finish = int(size_) - 1; + for (; finish >= start && isspace(text_[finish]); finish--); + + // intermediate storage: + the_text_t tmp; + tmp.fill('\0', size_t((finish + 1) - start)); + + size_t j = 0; + bool prev_ws = false; + for (int i = start; i <= finish; i++) + { + char c = isspace(text_[i]) ? ' ' : text_[i]; + if (c == ' ' && prev_ws) continue; + + prev_ws = (c == ' '); + tmp.text_[j] = c; + j++; + } + + the_text_t out(tmp.text(), j); + return out; +} + +//---------------------------------------------------------------- +// the_text_t::split +// +size_t +the_text_t::split(std::vector & tokens, + const char & separator, + const bool & empty_ok) const +{ + if (size_ == 0) + { + tokens.resize(0); + return 0; + } + + // find the separators: + typedef std::list list_t; + list_t separators; + for (size_t i = 0; i < size_; i++) + { + if (text_[i] != separator) continue; + separators.push_back(i); + } + separators.push_back(size_); + + std::list tmp; + + typedef std::list::iterator iter_t; + size_t a = 0; + for (iter_t i = separators.begin(); i != separators.end(); ++i) + { + size_t b = *i; + + if (b - a == 0 && empty_ok) + { + tmp.push_back(the_text_t("")); + } + else if (b - a > 0) + { + tmp.push_back(the_text_t(&text_[a], b - a)); + } + + a = b + 1; + } + + tokens.resize(tmp.size()); + a = 0; + for (std::list::iterator i = tmp.begin(); i != tmp.end(); ++i) + { + tokens[a] = *i; + a++; + } + + return tmp.size(); +} + +//---------------------------------------------------------------- +// the_text_t::splitAt +// Splits the string around the n'th occurance of split_char +// if n > size_, it is split around the last occurance. +std::vector the_text_t::splitAt(const char split_char, + unsigned int n) const +{ + char * cursor = text_; + + // Find the desired character. + unsigned int position = 0; + unsigned int found_pos = std::numeric_limits::max(); + while (*cursor != '\0') + { + if ( *cursor == split_char ) + { + n--; + found_pos = position; + if ( n == 0 ) + break; + } + cursor++; + position++; + } + + std::vector str_parts; + str_parts.resize(2); + if ( found_pos != std::numeric_limits::max() ) + { + // Split point found, do split. + the_text_t left(&text_[0], found_pos); + the_text_t right(&text_[found_pos + 1], size_ - found_pos - 1); + + str_parts[0] = left; + str_parts[1] = right; + } + else + { + // No split point found. + str_parts[0] = ""; + str_parts[1] = the_text_t( text_ ); + } + + return str_parts; +} + +//---------------------------------------------------------------- +// the_text_t::contains +// +size_t +the_text_t::contains(const char & symbol) const +{ + size_t count = 0; + for (size_t i = 0; i < size_; i++) + { + if (text_[i] != symbol) continue; + count++; + } + + return count; +} + + +//---------------------------------------------------------------- +// operator << +// +std::ostream & +operator << (std::ostream & out, const the_text_t & text) +{ + std::string tmp(text.text(), text.size()); + out << tmp; + return out; +} + +//---------------------------------------------------------------- +// operator >> +// +std::istream & +operator >> (std::istream & in, the_text_t & text) +{ + std::string tmp; + in >> tmp; + text.assign(tmp.data(), tmp.size()); + return in; +} + +//---------------------------------------------------------------- +// getline +// +std::istream & +getline(std::istream & in, the_text_t & text) +{ + std::string tmp; + getline(in, tmp); + if (!tmp.empty()) + { + std::size_t len = tmp.size(); + if (tmp[len - 1] == '\r') + { + // truncate the \r character + tmp.resize(len - 1); + } + } + + text.assign(tmp.data(), tmp.size()); + return in; +} + +//---------------------------------------------------------------- +// to_binary +// +the_text_t +to_binary(const unsigned char & byte, bool lsb_first) +{ + the_text_t str; + str.fill('0', 8); + + unsigned char mask = 1; + if (lsb_first) + { + for (int i = 0; i < 8; i++, mask *= 2) + { + str[i] = (byte & mask) ? '1' : '0'; + } + } + else + { + for (int i = 7; i > -1; i--, mask *= 2) + { + str[i] = (byte & mask) ? '1' : '0'; + } + } + + return str; +} diff --git a/src/itkIRUtils.cxx b/src/itkIRUtils.cxx new file mode 100644 index 0000000..db4e4ab --- /dev/null +++ b/src/itkIRUtils.cxx @@ -0,0 +1,314 @@ +// -*- Mode: c++; tab-width: 8; c-basic-offset: 2; indent-tabs-mode: t -*- +// NOTE: the first line of this file sets up source code indentation rules +// for Emacs; it is also a hint to anyone modifying this file. + +/* + 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 2 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, write to the Free Software + Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA +*/ + + +// File : the_utils.cxx +// Author : Pavel A. Koshevoy +// Created : Tue Nov 4 20:32:04 MST 2008 +// Copyright : (C) 2004-2008 University of Utah +// License : GPLv2 +// Description : utility functions + +// system includes: +#ifndef _USE_MATH_DEFINES +#define _USE_MATH_DEFINES +#endif + +// system includes: +#ifdef _WIN32 +#ifndef NOMINMAX +#define NOMINMAX +#endif +#include +#include +#else +#include +#include +#include +#include +#endif +#include + +// local includes: +#include "itkIRUtils.h" + + +//---------------------------------------------------------------- +// sleep_msec +// +void +sleep_msec(size_t msec) +{ +#ifdef WIN32 + Sleep((DWORD)(msec)); +#else + usleep(msec * 1000); +#endif +} + +//---------------------------------------------------------------- +// restore_console_stdio +// +bool +restore_console_stdio() +{ +#ifdef _WIN32 + AllocConsole(); + +#pragma warning(push) +#pragma warning(disable: 4996) + + freopen("conin$", "r", stdin); + freopen("conout$", "w", stdout); + freopen("conout$", "w", stderr); + +#pragma warning(pop) + + HANDLE std_out_handle = GetStdHandle(STD_OUTPUT_HANDLE); + if (std_out_handle == INVALID_HANDLE_VALUE) + { + return false; + } + + COORD console_buffer_size; + console_buffer_size.X = 80; + console_buffer_size.Y = 9999; + SetConsoleScreenBufferSize(std_out_handle, + console_buffer_size); +#endif + + return true; +} + + +namespace the +{ +#ifdef _WIN32 + //---------------------------------------------------------------- + // utf8_to_utf16 + // + static void + utf8_to_utf16(const char * utf8, wchar_t *& utf16) + { + int wcs_size = + MultiByteToWideChar(CP_UTF8, // encoding (ansi, utf, etc...) + 0, // flags (precomposed, composite,... ) + utf8, // source multi-byte character string + -1, // number of bytes in the source string + NULL, // wide-character destination + 0); // destination buffer size + + utf16 = new wchar_t[wcs_size + 1]; + MultiByteToWideChar(CP_UTF8, + 0, + utf8, + -1, + utf16, + wcs_size); + } +#endif + + //---------------------------------------------------------------- + // open_utf8 + // + int + open_utf8(const char * filename_utf8, int oflag, int pmode) + { + int fd = -1; + +#ifdef _WIN32 + // on windows utf-8 has to be converted to utf-16 + wchar_t * filename_utf16 = 0; + utf8_to_utf16(filename_utf8, filename_utf16); + + int sflag = _SH_DENYNO; + _wsopen_s(&fd, filename_utf16, oflag, sflag, pmode); + delete [] filename_utf16; + +#else + // assume utf-8 is supported natively: + fd = open(filename_utf8, oflag, pmode); +#endif + + return fd; + } + + //---------------------------------------------------------------- + // open_utf8 + // + void + open_utf8(std::fstream & fstream_to_open, + const char * filename_utf8, + std::ios_base::openmode mode) + { +#ifdef _WIN32 + // on windows utf-8 has to be converted to utf-16 + wchar_t * filename_utf16 = 0; + utf8_to_utf16(filename_utf8, filename_utf16); + + fstream_to_open.open(filename_utf16, mode); + delete [] filename_utf16; + +#else + // assume utf-8 is supported natively: + fstream_to_open.open(filename_utf8, mode); +#endif + } + + //---------------------------------------------------------------- + // fopen_utf8 + // + FILE * + fopen_utf8(const char * filename_utf8, const char * mode) + { + FILE * file = NULL; + +#ifdef _WIN32 + wchar_t * filename_utf16 = NULL; + utf8_to_utf16(filename_utf8, filename_utf16); + + wchar_t * mode_utf16 = NULL; + utf8_to_utf16(mode, mode_utf16); + + _wfopen_s(&file, filename_utf16, mode_utf16); + delete [] filename_utf16; + delete [] mode_utf16; +#else + file = fopen(filename_utf8, mode); +#endif + + return file; + } + + //---------------------------------------------------------------- + // rename_utf8 + // + int + rename_utf8(const char * old_utf8, const char * new_utf8) + { +#ifdef _WIN32 + wchar_t * old_utf16 = NULL; + utf8_to_utf16(old_utf8, old_utf16); + + wchar_t * new_utf16 = NULL; + utf8_to_utf16(new_utf8, new_utf16); + + int ret = _wrename(old_utf16, new_utf16); + + delete [] old_utf16; + delete [] new_utf16; +#else + + int ret = rename(old_utf8, new_utf8); +#endif + + return ret; + } + + //---------------------------------------------------------------- + // remove_utf8 + // + int + remove_utf8(const char * filename_utf8) + { +#ifdef _WIN32 + wchar_t * filename_utf16 = NULL; + utf8_to_utf16(filename_utf8, filename_utf16); + + int ret = _wremove(filename_utf16); + delete [] filename_utf16; +#else + + int ret = remove(filename_utf8); +#endif + + return ret; + } + + //---------------------------------------------------------------- + // rmdir_utf8 + // + int + rmdir_utf8(const char * dir_utf8) + { +#ifdef _WIN32 + wchar_t * dir_utf16 = NULL; + utf8_to_utf16(dir_utf8, dir_utf16); + + int ret = _wrmdir(dir_utf16); + delete [] dir_utf16; +#else + + int ret = remove(dir_utf8); +#endif + + return ret; + } + + //---------------------------------------------------------------- + // mkdir_utf8 + // + int + mkdir_utf8(const char * path_utf8) + { +#ifdef _WIN32 + wchar_t * path_utf16 = NULL; + utf8_to_utf16(path_utf8, path_utf16); + + int ret = _wmkdir(path_utf16); + delete [] path_utf16; +#else + + int ret = mkdir(path_utf8, S_IRWXU); +#endif + + return ret; + } + + //---------------------------------------------------------------- + // fseek64 + // + int + fseek64(FILE * file, off_t offset, int whence) + { +#ifdef _WIN32 + int ret = _fseeki64(file, offset, whence); +#else + int ret = fseeko(file, offset, whence); +#endif + + return ret; + } + + //---------------------------------------------------------------- + // ftell64 + // + off_t + ftell64(const FILE * file) + { +#ifdef _WIN32 + off_t pos = _ftelli64(const_cast(file)); +#else + off_t pos = ftello(const_cast(file)); +#endif + + return pos; + } +}