All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Skeleton.cpp
Go to the documentation of this file.
1 /* Copyright (C) 2001-2009 National Institute For Space Research (INPE) - Brazil.
2 
3  This file is part of the TerraLib - a Framework for building GIS enabled applications.
4 
5  TerraLib is free software: you can redistribute it and/or modify
6  it under the terms of the GNU Lesser General Public License as published by
7  the Free Software Foundation, either version 3 of the License,
8  or (at your option) any later version.
9 
10  TerraLib is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU Lesser General Public License for more details.
14 
15  You should have received a copy of the GNU Lesser General Public License
16  along with TerraLib. See COPYING. If not, write to
17  TerraLib Team at <terralib-team@terralib.org>.
18  */
19 
20 /*!
21  \file terralib/rp/Skeleton.cpp
22  \brief Creation of skeleton imagems.
23 */
24 
25 #include "Skeleton.h"
26 #include "Macros.h"
27 #include "Functions.h"
28 
29 #include "../raster/Band.h"
30 #include "../raster/Utils.h"
31 #include "../common/PlatformUtils.h"
32 
33 #include <boost/lexical_cast.hpp>
34 #include <boost/scoped_array.hpp>
35 
36 #include <memory>
37 #include <cmath>
38 #include <cstring>
39 
40 #define DIFFUSENEIGHBOR( neighborRow, neighborCol ) \
41  if( bufferMag[ neighborRow ][ neighborCol ] >= centerMag ) \
42  { \
43  ++validNeigborsNumber; \
44  newCenterX += bufferX[ neighborRow ][ neighborCol ]; \
45  newCenterY += bufferY[ neighborRow ][ neighborCol ]; \
46  }
47 
48 #define SKELSTRENGTHNEIGHBOR( neighborRow, neighborCol ) \
49  neiX = inputX[ neighborRow ][ neighborCol ]; \
50  neiY = inputY[ neighborRow ][ neighborCol ]; \
51  diffX = neiX - centerX; \
52  diffY = neiY - centerY; \
53  diffMag = std::sqrt( ( diffX * diffX ) + ( diffY * diffY ) ); \
54  if( diffMag != 0.0 ) \
55  { \
56  strength += ( ( ( diffX * neiX ) + ( diffY * neiY ) ) / diffMag ); \
57  }
58 
59 namespace te
60 {
61  namespace rp
62  {
63 
65  {
66  reset();
67  }
68 
70  {
71  reset();
72  operator=( other );
73  }
74 
76  {
77  reset();
78  }
79 
80  void Skeleton::InputParameters::reset() throw( te::rp::Exception )
81  {
82  m_inputRasterPtr = 0;
83  m_inputRasterBand = 0;
84  m_inputMaskRasterPtr = 0;
85  m_diffusionThreshold = 0.5;
86  m_diffusionRegularization = 0.5;
87  m_diffusionMaxIterations = 0;
88  m_enableMultiThread = true;
89  m_skeletonThreshold = 0.75;
90  m_enableProgress = false;
91  }
92 
94  const Skeleton::InputParameters& params )
95  {
96  reset();
97 
98  m_inputRasterPtr = params.m_inputRasterPtr;
99  m_inputRasterBand = params.m_inputRasterBand;
100  m_inputMaskRasterPtr = params.m_inputMaskRasterPtr;
101  m_diffusionThreshold = params.m_diffusionThreshold;
102  m_diffusionRegularization = params.m_diffusionRegularization;
103  m_diffusionMaxIterations = params.m_diffusionMaxIterations;
104  m_enableMultiThread = params.m_enableMultiThread;
105  m_skeletonThreshold = params.m_skeletonThreshold;
106  m_enableProgress = params.m_enableProgress;
107 
108  return *this;
109  }
110 
112  {
113  return new InputParameters( *this );
114  }
115 
117  {
118  reset();
119  }
120 
122  {
123  reset();
124  operator=( other );
125  }
126 
128  {
129  reset();
130  }
131 
132  void Skeleton::OutputParameters::reset() throw( te::rp::Exception )
133  {
134  m_rType.clear();
135  m_rInfo.clear();
136  m_outputRasterPtr.reset();
137  }
138 
140  const Skeleton::OutputParameters& params )
141  {
142  reset();
143 
144  m_rType = params.m_rType;
145  m_rInfo = params.m_rInfo;
146 
147  return *this;
148  }
149 
151  {
152  return new OutputParameters( *this );
153  }
154 
156  {
157  reset();
158  }
159 
161  {
162  }
163 
165  throw( te::rp::Exception )
166  {
167  if( ! m_isInitialized ) return false;
168 
169  Skeleton::OutputParameters* outParamsPtr = dynamic_cast<
170  Skeleton::OutputParameters* >( &outputParams );
171  TERP_TRUE_OR_THROW( outParamsPtr, "Invalid paramters" );
172 
173  // progress
174 
175  std::auto_ptr< te::common::TaskProgress > progressPtr;
177  {
178  progressPtr.reset( new te::common::TaskProgress );
179 
180  progressPtr->setTotalSteps( 7 );
181 
182  progressPtr->setMessage( "Generating the skeleton" );
183  }
184 
185  // loading raster data
186 
187  std::auto_ptr< te::rp::Matrix< double > > rasterDataPtr(
188  new te::rp::Matrix< double > ); ;
189  rasterDataPtr->reset( te::rp::Matrix< double >::AutoMemPol );
190  TERP_TRUE_OR_RETURN_FALSE( loadData( *rasterDataPtr ),
191  "Gradient maps build error" );
192 // CreateRasterFileFromMatrix( *rasterDataPtr, true, "rasterData.tif" );
193 
195  {
196  progressPtr->pulse();
197  if( ! progressPtr->isActive() ) return false;
198  }
199 
200  // Applying the smooth filter
201 
202  {
203  std::auto_ptr< te::rp::Matrix< double > > smoothedRasterDataPtr(
204  new te::rp::Matrix< double > ); ;
205  smoothedRasterDataPtr->reset( te::rp::Matrix< double >::AutoMemPol );
206  TERP_TRUE_OR_RETURN_FALSE( applyMeanSmooth( *rasterDataPtr, *smoothedRasterDataPtr ),
207  "Raster data smoothing error" );
208 // CreateRasterFileFromMatrix( *smoothedRasterDataPtr, true, "smoothedRasterData.tif" );
209  rasterDataPtr.reset( smoothedRasterDataPtr.release() );
210  }
211 
213  {
214  progressPtr->pulse();
215  if( ! progressPtr->isActive() ) return false;
216  }
217 
218  // calculating the edge strength map
219 
220  te::rp::Matrix< double > edgeStrengthMap;
221  edgeStrengthMap.reset( te::rp::Matrix< double >::AutoMemPol );
223  edgeStrengthMap ), "Edge strength map build error" );
224 // CreateRasterFileFromMatrix( edgeStrengthMap, true, "edgeStrengthMap.tif" );
225 
226  rasterDataPtr.reset();
227 
229  {
230  progressPtr->pulse();
231  if( ! progressPtr->isActive() ) return false;
232  }
233 
234  // Generating the initial vector field
235 
236  te::rp::Matrix< double > vecXMap;
238  te::rp::Matrix< double > vecYMap;
240 // TERP_TRUE_OR_RETURN_FALSE( getEdgeVecField( edgeStrengthMap, true,
241 // vecXMap, vecYMap ),
242  TERP_TRUE_OR_RETURN_FALSE( getGradientMaps( edgeStrengthMap, true, vecXMap, vecYMap ),
243  "Vector maps build error" );
244 // {
245 // CreateRasterFileFromMatrix( vecXMap, true, "vecXMap.tif" );
246 // CreateRasterFileFromMatrix( vecYMap, true, "vecYMap.tif" );
247 // createTifFromVecField( vecXMap, vecYMap, &edgeStrengthMap, 3, "vecMap" );
248 // te::rp::Matrix< double > vecMagMap;
249 // getMagnitude( vecXMap, vecYMap, vecMagMap );
250 // CreateRasterFileFromMatrix( vecMagMap, true, "vecMagMap.tif" );
251 // }
252 
254  {
255  progressPtr->pulse();
256  if( ! progressPtr->isActive() ) return false;
257  }
258 
259  // Applying a vector diffusion
260 
261  te::rp::Matrix< double > diffusedVecXMap;
262  diffusedVecXMap.reset( te::rp::Matrix< double >::AutoMemPol );
263  te::rp::Matrix< double > diffusedVecYMap;
264  diffusedVecYMap.reset( te::rp::Matrix< double >::AutoMemPol );
265  TERP_TRUE_OR_RETURN_FALSE( applyVecDiffusion( vecXMap, vecYMap, 0,
266  progressPtr.get(), diffusedVecXMap, diffusedVecYMap ),
267  "Vector maps build error" );
268 // {
269 // CreateRasterFileFromMatrix( diffusedVecXMap, true, "diffusedVecXMap.tif" );
270 // CreateRasterFileFromMatrix( diffusedVecYMap, true, "diffusedVecYMap.tif" );
271 // createTifFromVecField( diffusedVecXMap, diffusedVecYMap, rasterDataPtr.get(),
272 // 4, "diffusedVecMap" );
273 // te::rp::Matrix< double > diffusedVecXMagMap;
274 // getMagnitude( diffusedVecXMap, diffusedVecYMap, diffusedVecXMagMap );
275 // CreateRasterFileFromMatrix( diffusedVecXMagMap, true, "diffusedVecXMagMap.tif" );
276 // }
277 
278  vecXMap.clear();
279  vecYMap.clear();
280 
282  {
283  progressPtr->pulse();
284  if( ! progressPtr->isActive() ) return false;
285  }
286 
287  // Calculating the skeleton strength map
288 
289  te::rp::Matrix< double > skelSMap;
292  diffusedVecYMap, edgeStrengthMap, skelSMap ),
293  "Skeleton strength map build error" );
294 // CreateRasterFileFromMatrix( skelSMap, true, "skelSMap.tif" );
295 
296  diffusedVecXMap.reset();
297  diffusedVecYMap.reset();
298 
300  {
301  progressPtr->pulse();
302  if( ! progressPtr->isActive() ) return false;
303  }
304 
305  // creating the output raster
306 
307  {
308  std::vector< te::rst::BandProperty* > bandsProperties;
309  bandsProperties.push_back( new te::rst::BandProperty(
312  bandsProperties[ 0 ]->m_noDataValue = 0;
313  bandsProperties[ 0 ]->m_colorInterp = te::rst::GrayIdxCInt;
314  bandsProperties[ 0 ]->m_type = te::dt::UCHAR_TYPE;
315 
316  outParamsPtr->m_outputRasterPtr.reset(
318  outParamsPtr->m_rType,
320  bandsProperties,
321  outParamsPtr->m_rInfo,
322  0,
323  0 ) );
324  TERP_TRUE_OR_RETURN_FALSE( outParamsPtr->m_outputRasterPtr.get(),
325  "Output raster creation error" );
326 
327  te::rst::Band& outBand = *( outParamsPtr->m_outputRasterPtr->getBand( 0 ) );
328  const unsigned int nRows = outParamsPtr->m_outputRasterPtr->getNumberOfRows();
329  const unsigned int nCols = outParamsPtr->m_outputRasterPtr->getNumberOfColumns();
330  const unsigned int rowsBound = nRows - 4;
331  const unsigned int colsBound = nCols - 4;
332  unsigned int row = 0;
333  unsigned int col = 0;
334 
335  for( row = 0 ; row < nRows ; ++row )
336  {
337  for( col = 0 ; col < nCols ; ++col )
338  {
339  if(
340  ( skelSMap[ row ][ col ] > m_inputParameters.m_skeletonThreshold )
341  && ( row > 0 ) && ( row < rowsBound )
342  && ( col > 0 ) && ( col < colsBound )
343  )
344  {
345  outBand.setValue( col, row, 255 );
346  }
347  else
348  {
349  outBand.setValue( col, row, 0 );
350  }
351  }
352  }
353  }
354 
356  {
357  progressPtr->pulse();
358  if( ! progressPtr->isActive() ) return false;
359  }
360 
361  return true;
362  }
363 
364  void Skeleton::reset() throw( te::rp::Exception )
365  {
367  m_isInitialized = false;
368  }
369 
371  throw( te::rp::Exception )
372  {
373  reset();
374 
375  Skeleton::InputParameters const* inputParamsPtr = dynamic_cast<
376  Skeleton::InputParameters const* >( &inputParams );
377  TERP_TRUE_OR_THROW( inputParamsPtr, "Invalid paramters pointer" );
378 
379  m_inputParameters = *inputParamsPtr;
380 
381  // Checking the input raster
382 
384  "Invalid m_inputRasterPtr" )
385 
388  "Invalid raster" );
389 
393  "Invalid raster band" );
394 
396  {
399  "Invalid mask raster" );
400 
403  "Invalid mask raster band" );
404 
408  "Invalid mask raster number of rows" );
409 
413  "Invalid mask raster number of columns" );
414  }
415 
419  "Invalid diffusion threshold" );
420 
424  "Invalid diffusion regularization" );
425 
429  "Invalid diffusion regularization" );
430 
431  m_isInitialized = true;
432 
433  return true;
434  }
435 
437  {
438  return m_isInitialized;
439  }
440 
442  {
443  const unsigned int nRows = m_inputParameters.m_inputRasterPtr->getNumberOfRows();
444  const unsigned int nCols = m_inputParameters.m_inputRasterPtr->getNumberOfColumns();
445 
446  if( ! rasterData.reset( nRows,nCols ) )
447  return false;
448 
449  unsigned int row = 0;
450  unsigned int col = 0;
453  te::rst::Band const * const maskBandPtr = m_inputParameters.m_inputMaskRasterPtr ?
455  double maskValue = 0;
456 
457  for( row = 0 ; row < nRows ; ++row )
458  {
459  for( col = 0 ; col < nCols ; ++col )
460  {
461  if( maskBandPtr )
462  {
463  maskBandPtr->getValue( col, row, maskValue );
464  if( maskValue == 0.0 )
465  {
466  rasterData[ row ][ col ] = 0;
467  continue;
468  }
469  }
470 
471  band.getValue( col, row , rasterData[ row ][ col ] );
472  }
473  }
474 
475  return true;
476  }
477 
479  const bool unitVectors,
480  te::rp::Matrix< double >& gradXMap,
481  te::rp::Matrix< double >& gradYMap ) const
482  {
483  const unsigned int nRows = inputData.getLinesNumber();
484  TERP_TRUE_OR_RETURN_FALSE( nRows > 2, "Invalid rows number" );
485 
486  const unsigned int nCols = inputData.getColumnsNumber();
487  TERP_TRUE_OR_RETURN_FALSE( nCols > 2, "Invalid columns number" );
488 
489  const unsigned int lastRowIdx = nRows - 1;
490  const unsigned int lastColIdx = nCols - 1;
491 
492  if( ! gradXMap.reset( nRows,nCols ) )
493  return false;
494  if( ! gradYMap.reset( nRows,nCols ) )
495  return false;
496 
497  unsigned int row = 0;
498  unsigned int col = 0;
499 
500  for( row = 0 ; row < nRows ; ++row )
501  {
502  for( col = 0 ; col < nCols ; ++col )
503  {
504  gradXMap[ row ][ col ] = 0;
505  gradYMap[ row ][ col ] = 0;
506  }
507  }
508 
509  double value1 = 0;
510  double value2 = 0;
511  double value3 = 0;
512  double value4 = 0;
513  double value5 = 0;
514  double value6 = 0;
515  double value7 = 0;
516  double value8 = 0;
517  unsigned int nextRow = 0;
518  unsigned int nextCol = 0;
519  unsigned int prevRow = 0;
520  unsigned int prevCol = 0;
521  double gradX = 0;
522  double gradY = 0;
523  double mag = 0;
524 
525  for( row = 1 ; row < lastRowIdx ; ++row )
526  {
527  prevRow = row - 1;
528  nextRow = row + 1;
529 
530  for( col = 1 ; col < lastColIdx ; ++col )
531  {
532  prevCol = col - 1;
533  nextCol = col + 1;
534 
535  value1 = inputData[ prevRow ][ prevCol ];
536  value2 = inputData[ prevRow ][ col ];
537  value3 = inputData[ prevRow ][ nextCol ];
538  value4 = inputData[ row ][ prevCol ];
539  value5 = inputData[ row ][ nextCol ];
540  value6 = inputData[ nextRow ][ prevCol ];
541  value7 = inputData[ nextRow ][ col ];
542  value8 = inputData[ nextRow ][ nextCol ];
543 
544  if(
545  ( inputData[ row ][ col ] == DBL_MAX ) ||
546  ( value1 == DBL_MAX ) || ( value2 == DBL_MAX ) ||
547  ( value3 == DBL_MAX ) || ( value4 == DBL_MAX ) ||
548  ( value5 == DBL_MAX ) || ( value6 == DBL_MAX ) ||
549  ( value7 == DBL_MAX ) || ( value8 == DBL_MAX )
550  )
551  {
552  continue;
553  }
554  else
555  {
556  gradX = ( ( value3 + value5 + value8 ) -
557  ( value1 + value4 + value6 ) );
558  gradY = ( ( value1 + value2 + value3 ) -
559  ( value6 + value7 + value8 ) );
560 
561  if( unitVectors && ( gradX != 0.0 ) && ( gradY != 0.0 ) )
562  {
563  mag = std::sqrt( ( gradX * gradX ) + ( gradY * gradY ) );
564  gradX /= mag;
565  gradY /= mag;
566  }
567 
568  gradXMap[ row ][ col ] = gradX;
569  gradYMap[ row ][ col ] = gradY;
570  }
571  }
572  }
573 
574  return true;
575  }
576 
578  const te::rp::Matrix< double >& inputMap,
579  te::rp::Matrix< double >& edgeStrengthMap ) const
580  {
581  const unsigned int nRows = inputMap.getLinesNumber();
582  TERP_TRUE_OR_RETURN_FALSE( nRows > 4, "Invalid rows number" );
583 
584  const unsigned int nCols = inputMap.getColumnsNumber();
585  TERP_TRUE_OR_RETURN_FALSE( nCols > 4, "Invalid columns number" );
586 
587  const unsigned int lastRowIdx = nRows - 1;
588  const unsigned int lastColIdx = nCols - 1;
589 
590  if( ! edgeStrengthMap.reset( nRows,nCols ) )
591  return false;
592 
593  te::rp::Matrix< double > gradXMap;
595  te::rp::Matrix< double > gradYMap;
597  TERP_TRUE_OR_RETURN_FALSE( getGradientMaps( inputMap, false, gradXMap, gradYMap ),
598  "Gradient maps build error" );
599 
600  unsigned int row = 0;
601  unsigned int col = 0;
602 
603  for( row = 0 ; row < nRows ; ++row )
604  {
605  edgeStrengthMap[ row ][ 0 ] = 0;
606  edgeStrengthMap[ row ][ lastColIdx ] = 0;
607  }
608 
609  for( col = 0 ; col < nCols ; ++col )
610  {
611  edgeStrengthMap[ 0 ][ col ] = 0;
612  edgeStrengthMap[ lastRowIdx ][ col ] = 0;
613  }
614 
615  double x1 = 0;
616  double y1 = 0;
617  double x2 = 0;
618  double y2 = 0;
619  double x3 = 0;
620  double y3 = 0;
621  double x4 = 0;
622  double y4 = 0;
623  double x5 = 0;
624  double y5 = 0;
625  double x6 = 0;
626  double y6 = 0;
627  double x7 = 0;
628  double y7 = 0;
629  double x8 = 0;
630  double y8 = 0;
631  unsigned int nextRow = 0;
632  unsigned int nextCol = 0;
633  unsigned int prevRow = 0;
634  unsigned int prevCol = 0;
635  double strengthYUp = 0;
636  double strengthYDown = 0;
637  double strengthXLeft = 0;
638  double strengthXRight = 0;
639  double strength = 0;
640  double diffX = 0;
641  double diffY = 0;
642  double minStrength = DBL_MAX;
643  double maxStrength = -1.0 * DBL_MAX;
644 
645  for( row = 1 ; row < lastRowIdx ; ++row )
646  {
647  prevRow = row - 1;
648  nextRow = row + 1;
649 
650  for( col = 1 ; col < lastColIdx ; ++col )
651  {
652  prevCol = col - 1;
653  nextCol = col + 1;
654 
655  x1 = gradXMap[ prevRow ][ prevCol ];
656  x2 = gradXMap[ prevRow ][ col ];
657  x3 = gradXMap[ prevRow ][ nextCol ];
658  x4 = gradXMap[ row ][ prevCol ];
659  x5 = gradXMap[ row ][ nextCol ];
660  x6 = gradXMap[ nextRow ][ prevCol ];
661  x7 = gradXMap[ nextRow ][ col ];
662  x8 = gradXMap[ nextRow ][ nextCol ];
663 
664  y1 = gradYMap[ prevRow ][ prevCol ];
665  y2 = gradYMap[ prevRow ][ col ];
666  y3 = gradYMap[ prevRow ][ nextCol ];
667  y4 = gradYMap[ row ][ prevCol ];
668  y5 = gradYMap[ row ][ nextCol ];
669  y6 = gradYMap[ nextRow ][ prevCol ];
670  y7 = gradYMap[ nextRow ][ col ];
671  y8 = gradYMap[ nextRow ][ nextCol ];
672 
673  strengthXRight = x1 + x4 + x6;
674  strengthXLeft = x3 + x5 + x8;
675  strengthYUp = y6 + y7 + y8;
676  strengthYDown = y1 + y2 + y3;
677 
678  diffX = std::abs( strengthXRight - strengthXLeft );
679  diffY = std::abs( strengthYUp - strengthYDown );
680 
681  strength = std::max( diffX, diffY );
682 
683  edgeStrengthMap[ row ][ col ] = strength;
684 
685  if( minStrength > strength ) minStrength = strength;
686  if( maxStrength < strength ) maxStrength = strength;
687  }
688  }
689 
690  const double gain = ( minStrength == maxStrength ) ? 0.0 :
691  ( 1.0 / ( maxStrength - minStrength ) );
692 
693  for( row = 1 ; row < lastRowIdx ; ++row )
694  {
695  for( col = 1 ; col < lastColIdx ; ++col )
696  {
697  strength = edgeStrengthMap[ row ][ col ] - minStrength;
698  strength *= gain;
699  edgeStrengthMap[ row ][ col ] = strength;
700  }
701  }
702 
703  return true;
704  }
705 
707  const te::rp::Matrix< double >& inputVecFieldX,
708  const te::rp::Matrix< double >& inputVecFieldY,
709  te::rp::Matrix< double > const * const backGroundMapPtr,
710  const unsigned int vecPixelStep,
711  const std::string& tifFileName ) const
712  {
713  assert( inputVecFieldX.getColumnsNumber() == inputVecFieldY.getColumnsNumber() );
714  assert( inputVecFieldX.getLinesNumber() == inputVecFieldY.getLinesNumber() );
715  assert( backGroundMapPtr ? inputVecFieldX.getColumnsNumber() ==
716  backGroundMapPtr->getColumnsNumber() : true );
717  assert( backGroundMapPtr ? inputVecFieldX.getLinesNumber() ==
718  backGroundMapPtr->getLinesNumber() : true );
719 
720  std::map<std::string, std::string> rInfo;
721  rInfo["URI"] = tifFileName + ".tif";
722 
723  std::vector<te::rst::BandProperty*> bandsProperties;
724  bandsProperties.push_back(new te::rst::BandProperty( 0, te::dt::UCHAR_TYPE, "" ));
725  bandsProperties[0]->m_colorInterp = te::rst::RedCInt;
726  bandsProperties[0]->m_noDataValue = -1.0 * DBL_MAX;
727  bandsProperties.push_back(new te::rst::BandProperty( *bandsProperties[0] ));
728  bandsProperties[0]->m_colorInterp = te::rst::GreenCInt;
729  bandsProperties.push_back(new te::rst::BandProperty( *bandsProperties[0] ));
730  bandsProperties[0]->m_colorInterp = te::rst::BlueCInt;
731 
732  te::rst::Grid* newgrid = new te::rst::Grid( inputVecFieldX.getColumnsNumber(),
733  inputVecFieldX.getLinesNumber(), 0, -1 );
734 
735  std::auto_ptr< te::rst::Raster > outputRasterPtr(
736  te::rst::RasterFactory::make( "GDAL", newgrid, bandsProperties, rInfo, 0, 0));
737  TERP_TRUE_OR_THROW( outputRasterPtr.get(), "Output raster creation error");
738 
739  unsigned int line = 0;
740  unsigned int col = 0;
741  const unsigned int nLines = inputVecFieldX.getLinesNumber();
742  const unsigned int nCols = inputVecFieldX.getColumnsNumber();
743  const unsigned int lastLineIdx = inputVecFieldX.getLinesNumber() - 1;
744  const unsigned int lastColIdx = inputVecFieldX.getColumnsNumber() - 1;
745  double x = 0.0;
746  double y = 0.0;
747 
748  if( backGroundMapPtr )
749  {
750  double minBGValue = DBL_MAX;
751  double maxBGValue = -1.0 * DBL_MAX;
752  double bGValue = 0;
753 
754  for( line = 0 ; line < nLines ; ++line )
755  {
756  for( col = 0 ; col < nCols ; ++col )
757  {
758  bGValue = backGroundMapPtr->operator[]( line )[ col ];
759 
760  if( minBGValue > bGValue )
761  {
762  minBGValue = bGValue;
763  }
764 
765  if( maxBGValue < bGValue )
766  {
767  maxBGValue = bGValue;
768  }
769  }
770  }
771 
772  for( line = 0 ; line < nLines ; ++line )
773  {
774  for( col = 0 ; col < nCols ; ++col )
775  {
776  bGValue = backGroundMapPtr->operator[]( line )[ col ];
777  bGValue -= minBGValue;
778  bGValue *= ( 128.0 / ( maxBGValue - minBGValue ) );
779 
780  outputRasterPtr->setValue( col, line, bGValue, 0 );
781  outputRasterPtr->setValue( col, line, bGValue, 1 );
782  outputRasterPtr->setValue( col, line, bGValue, 2 );
783  }
784  }
785  }
786  else
787  {
788  for( line = 0 ; line < nLines ; ++line )
789  {
790  for( col = 0 ; col < nCols ; ++col )
791  {
792  outputRasterPtr->setValue( col, line, 0, 0 );
793  outputRasterPtr->setValue( col, line, 0, 1 );
794  outputRasterPtr->setValue( col, line, 0, 2 );
795  }
796  }
797  }
798 
799  for( line = 1 ; line < lastLineIdx ; line += vecPixelStep )
800  {
801  for( col = 1 ; col < lastColIdx ; col += vecPixelStep )
802  {
803  x = inputVecFieldX[ line ][ col ];
804  y = inputVecFieldY[ line ][ col ];
805  outputRasterPtr->setValue( col, line, 255, 0 );
806 
807  if( x == 0 )
808  {
809  if( y > 0 )
810  {
811  outputRasterPtr->setValue( col, line - 1, 255, 1 );
812  }
813  else if ( y < 0 )
814  {
815  outputRasterPtr->setValue( col, line + 1, 255, 1 );
816  }
817  }
818  else if( x > 0 )
819  {
820  if( y == 0 )
821  {
822  outputRasterPtr->setValue( col + 1, line, 255, 1 );
823  }
824  else if( y > 0 )
825  {
826  outputRasterPtr->setValue( col + 1, line - 1, 255, 1 );
827  }
828  else
829  {
830  outputRasterPtr->setValue( col + 1, line + 1, 255, 1 );
831  }
832  }
833  else // x < 0
834  {
835  if( y == 0 )
836  {
837  outputRasterPtr->setValue( col - 1, line, 255, 1 );
838  }
839  else if( y > 0 )
840  {
841  outputRasterPtr->setValue( col - 1, line - 1, 255, 1 );
842  }
843  else
844  {
845  outputRasterPtr->setValue( col - 1, line + 1, 255, 1 );
846  }
847  }
848  }
849  }
850  }
851 
853  const te::rp::Matrix< double >& inputX,
854  const te::rp::Matrix< double >& inputY,
855  te::rp::Matrix< double > const * const backgroundDataPtr,
856  te::common::TaskProgress* progressPtr,
857  te::rp::Matrix< double >& outputX,
858  te::rp::Matrix< double >& outputY ) const
859  {
860  assert( inputX.getColumnsNumber() == inputY.getColumnsNumber() );
861  assert( inputX.getLinesNumber() == inputY.getLinesNumber() );
862 
863  const unsigned int nRows = inputX.getLinesNumber();
864  TERP_TRUE_OR_RETURN_FALSE( nRows > 8, "Invalid rows number" );
865 
866  const unsigned int nCols = inputX.getColumnsNumber();
867  TERP_TRUE_OR_RETURN_FALSE( nCols > 8, "Invalid columns number" );
868 
869  const unsigned int rowSizeBytes = sizeof( double ) * nCols;
870 
871  if( ! outputX.reset( nRows,nCols ) )
872  return false;
873  if( ! outputY.reset( nRows,nCols ) )
874  return false;
875 
876  std::auto_ptr< te::rp::Matrix< double > > xBuf1Ptr;
877  xBuf1Ptr.reset( new te::rp::Matrix< double >() );
878  if( ! xBuf1Ptr->reset( nRows,nCols ) )
879  return false;
880 
881  std::auto_ptr< te::rp::Matrix< double > > yBuf1Ptr;
882  yBuf1Ptr.reset( new te::rp::Matrix< double >() );
883  if( ! yBuf1Ptr->reset( nRows,nCols ) )
884  return false;
885 
886  std::auto_ptr< te::rp::Matrix< double > > xBuf2Ptr;
887  xBuf2Ptr.reset( new te::rp::Matrix< double >() );
888  if( ! xBuf2Ptr->reset( nRows,nCols ) )
889  return false;
890 
891  std::auto_ptr< te::rp::Matrix< double > > yBuf2Ptr;
892  yBuf2Ptr.reset( new te::rp::Matrix< double >() );
893  if( ! yBuf2Ptr->reset( nRows,nCols ) )
894  return false;
895 
896  std::auto_ptr< te::rp::Matrix< double > > magBuf1Ptr;
897  magBuf1Ptr.reset( new te::rp::Matrix< double >() );
898  if( ! magBuf1Ptr->reset( nRows,nCols ) )
899  return false;
900 
901  std::auto_ptr< te::rp::Matrix< double > > magBuf2Ptr;
902  magBuf2Ptr.reset( new te::rp::Matrix< double >() );
903  if( ! magBuf2Ptr->reset( nRows,nCols ) )
904  return false;
905  if( ! getMagnitude( inputX, inputY, *magBuf2Ptr ) ) return false;
906 
907  double currentIterationResidue = DBL_MAX;
908 
909  boost::mutex mutex;
910 
911  ApplyVecDiffusionThreadParams threadParams;
912  threadParams.m_initialXBufPtr = &inputX;
913  threadParams.m_initialYBufPtr = &inputY;
914  threadParams.m_currentIterationResiduePtr = &currentIterationResidue;
915  threadParams.m_mutexPtr = &mutex;
917 
918  const unsigned int threadsNumber = m_inputParameters.m_enableMultiThread ?
920  const unsigned int rowsPerThread = m_inputParameters.m_enableMultiThread ?
921  (const unsigned int)(std::ceil( ((double)(nRows - 8 )) / ((double)threadsNumber))) : nRows;
922 
923  unsigned int iteration = 0;
924 
925  while(
926  (
928  ?
930  :
931  true
932  )
933  &&
934  ( currentIterationResidue > m_inputParameters.m_diffusionThreshold )
935  )
936  {
937  if( progressPtr )
938  {
939  if( ! progressPtr->isActive() ) return false;
940  }
941 
942  currentIterationResidue = 0;
943 
944  if( iteration == 0 )
945  {
946  threadParams.m_inputMagBufPtr = magBuf2Ptr.get();
947  threadParams.m_inputBufXPtr = &inputX;
948  threadParams.m_inputBufYPtr = &inputY;
949  threadParams.m_outputMagBufPtr = magBuf1Ptr.get();
950  threadParams.m_outputBufXPtr = xBuf1Ptr.get();
951  threadParams.m_outputBufYPtr = yBuf1Ptr.get();
952  }
953  else
954  {
955  if( threadParams.m_outputBufXPtr == xBuf1Ptr.get() )
956  {
957  threadParams.m_inputMagBufPtr = magBuf1Ptr.get();
958  threadParams.m_inputBufXPtr = xBuf1Ptr.get();
959  threadParams.m_inputBufYPtr = yBuf1Ptr.get();
960  threadParams.m_outputMagBufPtr = magBuf2Ptr.get();
961  threadParams.m_outputBufXPtr = xBuf2Ptr.get();
962  threadParams.m_outputBufYPtr = yBuf2Ptr.get();
963  }
964  else
965  {
966  threadParams.m_inputMagBufPtr = magBuf2Ptr.get();
967  threadParams.m_inputBufXPtr = xBuf2Ptr.get();
968  threadParams.m_inputBufYPtr = yBuf2Ptr.get();
969  threadParams.m_outputMagBufPtr = magBuf1Ptr.get();
970  threadParams.m_outputBufXPtr = xBuf1Ptr.get();
971  threadParams.m_outputBufYPtr = yBuf1Ptr.get();
972  }
973  }
974 
975  memcpy( threadParams.m_outputMagBufPtr->operator[]( 0 ),
976  threadParams.m_inputMagBufPtr->operator[]( 0 ), rowSizeBytes );
977  memcpy( threadParams.m_outputMagBufPtr->operator[]( 1 ),
978  threadParams.m_inputMagBufPtr->operator[]( 1 ), rowSizeBytes );
979  memcpy( threadParams.m_outputMagBufPtr->operator[]( 2 ),
980  threadParams.m_inputMagBufPtr->operator[]( 2 ), rowSizeBytes );
981  memcpy( threadParams.m_outputMagBufPtr->operator[]( 3 ),
982  threadParams.m_inputMagBufPtr->operator[]( 3 ), rowSizeBytes );
983 
984  memcpy( threadParams.m_outputMagBufPtr->operator[]( nRows - 4 ),
985  threadParams.m_inputMagBufPtr->operator[]( nRows - 4 ), rowSizeBytes );
986  memcpy( threadParams.m_outputMagBufPtr->operator[]( nRows - 3 ),
987  threadParams.m_inputMagBufPtr->operator[]( nRows - 3 ), rowSizeBytes );
988  memcpy( threadParams.m_outputMagBufPtr->operator[]( nRows - 2 ),
989  threadParams.m_inputMagBufPtr->operator[]( nRows - 2 ), rowSizeBytes );
990  memcpy( threadParams.m_outputMagBufPtr->operator[]( nRows - 1 ),
991  threadParams.m_inputMagBufPtr->operator[]( nRows - 1 ), rowSizeBytes );
992 
993  memcpy( threadParams.m_outputBufXPtr->operator[]( 0 ),
994  threadParams.m_inputBufXPtr->operator[]( 0 ), rowSizeBytes );
995  memcpy( threadParams.m_outputBufXPtr->operator[]( 1 ),
996  threadParams.m_inputBufXPtr->operator[]( 1 ), rowSizeBytes );
997  memcpy( threadParams.m_outputBufXPtr->operator[]( 2 ),
998  threadParams.m_inputBufXPtr->operator[]( 2 ), rowSizeBytes );
999  memcpy( threadParams.m_outputBufXPtr->operator[]( 3 ),
1000  threadParams.m_inputBufXPtr->operator[]( 3 ), rowSizeBytes );
1001 
1002  memcpy( threadParams.m_outputBufXPtr->operator[]( nRows - 4 ),
1003  threadParams.m_inputBufXPtr->operator[]( nRows - 4 ), rowSizeBytes );
1004  memcpy( threadParams.m_outputBufXPtr->operator[]( nRows - 3 ),
1005  threadParams.m_inputBufXPtr->operator[]( nRows - 3 ), rowSizeBytes );
1006  memcpy( threadParams.m_outputBufXPtr->operator[]( nRows - 2 ),
1007  threadParams.m_inputBufXPtr->operator[]( nRows - 2 ), rowSizeBytes );
1008  memcpy( threadParams.m_outputBufXPtr->operator[]( nRows - 1 ),
1009  threadParams.m_inputBufXPtr->operator[]( nRows - 1 ), rowSizeBytes );
1010 
1011  memcpy( threadParams.m_outputBufYPtr->operator[]( 0 ),
1012  threadParams.m_inputBufYPtr->operator[]( 0 ), rowSizeBytes );
1013  memcpy( threadParams.m_outputBufYPtr->operator[]( 1 ),
1014  threadParams.m_inputBufYPtr->operator[]( 1 ), rowSizeBytes );
1015  memcpy( threadParams.m_outputBufYPtr->operator[]( 2 ),
1016  threadParams.m_inputBufYPtr->operator[]( 2 ), rowSizeBytes );
1017  memcpy( threadParams.m_outputBufYPtr->operator[]( 3 ),
1018  threadParams.m_inputBufYPtr->operator[]( 3 ), rowSizeBytes );
1019 
1020  memcpy( threadParams.m_outputBufYPtr->operator[]( nRows - 4 ),
1021  threadParams.m_inputBufYPtr->operator[]( nRows - 4 ), rowSizeBytes );
1022  memcpy( threadParams.m_outputBufYPtr->operator[]( nRows - 3 ),
1023  threadParams.m_inputBufYPtr->operator[]( nRows - 3 ), rowSizeBytes );
1024  memcpy( threadParams.m_outputBufYPtr->operator[]( nRows - 2 ),
1025  threadParams.m_inputBufYPtr->operator[]( nRows - 2 ), rowSizeBytes );
1026  memcpy( threadParams.m_outputBufYPtr->operator[]( nRows - 1 ),
1027  threadParams.m_inputBufYPtr->operator[]( nRows - 1 ), rowSizeBytes );
1028 
1029  if( threadsNumber )
1030  {
1031  unsigned int startingRow = 4;
1032  std::vector< ApplyVecDiffusionThreadParams > threadsParamesVec;
1033  threadsParamesVec.resize( threadsNumber, threadParams );
1034  boost::thread_group threads;
1035 
1036  for( unsigned int threadIdx = 0 ; threadIdx < threadsNumber ;
1037  ++threadIdx )
1038  {
1039  threadsParamesVec[ threadIdx ].m_firstRowIdx = startingRow;
1040  threadsParamesVec[ threadIdx ].m_lastRowIdx = std::min( nRows - 5,
1041  startingRow + rowsPerThread - 1 );
1042  threads.add_thread( new boost::thread( applyVecDiffusionThreadEntry,
1043  &threadsParamesVec[ threadIdx ] ) );
1044  startingRow += rowsPerThread;
1045  }
1046 
1047  threads.join_all();
1048  }
1049  else
1050  {
1051  threadParams.m_firstRowIdx = 4;
1052  threadParams.m_lastRowIdx = nRows - 5;
1053  applyVecDiffusionThreadEntry( &threadParams );
1054  };
1055 
1056 // if( ( backgroundDataPtr != 0 ) && ( ( iteration < 10 ) || ( iteration % 100 == 0 ) ) )
1057 // {
1058 // CreateRasterFileFromMatrix( *threadParams.m_inputBufXPtr, true,
1059 // boost::lexical_cast< std::string >( iteration ) + "_diffusedInX.tif" );
1060 // CreateRasterFileFromMatrix( *threadParams.m_inputBufYPtr, true,
1061 // boost::lexical_cast< std::string >( iteration ) + "_diffusedInY.tif");
1062 // createTifFromVecField( *threadParams.m_inputBufXPtr, *threadParams.m_inputBufYPtr,
1063 // backgroundDataPtr, 3,
1064 // boost::lexical_cast< std::string >( iteration ) + "_diffusedInVecs");
1065 // CreateRasterFileFromMatrix( *threadParams.m_inputMagBufPtr, true,
1066 // boost::lexical_cast< std::string >( iteration ) + "_diffusedInMag.tif");
1067 //
1068 // CreateRasterFileFromMatrix( *threadParams.m_outputBufXPtr, true,
1069 // boost::lexical_cast< std::string >( iteration ) + "_diffusedOutX.tif" );
1070 // CreateRasterFileFromMatrix( *threadParams.m_outputBufYPtr, true,
1071 // boost::lexical_cast< std::string >( iteration ) + "_diffusedOutY.tif");
1072 // createTifFromVecField( *threadParams.m_outputBufXPtr, *threadParams.m_outputBufYPtr,
1073 // backgroundDataPtr, 3,
1074 // boost::lexical_cast< std::string >( iteration ) + "_diffusedOutVecs");
1075 // CreateRasterFileFromMatrix( *threadParams.m_outputMagBufPtr, true,
1076 // boost::lexical_cast< std::string >( iteration ) + "_diffusedOutMag.tif");
1077 // }
1078 //
1079 // std::cout << std::endl << "currentIteration=" << iteration << std::endl;
1080 // std::cout << std::endl << "currentIterationResidue=" << currentIterationResidue << std::endl;
1081 
1082  ++iteration;
1083  }
1084 
1085  for( unsigned int row = 0 ; row < nRows ; ++row )
1086  {
1087  memcpy( outputX[ row ], threadParams.m_outputBufXPtr->operator[]( row ),
1088  rowSizeBytes );
1089  memcpy( outputY[ row ], threadParams.m_outputBufYPtr->operator[]( row ),
1090  rowSizeBytes );
1091  }
1092 
1093  return true;
1094  }
1095 
1097  {
1098 // const te::rp::Matrix< double >& initBufX = *(paramsPtr->m_initialXBufPtr);
1099 // const te::rp::Matrix< double >& initBufY = *(paramsPtr->m_initialYBufPtr);
1100  const te::rp::Matrix< double >& iBufMag = *(paramsPtr->m_inputMagBufPtr);
1101  const te::rp::Matrix< double >& iBufX = *(paramsPtr->m_inputBufXPtr);
1102  const te::rp::Matrix< double >& iBufY = *(paramsPtr->m_inputBufYPtr);
1103  te::rp::Matrix< double >& oMagBuf = *(paramsPtr->m_outputMagBufPtr);
1104  te::rp::Matrix< double >& oBufX = *(paramsPtr->m_outputBufXPtr);
1105  te::rp::Matrix< double >& oBufY = *(paramsPtr->m_outputBufYPtr);
1106  const double& diffusionRegularization = paramsPtr->m_diffusionRegularization;
1107  const double complementDiffReg = 1.0 - diffusionRegularization;
1108  const unsigned int nCols = iBufX.getColumnsNumber();
1109  const unsigned int rowSizeBytes = sizeof( double ) * nCols;
1110  unsigned int row = 0;
1111  unsigned int col = 0;
1112 
1113  boost::scoped_array< double > bufMagRow0Handler( new double[ nCols ] );
1114  boost::scoped_array< double > bufMagRow1Handler( new double[ nCols ] );
1115  boost::scoped_array< double > bufMagRow2Handler( new double[ nCols ] );
1116  boost::scoped_array< double > bufXRow0Handler( new double[ nCols ] );
1117  boost::scoped_array< double > bufXRow1Handler( new double[ nCols ] );
1118  boost::scoped_array< double > bufXRow2Handler( new double[ nCols ] );
1119  boost::scoped_array< double > bufYRow0Handler( new double[ nCols ] );
1120  boost::scoped_array< double > bufYRow1Handler( new double[ nCols ] );
1121  boost::scoped_array< double > bufYRow2Handler( new double[ nCols ] );
1122 /* boost::scoped_array< double > bufInitXRow0Handler( new double[ nCols ] );
1123  boost::scoped_array< double > bufInitXRow1Handler( new double[ nCols ] );
1124  boost::scoped_array< double > bufInitXRow2Handler( new double[ nCols ] );
1125  boost::scoped_array< double > bufInitYRow0Handler( new double[ nCols ] );
1126  boost::scoped_array< double > bufInitYRow1Handler( new double[ nCols ] );
1127  boost::scoped_array< double > bufInitYRow2Handler( new double[ nCols ] );
1128  double* bufferInitX[ 3 ];
1129  double* bufferInitY[ 3 ]; */
1130  double* bufferMag[ 3 ];
1131  double* bufferX[ 3 ];
1132  double* bufferY[ 3 ];
1133  bufferMag[ 0 ] = bufMagRow0Handler.get();
1134  bufferMag[ 1 ] = bufMagRow1Handler.get();
1135  bufferMag[ 2 ] = bufMagRow2Handler.get();
1136 /* bufferInitX[ 0 ] = bufInitXRow0Handler.get();
1137  bufferInitX[ 1 ] = bufInitXRow1Handler.get();
1138  bufferInitX[ 2 ] = bufInitXRow2Handler.get();
1139  bufferInitY[ 0 ] = bufInitYRow0Handler.get();
1140  bufferInitY[ 1 ] = bufInitYRow1Handler.get();
1141  bufferInitY[ 2 ] = bufInitYRow2Handler.get(); */
1142  bufferX[ 0 ] = bufXRow0Handler.get();
1143  bufferX[ 1 ] = bufXRow1Handler.get();
1144  bufferX[ 2 ] = bufXRow2Handler.get();
1145  bufferY[ 0 ] = bufYRow0Handler.get();
1146  bufferY[ 1 ] = bufYRow1Handler.get();
1147  bufferY[ 2 ] = bufYRow2Handler.get();
1148 
1149  // Loading the two initial rows
1150 
1151  paramsPtr->m_mutexPtr->lock();
1152 
1153  memcpy( bufferMag[ 1 ], iBufMag[ paramsPtr->m_firstRowIdx - 1 ],
1154  rowSizeBytes );
1155  memcpy( bufferMag[ 2 ], iBufMag[ paramsPtr->m_firstRowIdx ],
1156  rowSizeBytes );
1157 
1158 // memcpy( bufferInitX[ 1 ], initBufX[ paramsPtr->m_firstRowIdx - 1 ],
1159 // rowSizeBytes );
1160 // memcpy( bufferInitX[ 2 ], initBufX[ paramsPtr->m_firstRowIdx ],
1161 // rowSizeBytes );
1162 //
1163 // memcpy( bufferInitY[ 1 ], initBufY[ paramsPtr->m_firstRowIdx - 1 ],
1164 // rowSizeBytes );
1165 // memcpy( bufferInitY[ 2 ], initBufY[ paramsPtr->m_firstRowIdx ],
1166 // rowSizeBytes );
1167 
1168  memcpy( bufferX[ 1 ], iBufX[ paramsPtr->m_firstRowIdx - 1 ],
1169  rowSizeBytes );
1170  memcpy( bufferX[ 2 ], iBufX[ paramsPtr->m_firstRowIdx ],
1171  rowSizeBytes );
1172 
1173  memcpy( bufferY[ 1 ], iBufY[ paramsPtr->m_firstRowIdx - 1 ],
1174  rowSizeBytes );
1175  memcpy( bufferY[ 2 ], iBufY[ paramsPtr->m_firstRowIdx ],
1176  rowSizeBytes );
1177 
1178  paramsPtr->m_mutexPtr->unlock();
1179 
1180  const unsigned int rowsBound = paramsPtr->m_lastRowIdx + 2;
1181  const unsigned int colsBound = nCols - 4;
1182 /* double centerInitX = 0;
1183  double centerInitY = 0; */
1184  double centerDiffX = 0;
1185  double centerDiffY = 0;
1186  double centerX = 0;
1187  double centerY = 0;
1188  double centerMag = 0;
1189  double newCenterX = 0;
1190  double newCenterY = 0;
1191  double* rowPtr = 0;
1192  unsigned int prevRow = 0;
1193  unsigned int prevCol = 0;
1194  unsigned int nextCol = 0;
1195  double currentIterationResidue = 0;
1196  double validNeigborsNumber = 0.0;
1197 
1198  boost::scoped_array< double > outRowMagBuffHandler( new double[ nCols ] );
1199  double* outRowMagBuff = outRowMagBuffHandler.get();
1200 
1201  boost::scoped_array< double > outRowXBuffHandler( new double[ nCols ] );
1202  double* outRowXBuff = outRowXBuffHandler.get();
1203 
1204  boost::scoped_array< double > outRowYBuffHandler( new double[ nCols ] );
1205  double* outRowYBuff = outRowYBuffHandler.get();
1206 
1207  for( row = ( paramsPtr->m_firstRowIdx + 1 ) ; row < rowsBound ; ++row )
1208  {
1209  prevRow = row - 1;
1210 
1211  // Rolling-up the internal buffers
1212 
1213 /* rowPtr = bufferInitX[ 0 ];
1214  bufferInitX[ 0 ] = bufferInitX[ 1 ];
1215  bufferInitX[ 1 ] = bufferInitX[ 2 ];
1216  bufferInitX[ 2 ] = rowPtr;
1217 
1218  rowPtr = bufferInitY[ 0 ];
1219  bufferInitY[ 0 ] = bufferInitY[ 1 ];
1220  bufferInitY[ 1 ] = bufferInitY[ 2 ];
1221  bufferInitY[ 2 ] = rowPtr; */
1222 
1223  rowPtr = bufferMag[ 0 ];
1224  bufferMag[ 0 ] = bufferMag[ 1 ];
1225  bufferMag[ 1 ] = bufferMag[ 2 ];
1226  bufferMag[ 2 ] = rowPtr;
1227 
1228  rowPtr = bufferX[ 0 ];
1229  bufferX[ 0 ] = bufferX[ 1 ];
1230  bufferX[ 1 ] = bufferX[ 2 ];
1231  bufferX[ 2 ] = rowPtr;
1232 
1233  rowPtr = bufferY[ 0 ];
1234  bufferY[ 0 ] = bufferY[ 1 ];
1235  bufferY[ 1 ] = bufferY[ 2 ];
1236  bufferY[ 2 ] = rowPtr;
1237 
1238  // Getting a new line
1239 
1240  paramsPtr->m_mutexPtr->lock();
1241 /* memcpy( bufferInitX[ 2 ], initBufX[ row ], rowSizeBytes );
1242  memcpy( bufferInitY[ 2 ], initBufY[ row ], rowSizeBytes ); */
1243  memcpy( bufferMag[ 2 ], iBufMag[ row ], rowSizeBytes );
1244  memcpy( bufferX[ 2 ], iBufX[ row ], rowSizeBytes );
1245  memcpy( bufferY[ 2 ], iBufY[ row ], rowSizeBytes );
1246  paramsPtr->m_mutexPtr->unlock();
1247 
1248  // Diffusing the buffer center line
1249 
1250  for( col = 4 ; col < colsBound ; ++col )
1251  {
1252  prevCol = col - 1;
1253  nextCol = col + 1;
1254 
1255 // centerInitX = bufferInitX[ 1 ][ col ];
1256 // centerInitY = bufferInitY[ 1 ][ col ];
1257 
1258  centerMag = bufferMag[ 1 ][ col ];
1259  centerX = bufferX[ 1 ][ col ];
1260  centerY = bufferY[ 1 ][ col ];
1261 
1262  validNeigborsNumber = 0.0;
1263  newCenterX = 0;
1264  newCenterY = 0;
1265 
1266  DIFFUSENEIGHBOR( 0, prevCol );
1267  DIFFUSENEIGHBOR( 0, col );
1268  DIFFUSENEIGHBOR( 0, nextCol );
1269  DIFFUSENEIGHBOR( 1, prevCol );
1270  DIFFUSENEIGHBOR( 1, nextCol );
1271  DIFFUSENEIGHBOR( 2, prevCol );
1272  DIFFUSENEIGHBOR( 2, col );
1273  DIFFUSENEIGHBOR( 2, nextCol );
1274 
1275  if( validNeigborsNumber != 0.0 )
1276  {
1277  newCenterX /= validNeigborsNumber;
1278  newCenterY /= validNeigborsNumber;
1279 
1280  newCenterX = ( newCenterX * diffusionRegularization ) +
1281  ( centerX * complementDiffReg );
1282  newCenterY = ( newCenterY * diffusionRegularization ) +
1283  ( centerY * complementDiffReg );
1284 
1285  centerDiffX = centerX - newCenterX;
1286  centerDiffY = centerY - newCenterY;
1287 
1288  if( centerDiffX != 0.0 )
1289  {
1290  if( centerX == 0.0 )
1291  {
1292  currentIterationResidue = std::max( currentIterationResidue, 1.0 );
1293  }
1294  else
1295  {
1296  currentIterationResidue =
1297  std::max(
1298  currentIterationResidue
1299  ,
1300  std::abs( centerDiffX / centerX )
1301  );
1302  }
1303  }
1304 
1305  if( centerDiffY != 0.0 )
1306  {
1307  if( centerY == 0.0 )
1308  {
1309  currentIterationResidue = std::max( currentIterationResidue, 1.0 );
1310  }
1311  else
1312  {
1313  currentIterationResidue =
1314  std::max(
1315  currentIterationResidue
1316  ,
1317  std::abs( centerDiffY / centerY )
1318  );
1319  }
1320  }
1321 
1322  outRowXBuff[ col ] = newCenterX;
1323  outRowYBuff[ col ] = newCenterY;
1324  outRowMagBuff[ col ] = std::sqrt( ( newCenterX * newCenterX ) + ( newCenterY *
1325  newCenterY ) );
1326  }
1327  else
1328  {
1329  outRowXBuff[ col ] = centerX;
1330  outRowYBuff[ col ] = centerY;
1331  outRowMagBuff[ col ] = centerMag;
1332  }
1333  }
1334 
1335  // copying the result to output
1336 
1337  paramsPtr->m_mutexPtr->lock();
1338 
1339  memcpy( oMagBuf[ prevRow ], outRowMagBuff, rowSizeBytes );
1340  memcpy( oBufX[ prevRow ], outRowXBuff, rowSizeBytes );
1341  memcpy( oBufY[ prevRow ], outRowYBuff, rowSizeBytes );
1342 
1343  oMagBuf[ prevRow ][ 0 ] = iBufMag[ 1 ][ 0 ];
1344  oMagBuf[ prevRow ][ 1 ] = iBufMag[ 1 ][ 1 ];
1345  oMagBuf[ prevRow ][ 2 ] = iBufMag[ 1 ][ 2 ];
1346  oMagBuf[ prevRow ][ 3 ] = iBufMag[ 1 ][ 3 ];
1347  oMagBuf[ prevRow ][ colsBound ] = iBufMag[ 1 ][ colsBound ];
1348  oMagBuf[ prevRow ][ colsBound + 1 ] = iBufMag[ 1 ][ colsBound + 1 ];
1349  oMagBuf[ prevRow ][ colsBound + 2 ] = iBufMag[ 1 ][ colsBound + 2 ];
1350  oMagBuf[ prevRow ][ colsBound + 3 ] = iBufMag[ 1 ][ colsBound + 3 ];
1351 
1352  oBufX[ prevRow ][ 0 ] = bufferX[ 1 ][ 0 ];
1353  oBufX[ prevRow ][ 1 ] = bufferX[ 1 ][ 1 ];
1354  oBufX[ prevRow ][ 2 ] = bufferX[ 1 ][ 2 ];
1355  oBufX[ prevRow ][ 3 ] = bufferX[ 1 ][ 3 ];
1356  oBufX[ prevRow ][ colsBound ] = bufferX[ 1 ][ colsBound ];
1357  oBufX[ prevRow ][ colsBound + 1 ] = bufferX[ 1 ][ colsBound + 1 ];
1358  oBufX[ prevRow ][ colsBound + 2 ] = bufferX[ 1 ][ colsBound + 2 ];
1359  oBufX[ prevRow ][ colsBound + 3 ] = bufferX[ 1 ][ colsBound + 3 ];
1360 
1361  oBufY[ prevRow ][ 0 ] = bufferY[ 1 ][ 0 ];
1362  oBufY[ prevRow ][ 1 ] = bufferY[ 1 ][ 1 ];
1363  oBufY[ prevRow ][ 2 ] = bufferY[ 1 ][ 2 ];
1364  oBufY[ prevRow ][ 3 ] = bufferY[ 1 ][ 3 ];
1365  oBufY[ prevRow ][ colsBound ] = bufferY[ 1 ][ colsBound ];
1366  oBufY[ prevRow ][ colsBound + 1 ] = bufferY[ 1 ][ colsBound + 1 ];
1367  oBufY[ prevRow ][ colsBound + 2 ] = bufferY[ 1 ][ colsBound + 2 ];
1368  oBufY[ prevRow ][ colsBound + 3 ] = bufferY[ 1 ][ colsBound + 3 ];
1369 
1370  (*paramsPtr->m_currentIterationResiduePtr) =
1371  std::max( (*paramsPtr->m_currentIterationResiduePtr),
1372  currentIterationResidue );
1373 
1374  paramsPtr->m_mutexPtr->unlock();
1375  }
1376  }
1377 
1379  const te::rp::Matrix< double >& inputX,
1380  const te::rp::Matrix< double >& inputY,
1381  const te::rp::Matrix< double >& edgeStrengthMap,
1382  te::rp::Matrix< double >& skelMap ) const
1383  {
1384  assert( inputX.getColumnsNumber() == inputY.getColumnsNumber() );
1385  assert( inputX.getLinesNumber() == inputY.getLinesNumber() );
1386 
1387  const unsigned int nRows = inputX.getLinesNumber();
1388  TERP_TRUE_OR_RETURN_FALSE( nRows > 2, "Invalid rows number" );
1389 
1390  const unsigned int nCols = inputX.getColumnsNumber();
1391  TERP_TRUE_OR_RETURN_FALSE( nCols > 2, "Invalid columns number" );
1392 
1393  if( ! skelMap.reset( nRows,nCols ) )
1394  return false;
1395 
1396  const unsigned int lastRowIdx = nRows - 1;
1397  const unsigned int lastColIdx = nCols - 1;
1398  unsigned int row = 0;
1399  unsigned int col = 0;
1400  unsigned int nextRow = 0;
1401  unsigned int nextCol = 0;
1402  unsigned int prevRow = 0;
1403  unsigned int prevCol = 0;
1404  double centerX = 0;
1405  double centerY = 0;
1406  double neiX = 0;
1407  double neiY = 0;
1408  double diffX = 0;
1409  double diffY = 0;
1410  double diffMag = 0;
1411  double strength = 0;
1412  double minStrength = DBL_MAX;
1413  double maxStrength = -1.0 * DBL_MAX;
1414 
1415  for( row = 0 ; row < nRows ; ++row )
1416  {
1417  skelMap[ row ][ 0 ] = 0;
1418  skelMap[ row ][ lastColIdx ] = 0;
1419  }
1420 
1421  for( col = 0 ; col < nCols ; ++col )
1422  {
1423  skelMap[ 0 ][ col ] = 0;
1424  skelMap[ lastRowIdx ][ col ] = 0;
1425  }
1426 
1427  for( row = 1 ; row < lastRowIdx ; ++row )
1428  {
1429  prevRow = row - 1;
1430  nextRow = row + 1;
1431 
1432  for( col = 1 ; col < lastColIdx ; ++col )
1433  {
1434  prevCol = col - 1;
1435  nextCol = col + 1;
1436 
1437  centerX = inputX[ row ][ col ];
1438  centerY = inputY[ row ][ col ];
1439 
1440  strength = 0;
1441 
1442  SKELSTRENGTHNEIGHBOR( prevRow, prevCol )
1443  SKELSTRENGTHNEIGHBOR( prevRow, col )
1444  SKELSTRENGTHNEIGHBOR( prevRow, nextCol )
1445  SKELSTRENGTHNEIGHBOR( row, prevCol )
1446  SKELSTRENGTHNEIGHBOR( row, nextCol )
1447  SKELSTRENGTHNEIGHBOR( nextRow, prevCol )
1448  SKELSTRENGTHNEIGHBOR( nextRow, col )
1449  SKELSTRENGTHNEIGHBOR( nextRow, nextCol )
1450 
1451  strength /= 8.0;
1452  strength = std::max( 0.0, strength );
1453 
1454  skelMap[ row ][ col ] = strength;
1455 
1456  if( minStrength > strength ) minStrength = strength;
1457  if( maxStrength < strength ) maxStrength = strength;
1458  }
1459  }
1460 
1461  const double gain = ( minStrength == maxStrength ) ? 0.0 :
1462  ( 1.0 / ( maxStrength - minStrength ) );
1463 
1464  for( row = 1 ; row < lastRowIdx ; ++row )
1465  {
1466  for( col = 1 ; col < lastColIdx ; ++col )
1467  {
1468  strength = skelMap[ row ][ col ];
1469  strength -= minStrength;
1470  strength *= gain;
1471  strength = ( strength + ( 1.0 - edgeStrengthMap[ row ][ col ] ) ) / 2.0;
1472  skelMap[ row ][ col ] = strength;
1473  }
1474  }
1475 
1476  return true;
1477  }
1478 
1479  } // end namespace rp
1480 
1481 
1482 } // end namespace te
1483 
bool isActive() const
Verify if the task is active.
te::rp::Matrix< double > const * m_initialXBufPtr
A pointer to the input buffer initial X component.
Definition: Skeleton.h:170
const OutputParameters & operator=(const OutputParameters &params)
Definition: Skeleton.cpp:139
Skeleton::InputParameters m_inputParameters
Input execution parameters.
Definition: Skeleton.h:213
#define DIFFUSENEIGHBOR(neighborRow, neighborCol)
Definition: Skeleton.cpp:40
bool applyVecDiffusion(const te::rp::Matrix< double > &inputX, const te::rp::Matrix< double > &inputY, te::rp::Matrix< double > const *const backgroundDataPtr, te::common::TaskProgress *progressPtr, te::rp::Matrix< double > &outputX, te::rp::Matrix< double > &outputY) const
Apply a vector diffusion over the given vector field.
Definition: Skeleton.cpp:852
std::auto_ptr< te::rst::Raster > m_outputRasterPtr
The generated output registered raster.
Definition: Skeleton.h:126
double m_diffusionRegularization
The diffusion regularization parameter.
Definition: Skeleton.h:182
AbstractParameters * clone() const
Create a clone copy of this instance.
Definition: Skeleton.cpp:150
unsigned int m_firstRowIdx
First row to process.
Definition: Skeleton.h:179
static void applyVecDiffusionThreadEntry(ApplyVecDiffusionThreadParams *paramsPtr)
Vector diffusion thread entry.
Definition: Skeleton.cpp:1096
te::rst::Raster const * m_inputMaskRasterPtr
A pointer to an input raster (where pixels with zero velues will be ignored) or an null pointer if no...
Definition: Skeleton.h:84
Blue channel color interpretation.
Definition: Enums.h:61
unsigned int getNumberOfRows() const
Returns the raster number of rows.
Definition: Raster.cpp:208
unsigned int m_lastRowIdx
Last row to process.
Definition: Skeleton.h:180
bool applyMeanSmooth(const te::rp::Matrix< MatrixElementT > &input, te::rp::Matrix< MatrixElementT > &output) const
Apply a mean filter.
Definition: Skeleton.h:245
double m_diffusionRegularization
A regularization parameter to control the variation from one iteration to the next one (higher values...
Definition: Skeleton.h:88
Raster Processing functions.
AbstractParameters * clone() const
Create a clone copy of this instance.
Definition: Skeleton.cpp:111
te::rp::Matrix< double > const * m_inputMagBufPtr
A pointer to the input magnitude buffer.
Definition: Skeleton.h:172
void clear()
Clear all allocated resources and go back to the initial default parameters.
Definition: Matrix.h:636
const InputParameters & operator=(const InputParameters &params)
Definition: Skeleton.cpp:93
void reset()
Clear all internal allocated resources and reset the parameters instance to its initial state...
Definition: Skeleton.cpp:132
te::rp::Matrix< double > * m_outputMagBufPtr
A pointer to the output magnitude buffer.
Definition: Skeleton.h:175
void reset()
Clear all internal allocated objects and reset the algorithm to its initial state.
Definition: Skeleton.cpp:364
Grid * getGrid()
It returns the raster grid.
Definition: Raster.cpp:94
void createTifFromVecField(const te::rp::Matrix< double > &inputVecFieldX, const te::rp::Matrix< double > &inputVecFieldY, te::rp::Matrix< double > const *const backGroundMapPtr, const unsigned int vecPixelStep, const std::string &tifFileName) const
Create a tiff file from a vector field.
Definition: Skeleton.cpp:706
double * m_currentIterationResiduePtr
A pointer the the current iteration residue;.
Definition: Skeleton.h:181
static Raster * make()
It creates and returns an empty raster with default raster driver.
unsigned int getLinesNumber() const
The number of current matrix lines.
Definition: Matrix.h:665
te::rp::Matrix< double > const * m_inputBufYPtr
A pointer to the input buffer Y component.
Definition: Skeleton.h:174
virtual void getValue(unsigned int c, unsigned int r, double &value) const =0
Returns the cell attribute value.
unsigned int GetPhysProcNumber()
Returns the number of physical processors.
void reset()
Clear all internal allocated resources and reset the parameters instance to its initial state...
Definition: Skeleton.cpp:80
virtual std::size_t getNumberOfBands() const =0
Returns the number of bands (dimension of cells attribute values) in the raster.
#define SKELSTRENGTHNEIGHBOR(neighborRow, neighborCol)
Definition: Skeleton.cpp:48
Skeleton output parameters.
Definition: Skeleton.h:118
bool m_enableMultiThread
Enable (true) the use of threads.
Definition: Skeleton.h:92
#define TERP_TRUE_OR_RETURN_FALSE(value, message)
Checks if value is true. For false values a warning message will be logged and a return of context wi...
Definition: Macros.h:183
bool getGradientMaps(const te::rp::Matrix< double > &inputData, const bool unitVectors, te::rp::Matrix< double > &gradXMap, te::rp::Matrix< double > &gradYMap) const
Create an gradient maps from the input image.
Definition: Skeleton.cpp:478
Red channel color interpretation.
Definition: Enums.h:59
Abstract parameters base interface.
A raster band description.
Definition: Band.h:63
A rectified grid is the spatial support for raster data.
Definition: Grid.h:55
te::rst::Raster const * m_inputRasterPtr
Input raster.
Definition: Skeleton.h:80
boost::mutex * m_mutexPtr
A pointer to the sync mutex.
Definition: Skeleton.h:178
std::string m_rType
Output raster data source type (as described in te::raster::RasterFactory ).
Definition: Skeleton.h:122
te::rp::Matrix< double > * m_outputBufXPtr
A pointer to the output buffer X component.
Definition: Skeleton.h:176
bool createSkeletonStrengthMap(const te::rp::Matrix< double > &inputX, const te::rp::Matrix< double > &inputY, const te::rp::Matrix< double > &edgeStrengthMap, te::rp::Matrix< double > &skelMap) const
Create a skeleton strength map.
Definition: Skeleton.cpp:1378
te::rp::Matrix< double > const * m_initialYBufPtr
A pointer to the input buffer initial Y component.
Definition: Skeleton.h:171
te::rp::Matrix< double > const * m_inputBufXPtr
A pointer to the input buffer X component.
Definition: Skeleton.h:173
bool isInitialized() const
Returns true if the algorithm instance is initialized and ready for execution.
Definition: Skeleton.cpp:436
bool m_enableProgress
Enable/Disable the progress interface (default:false).
Definition: Skeleton.h:96
unsigned int getNumberOfColumns() const
Returns the raster number of columns.
Definition: Raster.cpp:213
const Algorithm & operator=(const Algorithm &)
Definition: Algorithm.cpp:43
unsigned int m_diffusionMaxIterations
The maximum number of iterations to perform (zero: no iterations limit).
Definition: Skeleton.h:90
Skeleton input parameters.
Definition: Skeleton.h:76
A raster band description.
Definition: BandProperty.h:61
#define TERP_TRUE_OR_THROW(value, message)
Checks if value is true and throws an exception if not.
Definition: Macros.h:149
unsigned int getColumnsNumber() const
The number of current matrix columns.
Definition: Matrix.h:672
double m_skeletonThreshold
A threshold to select those pixels as being part of the final skeleton - valid range [0...
Definition: Skeleton.h:94
Raster Processing algorithm output parameters base interface.
bool initialize(const AlgorithmInputParameters &inputParams)
Initialize the algorithm instance making it ready for execution.
Definition: Skeleton.cpp:370
void reset()
Reset (clear) the active instance data.
Definition: Matrix.h:474
bool m_isInitialized
Tells if this instance is initialized.
Definition: Skeleton.h:215
Creation of skeleton imagems.
BandProperty * getProperty()
Returns the band property.
Definition: Band.cpp:370
virtual const Band * getBand(std::size_t i) const =0
Returns the raster i-th band.
This class can be used to inform the progress of a task.
Definition: TaskProgress.h:53
unsigned int m_inputRasterBand
Bands to process from the input raster.
Definition: Skeleton.h:82
bool loadData(te::rp::Matrix< double > &rasterData) const
Load data from the input raster.
Definition: Skeleton.cpp:441
te::rp::Matrix< double > * m_outputBufYPtr
A pointer to the output buffer X component.
Definition: Skeleton.h:177
bool execute(AlgorithmOutputParameters &outputParams)
Executes the algorithm using the supplied parameters.
Definition: Skeleton.cpp:164
Index into a lookup table.
Definition: Enums.h:57
Green channel color interpretation.
Definition: Enums.h:60
double m_diffusionThreshold
A threshold over the residue from one iteration to another, if the residue drops below this value the...
Definition: Skeleton.h:86
std::map< std::string, std::string > m_rInfo
The necessary information to create the output rasters (as described in te::raster::RasterFactory).
Definition: Skeleton.h:124
te::common::AccessPolicy getAccessPolicy() const
Returns the raster access policy.
Definition: Raster.cpp:89
Raster Processing algorithm input parameters base interface.
bool getMagnitude(const te::rp::Matrix< double > &xMap, const te::rp::Matrix< double > &yMap, te::rp::Matrix< MatrixElementT > &magnitude) const
Generate the magnitude map from the input vectors.
Definition: Skeleton.h:311
bool getEdgeStrengthMap(const te::rp::Matrix< double > &inputMap, te::rp::Matrix< double > &edgeStrengthMap) const
Create an Edge strenght Map from the input data.
Definition: Skeleton.cpp:577