lsrmBaatzSegmenter.txx 5.27 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
#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]);
		}
60
		spect_cost *= this->m_Param.m_SpectralWeight;
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81

		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);
			
			const lp::BoundingBox merged_bbox = lp::ContourOperations::MergeBoundingBoxes(n1->m_Bbox, n2->m_Bbox);
			const float bb1_perimeter = static_cast<float>(2*n1->m_Bbox.m_W + 2*n1->m_Bbox.m_H);
			const float bb2_perimeter = static_cast<float>(2*n2->m_Bbox.m_W + 2*n2->m_Bbox.m_H);
			const float mbb_perimeter = static_cast<float>(2 * merged_bbox.m_W + 2 * merged_bbox.m_H);

			smooth_f = a_sum*p3/mbb_perimeter - a1*p1/bb1_perimeter - a2*p2/bb2_perimeter;
			compact_f = a_sum*p3/std::sqrt(a_sum) - a1*p1/std::sqrt(a1) - a2*p2/std::sqrt(a2);

			shape_cost = this->m_Param.m_ShapeWeight * compact_f + (1-this->m_Param.m_ShapeWeight) * smooth_f;

82
			return (spect_cost + (1-this->m_Param.m_SpectralWeight)*shape_cost);
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
		}
		else
			return spect_cost;
	}

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

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

	template<class TImage>
	void
	BaatzSegmenter<TImage>::Update()
	{
		GraphOperatorType::InitNodes(this->m_InputImage, *this, FOUR);
		bool prev_merged = false;

		if(this->m_NumberOfIterations > 0)
		{
			prev_merged = GraphOperatorType::PerfomAllIterationsWithLMBFAndConstThreshold(*this);
																			
			this->m_Complete = !prev_merged;
			
			if(prev_merged && this->m_DoBFSegmentation)
			{
				prev_merged = GraphOperatorType::PerfomAllIterationsWithBFAndConstThreshold(*this);
				this->m_Complete = !prev_merged;  
			}
		}
		else
		{
			assert(this->m_DoBFSegmentation == true);
			prev_merged = GraphOperatorType::PerfomAllIterationsWithBFAndConstThreshold(*this);
			this->m_Complete = !prev_merged;
		}
	}
132
133

    template<class TImage>
134
    long long unsigned int
135
    BaatzSegmenter<TImage>::GetNodeMemory(NodePointerType &node)
136
137
138
139
140
141
    {

      unsigned int nBands = this->m_InputImage->GetNumberOfComponentsPerPixel();

      long long unsigned int memory = 0;
      memory += sizeof(NodePointerType);    // size of the node pointer
142
      memory += sizeof(NodeType);           // size of the node (actually, size of the base node)
143
      memory += 4 * nBands * sizeof(float); // size of the 4 attributes, multiplied by the nb. of bands
144
      memory += node->m_Edges.size() * sizeof(EdgeType);  // size of the edges
145
146
      return memory;
    }
147
148
149
150
151
152
153
154
155
156
} // end of namespace lsrm

#endif