Posts tagged "otb":

Installing OTB has never been so easy

You've heard about the Orfeo Toolbox library and its wonders, but urban legends say that it is difficult to install. Don't believe that. Maybe it was difficult to install, but this is not the case anymore.

Thanks to the heroic work of the OTB core development team, installing OTB has never been so easy. In this post, you will find the step-by-step procedure to compile OTB from source on a Debian 8.0 Jessie GNU/Linux distribution.

su -

visudo

# User privilege specification
root    ALL=(ALL:ALL) ALL

otb     ALL=(ALL:ALL) ALL

sudo aptitude install mercurial \
     cmake-curses-gui build-essential \
     qt4-dev-tools libqt4-core \
     libqt4-dev libboost1.55-dev \
     zlib1g-dev libopencv-dev curl \
     libcurl4-openssl-dev swig \
     libpython-dev 

mkdir OTB
cd OTB

hg clone http://hg.orfeo-toolbox.org/OTB

cd OTB
hg update 5.0.0-rc1
cd ../

wget https://www.orfeo-toolbox.org/packages/SuperBuild-archives.tar.bz2
tar xvjf SuperBuild-archives.tar.bz2

mkdir -p SuperBuild/OTB

cd SuperBuild/OTB

cmake \
    -DCMAKE_CXX_FLAGS:STRING=-std=c++14 \
    -DOTB_WRAP_PYTHON:BOOL=ON \
    -DCMAKE_BUILD_TYPE:STRING=Release \
    -DCMAKE_INSTALL_PREFIX:PATH=/home/otb/OTB/SuperBuild/ \
    -DDOWNLOAD_LOCATION:PATH=/home/otb/OTB/SuperBuild-archives/ \
    -DOTB_USE_OPENJPEG:BOOL=ON \
    -DUSE_SYSTEM_OPENJPEG:BOOL=OFF \
    ../../OTB/SuperBuild/

make -j4

export OTB_HOME=${HOME}/OTB/SuperBuild
export PATH=${OTB_HOME}/bin:$PATH
export LD_LIBRARY_PATH=${OTB_HOME}/lib

cd
source .bashrc 

otbcli_BandMath

otbApplicationLauncherCommandLine BandMath $OTB_HOME/lib/otb/applications/

A simple thread pool to run parallel PROSail simulations

In the otb-bv we use the OTB versions of the Prospect and Sail models to perform satellite reflectance simulations of vegetation.

The code for the simulation of a single sample uses the ProSail simulator configured with the satellite relative spectral responses, the acquisition parameters (angles) and the biophysical variables (leaf pigments, LAI, etc.):

ProSailType prosail;
prosail.SetRSR(satRSR);
prosail.SetParameters(prosailPars);
prosail.SetBVs(bio_vars);
auto result = prosail();

A simulation is computationally expensive and it would be difficult to parallelize the code. However, if many simulations are going to be run, it is worth using all the available CPU cores in the machine.

I show below how using C++11 standard support for threads allows to easily run many simulations in parallel.

Each simulation uses a set of variables given in an input file. We parse the sample file in order to get the input parameters for each sample and we construct a vector of simulations with the appropriate size to store the results.

otbAppLogINFO("Processing simulations ..." << std::endl);
auto bv_vec = parse_bv_sample_file(m_SampleFile);
auto sampleCount = bv_vec.size();
otbAppLogINFO("" << sampleCount << " samples read."<< std::endl);
std::vector<SimulationType> simus{sampleCount};

The simulation function is actually a lambda which will sequentially process a sequence of samples and store the results into the simus vector. We capture by reference the parameters which are the same for all simulations (the satellite relative spectral responses satRSR and the acquisition angles in prosailPars):

auto simulator = [&](std::vector<BVType>::const_iterator sample_first,
                     std::vector<BVType>::const_iterator sample_last,
                     std::vector<SimulationType>::iterator simu_first){
  ProSailType prosail;
  prosail.SetRSR(satRSR);
  prosail.SetParameters(prosailPars);
  while(sample_first != sample_last)
    {
    prosail.SetBVs(*sample_first);
    *simu_first = prosail();
    ++sample_first;
    ++simu_first;
    }
};    

We start by figuring out how to split the simulation into concurrent threads. How many cores are there?

auto num_threads = std::thread::hardware_concurrency();
otbAppLogINFO("" << num_threads << " CPUs available."<< std::endl);

So we define the size of the chunks we are going to run in parallel and we prepare a vector to store the threads:

auto block_size = sampleCount/num_threads;
if(num_threads>=sampleCount) block_size = sampleCount;
std::vector<std::thread> threads(num_threads);

Here, I choose to use as many threads as cores available, but this could be changed by a multiplicative factor if we know, for instance that disk I/O will introduce some idle time for each thread.

An now we can fill the vector with the threads that will process every block of simulations :

auto input_start = std::begin(bv_vec);
auto output_start = std::begin(simus);
for(size_t t=0; t<num_threads; ++t)
  {
  auto input_end = input_start;
  std::advance(input_end, block_size);
  threads[t] = std::thread(simulator,
                           input_start,
                           input_end,
                           output_start);
  input_start = input_end;
  std::advance(output_start, block_size);
  }

The std::thread takes the name of the function object to be called, followed by the arguments of this function, which in our case are the iterators to the beginning and the end of the block of samples to be processed and the iterator of the output vector where the results will be stored. We use std::advance to update the iterator positions.

After that, we have a vector of threads which have been started concurrently. In order to make sure that they have finished before trying to write the results to disk, we call join on each thread, which results in waiting for each thread to end:

std::for_each(threads.begin(),threads.end(),
              std::mem_fn(&std::thread::join));
otbAppLogINFO("" << sampleCount << " samples processed."<< std::endl);
for(const auto& s : simus)
  this->WriteSimulation(s);

This may no be the most efficient solution, nor the most general one. Using std::async and std::future would have allowed not to have to deal with block sizes, but in this solution we can easily decide the number of parallel threads that we want to use, which may be useful in a server shared with other users.

Should we always normalize image features before classification?

The quick answer is of course not. But if I am writing this post is because I have been myself giving a rule of thumb to students for a while without explaining the details.

When performing image classification (supervised or unsupervised) I always tell students to normalize the features. This normalization can be for instance just standardization (subtract the mean and divide by the standard deviation), although other approaches are possible.

The goal of this normalization is being able to compare apples and oranges: many classification algorithms are based on a distance in the feature space, and in order to give equal weight to all features, they have to have similar ranges. For instance, if we have 2 features like the NDVI, which ranges from -1 to 1, and the slope of a topographic surface which could be between \(0^\circ\) and \(90^\circ\), the Euclidean distance would give much more weight to the slope and the NDVI values would not have much influence in the classification.

This is why the TrainImagesClassifier and the ImageClassifier applications in the ORFEO Toolbox have the option to provide a statistics file with the mean and standard deviation of the features so they samples can be normalized.

This is needed for classifiers like SVM, unless custom kernels suited to particular sets of features are used.

Most OTB users have started using the Random Forest classifier since it was made available by the integration of the OpenCV Machine Learning module. And not long ago, someone told me that she did not notice any difference with and without feature normalization. Let's see why.

A Random Forest classifier uses a committee of decision trees. The decision trees split the samples at each node of the tree by thresholding one single feature. The learning of the classifier amounts basically at finding the value of the threshold at each node which is the optimum to split the samples.

The decision trees used by the Random Forest implementation in OpenCV use the Gini impurity in order to find the best split¹. At any node of the tree during the learning, the samples are sorted in increasing order of the feature. Then all possible threshold values are used to split the samples into 2 sets and the Gini impurity is computed for every possible split. The impurity index does not use the values of the features, but only the labels of the samples. The values of the features are only used in order to sort the samples.

Any pre-processing of the features which is monotonic will not change the value of the Gini impurity and therefore will have no effect on the training of a Random Forest classifier.

Data science in C?

Coursera is offering a Data Science specialization taught by professors from Johns Hopkins. I discovered it via one of their courses which is about reproducible research. I have been watching some of the video lectures and they are very interesting, since they combine data processing, programming and statistics.

The courses use the R language which is free software and is one the reference tools in the field of statistics, but it is also very much used for other kinds of data analysis.

While I was watching some of the lectures, I had some ideas to be implemented in some of my tools. Although I have used GSL and VXL¹ for linear algebra and optimization, I have never really needed statistics libraries in C or C++, so I ducked a bit and found the apophenia library², which builds on top of GSL, SQLite and Glib to provide very useful tools and data structures to do statistics in C.

Browsing a little bit more, I found that the author of apophenia has written a book "Modeling with data"³, which teaches you statistics like many books about R, but using C, SQLite and gnuplot.

This is the kind of technical book that I like most: good math and code, no fluff, just stuff!

The author sets the stage from the foreword. An example from page xii (the emphasis is mine):

" The politics of software

All of the software in this book is free software, meaning that it may be freely downloaded and distributed. This is because the book focuses on portability and replicability, and if you need to purchase a license every time you switch computers, then the code is not portable. If you redistribute a functioning program that you wrote based on the GSL or Apophenia, then you need to redistribute both the compiled final program and the source code you used to write the program. If you are publishing an academic work, you should be doing this anyway. If you are in a situation where you will distribute only the output of an analysis, there are no obligations at all. This book is also reliant on POSIX-compliant systems, because such systems were built from the ground up for writing and running replicable and portable projects. This does not exclude any current operating system (OS): current members of the Microsoft Windows family of OSes claim POSIX compliance, as do all OSes ending in X (Mac OS X, Linux, UNIX,…)."

Of course, the author knows the usual complaints about programming in C (or C++ for that matter) and spends many pages explaining his choice:

"I spent much of my life ignoring the fundamentals of computing and just hacking together projects using the package or language of the month: C++, Mathematica, Octave, Perl, Python, Java, Scheme, S-PLUS, Stata, R, and probably a few others that I’ve forgotten. Albee (1960, p 30)⁴ explains that “sometimes it’s necessary to go a long distance out of the way in order to come back a short distance correctly;” this is the distance I’ve gone to arrive at writing a book on data-oriented computing using a general and basic computing language. For the purpose of modeling with data, I have found C to be an easier and more pleasant language than the purpose-built alternatives—especially after I worked out that I could ignore much of the advice from books written in the 1980s and apply the techniques I learned from the scripting languages."

The author explains that C is a very simple language:

" Simplicity

C is a super-simple language. Its syntax has no special tricks for polymorphic operators, abstract classes, virtual inheritance, lexical scoping, lambda expressions, or other such arcana, meaning that you have less to learn. Those features are certainly helpful in their place, but without them C has already proven to be sufficient for writing some impressive programs, like the Mac and Linux operating systems and most of the stats packages listed above."

And he makes it really simple, since he actually teaches you C in one chapter of 50 pages (and 50 pages counting source code is not that much!). He does not teach you all C, though:

"As for the syntax of C, this chapter will cover only a subset. C has 32 keywords and this book will only use 18 of them."

At one point in the introduction I worried about the author bashing C++, which I like very much, but he actually gives a good explanation of the differences between C and C++ (emphasis and footnotes are mine):

"This is the appropriate time to answer a common intro-to-C question: What is the difference between C and C++? There is much confusion due to the almost-compatible syntax and similar name—when explaining the name C-double-plus⁵, the language’s author references the Newspeak language used in George Orwell’s 1984 (Orwell, 1949⁶; Stroustrup, 1986, p 4⁷). The key difference is that C++ adds a second scope paradigm on top of C’s file- and function-based scope: object-oriented scope. In this system, functions are bound to objects, where an object is effectively a struct holding several variables and functions. Variables that are private to the object are in scope only for functions bound to the object, while those that are public are in scope whenever the object itself is in scope. In C, think of one file as an object: all variables declared inside the file are private, and all those declared in a header file are public. Only those functions that have a declaration in the header file can be called outside of the file. But the real difference between C and C++ is in philosophy: C++ is intended to allow for the mixing of various styles of programming, of which object-oriented coding is one. C++ therefore includes a number of other features, such as yet another type of scope called namespaces, templates and other tools for representing more abstract structures, and a large standard library of templates. Thus, C represents a philosophy of keeping the language as simple and unchanging as possible, even if it means passing up on useful additions; C++ represents an all-inclusive philosophy, choosing additional features and conveniences over parsimony."

It is actually funny that I find myself using less and less class inheritance and leaning towards small functions (often templates) and when I use classes, it is usually to create functors. This is certainly due to the influence of the Algorithmic journeys of Alex Stepanov.

Is open and free global land cover mapping possible?

Short answer: yes.

Mid November took place in Toulouse "Le Capitole du libre", a conference on Free Software and Free Culture. The program this year was again full of interesting talks and workshops.

This year, I attended a workshop about contributing to Openstreetmap (OSM) using the JOSM software. The workshop was organised by Sébastien Dinot who is a massive contributor to OSM, and more importantly a very nice and passionate fellow.

I was very happy to learn to use JOSM and did 2 minor contributions right there.

During the workshop I learned that, over the past, OSM has been enriched using massive imports from open data sources, like for instance cadastral data bases from different countries or the Corine Land Cover data base. This has been possible thanks to the policies of many countries which have understood that the commons are important for the advancement of society. One example of this is the European INSPIRE initiative.

I was also interested to learn that what could be considered niche data, like agricultural land parcel data bases as for instance the French RPG have also been imported into OSM. Since I have been using the RPG at work for the last 4 years (see for example here or here), I was sympathetic with the difficulties of OSM contributors to efficiently exploit these data. I understood that the Corine Land Cover import was also difficult and the results were not fully satisfactory.

As a matter of fact, roads, buildings and other cartographic objects are easier to map than land cover, since they are discrete and sparse. They can be pointed, defined and characterised more easily than natural and semi-natural areas.

After that, I could not avoid making the link with what we do at work in terms of preparing the exploitation of upcoming satellite missions for automatic land cover map production.

One of our main interests is the use of Sentinel-2 images. It is the week end while I am writing this, so I will not use my free time to explain how land cover map production from multi-temporal satellite images work: I already did it in my day job.

What is therefore the link between what we do at work and OSM? The revolutionary thing from my point of view is the fact that Sentinel-2 data will be open and free, which means that the OSM project could use it to have a constantly up to date land cover layer.

Of course, Sentinel-2 data will come in huge volumes and a good amount of expertise will be needed to use them. However, several public agencies are paving the road in order to deliver data which is easy to use. For instance, the THEIA Land Data Centre will provide Sentinel-2 data which is ready to use for mapping. The data will be available with all the geometric and radiometric corrections of the best quality.

Actually, right now this is being done, for instance, for Landsat imagery. Of course, all these data is and will be available under open and free licences, which means that anyone can start right now learning how to use them.

However, going from images to land cover maps is not straightforward. Again, a good deal of expertise and efficient tools are needed in order to convert pixels into maps. This is what I have the chance to do at work: building tools to convert pixels into maps which are useful for real world applications.

Applying the same philosophy to tools as for data, the tools we produce are free and open. The core of all these tools is of course the Orfeo Toolbox, the Free Remote Sensing Image Processing Library from CNES. We have several times demonstrated that the tools are ready to efficiently exploit satellite imagery to produce maps. For instance, in this post here you even have the sequence of commands to generate land cover maps using satellite image time series.

This means that we have free data and free tools. Therefore, the complete pipeline is available for projects like OSM. OSM contributors could start right now getting familiar with these data and these tools.

Head over to CNES servers to get some Landsat data, install OTB, get familiar with the classification framework and see what could be done for OSM.

It is likely that some pieces may still be missing. For instance, the main approach for the map production is supervised classification. This means that we use machine learning algorithms to infer which land cover class is present at every given site using the images as input data. For these machine learning algorithms to work, we need training data, that is, we need to know before hand the correct land cover class in some places so the algorithm can be calibrated.

This training data is usually called ground truth and it is expensive and difficult to get. In a global mapping context, this can be a major drawback. However, there are interesting initiatives which could be leveraged to help here. For instance, Geo-Wiki comes to mind as a possible source of training data.

As always, talk is cheap, but it seems to me that exciting opportunities are available for open and free quality global mapping. This does not mean that the task is easy. It is not. There are many issues to be solved yet and some of them are at the research stage. But this should not stop motivated mappers and hackers to start learning to use the data and the tools.

xdg-open or "computer, do what I mean"

Command line nerds¹ say that graphical user interfaces are less efficient since the former use the keyboard and the latter tend to need the use of the mouse.

Graphical user interface fans say that command line based work-flows need the user to know more things about the system. Although I don't see why this would be an issue (what's the problem with getting to know better the system you use every day?), I can see one clear advantage to the point-and-click work-flow.

When the user wants to open a file, double-clicking on it just works. On the other hand, on the command line, the user needs to call the appropriate application and pass the file name as argument:

evince my_document.pdf

That means that the user has to know which is the application that is going to correctly deal with the given file. For example, this won't work:

evince my_movie.mp4

because evince is a document viewer and it only understands formats as pdf, postscript, dejavu or dvi, but it is not able to understand MPEG-4 videos.

So how can one replicate the point-and-click behavior on the command line? It would be nice to have something like:

open my_file

and that auto-magically the appropriate application was chosen.

It is of course possible, and it can be done in the same way as the graphical desktop environment does it: using xdg-utils.

Inside a desktop environment (e.g. GNOME, KDE, or Xfce), xdg-open simply passes the arguments to that desktop environment's file-opener application (gvfs-open, kde-open, or exo-open, respectively), which means that the associations are left up to the desktop environment. Therefore, on the command line, one can do:

xdg-open my_file.pdf

and the appropriate application will be used (evince on GNOME, okular on KDE, etc.).

When no desktop environment is detected (for example when one runs a standalone window manager, e.g. stumpwm), xdg-open will use its own configuration files.

Sometimes, the default associations between applications and file types may not suit the user. For instance, in my case, MPEG-4 videos were open with mplayer, but I prefer vlc. The xdg-mime tool is meant to help you change the default associations:

xdg-mime default vlc.desktop video/mp4

The vlc.desktop parameter is an xdg desktop file which describes the vlc applications. On my Debian GNU/Linux system, this files are located in /usr/share/applications. So I was able to look in there to see what applications xdg-open could use.

The parameter video/mp4 passed to xdg-mime is the type of file we want to associate the application with. In order to know what type of file we are dealing with, xdg-mime can query it for you:

xdg-mime query my_movie.mp4 

And the answer is:

video/mp4

Now, what happens if there is no desktop file for an application I want to use? For instance, remote sensing image visualization is not possible with the standard image viewers available on GNU/Linux. I personally use the OTB IceViewer, which is a lightweight application using the rendering engine developed for Monteverdi2.

In this case, I have to create a desktop file for the application. I my case, I have created the otbice.desktop file in /home/inglada/.local/share/applications/ with the following contents:

[Desktop Entry]
Type=Application
Name=OTB Ice Viewer
Exec=/home/inglada/OTB/builds/Ice/bin/otbiceviewer %U
Icon=otb-logo
Terminal=false
Categories=Graphics;2DGraphics;RasterGraphics;
MimeType=image/bmp;image/gif;image/x-portable-bitmap;image/x-portable-graymap;image/x-portable-pixmap;image/x-xbitmap;image/tiff;image/jpeg;image/png;image/x-xpixmap;image/jp2;image/jpeg2000;

The file has to be executable, so I have to do:

chmod +x /home/inglada/.local/share/applications/otbice.desktop

After that, I can associate the IceViewer with the image formats I want:

xdg-mime default otbice.desktop image/tiff

And it just works.

Why C++?

Recently, on the otb-users mailing list somebody asked why C++ was chosen to implement OTB.

The main reasons why I chose C++ for OTB were the following (by order of relevance):

1. Most of the existing libraries on which OTB was going to rely on were written in C++. Not only ITK, which is at the core of OTB, but also OSSIM and gdal.
2. C++ can be as efficient as C and has higher level abstraction mechanisms allowing for cleaner code and architecture (templates, classes, name spaces, etc.).
3. C++ was the language I was most familiar with.

C++ was first introduced in the late 70's and was first called "C with classes" since C++ was designed to provide Simula's facilities for program organization together with C's efficiency and flexibility for systems programming. The class concept (with derived classes and virtual functions) was borrowed from the Simula programming language.

C++ supports different programming styles:

• procedural and imperative styles due to its closeness to C,
• object orientation with classes, virtual functions and inheritance,
• generic programming through templates,
• functional programming through function objects.

Let's have a look at these different items.

OTB++11

auto x = sqrt(2.0);

auto x=0;
auto y=0.0;

typedef otb::ImageFileReader<ImageType> ReaderType;
ReaderType::Pointer reader = ReaderType::New();

typedef otb::ImageFileWriter<ImageType> WriterType;
WriterType::Pointer writer = WriterType::New();

auto reader = otb::ImageFileReader<ImageType>::New();
auto writer = otb::ImageFileWriter<ImageType>::New();

typedef itk::RescaleIntensityImageFilter<ImageType,
                                         OutputImageType> RescalerType;
RescalerType::Pointer rescaler = RescalerType::New();

auto rescaler = itk::RescaleIntensityImageFilter<ImageType,
                                                 OutputImageType>::New();

for (PolygonListType::Iterator it = polygonList->Begin();
     it != polygonList->End(); ++it)
  {
  DataNodeType::Pointer newPolygon = DataNodeType::New();
  newPolygon->SetPolygonExteriorRing(it.Get());
  tree->Add(newPolygon, multiPolygon);
  }

for (auto it = polygonList->Begin();
     it != polygonList->End(); ++it)
  {
  auto newPolygon = DataNodeType::New();
  newPolygon->SetPolygonExteriorRing(it.Get());
  tree->Add(newPolygon, multiPolygon);
  }

typedef otb::Image<unsigned int, 2> ImageType;
typedef otb::VectorImage<unsigned int, 2> VectorImageType;
typedef otb::Image<double, 2> DoubleImageType;
typedef otb::VectorImage<double, 2> DoubleVectorImageType;

using ImageType = otb::Image<unsigned int, 2>;
using VectorImageType = otb::VectorImage<unsigned int, 2>;
using DoubleImageType = otb::Image<double, 2>;
using DoubleVectorImageType = otb::VectorImage<double, 2>;

template<typename T> 
using Image2D = otb::Image<T, 2>;
auto im = Image2D<int>::New();  

ImageType::IndexType start = {0, 0};
decltype(start) start2 = {0, 0};

typedef otb::GenericMapProjection<otb::TransformDirection::INVERSE>
        InverseProjectionType;

InverseProjectionType::Pointer utmMapProjection =
                               InverseProjectionType::New();

typedef otb::OrthoRectificationFilter<ImageType, DoubleImageType,
                                      InverseProjectionType> 
             OrthoRectifFilterType;
OrthoRectifFilterType::Pointer orthoRectifPAN =
                               OrthoRectifFilterType::New();

auto utmMapProjection =
     otb::GenericMapProjection<otb::TransformDirection::INVERSE>::New();

auto orthoRectifPAN =
      otb::OrthoRectificationFilter<ImageType,
                                    DoubleImageType,
                                    decltype(utmMapProjection)::ObjectType>::New();

template<class TInput1, class TInput2, class TOutput>
class MyChangeDetector
{
public:
  MyChangeDetector() {}
  ~MyChangeDetector() {}
  inline TOutput operator ()(const TInput1& itA,
                             const TInput2& itB)
  {
    TOutput result = 0.0;

    for (unsigned long pos = 0; pos < itA.Size(); ++pos)
      {

      result += static_cast<TOutput>(itA.GetPixel(pos) - itB.GetPixel(pos));

      }
    return static_cast<TOutput>(result / itA.Size());
  }
};

template <class TInputImage1, class TInputImage2, class TOutputImage>
class ITK_EXPORT MyChangeDetectorImageFilter :
  public otb::BinaryFunctorNeighborhoodImageFilter<
      TInputImage1, TInputImage2, TOutputImage,
      MyChangeDetector<
          typename itk::ConstNeighborhoodIterator<TInputImage1>,
          typename itk::ConstNeighborhoodIterator<TInputImage2>,
          typename TOutputImage::PixelType> >
{
public:
  typedef MyChangeDetectorImageFilter Self;

  typedef typename otb::BinaryFunctorNeighborhoodImageFilter<
      TInputImage1, TInputImage2, TOutputImage,
      MyChangeDetector<
          typename itk::ConstNeighborhoodIterator<TInputImage1>,
          typename itk::ConstNeighborhoodIterator<TInputImage2>,
          typename TOutputImage::PixelType>
      >  Superclass;

  typedef itk::SmartPointer<Self>       Pointer;
  typedef itk::SmartPointer<const Self> ConstPointer;

  itkNewMacro(Self);

protected:
  MyChangeDetectorImageFilter() {}
  virtual ~MyChangeDetectorImageFilter() {}

private:
  MyChangeDetectorImageFilter(const Self &);
  void operator =(const Self&);

};

typedef MyChangeDetectorImageFilter<InputImageType1, InputImageType2,
    ChangeImageType>      FilterType;
FilterType::Pointer   filter = FilterType::New();

using CPT = ChangeImageType::PixelType;
using CNIType1 = itk::ConstNeighborhoodIterator<InputImageType1>;
using CNIType2 = itk::ConstNeighborhoodIterator<InputImageType2>;

std::function<CPT (CNIType1, CNIType2)>
  cdt = [](CNIType1 itA, CNIType2 itB){
  CPT result{0};
  for (auto pos = 0; pos < itA.Size(); ++pos)
    result += static_cast<CPT>(itA.GetPixel(pos) - itB.GetPixel(pos));
  return static_cast<CPT>(result / itA.Size());
};                                            

auto filter = otb::BinaryFunctorNeighborhoodImageFilter<InputImageType1,
                                                        InputImageType2,
                                                        ChangeImageType,
                                                        decltype(cdt)>::New();

Hack OTB on the Clojure REPL

    pacman -Suy

    useradd -m jordi
    passwd jordi

    pacman -S sudo
    vi /etc/sudoers

    sudo pacman -S mercurial make cmake gcc gdal openjdk6 \
                   mesa swig cvs bison links unzip rlwrap

    hg clone http://hg.orfeo-toolbox.org/OTB
    hg clone http://hg.orfeo-toolbox.org/OTB-Wrapping

    cvs -d :pserver:anonymous@pubilc.kitware.com:/cvsroot/CableSwig \
            co CableSwig

    mkdir builds
    cd builds/

    mkdir OTB
    cd OTB
    ccmake ../../OTB

    make

    cd ../
    mkdir CableSwig
    cd CableSwig/
    ccmake ../../CableSwig/
    make

    cd ../
    mkdir OTB-Wrapping
    cd OTB-Wrapping/
    ccmake ../../OTB-Wrapping/

    make
    cd

    mkdir Dev
    cd Dev/
    hg clone http://bitbucket.org/inglada/pwob

    links http://clojure.org

    mkdir src
    mv clojure-1.3.0.zip src/
    cd src/
    unzip clojure-1.3.0.zip

    cd clojure-1.3.0
    mkdir ~/.clojure
    mv clojure-1.3.0.jar ~/.clojure/

    ln -s /home/inglada/.clojure/clojure-1.3.0.jar /home/inglada/.clojure/clojure.jar

    cd ..
    rm -rf clojure-1.3.0*

    touch ~/.clj_completions

    export LD_LIBRARY_PATH="~/builds/OTB-Wrapping/lib/:~/builds/OTB/bin/"
    export ITK_AUTOLOAD_PATH="~/builds/OTB/bin/"

    cd ~/Dev/pwob/Clojure/src/Pwob
    ./otb-repl.sh

    Clojure 1.3.0
    user=>

    (import '(org.otb.application Registry))
    (def available-applications (Registry/GetAvailableApplications))

Lambda Omega Lambda

In my current research work, I am investigating the possibility of integrating domain expert knowledge together with image processing techniques (segmentation, classification) in order to provide accurate land cover maps from remote sensing image time series. When we look at this in terms of tools (software) there are a set of requirements which have to be met in order to be able to go from toy problems to operational processing chains. These operational chains are needed by scientists for carbon and water cycle studies, climate change analysis, etc. In the coming years, the image data needed for this kind of near-real-time mapping will be available through Earth observation missions like ESA's Sentinel program. Briefly, the constraints that the software tools for this task have are the following.

1. Interactive environment for exploratory analysis
2. Concurrent/parallel processing
3. Availability of state of the art, validated image processing algorithms
4. Symbolic, semantic information processing

This long post / short article describes some of the thoughts I have and some of the conclusions I have come to after a rather long analysis of the tools available.

    (make-otb-pipeline reader gradientFilter thresholdFilter writer)

    (def filter (RescaleIntensityImageFilter :itk (Image. : otb :Float 2)
                                                  (Image. : otb :UnsignedChar 2)))

    (def filter (itkRescaleIntensityImageFilterIF2IUC2.))

Geospatial analysis and Object-Based image analysis

I was searching in the web about the use of PostGIS data bases for object based image analysis (OBIA) and Google sent me to the OTB Wiki to a page that I wrote 6 months ago (first hit for "postgis object based image analysis").

It seems that this is a subject for which no much work is ongoing. Julien Michel will be presenting some results about how to put together OBIA and geospatial analysis (+ active learning and some other cool things) in a workshop in Paris next May.