All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
Skeleton.cpp
Go to the documentation of this file.
1 /* Copyright (C) 2008 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 
void clear()
Clear all allocated resources and go back to the initial default parameters.
Definition: Matrix.h:642
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:313
te::rp::Matrix< double > const * m_initialYBufPtr
A pointer to the input buffer initial Y component.
Definition: Skeleton.h:173
Index into a lookup table.
Definition: Enums.h:57
virtual void getValue(unsigned int c, unsigned int r, double &value) const =0
Returns the cell attribute value.
#define SKELSTRENGTHNEIGHBOR(neighborRow, neighborCol)
Definition: Skeleton.cpp:48
void reset()
Clear all internal allocated objects and reset the algorithm to its initial state.
Definition: Skeleton.cpp:364
double * m_currentIterationResiduePtr
A pointer the the current iteration residue;.
Definition: Skeleton.h:183
A raster band description.
Definition: BandProperty.h:61
Skeleton::InputParameters m_inputParameters
Input execution parameters.
Definition: Skeleton.h:215
Skeleton output parameters.
Definition: Skeleton.h:120
unsigned int getNumberOfColumns() const
Returns the raster number of columns.
Definition: Raster.cpp:213
virtual const Band * getBand(std::size_t i) const =0
Returns the raster i-th band.
te::rp::Matrix< double > const * m_inputMagBufPtr
A pointer to the input magnitude buffer.
Definition: Skeleton.h:174
bool m_isInitialized
Tells if this instance is initialized.
Definition: Skeleton.h:217
This class can be used to inform the progress of a task.
Definition: TaskProgress.h:53
Raster Processing algorithm output parameters base interface.
double m_diffusionRegularization
A regularization parameter to control the variation from one iteration to the next one (higher values...
Definition: Skeleton.h:90
bool loadData(te::rp::Matrix< double > &rasterData) const
Load data from the input raster.
Definition: Skeleton.cpp:441
Red channel color interpretation.
Definition: Enums.h:59
std::string m_rType
Output raster data source type (as described in te::raster::RasterFactory ).
Definition: Skeleton.h:124
#define DIFFUSENEIGHBOR(neighborRow, neighborCol)
Definition: Skeleton.cpp:40
double m_diffusionRegularization
The diffusion regularization parameter.
Definition: Skeleton.h:184
double m_diffusionThreshold
A threshold over the residue from one iteration to another, if the residue drops below this value the...
Definition: Skeleton.h:88
bool initialize(const AlgorithmInputParameters &inputParams)
Initialize the algorithm instance making it ready for execution.
Definition: Skeleton.cpp:370
bool isActive() const
Verify if the task is active.
void reset()
Clear all internal allocated resources and reset the parameters instance to its initial state...
Definition: Skeleton.cpp:80
te::rp::Matrix< double > const * m_initialXBufPtr
A pointer to the input buffer initial X component.
Definition: Skeleton.h:172
boost::mutex * m_mutexPtr
A pointer to the sync mutex.
Definition: Skeleton.h:180
Raster Processing functions.
bool applyMeanSmooth(const te::rp::Matrix< MatrixElementT > &input, te::rp::Matrix< MatrixElementT > &output) const
Apply a mean filter.
Definition: Skeleton.h:247
#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
const Algorithm & operator=(const Algorithm &)
Definition: Algorithm.cpp:43
te::common::AccessPolicy getAccessPolicy() const
Returns the raster access policy.
Definition: Raster.cpp:89
const InputParameters & operator=(const InputParameters &params)
Definition: Skeleton.cpp:93
unsigned int m_inputRasterBand
Bands to process from the input raster.
Definition: Skeleton.h:84
TECOMMONEXPORT unsigned int GetPhysProcNumber()
Returns the number of physical processors.
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:86
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
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
unsigned int getNumberOfRows() const
Returns the raster number of rows.
Definition: Raster.cpp:208
unsigned int getColumnsNumber() const
The number of current matrix columns.
Definition: Matrix.h:678
unsigned int m_firstRowIdx
First row to process.
Definition: Skeleton.h:181
BandProperty * getProperty()
Returns the band property.
Definition: Band.cpp:428
AbstractParameters * clone() const
Create a clone copy of this instance.
Definition: Skeleton.cpp:150
te::rp::Matrix< double > const * m_inputBufXPtr
A pointer to the input buffer X component.
Definition: Skeleton.h:175
virtual std::size_t getNumberOfBands() const =0
Returns the number of bands (dimension of cells attribute values) in the raster.
te::rp::Matrix< double > * m_outputBufXPtr
A pointer to the output buffer X component.
Definition: Skeleton.h:178
A raster band description.
Definition: Band.h:63
unsigned int m_diffusionMaxIterations
The maximum number of iterations to perform (zero: no iterations limit).
Definition: Skeleton.h:92
te::rp::Matrix< double > const * m_inputBufYPtr
A pointer to the input buffer Y component.
Definition: Skeleton.h:176
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:126
Grid * getGrid()
It returns the raster grid.
Definition: Raster.cpp:94
bool m_enableProgress
Enable/Disable the progress interface (default:false).
Definition: Skeleton.h:98
unsigned int m_lastRowIdx
Last row to process.
Definition: Skeleton.h:182
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
const OutputParameters & operator=(const OutputParameters &params)
Definition: Skeleton.cpp:139
void reset()
Reset (clear) the active instance data.
Definition: Matrix.h:480
static void applyVecDiffusionThreadEntry(ApplyVecDiffusionThreadParams *paramsPtr)
Vector diffusion thread entry.
Definition: Skeleton.cpp:1096
bool isInitialized() const
Returns true if the algorithm instance is initialized and ready for execution.
Definition: Skeleton.cpp:436
Abstract parameters base interface.
te::rp::Matrix< double > * m_outputMagBufPtr
A pointer to the output magnitude buffer.
Definition: Skeleton.h:177
double m_skeletonThreshold
A threshold to select those pixels as being part of the final skeleton - valid range [0...
Definition: Skeleton.h:96
te::rst::Raster const * m_inputRasterPtr
Input raster.
Definition: Skeleton.h:82
AbstractParameters * clone() const
Create a clone copy of this instance.
Definition: Skeleton.cpp:111
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
bool m_enableMultiThread
Enable (true) the use of threads.
Definition: Skeleton.h:94
static Raster * make()
It creates and returns an empty raster with default raster driver.
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
bool execute(AlgorithmOutputParameters &outputParams)
Executes the algorithm using the supplied parameters.
Definition: Skeleton.cpp:164
Raster Processing algorithm input parameters base interface.
void reset()
Clear all internal allocated resources and reset the parameters instance to its initial state...
Definition: Skeleton.cpp:132
Blue channel color interpretation.
Definition: Enums.h:61
A rectified grid is the spatial support for raster data.
Definition: Grid.h:68
Green channel color interpretation.
Definition: Enums.h:60
unsigned int getLinesNumber() const
The number of current matrix lines.
Definition: Matrix.h:671
#define TERP_TRUE_OR_THROW(value, message)
Checks if value is true and throws an exception if not.
Definition: Macros.h:149
Skeleton input parameters.
Definition: Skeleton.h:78
Creation of skeleton imagems.
std::auto_ptr< te::rst::Raster > m_outputRasterPtr
The generated output registered raster.
Definition: Skeleton.h:128
te::rp::Matrix< double > * m_outputBufYPtr
A pointer to the output buffer X component.
Definition: Skeleton.h:179