Commit 8f495e21f0604882f29da07e6848f3fc09d65507

Authored by Pierre Lassalle
1 parent a1903882
Exists in master

accept modif Jordi

code/src/algo/BaatzAlgorithm.txx
... ... @@ -165,26 +165,26 @@ BaatzAlgorithm<TInputImage>::Segmentation()
165 165 RegionPointerType bestneigh = nullptr;
166 166 bool prev_merged = true;
167 167  
168   - auto cw = m_Parameters.m_ColorWeight;
169   - auto sc = m_Parameters.m_Scale;
170   - m_ComputeCostFunction = [&](RegionPointerType r1, RegionPointerType r2)->double
171   - {
172   - double spectral_h;
173   - double spatial_h;
174   - double cost = 0;
175   -
176   - spectral_h = ColorComponentCostFusion(r1, r2);
177   - cost += cw * spectral_h;
178   -
179   - if(cost < sc)
180   - {
181   - spatial_h = CompactnessComponentCostFusion(r1, r2);
182   - cost += (1-cw)*spatial_h;
183   - return cost;
184   - }
185   - else
186   - return cost;
187   - };
  168 + auto cw = m_Parameters.m_ColorWeight;
  169 + auto sc = m_Parameters.m_Scale;
  170 + m_ComputeCostFunction = [&](RegionPointerType r1, RegionPointerType r2)->double
  171 + {
  172 + double spectral_h;
  173 + double spatial_h;
  174 + double cost = 0;
  175 +
  176 + spectral_h = ColorComponentCostFusion(r1, r2);
  177 + cost += cw * spectral_h;
  178 +
  179 + if(cost < sc)
  180 + {
  181 + spatial_h = CompactnessComponentCostFusion(r1, r2);
  182 + cost += (1-cw)*spatial_h;
  183 + return cost;
  184 + }
  185 + else
  186 + return cost;
  187 + };
188 188  
189 189 while(curr_step < max_iter && prev_merged)
190 190 {
... ...
code/src/algo/BaatzAlgorithm.txx~ 0 → 100644
... ... @@ -0,0 +1,338 @@
  1 +/*=========================================================================
  2 +
  3 + Program: Large Scale Segmentation (LSS)
  4 + Language: C++
  5 + author: Lassalle Pierre
  6 +
  7 +
  8 +
  9 + Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
  10 +
  11 +
  12 + This software is distributed WITHOUT ANY WARRANTY; without even
  13 + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  14 + PURPOSE. See the above copyright notices for more information.
  15 +
  16 +=========================================================================*/
  17 +#ifndef __BaatzAlgorithm_txx
  18 +#define __BaatzAlgorithm_txx
  19 +
  20 +template<class TInputImage>
  21 +void
  22 +BaatzAlgorithm<TInputImage>::SetDimensionForEncoder(unsigned int cols, unsigned int rows)
  23 +{
  24 + m_RMHandler.m_Encoder.SetCols(cols);
  25 + m_RMHandler.m_Encoder.SetRows(rows);
  26 +}
  27 +
  28 +template<class TInputImage>
  29 +void
  30 +BaatzAlgorithm<TInputImage>::InitFromTiffImage(const std::string& input_file)
  31 +{
  32 + auto reader = ReaderType::New();
  33 + reader->SetFileName(input_file);
  34 + reader->Update();
  35 +
  36 + auto regionToProcess = reader->GetOutput()->GetLargestPossibleRegion();
  37 + unsigned int cols = regionToProcess.GetSize()[0];
  38 + unsigned int rows = regionToProcess.GetSize()[1];
  39 + unsigned int bands = reader->GetOutput()->GetNumberOfComponentsPerPixel();
  40 +
  41 + m_RMHandler.m_Encoder.SetRows(rows); // Penser à n'utiliser qu'une seule méthode (SetRegionToProcess)
  42 + m_RMHandler.m_Encoder.SetCols(cols);
  43 +
  44 + m_RegionList.reserve(cols*rows);
  45 +
  46 + IteratorType it(reader->GetOutput(), regionToProcess);
  47 +
  48 + // Create the vertices
  49 + long unsigned int idx = 0;
  50 + for(it.GoToBegin(); !it.IsAtEnd(); ++it)
  51 + {
  52 + auto r = std::make_shared<BaatzRegion>();
  53 + m_RMHandler.InitGenericRegion(r, idx);
  54 + r->m_Area = 1;
  55 + r->m_Perimeter = 4;
  56 + r->m_Avg_Color.reserve(bands);
  57 + r->m_Avg_Color_Square.reserve(bands);
  58 + r->m_Color_Sum.reserve(bands);
  59 + r->m_Std_Color.reserve(bands);
  60 +
  61 + for(unsigned int b = 0; b < bands; b++)
  62 + {
  63 + r->m_Avg_Color.push_back(it.Get()[b]);
  64 + r->m_Avg_Color_Square.push_back(pow(it.Get()[b], 2));
  65 + r->m_Color_Sum.push_back(it.Get()[b]);
  66 + r->m_Std_Color.push_back(0);
  67 + }
  68 +
  69 + m_RegionList.push_back(r);
  70 + ++idx;
  71 + }
  72 + m_RMHandler.InitNeighborhood(m_RegionList);
  73 +}
  74 +
  75 +template<class TInputImage>
  76 +double
  77 +BaatzAlgorithm<TInputImage>::ColorComponentCostFusion(RegionPointerType r1, RegionPointerType r2)
  78 +{
  79 + unsigned int bands = r1->m_Avg_Color.size();
  80 + float mean[bands], colorSum[bands];
  81 + float squarePixels[bands];
  82 + float stddev[bands];
  83 + float stddevNew[bands];
  84 + double color_f[bands];
  85 + double color_h;
  86 + float a_current, a_neighbor, a_sum;
  87 + a_current = r1->m_Area;
  88 + a_neighbor = r2->m_Area;
  89 + a_sum = a_current+a_neighbor;
  90 +
  91 + for (unsigned int b = 0; b < bands; b++)
  92 + {
  93 + mean[b] = ((r1->m_Avg_Color[b]*a_current)+(r2->m_Avg_Color[b]*a_neighbor))/a_sum;
  94 + squarePixels[b] = (r1->m_Avg_Color_Square[b])+(r2->m_Avg_Color_Square[b]);
  95 + colorSum[b] = r1->m_Color_Sum[b] + r2->m_Color_Sum[b];
  96 + stddev[b] = 0;
  97 + stddevNew[b] = 0;
  98 + }
  99 +
  100 + for(unsigned int b = 0; b < bands; b++)
  101 + {
  102 + stddevNew[b] = squarePixels[b] - 2*mean[b]*colorSum[b] + a_sum*mean[b]*mean[b];
  103 + }
  104 +
  105 + /* calculates color factor per band and total */
  106 + color_h = 0;
  107 + for (unsigned int b = 0; b < bands; b++)
  108 + {
  109 + stddev[b] = sqrt(stddevNew[b]/a_sum);
  110 + color_f[b] = (a_current*r1->m_Std_Color[b]) + (a_neighbor*r2->m_Std_Color[b]);
  111 + color_f[b] = (a_sum*stddev[b])- color_f[b];
  112 + color_h += color_f[b];
  113 + }
  114 + return color_h;
  115 +}
  116 +
  117 +template<class TInputImage>
  118 +double
  119 +BaatzAlgorithm<TInputImage>::CompactnessComponentCostFusion(RegionPointerType r1, RegionPointerType r2)
  120 +{
  121 + double spatial_h, smooth_f, compact_f;
  122 + float area[3], perimeter[3], b_box_len[3]; /* 0-current segment; 1-neighbor segment; 2-merged (new) segment */
  123 +
  124 + /* area */
  125 + area[0] = r1->m_Area;
  126 + area[1] = r2->m_Area;
  127 + area[2] = area[0]+area[1];
  128 +
  129 + /* perimeter */
  130 + perimeter[0] = r1->m_Perimeter;
  131 + perimeter[1] = r2->m_Perimeter;
  132 + perimeter[2] = r1->m_Perimeter + r2->m_Perimeter - 2*m_RMHandler.GetConnections(r1, r2);
  133 +
  134 + /* bounding box lenght */
  135 + auto mbbox = m_RMHandler.GetResultingBbox(r1, r2);
  136 + b_box_len[0] = (r1->m_Bbox.GetSize(0))*2 + (r1->m_Bbox.GetSize(1))*2;
  137 + b_box_len[1] = (r2->m_Bbox.GetSize(0))*2 + (r2->m_Bbox.GetSize(1))*2;
  138 + b_box_len[2] = (mbbox.GetSize(0))*2 + (mbbox.GetSize(1))*2;
  139 +
  140 + /* smoothness factor */
  141 + smooth_f = (area[2]*perimeter[2]/b_box_len[2] -
  142 + (area[1]*perimeter[1]/b_box_len[1] + area[0]*perimeter[0]/b_box_len[0]));
  143 +
  144 + /* compactness factor */
  145 + compact_f = (area[2]*perimeter[2]/sqrt(area[2]) -
  146 + (area[1]*perimeter[1]/sqrt(area[1]) + area[0]*perimeter[0]/sqrt(area[0])));
  147 +
  148 + /* spatial heterogeneity */
  149 + spatial_h = m_Parameters.m_CompactnessWeight*compact_f + (1-m_Parameters.m_CompactnessWeight)*smooth_f;
  150 +
  151 + return spatial_h;
  152 +}
  153 +
  154 +template<class TInputImage>
  155 +void
  156 +BaatzAlgorithm<TInputImage>::Segmentation()
  157 +{
  158 + m_PreviousMerge = false;
  159 + unsigned int max_iter = 500;
  160 + if(m_NumberOfIterations > 0)
  161 + max_iter = m_NumberOfIterations;
  162 +
  163 + unsigned int curr_step = 0;
  164 + double cost;
  165 + RegionPointerType bestneigh = nullptr;
  166 + bool prev_merged = true;
  167 +
  168 + auto cw = m_Parameters.m_ColorWeight;
  169 + auto sc = m_Parameters.m_Scale;
  170 + m_ComputeCostFunction = [&](RegionPointerType r1, RegionPointerType r2)->double
  171 + {
  172 + double spectral_h;
  173 + double spatial_h;
  174 + double cost = 0;
  175 +
  176 + spectral_h = ColorComponentCostFusion(r1, r2);
  177 + cost += cw * spectral_h;
  178 +
  179 + if(cost < sc)
  180 + {
  181 + spatial_h = CompactnessComponentCostFusion(r1, r2);
  182 + cost += (1-cw)*spatial_h;
  183 + return cost;
  184 + }
  185 + else
  186 + return cost;
  187 + };
  188 +
  189 + while(curr_step < max_iter && prev_merged)
  190 + {
  191 + prev_merged = false;
  192 + std::cout << curr_step << std::endl;
  193 + m_RMHandler.UpdateMergingCost(m_RegionList, m_ComputeCostFunction);
  194 +
  195 + for(auto& r: m_RegionList)
  196 + {
  197 + if(m_RMHandler.IsLMBF(r, m_Parameters.m_Scale))
  198 + {
  199 + auto res = m_RMHandler.GetRegionsToMerge(r);
  200 + //std::cout << "Merge " << res.first->m_Id << " with " << res.second->m_Id << std::endl;
  201 + // User contribution
  202 + UpdateAttribute(res.first, res.second);
  203 + m_RMHandler.Update(res);
  204 + prev_merged = true;
  205 + }
  206 + }
  207 + m_RMHandler.RemoveExpiredVertices(m_RegionList);
  208 + curr_step++;
  209 + m_PreviousMerge = prev_merged;
  210 + }
  211 +}
  212 +
  213 +template<class TInputImage>
  214 +void
  215 +BaatzAlgorithm<TInputImage>::UpdateAttribute(RegionPointerType r1, RegionPointerType r2)
  216 +{
  217 + // Update the spectral attributes
  218 +
  219 + unsigned int bands = r1->m_Avg_Color.size();
  220 + float mean[bands], colorSum[bands];
  221 + float squarePixels[bands];
  222 + float a_current, a_neighbor, a_sum;
  223 + a_current = r1->m_Area;
  224 + a_neighbor = r2->m_Area;
  225 + a_sum = a_current+a_neighbor;
  226 +
  227 + for (unsigned int b = 0; b < bands; b++)
  228 + {
  229 + mean[b] = ((r1->m_Avg_Color[b]*a_current)+(r2->m_Avg_Color[b]*a_neighbor))/a_sum;
  230 + squarePixels[b] = (r1->m_Avg_Color_Square[b])+(r2->m_Avg_Color_Square[b]);
  231 + colorSum[b] = r1->m_Color_Sum[b] + r2->m_Color_Sum[b];
  232 + }
  233 +
  234 + for(unsigned int b = 0; b < bands; b++)
  235 + {
  236 + r1->m_Avg_Color[b] = mean[b];
  237 + r1->m_Avg_Color_Square[b] = squarePixels[b];
  238 + r1->m_Color_Sum[b] = colorSum[b];
  239 + r1->m_Std_Color[b] = sqrt((squarePixels[b] - 2*mean[b]*colorSum[b] + a_sum*mean[b]*mean[b])/a_sum);
  240 + }
  241 +
  242 + // Update spatial attributes
  243 + r1->m_Area += r2->m_Area;
  244 + r1->m_Perimeter += r2->m_Perimeter - 2*m_RMHandler.GetConnections(r1,r2);
  245 +}
  246 +
  247 +template<class TInputImage>
  248 +void
  249 +BaatzAlgorithm<TInputImage>::WriteLabelImage(const std::string& ofname)
  250 +{
  251 + typename InputImageType::IndexType index;
  252 + typename InputImageType::SizeType size;
  253 + typename InputImageType::RegionType region;
  254 + auto label_image = LabelImageType::New();
  255 +
  256 + index[0] = 0;
  257 + index[1] = 0;
  258 + size[0] = m_RMHandler.m_Encoder.GetCols();
  259 + size[1] = m_RMHandler.m_Encoder.GetRows();
  260 + region.SetIndex(index);
  261 + region.SetSize(size);
  262 +
  263 + label_image->SetRegions(region);
  264 + label_image->Allocate();
  265 +
  266 + unsigned int label = 1;
  267 + for(auto& v : m_RegionList)
  268 + {
  269 + auto pixels = m_RMHandler.m_Encoder.GenerateAllPixels(v->m_Id, v->m_Contour, v->m_Bbox);
  270 + for(auto& pix: pixels)
  271 + {
  272 + index[0] = pix % size[0];
  273 + index[1] = pix / size[0];
  274 + auto label_pixel = label_image->GetPixel(index);
  275 + label_image->SetPixel(index, label);
  276 + }
  277 + label++;
  278 + }
  279 +
  280 + auto label_writer = LabelWriterType::New();
  281 + label_writer->SetFileName(ofname);
  282 + label_writer->SetInput(label_image);
  283 + label_writer->Update();
  284 +}
  285 +
  286 +template<class TInputImage>
  287 +void
  288 +BaatzAlgorithm<TInputImage>::WriteResultingImage(const std::string& ofname)
  289 +{
  290 + auto out_image = InputImageType::New();
  291 + typename InputImageType::IndexType index;
  292 + typename InputImageType::SizeType size;
  293 + typename InputImageType::RegionType region;
  294 +
  295 + index[0] = 0;
  296 + index[1] = 0;
  297 + size[0] = m_RMHandler.m_Encoder.GetCols();
  298 + size[1] = m_RMHandler.m_Encoder.GetRows();
  299 + region.SetIndex(index);
  300 + region.SetSize(size);
  301 + out_image->SetRegions(region);
  302 + out_image->SetNumberOfComponentsPerPixel(3);
  303 + out_image->Allocate();
  304 +
  305 + for(unsigned int r = 0; r< size[1]; r++)
  306 + {
  307 + for(unsigned int c = 0; c<size[0]; c++)
  308 + {
  309 + index[0] = c;
  310 + index[1] = r;
  311 + auto pixel_value = out_image->GetPixel(index);
  312 + for(int b=0; b<3; b++)
  313 + pixel_value[b] = 0;
  314 + out_image->SetPixel(index, pixel_value);
  315 + }
  316 + }
  317 + unsigned int label = 1;
  318 + for(auto& v : m_RegionList)
  319 + {
  320 + auto pixels = m_RMHandler.m_Encoder.GeneratePixels(v->m_Id, v->m_Contour);
  321 + for(auto& pix: pixels)
  322 + {
  323 + index[0] = pix % size[0];
  324 + index[1] = pix / size[0];
  325 + auto pixel_value = out_image->GetPixel(index);
  326 + for(int b=0; b<3; b++)
  327 + pixel_value[b] = v->m_Avg_Color[b];
  328 + out_image->SetPixel(index, pixel_value);
  329 + }
  330 + }
  331 +
  332 + auto writer = WriterType::New();
  333 + writer->SetFileName(ofname);
  334 + writer->SetInput(out_image);
  335 + writer->Update();
  336 +}
  337 +
  338 +#endif
... ...