lsrmBaatzSegmenter.txx

#ifndef __LSRM_BAATZ_SEGMENTER_TXX
#define __LSRM_BAATZ_SEGMENTER_TXX
#include <otbImageFileReader.h>
#include <itkImageRegionIterator.h>

namespace lsrm
{

	template<class TImage>
	BaatzSegmenter<TImage>::BaatzSegmenter() : Superclass()
	{
	}
	
	template<class TImage>
	void
	BaatzSegmenter<TImage>::InitFromImage()
	{
		typedef itk::ImageRegionIterator<TImage> ImageIterator;

		this->m_ImageWidth = this->m_InputImage->GetLargestPossibleRegion().GetSize()[0];
		this->m_ImageHeight =this->m_InputImage->GetLargestPossibleRegion().GetSize()[1];
		this->m_NumberOfComponentsPerPixel = this->m_InputImage->GetNumberOfComponentsPerPixel();

		std::size_t idx = 0;
		ImageIterator it(this->m_InputImage, this->m_InputImage->GetLargestPossibleRegion());
		for(it.GoToBegin(); !it.IsAtEnd(); ++it)
		{
			this->m_Graph.m_Nodes[idx]->m_Means.reserve(this->m_NumberOfComponentsPerPixel);
			this->m_Graph.m_Nodes[idx]->m_SquareMeans.reserve(this->m_NumberOfComponentsPerPixel);
			this->m_Graph.m_Nodes[idx]->m_SpectralSum.reserve(this->m_NumberOfComponentsPerPixel);
			this->m_Graph.m_Nodes[idx]->m_Std.assign(this->m_NumberOfComponentsPerPixel, 0.0f);

			for(std::size_t b = 0; b < this->m_NumberOfComponentsPerPixel; ++b)
			{
				this->m_Graph.m_Nodes[idx]->m_Means.push_back(it.Get()[b]);
				this->m_Graph.m_Nodes[idx]->m_SquareMeans.push_back((it.Get()[b])*(it.Get()[b]));
				this->m_Graph.m_Nodes[idx]->m_SpectralSum.push_back(it.Get()[b]);
			}	
			++idx;
		}
	}

	template<class TImage>
	float
	BaatzSegmenter<TImage>::ComputeMergingCost(NodePointerType n1, NodePointerType n2)
	{
		const unsigned int a1 = n1->m_Area, a2 = n2->m_Area, a_sum = a1 + a2;

		float spect_cost = 0.0f;
		float mean, square_mean, sum, std;

		for (unsigned int b = 0; b < this->m_NumberOfComponentsPerPixel; b++)
		{
			mean = ((a1 * n1->m_Means[b]) + (a2 * n2->m_Means[b])) / a_sum;
			square_mean = n1->m_SquareMeans[b] + n2->m_SquareMeans[b];
			sum = n1->m_SpectralSum[b] + n2->m_SpectralSum[b];
			std = std::sqrt((square_mean - 2*mean*sum + a_sum * mean* mean) / a_sum);
			spect_cost += (a_sum * std - a1 * n1->m_Std[b] - a2 * n2->m_Std[b]);
		}
		spect_cost *= this->m_Param.m_SpectralWeight;

		if(spect_cost < this->m_Threshold)
		{
			float shape_cost, smooth_f, compact_f;

 			// Compute the shape merging cost
			const float p1 = static_cast<float>(n1->m_Perimeter);
			const float p2 = static_cast<float>(n2->m_Perimeter);
			const unsigned int boundary = (GraphOperatorType::FindEdge(n1, n2))->m_Boundary;
			const float p3 = p1 + p2 - 2 * static_cast<float>(boundary);