diff --git a/Code/Vahine/CMakeLists.txt b/Code/Vahine/CMakeLists.txt
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/Code/Vahine/COPYING b/Code/Vahine/COPYING
deleted file mode 100644
index 9912d3e97c69cc78f5ea74ccb83d861818b80263..0000000000000000000000000000000000000000
--- a/Code/Vahine/COPYING
+++ /dev/null
@@ -1,491 +0,0 @@
-CeCILL FREE SOFTWARE LICENSE AGREEMENT
-
-Notice
-
-This Agreement is a Free Software license agreement that is the result
-of discussions between its authors in order to ensure compliance with
-the two main principles guiding its drafting :
-
-* firstly, compliance with the principles governing the distribution
-of Free Software : access to source code, broad rights granted to
- users,
-* secondly, the election of a governing law, French law, with which
-it is conformant, both as regards the law of torts and
-intellectual property law, and the protection that it offers to
-both authors and holders of the economic rights over software.
-
-The authors of the CeCILL (for Ce[a] C[nrs] I[nria] L[ogiciel] L[ibre])
-license are :
-
- Commissariat l
- 'Energie Atomique - CEA, a public scientific, technical
-and industrial research establishment, having its principal place of
-business at 25 rue Leblanc, immeuble Le Ponant D, 75015 Paris, France.
-
-Centre National de la Recherche Scientifique - CNRS, a public scientific
-and technological establishment, having its principal place of business
-at 3 rue Michel-Ange, 75794 Paris cedex 16, France.
-
-Institut National de Recherche en Informatique et en Automatique -
-INRIA, a public scientific and technological establishment, having its
-principal place of business at Domaine de Voluceau, Rocquencourt, BP
-105, 78153 Le Chesnay cedex, France.
-
-
- Preamble
-
-The purpose of this Free Software license agreement is to grant users
-the right to modify and redistribute the software governed by this
-license within the framework of an open source distribution model.
-
-The exercising of these rights is conditional upon certain obligations
-for users so as to preserve this status for all subsequent redistributions.
-
-In consideration of access to the source code and the rights to copy,
-modify and redistribute granted by the license, users are provided only
-with a limited warranty and the software's
- author, the holder of the
- economic rights, and the successive licensors only have limited liability.
-
- In this respect, the risks associated with loading, using, modifying
- and/or developing or reproducing the software by the user are brought to
- the user
- 's attention, given its Free Software status, which may make it
-complicated to use, with the result that its use is reserved for
-developers and experienced professionals having in-depth computer
-knowledge. Users are therefore encouraged to load and test the
-suitability of the software as regards their requirements in conditions
-enabling the security of their systems and/or data to be ensured and,
-more generally, to use and operate it in the same conditions of
-security. This Agreement may be freely reproduced and published,
-provided it is not altered, and that no provisions are either added or
-removed herefrom.
-
-This Agreement may apply to any or all software for which the holder of
-the economic rights decides to submit the use thereof to its provisions.
-
-
- Article 1 - DEFINITIONS
-
-For the purpose of this Agreement, when the following expressions
-commence with a capital letter, they shall have the following meaning:
-
-Agreement: means this license agreement, and its possible subsequent
-versions and annexes.
-
-Software: means the software in its Object Code and/or Source Code form
-and, where applicable, its documentation, "as is" when the Licensee
-accepts the Agreement.
-
-Initial Software: means the Software in its Source Code and possibly its
-Object Code form and, where applicable, its documentation, "as is" when
-it is first distributed under the terms and conditions of the Agreement.
-
-Modified Software: means the Software modified by at least one
-Contribution.
-
-Source Code: means all the Software's
- instructions and program lines to
- which access is required so as to modify the Software.
-
- Object Code : means the binary files originating from the compilation of
- the Source Code.
-
- Holder : means the holder(s) of the economic rights over the Initial
- Software.
-
- Licensee : means the Software user(s) having accepted the Agreement.
-
- Contributor : means a Licensee having made at least one Contribution.
-
- Licensor : means the Holder, or any other individual or legal entity, who
- distributes the Software under the Agreement.
-
- Contribution : means any or all modifications, corrections, translations,
- adaptations and/or new functions integrated into the Software by any or
- all Contributors, as well as any or all Internal Modules.
-
- Module : means a set of sources files including their documentation that
- enables supplementary functions or services in addition to those offered
- by the Software.
-
- External Module : means any or all Modules, not derived from the
- Software, so that this Module and the Software run in separate address
- spaces, with one calling the other when they are run.
-
- Internal Module : means any or all Module, connected to the Software so
- that they both execute in the same address space.
-
- GNU GPL : means the GNU General Public License version 2 or any
- subsequent version, as published by the Free Software Foundation Inc.
-
- Parties : mean both the Licensee and the Licensor.
-
- These expressions may be used both in singular and plural form.
-
- Article 2 - PURPOSE
-
- The purpose of the Agreement is the grant by the Licensor to the
- Licensee of a non-exclusive, transferable and worldwide license for the
- Software as set forth in Article 5 hereinafter for the whole term of the
- protection granted by the rights over said Software.
-
- Article 3 - ACCEPTANCE
-
- 3.1 The Licensee shall be deemed as having accepted the terms and
- conditions of this Agreement upon the occurrence of the first of the
- following events :
-
- *(i) loading the Software by any or all means, notably, by
- downloading from a remote server, or by loading from a physical
- medium;
-*(ii) the first time the Licensee exercises any of the rights
-granted hereunder.
-
-3.2 One copy of the Agreement, containing a notice relating to the
-characteristics of the Software, to the limited warranty, and to the
-fact that its use is restricted to experienced users has been provided
-to the Licensee prior to its acceptance as set forth in Article 3.1
-hereinabove, and the Licensee hereby acknowledges that it has read and
-understood it.
-
-Article 4 - EFFECTIVE DATE AND TERM
-
-4.1 EFFECTIVE DATE
-
-The Agreement shall become effective on the date when it is accepted by
-the Licensee as set forth in Article 3.1.
-
-4.2 TERM
-
-The Agreement shall remain in force for the entire legal term of
-protection of the economic rights over the Software.
-
-Article 5 - SCOPE OF RIGHTS GRANTED
-
-The Licensor hereby grants to the Licensee, who accepts, the following
-rights over the Software for any or all use, and for the term of the
-Agreement, on the basis of the terms and conditions set forth hereinafter.
-
-Besides, if the Licensor owns or comes to own one or more patents
-protecting all or part of the functions of the Software or of its
-components, the Licensor undertakes not to enforce the rights granted by
-these patents against successive Licensees using, exploiting or
- modifying the Software.If these patents are transferred, the Licensor
- undertakes to have the transferees subscribe to the obligations set
- forth in this paragraph.
-
- 5.1 RIGHT OF USE
-
- The Licensee is authorized to use the Software, without any limitation
- as to its fields of application, with it being hereinafter specified
- that this comprises :
-
- 1. permanent or temporary reproduction of all or part of the Software
- by any or all means and in any or all form.
-
- 2. loading, displaying, running, or storing the Software on any or
- all medium.
-
- 3. entitlement to observe, study or test its operation so as to
- determine the ideas and principles behind any or all constituent
- elements of said Software.This shall apply when the Licensee
- carries out any or all loading, displaying, running, transmission
- or storage operation as regards the Software, that it is entitled
- to carry out hereunder.
-
- 5.2 ENTITLEMENT TO MAKE CONTRIBUTIONS
-
- The right to make Contributions includes the right to translate, adapt,
- arrange, or make any or all modifications to the Software, and the right
- to reproduce the resulting software.
-
- The Licensee is authorized to make any or all Contributions to the
- Software provided that it includes an explicit notice that it is the
- author of said Contribution and indicates the date of the creation thereof.
-
- 5.3 RIGHT OF DISTRIBUTION
-
- In particular, the right of distribution includes the right to publish,
- transmit and communicate the Software to the general public on any or
- all medium, and by any or all means, and the right to market, either in
- consideration of a fee, or free of charge, one or more copies of the
- Software by any means.
-
- The Licensee is further authorized to distribute copies of the modified
- or unmodified Software to third parties according to the terms and
- conditions set forth hereinafter.
-
- 5.3 .1 DISTRIBUTION OF SOFTWARE WITHOUT MODIFICATION
-
- The Licensee is authorized to distribute true copies of the Software in
- Source Code or Object Code form, provided that said distribution
- complies with all the provisions of the Agreement and is accompanied by :
-
- 1. a copy of the Agreement,
-
- 2. a notice relating to the limitation of both the Licensor
- 's
- warranty and liability as set forth in Articles 8 and 9,
-
-and that, in the event that only the Object Code of the Software is
-redistributed, the Licensee allows future Licensees unhindered access to
-the full Source Code of the Software by indicating how to access it, it
-being understood that the additional cost of acquiring the Source Code
-shall not exceed the cost of transferring the data.
-
-
- 5.3.2 DISTRIBUTION OF MODIFIED SOFTWARE
-
-When the Licensee makes a Contribution to the Software, the terms and
-conditions for the distribution of the resulting Modified Software
-become subject to all the provisions of this Agreement.
-
-The Licensee is authorized to distribute the Modified Software, in
-source code or object code form, provided that said distribution
-complies with all the provisions of the Agreement and is accompanied by:
-
- 1. a copy of the Agreement,
-
- 2. a notice relating to the limitation of both the Licensor's
- warranty and liability as set forth in Articles 8 and 9,
-
- and that, in the event that only the object code of the Modified
- Software is redistributed, the Licensee allows future Licensees
- unhindered access to the full source code of the Modified Software by
- indicating how to access it, it being understood that the additional
- cost of acquiring the source code shall not exceed the cost of
- transferring the data.
-
- 5.3 .3 DISTRIBUTION OF EXTERNAL MODULES
-
- When the Licensee has developed an External Module, the terms and
- conditions of this Agreement do not apply to said External Module, that
- may be distributed under a separate license agreement.
-
- 5.3 .4 COMPATIBILITY WITH THE GNU GPL
-
- The Licensee can include a code that is subject to the provisions of one
- of the versions of the GNU GPL in the Modified or unmodified Software,
- and distribute that entire code under the terms of the same version of
- the GNU GPL.
-
- The Licensee can include the Modified or unmodified Software in a code
- that is subject to the provisions of one of the versions of the GNU GPL,
- and distribute that entire code under the terms of the same version of
- the GNU GPL.
-
- Article 6 - INTELLECTUAL PROPERTY
-
- 6.1 OVER THE INITIAL SOFTWARE
-
- The Holder owns the economic rights over the Initial Software.Any or
- all use of the Initial Software is subject to compliance with the terms
- and conditions under which the Holder has elected to distribute its work
- and no one shall be entitled to modify the terms and conditions for the
- distribution of said Initial Software.
-
- The Holder undertakes that the Initial Software will remain ruled at
- least by this Agreement, for the duration set forth in Article 4.2.
-
- 6.2 OVER THE CONTRIBUTIONS
-
- The Licensee who develops a Contribution is the owner of the
- intellectual property rights over this Contribution as defined by
- applicable law.
-
- 6.3 OVER THE EXTERNAL MODULES
-
- The Licensee who develops an External Module is the owner of the
- intellectual property rights over this External Module as defined by
- applicable law and is free to choose the type of agreement that shall
- govern its distribution.
-
- 6.4 JOINT PROVISIONS
-
- The Licensee expressly undertakes :
-
- 1. not to remove, or modify, in any manner, the intellectual property
- notices attached to the Software;
-
- 2. to reproduce said notices, in an identical manner, in the copies
- of the Software modified or not.
-
- The Licensee undertakes not to directly or indirectly infringe the
- intellectual property rights of the Holder and/or Contributors on the
- Software and to take, where applicable, vis- -vis its staff, any and all
- measures required to ensure respect of said intellectual property rights
- of the Holder and/or Contributors.
-
- Article 7 - RELATED SERVICES
-
- 7.1 Under no circumstances shall the Agreement oblige the Licensor to
- provide technical assistance or maintenance services for the Software.
-
- However, the Licensor is entitled to offer this type of services.The
- terms and conditions of such technical assistance, and/or such
- maintenance, shall be set forth in a separate instrument.Only the
- Licensor offering said maintenance and/or technical assistance services
- shall incur liability therefor.
-
- 7.2 Similarly, any Licensor is entitled to offer to its licensees, under
- its sole responsibility, a warranty, that shall only be binding upon
- itself, for the redistribution of the Software and/or the Modified
- Software, under terms and conditions that it is free to decide.Said
- warranty, and the financial terms and conditions of its application,
- shall be subject of a separate instrument executed between the Licensor
- and the Licensee.
-
- Article 8 - LIABILITY
-
- 8.1 Subject to the provisions of Article 8.2, the Licensee shall be
- entitled to claim compensation for any direct loss it may have suffered
- from the Software as a result of a fault on the part of the relevant
- Licensor, subject to providing evidence thereof.
-
- 8.2 The Licensor
- 's liability is limited to the commitments made under
-this Agreement and shall not be incurred as a result of in particular:
-(i) loss due the Licensee's
- total or partial failure to fulfill its
- obligations, (ii) direct or consequential loss that is suffered by the
- Licensee due to the use or performance of the Software, and (iii) more
- generally, any consequential loss.In particular the Parties expressly
- agree that any or all pecuniary or business loss(i.e.loss of data,
- loss of profits, operating loss, loss of customers or orders,
- opportunity cost,
- any disturbance to business activities) or any or all
- legal proceedings instituted against the Licensee by a third party,
- shall constitute consequential loss and shall not provide entitlement to
- any or all compensation from the Licensor.
-
- Article 9 - WARRANTY
-
- 9.1 The Licensee acknowledges that the scientific and technical
- state-of-the-art when the Software was distributed did not enable all
- possible uses to be tested and verified, nor for the presence of
- possible defects to be detected.In this respect,
- the Licensee
- 's
-attention has been drawn to the risks associated with loading, using,
-modifying and/or developing and reproducing the Software which are
-reserved for experienced users.
-
-The Licensee shall be responsible for verifying, by any or all means,
-the suitability of the product for its requirements, its good working
-order, and for ensuring that it shall not cause damage to either persons
-or properties.
-
-9.2 The Licensor hereby represents, in good faith, that it is entitled
-to grant all the rights over the Software (including in particular the
-rights set forth in Article 5).
-
-9.3 The Licensee acknowledges that the Software is supplied "as is" by
-the Licensor without any other express or tacit warranty, other than
-that provided for in Article 9.2 and, in particular, without any warranty
-as to its commercial value, its secured, safe, innovative or relevant
-nature.
-
-Specifically, the Licensor does not warrant that the Software is free
-from any error, that it will operate without interruption, that it will
-be compatible with the Licensee's
- own equipment and software
- configuration,
- nor that it will meet the Licensee
- 's requirements.
-
-9.4 The Licensor does not either expressly or tacitly warrant that the
-Software does not infringe any third party intellectual property right
-relating to a patent, software or any other property right. Therefore,
-the Licensor disclaims any and all liability towards the Licensee
-arising out of any or all proceedings for infringement that may be
-instituted in respect of the use, modification and redistribution of the
-Software. Nevertheless, should such proceedings be instituted against
-the Licensee, the Licensor shall provide it with technical and legal
-assistance for its defense. Such technical and legal assistance shall be
-decided on a case-by-case basis between the relevant Licensor and the
-Licensee pursuant to a memorandum of understanding. The Licensor
-disclaims any and all liability as regards the Licensee's
- use of the
- name of the Software.No warranty is given as regards the existence of
- prior rights over the name of the Software or as regards the existence
- of a trademark.
-
- Article 10 - TERMINATION
-
- 10.1 In the event of a breach by the Licensee of its obligations
- hereunder, the Licensor may automatically terminate this Agreement
- thirty(30) days after notice has been sent to the Licensee and has
- remained ineffective.
-
- 10.2 A Licensee whose Agreement is terminated shall no longer be
- authorized to use, modify or distribute the Software.However, any
- licenses that it may have granted prior to termination of the Agreement
- shall remain valid subject to their having been granted in compliance
- with the terms and conditions hereof.
-
- Article 11 - MISCELLANEOUS
-
- 11.1 EXCUSABLE EVENTS
-
- Neither Party shall be liable for any or all delay, or failure to
- perform the Agreement, that may be attributable to an event of force
- majeure, an act of God or an outside cause, such as defective
- functioning or interruptions of the electricity or telecommunications
- networks, network paralysis following a virus attack, intervention by
- government authorities, natural disasters, water damage, earthquakes,
- fire, explosions, strikes and labor unrest, war, etc.
-
- 11.2 Any failure by either Party, on one or more occasions, to invoke
- one or more of the provisions hereof, shall under no circumstances be
- interpreted as being a waiver by the interested Party of its right to
- invoke said provision(s) subsequently.
-
- 11.3 The Agreement cancels and replaces any or all previous agreements,
- whether written or oral, between the Parties and having the same
- purpose, and constitutes the entirety of the agreement between said
- Parties concerning said purpose.No supplement or modification to the
- terms and conditions hereof shall be effective as between the Parties
- unless it is made in writing and signed by their duly authorized
- representatives.
-
- 11.4 In the event that one or more of the provisions hereof were to
- conflict with a current or future applicable act or legislative text,
- said act or legislative text shall prevail, and the Parties shall make
- the necessary amendments so as to comply with said act or legislative
- text.All other provisions shall remain effective.Similarly, invalidity
- of a provision of the Agreement, for any reason whatsoever, shall not
- cause the Agreement as a whole to be invalid.
-
- 11.5 LANGUAGE
-
- The Agreement is drafted in both French and English and both versions
- are deemed authentic.
-
- Article 12 - NEW VERSIONS OF THE AGREEMENT
-
- 12.1 Any person is authorized to duplicate and distribute copies of this
- Agreement.
-
- 12.2 So as to ensure coherence, the wording of this Agreement is
- protected and may only be modified by the authors of the License, who
- reserve the right to periodically publish updates or new versions of the
- Agreement, each with a separate number.These subsequent versions may
- address new issues encountered by Free Software.
-
- 12.3 Any Software distributed under a given version of the Agreement may
- only be subsequently distributed under the same version of the Agreement
- or a subsequent version, subject to the provisions of Article 5.3 .4.
-
- Article 13 - GOVERNING LAW AND JURISDICTION
-
- 13.1 The Agreement is governed by French law.The Parties agree to
- endeavor to seek an amicable solution to any disagreements or disputes
- that may arise during the performance of the Agreement.
-
- 13.2 Failing an amicable solution within two(2) months as from their
- occurrence, and unless emergency proceedings are necessary, the
- disagreements or disputes shall be referred to the Paris Courts having
- jurisdiction, by the more diligent Party.
-
- Version 2.0 dated 2006-0 9-05.
\ No newline at end of file
diff --git a/Code/Vahine/CenteringImageFilter.h b/Code/Vahine/CenteringImageFilter.h
deleted file mode 100644
index e7469835c019df20158ca3c2f0d28ee96a460270..0000000000000000000000000000000000000000
--- a/Code/Vahine/CenteringImageFilter.h
+++ /dev/null
@@ -1,74 +0,0 @@
-/**
- * CenteringImageFilter.h
- *
- * Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-#ifndef __VAHINECENTERINGIMAGE_H
-#define __VAHINECENTERINGIMAGE_H
-
-#include "VahineMacro.h"
-
-#include "itkImageToImageFilter.h"
-#include "itkStatisticsImageFilter.h"
-#include "itkShiftScaleImageFilter.h"
-
-namespace itk {
-
-template<class TInputImage,class TOutputImage> class ITK_EXPORT CenteringImageFilter :
- public ImageToImageFilter<TInputImage, TOutputImage> {
-public:
- /** Extract dimension from input and output image. */
- //itkStaticConstMacro(InputImageDimension, unsigned int,
- // TInputImage::ImageDimension);
- //itkStaticConstMacro(OutputImageDimension, unsigned int,
- // TOutputImage::ImageDimension);
-
- /** Standard Self typedef */
- typedef CenteringImageFilter Self;
- typedef ImageToImageFilter<TInputImage, TOutputImage> Superclass;
- typedef SmartPointer<Self> Pointer;
- typedef SmartPointer<const Self> ConstPointer;
-
- /** Method for creation through the object factory. */
- itkNewMacro(Self);
-
- /** Runtime information support. */
- itkTypeMacro(CenteringImageFilter, ImageToImageFilter);
-
- /** Image related typedefs. */
- typedef typename TInputImage::Pointer InputImagePointer;
- typedef typename TOutputImage::Pointer OutputImagePointer;
-
- /** NormalizeImageFilter must call modified on its internal filters */
- virtual void Modified() const;
-
-protected:
- CenteringImageFilter();
-
- /** GenerateData. */
- void GenerateData();
-
- // Override since the filter needs all the data for the algorithm
- void GenerateInputRequestedRegion();
-
-private:
- CenteringImageFilter(const Self&); //purposely not implemented
- void operator=(const Self&); //purposely not implemented
-
- typename StatisticsImageFilter<TInputImage>::Pointer m_StatisticsFilter;
- typename ShiftScaleImageFilter<TInputImage,TOutputImage>::Pointer m_ShiftScaleFilter;
-}; // end of class
-
-} // end namespace itk
-
-#include "CenteringImageFilter.txx"
-
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/CenteringImageFilter.txx b/Code/Vahine/CenteringImageFilter.txx
deleted file mode 100644
index 2d73831c6d9cf99904b57bfc5984e0a2befd62b2..0000000000000000000000000000000000000000
--- a/Code/Vahine/CenteringImageFilter.txx
+++ /dev/null
@@ -1,82 +0,0 @@
-/**
- * CenteringImageFilter.txx
- *
- * Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-#ifndef __VAHINECENTERING_TXX
-#define __VAHINECENTERING_TXX
-
-#include "CenteringImageFilter.h"
-#include "itkImageRegionIterator.h"
-#include "itkStatisticsImageFilter.h"
-#include "itkShiftScaleImageFilter.h"
-#include "itkCommand.h"
-#include "itkProgressAccumulator.h"
-
-namespace itk {
-
-template <class TInputImage, class TOutputImage>
-CenteringImageFilter<TInputImage, TOutputImage>::CenteringImageFilter() {
- m_StatisticsFilter = 0;
- m_StatisticsFilter = StatisticsImageFilter<TInputImage>::New();
- m_ShiftScaleFilter = ShiftScaleImageFilter<TInputImage,TOutputImage>::New();
-}
-
-template <class TInputImage, class TOutputImage>
-void CenteringImageFilter<TInputImage, TOutputImage>::GenerateInputRequestedRegion() {
- Superclass::GenerateInputRequestedRegion();
- if ( this->GetInput() ) {
- InputImagePointer image = const_cast< typename Superclass::InputImageType * >( this->GetInput() );
- image->SetRequestedRegionToLargestPossibleRegion();
- }
-}
-
-template <class TInputImage, class TOutputImage>
-void CenteringImageFilter<TInputImage, TOutputImage>::Modified() const {
- Superclass::Modified();
- m_StatisticsFilter->Modified();
- m_ShiftScaleFilter->Modified();
-}
-
-template <class TInputImage, class TOutputImage>
-void CenteringImageFilter<TInputImage, TOutputImage>::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->Update();
-
- // Set the parameters for Shift
- vahineDebug("CenteringImageFilter : mean = "<< m_StatisticsFilter->GetMean() );
- m_ShiftScaleFilter->SetShift(-m_StatisticsFilter->GetMean() );
- m_ShiftScaleFilter->SetInput(this->GetInput() );
-
- m_ShiftScaleFilter->GraftOutput(this->GetOutput() );
- //m_ShiftScaleFilter->GetOutput()->SetRequestedRegion(this->GetOutput()->GetRequestedRegion());
- m_ShiftScaleFilter->Update();
-
- //std::cout<<"shiftscale = "<<m_ShiftScaleFilter->GetOutput()<<std::endl;
-
- // Graft the mini pipeline output to this filters output
- this->GraftOutput(m_ShiftScaleFilter->GetOutput() );
-}
-
-} // end namespace itk
-
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/CovPixelSelectFilter.h b/Code/Vahine/CovPixelSelectFilter.h
deleted file mode 100644
index aecc3555ad9738d17776c4e611d876a77183c7f7..0000000000000000000000000000000000000000
--- a/Code/Vahine/CovPixelSelectFilter.h
+++ /dev/null
@@ -1,148 +0,0 @@
-/*
- * CovPixelSelectFilter.h
- *
- * Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINEFILTRAGEFILTER_H_
-#define __VAHINEFILTRAGEFILTER_H_
-
-#include <otb/Utilities/ITK/Utilities/vxl/core/vnl/algo/vnl_symmetric_eigensystem.h>
-#include <otb/Utilities/ITK/Utilities/vxl/core/vnl/algo/vnl_real_eigensystem.h>
-
-#include <otb/IO/otbImage.h>
-#include <otb/IO/otbVectorImage.h>
-#include <otb/Utilities/ITK/Common/itkImageToImageFilter.h>
-
-#include <otb/BasicFilters/otbStreamingStatisticsVectorImageFilter.h>
-#include <otb/BasicFilters/otbImportImageFilter.h>
-#include "itkVariableSizeMatrix.h"
-#include "itkArray.h"
-#include <otb/Utilities/ITK/Common/itkVector.h>
-#include <otb/Utilities/ITK/Numerics/Statistics/itkListSample.h>
-#include <otb/Utilities/ITK/Numerics/Statistics/itkGaussianMixtureModelComponent.h>
-#include <otb/Utilities/ITK/Numerics/Statistics/itkExpectationMaximizationMixtureModelEstimator.h>
-
-#include "VahineMacro.h"
-
-namespace otb {
-
-template <class TInputImage, class TOutputImage>
-class ITK_EXPORT VahineCovPixelSelectFilter
- : public itk::ImageToImageFilter<TInputImage, TOutputImage> {
-public:
- typedef TInputImage InputImageType;
- typedef TOutputImage OutputImageType;
-
- typedef VahineCovPixelSelectFilter Self;
- typedef itk::ImageToImageFilter<InputImageType, OutputImageType> Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
-
- // creation of SmartPointer
- itkNewMacro(Self);
- // runtime type information
- itkTypeMacro(VahineCovPixelSelectFilter, itk::ImageToImageFilter);
-
- // Input Image
- typedef typename InputImageType::Pointer ImagePointer;
- typedef typename InputImageType::PixelType PixelType;
- typedef typename InputImageType::RegionType RegionType;
- typedef typename InputImageType::InternalPixelType InternalPixelType;
- typedef typename itk::NumericTraits<InternalPixelType>::RealType RealType;
- typedef typename itk::ImageRegionConstIterator<InputImageType> ImageRegionConstIterator;
-
- // Output Image
- typedef typename OutputImageType::Pointer OutputImagePointer;
- typedef typename OutputImageType::PixelType OutputPixelType;
- typedef typename OutputImageType::IndexType OutputIndexType;
- typedef typename OutputImageType::SizeType OutputSizeType;
- typedef typename OutputImageType::RegionType OutputRegionType;
- typedef typename itk::ImageRegionIterator<OutputImageType> OutputRegionIterator;
-
- typedef StreamingStatisticsVectorImageFilter<InputImageType> StatsFilterType;
- typedef typename StatsFilterType::Pointer StatsFilterPointer;
-
- typedef typename itk::VariableSizeMatrix<RealType> MatrixType;
- typedef itk::Array<double> ItkArrayType;
- typedef vnl_vector<RealType> ArrayType;
- typedef typename std::vector<MatrixType> ArrayMatrixType;
-
- typedef itk::Vector<double, 1> MeasurementVectorType;
- typedef itk::Statistics::ListSample<MeasurementVectorType> SampleType;
- typedef itk::Statistics::GaussianMixtureModelComponent<SampleType> ComponentType;
- typedef itk::Statistics::ExpectationMaximizationMixtureModelEstimator<SampleType> EstimatorType;
-
- bool GetDebug(){
- return true;
- };
-
- itkSetMacro( ReferenceImage, ImagePointer );
- itkGetMacro( PCAVectors, MatrixType );
- itkGetMacro( PCAValues, ArrayType*);
- itkSetMacro( Explanation, double );
- itkGetMacro( Explanation, double );
- itkGetMacro( NumberOfPCAAxes, unsigned int );
- itkGetMacro( PCAReferenceProjections, ArrayMatrixType* );
- itkGetMacro( PCAInputProjections, ArrayMatrixType* );
- itkGetMacro( Distances, MatrixType* );
- itkGetMacro( GMMMean, ItkArrayType );
- itkGetMacro( GMMCovariance, ItkArrayType );
- itkGetMacro( GMMProportion, ItkArrayType );
- itkSetMacro( BoundDistance, ItkArrayType );
- itkSetMacro( DistanceWeight, ItkArrayType );
-protected:
- VahineCovPixelSelectFilter();
- virtual ~VahineCovPixelSelectFilter();
- void PrintSelf(std::ostream& os, itk::Indent indent) const;
-
- void GenerateData();
-
- unsigned int m_NumberOfPCAAxes;
- virtual void CalculatePCAVectors();
-
- virtual void CalculateNNS();
-
- virtual void UpdateNumberOfPCAAxes();
-
- virtual void CalculatePCAProjection(const unsigned int ImageNumberOfBands, ImagePointer image,
- ArrayMatrixType *matrix);
-
- virtual void ClassificationEM();
-
- virtual void WeightFilter();
-
-private:
- ImagePointer m_ReferenceImage;
- MatrixType m_PCAVectors;
- ArrayType* m_PCAValues;
- /* pourcent of explanation */
- double m_Explanation;
-
- ArrayMatrixType* m_PCAReferenceProjections;
- ArrayMatrixType* m_PCAInputProjections;
- MatrixType* m_Distances;
- ItkArrayType m_GMMMean;
- ItkArrayType m_GMMCovariance;
- ItkArrayType m_GMMProportion;
- ItkArrayType m_BoundDistance;
- ItkArrayType m_DistanceWeight;
-
- VahineCovPixelSelectFilter(const Self&); //not implemented
- void operator=(const Self&); //not implemented
-
-};
-
-} // end namespace otb
-
-#include "CovPixelSelectFilter.txx"
-
-#endif /* VAHINEFILTRAGEFILTER_H_ */
\ No newline at end of file
diff --git a/Code/Vahine/CovPixelSelectFilter.txx b/Code/Vahine/CovPixelSelectFilter.txx
deleted file mode 100644
index 6afb8e29f3c515187d7c2f81766a8aa08ec3a03a..0000000000000000000000000000000000000000
--- a/Code/Vahine/CovPixelSelectFilter.txx
+++ /dev/null
@@ -1,305 +0,0 @@
-/*
- * CovPixelSelectFilter.txx
- *
- * Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-#ifndef __VAHINEFILTRAGEFILTER_TXX
-#define __VAHINEFILTRAGEFILTER_TXX
-#include "CovPixelSelectFilter.h"
-#include "limits"
-
-namespace otb {
-
-template<class TInputImage, class TOutputImage>
-VahineCovPixelSelectFilter<TInputImage,TOutputImage>::VahineCovPixelSelectFilter() :
- m_NumberOfPCAAxes(0),
- m_Explanation(95.0) {
- vahineDebug("VahineCovPixelSelectFilter constructor");
-}
-
-template<class TInputImage, class TOutputImage>
-VahineCovPixelSelectFilter<TInputImage,TOutputImage>::~VahineCovPixelSelectFilter(){
- /*if(m_PCAReferenceProjections != NULL){
- delete m_PCAReferenceProjections;
- }
-delete m_PCAInputProjections;
-delete m_Distances;
-delete m_PCAValues;*/
-}
-
-/**
- * vnl_symmetric_eigensystem :
- * The columns of V are the eigenvectors, sorted by increasing eigenvalue,
- * from most negative to most positive.
- */
-template<class TInputImage, class TOutputImage>
-void VahineCovPixelSelectFilter<TInputImage, TOutputImage>::CalculatePCAVectors(){
- vahineDebug("CalculatePCAVectors");
- StatsFilterPointer statsFilter = StatsFilterType::New();
- statsFilter->SetInput(m_ReferenceImage);
- statsFilter->Update();
- MatrixType covarianceMatrix = statsFilter->GetCovariance();
- unsigned int N = covarianceMatrix.Cols();
- covarianceMatrix *= static_cast<double>(N)/(N-1);
- vnl_symmetric_eigensystem<RealType> eigen(covarianceMatrix.GetVnlMatrix() );
- m_PCAVectors = eigen.V;
- m_PCAValues = new ArrayType();
- (*m_PCAValues) = eigen.D.diagonal();
- double sum = static_cast<double>(m_PCAValues->sum() );
- (*m_PCAValues) *= 100.0/sum;
-}
-
-/**
- * Update the number of interresting PCA axes
- * for an explanation pourcentage
- */
-template<class TInputImage, class TOutputImage>
-void VahineCovPixelSelectFilter<TInputImage, TOutputImage>::UpdateNumberOfPCAAxes(){
- vahineDebug("UpdateNumberOfPCAAxes");
- double explanationSum = 0.0;
- unsigned int numberOfVectors = 0;
- for(unsigned int i = m_PCAValues->size(); i>0; i--) {
- explanationSum += (*m_PCAValues)(i-1);
- numberOfVectors++;
- vahineDebug("### explanation "<<explanationSum<<" "<<m_Explanation);
- if(explanationSum >= m_Explanation) {
- vahineDebug("Break number of vector = "<<numberOfVectors);
- break;
- }
- }
- m_NumberOfPCAAxes = numberOfVectors;
-}
-
-template<class TInputImage, class TOutputImage>
-void VahineCovPixelSelectFilter<TInputImage, TOutputImage>::CalculatePCAProjection(
- const unsigned int ImageNumberOfBands, ImagePointer image, ArrayMatrixType * matrix){
- vahineDebug("CalculatePCAProjection");
- MatrixType tempMatrix;
- RegionType region = image->GetLargestPossibleRegion();
- // size along X
- unsigned int cols = region.GetSize(0);
- // size along Y
- unsigned int rows = region.GetSize(1);
- vahineDebug("### size "<<rows<<","<<cols);
- tempMatrix.SetSize(rows, cols);
- tempMatrix.Fill(itk::NumericTraits<InternalPixelType>::Zero);
- matrix->assign(m_NumberOfPCAAxes, tempMatrix);
-
- ImageRegionConstIterator imageIt(image, region);
- // number of total PCA axes
- //const unsigned int nbVectors = m_PCAVectors.Cols();
- unsigned int i = 0;
- for(imageIt.GoToBegin(); !imageIt.IsAtEnd(); ++imageIt) {
- PixelType pixel = imageIt.Get();
- const InternalPixelType * buffer = pixel.GetDataPointer();
- vnl_vector<InternalPixelType> pix = vnl_vector<InternalPixelType>(buffer, ImageNumberOfBands);
- // convert from InternalPixelType to RealType
- ArrayType mypix( pix.size() );
- for ( unsigned int j=0; j<pix.size(); ++j ) {
- mypix[j] = static_cast<RealType>(pix[j]);
- }
- for(unsigned int n=0; n<m_NumberOfPCAAxes; n++) {
- //cout<<"### "<<i<<" : "<<floor(i/cols)<<","<< i%cols<<endl;
- ArrayType currentPCAVector = m_PCAVectors.GetVnlMatrix().get_column(ImageNumberOfBands-1-n);
- //cout<<"# PCA vect "<<currentPCAVector<<endl;
- matrix->at(n) (static_cast<unsigned int>(floor(i/cols) ), i%cols) = dot_product(mypix, currentPCAVector);
- }
- i++;
- }
-
-}
-
-/**
- * Nearest neighbor search : euclidian distance
- * //TODO use a another algo like kd-tree for decrease time calculation
- */
-template<class TInputImage, class TOutputImage>
-void VahineCovPixelSelectFilter<TInputImage, TOutputImage>::CalculateNNS(){
- vahineDebug("CalculateNNS");
- m_Distances = new MatrixType();
- // size of input image
- unsigned int icols = m_PCAInputProjections->at(1).Cols();
- unsigned int irows = m_PCAInputProjections->at(1).Rows();
- m_Distances->SetSize(irows, icols);
- // size of reference image
- unsigned int rcols = m_PCAReferenceProjections->at(1).Cols();
- unsigned int rrows = m_PCAReferenceProjections->at(1).Rows();
- double inCoord, refCoord, sum, distance, curDistance;
- std::numeric_limits<double> doubleLimits;
- m_BoundDistance = ItkArrayType(2);
- m_BoundDistance[0] = doubleLimits.max();
- m_BoundDistance[1] = doubleLimits.min();
- for(unsigned int i=0; i<irows; ++i) {
- for(unsigned int j=0; j<icols; ++j) {
- (*m_Distances)(i,j) = doubleLimits.max();
- for(unsigned int k=0; k<rrows; ++k) {
- for(unsigned int h=0; h<rcols; ++h) {
- sum = 0.0;
- for(unsigned int axe=0; axe<m_NumberOfPCAAxes; ++axe) {
- inCoord = m_PCAInputProjections->at(axe) (i,j);
- refCoord = m_PCAReferenceProjections->at(axe) (k,h);
- sum += pow(inCoord-refCoord,2);
- }
- curDistance = pow(sum, 0.5);
- if(curDistance < (*m_Distances)(i,j) ) {
- (*m_Distances)(i,j) = curDistance;
- if(m_BoundDistance[0]> curDistance) {
- m_BoundDistance[0] = curDistance;
- }
- if(m_BoundDistance[1]<curDistance) {
- m_BoundDistance[1] = curDistance;
- }
- }
- }
- }
- }
- }
-}
-
-template<class TInputImage, class TOutputImage>
-void VahineCovPixelSelectFilter<TInputImage, TOutputImage>::ClassificationEM(){
- vahineDebug("ClassificationEM");
- unsigned int numberOfClasses = 2;
- m_GMMMean = ItkArrayType();
- m_GMMCovariance = ItkArrayType();
- m_GMMProportion = ItkArrayType();
- m_GMMMean.SetSize(numberOfClasses);
- m_GMMCovariance.SetSize(numberOfClasses);
- m_GMMProportion.SetSize(numberOfClasses);
-
- SampleType::Pointer sample = SampleType::New();
- sample->SetMeasurementVectorSize(1);
- MeasurementVectorType mv;
- unsigned int cols = m_Distances->Cols();
- unsigned int rows = m_Distances->Rows();
- std::cout<<"rows = "<<rows<<" cols "<<cols<<std::endl;
- for(unsigned int i=0; i<rows; ++i) {
- for(unsigned int j=0; j<cols; ++j) {
- mv[0] = (*m_Distances)(i,j);
- sample->PushBack(mv);
- }
- }
-
- std::vector<ItkArrayType> initialParameters(numberOfClasses);
- double delta = (m_BoundDistance[1]-m_BoundDistance[0])/3;
- // parameters of first gaussian
- initialParameters[0] = ItkArrayType(2);
- initialParameters[0][0] = m_BoundDistance[0] + delta;
- initialParameters[0][1] = delta;
- // parameters of second gaussian
- initialParameters[1] = ItkArrayType(2);
- initialParameters[1][0] = m_BoundDistance[1] - delta;
- initialParameters[1][1] = delta;
-
- ItkArrayType initialProportions(numberOfClasses);
-
- EstimatorType::Pointer estimator = EstimatorType::New();
- estimator->SetSample(sample);
- estimator->SetMaximumIteration(1000);
-
- std::vector<ComponentType::Pointer> components;
- for(unsigned int i=0; i<numberOfClasses; ++i) {
- components.push_back(ComponentType::New() );
- (components[i])->SetSample(sample);
- (components[i])->SetParameters(initialParameters[i]);
- estimator->AddComponent( (ComponentType::Superclass*)( (components[i]).GetPointer() ) );
- initialProportions[i] = 1.0/numberOfClasses;
- }
-
- estimator->SetInitialProportions(initialProportions);
- estimator->Update();
-
- for(unsigned int i=0; i<numberOfClasses; ++i) {
- m_GMMMean[i] = (components[i])->GetFullParameters()[0];
- m_GMMCovariance[i] = (components[i])->GetFullParameters()[1];
- //std::cout<<" params "<<(components[i])->GetFullParameters()<<std::endl;
- //std::cout<<" prop "<<(*estimator->GetProportions())[i]<<std::endl;
- //std::cout<<" compo "<<(components[i])<<std::endl;
- }
- ComponentType::WeightArrayType *firstCompoWeight = (components[0])->GetWeights();
- m_DistanceWeight = ItkArrayType();
- m_DistanceWeight.SetSize(firstCompoWeight->GetSize() );
- m_DistanceWeight = (*firstCompoWeight);
-}
-
-template<class TInputImage, class TOutputImage>
-void VahineCovPixelSelectFilter<TInputImage, TOutputImage>::WeightFilter(){
- vahineDebug("WeightFilter");
- // size along X
- unsigned int cols = this->GetInput()->GetLargestPossibleRegion().GetSize(0);
- // size along Y
- unsigned int rows = this->GetInput()->GetLargestPossibleRegion().GetSize(1);
- //std::cout<<"cols "<<cols<<" "<<"rows "<<rows<<endl;
- // create mask
- OutputImagePointer mask = OutputImageType::New();
- OutputIndexType start;
- start[0] = 0;
- start[1] = 0;
- OutputSizeType size;
- size[0] = cols;
- size[1] = rows;
- OutputRegionType region;
- region.SetSize(size);
- region.SetIndex(start);
- mask->SetRegions(region);
- mask->Allocate();
- vahineDebug("weight size "<<m_DistanceWeight.Size() );
- vahineDebug("mask size "<<cols*rows);
- OutputRegionIterator maskIt(mask, mask->GetLargestPossibleRegion() );
- ItkArrayType::iterator weightIt = m_DistanceWeight.begin();
- for(maskIt.GoToBegin(); !maskIt.IsAtEnd(); ++maskIt, ++weightIt) {
- if( (*weightIt)>0.5) { // we have two class
- maskIt.Set(static_cast<OutputPixelType>(1.0) );
- } else {
- maskIt.Set(static_cast<OutputPixelType>(0.0) );
- }
-
- }
- this->GraftOutput(mask);
- vahineDebug("end weightfilter");
-}
-
-template<class TInputImage, class TOutputImage>
-void VahineCovPixelSelectFilter<TInputImage, TOutputImage>::GenerateData(){
- m_ReferenceImage->UpdateOutputInformation();
- const unsigned int ReferenceImageNumberOfBands = m_ReferenceImage->GetNumberOfComponentsPerPixel();
- const unsigned int InputNumberOfBands = this->GetInput()->GetNumberOfComponentsPerPixel();
- if( ReferenceImageNumberOfBands != InputNumberOfBands ) {
- itkExceptionMacro(
- "The number of bands of the reference image ("<<ReferenceImageNumberOfBands<<
- ") is different from the number of bands of the input image ("<<InputNumberOfBands<<"). ");
- }
- CalculatePCAVectors();
- UpdateNumberOfPCAAxes();
- m_PCAReferenceProjections = new ArrayMatrixType();
- vahineDebug("First PCA");
- m_ReferenceImage->Update();
- CalculatePCAProjection(ReferenceImageNumberOfBands, m_ReferenceImage, m_PCAReferenceProjections);
- m_PCAInputProjections = new ArrayMatrixType();
- vahineDebug("Second PCA");
- CalculatePCAProjection(ReferenceImageNumberOfBands,
- const_cast<InputImageType *>(this->GetInput() ), m_PCAInputProjections);
- CalculateNNS();
- ClassificationEM();
- WeightFilter();
-}
-
-/**
- * Standard "PrintSelf" method
- */
-template<class TInputImage, class TOutputImage>
-void VahineCovPixelSelectFilter<TInputImage, TOutputImage>::PrintSelf(std::ostream& os, itk::Indent indent) const {
- Superclass::PrintSelf( os, indent );
-}
-
-} // end namespace otb
-
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/ElmVcaFilter.h b/Code/Vahine/ElmVcaFilter.h
deleted file mode 100644
index 118bbb17d8877a94f5b2ea2c58e43935329194c8..0000000000000000000000000000000000000000
--- a/Code/Vahine/ElmVcaFilter.h
+++ /dev/null
@@ -1,170 +0,0 @@
-/*
- * ElmVcaFilter.h
- * Date: $Date$
- * Version: $Revision$
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2010 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef ELMVCAFILTER_H_
-#define ELMVCAFILTER_H_
-
-#include <limits>
-#include <cmath>
-#include <ctime>
-#include <vector>
-#include <list>
-#include "vnl/algo/vnl_symmetric_eigensystem.h"
-#include <vcl_algorithm.h>
-#include <vcl_vector.h>
-
-#include <otbStreamingStatisticsVectorImageFilter.h>
-#include <otbStreamingStatisticsImageFilter.h>
-#include <otbVectorImageToImageListFilter.h>
-#include <otbImageList.h>
-#include <otbImage.h>
-#include <otbVectorImage.h>
-
-#include <itkImageToImageFilter.h>
-#include <itkImageRegionConstIterator.h>
-#include <itkVariableLengthVector.h>
-#include <itkCovarianceCalculator.h>
-#include <itkImageToListAdaptor.h>
-#include <itkListSample.h>
-
-#include "VahineMacro.h"
-#include "VcaFilter.h"
-
-namespace otb {
-template <class TImage>
-class ITK_EXPORT VahineElmVcaFilter
- : public itk::ImageToImageFilter<TImage, TImage> {
-public:
- typedef TImage VectorImageType;
-
- typedef VahineElmVcaFilter Self;
- typedef itk::ImageToImageFilter<VectorImageType, VectorImageType> Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
-
- // creation of SmartPointer
- itkNewMacro(Self);
- // runtime type information
- itkTypeMacro(VahineElmVcaFilter, itk::ImageToImageFilter);
-
- // Input/Output Image
- typedef typename VectorImageType::Pointer VectorImagePointer;
- typedef typename VectorImageType::ConstPointer VectorImageConstPointer;
- typedef typename VectorImageType::PixelType PixelType;
- typedef typename VectorImageType::InternalPixelType InternalPixelType;
- typedef typename VectorImageType::RegionType RegionType;
- typedef typename VectorImageType::SizeType SizeType;
-
- // type for computation
- const static unsigned int Dimension = 2;
- typedef Image<InternalPixelType, Dimension> BandType;
-
- typedef VahineVCAFilter<VectorImageType> VCAFilterType;
- typedef typename VCAFilterType::Pointer VCAFilterPointer;
- typedef typename std::vector<unsigned int> IntegerList;
-
- bool GetDebug(){
- return true;
- };
- itkGetMacro( NbComponents, unsigned int );
- itkSetMacro( NbComponents, unsigned int );
- IntegerList GetMaxLikeHood();
-
- /**
- * MVSAFilter produces an image with same resolution that input image but
- * with a different number of bands. As such, this filter needs to provide an
- * implementation for GenerateOutputInformation() in order
- * to inform the pipeline execution model.
- */
- virtual void GenerateOutputInformation();
-
- // for testing
- friend class ElmVcaFilterFixture;
-protected:
- VahineElmVcaFilter();
- virtual ~VahineElmVcaFilter();
-
- void PrintSelf(std::ostream& os, itk::Indent indent) const;
-
- void GenerateData();
-
- unsigned int m_NumberOfBand; // protected for testing
-private:
- // internal type for computation
- typedef vnl_matrix<double> VnlMatrix;
- typedef vnl_vector<double> VnlVector;
- typedef vnl_vector<InternalPixelType> InternalVnlVector;
- typedef vnl_matrix<InternalPixelType> InternalVnlMatrix;
-
- typedef itk::ImageRegionConstIterator<VectorImageType> ConstIteratorType;
- typedef itk::ImageRegionIterator<VectorImageType> IteratorType;
- typedef otb::StreamingStatisticsVectorImageFilter<VectorImageType> StatsFilterType;
- typedef typename StatsFilterType::Pointer StatsFilterPointer;
-
- typedef ImageList<BandType> ImageListType;
- typedef typename ImageListType::Pointer ImageListPointer;
- typedef typename ImageListType::ConstIterator ConstListIteratorType;
- typedef VectorImageToImageListFilter<VectorImageType, ImageListType> DecomposeFilterType;
- typedef typename DecomposeFilterType::Pointer DecomposeFilterPointer;
- typedef StreamingStatisticsImageFilter<BandType> BandStatsFilterType;
- typedef typename BandStatsFilterType::Pointer BandStatsFilterPointer;
-
- /** Typedef for statistic computing. */
- typedef StreamingStatisticsVectorImageFilter<VectorImageType> StreamingStatisticsVectorImageFilterType;
- typedef typename StreamingStatisticsVectorImageFilterType::MatrixType MatrixType;
-
- typedef typename itk::Statistics::ListSample< itk::VariableLengthVector< InternalPixelType > > SampleType;
- typedef typename itk::Statistics::CovarianceCalculator< SampleType >
- CovarianceAlgorithmType;
- typedef typename CovarianceAlgorithmType::Pointer
- CovarianceAlgorithmPointer;
-
- std::numeric_limits<InternalPixelType> m_RealLimits;
- unsigned int m_NbComponents;
- IntegerList m_MaxLikeHood;
- VnlVector m_Mean;
-
- VnlMatrix myCovariance(VectorImagePointer image, unsigned int nbObserv);
-
- VnlMatrix streamingCovariance(VectorImagePointer image, unsigned int nbObserv);
-
- unsigned int Demelange(unsigned int);
-
- void LocalMaximum(VnlVector & likeHood);
-
- void SpectralMean(VnlVector & mean);
-
- void LikelihoodLog(VnlVector & eigenR, VnlVector & eigenK, unsigned int N, VnlVector & lH);
-
- VahineElmVcaFilter(const Self&); //not implemented
- void UpdateNumberOfComponents();
-
- void operator=(const Self&); //not implemented
-
- static double mysquare(const double x){
- return x*x;
- }
- static double mylog(const double x){
- return log(x);
- }
-
-};
-
-} // end namesapce otb
-
-#include "ElmVcaFilter.txx"
-
-#endif /* MVSAFILTER_H_ */
diff --git a/Code/Vahine/ElmVcaFilter.txx b/Code/Vahine/ElmVcaFilter.txx
deleted file mode 100644
index da146482cf685f1efab61d275931b5f7f3daeaf8..0000000000000000000000000000000000000000
--- a/Code/Vahine/ElmVcaFilter.txx
+++ /dev/null
@@ -1,324 +0,0 @@
-/*
- * MVSAFilter.txx
- * Date: $Date$
- * Version: $Revision$
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2010 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINEELMVCAFILTER_TXX
-#define __VAHINEELMVCAFILTER_TXX
-
-#include "ElmVcaFilter.h"
-
-namespace otb {
-
-template<class TImage>
-VahineElmVcaFilter<TImage>::VahineElmVcaFilter() :
- m_NbComponents(0) {
-
-}
-
-template<class TImage>
-VahineElmVcaFilter<TImage>::~VahineElmVcaFilter(){
-
-}
-
-template<class TImage>
-void VahineElmVcaFilter<TImage>::UpdateNumberOfComponents(){
- RegionType imageRegion = this->GetInput()->GetLargestPossibleRegion();
- SizeType size = imageRegion.GetSize();
- unsigned int NObserv = size[0]*size[1];
-
- m_NumberOfBand = this->GetInput()->GetNumberOfComponentsPerPixel();
- if(m_NbComponents == 0) {
- m_NbComponents = Demelange(NObserv);
- std::cout << "Using " << m_NbComponents << " endmembers" << std::endl;
- }
-}
-
-template<class TImage>
-void VahineElmVcaFilter<TImage>::GenerateData(){
- RegionType imageRegion = this->GetInput()->GetLargestPossibleRegion();
-
- VCAFilterPointer vcaFilter = VCAFilterType::New();
-
- vcaFilter->SetInput(this->GetInput() );
- vcaFilter->SetNumberEndMember(m_NbComponents);
- vcaFilter->SetDebug(false);
- vcaFilter->Update();
-
- VnlMatrix mixingMatrix = vcaFilter->GetMixingMatrix();
- vahineDebug("mixingMatrix=\n"); //<<mixingMatrix);
- vnl_matrix_inverse<double> PseudoInv(mixingMatrix);
- VnlMatrix pinv = PseudoInv.pinverse();
- vahineDebug("pinv=\n"); //<<pinv);
-
- VectorImagePointer output = VectorImageType::New();
- output->SetRegions(imageRegion);
- output->SetNumberOfComponentsPerPixel(m_NbComponents);
- output->Allocate();
-
- // fill output
- ConstIteratorType itInput( this->GetInput(), imageRegion );
- IteratorType itOutput( output, imageRegion );
-
- InternalVnlMatrix realPinv(pinv.rows(), pinv.cols() );
- VnlMatrix::iterator itPinv = pinv.begin();
- typename InternalVnlMatrix::iterator itRealPinv = realPinv.begin();
-
- while(itPinv != pinv.end() ) {
- (*itRealPinv) = static_cast<InternalPixelType>(*itPinv);
- ++itRealPinv; ++itPinv;
- }
-
- for(itInput.GoToBegin(), itOutput.GoToBegin(); !itInput.IsAtEnd(); ++itInput, ++itOutput) {
- PixelType pixelIn = itInput.Get();
- PixelType pixelOut = itOutput.Get();
- InternalVnlVector pixIn(pixelIn.GetDataPointer(), pixelIn.GetNumberOfElements() );
- InternalVnlVector u = realPinv * pixIn;
- //vahineDebug("u = "<<u);
- // TODO optimize with mem copy ?
- for(unsigned int i=0; i<u.size(); ++i) {
- pixelOut[i] = u(i);
- }
- itOutput.Set(pixelOut);
- }
- vahineDebug("end output");
- this->GraftOutput(output);
-}
-
-/**
- * Compute the mean value of each spectral band
- */
-template<class TImage>
-void VahineElmVcaFilter<TImage>::SpectralMean(VnlVector & mean){
- mean.set_size(m_NumberOfBand);
- VectorImagePointer image = const_cast<VectorImageType *>(this->GetInput() );
- DecomposeFilterPointer decomposer = DecomposeFilterType::New();
- decomposer->SetInput(image);
- ImageListPointer imageList = decomposer->GetOutput();
- imageList->Update();
- unsigned int i = 0;
- for (ConstListIteratorType it = imageList->Begin(); it!=imageList->End(); ++it) {
- // Compute the mean for each band with the statistical filter
- BandStatsFilterPointer meanFilter = BandStatsFilterType::New();
- meanFilter->SetInput(it.Get() );
- meanFilter->Update();
- mean(i) = static_cast<double>(meanFilter->GetMean() );
- ++i;
- }
- vahineDebug("mean = "<<mean);
-}
-
-template<class TImage>
-typename VahineElmVcaFilter<TImage>::VnlMatrix VahineElmVcaFilter<TImage>::myCovariance(VectorImagePointer image,
- unsigned int nbObserv){
- image->Update();
- typename SampleType::Pointer sample = SampleType::New();
- RegionType imageRegion = this->GetInput()->GetLargestPossibleRegion();
- ConstIteratorType it( this->GetInput(), imageRegion );
- for(it.GoToBegin(); !it.IsAtEnd(); ++it) {
- sample->PushBack( it.Get() );
- //vahineDebug("it "<<it.Get() );
- }
-
- CovarianceAlgorithmPointer covarianceAlgorithm = CovarianceAlgorithmType::New();
- covarianceAlgorithm->SetInputSample(sample);
- covarianceAlgorithm->SetMean(0);
- covarianceAlgorithm->Update();
- const itk::VariableSizeMatrix<double> *covarianceVSMat = covarianceAlgorithm->GetOutput();
- VnlMatrix covarianceMatrix = covarianceVSMat->GetVnlMatrix();
- typename CovarianceAlgorithmType::MeanType *mean = covarianceAlgorithm->GetMean();
- m_Mean.set_size(mean->size() );
- m_Mean.copy_in(mean->data_block() );
- return covarianceMatrix;
-}
-
-template<class TImage>
-typename VahineElmVcaFilter<TImage>::VnlMatrix VahineElmVcaFilter<TImage>::streamingCovariance(VectorImagePointer image,
- unsigned int nbObserv){
- image->Update();
- typename StreamingStatisticsVectorImageFilterType::Pointer covComputefilter =
- StreamingStatisticsVectorImageFilterType::New();
- covComputefilter->SetInput(image);
- covComputefilter->Update();
- VnlMatrix covarianceMatrix = (covComputefilter->GetCovariance() ).GetVnlMatrix();
- PixelType mean = covComputefilter->GetMean();
- m_Mean.set_size(mean.Size() );
- for(unsigned int i=0; i< mean.Size(); ++i) {
- m_Mean[i] = mean[i];
- }
- double regul = static_cast<double>(nbObserv)/(nbObserv-1);
- covarianceMatrix *= regul;
- return covarianceMatrix;
-}
-
-/**
- * Find the number of element by calculation of :
- * covariance and correlation matrix
- * take care of correlation calculation, the matrix is normalize
- * with the number of observations.
- */
-template<class TImage>
-unsigned int VahineElmVcaFilter<TImage>::Demelange(unsigned int nbObserv){
- std::time_t start = time(NULL);
-
- vahineDebug("statistic calculation start");
- VectorImagePointer image = const_cast<VectorImageType *>(this->GetInput() );
- VnlMatrix covarianceMatrix;
- //covarianceMatrix = myCovariance(image, nbObserv);
- covarianceMatrix = streamingCovariance(image, nbObserv);
- //covarianceMatrix = threadedCovariance(image, nbObserv);
- std::time_t end = time(NULL);
- vahineDebug("statistic calculation ended");
- vahineDebug("covariance calculation : "<< end - start << " seconds.");
-
- vahineDebug("covariance matrix\n"); //<<covarianceMatrix);
-
- unsigned int dim = covarianceMatrix.cols();
-
- VnlMatrix meanProduct(dim, dim);
- for(unsigned int j=0; j<dim; ++j) {
- VnlVector currentColumn(dim, m_Mean[j]);
- meanProduct.set_column(j, currentColumn);
- }
-
- VnlMatrix correlationMatrix(dim, dim);
- for(unsigned int i=0; i<dim; ++i) {
- meanProduct.set_row(i,meanProduct.get_row(i)*m_Mean[i]);
- }
- vahineDebug("meanProduct\n"); //<<meanProduct);
- correlationMatrix = covarianceMatrix + meanProduct;
- vahineDebug("correlation matrix\n"); //<<correlationMatrix);
-
- vnl_symmetric_eigensystem<double> eigenK(covarianceMatrix);
- VnlVector eigenCovariance = eigenK.D.diagonal();
- vcl_sort(eigenCovariance.begin(),eigenCovariance.end() );
- vnl_symmetric_eigensystem<double> eigenR(correlationMatrix);
- VnlVector eigenCorrelation = eigenR.D.diagonal();
- vcl_sort(eigenCorrelation.begin(),eigenCorrelation.end() );
- eigenCorrelation.flip();
- eigenCovariance.flip();
-
- std::cout << "eigen covariance = " << eigenCovariance << std::endl;
- std::cout << "eigen correlation = " << eigenCorrelation << std::endl;
- VnlVector likeHood;
- LikelihoodLog(eigenCorrelation, eigenCovariance, nbObserv, likeHood);
- std::cout << "likeHood = " << likeHood << std::endl;
- LocalMaximum(likeHood);
- for(unsigned int i=0; i<m_MaxLikeHood.size(); ++i) {
- if(m_MaxLikeHood[i]>=2)
- return m_MaxLikeHood[i];
- }
- // if maxlikehood is empty, normally never in real case
- // we search maximum value
- // this is only for test with a small image
- double max = likeHood.max_value();
- for(unsigned int i=0; i<likeHood.size(); ++i) {
- if(abs(likeHood[i] - max) < m_RealLimits.epsilon() ) {
- return i;
- }
- }
- // if we return 0 it's a bug !
- // TODO throw exception
- return 0;
-}
-
-/**
- * save all local maximum index
- * min that the index start to 0
- */
-template<class TImage>
-void VahineElmVcaFilter<TImage>::LocalMaximum(VnlVector & likeHood){
- unsigned int n = likeHood.size()-1;
- VnlVector df(n), fi(n), fii(n);
-
- fi.copy_in(likeHood.data_block() );
- fii.copy_in(likeHood.data_block()+1);
- df = fii-fi;
- for(unsigned int i=0; i<df.size()-1; ++i) {
- if(df(i)>=0 && df(i+1)<0) {
- m_MaxLikeHood.push_back(i+1);
- }
- }
-}
-
-/**
- * calculate the likelihodd function of eigen vector
- * and the logarithme of this function.
- */
-template<class TImage>
-void VahineElmVcaFilter<TImage>::LikelihoodLog(VnlVector & eigenR, VnlVector & eigenK, unsigned int NObserv,
- VnlVector & lH){
- unsigned int Ns = eigenR.size();
-
- lH.set_size(Ns);
- double coef = 2.0/NObserv;
- for(unsigned int i=0; i<Ns; ++i) {
- unsigned int Nl = Ns - i;
- VnlVector r(Nl), k(Nl), ZSquare(Nl), Sigma(Nl), T(Nl);
- r.copy_in(&eigenR.data_block()[i]);
- k.copy_in(&eigenK.data_block()[i]);
- ZSquare = (r-k).apply(&mysquare);
- Sigma = r.apply(&mysquare) + k.apply(&mysquare);
- Sigma *= coef;
- VnlVector::iterator itSigma = Sigma.begin();
- VnlVector::iterator itZSquare = ZSquare.begin();
- VnlVector::iterator itT = T.begin();
- while(itSigma != Sigma.end() ) {
- *itT = (*itZSquare)/(*itSigma);
- itSigma++;
- itZSquare++;
- itT++;
- }
- Sigma=Sigma.apply(&log);
- lH(i) = -0.5*T.sum() - 0.5*Sigma.sum();
- }
- std::cout << "likelihood " << lH << std::endl;
-}
-
-template<class TImage>
-typename VahineElmVcaFilter<TImage>::IntegerList VahineElmVcaFilter<TImage>::GetMaxLikeHood(){
- return m_MaxLikeHood;
-}
-
-/**
- * Standard "PrintSelf" method
- */
-template<class TImage>
-void VahineElmVcaFilter<TImage>::PrintSelf(std::ostream& os, itk::Indent indent) const {
- Superclass::PrintSelf( os, indent );
-}
-
-template<class TImage>
-void VahineElmVcaFilter<TImage>::GenerateOutputInformation(){
- Superclass::GenerateOutputInformation();
- // get pointers to the input and output
- VectorImageConstPointer inputPtr = this->GetInput();
- VectorImagePointer outputPtr = this->GetOutput();
-
- if ( !inputPtr || !outputPtr ) {
- return;
- }
- // we need to compute the output spacing, the output image size, and the
- // output image start index
- const SizeType& inputSize = inputPtr->GetLargestPossibleRegion().GetSize();
- RegionType outputLargestPossibleRegion;
- outputLargestPossibleRegion.SetSize( inputSize );
- outputPtr->SetLargestPossibleRegion( outputLargestPossibleRegion );
- UpdateNumberOfComponents();
- outputPtr->SetNumberOfComponentsPerPixel(m_NbComponents);
-}
-
-} // end namespace otb
-
-#endif
diff --git a/Code/Vahine/FlatReshapeFilter.h b/Code/Vahine/FlatReshapeFilter.h
deleted file mode 100644
index 1eae46be1b59bf093cf6b2bc48906a4506d9328f..0000000000000000000000000000000000000000
--- a/Code/Vahine/FlatReshapeFilter.h
+++ /dev/null
@@ -1,115 +0,0 @@
-/*
- * FlatReshapeFilter.h
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINEFLATRESHAPEFILTER_H_
-#define __VAHINEFLATRESHAPEFILTER_H_
-
-#include <otb/IO/otbImage.h>
-#include <otb/IO/otbVectorImage.h>
-#include <itkImageToImageFilter.h>
-#include <itkImageRegionIterator.h>
-#include <itkImageLinearConstIteratorWithIndex.h>
-#include <itkImageLinearIteratorWithIndex.h>
-#include <itkImageRegionConstIteratorWithIndex.h>
-#include <itkExceptionObject.h>
-
-#include "VahineMacro.h"
-
-namespace otb {
-/**
- * \class VahineFlatReshapeFilter
- * \brief Filter to apply a mask to a hyperspectral image.
- */
-template <class TBaseImage, class TMaskImage>
-class ITK_EXPORT VahineFlatReshapeFilter
- : public itk::ImageToImageFilter<TBaseImage, TBaseImage> {
-public:
- typedef TBaseImage BaseImageType;
- typedef TMaskImage MaskImageType;
-
- typedef VahineFlatReshapeFilter Self;
- typedef itk::ImageToImageFilter<BaseImageType, BaseImageType> Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
- typedef typename MaskImageType::Pointer MaskImagePointer;
-
- // creation of SmartPointer
- itkNewMacro(Self);
- // runtime type information
- itkTypeMacro(VahineFlatReshapeFilter, itk::ImageToImageFilter);
-
- typedef typename BaseImageType::Pointer BaseImagePointer;
- typedef typename BaseImageType::ConstPointer BaseImageConstPointer;
- typedef typename BaseImageType::IndexType BaseImageIndex;
- typedef typename BaseImageType::SizeType BaseImageSize;
- typedef typename BaseImageType::RegionType BaseImageRegion;
- typedef typename BaseImageType::PixelType PixelType;
- typedef typename BaseImageType::InternalPixelType InternalPixelType;
- typedef typename MaskImageType::RegionType MaskImageRegion;
- typedef typename MaskImageType::IndexType MaskImageIndex;
- typedef typename MaskImageType::SizeType MaskImageSize;
-
- typedef itk::ImageRegionConstIterator<BaseImageType> ConstBaseIteratorType;
- typedef itk::ImageRegionConstIterator<MaskImageType> ConstMaskIteratorType;
- typedef itk::ImageRegionIterator<BaseImageType> OutputIteratorType;
-
- typedef itk::ImageRegionConstIteratorWithIndex<MaskImageType> ConstMaskIteratorIndexType;
-
- itkSetMacro( Mask, MaskImagePointer );
- itkSetMacro( HighThreshold, double );
- itkGetMacro( HighThreshold, double );
- itkSetMacro( LowThreshold, double );
- itkGetMacro( LowThreshold, double );
- itkSetMacro( FlatReshape, bool );
- virtual void GenerateOutputInformation();
-
- virtual void GenerateInputRequestedRegion();
-
- bool GetDebug(){
- return true;
- };
-protected:
- VahineFlatReshapeFilter();
- virtual ~VahineFlatReshapeFilter();
- void PrintSelf(std::ostream& os, itk::Indent indent) const;
-
- void GenerateData();
-
- //void ThreadedGenerateData(const BaseImageRegion& outputRegionForThread, int
- // threadId );
-private:
- VahineFlatReshapeFilter(const Self&); //not implemented
- void operator=(const Self&); //not implemented
-
- //void GetMaskRegion(BaseImageRegion outputRegion, MaskImageIndex &
- // maskImageStartIndex, MaskImageRegion & maskImageRegion);
- unsigned int GetNumberOfActivPixels();
-
- void CheckInputBaseAndMaskSize();
-
- MaskImagePointer m_Mask;
- MaskImageRegion m_MaskLargestPossibleRegion;
- double m_HighThreshold;
- double m_LowThreshold;
- BaseImageRegion m_OutRegion;
- bool m_FlatReshape;
- InternalPixelType m_NullValue;
-
-};
-
-} // end namespace otb
-
-#include "FlatReshapeFilter.txx"
-
-#endif /* FLATRESHAPEFILTER_H_ */
\ No newline at end of file
diff --git a/Code/Vahine/FlatReshapeFilter.txx b/Code/Vahine/FlatReshapeFilter.txx
deleted file mode 100644
index 138fd995cf6d3ac03d9de67d7959f9d8e4ba9fc1..0000000000000000000000000000000000000000
--- a/Code/Vahine/FlatReshapeFilter.txx
+++ /dev/null
@@ -1,280 +0,0 @@
-/*
- * FlatReshapeFilter.txx
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-#ifndef __VAHINEFLATRESHAPEFILTER_TXX
-#define __VAHINEFLATRESHAPEFILTER_TXX
-#include <itkProgressReporter.h>
-#include "FlatReshapeFilter.h"
-
-namespace otb {
-
-template<class TBaseImage, class TMaskImage>
-VahineFlatReshapeFilter<TBaseImage, TMaskImage>::VahineFlatReshapeFilter() :
- m_LowThreshold(1.0), m_HighThreshold(1.0),
- m_FlatReshape(true),
- m_NullValue(0.0){
-
-}
-
-template<class TBaseImage, class TMaskImage>
-VahineFlatReshapeFilter<TBaseImage, TMaskImage>::~VahineFlatReshapeFilter(){
-}
-
-template<class TBaseImage, class TMaskImage>
-void VahineFlatReshapeFilter<TBaseImage, TMaskImage>::CheckInputBaseAndMaskSize(){
- vahineDebug("CheckInputBaseAndMaskSize");
- m_Mask->UpdateOutputInformation();
- BaseImageRegion imageRegion = this->GetInput()->GetLargestPossibleRegion();
- m_MaskLargestPossibleRegion = m_Mask->GetLargestPossibleRegion();
- if(m_FlatReshape) {
- if(m_MaskLargestPossibleRegion != imageRegion) {
- itkExceptionMacro("The size of base image and mask image is different.");
- }
- } else {
- // unflat
- vahineDebug("unflat");
- unsigned int nbPixelMask = GetNumberOfActivPixels();
- BaseImageSize imageSize = imageRegion.GetSize();
- unsigned int nbPixelImage = imageSize[0]*imageSize[1];
- if(nbPixelMask != nbPixelImage) {
- itkExceptionMacro("The number of pixel of base image and the number of valid pixel in mask image is different.");
- }
- }
-}
-
-/*
-template<class TBaseImage, class TMaskImage>
- void VahineFlatReshapeFilter<TBaseImage, TMaskImage>::GetMaskRegion(BaseImageRegion outputRegion, MaskImageIndex & maskImageStartIndex, MaskImageRegion & maskImageRegion){
- vahineDebug("GetMaskRegion");
- ConstMaskIteratorIndexType itMask(m_Mask, m_MaskLargestPossibleRegion);
- itMask.GoToBegin();
- //MaskImageRegion maskRegion;
- MaskImageIndex currentIndex, maskStartIndex;
- MaskImageSize maskSize;
- unsigned int count = 0;
- while(count<outputRegion.GetIndex()[1]+1){
- if(static_cast<double>(itMask.Get()) > m_Threshold){
- ++count;
- maskStartIndex = itMask.GetIndex();
- }
- ++itMask;
- }
- while(count<outputRegion.GetSize()[1]){
- if(static_cast<double>(itMask.Get()) > m_Threshold){
- ++count;
- currentIndex = itMask.GetIndex();
- }
- ++itMask;
- }
- maskImageStartIndex = maskStartIndex;
- maskStartIndex[0] = 0; // force first column
- maskSize[0] = m_MaskLargestPossibleRegion.GetSize()[0]; // force width of image
- maskSize[1] = currentIndex[1] - maskStartIndex[1] + 1;
- maskImageRegion.SetIndex( maskStartIndex );
- maskImageRegion.SetSize( maskSize );
-}
-*/
-
-/**
- * Get number of activ pixel
- * We need to update mask image for accessing to mask value
- */
-template<class TBaseImage, class TMaskImage>
-unsigned int VahineFlatReshapeFilter<TBaseImage, TMaskImage>::GetNumberOfActivPixels(){
- vahineDebug("GetNumberOfActivPixels");
- m_Mask->Update();
- ConstMaskIteratorType itMask(m_Mask, m_Mask->GetLargestPossibleRegion() );
- itMask.GoToBegin();
- unsigned int count = 0;
- unsigned int nbpixel = 0;
- while(!itMask.IsAtEnd() ) {
- double maskValue = static_cast<double>(itMask.Get() );
- if( m_LowThreshold <= maskValue && maskValue <= m_HighThreshold ) {
- ++count;
- }
- ++itMask;
- ++nbpixel;
- }
- vahineDebug("count number of activ pixel = "<<count);
- vahineDebug("low threshold ="<<m_LowThreshold<<" high threshold ="<<m_HighThreshold);
- return count;
-}
-
-/**
- * Output image have one columns
- * We update ouput size image
- */
-template<class TBaseImage, class TMaskImage>
-void VahineFlatReshapeFilter<TBaseImage, TMaskImage>::GenerateOutputInformation(){
- vahineDebug("GenerateOutputInformation");
- // call the superclass' implementation of this method
- Superclass::GenerateOutputInformation();
-
- // get pointers to the input and output
- BaseImageConstPointer inputPtr = this->GetInput();
- BaseImagePointer outputPtr = this->GetOutput();
-
- if ( !inputPtr || !outputPtr ) {
- return;
- }
-
- CheckInputBaseAndMaskSize();
-
- BaseImageIndex outputStartIndex;
- outputStartIndex.Fill(0);
- BaseImageSize outputSize;
- if( m_FlatReshape ) {
- // size along X
- outputSize[0] = 1;
- outputSize[1] = GetNumberOfActivPixels();
- } else {
- // unflat
- outputSize = m_Mask->GetLargestPossibleRegion().GetSize();
- }
- m_OutRegion.SetSize( outputSize );
- m_OutRegion.SetIndex( outputStartIndex );
- outputPtr->SetLargestPossibleRegion( m_OutRegion );
- outputPtr->SetNumberOfComponentsPerPixel(inputPtr->GetNumberOfComponentsPerPixel() );
-}
-
-template<class TBaseImage, class TMaskImage>
-void VahineFlatReshapeFilter<TBaseImage, TMaskImage>::GenerateInputRequestedRegion(){
- vahineDebug("GenerateInputRequestedRegion");
- // call the superclass' implementation of this method
- Superclass::GenerateInputRequestedRegion();
-
- // get pointers to the input and output
- BaseImagePointer inputPtr = const_cast<BaseImageType *>(this->GetInput() );
- BaseImagePointer outputPtr = this->GetOutput();
-
- if ( !inputPtr || !outputPtr ) {
- return;
- }
-
- // we need to compute the input requested region (size and start index)
-
- BaseImageRegion inputRequestedRegion = inputPtr->GetLargestPossibleRegion();
- /*
- const BaseImageRegion & outputRequestedRegion = outputPtr->GetRequestedRegion();
- BaseImageRegion inputRequestedRegion = inputPtr->GetRequestedRegion();
- vahineDebug("output request "<<outputRequestedRegion);
- BaseImageIndex dummy;
- GetMaskRegion(outputRequestedRegion, dummy, inputRequestedRegion);
-
- inputRequestedRegion.Crop( inputPtr->GetLargestPossibleRegion() );
- */
- inputPtr->SetRequestedRegion( inputRequestedRegion );
- vahineDebug("input request "<<inputRequestedRegion);
-}
-
-/**
- * Each activ paixel is an row with one pixel wise in new image
- * z-size (number of band) is unchanged
- */
-template<class TBaseImage, class TMaskImage>
-void VahineFlatReshapeFilter<TBaseImage, TMaskImage>::GenerateData(){
- CheckInputBaseAndMaskSize();
- MaskImageRegion maskRegion = m_Mask->GetLargestPossibleRegion();
- BaseImageRegion baseRegion = this->GetInput()->GetLargestPossibleRegion();
- vahineDebug("base region "<<baseRegion<<" mask region "<<maskRegion);
- vahineDebug("out region "<<m_OutRegion);
- ConstBaseIteratorType itBaseImage(this->GetInput(), baseRegion);
- ConstMaskIteratorType itMask(m_Mask, maskRegion);
-
- BaseImagePointer outputImage = BaseImageType::New();
- outputImage->SetRegions(m_OutRegion);
- unsigned int componentsPerPixel = this->GetInput()->GetNumberOfComponentsPerPixel();
- vahineDebug("components per pixel = "<<componentsPerPixel);
- outputImage->SetNumberOfComponentsPerPixel(componentsPerPixel);
- outputImage->Allocate();
- OutputIteratorType itOutput(outputImage, m_OutRegion);
- itOutput.GoToBegin();
- if( m_FlatReshape ) {
- for(itBaseImage.GoToBegin(), itMask.GoToBegin(); !itBaseImage.IsAtEnd(); ++itBaseImage, ++itMask) {
- double maskValue = static_cast<double>(itMask.Get() );
- if(m_LowThreshold <= maskValue && maskValue <= m_HighThreshold ) {
- itOutput.Set(itBaseImage.Get() );
- ++itOutput;
- }
- }
- } else {
- // unflat operation
- itBaseImage.GoToBegin();
- for(itMask.GoToBegin(), itOutput.GoToBegin(); !itMask.IsAtEnd(); ++itMask, ++itOutput) {
- double maskValue = static_cast<double>(itMask.Get() );
- if(m_LowThreshold <= maskValue && maskValue <= m_HighThreshold ) {
- itOutput.Set(itBaseImage.Get() );
- ++itBaseImage;
- } else {
- PixelType nullPixel;
- nullPixel.Reserve(componentsPerPixel);
- nullPixel.Fill(m_NullValue);
- itOutput.Set(nullPixel);
- }
- }
- }
- this->GraftOutput(outputImage);
-}
-
-/**
- * Each activ pixel is an row with one pixel wise in new image
- * z-size (number of band) is unchanged
- */
-/*
-template<class TBaseImage, class TMaskImage>
- void VahineFlatReshapeFilter<TBaseImage, TMaskImage>::ThreadedGenerateData(const BaseImageRegion& outputRegionForThread, int threadId ){
- vahineDebug("ThreadedGenerateData");
- // Get the input and output pointers
- BaseImageConstPointer inputPtr = this->GetInput();
- BaseImagePointer outputPtr = this->GetOutput();
-
- // Define/declare an iterator that will walk the output region for this
- // thread.
- OutputIteratorType outIt(outputPtr, outputRegionForThread);
-
- // support progress methods/callbacks
- itk::ProgressReporter progress(this, threadId, outputRegionForThread.GetNumberOfPixels());
-
- MaskImageRegion maskImageRegion = MaskImageRegion();
- MaskImageIndex maskStartIndex = MaskImageIndex();
- GetMaskRegion(outputRegionForThread, maskStartIndex, maskImageRegion);
- ConstMaskIteratorIndexType maskIt(m_Mask, maskImageRegion);
-
- maskIt.SetIndex(maskStartIndex);
-
- // walk the output region, and sample the input image
- outIt.GoToBegin();
- while ( !outIt.IsAtEnd() ) {
- // copy the input pixel to the output
- vahineDebug("maskIt threadId "<<threadId<<" index "<<maskIt.GetIndex());
- while(static_cast<double>(m_Mask->GetPixel(maskIt.GetIndex()))<=m_Threshold){
- ++maskIt;
- }
- outIt.Set(inputPtr->GetPixel(maskIt.GetIndex()));
- ++outIt;
- ++maskIt;
- progress.CompletedPixel();
- }
-}
-*/
-/**
- * Standard "PrintSelf" method
- */
-template<class TBaseImage, class TMaskImage>
-void VahineFlatReshapeFilter<TBaseImage, TMaskImage>::PrintSelf(std::ostream& os, itk::Indent indent) const {
- Superclass::PrintSelf( os, indent );
-}
-
-} // end namespace otb
-
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/GrsirFilter.h b/Code/Vahine/GrsirFilter.h
deleted file mode 100644
index ccf1890913bbf5c35a8b25fdcdab8c38397c8674..0000000000000000000000000000000000000000
--- a/Code/Vahine/GrsirFilter.h
+++ /dev/null
@@ -1,269 +0,0 @@
-/*
- * GrsirFilter.h
- *
- * Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINEGRSIRFILTER_H
-#define __VAHINEGRSIRFILTER_H
-
-#include "VahineMacro.h"
-
-#include <utility>
-#include <set>
-#include <list>
-#include <vnl/vnl_vector_ref.h>
-#include <vnl/algo/vnl_real_eigensystem.h>
-#include <boost/lexical_cast.hpp>
-
-#include "gdal.h"
-#include "otbImage.h"
-
-#include "otbPerBandVectorImageFilter.h"
-#include "itkMeanImageFilter.h"
-#include "CenteringImageFilter.h"
-#include "otbPerBandVectorImageFilter.h"
-#include "itkMatrix.h"
-#include "itkArray.h"
-#include "itkVariableSizeMatrix.h"
-
-#include "VahineInterpol.h"
-#include "VahineSlice.h"
-#include "TikonovFilter.h"
-#include "otbStreamingStatisticsImageFilter.h"
-#include "otbMultiChannelExtractROI.h"
-#include "otbImageListToVectorImageFilter.h"
-#include "otbMultiToMonoChannelExtractROI.h"
-#include "InverseFilter.h"
-#include "VahineXml.h"
-
-namespace otb {
-
-class GrsirException { };
-
-// functor predicat for order our custom pair
-// don't use the internal pair operator< because it use the two elements for
-// comparaison
-// and we don't have operator< in the second type
-template < typename PairType >
-class SortPairIncrease {
-public:
- bool operator()( const PairType & left , const PairType & right ){
- return left.first < right.first;
- }
-
-};
-
-template < typename PairType >
-class SortPairDecrease {
-public:
- bool operator()( const PairType & left , const PairType & right ){
- return left.first > right.first;
- }
-
-};
-
-class ManagedTypes {
-public:
- // Convert a VariableLengthVector to an Array.
- template < class TElement >
- static itk::Array<TElement> ToManagedVariableLengthVector( itk::VariableLengthVector<TElement> vector ){
- typedef typename itk::Array<TElement> ArrayType;
- ArrayType result;
- result.SetSize(vector.Size() );
- for (unsigned int i=0; i<vector.Size(); ++i)
- result[i] = static_cast<TElement>(vector[i]);
- return result;
- }
-
-};
-
-template <class TInputImage, class TOutputImage>
-class ITK_EXPORT VahineGrsirFilter
- : public itk::ImageToImageFilter<TInputImage, TOutputImage> {
-public:
-
- typedef VahineGrsirFilter Self;
- typedef itk::ImageToImageFilter<TInputImage, TOutputImage> Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
-
- // creation of SmartPointer
- itkNewMacro(Self);
- // runtime type information
- itkTypeMacro(VahineGrsirFilter, itk::ImageToImageFilter);
-
- // Type to use for computation
- const static unsigned int Dimension = 2;
- typedef TInputImage InputImageType;
- typedef typename InputImageType::Pointer InputImagePointer;
- typedef typename InputImageType::PixelType PixelType;
- typedef typename InputImageType::InternalPixelType InternalPixelType;
- typedef typename InputImageType::IndexType IndexType;
- typedef typename itk::NumericTraits<InternalPixelType>::RealType RealType;
- typedef typename itk::VariableSizeMatrix<RealType> MatrixType;
-
- typedef Image<InternalPixelType, Dimension> BandType;
- typedef typename BandType::PixelType BandPixelType;
- typedef typename BandType::IndexType BandIndexType;
- typedef typename BandType::Pointer BandPointer;
- typedef itk::ImageRegionIterator<BandType> BandIterator;
- typedef ImageList<BandType> ImageBandListType;
- typedef typename ImageBandListType::Pointer ImageBandListPointer;
-
- typedef InputImageType VectorImageType;
- typedef typename VectorImageType::Pointer VectorImagePointer;
- typedef itk::ImageRegionConstIterator<VectorImageType> ConstVectorIteratorType;
-
- typedef vnl_matrix<RealType> VnlMatrixType;
- typedef vnl_vector<RealType> VnlVectorType;
-
- typedef std::pair<double, VnlVectorType> EigenPairType;
- typedef std::list<EigenPairType> EigenListType;
-
- typedef std::vector<double> RegulVectorType;
- typedef RegulVectorType* RegulVectorPointer;
-
- typedef Image<InternalPixelType, 2> ProjType;
- typedef typename ProjType::Pointer ProjPointer;
- typedef itk::ImageRegionIterator<ProjType> IteratorProjType;
-
- // used internal filter
- typedef otb::StreamingStatisticsImageFilter<BandType> StatsFilterType;
- typedef otb::VahineTikonovFilter<VectorImageType> TikonovFilterType;
-
- typedef typename std::list<vahine::Slice<InternalPixelType> > SliceListType;
- typedef typename SliceListType::iterator IteratorSliceListType;
-
- typedef XmlLog* XmlLogPointer;
-
- // members accessor
- itkGetMacro( RegulParams, RegulVectorPointer);
- itkSetMacro( RegulParams, RegulVectorPointer);
- itkGetMacro( NSlice, int);
- itkSetMacro( NSlice, int);
-
- itkSetMacro( Noise, VectorImagePointer );
- itkSetMacro( Learning, VectorImagePointer );
- itkSetMacro( LearningParam, VectorImagePointer );
-
- itkGetMacro( Epsilon, float);
- itkSetMacro( Epsilon, float);
-
- itkSetMacro( SearchBestRegulParam, bool );
-
- void SetTikonovFilter();
-
- void SetXmlLog(XmlLogPointer xmlLog);
-
- void PushBackLearningBand(unsigned int);
-
- bool GetDebug(){
- return true;
- };
-
- void PrintSelf(std::ostream& os, itk::Indent indent) const {
- Superclass::PrintSelf(os, indent);
- };
-
- SliceListType &GetSliceList();
-
-protected:
- VahineGrsirFilter();
- typedef itk::CenteringImageFilter<BandType, BandType> CenteringFilterType;
- typedef PerBandVectorImageFilter<VectorImageType, VectorImageType, CenteringFilterType> PerBandVectorImageFilterType;
- void GenerateData();
-
-private:
-
- /*** INTERNAL TYPE ***/
-
- // ParameterPairType :
- // * first element is Y value (float)
- // * second element is its index
- typedef std::pair<RealType, BandIndexType> ParameterPairType;
- typedef std::list<ParameterPairType> ParameterListType;
- typedef typename ParameterListType::const_iterator ParameterListIteratorType;
- typedef std::vector<unsigned int> BandNumberListType;
- typedef BandNumberListType::iterator BandNumberListIterator;
- typedef itk::ImageRegionConstIterator<BandType> ConstIteratorType;
- typedef MultiToMonoChannelExtractROI<BandPixelType, BandPixelType> ExtractChannelType;
- typedef ImageListToVectorImageFilter<ImageBandListType, VectorImageType> ImageListToVectorImageFilterType;
- typedef typename ImageListToVectorImageFilterType::Pointer ImageListToVectorImageFilterPointer;
-
- typedef VahineInverseFilter<InputImageType, BandType> VahineInverseType;
- typedef typename VahineInverseType::Pointer VahineInversePointer;
-
- itkSetMacro( Parameters, BandPointer );
-
- struct EigenStruct {
- double nrmse;
- VnlVectorType eigenVector;
- VahineInterpol<RealType> * interpolation;
- };
-
- /*** END INTERNAL TYPE ***/
-
- RegulVectorPointer m_RegulParams;
- MatrixType m_RegularizationMatrix;
- MatrixType m_CovarianceMatrix;
- int m_NSlice;
- BandPointer m_Parameters;
- VectorImagePointer m_Learning;
- VectorImagePointer m_Noise;
- VectorImagePointer m_LearningParam;
- bool m_SearchBestRegulParam;
-
- float m_Epsilon;
-
- SliceListType m_SliceList;
- BandNumberListType m_LearningParamBand;
-
- typename CenteringFilterType::Pointer m_CenteringFilter;
- typename PerBandVectorImageFilterType::Pointer m_CenteringVectorFilter;
- typename TikonovFilterType::Pointer m_TikonovFilter;
- BandPointer m_currentImageInv;
- ImageBandListPointer m_ImageBandList;
- VectorImagePointer m_LearningCenteringImage;
-
- VectorImagePointer m_OutputImage;
-
- XmlLogPointer m_XmlLog;
- //VahineInversePointer m_Inverse;
-
- VahineGrsirFilter(const Self&); //not implemented
- void operator=(const Self&); //not implemented
-
- void Slicing(ParameterListType parameterList, VectorImagePointer image);
-
- void MakeMean(unsigned int nObserv);
-
- void UpdateTikonovFilter(double regulParam);
-
- void EstimateEDR(VnlMatrixType &gamma, unsigned int numberOfBand, EigenStruct &eigenStruct);
-
- double CalculateSIRC(VnlMatrixType &gamma, VnlVectorType &axe);
-
- void Inverse(unsigned int numberOfBand, unsigned int nObserv);
-
- void InverseImage(EigenStruct eigenStruct, InputImagePointer image);
-
- void CalculateNRMSE(EigenStruct &eigenStruct, IteratorProjType &itLearningNoiseEstim);
-
- void CalculateEigenVector(VnlMatrixType &gamma, EigenStruct &eigenStruct);
-
-};
-
-} // End namespace otb
-
-#include "GrsirFilter.txx"
-
-#endif /* GRSIRFILTER_H_ */
\ No newline at end of file
diff --git a/Code/Vahine/GrsirFilter.txx b/Code/Vahine/GrsirFilter.txx
deleted file mode 100644
index 252855b14beb50d8a3d73848c031a316d38853a0..0000000000000000000000000000000000000000
--- a/Code/Vahine/GrsirFilter.txx
+++ /dev/null
@@ -1,487 +0,0 @@
-/**
- * Grsir.txx
- *
- * Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINEGRSIRFILTER_TXX
-#define __VAHINEGRSIRFILTER_TXX
-
-#include "GrsirFilter.h"
-
-namespace otb {
-
-template<class TInputImage, class TOutputImage>
-VahineGrsirFilter<TInputImage, TOutputImage>::VahineGrsirFilter(){
- m_Epsilon = 0.000001;
- m_CenteringFilter = CenteringFilterType::New();
- m_CenteringVectorFilter = PerBandVectorImageFilterType::New();
- m_CenteringVectorFilter->SetFilter(m_CenteringFilter);
- m_ImageBandList = ImageBandListType::New();
- m_SearchBestRegulParam = true;
- //m_ImageBandList->DisconnectPipeline();
- //m_Inverse = VahineInverseType::New();
-}
-
-template<class TInputImage, class TOutputImage>
-void VahineGrsirFilter<TInputImage, TOutputImage>::Slicing(ParameterListType parameterList, VectorImagePointer image){
- unsigned int nObserv = image->GetLargestPossibleRegion().GetNumberOfPixels();
-
- // distribute value in several class
-
- vahineDebug("VahineGrsirFilter : distribute data into class param="<<parameterList.size()<<" observ="<<nObserv);
- ParameterListIteratorType itParameterList = parameterList.begin();
- ParameterListIteratorType itLastParameter = parameterList.begin();
- if(!m_SliceList.empty() ) {
- vahineDebug("**** ERROR ****");
- throw GrsirException();
- }
- int i = 0;
- while(itParameterList != parameterList.end() ) {
- i++;
- BandIndexType index = (*itParameterList).second;
- PixelType pixel = image->GetPixel(index);
- if(m_SliceList.empty() ) {
- // push into list the first slice
- m_SliceList.push_back(vahine::Slice<InternalPixelType>(index, pixel, itParameterList->first) );
- } else {
- if( (itParameterList->first - itLastParameter->first)<m_Epsilon) {
- //vahineDebug("value closer : same slice ... current =
- // "<<itParameterList->first<<" last = "<<itLastParameter->first);
- // sum current pixel value to the current slice
- try{
- m_SliceList.back().add(index, pixel, itParameterList->first);
- }catch(...) {
- vahineDebug("add element failed");
- }
- } else {
- //vahineDebug("current pixel is in next class");
- // the current pixel is in next class
- // create the last mean object
- MakeMean(nObserv);
- // create a new slice with current pixel
- //vahineDebug("create a new slice with current pixel");
- m_SliceList.push_back(vahine::Slice<InternalPixelType>(index, pixel, itParameterList->first) );
- }
- itLastParameter++;
- }
- itParameterList++;
- }
- if(!parameterList.empty() ) {
- // create the last mean object
- MakeMean(nObserv);
- }
-}
-
-template<class TInputImage, class TOutputImage>
-void VahineGrsirFilter<TInputImage, TOutputImage>::MakeMean(unsigned int nObserv){
- m_SliceList.back().setMeanPixel();
- m_SliceList.back().setPropOfPixel(nObserv);
-}
-
-template<class TInputImage, class TOutputImage>
-typename VahineGrsirFilter<TInputImage, TOutputImage>::SliceListType &
-VahineGrsirFilter<TInputImage, TOutputImage>::GetSliceList(){
- return m_SliceList;
-}
-
-template<class TInputImage, class TOutputImage>
-void VahineGrsirFilter<TInputImage, TOutputImage>::SetTikonovFilter(){
- m_TikonovFilter = TikonovFilterType::New();
- m_TikonovFilter->SetInput(m_Learning);
- m_TikonovFilter->UpdateCov();
- m_CovarianceMatrix = m_TikonovFilter->GetCovMatrix();
- VnlMatrixType m = m_CovarianceMatrix.GetVnlMatrix();
- for(int i=0; i<1; i++) {
- for(int j=0; j<10; j++) {
- vahineDebug("cov = "<<m.get(i,j) );
- }
- }
-}
-
-template<class TInputImage, class TOutputImage>
-void VahineGrsirFilter<TInputImage, TOutputImage>::SetXmlLog(XmlLogPointer xmlLog){
- m_XmlLog = xmlLog;
-}
-
-template<class TInputImage, class TOutputImage>
-void VahineGrsirFilter<TInputImage, TOutputImage>::UpdateTikonovFilter(double regulParam){
- m_TikonovFilter->SetRegularizationParam(regulParam);
- m_TikonovFilter->Update();
- m_RegularizationMatrix = m_TikonovFilter->GetTikonovMatrix();
- VnlMatrixType m = m_RegularizationMatrix.GetVnlMatrix();
- for(int i=0; i<1; i++) {
- for(int j=0; j<10; j++) {
- vahineDebug("tiko = "<<m.get(i,j) );
- }
- }
-}
-
-template<class TInputImage, class TOutputImage>
-void VahineGrsirFilter<TInputImage, TOutputImage>::GenerateData(){
- vahineDebug("VahineGrsirFilter : GenerateData");
- m_XmlLog->CreateElement("filter-grsir");
-
- // centering data
- m_CenteringVectorFilter->SetInput(m_Learning);
- m_CenteringVectorFilter->Update();
- m_LearningCenteringImage = m_CenteringVectorFilter->GetOutput();
- unsigned int numberOfBand = m_LearningCenteringImage->GetNumberOfComponentsPerPixel();
- unsigned int nObserv = m_LearningCenteringImage->GetLargestPossibleRegion().GetNumberOfPixels();
-
- // Extract the correct parameter band
- typename ExtractChannelType::Pointer extractParamBand = ExtractChannelType::New();
- extractParamBand->SetExtractionRegion(m_LearningParam->GetLargestPossibleRegion() );
- extractParamBand->SetInput(m_LearningParam);
-
- // loop on parameter
- BandNumberListIterator it = m_LearningParamBand.begin();
- while(it != m_LearningParamBand.end() ) {
- vahineDebug("Inversion for parameter : "<<(*it) );
- extractParamBand->SetChannel(*it);
- m_Parameters = extractParamBand->GetOutput();
- extractParamBand->Update();
-
- vahineDebug("Imagebandlist before size ="<<m_ImageBandList->Size() );
- //******************* Inverse call *********************/
- m_XmlLog->CreateElement("optimization");
- m_XmlLog->CreateAttribute("parameter",boost::lexical_cast<std::string>( (*it) ) );
- Inverse(numberOfBand, nObserv);
- m_XmlLog->AppendToSup();
- //***************************************************/
- m_ImageBandList->Update();
-
- //TODO cleanup internal variable ?
- ++it;
- //break; // for debug one parameter
- }
-
- // convert image band list to vector image
- ImageListToVectorImageFilterPointer imageBandList2VectorImage = ImageListToVectorImageFilterType::New();
- m_ImageBandList->Update();
- imageBandList2VectorImage->SetInput(m_ImageBandList);
- vahineDebug("image size = "<<m_ImageBandList->Size() );
- imageBandList2VectorImage->UpdateOutputInformation();
- imageBandList2VectorImage->Update();
- vahineDebug("get ouput");
- m_OutputImage = imageBandList2VectorImage->GetOutput();
-
- ConstVectorIteratorType itt(m_OutputImage, m_OutputImage->GetRequestedRegion() );
- // for debugging
- int h = 1;
- for(itt.GoToBegin(); !itt.IsAtEnd(); ++itt) {
- if(h>=120 && h<=128) {
- vahineDebug("OutputImage ="<<itt.Get() );
- }
- h++;
- }
-
- this->GraftOutput(m_OutputImage);
-}
-
-template<class TInputImage, class TOutputImage>
-void VahineGrsirFilter<TInputImage, TOutputImage>::PushBackLearningBand(unsigned int band){
- m_LearningParamBand.push_back(band);
-}
-
-/**
- * Sortering Y parameters
- * Calculate the best eigen vector and apply it to input image
- * for inversion
- */
-template<class TInputImage, class TOutputImage>
-void
-VahineGrsirFilter<TInputImage, TOutputImage>::Inverse(unsigned int numberOfBand, unsigned int nObserv){
-
- // Sortering Y !!
- vahineDebug("VahineGrsirFilter : Sortering Y");
- vahineDebug("Number of band learning = "<<numberOfBand);
- vahineDebug("Number of observation learning = "<<nObserv);
- m_Parameters->Update();
- vahineDebug("Number of parameter values = "<<m_Parameters->GetLargestPossibleRegion().GetNumberOfPixels() );
-
- ConstIteratorType itParameters( m_Parameters, m_Parameters->GetLargestPossibleRegion() );
- ParameterListType parameterList;
-
- for(itParameters.GoToBegin(); !itParameters.IsAtEnd(); ++itParameters) {
- BandIndexType index = itParameters.GetIndex();
- typename BandType::PixelType value = itParameters.Get();
- //vahineDebug("VahineGrsirFilter : index = "<<index<<" param v = "<< value);
- ParameterPairType element = make_pair(value, index);
- parameterList.push_back(element);
- }
-
- // sort pair by Y value
- vahineDebug("VahineGrsirFilter : Sort pair by Y value");
- parameterList.sort(SortPairIncrease<ParameterPairType>() );
-
- Slicing(parameterList, m_LearningCenteringImage);
-
- // calculate the between slice mean covariance matrix
- vahineDebug("VahineGrsirFilter : Calculate slice mean");
- vahineDebug("VahineGrsirFilter : number of slice = "<<m_SliceList.size() );
- vahineDebug("VahineGrsirFilter : Gamma size = "<<numberOfBand<<"*"<<numberOfBand);
- VnlMatrixType gamma = VnlMatrixType(numberOfBand, numberOfBand, 0);
- IteratorSliceListType itSlice = m_SliceList.begin();
- while(itSlice != m_SliceList.end() ) {
- VnlVectorType meanH = ManagedTypes::ToManagedVariableLengthVector<RealType>(itSlice->getMeanPixel() );
- //vahineDebug("VahineGrsirFilter : meanH = "<<meanH);
- //vahineDebug("propOfPixel = "<<itSlice->getPropOfPixel());
- gamma += itSlice->getPropOfPixel() * outer_product(meanH, meanH);
- ++itSlice;
- }
-
- // Grsir core : estimate GRSIR axes
- std::vector<double> nrmseVector;
- // eigenStructVEctor : used with no auto regularization option
- // all regularization and projection are used
- std::vector<EigenStruct> eigenStructVector;
- double regulMinNrmse;
- int regulIndex;
- EigenStruct bestEigenStruct;
- RegulVectorType::iterator itRegulParam = m_RegulParams->begin();
- while( itRegulParam != m_RegulParams->end() ) {
- // set regularization matrix
- double regul = *itRegulParam;
- UpdateTikonovFilter(regul);
- EigenStruct currentEigenStruct;
- // ********************** ESTIMATE EDR CALL**************************
- EstimateEDR(gamma, numberOfBand, currentEigenStruct);
- //vahineDebug("Regul = "<<regul<<" Nrmse = "<<currentEigenStruct.nrmse);
- nrmseVector.push_back(currentEigenStruct.nrmse);
- if(!m_SearchBestRegulParam) {
- eigenStructVector.push_back(currentEigenStruct);
- }
- if(itRegulParam == m_RegulParams->begin() ) {
- // first
- bestEigenStruct = currentEigenStruct;
- regulMinNrmse = regul;
- regulIndex = itRegulParam-m_RegulParams->begin();
- } else {
- if(bestEigenStruct.nrmse>currentEigenStruct.nrmse) {
- bestEigenStruct = currentEigenStruct;
- regulMinNrmse = regul;
- regulIndex = itRegulParam-m_RegulParams->begin();
- }
- }
- ++itRegulParam;
- }
- std::vector<double>::iterator it = nrmseVector.begin();
- while(it != nrmseVector.end() ) {
- vahineDebug("all nrmse = "<<*it);
- ++it;
- }
-
- vahineDebug("min nrmse is "<<bestEigenStruct.nrmse<<" regul value is "<<regulMinNrmse);
- //vahineDebug("gamma "<<gamma);
- double cgrsir = CalculateSIRC(gamma, bestEigenStruct.eigenVector);
- vahineDebug("cgrsir = "<<cgrsir);
- m_XmlLog->CreateElementTextApp("nrmse", boost::lexical_cast<std::string>(bestEigenStruct.nrmse) );
- m_XmlLog->CreateElement("regularisation");
- m_XmlLog->CreateAttribute("index",boost::lexical_cast<std::string>(regulIndex) );
- m_XmlLog->CreateTextNode(boost::lexical_cast<std::string>(regulMinNrmse) );
- m_XmlLog->AppendToSup();
- m_XmlLog->CreateElementTextApp("cgrsir",boost::lexical_cast<std::string>(cgrsir) );
-
- // ************************************************************************
- // *********************** inverse image
- // **************************
- // ************************************************************************
- InputImagePointer image = const_cast<InputImageType *>(this->GetInput() );
-
- if(m_SearchBestRegulParam) {
- vahineDebug("InverseImage for best regularization parameter");
- InverseImage(bestEigenStruct, image);
- } else {
- vahineDebug("InverseImage for all regularization parameter");
- typename std::vector<EigenStruct>::iterator it = eigenStructVector.begin();
- while(it != eigenStructVector.end() ) {
- vahineDebug("inverse...");
- InverseImage(*it, image);
- ++it;
- }
- }
- // clear list of slice for next parameter @see Slicing
- m_SliceList.clear();
-}
-
-template<class TInputImage, class TOutputImage>
-void
-VahineGrsirFilter<TInputImage, TOutputImage>::InverseImage(EigenStruct eigenStruct, InputImagePointer image){
- VahineInversePointer m_Inverse = VahineInverseType::New();
-
- m_Inverse->SetInput(image);
- //m_Inverse->SetReleaseDataFlag(true);
- //m_Inverse->SetReleaseDataBeforeUpdateFlag(true);
- m_Inverse->SetGrsirAxe(eigenStruct.eigenVector);
- m_Inverse->SetInterpolation(eigenStruct.interpolation);
- m_Inverse->UpdateOutputInformation();
- m_ImageBandList->PushBack(m_Inverse->GetOutput() );
- m_Inverse->Update();
- m_ImageBandList->Update();
-}
-
-template<class TInputImage, class TOutputImage>
-void
-VahineGrsirFilter<TInputImage, TOutputImage>::EstimateEDR(VnlMatrixType &gamma, unsigned int numberOfBand,
- EigenStruct &eigenStruct){
- // calculate eigen vector
- CalculateEigenVector(gamma, eigenStruct);
-
- // project learning value to first eigen vector
- // project learning noise value to first eigen vector
- vahineDebug("VahineGrsirFilter : projection of learning data");
- ConstVectorIteratorType itLearningData(m_Learning, m_Learning->GetLargestPossibleRegion() );
- ConstVectorIteratorType itLearningNoise(m_Noise, m_Noise->GetLargestPossibleRegion() );
-
- ProjPointer learningDataProj = ProjType::New();
- learningDataProj->SetRegions( m_Learning->GetLargestPossibleRegion() );
- learningDataProj->Allocate();
- IteratorProjType itLearningDataProj(learningDataProj, learningDataProj->GetLargestPossibleRegion() );
-
- ProjPointer learningNoiseProj = ProjType::New();
- learningNoiseProj->SetRegions( m_Noise->GetLargestPossibleRegion() );
- learningNoiseProj->Allocate();
- IteratorProjType itLearningNoiseProj(learningNoiseProj, learningNoiseProj->GetLargestPossibleRegion() );
-
- for(itLearningData.GoToBegin(), itLearningDataProj.GoToBegin(); !itLearningData.IsAtEnd();
- ++itLearningData, ++itLearningDataProj) {
- PixelType currentPixel = itLearningData.Get();
- VnlVectorType pixel = ManagedTypes::ToManagedVariableLengthVector<RealType>(currentPixel);
- RealType tmp = static_cast<RealType>(dot_product(pixel, eigenStruct.eigenVector) );
- itLearningDataProj.Set(tmp);
- }
-
- for(itLearningNoise.GoToBegin(), itLearningNoiseProj.GoToBegin(); !itLearningNoise.IsAtEnd();
- ++itLearningNoise, ++itLearningNoiseProj) {
- PixelType currentPixel = itLearningNoise.Get();
- VnlVectorType pixel = ManagedTypes::ToManagedVariableLengthVector<RealType>(currentPixel);
- RealType tmp = static_cast<RealType>(dot_product(pixel, eigenStruct.eigenVector) );
- itLearningNoiseProj.Set(tmp);
- }
-
- // calculate the mean values of the current functionnal
- IteratorSliceListType itSlice = m_SliceList.begin();
- unsigned int sizeOfInterpol = m_SliceList.size();
- while( itSlice != m_SliceList.end() ) {
- itSlice->calculateMean(learningDataProj);
- ++itSlice;
- }
- m_SliceList.sort();
-
- vahineDebug("get mean and param sizeOfInterpol = "<<sizeOfInterpol);
- VnlVectorType dataVector = VnlVectorType(sizeOfInterpol);
- VnlVectorType paramVector = VnlVectorType(sizeOfInterpol);
- unsigned int i = 0;
- itSlice = m_SliceList.begin();
- while( itSlice != m_SliceList.end() ) {
- dataVector[i] = (*itSlice).getMeanData();
- paramVector[i] = (*itSlice).getMeanParam();
- //vahineDebug("data "<<dataVector[i]<<" param "<<paramVector[i]);
- ++itSlice;
- i++;
- }
-
- // interpolation for learning noise
- vahineDebug("interpolation for learning noise data");
- ProjPointer learningNoiseEstim = ProjType::New();
- learningNoiseEstim->SetRegions( m_Learning->GetLargestPossibleRegion() );
- learningNoiseEstim->Allocate();
- IteratorProjType itLearningNoiseEstim(learningNoiseEstim, learningNoiseEstim->GetLargestPossibleRegion() );
-
- eigenStruct.interpolation = new VahineInterpol<RealType>(dataVector.size() );
- eigenStruct.interpolation->SetVx(dataVector);
- eigenStruct.interpolation->SetVy(paramVector);
- eigenStruct.interpolation->Update();
- vahineDebug("after allocation of vahine interpol");
- for(itLearningNoiseEstim.GoToBegin(), itLearningNoiseProj.GoToBegin(); !itLearningNoiseEstim.IsAtEnd();
- ++itLearningNoiseEstim, ++itLearningNoiseProj) {
- RealType data = eigenStruct.interpolation->get(itLearningNoiseProj.Get() );
- itLearningNoiseEstim.Set(data);
- }
- // calculate nrmse
- CalculateNRMSE(eigenStruct, itLearningNoiseEstim);
-}
-
-/**
- * Calcul and sort the eigen vector of the gamma matrix
- * Sortering is in decrease order of eigen value
- */
-template<class TInputImage, class TOutputImage>
-void
-VahineGrsirFilter<TInputImage, TOutputImage>::CalculateEigenVector(VnlMatrixType &gamma, EigenStruct &eigenStruct){
- EigenListType eigenList;
-
- // calculation of eigenvector
- vahineDebug("VahineGrsirFilter : calculate eigenvector");
- // construct the eigenvector problem
- vnl_real_eigensystem eigen( m_RegularizationMatrix.GetVnlMatrix() * gamma );
- // extract eigenvector and ordering in decrease order
- //EigenListType eigenList;
- vnl_diag_matrix<vcl_complex<double> >::iterator it;
- for (it = eigen.D.begin(); it != eigen.D.end(); ++it) {
- VnlVectorType vector;
- vector.set_size( eigen.V.rows() );
- double value = (*it).real();
- for (int i=0; i < eigen.V.rows(); ++i) {
- vector(i) = eigen.V(i,it-eigen.D.begin() ).real();
- }
- EigenPairType currentEigenPair = std::make_pair(value, vector);
- eigenList.push_back(currentEigenPair);
- //vahineDebug("Eigen value = "<<currentEigenPair.first);
- }
-
- // sorting eigen Vector by decrease eigen value
- vahineDebug("VahineGrsirFilter : sortering eigen vector");
- eigenList.sort(SortPairDecrease<EigenPairType>() );
- eigenStruct.eigenVector = eigenList.front().second;
- vahineDebug("FirstEigenValue :"<<eigenList.front().first);
- vahineDebug("FirstEigenVector : "<<eigenStruct.eigenVector);
-}
-
-/**
- * Normelize Root Min Square Error calculation
- */
-template<class TInputImage, class TOutputImage>
-void VahineGrsirFilter<TInputImage, TOutputImage>::CalculateNRMSE(EigenStruct &eigenStruct,
- IteratorProjType &itLearningNoiseEstim) {
- vahineDebug("calculate nmrse");
- typename StatsFilterType::Pointer statsFilter = StatsFilterType::New();
- statsFilter->SetInput(m_Parameters);
- statsFilter->Update();
- RealType parameterMean = statsFilter->GetMean();
- RealType diffSum = 0;
- RealType diffParamSum = 0;
- ConstIteratorType itParameters( m_Parameters, m_Parameters->GetLargestPossibleRegion() );
- for(itParameters.GoToBegin(), itLearningNoiseEstim.GoToBegin();
- !itParameters.IsAtEnd();
- ++itParameters, ++itLearningNoiseEstim) {
- diffSum += pow(itParameters.Get() - itLearningNoiseEstim.Get(), 2);
- diffParamSum += pow(itParameters.Get() - parameterMean, 2);
- }
- eigenStruct.nrmse = vcl_sqrt(diffSum/diffParamSum);
-}
-
-/**
- * Calcul SIR criterion
- * ratio beetween the beetween-slices variance and intra-slice variance
- */
-template<class TInputImage, class TOutputImage>
-double VahineGrsirFilter<TInputImage, TOutputImage>::CalculateSIRC(VnlMatrixType &gamma, VnlVectorType &axe){
- RealType sGamma = dot_product(axe, gamma*axe);
- RealType sSigma = dot_product(axe, m_CovarianceMatrix*axe);
-
- return sGamma/sSigma;
-}
-
-} // end namspace otb
-
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/GrsirInverseFilter.h b/Code/Vahine/GrsirInverseFilter.h
deleted file mode 100644
index 4b8a970c66a0aa34d95d3b05f4785a07bea035b1..0000000000000000000000000000000000000000
--- a/Code/Vahine/GrsirInverseFilter.h
+++ /dev/null
@@ -1,162 +0,0 @@
-/*
- * GrsirInverseFilter.h
- * Date: $Date$
- * Version: $Revision$
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2010 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __GRSIRINVERSEFILTER_H_
-#define __GRSIRINVERSEFILTER_H_
-
-#include <bitset>
-#include <itkImageToImageFilter.h>
-#include <itkImageLinearIteratorWithIndex.h>
-#include <otb/BasicFilters/otbImageListToVectorImageFilter.h>
-#include <otb/Common/otbMultiToMonoChannelExtractROI.h>
-#include <otb/Common/otbImageList.h>
-#include <otb/IO/otbImage.h>
-
-#include "VahineMacro.h"
-#include "VahineInterpol.h"
-#include "VahineException.h"
-#include "InverseFilter.h"
-
-namespace otb {
-
-template <class TImage>
-class ITK_EXPORT VahineGrsirInverseFilter
- : public itk::ImageToImageFilter<TImage, TImage> {
-public:
-
- typedef VahineGrsirInverseFilter Self;
- typedef itk::ImageToImageFilter<TImage, TImage> Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
-
- // creation of SmartPointer
- itkNewMacro(Self);
- // runtime type information
- itkTypeMacro(VahineGrsirInverseFilter, itk::ImageToImageFilter);
-
- // Types to used for computation
- const static unsigned int Dimension = 2;
- typedef TImage ImageType;
- typedef typename ImageType::Pointer ImagePointer;
- typedef typename Superclass::InputImageConstPointer ImageConstPointer;
- typedef typename ImageType::InternalPixelType InternalPixelType;
- typedef typename itk::NumericTraits<InternalPixelType>::RealType RealType;
- typedef vnl_vector<InternalPixelType> VnlVectorType;
-
- // type for managing learning parameter
- typedef Image<InternalPixelType, Dimension> BandType;
- typedef typename BandType::Pointer BandPointer;
- typedef MultiToMonoChannelExtractROI<InternalPixelType,InternalPixelType> ExtractROImonoFilterType;
- typedef typename ExtractROImonoFilterType::Pointer ExtractROImonoFilterPointer;
- typedef VahineInverseFilter<ImageType, BandType> VahineInverseType;
- typedef typename VahineInverseType::Pointer VahineInversePointer;
-
- typedef itk::ImageLinearIteratorWithIndex<BandType> BandLinearIterator;
- // type for output
- typedef ImageList<BandType> InverseResultListType;
- typedef typename InverseResultListType::Pointer InverseResultListPointer;
- typedef ImageListToVectorImageFilter<InverseResultListType, ImageType> ImageListToVectorImageFilterType;
- typedef typename ImageListToVectorImageFilterType::Pointer ImageListToVectorImageFilterPointer;
-
- /**
- * 1 : 1 == inverse this parameter
- * 2-3 : regularization parameter
- * 00 : not use
- * 01 : optimal regularization parameter
- * 10 : custom regularization parameter
- * 11 : inverse all regularization parameter
- * 4-7 : inverse one value : index between 0 and 2^4 (16)
- * 8 : not use
- * exemple :
- * 0000 0011 inverse with optimal regularization value
- * 0001 0101 inverse with the value with index 1 (second value)
- * 0000 0111 inverse for all regularization values
- */
- typedef std::bitset<8> BandParameterType;
- const BandParameterType INVERSE_MASK;
- const BandParameterType OPTREGUL_MASK;
- const BandParameterType ALLREGUL_MASK;
- const BandParameterType CUSREGUL_MASK;
- const BandParameterType REGULVALUE_MASK;
-
- const unsigned int NbEigenVector;
-
- void PushBackParameterBand(BandParameterType);
-
- itkSetMacro( LearningParam, ImagePointer );
-
- bool GetDebug(){
- return true;
- };
-
- void PrintSelf(std::ostream& os, itk::Indent indent) const {
- Superclass::PrintSelf(os, indent);
- };
-
- void GenerateOutputInformation();
-
- void GenerateInputRequestedRegion();
-
-protected:
- VahineGrsirInverseFilter();
- void GenerateData();
-
-private:
-
- struct EigenStruct {
- VnlVectorType eigenVector;
- vahine::VahineInterpol<InternalPixelType> * interpolation;
- };
-
- // learning parameters : one band by parameter
- ImagePointer m_LearningParam;
- // result list
- InverseResultListPointer m_ResultList;
-
- typedef std::list<BandParameterType> ParameterListType;
- typedef ParameterListType::iterator ItParameterList;
- ParameterListType m_parameter;
-
- VahineGrsirInverseFilter(const Self&); //not implemented
- void operator=(const Self&); //not implemented
-
- inline bool checkParameter(BandParameterType b, BandParameterType mask){
- return ( (b &= mask) == mask);
- }
-
- unsigned int getNumberOfRegulParameters(BandPointer learnedParameters);
-
- BandPointer getLearningParameters(unsigned int bandNumber);
-
- unsigned int getIndexForBestRegul(BandPointer learnedParameters);
-
- unsigned int getIndexForCustRegul(BandParameterType param);
-
- unsigned int getNbRegularization(BandPointer learnedParameters);
-
- VnlVectorType getDataFromLine(BandPointer band, unsigned int, bool);
-
- void UpdateEigenStruct(BandPointer band, EigenStruct & eigenStruct, unsigned int index);
-
- void InverseImage(EigenStruct eigenStruct);
-
-};
-
-} // end namespace otb
-
-#include "GrsirInverseFilter.txx"
-
-#endif /* GRSIRINVERSEFILTER_H_ */
\ No newline at end of file
diff --git a/Code/Vahine/GrsirInverseFilter.txx b/Code/Vahine/GrsirInverseFilter.txx
deleted file mode 100644
index 356d54d8195c45e249cb42a65c8b2a15a146508e..0000000000000000000000000000000000000000
--- a/Code/Vahine/GrsirInverseFilter.txx
+++ /dev/null
@@ -1,269 +0,0 @@
-/*
- * GrsirInverseFilter.txx
- * Date: $Date$
- * Version: $Revision$
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2010 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINEGRSIRINVERSEFILTER_TXX
-#define __VAHINEGRSIRINVERSEFILTER_TXX
-
-#include "GrsirInverseFilter.h"
-
-namespace otb {
-
-template<class TImage>
-VahineGrsirInverseFilter<TImage>::VahineGrsirInverseFilter() :
- INVERSE_MASK(0x01),
- OPTREGUL_MASK(0x02),
- CUSREGUL_MASK(0x04),
- ALLREGUL_MASK(0x06),
- REGULVALUE_MASK(0x78),
- NbEigenVector(1)
-{
- m_ResultList = InverseResultListType::New();
-}
-
-template<class TImage>
-void VahineGrsirInverseFilter<TImage>::PushBackParameterBand(BandParameterType param){
- m_parameter.push_back(param);
-}
-
-template<class TImage>
-void VahineGrsirInverseFilter<TImage>::GenerateOutputInformation(){
- vahineDebug("VahineGrsirInverseFilter::GenerateOutputInformation");
- Superclass::GenerateOutputInformation();
-
- // get pointers to the input and output
- ImageConstPointer inputPtr = this->GetInput();
- ImagePointer outputPtr = this->GetOutput();
-
- if ( !inputPtr || !outputPtr ) {
- return;
- }
-
- unsigned int zSize = 0;
- unsigned int physicalNumber = 1;
- ItParameterList itParam = m_parameter.begin();
- while(itParam != m_parameter.end() ) {
- if(checkParameter(*itParam, INVERSE_MASK) ) {
- if(checkParameter(*itParam, ALLREGUL_MASK) ) {
- zSize += getNumberOfRegulParameters(getLearningParameters(physicalNumber) );
- } else {
- // only one band for this physical parameter
- ++zSize;
- }
- }
- ++physicalNumber;
- ++itParam;
- }
- vahineDebug("Number of components per pixel for output = "<<zSize);
- //inputPtr->UpdateOutputInformation();
- outputPtr->SetLargestPossibleRegion( inputPtr->GetLargestPossibleRegion() );
- outputPtr->SetNumberOfComponentsPerPixel( zSize );
-}
-
-template<class TImage>
-void VahineGrsirInverseFilter<TImage>::GenerateInputRequestedRegion(){
- vahineDebug("VahineGrsirLearningFilter::GenerateInputRequestedRegion");
- ImagePointer inputPtr = const_cast<ImageType * >( this->GetInput() );
- typename ImageType::RegionType region = inputPtr->GetLargestPossibleRegion();
- vahineDebug("###################### region = "<<region);
- inputPtr->SetRequestedRegion(region);
-}
-
-/**
-* extract all learning parameters for one physical parameter
-*/
-template<class TImage>
-typename VahineGrsirInverseFilter<TImage>::BandPointer
-VahineGrsirInverseFilter<TImage>::getLearningParameters(unsigned int bandNumber){
- vahineDebug("getLearningParameters bandNumber="<<bandNumber);
- ExtractROImonoFilterPointer extractFilter = ExtractROImonoFilterType::New();
- extractFilter->SetChannel(bandNumber);
- extractFilter->SetInput(m_LearningParam);
- BandPointer learningParmeters = extractFilter->GetOutput();
- learningParmeters->Update();
- return learningParmeters;
-}
-
-/*
- * first value of learned parameters is the number of regul parameters
- */
-template<class TImage>
-unsigned int VahineGrsirInverseFilter<TImage>::getNumberOfRegulParameters(BandPointer learnedParameters){
- BandLinearIterator itLearnedParameter(learnedParameters, learnedParameters->GetLargestPossibleRegion() );
-
- itLearnedParameter.SetDirection(0); // x dimension
- itLearnedParameter.GoToBegin();
- return static_cast<unsigned int>(itLearnedParameter.Get() );
-}
-
-/*
- * second value of learned parameter is the best regul parameters
- */
-template<class TImage>
-unsigned int VahineGrsirInverseFilter<TImage>::getIndexForBestRegul(BandPointer learnedParameters){
- BandLinearIterator itLearnedParameter(learnedParameters, learnedParameters->GetLargestPossibleRegion() );
-
- itLearnedParameter.SetDirection(0); // x dimension
- itLearnedParameter.GoToBegin();
- ++itLearnedParameter; // skip first value
- return static_cast<unsigned int>(itLearnedParameter.Get() );
-}
-
-/** convert bitset parameter to unsigend int index
- */
-template<class TImage>
-unsigned int VahineGrsirInverseFilter<TImage>::getIndexForCustRegul(BandParameterType param){
- BandParameterType number = param &= REGULVALUE_MASK;
-
- number>>=3;
- unsigned int index = number.to_ulong();
- vahineDebug("getIndexForCustRegul "<< index);
- return index;
-}
-
-template<class TImage>
-unsigned int VahineGrsirInverseFilter<TImage>::getNbRegularization(BandPointer learnedParameters){
- BandLinearIterator itLearnedParameter(learnedParameters, learnedParameters->GetLargestPossibleRegion() );
-
- itLearnedParameter.SetDirection(0); // x dimension
- itLearnedParameter.GoToBegin();
- return static_cast<unsigned int>(itLearnedParameter.Get() );
-}
-
-/**
- * line start from 0
- */
-template<class TImage>
-typename VahineGrsirInverseFilter<TImage>::VnlVectorType
-VahineGrsirInverseFilter<TImage>::getDataFromLine(BandPointer band, unsigned int index, bool allLine){
- VnlVectorType vector;
- BandLinearIterator itLearnedParameter(band, band->GetLargestPossibleRegion() );
-
- itLearnedParameter.SetDirection(0);
- itLearnedParameter.GoToBegin();
- vahineDebug("getDataFromLine number = "<<index);
- for(unsigned int i=0; i<index; ++i) {
- itLearnedParameter.NextLine();
- }
- unsigned int size = band->GetLargestPossibleRegion().GetSize()[0];
- if(!allLine) {
- // first element is size
- size = static_cast<unsigned int>(itLearnedParameter.Get() );
- } else {
- // no size if all line requested
- --itLearnedParameter;
- }
- vahineDebug("getDataFromLine size = "<<size);
- vector.set_size(size);
- for(unsigned int j=0; j<vector.size(); ++j) {
- ++itLearnedParameter;
- vector[j] = static_cast<RealType>(itLearnedParameter.Get() );
- }
- return vector;
-}
-
-template<class TImage>
-void VahineGrsirInverseFilter<TImage>::UpdateEigenStruct(BandPointer band, EigenStruct & eigenStruct,
- unsigned int index){
- vahineDebug("UpdateEigenStruct index="<<index);
- // index start from 0
- unsigned int offset = index*(2+NbEigenVector);
- VnlVectorType xVector = getDataFromLine(band, 1+offset, false);
- VnlVectorType yVector = getDataFromLine(band, 2+offset, false);
- eigenStruct.eigenVector = getDataFromLine(band, 3+offset, true);
- eigenStruct.interpolation = new vahine::VahineInterpol<InternalPixelType>(xVector.size() );
- eigenStruct.interpolation->SetVx(xVector);
- eigenStruct.interpolation->SetVy(yVector);
- eigenStruct.interpolation->Update();
-}
-
-template<class TImage>
-void VahineGrsirInverseFilter<TImage>::GenerateData(){
- vahineDebug("VahineGrsirInverseFilter : GenerateData");
-
- // check if we have the same number of physical parameter
- m_LearningParam->Update();
- unsigned int nbBand = m_LearningParam->GetNumberOfComponentsPerPixel();
- unsigned int nbParam = m_parameter.size();
- if( nbBand != nbParam ) {
- vahineDebug(
- "Error parameter bitset list must be have the same number of parameter than the learning parameter cube");
- vahineDebug("nbBand="<<nbBand<<" nbParam="<<nbParam);
- throw GrsirInverseException();
- }
-
- ItParameterList itParam = m_parameter.begin();
- // band number start with index 1
- unsigned int bandNumber = 1;
- BandPointer currentLearningBand;
- while(itParam != m_parameter.end() ) {
- if(checkParameter(*itParam, INVERSE_MASK) ) {
- currentLearningBand = getLearningParameters(bandNumber);
- if( checkParameter(*itParam, ALLREGUL_MASK) ) {
- vahineDebug("InverseImage for all regularization parameter");
- unsigned int nbRegularization = getNbRegularization(currentLearningBand);
- vahineDebug("nbRegularization = "<<nbRegularization);
- for(unsigned int i=0; i<nbRegularization; ++i) {
- EigenStruct currentEigenStruct;
- UpdateEigenStruct(currentLearningBand, currentEigenStruct, i);
- InverseImage(currentEigenStruct);
- }
- } else {
- if(checkParameter(*itParam, CUSREGUL_MASK) ) {
- vahineDebug("InverseImage for custom regularization parameter");
- EigenStruct customEigenStruct;
- unsigned int index = getIndexForCustRegul(*itParam);
- UpdateEigenStruct(currentLearningBand, customEigenStruct, index);
- InverseImage(customEigenStruct);
- } else {
- // default choice
- vahineDebug("InverseImage for best regularization parameter");
- EigenStruct bestEigenStruct;
- // index = [1 ; n-1]
- unsigned int index = getIndexForBestRegul(currentLearningBand);
- UpdateEigenStruct(currentLearningBand, bestEigenStruct, index);
- InverseImage(bestEigenStruct);
- }
- }
- }
- ++bandNumber;
- ++itParam;
- } // end while
- // convert image band list to vector image
- ImageListToVectorImageFilterPointer imageBandList2VectorImage = ImageListToVectorImageFilterType::New();
- imageBandList2VectorImage->SetInput(m_ResultList);
- vahineDebug("image size = "<<m_ResultList->Size() );
- imageBandList2VectorImage->Update();
- vahineDebug("get ouput");
- this->GraftOutput(imageBandList2VectorImage->GetOutput() );
-}
-
-/**
- * inverse method for one hyperspectral image and eigenStructure
- */
-template<class TImage>
-void VahineGrsirInverseFilter<TImage>::InverseImage(EigenStruct eigenStruct){
- vahineDebug("InverseImage");
- VahineInversePointer m_Inverse = VahineInverseType::New();
- m_Inverse->SetInput(this->GetInput() );
- m_Inverse->SetGrsirAxe(eigenStruct.eigenVector);
- m_Inverse->SetInterpolation(eigenStruct.interpolation);
- m_ResultList->PushBack(m_Inverse->GetOutput() );
- m_Inverse->Update();
-}
-
-} // end namespace otb
-
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/GrsirLearningFilter.h b/Code/Vahine/GrsirLearningFilter.h
deleted file mode 100644
index 06348b350bd0f7dc52e2133b1870481585fffe2f..0000000000000000000000000000000000000000
--- a/Code/Vahine/GrsirLearningFilter.h
+++ /dev/null
@@ -1,311 +0,0 @@
-/*
- * GrsirLearningFilter.h
- * Date: $Date$
- * Version: $Revision$
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2010 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __GRSIRLEARNINGFILTER_H_
-#define __GRSIRLEARNINGFILTER_H_
-
-#include <list>
-#include <set>
-#include <utility>
-#include <math.h>
-
-#include <boost/lexical_cast.hpp>
-#include <gdal.h>
-
-#include <otb/IO/otbImage.h>
-#include <otb/BasicFilters/otbImageListToVectorImageFilter.h>
-#include <otb/BasicFilters/otbPerBandVectorImageFilter.h>
-#include <otb/BasicFilters/otbStreamingStatisticsImageFilter.h>
-#include <otb/Common/otbMultiToMonoChannelExtractROI.h>
-
-#include <itkArray.h>
-#include <itkImageToImageFilter.h>
-#include <itkImageRegionConstIterator.h>
-#include <itkImageLinearIteratorWithIndex.h>
-#include <itkMatrix.h>
-
-#include <vnl/vnl_vector_ref.h>
-#include <vnl/algo/vnl_real_eigensystem.h>
-
-#include "VahineMacro.h"
-#include "VahineInterpol.h"
-#include "VahineSlice.h"
-#include "VahineXml.h"
-
-#include "CenteringImageFilter.h"
-#include "TikonovFilter.h"
-
-namespace otb {
-
-class GrsirLearningException { };
-
-// functor predicat for order our custom pair
-// don't use the internal pair operator< because it use the two elements for
-// comparaison
-// and we don't have operator< in the second type
-template < typename PairType >
-class SortPairIncrease {
-public:
- bool operator()( const PairType & left , const PairType & right ){
- return left.first < right.first;
- }
-
-};
-
-template < typename PairType >
-class SortPairDecrease {
-public:
- bool operator()( const PairType & left , const PairType & right ){
- return left.first > right.first;
- }
-
-};
-
-class ManagedTypes {
-public:
- // Convert a VariableLengthVector to an Array.
- template < class TElement >
- static itk::Array<TElement> ToManagedVariableLengthVector( itk::VariableLengthVector<TElement> vector ){
- typedef typename itk::Array<TElement> ArrayType;
- ArrayType result;
- result.SetSize(vector.Size() );
- for (unsigned int i=0; i<vector.Size(); ++i)
- result[i] = static_cast<TElement>(vector[i]);
- return result;
- }
-
-};
-
-template <class TImage>
-class ITK_EXPORT VahineGrsirLearningFilter
- : public itk::ImageToImageFilter<TImage, TImage> {
-public:
-
- typedef VahineGrsirLearningFilter Self;
- typedef itk::ImageToImageFilter<TImage, TImage> Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
-
- // creation of SmartPointer
- itkNewMacro(Self);
- // runtime type information
- itkTypeMacro(VahineGrsirLearningFilter, itk::ImageToImageFilter);
-
- // Type to use for computation
- const static unsigned int Dimension = 2;
- typedef TImage ImageType;
- typedef typename ImageType::Pointer ImagePointer;
- typedef typename Superclass::InputImageConstPointer ImageConstPointer;
- typedef typename ImageType::PixelType PixelType;
- typedef typename ImageType::InternalPixelType InternalPixelType;
- typedef typename ImageType::IndexType IndexType;
- typedef typename ImageType::RegionType RegionType;
-
- typedef typename itk::NumericTraits<InternalPixelType>::RealType RealType;
- typedef typename itk::VariableSizeMatrix<RealType> MatrixType;
-
- typedef Image<InternalPixelType, Dimension> BandType;
- typedef typename BandType::PixelType BandPixelType;
- typedef typename BandType::IndexType BandIndexType;
- typedef typename BandType::Pointer BandPointer;
- typedef itk::ImageRegionIterator<BandType> BandIterator;
- typedef itk::ImageLinearIteratorWithIndex<BandType> BandLinearIterator;
- typedef ImageList<BandType> ImageBandListType;
- typedef typename ImageBandListType::Pointer ImageBandListPointer;
-
- typedef itk::ImageRegionConstIterator<ImageType> ConstImageIteratorType;
-
- typedef vnl_matrix<RealType> VnlMatrixType;
- typedef vnl_vector<RealType> VnlVectorType;
-
- typedef std::pair<double, VnlVectorType> EigenPairType;
- typedef std::list<EigenPairType> EigenListType;
-
- typedef std::vector<double> RegulVectorType;
- typedef RegulVectorType* RegulVectorPointer;
-
- typedef Image<InternalPixelType, 2> ProjType;
- typedef typename ProjType::Pointer ProjPointer;
- typedef itk::ImageRegionIterator<ProjType> IteratorProjType;
-
- // used internal filter
- typedef otb::StreamingStatisticsImageFilter<BandType> StatsFilterType;
- typedef otb::VahineTikonovFilter<ImageType> TikonovFilterType;
-
- typedef typename std::list<vahine::Slice<InternalPixelType> > SliceListType;
- typedef typename SliceListType::iterator IteratorSliceListType;
-
- typedef XmlLog* XmlLogPointer;
-
- /*** members accessor ***/
- /** tikonov regulation parameters */
- itkGetMacro( RegulParameters, RegulVectorPointer);
- void SetRegulParameters(RegulVectorPointer);
-
- itkGetMacro( NSlice, int);
- itkSetMacro( NSlice, int);
-
- /** learning spectrums database image and
- * learning noisy spectrums database image
- */
- itkSetMacro( LearningNoiseData, ImagePointer );
- itkSetMacro( LearningData, ImagePointer );
-
- /** mesurement of difference for two closest parameter value
- * if differrence is smaller than epsilon the two values are
- * in the same slice.
- */
- itkGetMacro( Epsilon, float);
- itkSetMacro( Epsilon, float);
-
- /** if true calculate and save only the best projection
- * use the minimal nrmse criteria */
- itkSetMacro( SaveOnlyBestRegulParam, bool );
-
- /** Number of eigen vector that must be calculated */
- itkGetMacro( NumberOfEigenVector, unsigned int );
- itkSetMacro( NumberOfEigenVector, unsigned int );
-
- void SetTikonovFilter();
-
- void SetXmlLog(XmlLogPointer xmlLog);
-
- void PushBackLearningBand(unsigned int);
-
- bool GetDebug(){
- return true;
- };
-
- void PrintSelf(std::ostream& os, itk::Indent indent) const {
- Superclass::PrintSelf(os, indent);
- };
-
- SliceListType &GetSliceList();
-
- void GenerateOutputInformation();
-
- void GenerateInputRequestedRegion();
-
-protected:
- VahineGrsirLearningFilter();
- typedef itk::CenteringImageFilter<BandType, BandType> CenteringFilterType;
- typedef PerBandVectorImageFilter<ImageType, ImageType, CenteringFilterType> PerBandVectorImageFilterType;
- void GenerateData();
-
-private:
-
- /*** INTERNAL TYPE ***/
-
- // ParameterPairType :
- // * first element is Y value (float)
- // * second element is its index
- typedef std::pair<RealType, BandIndexType> ParameterPairType;
- typedef std::list<ParameterPairType> ParameterListType;
- typedef typename ParameterListType::const_iterator ParameterListIteratorType;
- typedef std::vector<unsigned int> BandNumberListType;
- typedef BandNumberListType::iterator BandNumberListIterator;
- typedef itk::ImageRegionConstIterator<BandType> ConstBandIteratorType;
- typedef MultiToMonoChannelExtractROI<BandPixelType, BandPixelType> ExtractChannelType;
- typedef ImageListToVectorImageFilter<ImageBandListType, ImageType> ImageListToVectorImageFilterType;
- typedef typename ImageListToVectorImageFilterType::Pointer ImageListToVectorImageFilterPointer;
-
- itkSetMacro( Parameters, BandPointer );
-
- struct EigenStruct {
- double nrmse;
- VnlVectorType eigenVector;
- VahineInterpol<RealType> * interpolation;
- };
-
- /*** END INTERNAL TYPE ***/
-
- /** tikonov regulation parameters */
- RegulVectorPointer m_RegulParameters;
-
- MatrixType m_RegularizationMatrix;
- MatrixType m_CovarianceMatrix;
- int m_NSlice;
- BandPointer m_Parameters;
-
- /** learning spectrums database image and
- * learning noisy spectrums database image
- */
- ImagePointer m_LearningData;
- ImagePointer m_LearningNoiseData;
-
- /** mesurement of difference for two closest parameter value
- * if differrence is smaller than epsilon the two values are
- * in the same slice.
- */
- float m_Epsilon;
-
- /** if true calculate and save only the best projection
-* use the minimal nrmse criteria */
- bool m_SaveOnlyBestRegulParam;
-
- /** Number of eigen vector that must be calculated */
- unsigned int m_NumberOfEigenVector;
-
- SliceListType m_SliceList;
- /** list of number of parameter band for grsisr calculation */
- BandNumberListType m_LearningParamBand;
-
- typename CenteringFilterType::Pointer m_CenteringFilter;
- typename PerBandVectorImageFilterType::Pointer m_CenteringVectorFilter;
- typename TikonovFilterType::Pointer m_TikonovFilter;
- BandPointer m_currentImageInv;
- ImageBandListPointer m_ImageBandList;
- ImagePointer m_LearningCenteringImage;
-
- ImagePointer m_OutputImage;
-
- RegionType m_OutputLargestPossibleRegion;
-
- XmlLogPointer m_XmlLog;
-
- VahineGrsirLearningFilter(const Self&); //not implemented
- void operator=(const Self&); //not implemented
-
- void Slicing(ParameterListType parameterList, ImagePointer image);
-
- void MakeMean(unsigned int nObserv);
-
- void UpdateTikonovFilter(double regulParam);
-
- void EstimateEDR(VnlMatrixType &gamma, unsigned int numberOfBand, EigenStruct &eigenStruct);
-
- double CalculateSIRC(VnlMatrixType &gamma, VnlVectorType &axe);
-
- void CalculateNRMSE(EigenStruct &eigenStruct, IteratorProjType &itLearningNoiseEstim);
-
- void CalculateEigenVector(VnlMatrixType &gamma, EigenStruct &eigenStruct);
-
- void Learning(unsigned int numberOfBand, unsigned int nObserv);
-
- void AddRegularizationParameters(BandLinearIterator &itLearnedFunctionnal, unsigned int bestRegul,
- RegulVectorPointer regulVector);
-
- void AddLearnedFunctionnal(EigenStruct eigenStruct, BandLinearIterator &itLearnedFunctionnal);
-
- void CompleteLine(BandLinearIterator &itLearnedFunctionnal);
-
-};
-
-} // End namespace otb
-
-#include "GrsirLearningFilter.txx"
-
-#endif /* GRSIRLEARNINGFILTER_H_ */
\ No newline at end of file
diff --git a/Code/Vahine/GrsirLearningFilter.txx b/Code/Vahine/GrsirLearningFilter.txx
deleted file mode 100644
index ba10d689caf480250f033eda590a8557aed60d4c..0000000000000000000000000000000000000000
--- a/Code/Vahine/GrsirLearningFilter.txx
+++ /dev/null
@@ -1,646 +0,0 @@
-/*
- * GrsirLearningFilter.txx
- * Date: $Date$
- * Version: $Revision$
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2010 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINEGRSIRLEARNINGFILTER_TXX
-#define __VAHINEGRSIRLEARNINGFILTER_TXX
-
-#include "GrsirLearningFilter.h"
-
-namespace otb {
-
-template<class TImage>
-VahineGrsirLearningFilter<TImage>::VahineGrsirLearningFilter(){
- m_Epsilon = 1E-6;
- m_SaveOnlyBestRegulParam = false;
- m_NumberOfEigenVector = 1; // several eigen vector not completly implemented
- // !
- /** default regul parameter is range [1E-15, 1E-2]
- */
- m_RegulParameters = new RegulVectorType();
- for(int i=-15; i<-2; ++i) {
- m_RegulParameters->push_back(pow(10.0,i) );
- }
- m_CenteringFilter = CenteringFilterType::New();
- m_CenteringVectorFilter = PerBandVectorImageFilterType::New();
- m_CenteringVectorFilter->SetFilter(m_CenteringFilter);
- m_ImageBandList = ImageBandListType::New();
- m_XmlLog = new XmlLog("VahineGrsirLearningFilter");
-}
-
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::SetRegulParameters(RegulVectorPointer p){
- delete m_RegulParameters;
- m_RegulParameters = p;
-}
-
-/**
- * X size is the number of component of one eigen vector : it's the number of band of
- * a learning database image
- * Y size is the sum of :
- * one : regularization parameters list
- * one : Xvalues of functionnal
- * one : Yvalues of functionnal
- * the number of eigen vector requested
- * multiply by the number of regul parameters requested, one if only the best is requested
- */
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::GenerateOutputInformation(){
- vahineDebug("VahineGrsirLearningFilter::GenerateOutputInformation");
- // call the superclass' implementation of this method
- Superclass::GenerateOutputInformation();
-
- // get pointers to the input and output
- ImageConstPointer inputPtr = this->GetInput();
- ImagePointer outputPtr = this->GetOutput();
-
- if ( !inputPtr || !outputPtr ) {
- return;
- }
-
- typename TImage::SizeType outputSize;
- m_LearningData->UpdateOutputInformation();
- outputSize[0] = m_LearningData->GetNumberOfComponentsPerPixel();
- outputSize[1] = 1+(2+m_NumberOfEigenVector)*m_RegulParameters->size();
- vahineDebug("outputSize "<<outputSize<<","<<m_NumberOfEigenVector<<","<<m_RegulParameters->size() );
- if( m_SaveOnlyBestRegulParam ) {
- outputSize[1] = 3+m_NumberOfEigenVector;
- }
-
- m_OutputLargestPossibleRegion.SetSize( outputSize );
- outputPtr->SetLargestPossibleRegion( m_OutputLargestPossibleRegion );
- unsigned int nParam = inputPtr->GetNumberOfComponentsPerPixel();
-
- // paramBand not set then all param are calculated [1, 2, ...]
- if(m_LearningParamBand.empty() ) {
- vahineDebug("learningParamBand is not set ! Using all bands.");
- for(unsigned int i=0; i<nParam; ++i) {
- m_LearningParamBand.push_back(i+1);
- }
- }
- unsigned int numberOfBands = m_LearningParamBand.size();
- vahineDebug("Output number of bands = "<<numberOfBands);
- outputPtr->SetNumberOfComponentsPerPixel( numberOfBands );
-}
-
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::GenerateInputRequestedRegion(){
- vahineDebug("VahineGrsirLearningFilter::GenerateInputRequestedRegion");
- ImagePointer inputPtr = const_cast<ImageType * >( this->GetInput() );
- typename ImageType::IndexType index = inputPtr->GetLargestPossibleRegion().GetIndex();
- typename ImageType::SizeType size;
- size.Fill(0);
- typename ImageType::RegionType region;
- region.SetSize(size);
- region.SetIndex(index);
- inputPtr->SetRequestedRegion(region);
-}
-
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::Slicing(ParameterListType parameterList, ImagePointer image){
- unsigned int nObserv = image->GetLargestPossibleRegion().GetNumberOfPixels();
-
- // distribute value in several class
-
- vahineDebug("VahineGrsirFilter : distribute data into class param="<<parameterList.size()<<" observ="<<nObserv);
- ParameterListIteratorType itParameterList = parameterList.begin();
- ParameterListIteratorType itLastParameter = parameterList.begin();
- if(!m_SliceList.empty() ) {
- vahineDebug("**** ERROR Slice list not empty ****");
- throw GrsirLearningException();
- }
- int i = 0;
- while(itParameterList != parameterList.end() ) {
- i++;
- BandIndexType index = (*itParameterList).second;
- PixelType pixel = image->GetPixel(index);
- if(m_SliceList.empty() ) {
- // push into list the first slice
- m_SliceList.push_back(vahine::Slice<InternalPixelType>(index, pixel, itParameterList->first) );
- } else {
- if( (itParameterList->first - itLastParameter->first)<m_Epsilon) {
- //vahineDebug("value closer : same slice ... current =
- // "<<itParameterList->first<<" last = "<<itLastParameter->first);
- // sum current pixel value to the current slice
- try{
- m_SliceList.back().add(index, pixel, itParameterList->first);
- }catch(...) {
- vahineDebug("add element failed");
- }
- } else {
- //vahineDebug("current pixel is in next class");
- // the current pixel is in next class
- // create the last mean object
- MakeMean(nObserv);
- // create a new slice with current pixel
- //vahineDebug("create a new slice with current pixel");
- m_SliceList.push_back(vahine::Slice<InternalPixelType>(index, pixel, itParameterList->first) );
- }
- itLastParameter++;
- }
- itParameterList++;
- }
- if(!parameterList.empty() ) {
- // create the last mean object
- MakeMean(nObserv);
- }
- if(m_SliceList.size() <= 1) {
- vahineDebug("*** ERROR Slice list must be have more than one element ***");
- throw GrsirLearningException();
- }
-}
-
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::MakeMean(unsigned int nObserv){
- m_SliceList.back().setMeanPixel();
- m_SliceList.back().setPropOfPixel(nObserv);
-}
-
-template<class TImage>
-typename VahineGrsirLearningFilter<TImage>::SliceListType &
-VahineGrsirLearningFilter<TImage>::GetSliceList(){
- return m_SliceList;
-}
-
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::SetTikonovFilter(){
- m_TikonovFilter = TikonovFilterType::New();
- m_TikonovFilter->SetInput(m_LearningData);
- //m_LearningData->Update();
- m_TikonovFilter->UpdateCov();
- m_CovarianceMatrix = m_TikonovFilter->GetCovMatrix();
- VnlMatrixType m = m_CovarianceMatrix.GetVnlMatrix();
- for(int i=0; i<1; i++) {
- for(int j=0; j<10; j++) {
- vahineDebug("cov = "<<m.get(i,j) );
- }
- }
-}
-
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::SetXmlLog(XmlLogPointer xmlLog){
- delete m_XmlLog;
- m_XmlLog = xmlLog;
-}
-
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::UpdateTikonovFilter(double regulParam){
- m_TikonovFilter->SetRegularizationParam(regulParam);
- m_TikonovFilter->Update();
- m_RegularizationMatrix = m_TikonovFilter->GetTikonovMatrix();
- VnlMatrixType m = m_RegularizationMatrix.GetVnlMatrix();
- for(int i=0; i<1; i++) {
- for(int j=0; j<10; j++) {
- vahineDebug("tiko = "<<m.get(i,j) );
- }
- }
-}
-
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::GenerateData(){
- vahineDebug("VahineGrsirLearningFilter : GenerateData");
- m_XmlLog->CreateElement("filter-grsir");
-
- SetTikonovFilter();
-
- // centering data
- m_CenteringVectorFilter->SetInput(m_LearningData);
- m_CenteringVectorFilter->Update();
- m_LearningCenteringImage = m_CenteringVectorFilter->GetOutput();
- unsigned int numberOfBand = m_LearningCenteringImage->GetNumberOfComponentsPerPixel();
- unsigned int nObserv = m_LearningCenteringImage->GetLargestPossibleRegion().GetNumberOfPixels();
-
- // Extract the correct parameter band
- // parameter is input image
- typename ExtractChannelType::Pointer extractParamBand = ExtractChannelType::New();
- extractParamBand->SetExtractionRegion(this->GetInput()->GetLargestPossibleRegion() );
- extractParamBand->SetInput(this->GetInput() );
-
- // loop on parameter
- BandNumberListIterator it = m_LearningParamBand.begin();
- while(it != m_LearningParamBand.end() ) {
- vahineDebug("Inversion for parameter : "<<(*it) );
- extractParamBand->SetChannel(*it);
- m_Parameters = extractParamBand->GetOutput();
- extractParamBand->Update();
-
- vahineDebug("Imagebandlist before size ="<<m_ImageBandList->Size() );
- //******************* Inverse call *********************/
- m_XmlLog->CreateElement("optimization");
- m_XmlLog->CreateAttribute("parameter",boost::lexical_cast<std::string>( (*it) ) );
- Learning(numberOfBand, nObserv);
- m_XmlLog->AppendToSup();
- //***************************************************/
- m_ImageBandList->Update();
-
- //TODO cleanup internal variable ?
- ++it;
- //break; // for debug one parameter
- }
-
- // convert image band list to vector image
- ImageListToVectorImageFilterPointer imageBandList2VectorImage = ImageListToVectorImageFilterType::New();
- m_ImageBandList->Update();
- imageBandList2VectorImage->SetInput(m_ImageBandList);
- vahineDebug("image size = "<<m_ImageBandList->Size() );
- imageBandList2VectorImage->UpdateOutputInformation();
- imageBandList2VectorImage->Update();
- vahineDebug("get ouput");
- m_OutputImage = imageBandList2VectorImage->GetOutput();
-
- // for debugging
- /*
- ConstVectorIteratorType itt(m_OutputImage, m_OutputImage->GetRequestedRegion());
- int h = 1;
- for(itt.GoToBegin();!itt.IsAtEnd(); ++itt){
- if(h>=120 && h<=128){
- vahineDebug("OutputImage ="<<itt.Get());
- }
- h++;
- }
- */
-
- this->GraftOutput(m_OutputImage);
-}
-
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::PushBackLearningBand(unsigned int band){
- vahineDebug("learned band number "<<band);
- m_LearningParamBand.push_back(band);
-}
-
-/**
- * Sortering Y parameters
- * Calculate the best eigen vector and apply it to input image
- * for inversion
- */
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::Learning(unsigned int numberOfBand, unsigned int nObserv){
-
- // Sortering Y !!
- vahineDebug("VahineGrsirFilter : Sortering Y");
- vahineDebug("Number of band learning = "<<numberOfBand);
- vahineDebug("Number of observation learning = "<<nObserv);
- m_Parameters->Update();
- vahineDebug("Number of parameter values = "<<m_Parameters->GetLargestPossibleRegion().GetNumberOfPixels() );
-
- ConstBandIteratorType itParameters( m_Parameters, m_Parameters->GetLargestPossibleRegion() );
- ParameterListType parameterList;
-
- for(itParameters.GoToBegin(); !itParameters.IsAtEnd(); ++itParameters) {
- BandIndexType index = itParameters.GetIndex();
- typename BandType::PixelType value = itParameters.Get();
- //vahineDebug("VahineGrsirFilter : index = "<<index<<" param v = "<< value);
- ParameterPairType element = make_pair(value, index);
- parameterList.push_back(element);
- }
-
- // sort pair by Y value
- vahineDebug("VahineGrsirFilter : Sort pair by Y value");
- parameterList.sort(SortPairIncrease<ParameterPairType>() );
-
- Slicing(parameterList, m_LearningCenteringImage);
-
- // calculate the between slice mean covariance matrix
- vahineDebug("VahineGrsirFilter : Calculate slice mean");
- vahineDebug("VahineGrsirFilter : number of slice = "<<m_SliceList.size() );
- vahineDebug("VahineGrsirFilter : Gamma size = "<<numberOfBand<<"*"<<numberOfBand);
- VnlMatrixType gamma = VnlMatrixType(numberOfBand, numberOfBand, 0);
- IteratorSliceListType itSlice = m_SliceList.begin();
- while(itSlice != m_SliceList.end() ) {
- VnlVectorType meanH = ManagedTypes::ToManagedVariableLengthVector<RealType>(itSlice->getMeanPixel() );
- //vahineDebug("VahineGrsirFilter : meanH = "<<meanH);
- //vahineDebug("propOfPixel = "<<itSlice->getPropOfPixel());
- gamma += itSlice->getPropOfPixel() * outer_product(meanH, meanH);
- ++itSlice;
- }
-
- // Grsir core : estimate GRSIR axes
- std::vector<double> nrmseVector;
- // eigenStructVEctor : used with no auto regularization option
- // all regularization and projection are used
- std::vector<EigenStruct> eigenStructVector;
- double regulMinNrmse;
- int regulIndex;
- EigenStruct bestEigenStruct;
- RegulVectorType::iterator itRegulParam = m_RegulParameters->begin();
- while( itRegulParam != m_RegulParameters->end() ) {
- // set regularization matrix
- double regul = *itRegulParam;
- UpdateTikonovFilter(regul);
- EigenStruct currentEigenStruct;
- // ********************** ESTIMATE EDR CALL**************************
- EstimateEDR(gamma, numberOfBand, currentEigenStruct);
- vahineDebug("Regul = "<<regul<<" Nrmse = "<<currentEigenStruct.nrmse);
- nrmseVector.push_back(currentEigenStruct.nrmse);
- if(!m_SaveOnlyBestRegulParam) {
- eigenStructVector.push_back(currentEigenStruct);
- }
- if(itRegulParam == m_RegulParameters->begin() ) {
- // first
- bestEigenStruct = currentEigenStruct;
- regulMinNrmse = regul;
- regulIndex = itRegulParam-m_RegulParameters->begin();
- } else {
- if(bestEigenStruct.nrmse>currentEigenStruct.nrmse) {
- bestEigenStruct = currentEigenStruct;
- regulMinNrmse = regul;
- regulIndex = itRegulParam-m_RegulParameters->begin();
- }
- }
- ++itRegulParam;
- }
- std::vector<double>::iterator it = nrmseVector.begin();
- while(it != nrmseVector.end() ) {
- vahineDebug("all nrmse = "<<*it);
- ++it;
- }
-
- vahineDebug("min nrmse is "<<bestEigenStruct.nrmse<<" regul value is "<<regulMinNrmse);
- //vahineDebug("gamma "<<gamma);
- double cgrsir = CalculateSIRC(gamma, bestEigenStruct.eigenVector);
- vahineDebug("cgrsir = "<<cgrsir);
- m_XmlLog->CreateElementTextApp("nrmse", boost::lexical_cast<std::string>(bestEigenStruct.nrmse) );
- m_XmlLog->CreateElement("regularisation");
- m_XmlLog->CreateAttribute("index",boost::lexical_cast<std::string>(regulIndex) );
- m_XmlLog->CreateTextNode(boost::lexical_cast<std::string>(regulMinNrmse) );
- m_XmlLog->AppendToSup();
- m_XmlLog->CreateElementTextApp("cgrsir",boost::lexical_cast<std::string>(cgrsir) );
-
- /*** Output fonctionnal interpolation data ***/
-
- BandPointer pLearnedFunctionnal = BandType::New();
- pLearnedFunctionnal->SetRegions(m_OutputLargestPossibleRegion);
- pLearnedFunctionnal->Allocate();
- BandLinearIterator itLearnedFunctionnal(pLearnedFunctionnal, pLearnedFunctionnal->GetRequestedRegion() );
- itLearnedFunctionnal.SetDirection(0); // x dimension
-
- AddRegularizationParameters(itLearnedFunctionnal, regulIndex, m_RegulParameters);
-
- if(m_SaveOnlyBestRegulParam) {
- vahineDebug("Functionnal for best regularization parameter");
- AddLearnedFunctionnal(bestEigenStruct, itLearnedFunctionnal);
- } else {
- vahineDebug("Functionnal for all regularization parameter");
- typename std::vector<EigenStruct>::iterator it = eigenStructVector.begin();
- while(it != eigenStructVector.end() ) {
- AddLearnedFunctionnal(*it, itLearnedFunctionnal);
- ++it;
- }
- }
-
- m_ImageBandList->PushBack(pLearnedFunctionnal);
-
- // clear list of slice for next parameter @see Slicing
- m_SliceList.clear();
-}
-
-/**
- * this method set several regularization param
- * line is list of regularization params
- * first value is optimal value
- */
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::AddRegularizationParameters(BandLinearIterator &itLearnedFunctionnal,
- unsigned int bestRegul,
- RegulVectorPointer regulVector){
- itLearnedFunctionnal.GoToBegin();
- itLearnedFunctionnal.Set(m_RegulParameters->size() );
- ++itLearnedFunctionnal;
- itLearnedFunctionnal.Set(bestRegul);
- ++itLearnedFunctionnal;
- RegulVectorType::iterator itRegulParam = m_RegulParameters->begin();
- // regularization not use in inverse filter but we save the logarithm value
- while(itRegulParam != m_RegulParameters->end() ) {
- itLearnedFunctionnal.Set(static_cast<RealType>(log10(*itRegulParam) ) );
- ++itRegulParam;
- ++itLearnedFunctionnal;
- }
- CompleteLine(itLearnedFunctionnal);
-
-}
-
-/**
- * Add eigen vector and functionnal interpolation values to a band pointer for ouput
- */
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::AddLearnedFunctionnal(EigenStruct eigenStruct,
- BandLinearIterator &itLearnedFunctionnal){
- itLearnedFunctionnal.NextLine();
- itLearnedFunctionnal.GoToBeginOfLine();
- VnlVectorType xvalues = eigenStruct.interpolation->GetVx();
- VnlVectorType yvalues = eigenStruct.interpolation->GetVy();
- // next line is xvalue
- vahineDebug("size="<<xvalues.size() );
- itLearnedFunctionnal.Set(static_cast<RealType>(xvalues.size() ) );
- ++itLearnedFunctionnal;
- for(int i=0; i<xvalues.size(); ++i) {
- RealType xval = static_cast<RealType>(xvalues[i]);
- itLearnedFunctionnal.Set(xval);
- vahineDebug("index="<<itLearnedFunctionnal.GetIndex()<<" xval="<<xval);
- ++itLearnedFunctionnal;
- }
- CompleteLine(itLearnedFunctionnal);
- itLearnedFunctionnal.NextLine();
- itLearnedFunctionnal.GoToBeginOfLine();
- // next line is yvalue
- itLearnedFunctionnal.Set(static_cast<RealType>(yvalues.size() ) );
- ++itLearnedFunctionnal;
- for(int i=0; i<yvalues.size(); ++i) {
- RealType yval = static_cast<RealType>(yvalues[i]);
- itLearnedFunctionnal.Set(yval);
- vahineDebug("index="<<itLearnedFunctionnal.GetIndex()<<" yval="<<yval);
- ++itLearnedFunctionnal;
- }
- CompleteLine(itLearnedFunctionnal);
- itLearnedFunctionnal.NextLine();
- itLearnedFunctionnal.GoToBeginOfLine();
- // next line is best eigen vector (first eigen vector)
- for(int i=0; i<eigenStruct.eigenVector.size(); ++i) {
- itLearnedFunctionnal.Set(eigenStruct.eigenVector[i]);
- ++itLearnedFunctionnal;
- }
-}
-
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::EstimateEDR(VnlMatrixType &gamma, unsigned int numberOfBand,
- EigenStruct &eigenStruct){
- // calculate eigen vector
- CalculateEigenVector(gamma, eigenStruct);
-
- vahineDebug("VahineGrsirFilter : projection of learning data");
-
- //TODO refactor this two loop
- // project learning value to first eigen vector
- ConstImageIteratorType itLearningData(m_LearningData, m_LearningData->GetLargestPossibleRegion() );
- m_LearningData->Update();
- ProjPointer learningDataProj = ProjType::New();
- learningDataProj->SetRegions( m_LearningData->GetLargestPossibleRegion() );
- learningDataProj->Allocate();
- IteratorProjType itLearningDataProj(learningDataProj, learningDataProj->GetLargestPossibleRegion() );
-
- for(itLearningData.GoToBegin(), itLearningDataProj.GoToBegin(); !itLearningData.IsAtEnd();
- ++itLearningData, ++itLearningDataProj) {
- PixelType currentPixel = itLearningData.Get();
- VnlVectorType pixel = ManagedTypes::ToManagedVariableLengthVector<RealType>(currentPixel);
- RealType tmp = static_cast<RealType>(dot_product(pixel, eigenStruct.eigenVector) );
- itLearningDataProj.Set(tmp);
- }
-
- // project learning noise value to first eigen vector
- ConstImageIteratorType itLearningNoise(m_LearningNoiseData, m_LearningNoiseData->GetLargestPossibleRegion() );
- m_LearningNoiseData->Update();
- ProjPointer learningNoiseProj = ProjType::New();
- learningNoiseProj->SetRegions( m_LearningNoiseData->GetLargestPossibleRegion() );
- learningNoiseProj->Allocate();
- IteratorProjType itLearningNoiseProj(learningNoiseProj, learningNoiseProj->GetLargestPossibleRegion() );
-
- for(itLearningNoise.GoToBegin(), itLearningNoiseProj.GoToBegin(); !itLearningNoise.IsAtEnd();
- ++itLearningNoise, ++itLearningNoiseProj) {
- PixelType currentPixel = itLearningNoise.Get();
- VnlVectorType pixel = ManagedTypes::ToManagedVariableLengthVector<RealType>(currentPixel);
- RealType tmp = static_cast<RealType>(dot_product(pixel, eigenStruct.eigenVector) );
- itLearningNoiseProj.Set(tmp);
- }
- //end TODO refactor
-
- // calculate the mean values of the current functionnal
- IteratorSliceListType itSlice = m_SliceList.begin();
- unsigned int sizeOfInterpol = m_SliceList.size();
- while( itSlice != m_SliceList.end() ) {
- itSlice->calculateMean(learningDataProj);
- ++itSlice;
- }
- m_SliceList.sort();
-
- vahineDebug("get mean and param sizeOfInterpol = "<<sizeOfInterpol);
- VnlVectorType dataVector = VnlVectorType(sizeOfInterpol);
- VnlVectorType paramVector = VnlVectorType(sizeOfInterpol);
- unsigned int i = 0;
- itSlice = m_SliceList.begin();
- while( itSlice != m_SliceList.end() ) {
- dataVector[i] = (*itSlice).getMeanData();
- paramVector[i] = (*itSlice).getMeanParam();
- vahineDebug("data "<<dataVector[i]<<" param "<<paramVector[i]);
- ++itSlice;
- i++;
- }
-
- // interpolation for learning noise
- vahineDebug("interpolation for learning noise data "<<m_LearningNoiseData->GetLargestPossibleRegion() );
- ProjPointer learningNoiseEstim = ProjType::New();
- learningNoiseEstim->SetRegions( m_LearningNoiseData->GetLargestPossibleRegion() );
- learningNoiseEstim->Allocate();
- IteratorProjType itLearningNoiseEstim(learningNoiseEstim, learningNoiseEstim->GetLargestPossibleRegion() );
-
- eigenStruct.interpolation = new VahineInterpol<RealType>(dataVector.size() );
- eigenStruct.interpolation->SetVx(dataVector);
- eigenStruct.interpolation->SetVy(paramVector);
- eigenStruct.interpolation->Update();
- vahineDebug("after allocation of vahine interpol");
- for(itLearningNoiseEstim.GoToBegin(), itLearningNoiseProj.GoToBegin(); !itLearningNoiseProj.IsAtEnd();
- ++itLearningNoiseEstim, ++itLearningNoiseProj) {
- RealType data = eigenStruct.interpolation->get(itLearningNoiseProj.Get() );
- itLearningNoiseEstim.Set(data);
- }
- // calculate nrmse
- CalculateNRMSE(eigenStruct, itLearningNoiseEstim);
-}
-
-/**
- * Calcul and sort the eigen vector of the gamma matrix
- * Sortering is in decrease order of eigen value
- */
-template<class TImage>
-void
-VahineGrsirLearningFilter<TImage>::CalculateEigenVector(VnlMatrixType &gamma, EigenStruct &eigenStruct){
- EigenListType eigenList;
-
- // calculation of eigenvector
- vahineDebug("VahineGrsirFilter : calculate eigenvector");
- // construct the eigenvector problem
- vnl_real_eigensystem eigen( m_RegularizationMatrix.GetVnlMatrix() * gamma );
- // extract eigenvector and ordering in decrease order
- //EigenListType eigenList;
- vnl_diag_matrix<vcl_complex<double> >::iterator it;
- for (it = eigen.D.begin(); it != eigen.D.end(); ++it) {
- VnlVectorType vector;
- vector.set_size( eigen.V.rows() );
- double value = (*it).real();
- for (int i=0; i < eigen.V.rows(); ++i) {
- vector(i) = eigen.V(i,it-eigen.D.begin() ).real();
- }
- EigenPairType currentEigenPair = std::make_pair(value, vector);
- eigenList.push_back(currentEigenPair);
- //vahineDebug("Eigen value = "<<currentEigenPair.first);
- }
-
- // sorting eigen Vector by decrease eigen value
- vahineDebug("VahineGrsirFilter : sortering eigen vector");
- eigenList.sort(SortPairDecrease<EigenPairType>() );
- eigenStruct.eigenVector = eigenList.front().second;
- vahineDebug("FirstEigenValue :"<<eigenList.front().first);
- vahineDebug("FirstEigenVector : "<<eigenStruct.eigenVector);
-}
-
-/**
- * Normelize Root Min Square Error calculation
- */
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::CalculateNRMSE(EigenStruct &eigenStruct,
- IteratorProjType &itLearningNoiseEstim) {
- vahineDebug("calculate nmrse");
- typename StatsFilterType::Pointer statsFilter = StatsFilterType::New();
- statsFilter->SetInput(m_Parameters);
- statsFilter->Update();
- RealType parameterMean = statsFilter->GetMean();
- RealType diffSum = 0;
- RealType diffParamSum = 0;
- ConstBandIteratorType itParameters( m_Parameters, m_Parameters->GetLargestPossibleRegion() );
- for(itParameters.GoToBegin(), itLearningNoiseEstim.GoToBegin();
- !itParameters.IsAtEnd();
- ++itParameters, ++itLearningNoiseEstim) {
- diffSum += pow(itParameters.Get() - itLearningNoiseEstim.Get(), 2);
- diffParamSum += pow(itParameters.Get() - parameterMean, 2);
- }
- eigenStruct.nrmse = vcl_sqrt(diffSum/diffParamSum);
-}
-
-/**
- * Calcul SIR criterion
- * ratio beetween the beetween-slices variance and intra-slice variance
- */
-template<class TImage>
-double VahineGrsirLearningFilter<TImage>::CalculateSIRC(VnlMatrixType &gamma, VnlVectorType &axe){
- RealType sGamma = dot_product(axe, gamma*axe);
- RealType sSigma = dot_product(axe, m_CovarianceMatrix*axe);
-
- return sGamma/sSigma;
-}
-
-/**
- * complete with null value
- */
-template<class TImage>
-void VahineGrsirLearningFilter<TImage>::CompleteLine(BandLinearIterator &itLearnedFunctionnal){
- while(!itLearnedFunctionnal.IsAtEndOfLine() ) {
- itLearnedFunctionnal.Set(0); //TODO replace with undefined value
- ++itLearnedFunctionnal;
- }
-}
-
-} // end namspace otb
-
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/InverseFilter.h b/Code/Vahine/InverseFilter.h
deleted file mode 100644
index d1f68fe227cca020aa0ee1a1097f1b1dd3b069c8..0000000000000000000000000000000000000000
--- a/Code/Vahine/InverseFilter.h
+++ /dev/null
@@ -1,75 +0,0 @@
-/*
- * InverseFilter.h
- *
- * Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINEINVERSEFILTER_H
-#define __VAHINEINVERSEFILTER_H
-
-#include <itkImageToImageFilter.h>
-#include <vnl/vnl_vector_ref.h>
-
-#include "VahineMacro.h"
-#include "VahineInterpol.h"
-
-namespace otb {
-template <class TInputImage, class TOutputImage> class ITK_EXPORT VahineInverseFilter :
- public itk::ImageToImageFilter<TInputImage, TOutputImage > {
-public:
- typedef VahineInverseFilter Self;
- typedef itk::ImageToImageFilter<TInputImage, TOutputImage> Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
- // creation of SmartPointer
- itkNewMacro(Self);
- // runtime type information
- itkTypeMacro(VahineInverseFilter, ImageToImageFilter);
-
- typedef TInputImage InputImageType;
- typedef typename TInputImage::Pointer InputImagePointer;
- typedef TOutputImage OutputImageType;
- typedef typename TOutputImage::Pointer OutputImagePointer;
- typedef typename TInputImage::PixelType PixelType;
- typedef typename TInputImage::InternalPixelType InternalPixelType;
- typedef typename itk::NumericTraits<InternalPixelType>::RealType RealType;
- typedef vnl_vector<InternalPixelType> VnlVectorType;
- typedef itk::ImageRegionConstIterator<InputImageType> ConstVectorIteratorType;
- typedef itk::ImageRegionIterator<OutputImageType> OutIterator;
-
- typedef typename vahine::VahineInterpol<InternalPixelType> * InterpolPointer;
-
- itkSetMacro( GrsirAxe, VnlVectorType );
- itkSetMacro( Interpolation, InterpolPointer)
- bool GetDebug(){
- return true;
- };
-protected:
- VahineInverseFilter();
- virtual ~VahineInverseFilter(){
- };
- virtual void PrintSelf(std::ostream& os, itk::Indent indent) const {
- Superclass::PrintSelf(os, indent);
- };
- void GenerateData();
-
-private:
- VahineInverseFilter( const Self& ); //purposely not implemented
- void operator=( const Self& ); //purposely not implemented
-
- VnlVectorType m_GrsirAxe;
- InterpolPointer m_Interpolation;
-};
-} //end namespace otb
-
-#include "InverseFilter.txx"
-
-#endif /* INVERSEFILTER_H_ */
\ No newline at end of file
diff --git a/Code/Vahine/InverseFilter.txx b/Code/Vahine/InverseFilter.txx
deleted file mode 100644
index e71491f98253b8c257687c32bdff0e0c3f12cb73..0000000000000000000000000000000000000000
--- a/Code/Vahine/InverseFilter.txx
+++ /dev/null
@@ -1,67 +0,0 @@
-/*
- * InverseFilter.txx
- *
- * Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-#ifndef __VAHINEINVERSEFILTER_TXX
-#define __VAHINEINVERSEFILTER_TXX
-#include "InverseFilter.h"
-
-namespace otb {
-template<class TInputImage, class TOutputImage>
-VahineInverseFilter<TInputImage,TOutputImage>::VahineInverseFilter(){
- m_GrsirAxe = VnlVectorType(0);
- m_Interpolation = 0;
-}
-
-template<class TInputImage, class TOutputImage>
-void VahineInverseFilter<TInputImage, TOutputImage>::GenerateData(){
- InputImagePointer pInputImage = const_cast<InputImageType *>(this->GetInput() );
-
- pInputImage->UpdateOutputInformation();
- OutputImagePointer pOutputImage = OutputImageType::New();
-
- unsigned int size = m_GrsirAxe.size();
- vahineDebug("size of grsir axe = "<<size);
- vahineDebug("input image largest possible region ="<<pInputImage->GetLargestPossibleRegion() );
- pOutputImage->SetRegions(pInputImage->GetLargestPossibleRegion() );
- pOutputImage->Allocate();
-
- ConstVectorIteratorType itConstInputImage(pInputImage, pInputImage->GetLargestPossibleRegion() );
- OutIterator itOutputImage(pOutputImage,pOutputImage->GetLargestPossibleRegion() );
-
- unsigned int nbInversedPixel = 0;
- for(itConstInputImage.GoToBegin(), itOutputImage.GoToBegin();
- !itConstInputImage.IsAtEnd();
- ++itConstInputImage, ++itOutputImage) {
- PixelType pixel = itConstInputImage.Get(); // input pixel
- RealType value = 0; // output value
- try{
- VnlVectorType p(pixel.GetSize() );
- p.copy_in(pixel.GetDataPointer() );
- value = static_cast<RealType>(dot_product(p,m_GrsirAxe) );
- if(value != 0) {
- value = m_Interpolation->get(value);
- ++nbInversedPixel;
- }
- }catch(...) {
- vahineDebug("Error in generate data");
- throw OVERFLOW; // for debugging
- }
- itOutputImage.Set(value);
- }
- vahineDebug("Number of pixel inversed = "<<nbInversedPixel<<" pointer output information"<<pOutputImage);
- this->GraftOutput(pOutputImage);
-}
-
-}
-
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/LICENCE b/Code/Vahine/LICENCE
deleted file mode 100644
index 8be4ec9a172905da8a570deb6427673ddfd01896..0000000000000000000000000000000000000000
--- a/Code/Vahine/LICENCE
+++ /dev/null
@@ -1,579 +0,0 @@
-CONTRAT DE LICENCE DE LOGICIEL LIBRE CeCILL
-
-Avertissement
-
-Ce contrat est une licence de logiciel libre issue d
-'une concertation
-entre ses auteurs afin que le respect de deux grands principes préside Ã
-sa rédaction:
-
- * d'une
-part, le respect des principes de diffusion des logiciels
-libres: accàs au code source, droits tendus confà ràs aux
- utilisateurs,
-* d 'autre part, la désignation d'un droit applicable, le droit
-franà ais, auquel elle est conforme, tant au regard du droit de la
-responsabilità civile que du droit de la proprià tà intellectuelle
-et de la protection qu
-'il offre aux auteurs et titulaires des
- droits patrimoniaux sur un logiciel.
-
-Les auteurs de la licence CeCILL (pour Ce[a] C[nrs] I[nria] L[ogiciel]
-L[ibre]) sont:
-
-Commissariat à l'Energie
-Atomique - CEA, tablissement public de
-recherche caractà re scientifique, technique et industriel, dont le
-sià ge est situà 25 rue Leblanc, immeuble Le Ponant D, 75015 Paris.
-
-Centre National de la Recherche Scientifique - CNRS, tablissement
-public caractà re scientifique et technologique, dont le sià ge est
- situà 3 rue Michel-Ange, 75794 Paris cedex 16.
-
-Institut National de Recherche en Informatique et en Automatique -
-INRIA, tablissement public caractà re scientifique et technologique,
-dont le sià ge est situà Domaine de Voluceau, Rocquencourt, BP 105, 78153
-Le Chesnay cedex.
-
-Prà ambule
-
-Ce contrat est une licence de logiciel libre dont l
-'objectif est de
-conférer aux utilisateurs la liberté de modification et de
-redistribution du logiciel régi par cette licence dans le cadre d'un
-modà le de diffusion en logiciel libre.
-
-L'exercice de ces libertés est assorti de certains devoirs à la charge
-des utilisateurs afin de préserver ce statut au cours des
-redistributions ultérieures.
-
-L'accessibilità au code source et les droits de copie, de modification
-et de redistribution qui en dà coulent ont pour contrepartie de n 'offrir
-aux utilisateurs qu'une garantie limità e et de ne faire peser sur
- l 'auteur du logiciel, le titulaire des droits patrimoniaux et les
-concédants successifs qu'une responsabilità restreinte.
-
-A cet gard l 'attention de l'utilisateur est attirà e sur les risques
-associàs au chargement,
- l 'utilisation, Ã la modification et/ou au
-développement et à la reproduction du logiciel par l'utilisateur tant
-donnà sa spà cificità de logiciel libre, qui peut le rendre complexe
-manipuler et qui le rà serve donc des dà veloppeurs ou des
-professionnels avertis possà dant des connaissances informatiques
- approfondies.Les utilisateurs sont donc invitàs charger et tester
- l 'adéquation du logiciel à leurs besoins dans des conditions permettant
-d'assurer la sà curità de leurs systà mes et/ou de leurs donnà es et, plus
-gà nà ralement, l 'utiliser et l'exploiter dans les mà mes conditions de
-sà curitÃ.Ce contrat peut tre reproduit et diffusà librement, sous
-rà serve de le conserver en l 'état, sans ajout ni suppression de clauses.
-
-Ce contrat est susceptible de s'appliquer tout logiciel dont le
-titulaire des droits patrimoniaux dà cide de soumettre l 'exploitation aux
-dispositions qu'il contient.
-
-Article 1 - DEFINITIONS
-
-Dans ce contrat, les termes suivants,
- lorsqu
-'ils seront écrits avec une
-lettre capitale, auront la signification suivante:
-
-Contrat: désigne le présent contrat de licence, ses éventuelles versions
-postérieures et annexes.
-
-Logiciel: désigne le logiciel sous sa forme de Code Objet et/ou de Code
-Source et le cas échéant sa documentation, dans leur état au moment de
-l'acceptation
-du Contrat par le LicenciÃ.
-
-Logiciel Initial : dà signe le Logiciel sous sa forme de Code Source et
-ventuellement de Code Objet et le cas chà ant sa documentation, dans
-leur tat au moment de leur premià re diffusion sous les termes du Contrat.
-
-Logiciel Modifià : dà signe le Logiciel modifià par au moins une
- Contribution.
-
-Code Source : dà signe l 'ensemble des instructions et des lignes de
-programme du Logiciel et auquel l'accàs est nà cessaire en vue de
-modifier le Logiciel.
-
-Code Objet : dà signe les fichiers binaires issus de la compilation du
-Code Source.
-
-Titulaire : dà signe le ou les dà tenteurs des droits patrimoniaux d
-'auteur
-sur le Logiciel Initial.
-
-Licencié: désigne le ou les utilisateurs du Logiciel ayant accepté le
-Contrat.
-
-Contributeur: désigne le Licencié auteur d'au
-moins une Contribution.
-
-Concà dant : dà signe le Titulaire ou toute personne physique ou morale
-distribuant le Logiciel sous le Contrat.
-
-Contribution : dà signe l
-'ensemble des modifications, corrections,
-traductions, adaptations et/ou nouvelles fonctionnalités intégrées dans
-le Logiciel par tout Contributeur, ainsi que tout Module Interne.
-
-Module: désigne un ensemble de fichiers sources y compris leur
-documentation qui permet de réaliser des fonctionnalités ou services
-supplémentaires à ceux fournis par le Logiciel.
-
-Module Externe: désigne tout Module, non dérivé du Logiciel, tel que ce
-Module et le Logiciel s'exÃ
-cutent dans des espaces d 'adressage
-différents, l'un appelant l
-'autre au moment de leur exécution.
-
-Module Interne: désigne tout Module lié au Logiciel de telle sorte
-qu'ils s
-'exécutent dans le même espace d'adressage.
-
-GNU GPL : dà signe la GNU General Public License dans sa version 2 ou
-toute version ultà rieure, telle que publià e par Free Software Foundation
- Inc.
-
-Parties : dà signe collectivement le Licencià et le Concà dant.
-
-Ces termes s
-'entendent au singulier comme au pluriel.
-
-
- Article 2 - OBJET
-
-Le Contrat a pour objet la concession par le Concédant au Licencié d'une
-licence non exclusive, cessible et mondiale du Logiciel telle que
-dà finie ci-
-apràs l
-'article 5 pour toute la durée de protection des droits
-portant sur ce Logiciel.
-
-
- Article 3 - ACCEPTATION
-
-3.1 L'acceptation
-par le Licencià des termes du Contrat est rà putà e
-acquise du fait du premier des faits suivants :
-
-* (i) le chargement du Logiciel par tout moyen notamment par
-tà là chargement partir d 'un serveur distant ou par chargement Ã
- partir d'un support physique;
-*(ii) le premier exercice par le Licencià de l
-'un quelconque des
- droits concédés par le Contrat.
-
-3.2 Un exemplaire du Contrat, contenant notamment un avertissement
-relatif aux spécificités du Logiciel, à la restriction de garantie et Ã
-la limitation à un usage par des utilisateurs expérimentés a été mis Ã
-disposition du Licencié préalablement à son acceptation telle que
-définie à l'article
-3.1 ci dessus et le Licencià reconnaà t en avoir pris
-connaissance.
-
-Article 4 - ENTREE EN VIGUEUR ET DUREE
-
-4.1 ENTREE EN VIGUEUR
-
-Le Contrat entre en vigueur la date de son acceptation par le LicenciÃ
-telle que dà finie en 3.1.
-
-4.2 DUREE
-
-Le Contrat produira ses effets pendant toute la durà e là gale de
-protection des droits patrimoniaux portant sur le Logiciel.
-
-Article 5 - ETENDUE DES DROITS CONCEDES
-
-Le Concà dant concà de au LicenciÃ, qui accepte, les droits suivants sur
-le Logiciel pour toutes destinations et pour la durà e du Contrat dans
-les conditions ci-apràs dà taillà es.
-
-Par ailleurs, si le Concà dant dà tient ou venait dà tenir un ou
-plusieurs brevets d
-'invention protégeant tout ou partie des
-fonctionnalités du Logiciel ou de ses composants, il s'engage ne pas
-opposer les ventuels droits confà ràs par ces brevets aux Licencià s
-successifs qui utiliseraient, exploiteraient ou modifieraient le
-Logiciel.En cas de cession de ces brevets,
-le Concà dant s
-'engage Ã
-faire reprendre les obligations du présent alinéa aux cessionnaires.
-
-
- 5.1 DROIT D'UTILISATION
-
-Le Licencià est autorisà utiliser le Logiciel, sans restriction quant
-aux domaines d
-'application, étant ci-après précisé que cela comporte:
-
- 1. la reproduction permanente ou provisoire du Logiciel en tout ou
- partie par tout moyen et sous toute forme.
-
- 2. le chargement, l'affichage ,
-l 'exécution, ou le stockage du
- Logiciel sur tout support.
-
- 3. la possibilité d'en observer, d 'en étudier, ou d'en tester le
-fonctionnement afin de dà terminer les idà es et principes qui sont
-la base de n
-'importe quel élément de ce Logiciel; et ceci,
- lorsque le Licencié effectue toute opération de chargement,
- d'affichage ,
-d 'exécution, de transmission ou de stockage du
- Logiciel qu'il est en droit d 'effectuer en vertu du Contrat.
-
-
- 5.2 DROIT D'APPORTER DES CONTRIBUTIONS
-
-Le droit d 'apporter des Contributions comporte le droit de traduire,
-d'adapter, d 'arranger ou d'apporter toute autre modification au Logiciel
-et le droit de reproduire le logiciel en rà sultant.
-
-Le Licencià est autorisà apporter toute Contribution au Logiciel sous
-rà serve de mentionner, de faà on explicite,
-son nom en tant qu
-'auteur de
-cette Contribution et la date de création de celle-ci.
-
-
- 5.3 DROIT DE DISTRIBUTION
-
-Le droit de distribution comporte notamment le droit de diffuser, de
-transmettre et de communiquer le Logiciel au public sur tout support et
-par tout moyen ainsi que le droit de mettre sur le marché à titre
-onéreux ou gratuit, un ou des exemplaires du Logiciel par tout procédé.
-
-Le Licencié est autorisé à distribuer des copies du Logiciel, modifié ou
-non, à des tiers dans les conditions ci-après détaillées.
-
-
- 5.3.1 DISTRIBUTION DU LOGICIEL SANS MODIFICATION
-
-Le Licencié est autorisé à distribuer des copies conformes du Logiciel,
-sous forme de Code Source ou de Code Objet, Ã condition que cette
-distribution respecte les dispositions du Contrat dans leur totalité et
-soit accompagnée:
-
- 1. d'un
-exemplaire du Contrat,
-
-2. d
-'un avertissement relatif à la restriction de garantie et de
- responsabilité du Concédant telle que prévue aux articles 8
- et 9,
-
-et que, dans le cas où seul le Code Objet du Logiciel est redistribué,
-le Licencié permette aux futurs Licenciés d'accÃ
-der facilement au Code
-Source complet du Logiciel en indiquant les modalitàs d 'accès, étant
-entendu que le coût additionnel d'acquisition du Code Source ne devra
-pas excà der le simple coà t de transfert des donnà es.
-
-5.3 .2 DISTRIBUTION DU LOGICIEL MODIFIE
-
-Lorsque le Licencià apporte une Contribution au Logiciel, les conditions
-de distribution du Logiciel Modifià en rà sultant sont alors soumises
-l
-'intégralité des dispositions du Contrat.
-
-Le Licencié est autorisé à distribuer le Logiciel Modifié, sous forme de
-code source ou de code objet, Ã condition que cette distribution
-respecte les dispositions du Contrat dans leur totalité et soit
-accompagnée:
-
- 1. d'un
-exemplaire du Contrat,
-
-2. d
-'un avertissement relatif à la restriction de garantie et de
- responsabilité du Concédant telle que prévue aux articles 8
- et 9,
-
-et que, dans le cas où seul le code objet du Logiciel Modifié est
-redistribué, le Licencié permette aux futurs Licenciés d'accÃ
-der
-facilement au code source complet du Logiciel Modifià en indiquant les
-modalitàs d 'accès, étant entendu que le coût additionnel d'acquisition
-du code source ne devra pas excà der le simple coà t de transfert des donnà es.
-
-5.3 .3 DISTRIBUTION DES MODULES EXTERNES
-
-Lorsque le Licencià a dà veloppà un Module Externe les conditions du
-Contrat ne s
-'appliquent pas à ce Module Externe, qui peut être distribué
-sous un contrat de licence différent.
-
-
- 5.3.4 COMPATIBILITE AVEC LA LICENCE GNU GPL
-
-Le Licencié peut inclure un code soumis aux dispositions d'une
-des
-versions de la licence GNU GPL dans le Logiciel modifià ou non et
-distribuer l
-'ensemble sous les conditions de la même version de la
-licence GNU GPL.
-
-Le Licencié peut inclure le Logiciel modifié ou non dans un code soumis
-aux dispositions d'une
-des versions de la licence GNU GPL et distribuer
-l
-'ensemble sous les conditions de la même version de la licence GNU GPL.
-
-
- Article 6 - PROPRIETE INTELLECTUELLE
-
-
- 6.1 SUR LE LOGICIEL INITIAL
-
-Le Titulaire est détenteur des droits patrimoniaux sur le Logiciel
-Initial. Toute utilisation du Logiciel Initial est soumise au respect
-des conditions dans lesquelles le Titulaire a choisi de diffuser son
-oeuvre et nul autre n'a
-la facultà de modifier les conditions de
-diffusion de ce Logiciel Initial.
-
-Le Titulaire s
-'engage à ce que le Logiciel Initial reste au moins régi
-par le Contrat et ce, pour la durée visée à l'article 4.2
-.
-
-6.2 SUR LES CONTRIBUTIONS
-
-Le Licencià qui a dà veloppà une Contribution est titulaire sur celle-ci
-des droits de proprià tà intellectuelle dans les conditions dà finies par
-la là gislation applicable.
-
-6.3 SUR LES MODULES EXTERNES
-
-Le Licencià qui a dà veloppà un Module Externe est titulaire sur celui-ci
-des droits de proprià tà intellectuelle dans les conditions dà finies par
-la là gislation applicable et reste libre du choix du contrat rà gissant
-sa diffusion.
-
-6.4 DISPOSITIONS COMMUNES
-
-Le Licencià s
-'engage expressément:
-
- 1. à ne pas supprimer ou modifier de quelque manière que ce soit les
- mentions de propriété intellectuelle apposées sur le Logiciel;
-
- 2. à reproduire à l'identique
-lesdites mentions de proprià tÃ
-intellectuelle sur les copies du Logiciel modifià ou non.
-
-Le Licencià s
-'engage à ne pas porter atteinte, directement ou
-indirectement, aux droits de propriété intellectuelle du Titulaire et/ou
-des Contributeurs sur le Logiciel et à prendre, le cas échéant, Ã
-l'Ã
-gard de son personnel toutes les mesures nà cessaires pour assurer le
-respect des dits droits de proprià tà intellectuelle du Titulaire et/ou
-des Contributeurs.
-
-Article 7 - SERVICES ASSOCIES
-
-7.1 Le Contrat n 'oblige en aucun cas le Concédant à la réalisation de
-prestations d'assistance technique ou de maintenance du Logiciel.
-
-Cependant le Concà dant reste libre de proposer ce type de services.Les
-termes et conditions d 'une telle assistance technique et/ou d'une telle
-maintenance seront alors dà terminà s dans un acte sà parÃ.Ces actes de
-maintenance et/
-ou assistance technique n
-'engageront que la seule
-responsabilité du Concédant qui les propose.
-
-7.2 De même, tout Concédant est libre de proposer, sous sa seule
-responsabilité, à ses licenciés une garantie, qui n'engagera
-que lui,
-lors de la redistribution du Logiciel et/ou du Logiciel Modifià et ce,
-dans les conditions qu 'il souhaite. Cette garantie et les modalités
-financières de son application feront l'objet d
-'un acte séparé entre le
-Concédant et le Licencié.
-
-
- Article 8 - RESPONSABILITE
-
-8.1 Sous réserve des dispositions de l'article
-8.2, le Licencià a la
-facultÃ, sous rà serve de prouver la faute du Concà dant concernÃ, de
-solliciter la rà paration du prà judice direct qu
-'il subirait du fait du
-Logiciel et dont il apportera la preuve.
-
-8.2 La responsabilité du Concédant est limitée aux engagements pris en
-application du Contrat et ne saurait être engagée en raison notamment:
-(i) des dommages dus à l'inexÃ
-cution, totale ou partielle, de ses
-obligations par le LicenciÃ, (ii) des dommages directs ou indirects
-dà coulant de l 'utilisation ou des performances du Logiciel subis par le
-Licencié et (iii) plus généralement d'un quelconque dommage indirect.En
-particulier, les Parties conviennent expressà ment que tout prà judice
-financier ou commercial(
- par exemple perte de donnà es, perte de
- bà nà fices,
- perte d
- 'exploitation, perte de clientèle ou de commandes,
-manque à gagner, trouble commercial quelconque) ou toute action dirigée
-contre le Licencié par un tiers, constitue un dommage indirect et
-n'ouvre
- pas droit rà paration par le Concà dant.
-
- Article 9 - GARANTIE
-
- 9.1 Le Licencià reconnaà t que l
- 'état actuel des connaissances
-scientifiques et techniques au moment de la mise en circulation du
-Logiciel ne permet pas d'en
- tester et d 'en vérifier toutes les
-utilisations ni de détecter l'existence d 'éventuels défauts. L'attention
- du Licencià a tà attirà e sur ce point sur les risques associàs au
- chargement,
- l
- 'utilisation, la modification et/ou au développement et Ã
-la reproduction du Logiciel qui sont réservés à des utilisateurs avertis.
-
-Il relève de la responsabilité du Licencié de contrôler, par tous
-moyens, l'adÃ
- quation du produit ses besoins, son bon fonctionnement et
- de s 'assurer qu'il ne causera pas de dommages aux personnes et aux biens.
-
- 9.2 Le Concà dant dà clare de bonne foi tre en droit de concà der
- l
- 'ensemble des droits attachés au Logiciel (comprenant notamment les
-droits visés à l'article 5).
-
-9.3 Le Licencià reconnaà t que le Logiciel est fourni "en l'état" par le
-Concà dant sans autre garantie, expresse ou tacite, que celle prà vue
-l
-'article 9.2 et notamment sans aucune garantie sur sa valeur commerciale,
-son caractère sécurisé, innovant ou pertinent.
-
-En particulier, le Concédant ne garantit pas que le Logiciel est exempt
-d'erreur ,
-qu 'il fonctionnera sans interruption, qu'il sera compatible
-avec l 'équipement du Licencié et sa configuration logicielle ni qu'il
-remplira les besoins du LicenciÃ.
-
-9.4 Le Concà dant ne garantit pas, de manià re expresse ou tacite, que le
-Logiciel ne porte pas atteinte un quelconque droit de proprià tÃ
-intellectuelle d
-'un tiers portant sur un brevet, un logiciel ou sur tout
-autre droit de propriété. Ainsi, le Concédant exclut toute garantie au
-profit du Licencié contre les actions en contrefaçon qui pourraient être
-diligentées au titre de l'utilisation ,
-de la modification, et de la
-redistribution du Logiciel.NÃ anmoins, si de telles actions sont
-exercà es contre le LicenciÃ, le Concà dant lui apportera son aide
-technique et juridique pour sa dà fense.Cette aide technique et
-juridique est dà terminà e au cas par cas entre le Concà dant concernà et
-le Licencià dans le cadre d 'un protocole d'accord.Le Concà dant dà gage
-toute responsabilità quant l 'utilisation de la dénomination du
-Logiciel par le Licencié. Aucune garantie n'est apportà e quant
-l 'existence de droits antérieurs sur le nom du Logiciel et sur
-l'existence d
-'une marque.
-
-
- Article 10 - RESILIATION
-
-10.1 En cas de manquement par le Licencié aux obligations mises à sa
-charge par le Contrat, le Concédant pourra résilier de plein droit le
-Contrat trente (30) jours après notification adressée au Licencié et
-restée sans effet.
-
-10.2 Le Licencié dont le Contrat est résilié n'est
-plus autorisÃ
-utiliser, modifier ou distribuer le Logiciel.Cependant, toutes les
-licences qu 'il aura concédées antérieurement à la résiliation du Contrat
-resteront valides sous réserve qu'elles aient tà effectuà es en
-conformità avec le Contrat.
-
-Article 11 - DISPOSITIONS DIVERSES
-
-11.1 CAUSE EXTERIEURE
-
-Aucune des Parties ne sera responsable d 'un retard ou d'une dà faillance
-d
-'exécution du Contrat qui serait dû à un cas de force majeure, un cas
-fortuit ou une cause extérieure, telle que, notamment, le mauvais
-fonctionnement ou les interruptions du réseau électrique ou de
-télécommunication, la paralysie du réseau liée à une attaque
-informatique, l'intervention
-des autoritàs gouvernementales, les
-catastrophes naturelles, les dà gà ts des eaux, les tremblements de terre,
-le feu, les explosions, les grà ves et les conflits sociaux, l 'état de
-guerre...
-
-11.2 Le fait, par l'une ou l 'autre des Parties, d'omettre en une ou
-plusieurs occasions de se prà valoir d
-'une ou plusieurs dispositions du
-Contrat, ne pourra en aucun cas impliquer renonciation par la Partie
-intéressée à s'en
-prà valoir ultà rieurement.
-
-11.3 Le Contrat annule et remplace toute convention antà rieure, crite
-ou orale,
-entre les Parties sur le mà me objet et constitue l
-'accord
-entier entre les Parties sur cet objet. Aucune addition ou modification
-aux termes du Contrat n'aura d
-'effet à l'à gard des Parties moins
-d 'être faite par écrit et signée par leurs représentants dûment habilités.
-
-11.4 Dans l'hypothà se où une ou plusieurs des dispositions du Contrat
-s
-'avèrerait contraire à une loi ou à un texte applicable, existants ou
-futurs, cette loi ou ce texte prévaudrait, et les Parties feraient les
-amendements nécessaires pour se conformer à cette loi ou à ce texte.
-Toutes les autres dispositions resteront en vigueur. De même, la
-nullité, pour quelque raison que ce soit, d'une
-des dispositions du
-Contrat ne saurait entraà ner la nullità de l
-'ensemble du Contrat.
-
-
- 11.5 LANGUE
-
-Le Contrat est rédigé en langue française et en langue anglaise, ces
-deux versions faisant également foi.
-
-
- Article 12 - NOUVELLES VERSIONS DU CONTRAT
-
-12.1 Toute personne est autorisée à copier et distribuer des copies de
-ce Contrat.
-
-12.2 Afin d'en
-prà server la cohà rence, le texte du Contrat est protà gÃ
-et ne peut tre modifià que par les auteurs de la licence, lesquels se
-rà servent le droit de publier pà riodiquement des mises jour ou de
-nouvelles versions du Contrat, qui possà deront chacune un numà ro
-distinct.Ces versions ultà rieures seront susceptibles de prendre en
-compte de nouvelles problà matiques rencontrà es par les logiciels libres.
-
-12.3 Tout Logiciel diffusà sous une version donnà e du Contrat ne pourra
-faire l 'objet d'une diffusion ultà rieure que sous la mà me version du
-Contrat ou une version postà rieure, sous rà serve des dispositions de
-l
-'article 5.3.4.
-
-
- Article 13 - LOI APPLICABLE ET COMPETENCE TERRITORIALE
-
-13.1 Le Contrat est régi par la loi française. Les Parties conviennent
-de tenter de régler à l'amiable
-les diffà rends ou litiges qui
-viendraient se produire par suite ou l 'occasion du Contrat.
-
-13.2 A défaut d'accord amiable dans un dà lai de deux(2) mois compter
-de leur survenance et sauf situation relevant d 'une procédure d'urgence,
-les diffà rends ou litiges seront portàs par la Partie la plus diligente
-devant les Tribunaux compà tents de Paris.
-
-Version 2.0 du 2006-0 9-05.
\ No newline at end of file
diff --git a/Code/Vahine/MaskUpdateFilter.h b/Code/Vahine/MaskUpdateFilter.h
deleted file mode 100644
index 4e8c11d515ec2601de81e7728433b7851dcc16dc..0000000000000000000000000000000000000000
--- a/Code/Vahine/MaskUpdateFilter.h
+++ /dev/null
@@ -1,80 +0,0 @@
-/*
- * MaskUpdateFilter.h
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2010 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __MASKUPDATEFILTER_H_
-#define __MASKUPDATEFILTER_H_
-
-#include <otb/IO/otbImage.h>
-#include <itkImageToImageFilter.h>
-#include <itkImageRegionIterator.h>
-#include <itkImageRegionConstIterator.h>
-#include "VahineMacro.h"
-#include "VahineException.h"
-
-namespace otb {
-template <class TMaskImage>
-class ITK_EXPORT VahineMaskUpdateFilter : public itk::ImageToImageFilter<TMaskImage, TMaskImage> {
-public:
- const static unsigned int Dimension = 2;
- typedef TMaskImage MaskImageType;
- typedef VahineMaskUpdateFilter Self;
- typedef itk::ImageToImageFilter<MaskImageType, MaskImageType> Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
- typedef typename MaskImageType::PixelType MaskImagePixelType;
- typedef typename MaskImageType::Pointer MaskImagePointer;
- typedef itk::ImageRegionConstIterator<MaskImageType> ConstMaskIteratorType;
- typedef itk::ImageRegionIterator<MaskImageType> MaskIteratorType;
- typedef typename MaskImageType::RegionType MaskRegionType;
- typedef typename MaskImageType::InternalPixelType MaskInternalPixelType;
- typedef Image<MaskInternalPixelType, Dimension> BandType;
- typedef typename BandType::PixelType BandPixelType;
- typedef itk::ImageRegionConstIterator<BandType> ConstBandIteratorType;
-
- // creation of SmartPointer
- itkNewMacro(Self);
- // runtime type information
- itkTypeMacro(VahineFlatReshapeFilter, itk::ImageToImageFilter);
-
- itkSetMacro( RemedialMask, MaskImagePointer );
- itkSetMacro( InputBand, unsigned int);
- itkSetMacro( RemedialBand, unsigned int);
- bool GetDebug(){
- return true;
- };
-protected:
- VahineMaskUpdateFilter();
- virtual ~VahineMaskUpdateFilter();
- void PrintSelf(std::ostream& os, itk::Indent indent) const;
-
- void GenerateData();
-
- //void ThreadedGenerateData(const BaseImageRegion& outputRegionForThread, int
- // threadId );
-private:
- VahineMaskUpdateFilter(const Self&); // not implemented
- void operator=(const Self&); //not implemented
-
- MaskImagePointer m_RemedialMask;
- unsigned int m_InputBand;
- unsigned int m_RemedialBand;
- double m_LowThreshold;
- double m_HighThreshold;
-};
-
-} // end namesapce otb
-
-#include "MaskUpdateFilter.txx"
-
-#endif /* MASKUPDATEFILTER_H_ */
\ No newline at end of file
diff --git a/Code/Vahine/MaskUpdateFilter.txx b/Code/Vahine/MaskUpdateFilter.txx
deleted file mode 100644
index 35b5f19e77cd79e25989b0da065f6a1d426b0498..0000000000000000000000000000000000000000
--- a/Code/Vahine/MaskUpdateFilter.txx
+++ /dev/null
@@ -1,81 +0,0 @@
-/*
- * MaskUpdateFilter.txx
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2010 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINEMASKUPDATEFILTER_TXX
-#define __VAHINEMASKUPDATEFILTER_TXX
-#include "MaskUpdateFilter.h"
-
-namespace otb {
-
-template<class TMaskImage>
-VahineMaskUpdateFilter<TMaskImage>::VahineMaskUpdateFilter() :
- m_LowThreshold(1.0), m_HighThreshold(1.0){
-}
-
-template<class TMaskImage>
-VahineMaskUpdateFilter<TMaskImage>::~VahineMaskUpdateFilter(){
-}
-
-template<class TMaskImage>
-void VahineMaskUpdateFilter<TMaskImage>::GenerateData(){
- const unsigned int inputMaskNumberOfbands = this->GetInput()->GetNumberOfComponentsPerPixel();
-
- m_RemedialMask->Update();
- MaskRegionType inputMaskRegion = this->GetInput()->GetLargestPossibleRegion();
- MaskRegionType remedialMaskRegion = m_RemedialMask->GetLargestPossibleRegion();
-
- ConstMaskIteratorType itRemedialMask(m_RemedialMask, remedialMaskRegion);
- ConstMaskIteratorType itInputMask(this->GetInput(), inputMaskRegion);
-
- MaskImagePointer maskOutput = MaskImageType::New();
- maskOutput->SetRegions(inputMaskRegion);
- maskOutput->SetNumberOfComponentsPerPixel(inputMaskNumberOfbands);
- if(m_InputBand>inputMaskNumberOfbands) {
- vahineDebug("Error on input band number");
- throw VahineException();
- }
- maskOutput->Allocate();
-
- MaskIteratorType itOutputMask(maskOutput, inputMaskRegion);
- itOutputMask.Begin();
- itRemedialMask.Begin();
- itInputMask.Begin();
- while(!itOutputMask.IsAtEnd() ) {
- MaskImagePixelType inputPixel = itInputMask.Get();
- MaskImagePixelType remedialValue = itRemedialMask.Get();
- //vahineDebug("in "<<inputPixel[m_InputBand]<<" rem
- // "<<remedialValue[m_RemedialBand]);
- if(m_LowThreshold <= inputPixel[m_InputBand] && inputPixel[m_InputBand] <= m_HighThreshold) {
- inputPixel[m_InputBand] = remedialValue[m_RemedialBand];
- ++itRemedialMask;
- }
- itOutputMask.Set(inputPixel);
- ++itInputMask;
- ++itOutputMask;
- }
-
- this->GraftOutput(maskOutput);
-}
-
-/**
- * Standard "PrintSelf" method
- */
-template<class TMaskImage>
-void VahineMaskUpdateFilter<TMaskImage>::PrintSelf(std::ostream& os, itk::Indent indent) const {
- Superclass::PrintSelf( os, indent );
-}
-
-} // end namespace otb
-
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/SyntheticBaseGeometricFilter.h b/Code/Vahine/SyntheticBaseGeometricFilter.h
deleted file mode 100644
index 8c51748f5a6ecfca970002a5b09caa06f82d0216..0000000000000000000000000000000000000000
--- a/Code/Vahine/SyntheticBaseGeometricFilter.h
+++ /dev/null
@@ -1,88 +0,0 @@
-/*
- * SyntheticBaseGeometricFilter.h
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINESYNTHETICBASEGEOMETRICFILTER_H_
-#define __VAHINESYNTHETICBASEGEOMETRICFILTER_H_
-
-#include "VahineMacro.h"
-#include "otbVectorImage.h"
-#include "itkImageToImageFilter.h"
-#include "itkImageRegionIterator.h"
-#include "itkImageLinearConstIteratorWithIndex.h"
-
-namespace otb {
-/**
- * \class VahineSyntheticBaseGeometricFilter
- * \brief Filter to select one or more synthetic database
- *
- * This class select one or more synthetic database.
- *
- * A synthetic database is a spectral cube were each column have
- * different geometry. With this filter you can select one geometry
- * i.e. one column.
- */
-template <class TImage>
-class ITK_EXPORT VahineSyntheticBaseGeometricFilter
- : public itk::ImageToImageFilter<TImage, TImage> {
-public:
- typedef TImage InputImageType;
- typedef TImage OutputImageType;
-
- typedef VahineSyntheticBaseGeometricFilter Self;
- typedef itk::ImageToImageFilter<InputImageType, InputImageType> Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
-
- // creation of SmartPointer
- itkNewMacro(Self);
- // runtime type information
- itkTypeMacro(VahineSyntheticFilter, itk::ImageToImageFilter);
-
- typedef itk::ImageLinearConstIteratorWithIndex<InputImageType> LinearIteratorType;
- typedef itk::ImageRegionIterator<OutputImageType> IteratorType;
- typedef typename OutputImageType::Pointer OutputImagePointer;
- typedef typename OutputImageType::IndexType OutputIndexType;
- typedef typename OutputImageType::SizeType OutputSizeType;
- typedef typename OutputImageType::RegionType OutputRegionType;
-
- itkSetMacro( ColumnIndex, unsigned int );
- itkGetMacro( ColumnIndex, unsigned int );
-
- bool GetDebug(){
- return m_Debug;
- };
- void SetDebug(bool debug){
- m_Debug = debug;
- };
-protected:
- VahineSyntheticBaseGeometricFilter();
- virtual ~VahineSyntheticBaseGeometricFilter();
- void PrintSelf(std::ostream& os, itk::Indent indent) const;
-
- void GenerateData();
-
-private:
- bool m_Debug;
- unsigned int m_ColumnIndex;
-
- VahineSyntheticBaseGeometricFilter(const Self&); //not implemented
- void operator=(const Self&); //not implemented
-
-};
-
-} // end namespace otb
-
-#include "SyntheticBaseGeometricFilter.txx"
-
-#endif /* __VAHINESYNTHETICBASEGEOMETRICFILTER_H_ */
\ No newline at end of file
diff --git a/Code/Vahine/SyntheticBaseGeometricFilter.txx b/Code/Vahine/SyntheticBaseGeometricFilter.txx
deleted file mode 100644
index a424ba84c795ce72a33193be247f7bb4bc104768..0000000000000000000000000000000000000000
--- a/Code/Vahine/SyntheticBaseGeometricFilter.txx
+++ /dev/null
@@ -1,80 +0,0 @@
-/*
- * SyntheticBaseGeometricFilter.txx
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-#ifndef __VAHINESYNTHETICBASEGEOMETRICFILTER_TXX
-#define __VAHINESYNTHETICBASEGEOMETRICFILTER_TXX
-#include "SyntheticBaseGeometricFilter.h"
-
-namespace otb {
-
-template<class TImage>
-VahineSyntheticBaseGeometricFilter<TImage>::VahineSyntheticBaseGeometricFilter() :
- m_Debug(false){
-}
-
-template<class TImage>
-VahineSyntheticBaseGeometricFilter<TImage>::~VahineSyntheticBaseGeometricFilter(){
-}
-
-template<class TImage>
-void VahineSyntheticBaseGeometricFilter<TImage>::GenerateData(){
- // size along X
- unsigned int cols = this->GetInput()->GetLargestPossibleRegion().GetSize(0);
- // size along Y
- unsigned int rows = this->GetInput()->GetLargestPossibleRegion().GetSize(1);
-
- if(m_ColumnIndex>=cols) {
- itkExceptionMacro("Invalid columns index");
- }
- const unsigned int InputNumberOfbands = this->GetInput()->GetNumberOfComponentsPerPixel();
-
- // create output image
- OutputImagePointer image = OutputImageType::New();
- OutputIndexType start;
- start[0] = 0;
- start[1] = 1;
- OutputSizeType size;
- size[0] = 1; // one columns
- size[1] = rows;
- OutputRegionType region;
- region.SetSize(size);
- region.SetIndex(start);
- image->SetRegions(region);
- image->SetNumberOfComponentsPerPixel(InputNumberOfbands);
- image->Allocate();
-
- LinearIteratorType inputIt(this->GetInput(), this->GetInput()->GetLargestPossibleRegion() );
- IteratorType outputIt(image, image->GetLargestPossibleRegion() );
- inputIt.SetDirection(1);
- unsigned i = 0;
- while(i<m_ColumnIndex) {
- inputIt.NextLine();
- ++i;
- }
- for(inputIt.GoToBeginOfLine(), outputIt.GoToBegin(); !inputIt.IsAtEndOfLine(); ++inputIt, ++outputIt) {
- outputIt.Set(inputIt.Get() );
- }
- this->GraftOutput(image);
-}
-
-/**
- * Standard "PrintSelf" method
- */
-template<class TImage>
-void VahineSyntheticBaseGeometricFilter<TImage>::PrintSelf(std::ostream& os, itk::Indent indent) const {
- Superclass::PrintSelf( os, indent );
-}
-
-} // end namespace otb
-
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/TikonovFilter.h b/Code/Vahine/TikonovFilter.h
deleted file mode 100644
index 191a03c385b973cad88944fbce1f7542b9800648..0000000000000000000000000000000000000000
--- a/Code/Vahine/TikonovFilter.h
+++ /dev/null
@@ -1,101 +0,0 @@
-/**
- * TikonovFilter.h
- *
- * Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-#ifndef __VAHINETIKONOVFILTER_H
-#define __VAHINETIKONOVFILTER_H
-
-#include "VahineMacro.h"
-
-#include "gdal.h"
-#include "otbImage.h"
-#include "otbVectorImage.h"
-#include "otbImageFileReader.h"
-#include "otbImageFileWriter.h"
-#include "otbStreamingImageFileWriter.h"
-#include <vnl/vnl_matrix.h>
-
-#include "itkVariableSizeMatrix.h"
-#include "itkImageToImageFilter.h"
-#include "otbStreamingStatisticsVectorImageFilter.h"
-
-namespace otb {
-template <class TInputImage>
-class ITK_EXPORT VahineTikonovFilter :
- public itk::ImageToImageFilter<TInputImage, TInputImage > {
-public:
- typedef VahineTikonovFilter Self;
- typedef itk::ImageToImageFilter<TInputImage, TInputImage> Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
-
- // creation of SmartPointer
- itkNewMacro(Self);
- // runtime type information
- itkTypeMacro(VahineTikonovFilter, ImageToImageFilter);
-
- /** Type to use for computation */
- typedef TInputImage InputImageType;
- typedef typename InputImageType::Pointer InputImagePointer;
- typedef itk::ImageRegionConstIterator<InputImageType> ConstIterator;
- typedef typename InputImageType::PixelType PixelType;
- typedef typename InputImageType::InternalPixelType InternalPixelType;
- typedef typename InputImageType::SizeType SizeType;
- typedef typename itk::NumericTraits<InternalPixelType>::RealType RealType;
- typedef itk::VariableLengthVector<RealType> RealPixelType;
-
- typedef typename itk::VariableSizeMatrix<RealType> MatrixType;
- typedef vnl_matrix<InternalPixelType> VnlMatrixType;
- typedef otb::StreamingStatisticsVectorImageFilter<TInputImage> StatsFilterType;
-
- itkStaticConstMacro(InputImageDimension, unsigned int, TInputImage::ImageDimension);
-
- itkSetMacro( RegularizationParam, double );
- MatrixType GetTikonovMatrix() {
- return this->m_tikonovMatrix;
- };
- MatrixType GetCovMatrix(){
- return this->m_covMatrix;
- };
- void UpdateCov();
-
- bool GetDebug(){
- return true;
- };
-protected:
- VahineTikonovFilter();
- virtual ~VahineTikonovFilter(){
- };
- virtual void PrintSelf(std::ostream& os, itk::Indent indent) const {
- Superclass::PrintSelf(os, indent);
- };
- void GenerateData();
-
- void SetRegularizationMatrix(MatrixType matrix);
-
- typename StatsFilterType::Pointer m_statsFilter;
-private:
- double m_RegularizationParam;
- MatrixType m_regulMatrix;
- MatrixType m_covMatrix;
- MatrixType m_tikonovMatrix;
- bool m_covUpdated;
-
- VahineTikonovFilter(const Self&); //not implemented
- void operator=(const Self&); //not implemented
-
-};
-} // End namespace otb
-
-#include "TikonovFilter.txx"
-
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/TikonovFilter.txx b/Code/Vahine/TikonovFilter.txx
deleted file mode 100644
index 75b3551cd81323645dff7838ebc9e7e6eb737057..0000000000000000000000000000000000000000
--- a/Code/Vahine/TikonovFilter.txx
+++ /dev/null
@@ -1,72 +0,0 @@
-/**
- * TikonovFilter.txx
- *
- * Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-#ifndef __VAHINETIKONOVFILTER_TXX
-#define __VAHINETIKONOVFILTER_TXX
-
-#include "TikonovFilter.h"
-
-namespace otb {
-
-template<class TInputImage>
-VahineTikonovFilter<TInputImage>::VahineTikonovFilter(){
- m_statsFilter = StatsFilterType::New();
- // default regularisation parameter
- m_RegularizationParam = 1E-1;
- m_covUpdated = false;
-}
-
-template<class TInputImage>
-void VahineTikonovFilter<TInputImage>::UpdateCov(){
- InputImagePointer image = const_cast<InputImageType *>(this->GetInput() );
-
- image->Update();
- SizeType size = image->GetLargestPossibleRegion().GetSize();
- //ConstIterator it(image, image->GetRequestedRegion());
-
- m_statsFilter->SetInput(image);
- m_statsFilter->Update();
- m_covMatrix = m_statsFilter->GetCovariance();
- unsigned int N = size[0]*size[1];
- m_covMatrix *= static_cast<double>(N)/(N-1);
- m_covUpdated = true;
-}
-
-template<class TInputImage>
-void VahineTikonovFilter<TInputImage>::GenerateData(){
- if(!m_covUpdated) {
- UpdateCov();
- }
- vahineDebug("regul = "<< m_RegularizationParam);
- vahineDebug("cov matrix size = " << m_covMatrix.Rows() <<"*"<< m_covMatrix.Cols() );
- int N = m_covMatrix.Cols();
- MatrixType reg, reg_t, cov_t;
- reg.SetSize(N, N);
- reg.SetIdentity();
- // TODO refactor trans(Identity) * Identity = Identity
- reg_t = reg.GetTranspose();
- reg = reg_t * reg;
- reg *= m_RegularizationParam;
-
- vahineDebug("reg size = "<<reg.Rows()<<"*"<<reg.Cols() );
-
- cov_t = m_covMatrix.GetTranspose();
-
- m_tikonovMatrix = (cov_t * m_covMatrix) + reg;
- m_tikonovMatrix = m_tikonovMatrix.GetInverse();
- m_tikonovMatrix = m_tikonovMatrix * cov_t;
- vahineDebug("sigma size = "<<m_tikonovMatrix.Rows()<<"*"<<m_tikonovMatrix.Cols() );
-}
-
-}
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/VahineException.h b/Code/Vahine/VahineException.h
deleted file mode 100644
index c604941913049c3287df67dc853608f3ab3f5930..0000000000000000000000000000000000000000
--- a/Code/Vahine/VahineException.h
+++ /dev/null
@@ -1,25 +0,0 @@
-/*
- * VahineException.h
- * Date: $Date$
- * Version: $Revision$
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2010 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef VAHINEEXCEPTION_H_
-#define VAHINEEXCEPTION_H_
-
-namespace otb {
-class GrsirInverseException { };
-class VahineException { };
-}
-
-#endif /* VAHINEEXCEPTION_H_ */
\ No newline at end of file
diff --git a/Code/Vahine/VahineInterpol.h b/Code/Vahine/VahineInterpol.h
deleted file mode 100644
index 1b470d09027a6723987ba4dbb1bfb911f679b5b4..0000000000000000000000000000000000000000
--- a/Code/Vahine/VahineInterpol.h
+++ /dev/null
@@ -1,73 +0,0 @@
-/*
- * VahineInterpol.h
- *
- * Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINEINTERPOL_H
-#define __VAHINEINTERPOL_H
-
-#include <vnl/vnl_vector_ref.h>
-#include "VahineMacro.h"
-
-namespace vahine {
-
-/**
- * Class of interpolation
- * For use it you can set \a Vx and \a Vy values and after call \a Update method.
- * You can get the interpolate value y of a input data x with the method \a Get(x)
- */
-template< typename RealType >
-class VahineInterpol {
-public:
- typedef vnl_vector<RealType> VectorType;
- VahineInterpol();
- VahineInterpol(unsigned int size);
- ~VahineInterpol();
- RealType get(RealType x);
-
- bool GetDebug(){
- return true;
- };
- std::string GetNameOfClass(){
- return "VahineInterpol";
- }
- void SetVx(VectorType vx);
-
- VectorType GetVx();
-
- void SetVy(VectorType vy);
-
- VectorType GetVy();
-
- void Update();
-
-private:
- unsigned int m_size;
- /**
- * vector of coefs for linear interpolation
- */
- VectorType m_aCoefs;
- VectorType m_vx;
- VectorType m_vy;
- void intervals(VectorType &v, VectorType &h);
-
- void linearCoeffs(VectorType &dx, VectorType &dy);
-
- void linear();
-
-};
-
-} // End namespace vahine
-
-#include "VahineInterpol.txx"
-
-#endif /* VAHINEINTERPOL_H_ */
\ No newline at end of file
diff --git a/Code/Vahine/VahineInterpol.txx b/Code/Vahine/VahineInterpol.txx
deleted file mode 100644
index 066ef55b9f2d2d8e4ddc18125feedcf761d53be5..0000000000000000000000000000000000000000
--- a/Code/Vahine/VahineInterpol.txx
+++ /dev/null
@@ -1,125 +0,0 @@
-/**
- * VahineInterpol.txx
- *
- * Interpolation class for vahine project
- * This class set in vector object several tuple of value (vx, vy)
- * and calculate a linear interpolation for each range
- *
- * Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINEINTERPOL_TXX
-#define __VAHINEINTERPOL_TXX
-
-#include "VahineInterpol.h"
-
-namespace vahine {
-template< typename RealType >
-VahineInterpol<RealType>::VahineInterpol(){
-}
-
-template< typename RealType >
-VahineInterpol<RealType>::VahineInterpol(unsigned int size) {
- //vahineDebug("Create interpol size = "<<size);
- m_size = size;
- m_vx = VectorType(size);
- m_vy = VectorType(size);
- m_aCoefs = VectorType(size);
-}
-
-template< typename RealType >
-VahineInterpol<RealType>::~VahineInterpol() {
-}
-
-template< typename RealType>
-void VahineInterpol<RealType>::SetVx(VectorType vx){
- m_vx = vx;
-}
-
-template< typename RealType>
-typename VahineInterpol<RealType>::VectorType VahineInterpol<RealType>::GetVx(){
- return m_vx;
-}
-
-template< typename RealType>
-void VahineInterpol<RealType>::SetVy(VectorType vy){
- m_vy = vy;
-}
-
-template< typename RealType>
-typename VahineInterpol<RealType>::VectorType VahineInterpol<RealType>::GetVy(){
- return m_vy;
-}
-
-/**
- * Calculate the interpolation with the current method
- * Actually only linear method are available
- */
-template < typename RealType>
-void VahineInterpol<RealType>::Update(){
- this->linear();
-}
-
-template< typename RealType >
-void
-VahineInterpol<RealType>::intervals(VectorType &v, VectorType &h){
- vahineDebug("intervals");
- for(int i=0; i<m_size-1; ++i)
- h[i] = v[i+1]-v[i];
- h[m_size-1] = 0.0; //to complete vector
-}
-
-template< typename RealType >
-void VahineInterpol<RealType>::linearCoeffs(VectorType &dx, VectorType &dy){
- for(int i = 0; i<m_size-1; i++)
- m_aCoefs[i] = dy[i]/dx[i];
-}
-
-/**
- * This method calculate all linear coef \a m_aCoefs for each range
- */
-template< typename RealType >
-void VahineInterpol<RealType>::linear(){
- VectorType dx = VectorType(m_size);
-
- this->intervals(m_vx, dx);
- VectorType dy = VectorType(m_size);
- this->intervals(m_vy, dy);
- this->linearCoeffs(dx, dy);
-}
-
-template<typename RealType>
-RealType VahineInterpol<RealType>::get(RealType x){
- // find appropriate segment
- //unsigned int n = m_vx.size();
- // interpolation outside the range is a constant value.
- if (x<m_vx[0]) {
- // TODO add warning
- //std::cout<<"interpolation inf range"<<m_vx[0]<<std::endl;
- return m_vy[0];
- }
- if (x>m_vx[m_size-1]) {
- // TODO add warning
- //std::cout<<"interpolation sup range "<<m_vx[n-1]<<std::endl;
- return m_vy[m_size-1];
- }
- // interpolation in the range
- int i;
- for(i=0; i<m_size; i++) if (x <m_vx[i]) break;
- // evaluate linear interpolation for i--
- i--;
- double t = m_aCoefs[i]*x + m_vy[i] - m_aCoefs[i]*m_vx[i];
- return t;
-}
-
-}
-
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/VahineMacro.h b/Code/Vahine/VahineMacro.h
deleted file mode 100644
index b6bb7c7183dcf661671cb2b8dff8f6acd3950a5b..0000000000000000000000000000000000000000
--- a/Code/Vahine/VahineMacro.h
+++ /dev/null
@@ -1,25 +0,0 @@
-/*
- * VahineMacro.h
- *
- * Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-#ifndef __VAHINE_MACRO_H
-#define __VAHINE_MACRO_H
-
-#define OTB_SHOW_ALL_MSG_DEBUG
-#include <otbMacro.h>
-
-#define vahineDebug(x) {otbDebugMacro(<<"[VAHINE] "<< x <<std::endl) }
-#define vahineDebugNel(x) {otbDebugMacro(<<"[VAHINE] "<< x) }
-
-#define vahineWarning(x) {itkWarningMacro(<<"[VAHINE] "<< x )}
-
-#endif /* __VAHINE_MACRO_H */
diff --git a/Code/Vahine/VahineSaxHandler.h b/Code/Vahine/VahineSaxHandler.h
deleted file mode 100644
index 22397afe8ac258a830ecbe8912f0fdc28e731348..0000000000000000000000000000000000000000
--- a/Code/Vahine/VahineSaxHandler.h
+++ /dev/null
@@ -1,181 +0,0 @@
-/*
- * VahineSaxHandler.h
- *
-* Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINESAXHANDLER_H
-#define __VAHINESAXHANDLER_H
-
-#include <xercesc/sax/HandlerBase.hpp>
-#include <xercesc/parsers/SAXParser.hpp>
-#include <xercesc/util/PlatformUtils.hpp>
-#include <xercesc/util/BinInputStream.hpp>
-#include <xercesc/sax/InputSource.hpp>
-
-#include <stack>
-#include <vector>
-#include <string>
-
-XERCES_CPP_NAMESPACE_USE
-
-namespace vahine {
-
-class XmlError {
-public:
- XmlError(){
- };
- XmlError(std::string m){
- msg =m;
- };
-private:
- std::string msg;
-}; // class exception
-
-class XmlTagWithChar {
-private:
- std::string value;
-public:
- XmlTagWithChar() {
- };
- void addChars(const char* buf, unsigned int len) {
- value += std::string(buf, len);
- }
-
- std::string getValue() {
- return value;
- }
-};
-
-class InverseParam : public XmlTagWithChar {
-private:
- bool activ;
-public:
- InverseParam(bool a){
- activ = a;
- }
- bool isActiv() {
- return activ;
- }
- std::string getDescription(){
- return getValue();
- }
-};
-
-class Band : public XmlTagWithChar {
-private:
- bool activ;
-public:
- Band(bool a) {
- activ = a;
- }
- bool isActiv() {
- return activ;
- }
- // convert string to double
- const double getLength() {
- return strtod(getValue().c_str(), NULL);
- }
- unsigned int number;
-};
-
-class WaveMask {
-public:
- typedef std::vector<Band> BandList;
- enum ApplyValue {
- UNKNOWN,
- IMAGE,
- LEARNDATA,
- };
-
- WaveMask() : isMaskDefaultOn(true) {
- }
- ~WaveMask() {
- };
- BandList bandList;
- bool isMaskDefaultOn;
- unsigned int numberOfBand;
- ApplyValue applyOn;
-
-};
-
-class SpatialMask : public XmlTagWithChar {
-public:
- unsigned int band;
- std::string getDescription(){
- return getValue();
- }
-};
-
-class GrsirParam {
-public:
- typedef std::vector<InverseParam> InverseParamList;
- typedef WaveMask::BandList BandList;
- typedef std::vector<double> RegulationList;
-
- GrsirParam() : regulation_auto(true){
- }
- ~GrsirParam(){
- }
-
- // file name list
- std::string image;
- std::string learning;
- std::string learning_noise;
- std::string learning_param;
- std::string image_param;
- std::string spatial_mask;
-
- InverseParamList inverse;
- RegulationList regulation;
- bool regulation_auto;
- SpatialMask spatialMask;
- WaveMask waveMask;
-};
-
-class VahineSAXHandler : public HandlerBase {
-private:
- bool charsValid;
- GrsirParam& grsir;
- std::string fileName;
- std::string tmp;
- // List of tags which can be found in the config XML file.
- enum ConfigTag {
- UNKNOWN,
- GRSIR,
- MASKWAVE,
- BAND,
- BANDS,
- FILE,
- INVERSE,
- PARAMETER,
- REGULATION,
- MASKSPATIAL
- };
- // The tag of the current read node.
- std::stack<ConfigTag> currentNodes;
- bool checkConfig();
-
-public:
- VahineSAXHandler(GrsirParam& g) : grsir(g), charsValid(false) {
- }
- void startElement(const XMLCh* const, AttributeList&);
-
- void endElement(const XMLCh* const name);
-
- void characters(const XMLCh* const buf, unsigned int len);
-
- void fatalError(const SAXParseException&);
-
-};
-
-} //end namespace vahine
-#endif /* VAHINESAXHANDLER_H_ */
\ No newline at end of file
diff --git a/Code/Vahine/VahineSlice.h b/Code/Vahine/VahineSlice.h
deleted file mode 100644
index 2e899fbdd413df93741990752498ea54f8db3ce7..0000000000000000000000000000000000000000
--- a/Code/Vahine/VahineSlice.h
+++ /dev/null
@@ -1,94 +0,0 @@
-/*
- * VahineSlice.h
- *
- * Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINESLICE_H
-#define __VAHINESLICE_H
-
-#include "otbImage.h"
-#include "itkVariableLengthVector.h"
-#include "itkImageRegionIterator.h"
-
-namespace vahine {
-/** Slice Class :
- * nPixel, count the number of observation in the slice
- * index, list of all image index in the current slice
- * sumPixel, sum of centering data foreach band
- * sumParam, sum of parameter value for all index in current slice
- * meanParam,
- * meanData,
- */
-template<class TReal>
-class Slice {
-public:
- typedef TReal RealType;
- typedef typename itk::VariableLengthVector<RealType> PixelType;
- typedef typename otb::Image<RealType, 2> ImageType;
- typedef typename ImageType::IndexType ImageIndexType;
- typedef std::list<ImageIndexType> ImageIndexListType;
- typedef typename ImageType::Pointer ProjPointer;
- typedef typename itk::ImageRegionIterator<ImageType> IteratorProjType;
-
- Slice();
- Slice(const Slice& s);
- Slice(ImageIndexType idx, PixelType pixel, RealType param);
- ~Slice(){
- };
- Slice& operator=(const Slice& s);
-
- bool operator<(const Slice& s);
-
- void add(ImageIndexType idx, PixelType pixel, RealType param);
-
- void calculateMean(ProjPointer projection);
-
- RealType getMeanParam();
-
- RealType getMeanData();
-
- PixelType getMeanPixel();
-
- bool GetDebug(){
- return true;
- }
- std::string GetNameOfClass(){
- return "Slice";
- }
- void setMeanPixel();
-
- void setPropOfPixel(unsigned int nObserv);
-
- double getPropOfPixel();
-
-private:
-
- int nPixel;
- PixelType sumPixel;
- PixelType meanPixel;
- ImageIndexListType index;
- RealType sumParam;
- RealType meanParam;
- RealType meanData;
- double propOfPixel;
-
- void setMeanParam();
-
- void setMeanData(ProjPointer projection);
-
-};
-
-} //end namespace vahine
-
-#include "VahineSlice.txx"
-
-#endif /* VAHINESLICE_H_ */
\ No newline at end of file
diff --git a/Code/Vahine/VahineSlice.txx b/Code/Vahine/VahineSlice.txx
deleted file mode 100644
index b726211d4bfae38d6ed887d743768aeb9b99977a..0000000000000000000000000000000000000000
--- a/Code/Vahine/VahineSlice.txx
+++ /dev/null
@@ -1,135 +0,0 @@
-/*
- * VahineSlice.txx
- *
-* Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-#ifndef __VAHINESLICE_TXX
-#define __VAHINESLICE_TXX
-
-#include "VahineSlice.h"
-
-using namespace vahine;
-using namespace std;
-
-template<class TReal>
-vahine::Slice<TReal>::Slice(){
- nPixel = 0;
- meanParam = 0;
- meanData = 0;
- sumPixel = PixelType();
- sumParam = 0;
- meanPixel = PixelType();
-}
-
-template<class TReal>
-vahine::Slice<TReal>::Slice(const Slice& s){
- nPixel = s.nPixel;
- index = s.index;
- sumPixel = s.sumPixel;
- sumParam = s.sumParam;
- meanParam = s.meanParam;
- meanData = s.meanData;
-}
-
-template<class TReal>
-vahine::Slice<TReal>::Slice(ImageIndexType idx, PixelType pixel, RealType param){
- nPixel = 1;
- index.push_back(idx);
- //sumPixel.SetSize(pixel.GetSize());
- sumPixel = pixel;
- sumParam = param;
- meanParam = 0;
- meanData = 0;
- propOfPixel = 0;
- //meanPixel.SetSize(pixel.GetSize());
-}
-
-template<class TReal>
-vahine::Slice<TReal>& vahine::Slice<TReal>::operator=(const Slice& s){
- nPixel = s.nPixel;
- sumPixel = s.sumPixel;
- meanPixel = s.meanPixel;
- index = s.index;
- propOfPixel = s.propOfPixel;
-
- sumParam = s.sumParam;
- meanParam = s.meanParam;
- meanData = s.meanData;
-}
-
-template<class TReal>
-bool vahine::Slice<TReal>::operator<(const Slice& s){
- return meanData < s.meanData;
-}
-
-template<class TReal>
-void vahine::Slice<TReal>::add(ImageIndexType idx, PixelType pixel, RealType param){
- nPixel++;
- index.push_back(idx);
- sumPixel += pixel;
- sumParam += param;
-}
-
-template<class TReal>
-void vahine::Slice<TReal>::calculateMean(ProjPointer projection){
- setMeanParam();
- setMeanData(projection);
-}
-
-template<class TReal>
-typename vahine::Slice<TReal>::RealType vahine::Slice<TReal>::getMeanParam(){
- return meanParam;
-}
-
-template<class TReal>
-typename Slice<TReal>::RealType Slice<TReal>::getMeanData(){
- return meanData;
-}
-
-template<class TReal>
-void Slice<TReal>::setMeanParam(){
- meanParam = sumParam / static_cast<RealType>(nPixel);
-}
-
-template<class TReal>
-void Slice<TReal>::setPropOfPixel(unsigned int nObserv){
- propOfPixel = nPixel / static_cast<double>(nObserv);
-}
-
-template<class TReal>
-double Slice<TReal>::getPropOfPixel(){
- return propOfPixel;
-}
-
-template<class TReal>
-void Slice<TReal>::setMeanPixel(){
- meanPixel = sumPixel / static_cast<RealType>(nPixel);
-}
-
-template<class TReal>
-typename Slice<TReal>::PixelType Slice<TReal>::getMeanPixel(){
- return meanPixel;
-}
-
-template<class TReal>
-void Slice<TReal>::setMeanData(ProjPointer projection){
- typename ImageIndexListType::iterator itImageIndex = index.begin();
- IteratorProjType itProj(projection, projection->GetLargestPossibleRegion() );
- RealType sumData = static_cast<RealType>(0);
- while( itImageIndex != index.end() ) {
- itProj.SetIndex(*itImageIndex);
- sumData += itProj.Get();
- ++itImageIndex;
- }
- meanData = sumData / static_cast<RealType>(nPixel);
-}
-
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/VahineXml.h b/Code/Vahine/VahineXml.h
deleted file mode 100644
index 7e9300384025593dc67b53eef95a35282a76869c..0000000000000000000000000000000000000000
--- a/Code/Vahine/VahineXml.h
+++ /dev/null
@@ -1,65 +0,0 @@
-/*
- * VahineXml.h
- *
-* Vahine Framework
- * Copyright (C) 2008, 2009 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINEXML_H
-#define __VAHINEXML_H
-
-#include <xercesc/dom/DOM.hpp>
-#include <xercesc/framework/StdOutFormatTarget.hpp>
-#include <xercesc/framework/LocalFileFormatTarget.hpp>
-
-#include <string>
-#include <list>
-
-XERCES_CPP_NAMESPACE_USE
-
-namespace vahine {
-
-class XmlLog {
-public:
- XmlLog(std::string moduleName);
- ~XmlLog();
-
- typedef DOMElement* DOMElementPointer;
- typedef DOMDocument* DOMDocumentPointer;
- typedef DOMImplementation* DOMImplementationPointer;
-
- void CreateElement(std::string name);
-
- void AppendAllToRoot();
-
- void AppendToSup();
-
- void CreateTextNode(std::string text);
-
- void CreateElementTextApp(std::string name, std::string text);
-
- void CreateAttribute(std::string name, std::string value);
-
- void SetFileName(std::string fileName);
-
- void Write();
-
-private:
-
- std::string m_XmlFileName;
- DOMImplementationPointer m_Implementation;
- DOMDocumentPointer m_Document;
- DOMElementPointer m_RootElement;
- std::list<DOMElementPointer> m_list;
-};
-
-} // end namespace vahine
-
-#endif
\ No newline at end of file
diff --git a/Code/Vahine/VcaFilter.h b/Code/Vahine/VcaFilter.h
deleted file mode 100644
index 0cbc5a0930f62bf37fb317bcdca84ce92187a797..0000000000000000000000000000000000000000
--- a/Code/Vahine/VcaFilter.h
+++ /dev/null
@@ -1,150 +0,0 @@
-/*
- * VCAFilter.h
- * Date: $Date$
- * Version: $Revision$
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2010 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef VCAFILTER_H_
-#define VCAFILTER_H_
-
-#include <vnl/vnl_matrix.h>
-#include <vnl/vnl_random.h>
-#include <vnl/algo/vnl_svd.h>
-#include <vnl/algo/vnl_matrix_inverse.h>
-
-#include <otbVectorImage.h>
-#include <itkArray.h>
-#include <itkImageToImageFilter.h>
-#include <itkImageRegionIterator.h>
-#include <itkImageRegionConstIterator.h>
-#include <itkVariableLengthVector.h>
-
-#include "VahineMacro.h"
-
-namespace otb {
-
-template <class TImage>
-class ITK_EXPORT VahineVCAFilter
- : public itk::ImageToImageFilter<TImage, TImage> {
-public:
- typedef TImage VectorImageType;
-
- typedef VahineVCAFilter Self;
- typedef itk::ImageToImageFilter<VectorImageType, VectorImageType> Superclass;
- typedef itk::SmartPointer<Self> Pointer;
- typedef itk::SmartPointer<const Self> ConstPointer;
-
- // creation of SmartPointer
- itkNewMacro(Self);
- // runtime type information
- itkTypeMacro(VahineVCAFilter, itk::ImageToImageFilter);
-
- // Input/Output Image
- typedef typename Superclass::InputImageConstPointer ImageConstPointer;
- typedef typename VectorImageType::Pointer VectorImagePointer;
- typedef typename VectorImageType::PixelType PixelType;
- typedef typename VectorImageType::RegionType RegionType;
- typedef typename VectorImageType::SizeType SizeType;
- typedef typename VectorImageType::IndexType IndexType;
- typedef itk::ImageRegionConstIterator<VectorImageType> ConstIteratorType;
- typedef itk::ImageRegionIterator<VectorImageType> IteratorType;
- typedef typename VectorImageType::InternalPixelType InternalPixelType;
-
- // type for computation
- typedef vnl_matrix<double> VnlMatrix;
- typedef vnl_vector<double> VnlVector;
- typedef typename std::vector<IndexType> IndexArray;
- typedef VnlMatrix::iterator VnlMatrixIterator;
-
- bool GetDebug(){
- return true;
- };
- itkSetMacro( NumberEndMember, unsigned int );
- itkGetMacro( NumberEndMember, unsigned int );
- itkGetMacro( NumberOfBand, unsigned int);
- itkSetMacro( SNR, double );
- itkGetMacro( SNR, double );
- itkGetMacro( MixingMatrix, VnlMatrix );
- itkGetMacro( ProjectedData, VnlMatrix );
- itkSetMacro( Y, VnlMatrix );
- itkGetMacro( Y, VnlMatrix );
- IndexArray GetPurePixelsIndex();
-
- void GenerateOutputInformation();
-
- // for testing
- friend class VCAFilterFixture;
-protected:
- VahineVCAFilter();
- virtual ~VahineVCAFilter();
-
- void PrintSelf(std::ostream& os, itk::Indent indent) const;
-
- void GenerateData();
-
- unsigned int m_NumberOfBand; // protected for testing
- SizeType m_InputSize; // protected for testing
-private:
- // type for computation
- typedef typename std::vector<unsigned int> IndiceArray;
-
- void FillRandom(VnlVector & w);
-
- void VCA(VnlMatrix data);
-
- double SNR(VnlMatrix & M);
-
- double EstimateSNR(VnlMatrix & M, VnlVector mean, VnlMatrix & X);
-
- double TheoricalSNR(unsigned int p);
-
- void MeanCalculation(VnlMatrix const& M, VnlMatrix & MMean, VnlMatrix & Centered, VnlMatrix & Result,
- VnlVector & mean);
-
- void ProjectivProjection(VnlMatrix const& M, VnlMatrix & Rp);
-
- void SubspaceProjection(VnlMatrix const& M, VnlMatrix & Rp);
-
- void SingularValueAndVector(VnlMatrix const& M, VnlMatrix & Ud, unsigned int p);
-
- void UpdateIndex();
-
- unsigned int arg_max(const VnlVector &v);
-
- unsigned int m_NumberEndMember;
- double m_SNR;
- VnlMatrix m_MixingMatrix, m_Y, m_ProjectedData;
- IndiceArray m_PurePixelsIndices;
- IndexArray m_PurePixelsIndex;
- vnl_random m_RandomGenerator;
-
- // TODO flag to speed calculation
- // if false only calculate mixing matrix.
- //bool m_ProjectData;
-
- VahineVCAFilter(const Self&); //not implemented
- void operator=(const Self&); //not implemented
-
- static double mysquare(const double x){
- return x*x;
- }
- static double myabs(const double x){
- return ( (x<0) ? x* -1 : x);
- }
-};
-
-} // end namesapce otb
-
-#include "VcaFilter.txx"
-
-#endif /* VCAFILTER_H_ */
diff --git a/Code/Vahine/VcaFilter.txx b/Code/Vahine/VcaFilter.txx
deleted file mode 100644
index 94de75f0a29c52d90dd560efb097ef4753e7856a..0000000000000000000000000000000000000000
--- a/Code/Vahine/VcaFilter.txx
+++ /dev/null
@@ -1,407 +0,0 @@
-/*
- * VCAFilter.txx
- * Date: $Date$
- * Version: $Revision$
- *
- * Vahine Framework
- * Copyright (C) 2008 to 2010 Ludovic Léau-Mercier and Laboratoire de Planétologie de Grenoble
- * See LICENCE and COPYING for details.
- *
- * 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
- * CeCILLl License for more details and http://www.cecill.info
- *
- */
-
-#ifndef __VAHINEVCAFILTER_TXX
-#define __VAHINEVCAFILTER_TXX
-#include "VcaFilter.h"
-
-namespace otb {
-
-template<class TImage>
-VahineVCAFilter<TImage>::VahineVCAFilter() :
- m_SNR(NAN),
- m_MixingMatrix(),
- m_PurePixelsIndices(),
- m_PurePixelsIndex(),
- m_ProjectedData(),
- m_Y(),
- m_NumberEndMember(0) {
- m_RandomGenerator.reseed();
-}
-
-template<class TImage>
-VahineVCAFilter<TImage>::~VahineVCAFilter(){
- /*m_Y.clear();
- m_ProjectedData.clear();
- m_MixingMatrix.clear();*/
-}
-
-template<class TImage>
-void VahineVCAFilter<TImage>::GenerateOutputInformation(){
- // call the superclass' implementation of this method
- Superclass::GenerateOutputInformation();
-
- // get pointers to the input and output
- ImageConstPointer inputPtr = this->GetInput();
- VectorImagePointer outputPtr = this->GetOutput();
-
- if ( !inputPtr || !outputPtr ) {
- return;
- }
- unsigned int numberOfBands = inputPtr->GetNumberOfComponentsPerPixel();
- RegionType imageRegion = this->GetInput()->GetLargestPossibleRegion();
- outputPtr->SetLargestPossibleRegion( imageRegion );
- outputPtr->SetNumberOfComponentsPerPixel( numberOfBands );
-}
-
-template<class TImage>
-void VahineVCAFilter<TImage>::GenerateData(){
-
- RegionType imageRegion = this->GetInput()->GetLargestPossibleRegion();
-
- m_InputSize = imageRegion.GetSize();
- m_NumberOfBand = this->GetInput()->GetNumberOfComponentsPerPixel();
-
- if(m_NumberEndMember < 1 || m_NumberEndMember > m_NumberOfBand) {
- itkExceptionMacro("Endmember parameter must be integer between 1 and number of band of input image");
- }
-
- VectorImagePointer output = VectorImageType::New();
-
- // fill the data matrix from the input image
- // TODO add threading for this
- ConstIteratorType itInput( this->GetInput(), imageRegion );
- vahineDebug("input band="<<m_NumberOfBand<<" size[0]="<<m_InputSize[0]<<" size[1]="<<m_InputSize[1]);
- VnlMatrix data(m_NumberOfBand, m_InputSize[0]*m_InputSize[1]);
-
- unsigned int idx(0);
- for(itInput.GoToBegin(); !itInput.IsAtEnd(); ++itInput) {
- PixelType pixel = itInput.Get();
- VnlVector column(m_NumberOfBand);
- for(unsigned int j(0); j<m_NumberOfBand; ++j) {
- column[j] = static_cast<double>(pixel[j]);
- }
- data.set_column(idx, column);
- ++idx;
- }
- //vahineDebug("data = "<<data);
-
- // run calculation
- VCA(data);
-
- // convert projected data to image for output
- // TODO add threading for this
- vahineDebug("start output");
- //VectorImagePointer output = VectorImageType::New();
- output->SetRegions(imageRegion);
- output->SetNumberOfComponentsPerPixel(m_NumberOfBand);
- output->Allocate();
- IteratorType itOutput( output, imageRegion );
- idx = 0;
- for(itOutput.GoToBegin(); !itOutput.IsAtEnd(); ++itOutput) {
- VnlVector column = m_ProjectedData.get_column(idx);
- PixelType pixel = itOutput.Get();
- for(unsigned int j(0); j<m_NumberOfBand; ++j) {
- pixel[j] = static_cast<InternalPixelType>(column(j) );
- }
- itOutput.Set(pixel);
- ++idx;
- }
- vahineDebug("end output");
- this->GraftOutput(output);
-}
-
-template<class TImage>
-void VahineVCAFilter<TImage>::FillRandom(VnlVector & w){
- VnlVector::iterator it = w.begin();
-
- while(it != w.end() ) {
- (*it) = m_RandomGenerator.drand32(1.0);
- ++it;
- }
-}
-
-template<class TImage>
-void VahineVCAFilter<TImage>::VCA(VnlMatrix data){
- double snr, snr_th;
-
- //VnlMatrix Y;
- //m_PurePixelsIndices.SetSize(m_NumberEndMember);
- m_MixingMatrix.set_size(m_NumberOfBand, m_NumberEndMember);
- //m_PurePixelsIndices.resize(m_NumberEndMember, 0);
-
- if(isnan(m_SNR) ) {
- // SNR Estimates
- snr = SNR(data);
- }
- snr_th = TheoricalSNR(m_NumberEndMember);
- vahineDebug("Snr estimated or set = "<<snr);
- vahineDebug("Snr theorical = "<<snr_th);
- if(snr < snr_th) {
- vahineDebug("Select the projectiv projection");
- ProjectivProjection(data, m_ProjectedData);
- } else {
- vahineDebug("Select projection to subspace");
- SubspaceProjection(data, m_ProjectedData);
- }
-
- // VCA Algorithm
- vahineDebug("start VCA algorithm");
- VnlMatrix A(m_NumberEndMember, m_NumberEndMember, 0);
- A(m_NumberEndMember-1, 0) = 1;
- VnlVector w(m_NumberEndMember);
- VnlVector f, fsquare, v, col;
- for(unsigned int i=0; i<m_NumberEndMember; i++) {
- FillRandom(w);
- vnl_svd<double> Asvd(A);
- VnlMatrix Ainv = Asvd.pinverse();
- f = w - A*Ainv*w;
- fsquare = f.apply(&mysquare);
- f = f / sqrt(fsquare.sum() );
- v = f*m_Y;
- v=v.apply(&myabs);
- // search pure pixel
- unsigned int j = arg_max(v);
- col = m_Y.get_column(j); // save Y columns values
- m_PurePixelsIndices.push_back(j); // save indice number
- A.set_column(i, col);
- }
-
- for(unsigned int i=0; i<m_NumberEndMember; ++i) {
- vahineDebug("pure pixel="<<m_PurePixelsIndices[i]);
- m_MixingMatrix.set_column(i, m_ProjectedData.get_column(m_PurePixelsIndices[i]) );
- }
-
- UpdateIndex();
- vahineDebug("end VCA method");
-}
-
-/**
- * deprecated with last vnl_vector : use vnl_vector::arg_max
- */
-template<class TImage>
-unsigned int VahineVCAFilter<TImage>::arg_max(const VnlVector &v){
- unsigned int n = v.size();
-
- if (n==0) return unsigned(-1); // the minimum of an empty set is undefined
- double tmp = v[0];
- unsigned int idx = 0;
- for (unsigned int i=1; i<n; ++i)
- if (v[i] > tmp) {
- tmp = v[i];
- idx = i;
- }
- return idx;
-}
-
-/**
- * update index of pure pixel vector with indice of each pure pixel
- * in matrix calculation
- */
-template<class TImage>
-void VahineVCAFilter<TImage>::UpdateIndex(){
- vahineDebug("Update index");
- //m_PurePixelsIndex.resize(2,IndexType());
- for(unsigned int i = 0; i < m_NumberEndMember; ++i) {
- unsigned int curIndice = m_PurePixelsIndices[i];
- IndexType pixelIndex;
- pixelIndex[0] = curIndice % m_InputSize[0];
- pixelIndex[1] = static_cast<unsigned int>(curIndice/m_InputSize[0]);
- vahineDebug("pixelIndex "<<pixelIndex);
- m_PurePixelsIndex.push_back(pixelIndex);
- }
-}
-
-template<class TImage>
-double VahineVCAFilter<TImage>::SNR(VnlMatrix & M){
- VnlMatrix Ud, Centered, Xp, MMean, Result;
- VnlVector mean;
-
- MeanCalculation(M, MMean, Centered, Result, mean);
- SingularValueAndVector(Result, Ud, m_NumberEndMember); // compute the
- // p-projection
- // matrix
- //vahineDebug("Ud="<<Ud);
- //vahineDebug("Centered="<<Centered);
- Xp = Ud.transpose() * Centered; // project thezeros mean data onto
- // p-subspace
- //vahineDebug("Xp="<<Xp);
- double snr = EstimateSNR(M, mean, Xp);
- //vahineDebug("snr = "<<snr);
- return snr;
-}
-
-/**
- * SNR Estimation
- * L number of bands
- * N number of pixels
- * p number of endmembers
- */
-template<class TImage>
-double VahineVCAFilter<TImage>::EstimateSNR(VnlMatrix & M, VnlVector mean, VnlMatrix & X){
- unsigned int L = M.rows();
- unsigned int N = M.cols();
- unsigned int p = X.rows();
- double px(0.0), py(0.0);
- VnlMatrixIterator itM = M.begin();
-
- while(itM != M.end() ) {
- py += pow( (*itM),2);
- ++itM;
- }
- py /= N;
- //TODO with apply() square and sum
- VnlMatrixIterator itX = X.begin();
- //X.apply(&mysquare);
- while(itX != X.end() ) {
- px += mysquare(*itX);
- ++itX;
- }
- px /= N;
- px += dot_product(mean,mean);
- //vahineDebug("px="<<px<<" py="<<py<<" p="<<p<<" L="<<L);
- return 10*log10( (px-p/static_cast<double>(L)*py)/(py-px) );
-}
-
-template<class TImage>
-double VahineVCAFilter<TImage>::TheoricalSNR(unsigned int p){
- return 15+10*log10(p);
-}
-
-template<class TImage>
-void VahineVCAFilter<TImage>::MeanCalculation(VnlMatrix const& M, VnlMatrix& MMean, VnlMatrix & Centered,
- VnlMatrix & Result,
- VnlVector & mean){
- unsigned int col = M.cols();
- unsigned int row = M.rows();
-
- // row mean calculation
- mean.set_size(row);
- for(unsigned int i = 0; i<row; ++i) {
- VnlVector currentRow = M.get_row(i);
- mean(i) = currentRow.mean();
- }
- MMean.set_size(row,col);
- for(unsigned int i = 0; i < col; ++i) {
- MMean.set_column(i,mean);
- }
- Centered.set_size(row, col);
- Centered = M;
- Centered -= MMean;
- Result = Centered * Centered.transpose();
- Result /= col;
-}
-
-template<class TImage>
-void VahineVCAFilter<TImage>::ProjectivProjection(VnlMatrix const& M, VnlMatrix & Rp){
- VnlMatrix Ud, Centered, Xp, MMean, Result;
- VnlVector mean;
- unsigned int d = m_NumberEndMember - 1;
- unsigned int cols = M.cols();
-
- MeanCalculation(M, MMean, Centered, Result, mean);
- SingularValueAndVector(Result, Ud, d); // compute the p-projection matrix
- Xp = Ud.transpose() * Centered; // project thezeros mean data onto
- // p-subspace
- Rp = Ud * Xp + MMean; // again in dimension L see matlab
- // script
- //double(*f)(double) = □
- VnlMatrix XpSquare = Xp.apply(&mysquare); // Square root all elements
- //unsigned int cols = ProjASq.cols();
- VnlVector Sum(cols);
- for(unsigned int i =0; i<cols; ++i) {
- Sum(i) = (XpSquare.get_column(i) ).sum();
- }
- double c = sqrt(Sum.max_value() );
- m_Y.set_size(m_NumberEndMember,cols);
- for(unsigned int i=0; i<d; ++i) {
- m_Y.set_row(i,VnlVector(Xp.get_row(i) ) );
- }
- m_Y.set_row(d,VnlVector(cols,c) );
-}
-
-template<class TImage>
-void VahineVCAFilter<TImage>::SubspaceProjection(VnlMatrix const& M, VnlMatrix & Rp){
- VnlMatrix Result, Ud, Xp, MMean, Mult, Temp;
- unsigned int d = m_NumberEndMember;
- unsigned int cols = M.cols();
- unsigned int rows = M.rows();
-
- Result = M*M.transpose();
- Result /= cols;
- SingularValueAndVector(Result, Ud, d); // compute the p-projection matrix
- Xp = Ud.transpose()*M;
- Rp = Ud * Xp; // again in dimension L see matlab script
- // row mean calculation
- VnlVector mean(d);
- for(unsigned int i = 0; i<d; ++i) {
- VnlVector currentRow = Xp.get_row(i);
- mean(i) = currentRow.mean();
- }
- MMean.set_size(d,cols);
- for(unsigned int i = 0; i < cols; ++i) {
- MMean.set_column(i,mean);
- }
- // element by element multiplication
- Mult.set_size(d, cols);
- VnlMatrixIterator itMMean = MMean.begin();
- VnlMatrixIterator itXp = Xp.begin();
- VnlMatrixIterator itMult = Mult.begin();
- while(itMMean != MMean.end() ) {
- (*itMult) = (*itMMean) * (*itXp);
- ++itMMean;
- ++itXp;
- ++itMult;
- }
- VnlVector sums(cols);
- for(unsigned int i=0; i< cols; ++i) {
- sums(i) = Mult.get_column(i).sum();
- }
- Temp.set_size(d,cols);
- for(unsigned int i=0; i<d; ++i) {
- Temp.set_row(i,sums);
- }
- m_Y.set_size(d,cols);
- VnlMatrixIterator itTemp = Temp.begin();
- VnlMatrixIterator itY = m_Y.begin();
- itXp = Xp.begin();
- while(itTemp != Temp.end() ) {
- *itY = (*itXp)/(*itTemp);
- ++itXp;
- ++itTemp;
- ++itY;
- }
-}
-
-template<class TImage>
-void VahineVCAFilter<TImage>::SingularValueAndVector(VnlMatrix const& M, VnlMatrix & Ud, unsigned int p){
- vnl_svd<double> svd(M);
- VnlMatrix fullUd = svd.U();
- unsigned int row = fullUd.rows();
- unsigned int col = fullUd.cols();
- // keep p columns
- Ud.set_size(row,p);
- for(unsigned int i=0; i<p; ++i) {
- Ud.set_column(i,fullUd.get_column(i) );
- }
-}
-
-template<class TImage>
-typename VahineVCAFilter<TImage>::IndexArray VahineVCAFilter<TImage>::GetPurePixelsIndex(){
- return m_PurePixelsIndex;
-}
-
-/**
- * Standard "PrintSelf" method
- */
-template<class TImage>
-void VahineVCAFilter<TImage>::PrintSelf(std::ostream& os, itk::Indent indent) const {
- Superclass::PrintSelf( os, indent );
-}
-
-} //end namespace otb
-
-#endif
\ No newline at end of file