00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef __itkBoxUtilities_h
00019 #define __itkBoxUtilities_h
00020 #include "itkProgressReporter.h"
00021 #include "itkShapedNeighborhoodIterator.h"
00022 #include "itkImageRegionConstIterator.h"
00023 #include "itkImageRegionIterator.h"
00024 #include "itkConstantBoundaryCondition.h"
00025 #include "itkImageRegionIteratorWithIndex.h"
00026 #include "itkImageRegionIterator.h"
00027 #include "itkImageRegionConstIterator.h"
00028 #include "itkOffset.h"
00029 #include "itkNumericTraits.h"
00030 #include "itkNeighborhoodAlgorithm.h"
00031 #include "itkShapedNeighborhoodIterator.h"
00032 #include "itkZeroFluxNeumannBoundaryCondition.h"
00033
00034 namespace itk {
00035
00036 template< class TIterator >
00037 TIterator*
00038 setConnectivityEarlyBox( TIterator* it, bool fullyConnected=false )
00039 {
00040
00041 typename TIterator::OffsetType offset;
00042 it->ClearActiveList();
00043 if( !fullyConnected)
00044 {
00045
00046
00047 offset.Fill( 0 );
00048 for( unsigned int d=0; d < TIterator::Dimension; ++d )
00049 {
00050 offset[d] = -1;
00051 it->ActivateOffset( offset );
00052 offset[d] = 0;
00053 }
00054 }
00055 else
00056 {
00057
00058
00059 unsigned int centerIndex = it->GetCenterNeighborhoodIndex();
00060 for( unsigned int d=0; d < centerIndex; d++ )
00061 {
00062 offset = it->GetOffset( d );
00063
00064 bool keep = true;
00065 for (unsigned i = 0; i < TIterator::Dimension; i++)
00066 {
00067 if (offset[i] > 0)
00068 {
00069 keep = false;
00070 break;
00071 }
00072 }
00073 if (keep)
00074 {
00075 it->ActivateOffset( offset );
00076 }
00077 }
00078 offset.Fill(0);
00079 it->DeactivateOffset( offset );
00080 }
00081 return it;
00082 }
00083
00084 template <class TInputImage, class TOutputImage>
00085 void
00086 BoxAccumulateFunction(const TInputImage * inputImage,
00087 const TOutputImage * outputImage,
00088 typename TInputImage::RegionType inputRegion,
00089 typename TOutputImage::RegionType outputRegion,
00090 ProgressReporter &progress)
00091 {
00092
00093 typedef TInputImage InputImageType;
00094 typedef typename TInputImage::RegionType RegionType;
00095 typedef typename TInputImage::SizeType SizeType;
00096 typedef typename TInputImage::IndexType IndexType;
00097 typedef typename TInputImage::PixelType PixelType;
00098 typedef typename TInputImage::OffsetType OffsetType;
00099 typedef TOutputImage OutputImageType;
00100 typedef typename TOutputImage::PixelType OutputPixelType;
00101 typedef typename TInputImage::PixelType InputPixelType;
00102
00103 typedef ImageRegionConstIterator<TInputImage> InputIterator;
00104 typedef ImageRegionIterator<TOutputImage> OutputIterator;
00105
00106 typedef ShapedNeighborhoodIterator<TOutputImage> NOutputIterator;
00107 InputIterator inIt( inputImage, inputRegion);
00108 typename TInputImage::SizeType kernelRadius;
00109 kernelRadius.Fill(1);
00110
00111 NOutputIterator noutIt( kernelRadius, outputImage, outputRegion);
00112
00113 setConnectivityEarlyBox( &noutIt, true );
00114
00115 ConstantBoundaryCondition<OutputImageType> oBC;
00116 oBC.SetConstant( NumericTraits< OutputPixelType >::Zero );
00117 noutIt.OverrideBoundaryCondition(&oBC);
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127 std::vector<int> Weights;
00128 typename NOutputIterator::ConstIterator sIt;
00129 for( typename NOutputIterator::IndexListType::const_iterator idxIt = noutIt.GetActiveIndexList().begin();
00130 idxIt != noutIt.GetActiveIndexList().end();
00131 idxIt++ )
00132 {
00133 OffsetType offset = noutIt.GetOffset(*idxIt);
00134 int w = -1;
00135 for (unsigned k = 0; k < InputImageType::ImageDimension; k++)
00136 {
00137 if (offset[k] != 0)
00138 w *= offset[k];
00139 }
00140
00141 Weights.push_back(w);
00142 }
00143
00144 for (inIt.GoToBegin(), noutIt.GoToBegin(); !noutIt.IsAtEnd(); ++inIt, ++noutIt )
00145 {
00146 OutputPixelType Sum = 0;
00147 int k;
00148 for (k = 0, sIt = noutIt.Begin(); !sIt.IsAtEnd();++sIt, ++k)
00149 {
00150 Sum += sIt.Get() * Weights[k];
00151 }
00152 noutIt.SetCenterPixel(Sum + inIt.Get());
00153 progress.CompletedPixel();
00154 }
00155 }
00156
00157
00158 template <class ImType>
00159 std::vector<typename ImType::OffsetType> CornerOffsets(const ImType * im)
00160 {
00161 typedef ShapedNeighborhoodIterator<ImType> NIterator;
00162 typename ImType::SizeType unitradius;
00163 unitradius.Fill(1);
00164 NIterator N1(unitradius, im, im->GetRequestedRegion());
00165 unsigned int centerIndex = N1.GetCenterNeighborhoodIndex();
00166 typename NIterator::OffsetType offset;
00167 std::vector<ITK_TYPENAME ImType::OffsetType> result;
00168 for( unsigned int d=0; d < centerIndex*2 + 1; d++ )
00169 {
00170 offset = N1.GetOffset( d );
00171
00172 bool corner = true;
00173 for (unsigned k = 0; k < ImType::ImageDimension; k++)
00174 {
00175 if (offset[k] == 0)
00176 {
00177 corner = false;
00178 break;
00179 }
00180 }
00181 if (corner)
00182 result.push_back(offset);
00183 }
00184 return(result);
00185 }
00186
00187 template <class TInputImage, class TOutputImage>
00188 void
00189 BoxMeanCalculatorFunction(const TInputImage * accImage,
00190 TOutputImage * outputImage,
00191 typename TInputImage::RegionType inputRegion,
00192 typename TOutputImage::RegionType outputRegion,
00193 typename TInputImage::SizeType Radius,
00194 ProgressReporter &progress)
00195 {
00196
00197 typedef TInputImage InputImageType;
00198 typedef typename TInputImage::RegionType RegionType;
00199 typedef typename TInputImage::SizeType SizeType;
00200 typedef typename TInputImage::IndexType IndexType;
00201 typedef typename TInputImage::PixelType PixelType;
00202 typedef typename TInputImage::OffsetType OffsetType;
00203 typedef typename TInputImage::OffsetValueType OffsetValueType;
00204 typedef TOutputImage OutputImageType;
00205 typedef typename TOutputImage::PixelType OutputPixelType;
00206 typedef typename TInputImage::PixelType InputPixelType;
00207
00208 typedef NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>
00209 FaceCalculatorType;
00210 typedef typename FaceCalculatorType::FaceListType FaceListType;
00211 typedef typename FaceCalculatorType::FaceListType::iterator FaceListTypeIt;
00212 FaceCalculatorType faceCalculator;
00213
00214 FaceListType faceList;
00215 FaceListTypeIt fit;
00216 ZeroFluxNeumannBoundaryCondition<TInputImage> nbc;
00217
00218
00219
00220 SizeType kernelSize;
00221 SizeType internalRadius;
00222 SizeType RegionLimit;
00223
00224 IndexType RegionStart = inputRegion.GetIndex();
00225 for( int i=0; i<TInputImage::ImageDimension; i++ )
00226 {
00227 kernelSize[i] = Radius[i] * 2 + 1;
00228 internalRadius[i] = Radius[i] + 1;
00229 RegionLimit[i] = inputRegion.GetSize()[i] + RegionStart[i] - 1;
00230 }
00231
00232 typedef typename NumericTraits<OutputPixelType>::RealType AccPixType;
00233
00234 std::vector<OffsetType> UnitCorners = CornerOffsets<TInputImage>(accImage);
00235 std::vector<OffsetType> RealCorners;
00236 std::vector<AccPixType> Weights;
00237
00238 for (unsigned k = 0; k < UnitCorners.size(); k++)
00239 {
00240 int prod = 1;
00241 OffsetType ThisCorner;
00242 for (unsigned i = 0; i < TInputImage::ImageDimension; i++)
00243 {
00244 prod *= UnitCorners[k][i];
00245 if (UnitCorners[k][i] > 0)
00246 {
00247 ThisCorner[i] = Radius[i];
00248 }
00249 else
00250 {
00251 ThisCorner[i] = -(Radius[i]+1);
00252 }
00253 }
00254 Weights.push_back((AccPixType)prod);
00255 RealCorners.push_back(ThisCorner);
00256 }
00257
00258
00259 faceList = faceCalculator(accImage, outputRegion, internalRadius);
00260
00261 for (fit = faceList.begin(); fit != faceList.end(); ++fit)
00262 {
00263 if (fit == faceList.begin())
00264 {
00265
00266
00267
00268 AccPixType pixelscount = 1;
00269 for (unsigned i = 0; i < TInputImage::ImageDimension; i++)
00270 {
00271 pixelscount *= (AccPixType)(2*Radius[i] + 1);
00272 }
00273
00274 typedef ImageRegionIterator<OutputImageType> OutputIteratorType;
00275 typedef ImageRegionConstIterator<InputImageType> InputIteratorType;
00276
00277 typedef std::vector<InputIteratorType> CornerItVecType;
00278 CornerItVecType CornerItVec;
00279
00280 for (unsigned k = 0; k < RealCorners.size(); k++)
00281 {
00282 typename InputImageType::RegionType tReg=(*fit);
00283 tReg.SetIndex(tReg.GetIndex() + RealCorners[k]);
00284 InputIteratorType tempIt(accImage, tReg);
00285 tempIt.GoToBegin();
00286 CornerItVec.push_back(tempIt);
00287 }
00288
00289 OutputIteratorType oIt(outputImage, *fit);
00290
00291 for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
00292 {
00293 AccPixType Sum = 0;
00294
00295 for (unsigned k = 0; k < CornerItVec.size(); k++)
00296 {
00297 Sum += Weights[k] * CornerItVec[k].Get();
00298
00299 ++(CornerItVec[k]);
00300 }
00301 oIt.Set(static_cast<OutputPixelType>(Sum/pixelscount));
00302 progress.CompletedPixel();
00303 }
00304 }
00305 else
00306 {
00307
00308 typedef ImageRegionIteratorWithIndex<OutputImageType> OutputIteratorType;
00309 OutputIteratorType oIt(outputImage, *fit);
00310
00311 for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
00312 {
00313
00314
00315
00316 RegionType currentKernelRegion;
00317 currentKernelRegion.SetSize( kernelSize );
00318
00319 IndexType kernelRegionIdx = oIt.GetIndex();
00320 IndexType CentIndex = kernelRegionIdx;
00321 for( int i=0; i<TInputImage::ImageDimension; i++ )
00322 {
00323 kernelRegionIdx[i] -= Radius[i];
00324 }
00325 currentKernelRegion.SetIndex( kernelRegionIdx );
00326 currentKernelRegion.Crop( inputRegion );
00327 long edgepixelscount = currentKernelRegion.GetNumberOfPixels();
00328 AccPixType Sum = 0;
00329
00330
00331
00332
00333
00334
00335
00336
00337 for (unsigned k = 0; k < RealCorners.size(); k++)
00338 {
00339 IndexType ThisCorner = CentIndex + RealCorners[k];
00340 bool IncludeCorner = true;
00341 for (unsigned j = 0; j < TInputImage::ImageDimension; j++)
00342 {
00343 if (UnitCorners[k][j] > 0)
00344 {
00345
00346 if( ThisCorner[j] > static_cast< OffsetValueType>( RegionLimit[j] ) )
00347 {
00348 ThisCorner[j] = static_cast< OffsetValueType>( RegionLimit[j] );
00349 }
00350 }
00351 else
00352 {
00353
00354 if (ThisCorner[j] < RegionStart[j])
00355 {
00356 IncludeCorner = false;
00357 break;
00358 }
00359 }
00360 }
00361 if (IncludeCorner)
00362 {
00363 Sum += accImage->GetPixel(ThisCorner) * Weights[k];
00364 }
00365 }
00366
00367 oIt.Set(static_cast<OutputPixelType>(Sum/(AccPixType)edgepixelscount));
00368 progress.CompletedPixel();
00369 }
00370 }
00371 }
00372 }
00373
00374
00375 template <class TInputImage, class TOutputImage>
00376 void
00377 BoxSigmaCalculatorFunction(const TInputImage * accImage,
00378 TOutputImage * outputImage,
00379 typename TInputImage::RegionType inputRegion,
00380 typename TOutputImage::RegionType outputRegion,
00381 typename TInputImage::SizeType Radius,
00382 ProgressReporter &progress)
00383 {
00384
00385 typedef TInputImage InputImageType;
00386 typedef typename TInputImage::RegionType RegionType;
00387 typedef typename TInputImage::SizeType SizeType;
00388 typedef typename TInputImage::IndexType IndexType;
00389 typedef typename TInputImage::PixelType PixelType;
00390 typedef typename TInputImage::OffsetType OffsetType;
00391 typedef typename TInputImage::OffsetValueType OffsetValueType;
00392 typedef TOutputImage OutputImageType;
00393 typedef typename TOutputImage::PixelType OutputPixelType;
00394 typedef typename TInputImage::PixelType InputPixelType;
00395 typedef typename InputPixelType::ValueType ValueType;
00396
00397 typedef typename NeighborhoodAlgorithm::ImageBoundaryFacesCalculator<InputImageType>
00398 FaceCalculatorType;
00399 typedef typename FaceCalculatorType::FaceListType FaceListType;
00400 typedef typename FaceCalculatorType::FaceListType::iterator FaceListTypeIt;
00401 FaceCalculatorType faceCalculator;
00402
00403 FaceListType faceList;
00404 FaceListTypeIt fit;
00405 ZeroFluxNeumannBoundaryCondition<TInputImage> nbc;
00406
00407
00408
00409 SizeType kernelSize;
00410 SizeType internalRadius;
00411 SizeType RegionLimit;
00412 IndexType RegionStart = inputRegion.GetIndex();
00413 for( int i=0; i<TInputImage::ImageDimension; i++ )
00414 {
00415 kernelSize[i] = Radius[i] * 2 + 1;
00416 internalRadius[i] = Radius[i] + 1;
00417 RegionLimit[i] = inputRegion.GetSize()[i] + RegionStart[i] - 1;
00418 }
00419
00420 typedef typename NumericTraits<OutputPixelType>::RealType AccPixType;
00421
00422 std::vector<OffsetType> UnitCorners = CornerOffsets<TInputImage>(accImage);
00423 std::vector<OffsetType> RealCorners;
00424 std::vector<AccPixType> Weights;
00425
00426 for (unsigned k = 0; k < UnitCorners.size(); k++)
00427 {
00428 int prod = 1;
00429 OffsetType ThisCorner;
00430 for (unsigned i = 0; i < TInputImage::ImageDimension; i++)
00431 {
00432 prod *= UnitCorners[k][i];
00433 if (UnitCorners[k][i] > 0)
00434 {
00435 ThisCorner[i] = Radius[i];
00436 }
00437 else
00438 {
00439 ThisCorner[i] = -(Radius[i]+1);
00440 }
00441 }
00442 Weights.push_back((AccPixType)prod);
00443 RealCorners.push_back(ThisCorner);
00444 }
00445
00446
00447 faceList = faceCalculator(accImage, outputRegion, internalRadius);
00448
00449 for (fit = faceList.begin(); fit != faceList.end(); ++fit)
00450 {
00451 if (fit == faceList.begin())
00452 {
00453
00454
00455
00456 AccPixType pixelscount = 1;
00457 for (unsigned i = 0; i < TInputImage::ImageDimension; i++)
00458 {
00459 pixelscount *= (AccPixType)(2*Radius[i] + 1);
00460 }
00461
00462 typedef ImageRegionIterator<OutputImageType> OutputIteratorType;
00463 typedef ImageRegionConstIterator<InputImageType> InputIteratorType;
00464
00465 typedef std::vector<InputIteratorType> CornerItVecType;
00466 CornerItVecType CornerItVec;
00467
00468 for (unsigned k = 0; k < RealCorners.size(); k++)
00469 {
00470 typename InputImageType::RegionType tReg=(*fit);
00471 tReg.SetIndex(tReg.GetIndex() + RealCorners[k]);
00472 InputIteratorType tempIt(accImage, tReg);
00473 tempIt.GoToBegin();
00474 CornerItVec.push_back(tempIt);
00475 }
00476
00477 OutputIteratorType oIt(outputImage, *fit);
00478
00479 for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
00480 {
00481 AccPixType Sum = 0;
00482 AccPixType SquareSum = 0;
00483
00484 for (unsigned k = 0; k < CornerItVec.size(); k++)
00485 {
00486 const InputPixelType & i = CornerItVec[k].Get();
00487 Sum += Weights[k] * i[0];
00488 SquareSum += Weights[k] * i[1];
00489
00490 ++(CornerItVec[k]);
00491 }
00492
00493 oIt.Set(static_cast<OutputPixelType>( vcl_sqrt( ( SquareSum - Sum*Sum/pixelscount ) / ( pixelscount -1 ) ) ) );
00494 progress.CompletedPixel();
00495 }
00496 }
00497 else
00498 {
00499
00500 typedef ImageRegionIteratorWithIndex<OutputImageType> OutputIteratorType;
00501 OutputIteratorType oIt(outputImage, *fit);
00502
00503 for (oIt.GoToBegin(); !oIt.IsAtEnd(); ++oIt)
00504 {
00505
00506
00507
00508 RegionType currentKernelRegion;
00509 currentKernelRegion.SetSize( kernelSize );
00510
00511 IndexType kernelRegionIdx = oIt.GetIndex();
00512 IndexType CentIndex = kernelRegionIdx;
00513 for( int i=0; i<TInputImage::ImageDimension; i++ )
00514 {
00515 kernelRegionIdx[i] -= Radius[i];
00516 }
00517 currentKernelRegion.SetIndex( kernelRegionIdx );
00518 currentKernelRegion.Crop( inputRegion );
00519 long edgepixelscount = currentKernelRegion.GetNumberOfPixels();
00520 AccPixType Sum = 0;
00521 AccPixType SquareSum = 0;
00522
00523
00524
00525
00526
00527
00528
00529
00530 for (unsigned k = 0; k < RealCorners.size(); k++)
00531 {
00532 IndexType ThisCorner = CentIndex + RealCorners[k];
00533 bool IncludeCorner = true;
00534 for (unsigned j = 0; j < TInputImage::ImageDimension; j++)
00535 {
00536 if (UnitCorners[k][j] > 0)
00537 {
00538
00539 if( ThisCorner[j] > static_cast< OffsetValueType >( RegionLimit[j]) )
00540 {
00541 ThisCorner[j] = static_cast< OffsetValueType>( RegionLimit[j] );
00542 }
00543 }
00544 else
00545 {
00546
00547 if (ThisCorner[j] < RegionStart[j])
00548 {
00549 IncludeCorner = false;
00550 break;
00551 }
00552 }
00553 }
00554 if (IncludeCorner)
00555 {
00556 const InputPixelType & i = accImage->GetPixel(ThisCorner);
00557 Sum += Weights[k] * i[0];
00558 SquareSum += Weights[k] * i[1];
00559 }
00560 }
00561
00562 oIt.Set(static_cast<OutputPixelType>( vcl_sqrt( ( SquareSum - Sum*Sum/edgepixelscount ) / ( edgepixelscount -1 ) ) ) );
00563 progress.CompletedPixel();
00564 }
00565 }
00566 }
00567 }
00568
00569
00570 template <class TInputImage, class TOutputImage>
00571 void
00572 BoxSquareAccumulateFunction(const TInputImage * inputImage,
00573 TOutputImage * outputImage,
00574 typename TInputImage::RegionType inputRegion,
00575 typename TOutputImage::RegionType outputRegion,
00576 ProgressReporter &progress)
00577 {
00578
00579 typedef TInputImage InputImageType;
00580 typedef typename TInputImage::RegionType RegionType;
00581 typedef typename TInputImage::SizeType SizeType;
00582 typedef typename TInputImage::IndexType IndexType;
00583 typedef typename TInputImage::PixelType PixelType;
00584 typedef typename TInputImage::OffsetType OffsetType;
00585 typedef TOutputImage OutputImageType;
00586 typedef typename TOutputImage::PixelType OutputPixelType;
00587 typedef typename OutputPixelType::ValueType ValueType;
00588 typedef typename TInputImage::PixelType InputPixelType;
00589
00590 typedef ImageRegionConstIterator<TInputImage> InputIterator;
00591 typedef ImageRegionIterator<TOutputImage> OutputIterator;
00592
00593 typedef ShapedNeighborhoodIterator<TOutputImage> NOutputIterator;
00594 InputIterator inIt( inputImage, inputRegion);
00595 typename TInputImage::SizeType kernelRadius;
00596 kernelRadius.Fill(1);
00597
00598 NOutputIterator noutIt( kernelRadius, outputImage, outputRegion);
00599
00600 setConnectivityEarlyBox( &noutIt, true );
00601
00602 ConstantBoundaryCondition<OutputImageType> oBC;
00603 oBC.SetConstant( NumericTraits< OutputPixelType >::Zero );
00604 noutIt.OverrideBoundaryCondition(&oBC);
00605
00606
00607
00608
00609
00610
00611
00612
00613
00614 std::vector<int> Weights;
00615 typename NOutputIterator::ConstIterator sIt;
00616 for( typename NOutputIterator::IndexListType::const_iterator idxIt = noutIt.GetActiveIndexList().begin();
00617 idxIt != noutIt.GetActiveIndexList().end();
00618 idxIt++ )
00619 {
00620 OffsetType offset = noutIt.GetOffset(*idxIt);
00621 int w = -1;
00622 for (unsigned k = 0; k < InputImageType::ImageDimension; k++)
00623 {
00624 if (offset[k] != 0)
00625 w *= offset[k];
00626 }
00627 Weights.push_back(w);
00628 }
00629
00630 for (inIt.GoToBegin(), noutIt.GoToBegin(); !noutIt.IsAtEnd(); ++inIt, ++noutIt )
00631 {
00632 ValueType Sum = 0;
00633 ValueType SquareSum = 0;
00634 int k;
00635 for (k = 0, sIt = noutIt.Begin(); !sIt.IsAtEnd();++sIt, ++k)
00636 {
00637 const OutputPixelType & v = sIt.Get();
00638 Sum += v[0] * Weights[k];
00639 SquareSum += v[1] * Weights[k];
00640 }
00641 OutputPixelType o;
00642 const InputPixelType & i = inIt.Get();
00643 o[0] = Sum + i;
00644 o[1] = SquareSum + i*i;
00645 noutIt.SetCenterPixel( o );
00646 progress.CompletedPixel();
00647 }
00648 }
00649
00650
00651 }
00652
00653 #endif
00654