vtkImageIslandRemoval2D.cxx 19.3 KB
Newer Older
1 2 3
/*=========================================================================

  Program:   Visualization Toolkit
Charles Law's avatar
Charles Law committed
4
  Module:    vtkImageIslandRemoval2D.cxx
5

6
  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 8
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
Will Schroeder's avatar
Will Schroeder committed
9

10 11
     This software is distributed WITHOUT ANY WARRANTY; without even
     the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12
     PURPOSE.  See the above copyright notice for more information.
13 14

=========================================================================*/
Charles Law's avatar
Charles Law committed
15
#include "vtkImageIslandRemoval2D.h"
16 17

#include "vtkImageData.h"
Amy Squillacote's avatar
Amy Squillacote committed
18 19
#include "vtkInformation.h"
#include "vtkInformationVector.h"
20
#include "vtkObjectFactory.h"
Amy Squillacote's avatar
Amy Squillacote committed
21
#include "vtkStreamingDemandDrivenPipeline.h"
22

Brad King's avatar
Brad King committed
23
vtkStandardNewMacro(vtkImageIslandRemoval2D);
24 25 26

//----------------------------------------------------------------------------
// Constructor: Sets default filter to be identity.
Charles Law's avatar
Charles Law committed
27
vtkImageIslandRemoval2D::vtkImageIslandRemoval2D()
28
{
29
  this->AreaThreshold = 0;
30
  this->SetAreaThreshold(4);
31
  this->SquareNeighborhood = 1;
32
  this->SquareNeighborhoodOff();
33
  this->ReplaceValue = 0;
34
  this->SetReplaceValue(255);
35
  this->IslandValue = 255;
36 37 38 39
  this->SetIslandValue(0);
}

//----------------------------------------------------------------------------
40
void vtkImageIslandRemoval2D::PrintSelf(ostream& os, vtkIndent indent)
41
{
Brad King's avatar
Brad King committed
42
  this->Superclass::PrintSelf(os,indent);
43 44 45 46 47 48 49 50 51 52 53
  os << indent << "AreaThreshold: " << this->AreaThreshold;
  if (this->SquareNeighborhood)
    {
    os << indent << "Neighborhood: Square";
    }
  else
    {
    os << indent << "Neighborhood: Cross";
    }
  os << indent << "IslandValue: " << this->IslandValue;
  os << indent << "ReplaceValue: " << this->ReplaceValue;
54

55 56
}

57

58 59 60 61


//----------------------------------------------------------------------------
// The templated execute function handles all the data types.
62
// Codes:  0 => unvisited. 1 => visted don't know.
63 64 65 66 67
//         2 => visted keep.  3 => visited replace.
// Please excuse the length of this function.  The easiest way to choose
// neighborhoods is to check neighbors one by one directly.  Also, I did
// not want to break the templated function into pieces.
template <class T>
Ken Martin's avatar
Ken Martin committed
68 69 70 71
void vtkImageIslandRemoval2DExecute(vtkImageIslandRemoval2D *self,
                                    vtkImageData *inData, T *inPtr,
                                    vtkImageData *outData, T *outPtr,
                                    int outExt[6])
72
{
73
  int outIdx0, outIdx1, outIdx2;
74
  vtkIdType outInc0, outInc1, outInc2;
75
  T *outPtr0, *outPtr1, *outPtr2;
76
  vtkIdType inInc0, inInc1, inInc2;
77
  T *inPtr0, *inPtr1, *inPtr2;
Charles Law's avatar
Charles Law committed
78
  vtkImage2DIslandPixel *pixels;  // All the pixels visited so far.
79
  int numPixels;      // The number of pixels visited so far.
Charles Law's avatar
Charles Law committed
80
  vtkImage2DIslandPixel *newPixel;   // The last pixel in the list.
81
  int nextPixelIdx;   // The index of the next pixel to grow.
Charles Law's avatar
Charles Law committed
82
  vtkImage2DIslandPixel *nextPixel;  // The next pixel to grow.
83 84 85 86 87 88
  int keepValue;
  int area;
  int squareNeighborhood;
  T islandValue;
  T replaceValue;
  T *inNeighborPtr, *outNeighborPtr;
89
  int idxC, maxC;
90 91
  unsigned long count = 0;
  unsigned long target;
92

93 94
  squareNeighborhood = self->GetSquareNeighborhood();
  area = self->GetAreaThreshold();
95 96
  islandValue = static_cast<T>(self->GetIslandValue());
  replaceValue = static_cast<T>(self->GetReplaceValue());
97

98 99 100
  outData->GetIncrements(outInc0, outInc1, outInc2);
  inData->GetIncrements(inInc0, inInc1, inInc2);
  maxC = outData->GetNumberOfScalarComponents();
101

102
  // Loop through pixels setting all output to 0 (unvisited).
103
  for (idxC = 0; idxC < maxC; idxC++)
104
    {
105 106
    outPtr2 = outPtr + idxC;
    for (outIdx2 = outExt[4]; outIdx2 <= outExt[5]; ++outIdx2)
107
      {
108 109
      outPtr1 = outPtr2;
      for (outIdx1 = outExt[2]; outIdx1 <= outExt[3]; ++outIdx1)
110 111 112 113
        {
        outPtr0 = outPtr1;
        for (outIdx0 = outExt[0]; outIdx0 <= outExt[1]; ++outIdx0)
          {
114
          *outPtr0 = static_cast<T>(0);  // Unvisited
115 116 117 118
          outPtr0 += outInc0;
          }
        outPtr1 += outInc1;
        }
119
      outPtr2 += outInc2;
120 121
      }
    }
122 123 124 125 126 127 128

  // update the progress and possibly abort
  self->UpdateProgress(0.1);
  if (self->AbortExecute)
    {
    return;
    }
129

130
  // In case all 8 neighbors get added before we test the number.
131 132
  pixels = new vtkImage2DIslandPixel [area + 8];

133 134
  target = static_cast<unsigned long>(maxC*(outExt[5]-outExt[4]+1)*
                                      (outExt[3]-outExt[2]+1)/50.0);
135 136
  target++;

137
  // Loop though all pixels looking for islands
138
  for (idxC = 0; idxC < maxC; idxC++)
139
    {
140 141
    outPtr2 = outPtr + idxC;
    inPtr2 = inPtr + idxC;
142
    for (outIdx2 = outExt[4];
143
         !self->AbortExecute && outIdx2 <= outExt[5]; ++outIdx2)
144
      {
145
      if (!(count%target))
146 147 148
        {
        self->UpdateProgress(0.1 + 0.8*count/(50.0*target));
        }
149
      count++;
150 151 152
      outPtr1 = outPtr2;
      inPtr1 = inPtr2;
      for (outIdx1 = outExt[2]; outIdx1 <= outExt[3]; ++outIdx1)
153 154 155 156 157 158 159 160 161 162 163 164 165 166
        {
        outPtr0 = outPtr1;
        inPtr0 = inPtr1;
        for (outIdx0 = outExt[0]; outIdx0 <= outExt[1]; ++outIdx0)
          {
          if (*outPtr0 == 0)
            {
            if (*inPtr0 != islandValue)
              {
              // not an island, keep and go on
              *outPtr0 = 2;
              }
            else
              {
167

168 169 170
              // Start an island search
              // Save first pixel.
              newPixel = pixels;
171 172
              newPixel->inPtr = static_cast<void *>(inPtr0);
              newPixel->outPtr = static_cast<void *>(outPtr0);
173 174 175 176 177 178 179 180 181 182 183 184 185 186
              newPixel->idx0 = outIdx0;
              newPixel->idx1 = outIdx1;
              numPixels = 1;
              nextPixelIdx = 0;
              nextPixel = pixels;
              *outPtr0 = 1;  // visited don't know
              keepValue = 1;
              // breadth first search
              while (keepValue == 1)  // don't know
                {
                // Check all the neighbors
                // left
                if (nextPixel->idx0 > outExt[0])
                  {
187
                  inNeighborPtr = static_cast<T *>(nextPixel->inPtr) - inInc0;
188 189
                  if ( *inNeighborPtr == islandValue)
                    {
190
                    outNeighborPtr = static_cast<T *>(nextPixel->outPtr) - outInc0;
191 192 193 194 195 196 197 198 199 200
                    if ( *outNeighborPtr == 2)
                      {
                      // This is part of a bigger island.
                      keepValue = 2;
                      }
                    if ( *outNeighborPtr == 0)
                      {
                      // New pixel to add
                      ++newPixel;
                      ++numPixels;
201 202
                      newPixel->inPtr = static_cast<void *>(inNeighborPtr);
                      newPixel->outPtr = static_cast<void *>(outNeighborPtr);
203 204 205 206 207 208 209 210 211
                      newPixel->idx0 = nextPixel->idx0 - 1;
                      newPixel->idx1 = nextPixel->idx1;
                      *outNeighborPtr = 1;  // visited don't know
                      }
                    }
                  }
                // right
                if (nextPixel->idx0 < outExt[1])
                  {
212
                  inNeighborPtr = static_cast<T *>(nextPixel->inPtr) + inInc0;
213 214
                  if ( *inNeighborPtr == islandValue)
                    {
215 216
                    outNeighborPtr = static_cast<T *>(nextPixel->outPtr)
                      + outInc0;
217 218 219 220 221 222 223 224 225 226
                    if ( *outNeighborPtr == 2)
                      {
                      // This is part of a bigger island.
                      keepValue = 2;
                      }
                    if ( *outNeighborPtr == 0)
                      {
                      // New pixel to add
                      ++newPixel;
                      ++numPixels;
227 228
                      newPixel->inPtr = static_cast<void *>(inNeighborPtr);
                      newPixel->outPtr = static_cast<void *>(outNeighborPtr);
229 230 231 232 233 234 235 236 237
                      newPixel->idx0 = nextPixel->idx0 + 1;
                      newPixel->idx1 = nextPixel->idx1;
                      *outNeighborPtr = 1;  // visited don't know
                      }
                    }
                  }
                // up
                if (nextPixel->idx1 > outExt[2])
                  {
238
                  inNeighborPtr = static_cast<T *>(nextPixel->inPtr) - inInc1;
239 240
                  if ( *inNeighborPtr == islandValue)
                    {
241
                    outNeighborPtr = static_cast<T *>(nextPixel->outPtr) - outInc1;
242 243 244 245 246 247 248 249 250 251
                    if ( *outNeighborPtr == 2)
                      {
                      // This is part of a bigger island.
                      keepValue = 2;
                      }
                    if ( *outNeighborPtr == 0)
                      {
                      // New pixel to add
                      ++newPixel;
                      ++numPixels;
252 253
                      newPixel->inPtr = static_cast<void *>(inNeighborPtr);
                      newPixel->outPtr = static_cast<void *>(outNeighborPtr);
254 255 256 257 258 259 260 261 262
                      newPixel->idx0 = nextPixel->idx0;
                      newPixel->idx1 = nextPixel->idx1 - 1;
                      *outNeighborPtr = 1;  // visited don't know
                      }
                    }
                  }
                // down
                if (nextPixel->idx1 < outExt[3])
                  {
263
                  inNeighborPtr = static_cast<T *>(nextPixel->inPtr) + inInc1;
264 265
                  if ( *inNeighborPtr == islandValue)
                    {
266
                    outNeighborPtr = static_cast<T *>(nextPixel->outPtr) + outInc1;
267 268 269 270 271 272 273 274 275 276
                    if ( *outNeighborPtr == 2)
                      {
                      // This is part of a bigger island.
                      keepValue = 2;
                      }
                    if ( *outNeighborPtr == 0)
                      {
                      // New pixel to add
                      ++newPixel;
                      ++numPixels;
277 278
                      newPixel->inPtr = static_cast<void *>(inNeighborPtr);
                      newPixel->outPtr = static_cast<void *>(outNeighborPtr);
279 280 281 282 283 284 285 286 287 288 289 290
                      newPixel->idx0 = nextPixel->idx0;
                      newPixel->idx1 = nextPixel->idx1 + 1;
                      *outNeighborPtr = 1;  // visited don't know
                      }
                    }
                  }
                // Corners
                if (squareNeighborhood)
                  {
                  // upper left
                  if (nextPixel->idx0 > outExt[0] && nextPixel->idx1 > outExt[2])
                    {
291
                    inNeighborPtr = static_cast<T *>(nextPixel->inPtr) - inInc0 - inInc1;
292 293
                    if ( *inNeighborPtr == islandValue)
                      {
294
                      outNeighborPtr = static_cast<T *>(nextPixel->outPtr) - outInc0 -outInc1;
295 296 297 298 299 300 301 302 303 304
                      if ( *outNeighborPtr == 2)
                        {
                        // This is part of a bigger island.
                        keepValue = 2;
                        }
                      if ( *outNeighborPtr == 0)
                        {
                        // New pixel to add
                        ++newPixel;
                        ++numPixels;
305 306
                        newPixel->inPtr = static_cast<void *>(inNeighborPtr);
                        newPixel->outPtr = static_cast<void *>(outNeighborPtr);
307 308 309 310 311 312 313 314 315
                        newPixel->idx0 = nextPixel->idx0 - 1;
                        newPixel->idx1 = nextPixel->idx1 - 1;
                        *outNeighborPtr = 1;  // visited don't know
                        }
                      }
                    }
                  // upper right
                  if (nextPixel->idx0 < outExt[1] && nextPixel->idx1 > outExt[2])
                    {
316
                    inNeighborPtr = static_cast<T *>(nextPixel->inPtr) + inInc0 - inInc1;
317 318
                    if ( *inNeighborPtr == islandValue)
                      {
319
                      outNeighborPtr = static_cast<T *>(nextPixel->outPtr) + outInc0 -outInc1;
320 321 322 323 324 325 326 327 328 329
                      if ( *outNeighborPtr == 2)
                        {
                        // This is part of a bigger island.
                        keepValue = 2;
                        }
                      if ( *outNeighborPtr == 0)
                        {
                        // New pixel to add
                        ++newPixel;
                        ++numPixels;
330 331
                        newPixel->inPtr = static_cast<void *>(inNeighborPtr);
                        newPixel->outPtr = static_cast<void *>(outNeighborPtr);
332 333 334 335 336 337 338 339 340
                        newPixel->idx0 = nextPixel->idx0 + 1;
                        newPixel->idx1 = nextPixel->idx1 - 1;
                        *outNeighborPtr = 1;  // visited don't know
                        }
                      }
                    }
                  // lower left
                  if (nextPixel->idx0 > outExt[0] && nextPixel->idx1 < outExt[3])
                    {
341
                    inNeighborPtr = static_cast<T *>(nextPixel->inPtr) - inInc0 + inInc1;
342 343
                    if ( *inNeighborPtr == islandValue)
                      {
344
                      outNeighborPtr = static_cast<T *>(nextPixel->outPtr) - outInc0 +outInc1;
345 346 347 348 349 350 351 352 353 354
                      if ( *outNeighborPtr == 2)
                        {
                        // This is part of a bigger island.
                        keepValue = 2;
                        }
                      if ( *outNeighborPtr == 0)
                        {
                        // New pixel to add
                        ++newPixel;
                        ++numPixels;
355 356
                        newPixel->inPtr = static_cast<void *>(inNeighborPtr);
                        newPixel->outPtr = static_cast<void *>(outNeighborPtr);
357 358 359 360 361 362 363 364 365
                        newPixel->idx0 = nextPixel->idx0 - 1;
                        newPixel->idx1 = nextPixel->idx1 + 1;
                        *outNeighborPtr = 1;  // visited don't know
                        }
                      }
                    }
                  // lower right
                  if (nextPixel->idx0 < outExt[1] && nextPixel->idx1 < outExt[3])
                    {
366
                    inNeighborPtr = static_cast<T *>(nextPixel->inPtr) + inInc0 + inInc1;
367 368
                    if ( *inNeighborPtr == islandValue)
                      {
369
                      outNeighborPtr = static_cast<T *>(nextPixel->outPtr) + outInc0 +outInc1;
370 371 372 373 374 375 376 377 378 379
                      if ( *outNeighborPtr == 2)
                        {
                        // This is part of a bigger island.
                        keepValue = 2;
                        }
                      if ( *outNeighborPtr == 0)
                        {
                        // New pixel to add
                        ++newPixel;
                        ++numPixels;
380 381
                        newPixel->inPtr = static_cast<void *>(inNeighborPtr);
                        newPixel->outPtr = static_cast<void *>(outNeighborPtr);
382 383 384 385 386 387 388
                        newPixel->idx0 = nextPixel->idx0 + 1;
                        newPixel->idx1 = nextPixel->idx1 + 1;
                        *outNeighborPtr = 1;  // visited don't know
                        }
                      }
                    }
                  }
389

390 391 392
                // Move to the next pixel to grow.
                ++nextPixel;
                ++nextPixelIdx;
393

394 395 396 397 398
                // Have we visted enogh pixels to determine this is a keeper?
                if (keepValue == 1 && numPixels >= area)
                  {
                  keepValue = 2;
                  }
399

400 401 402 403 404 405 406
                // Have we run out of pixels to grow?
                if (keepValue == 1 && nextPixelIdx >= numPixels)
                  {
                  // The island is too small. Set island values too replace.
                  keepValue = 3;
                  }
                }
407

408 409 410 411
              // Change "don't knows" to keep value
              nextPixel = pixels;
              for (nextPixelIdx = 0; nextPixelIdx < numPixels; ++nextPixelIdx)
                {
412
                *(static_cast<T *>(nextPixel->outPtr)) = keepValue;
413 414 415 416
                ++nextPixel;
                }
              }
            }
417

418 419 420 421 422 423
          outPtr0 += outInc0;
          inPtr0 += inInc0;
          }
        outPtr1 += outInc1;
        inPtr1 += inInc1;
        }
424 425
      outPtr2 += outInc2;
      inPtr2 += inInc2;
426 427
      }
    }
428

429
  delete [] pixels;
430

431 432 433 434 435 436 437
  // update the progress and possibly abort
  self->UpdateProgress(0.9);
  if (self->AbortExecute)
    {
    return;
    }

438
  // Loop though all pixels actually copying and replacing.
439
  for (idxC = 0; idxC < maxC; idxC++)
440
    {
441 442 443
    outPtr2 = outPtr + idxC;
    inPtr2 = inPtr + idxC;
    for (outIdx2 = outExt[4]; outIdx2 <= outExt[5]; ++outIdx2)
444
      {
445 446 447
      outPtr1 = outPtr2;
      inPtr1 = inPtr2;
      for (outIdx1 = outExt[2]; outIdx1 <= outExt[3]; ++outIdx1)
448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466
        {
        outPtr0 = outPtr1;
        inPtr0 = inPtr1;
        for (outIdx0 = outExt[0]; outIdx0 <= outExt[1]; ++outIdx0)
          {
          if (*outPtr0 == 3)
            {
            *outPtr0 = replaceValue;
            }
          else
            {
            *outPtr0 = *inPtr0;
            }
          inPtr0 += inInc0;
          outPtr0 += outInc0;
          }
        inPtr1 += inInc1;
        outPtr1 += outInc1;
        }
467 468
      inPtr2 += inInc2;
      outPtr2 += outInc2;
469 470 471 472
      }
    }
}

473

474
//----------------------------------------------------------------------------
475
// This method uses the input data to fill the output data.
476
// It can handle any type data, but the two datas must have the same
477
// data type.  Assumes that in and out have the same lower extent.
478
int vtkImageIslandRemoval2D::RequestData(
Amy Squillacote's avatar
Amy Squillacote committed
479 480 481
  vtkInformation *vtkNotUsed(request),
  vtkInformationVector **inputVector,
  vtkInformationVector *outputVector)
482
{
Amy Squillacote's avatar
Amy Squillacote committed
483
  int outExt[6];
484

Amy Squillacote's avatar
Amy Squillacote committed
485 486 487 488 489 490 491 492
  // get the data object
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkImageData *inData = vtkImageData::SafeDownCast(
    inInfo->Get(vtkDataObject::DATA_OBJECT()));
  vtkInformation *outInfo = outputVector->GetInformationObject(0);
  vtkImageData *outData = vtkImageData::SafeDownCast(
    outInfo->Get(vtkDataObject::DATA_OBJECT()));

493 494
  int wholeExtent[6];
  int extent[6];
495

496
  // We need to allocate our own scalars.
Amy Squillacote's avatar
Amy Squillacote committed
497 498
  outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), wholeExtent);
  outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), extent);
499 500 501 502 503
  extent[0] = wholeExtent[0];
  extent[1] = wholeExtent[1];
  extent[2] = wholeExtent[2];
  extent[3] = wholeExtent[3];
  outData->SetExtent(extent);
504
  outData->AllocateScalars(outInfo);
505

506
  // this filter expects that input is the same type as output.
507
  if (inData->GetScalarType() != outData->GetScalarType())
508
    {
509
    vtkErrorMacro(<< "Execute: input ScalarType, "
510
                  << vtkImageScalarTypeNameMacro(inData->GetScalarType())
Charles Law's avatar
Charles Law committed
511
                  << ", must match out ScalarType "
512
                  << vtkImageScalarTypeNameMacro(outData->GetScalarType()));
513
    return 1;
514 515
    }

Amy Squillacote's avatar
Amy Squillacote committed
516
  outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), outExt);
517 518
  void *inPtr = inData->GetScalarPointerForExtent(outExt);
  void *outPtr = outData->GetScalarPointerForExtent(outExt);
519

520
  switch (inData->GetScalarType())
521
    {
522
    vtkTemplateMacro(
523 524
      vtkImageIslandRemoval2DExecute(this, inData,
                                     static_cast<VTK_TT *>(inPtr), outData,
525
                                     static_cast<VTK_TT *>(outPtr), outExt));
526
    default:
Charles Law's avatar
Charles Law committed
527
      vtkErrorMacro(<< "Execute: Unknown ScalarType");
528 529 530 531
      return 1;
    }

  return 1;
532 533 534 535 536 537 538 539
}