PageRenderTime 4ms CodeModel.GetById 20ms app.highlight 68ms RepoModel.GetById 1ms app.codeStats 0ms

/Modules/Scripted/EditorLib/Logic/vtkImageConnectivity.cxx

https://github.com/lorensen/Slicer
C++ | 890 lines | 645 code | 76 blank | 169 comment | 136 complexity | 193ef064e7def273af510124cd1ccdad MD5 | raw file
  1/*=auto=========================================================================
  2
  3  Portions (c) Copyright 2005 Brigham and Women's Hospital (BWH) All Rights Reserved.
  4
  5  See COPYRIGHT.txt
  6  or http://www.slicer.org/copyright/copyright.txt for details.
  7
  8=========================================================================auto=*/
  9#include "vtkImageConnectivity.h"
 10
 11#include "vtkObjectFactory.h"
 12#include "vtkImageData.h"
 13#include <vtkInformation.h>
 14#include <vtkStreamingDemandDrivenPipeline.h>
 15
 16#include <limits.h>
 17#include <assert.h>
 18#include <stddef.h>
 19
 20//----------------------------------------------------------------------------
 21vtkStandardNewMacro(vtkImageConnectivity);
 22
 23//----------------------------------------------------------------------------
 24// Description:
 25// Constructor sets default values
 26vtkImageConnectivity::vtkImageConnectivity()
 27{
 28  this->Background = 0;
 29  this->MinForeground = VTK_SHORT_MIN;
 30  this->MaxForeground = VTK_SHORT_MAX;
 31  this->MinSize = 10000;
 32  this->Function = CONNECTIVITY_MEASURE;
 33  this->OutputLabel = 1;
 34  this->SliceBySlice = 0;
 35  this->LargestIslandSize = this->IslandSize = 0;
 36  this->Seed[0] = this->Seed[1] = this->Seed[2] = 0;
 37}
 38
 39//----------------------------------------------------------------------------
 40const char* vtkImageConnectivity::GetFunctionString()
 41{
 42  switch (this->Function)
 43    {
 44    case CONNECTIVITY_IDENTIFY:
 45      return "IdentifyIslands";
 46    case CONNECTIVITY_REMOVE:
 47      return "RemoveIslands";
 48    case CONNECTIVITY_CHANGE:
 49      return "ChangeIsland";
 50    case CONNECTIVITY_MEASURE:
 51      return "MeasureIsland";
 52    case CONNECTIVITY_SAVE:
 53      return "SaveIsland";
 54    default:
 55      return "ERROR: Unknown";
 56    }
 57}
 58
 59//************************************************************************
 60// The following procedures were written by Andre Robatino
 61// in November, 1999
 62//
 63// connect
 64// recursive_copy
 65//************************************************************************
 66
 67int connect(unsigned int, size_t *, char *, char, size_t *, size_t *);
 68static void recursive_copy(int, size_t);
 69
 70static size_t *g_axis_len, *g_outimagep, num_stride_index_bits;
 71static char *g_inimagep, g_inbackground;
 72
 73int connect(
 74     unsigned int rank,
 75     size_t *axis_len,
 76     char *inimage,
 77     char inbackground,
 78     size_t *outimage,
 79     size_t *num_components) /* set to NULL if not interested */
 80{
 81  unsigned int i;
 82  register unsigned int axisv;
 83  size_t data_len;
 84  register size_t label, *outimagep, *outimage_end, *imagep, *new_imagep, boundary_mask_start, component_mask, axis_mask;
 85  register ptrdiff_t *stride, stridev;
 86
 87  //assert(rank >= 0); // Note: Always true !
 88  if (rank == 0) {
 89    *outimage = (*inimage != inbackground);
 90    if (num_components) *num_components = *outimage;
 91    return 0;
 92  }
 93  label = 2*rank;
 94  num_stride_index_bits = 1;
 95  while (label >>= 1) num_stride_index_bits++;
 96  assert(num_stride_index_bits + 2*rank + 1 <= CHAR_BIT*sizeof(size_t));
 97  boundary_mask_start = 1<<num_stride_index_bits;
 98  axis_mask = boundary_mask_start - 1;
 99  stride = (ptrdiff_t *)malloc((2*rank + 1)*sizeof(ptrdiff_t));
100  assert(stride);
101  data_len = 1;
102  for (i=0; i<rank; i++) {
103    assert(axis_len[i] > 1);
104    stride[2*i + 1] = -(stride[2*i] = data_len);
105    data_len *= axis_len[i];
106  }
107  stride[2*rank] = 0;
108  g_axis_len = axis_len;
109  g_inimagep = inimage;
110  g_inbackground = inbackground;
111  g_outimagep = outimage;
112  component_mask = 0;
113  component_mask = ~(~component_mask>>1);
114  recursive_copy(rank-1, component_mask);
115  outimagep = outimage;
116  outimage_end = outimage + data_len;
117  label = 0;
118  do {
119    if (!(*outimagep & component_mask)) continue;
120    imagep = outimagep;
121    *imagep ^= component_mask;
122    axisv = 0;
123    label++;
124    while (1) {
125      while ( (stridev = stride[axisv]) ) {
126        if ((*imagep & boundary_mask_start<<axisv) && (*(new_imagep = imagep + stridev) & component_mask)) {
127          imagep = new_imagep;
128          *imagep ^= component_mask;
129          *imagep |= axisv;
130          axisv = 0;
131          continue;
132        }
133        axisv++;
134      }
135      if (imagep == outimagep) break;
136      axisv = *imagep & axis_mask;
137      *imagep = label;
138      imagep -= stride[axisv++];
139    }
140    *outimagep = label;
141  } while (++outimagep < outimage_end);
142  if (num_components) *num_components = label;
143  free(stride);
144  return 0;
145}
146
147static void recursive_copy(
148     int axis,
149     register size_t mask)
150{
151  register size_t len = g_axis_len[axis] - 2;
152
153  if (axis == 0) {
154    register size_t *outimagep;
155    register char *inimagep, inbackground;
156
157    inimagep = g_inimagep;
158    inbackground = g_inbackground;
159    outimagep = g_outimagep;
160    mask |= 1<<num_stride_index_bits;
161    *outimagep++ = (*inimagep++ == inbackground)? 0 : mask;
162    mask |= 2<<num_stride_index_bits;
163    while (len--) *outimagep++ = (*inimagep++ == inbackground)? 0 : mask;
164    mask ^= 1<<num_stride_index_bits;
165    *outimagep++ = (*inimagep++ == inbackground)? 0 : mask;
166    g_outimagep = outimagep;
167    g_inimagep = inimagep;
168  } else {
169    mask |= 1<<(num_stride_index_bits + 2*axis);
170    recursive_copy(axis-1, mask);
171    mask |= 2<<(num_stride_index_bits + 2*axis);
172    while (len--) recursive_copy(axis-1, mask);
173    mask ^= 1<<(num_stride_index_bits + 2*axis);
174    recursive_copy(axis-1, mask);
175  }
176  return;
177}
178
179//***********************************************************************
180// End Andre's cool routines.
181//************************************************************************
182
183
184static void vtkImageConnectivityExecute(vtkImageConnectivity *self,
185                     vtkImageData *inData, short *inPtr,
186                     vtkImageData *outData, short *outPtr,
187                     int outExt[6])
188{
189  // For looping though output (and input) pixels.
190  int outMin0, outMax0, outMin1, outMax1, outMin2, outMax2;
191  int outIdx0, outIdx1, outIdx2;
192  vtkIdType inInc0, inInc1, inInc2;
193  vtkIdType outInc0, outInc1, outInc2;
194  short *inPtr0, *outPtr0, *outPtr1;
195  short minForegnd = (short)self->GetMinForeground();
196  short maxForegnd = (short)self->GetMaxForeground();
197  short newLabel = (short)self->GetOutputLabel();
198  short seedLabel = 0;
199  int largest, len=1, nxy, z, nz = 0, j;
200  int *census = NULL;
201  int seed[3];
202  int minSize = self->GetMinSize();
203  short pix;
204  int identifyIslands = self->GetFunction() == CONNECTIVITY_IDENTIFY;
205  int removeIslands   = self->GetFunction() == CONNECTIVITY_REMOVE;
206  int changeIsland    = self->GetFunction() == CONNECTIVITY_CHANGE;
207  int saveIsland      = self->GetFunction() == CONNECTIVITY_SAVE;
208  int measureIsland   = self->GetFunction() == CONNECTIVITY_MEASURE;
209  int sliceBySlice    = self->GetSliceBySlice();
210
211  // connect
212  size_t conSeedLabel = 0, i, idx, dz;
213  int rank;
214  size_t *axis_len=NULL;
215  unsigned short bg = self->GetBackground();
216  short bgMask = 0;
217  short fgMask = 1;
218  char inbackground = (char)bgMask;
219  char *conInput=NULL;
220  size_t *conOutput=NULL;
221  size_t *numIslands=NULL;
222
223  // Image bounds
224  outMin0 = outExt[0];   outMax0 = outExt[1];
225  outMin1 = outExt[2];   outMax1 = outExt[3];
226  outMin2 = outExt[4];   outMax2 = outExt[5];
227
228  // Computer Parameters for connect().
229  rank = (outExt[5]==outExt[4]) ? 2 : 3;
230  axis_len = new size_t[rank+1];
231  axis_len[0] = outExt[1]-outExt[0]+1;
232  axis_len[1] = outExt[3]-outExt[2]+1;
233  axis_len[2] = outExt[5]-outExt[4]+1;
234  for (j=0; j<rank; j++)
235  {
236    len *= axis_len[j];
237  }
238  conInput = new char[len];
239  conOutput = new size_t[len];
240  numIslands = new size_t[axis_len[2]];
241
242  // Get increments to march through data continuously
243  outData->GetContinuousIncrements(outExt, outInc0, outInc1, outInc2);
244  inData->GetContinuousIncrements(outExt, inInc0, inInc1, inInc2);
245
246  ///////////////////////////////////////////////////////////////
247  // Save, Change, Measure:
248  // ----------------------
249  // Get the seed
250  //
251  //   seedLabel = inData[xSeed,ySeed,zSeed]
252  //
253  // If the seed is out of bounds, return the input
254  //
255  //   outData[i] = inData[i]
256  //
257  ///////////////////////////////////////////////////////////////
258
259  if (changeIsland || measureIsland || saveIsland)
260    {
261    self->GetSeed(seed);
262
263    if (seed[0] < outMin0 || seed[0] > outMax0 ||
264        seed[1] < outMin1 || seed[1] > outMax1 ||
265        seed[2] < outMin2 || seed[2] > outMax2)
266      {
267      //
268      // Out of bounds -- abort!
269      //
270      inPtr0 = inPtr;
271      outPtr0 = outPtr;
272      for (outIdx2 = outMin2; outIdx2 <= outMax2; outIdx2++)
273        {
274        for (outIdx1 = outMin1; outIdx1 <= outMax1; outIdx1++)
275          {
276          for (outIdx0 = outMin0; outIdx0 <= outMax0; outIdx0++)
277            {
278            *outPtr0 = *inPtr0;
279            outPtr0++;
280            inPtr0++;
281            }//for0
282          outPtr0 += outInc1;
283          inPtr0 += inInc1;
284          }//for1
285        outPtr0 += outInc2;
286        inPtr0 += inInc2;
287        }//for2
288
289      fprintf(stderr, "Seed %d,%d,%d out of bounds in CCA.\n",
290        seed[0], seed[1], seed[2]);
291      return;
292      }
293
294    //
295    // In bounds!
296    //
297    outPtr1 = (short*)inData->GetScalarPointer(seed[0], seed[1], seed[2]);
298    seedLabel = *outPtr1;
299    }
300
301  ///////////////////////////////////////////////////////////////
302  // Remove, Identify:
303  // ----------------------
304  // Create a mask where everything outside [min,max] or in the
305  // sea (bg) is in the background.
306  //
307  //   conInput[i] = fgMask,   inData[i] != bg
308  //               = bgMask,   else
309  //
310  //   conInput[i] = unchanged
311  //               = bgMask,   where inData[i] not on [min,max]
312  //
313  ///////////////////////////////////////////////////////////////
314
315  if (removeIslands || identifyIslands)
316    {
317    inPtr0 = inPtr;
318    i = 0;
319    for (outIdx2 = outMin2; outIdx2 <= outMax2; outIdx2++)
320      {
321      for (outIdx1 = outMin1; outIdx1 <= outMax1; outIdx1++)
322        {
323        for (outIdx0 = outMin0; outIdx0 <= outMax0; outIdx0++)
324          {
325          if (*inPtr0 != bg)
326            {
327            conInput[i] = fgMask;
328            }
329          else
330            {
331            conInput[i] = bgMask;
332            }
333          inPtr0++;
334          i++;
335          }//for0
336        inPtr0 += inInc1;
337        }//for1
338      inPtr0 += inInc2;
339      }//for2
340
341    // Optionally threshold [min,max]
342    if(minForegnd > VTK_SHORT_MIN || maxForegnd < VTK_SHORT_MAX)
343      {
344      inPtr0 = inPtr;
345      i=0;
346      for (outIdx2 = outMin2; outIdx2 <= outMax2; outIdx2++)
347        {
348        for (outIdx1 = outMin1; outIdx1 <= outMax1; outIdx1++)
349          {
350          for (outIdx0 = outMin0; outIdx0 <= outMax0; outIdx0++)
351            {
352            pix = *inPtr0;
353            if (pix < minForegnd || pix > maxForegnd)
354              {
355              conInput[i] = bgMask;
356              }
357            i++;
358            inPtr0++;
359            }//for0
360          inPtr0 += inInc1;
361          }//for1
362        inPtr0 += inInc2;
363        }//for2
364      }
365    }
366
367
368  ///////////////////////////////////////////////////////////////
369  // Save, Change, Measure:
370  // ----------------------
371  //
372  // Create a mask where everything not equal to seedLabel is
373  // in the background.
374  //
375  //     conInput[i] = fgMask,  inData[i] == seedLabel
376  //                 = bgMask,  else
377  //
378  ///////////////////////////////////////////////////////////////
379
380  if (saveIsland || changeIsland || measureIsland)
381    {
382    inPtr0 = inPtr;
383    i = 0;
384    for (outIdx2 = outMin2; outIdx2 <= outMax2; outIdx2++)
385      {
386      for (outIdx1 = outMin1; outIdx1 <= outMax1; outIdx1++)
387        {
388        for (outIdx0 = outMin0; outIdx0 <= outMax0; outIdx0++)
389          {
390          if (*inPtr0 == seedLabel)
391            {
392            conInput[i] = fgMask;
393            }
394          else
395            {
396            conInput[i] = bgMask;
397            }
398          inPtr0++;
399          i++;
400          }//for0
401        inPtr0 += inInc1;
402        }//for1
403      inPtr0 += inInc2;
404      }//for2
405    }
406
407
408  ///////////////////////////////////////////////////////////////
409  // Save, Change, Measure, Remove, Identify
410  // ---------------------------------------
411  // Run Connectivity
412  //
413  ///////////////////////////////////////////////////////////////
414
415  if (saveIsland || changeIsland || measureIsland || removeIslands || identifyIslands)
416    {
417    nz = 1;
418    if (sliceBySlice && removeIslands)
419      {
420      // If SliceBySlice, then call connect() for each slice
421      nxy = axis_len[0] * axis_len[1];
422      nz = axis_len[2];
423      rank = 2;
424      int axis_len2 = axis_len[2];
425      axis_len[2] = 1;
426
427      for (z=0; z < nz; z++)
428        {
429        connect(rank, axis_len, &conInput[nxy*z], inbackground,
430          &conOutput[nxy*z], &numIslands[z]);
431        }
432      axis_len[2] = axis_len2;
433      }
434    else
435      {
436      connect(rank, axis_len, conInput, inbackground, conOutput, &numIslands[0]);
437      }
438    }
439
440
441  ///////////////////////////////////////////////////////////////
442  // Save, Change, Measure
443  // -----------------------------
444  // Get conSeedLabel
445  //
446  //   conSeedLabel = conOutput[xSeed,ySeed,zSeed]
447  //
448  ///////////////////////////////////////////////////////////////
449
450  if (saveIsland || changeIsland || measureIsland)
451    {
452    i = seed[2]*axis_len[1]*axis_len[0] + seed[1]*axis_len[0] + seed[0];
453    conSeedLabel = conOutput[i];
454    }
455
456
457  ///////////////////////////////////////////////////////////////
458  // Measure, Remove
459  // -----------------------------
460  // Count size of each island in conOutput
461  //
462  //   census[c] = COUNT(conOutput[c]),  forall c on [0,numIslands]
463  //
464  ///////////////////////////////////////////////////////////////
465
466  if (removeIslands || measureIsland)
467    {
468    // For each label value, count the number of pixels with that label
469    // If SliceBySlice, then work on each slice one at a time
470    len = 0;
471    for (z=0; z<nz; z++)
472      {
473      len += numIslands[z] + 1;
474      }
475    census = new int[len];
476    memset(census, 0, len*sizeof(int));
477
478    if (nz == 1)
479      {
480      i = 0;
481      for (outIdx2 = outMin2; outIdx2 <= outMax2; outIdx2++)
482        {
483        for (outIdx1 = outMin1; outIdx1 <= outMax1; outIdx1++)
484          {
485          for (outIdx0 = outMin0; outIdx0 <= outMax0; outIdx0++)
486            {
487            idx = conOutput[i];
488            // Note: Since 'idx' is defined as 'size_t' it's of type 'unsigned'
489            // It means 'idx >= 0' is always true.
490            if (/*idx >= 0 && */idx <= numIslands[0])
491              {
492              census[idx] = census[idx] + 1;
493              }
494            i++;
495            }//for0
496          }//for1
497        }//for2
498      }
499    else
500      {
501      dz = 0;
502      i = 0;
503      for (z=0; z < nz; z++)
504        {
505        for (outIdx1 = outMin1; outIdx1 <= outMax1; outIdx1++)
506          {
507          for (outIdx0 = outMin0; outIdx0 <= outMax0; outIdx0++)
508            {
509            idx = conOutput[i];
510            // Note: Since 'idx' is defined as 'size_t' it's of type 'unsigned'
511            // It means 'idx >= 0' is always true.
512            if (/*idx >= 0 && */idx <= numIslands[z])
513              {
514              census[dz+idx] = census[dz+idx] + 1;
515              }
516            i++;
517            }//for0
518          }//for1
519
520        dz += numIslands[z]+1;
521        }//forz
522      }
523    }
524
525
526  ///////////////////////////////////////////////////////////////
527  // Remove
528  // -----------------------------
529  // Output gets input except where islands too small
530  //
531  //   outData[i] = inData[i],  census[conOutput[i]] >= minIslandSize
532  //              = bg,    else
533  //
534  ///////////////////////////////////////////////////////////////
535
536  if (removeIslands)
537    {
538    if (nz == 1)
539      {
540      inPtr0 = inPtr;
541      outPtr0 = outPtr;
542      i = 0;
543      for (outIdx2 = outMin2; outIdx2 <= outMax2; outIdx2++)
544        {
545        for (outIdx1 = outMin1; outIdx1 <= outMax1; outIdx1++)
546          {
547          for (outIdx0 = outMin0; outIdx0 <= outMax0; outIdx0++)
548            {
549            if (census[conOutput[i]] >= minSize)
550              {
551              *outPtr0 = *inPtr0;
552              }
553            else
554              {
555              *outPtr0 = bg;
556              }
557            i++;
558            outPtr0++;
559            inPtr0++;
560            }//for0
561          outPtr0 += outInc1;
562          inPtr0 += inInc1;
563          }//for1
564        outPtr0 += outInc2;
565        inPtr0 += inInc2;
566        }//for2
567      }
568    else
569      {
570      dz = 0;
571      i = 0;
572      inPtr0 = inPtr;
573      outPtr0 = outPtr;
574      for (z=0; z < nz; z++)
575        {
576        for (outIdx1 = outMin1; outIdx1 <= outMax1; outIdx1++)
577          {
578          for (outIdx0 = outMin0; outIdx0 <= outMax0; outIdx0++)
579            {
580            if (census[dz+conOutput[i]] >= minSize)
581              {
582              *outPtr0 = *inPtr0;
583              }
584            else
585              {
586              *outPtr0 = bg;
587              }
588            i++;
589            outPtr0++;
590            inPtr0++;
591            }//for0
592          outPtr0 += outInc1;
593          inPtr0 += inInc1;
594          }//for1
595        outPtr0 += outInc2;
596        inPtr0 += inInc2;
597
598        dz += numIslands[z] + 1;
599        }//z
600      }//else
601    }
602
603
604  ///////////////////////////////////////////////////////////////
605  // Measure
606  // -----------------------------
607  // Store statistics, and return output = input.
608  //
609  //   islandSize = census[conSeedLabel]
610  //   largest    = MAX(census[c])
611  //   outData[i] = inData[i]
612  //
613  ///////////////////////////////////////////////////////////////
614
615  if (measureIsland)
616    {
617    // Find largest island
618    largest = 0;
619    for (i=0; i<=numIslands[0]; i++)
620      {
621      if (i != bg)
622        {
623        if (census[i] > largest)
624          {
625          largest = census[i];
626          }
627        }
628      }
629    self->SetLargestIslandSize(largest);
630
631    // Measure island at seed
632    self->SetIslandSize(census[conSeedLabel]);
633
634    // Return output values to be the inputs
635    inPtr0 = inPtr;
636    outPtr0 = outPtr;
637    for (outIdx2 = outMin2; outIdx2 <= outMax2; outIdx2++)
638      {
639      for (outIdx1 = outMin1; outIdx1 <= outMax1; outIdx1++)
640        {
641        for (outIdx0 = outMin0; outIdx0 <= outMax0; outIdx0++)
642          {
643          *outPtr0 = *inPtr0;
644          outPtr0++;
645          inPtr0++;
646          }//for0
647        outPtr0 += outInc1;
648        inPtr0 += inInc1;
649        }//for1
650      outPtr0 += outInc2;
651      inPtr0 += inInc2;
652      }//for2
653    }
654
655
656  ///////////////////////////////////////////////////////////////
657  // Identify
658  // -----------------------------
659  // Output gets the output of connect()
660  //
661  //   outData[i] = conOutput[i]
662  //
663  ///////////////////////////////////////////////////////////////
664
665  if (identifyIslands)
666    {
667    outPtr0 = outPtr;
668    i = 0;
669    for (outIdx2 = outMin2; outIdx2 <= outMax2; outIdx2++)
670      {
671      for (outIdx1 = outMin1; outIdx1 <= outMax1; outIdx1++)
672        {
673        for (outIdx0 = outMin0; outIdx0 <= outMax0; outIdx0++)
674          {
675          *outPtr0 = (short)conOutput[i];
676          i++;
677          outPtr0++;
678          }//for0
679        outPtr0 += outInc1;
680        }//for1
681      outPtr0 += outInc2;
682      }//for2
683    }
684
685  ///////////////////////////////////////////////////////////////
686  // Remove, Identify
687  // -----------------------------
688  // Output gets input where the input was thresholded away
689  //
690  //   outData[i] = inData[i],  inData[i] outside [min,max]
691  //              = do nothing, else
692  //
693  ///////////////////////////////////////////////////////////////
694
695  if (removeIslands || identifyIslands)
696    {
697    if(minForegnd > VTK_SHORT_MIN || maxForegnd < VTK_SHORT_MAX)
698      {
699      inPtr0 = inPtr;
700      outPtr0 = outPtr;
701      for (outIdx2 = outMin2; outIdx2 <= outMax2; outIdx2++)
702        {
703        for (outIdx1 = outMin1; outIdx1 <= outMax1; outIdx1++)
704          {
705          for (outIdx0 = outMin0; outIdx0 <= outMax0; outIdx0++)
706            {
707            pix = *inPtr0;
708            if (pix < minForegnd || pix > maxForegnd)
709              {
710              *outPtr0 = pix;
711              }
712            inPtr0++;
713            outPtr0++;
714            }//for0
715          inPtr0 += inInc1;
716          outPtr0 += outInc1;
717          }//for1
718        inPtr0 += inInc2;
719        outPtr0 += outInc2;
720        }//for2
721      }
722    }
723
724  if (removeIslands || measureIsland)
725    {
726    delete [] census;
727    }
728
729
730  ///////////////////////////////////////////////////////////////
731  // Save
732  // -----------------------------
733  // Output gets input where seedLabel, else bg
734  //
735  //   outData[i] = inData[i],  conOutput[i] == conSeedLabel
736  //              = bg,         else
737  //
738  ///////////////////////////////////////////////////////////////
739
740  if (saveIsland)
741    {
742    inPtr0 = inPtr;
743    outPtr0 = outPtr;
744    i = 0;
745    for (outIdx2 = outMin2; outIdx2 <= outMax2; outIdx2++)
746      {
747      for (outIdx1 = outMin1; outIdx1 <= outMax1; outIdx1++)
748        {
749        for (outIdx0 = outMin0; outIdx0 <= outMax0; outIdx0++)
750          {
751          if (conOutput[i] == conSeedLabel)
752            {
753            *outPtr0 = *inPtr0;
754            }
755          else
756            {
757            *outPtr0 = bg;
758            }
759          i++;
760          outPtr0++;
761          inPtr0++;
762          }//for0
763        outPtr0 += outInc1;
764        inPtr0 += inInc1;
765        }//for1
766      outPtr0 += outInc2;
767      inPtr0 += inInc2;
768      }//for2
769    }
770
771
772  ///////////////////////////////////////////////////////////////
773  // Change
774  // -----------------------------
775  // Output gets newLabel where seedLabel, else input
776  //
777  //   outData[i] = newLabel,   conOutput[i] == conSeedLabel
778  //              = inData[i],  else
779  //
780  ///////////////////////////////////////////////////////////////
781
782  if (changeIsland)
783    {
784    inPtr0 = inPtr;
785    outPtr0 = outPtr;
786    i = 0;
787    for (outIdx2 = outMin2; outIdx2 <= outMax2; outIdx2++)
788      {
789      for (outIdx1 = outMin1; outIdx1 <= outMax1; outIdx1++)
790        {
791        for (outIdx0 = outMin0; outIdx0 <= outMax0; outIdx0++)
792          {
793          if (conOutput[i] == conSeedLabel)
794            {
795            *outPtr0 = newLabel;
796            }
797          else
798            {
799            *outPtr0 = *inPtr0;
800            }
801          i++;
802          outPtr0++;
803          inPtr0++;
804          }//for0
805        outPtr0 += outInc1;
806        inPtr0 += inInc1;
807        }//for1
808      outPtr0 += outInc2;
809      inPtr0 += inInc2;
810      }//for2
811    }
812
813  ///////////////////////////////////////////////////////////////
814  // Cleanup
815  ///////////////////////////////////////////////////////////////
816
817  delete [] axis_len;
818  delete [] numIslands;
819  delete [] conInput;
820  delete [] conOutput;
821}
822
823
824//----------------------------------------------------------------------------
825// Description:
826// This method is passed a input and output data, and executes the filter
827// algorithm to fill the output from the input.
828// It just executes a switch statement to call the correct function for
829// the datas data types.
830#if (VTK_MAJOR_VERSION <= 5)
831void vtkImageConnectivity::ExecuteData(vtkDataObject *)
832{
833  vtkImageData *inData = this->GetImageDataInput(0);
834  vtkImageData *outData = this->GetOutput();
835  outData->SetExtent(outData->GetWholeExtent());
836  outData->AllocateScalars();
837
838  int outExt[6], s;
839  outData->GetWholeExtent(outExt);
840#else
841void vtkImageConnectivity::ExecuteDataWithInformation(vtkDataObject *output, vtkInformation* outInfo)
842{
843  vtkImageData *inData = vtkImageData::SafeDownCast(this->GetInput());
844  vtkImageData *outData = this->AllocateOutputData(output, outInfo);
845
846  int outExt[6], s;
847  outInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), outExt);
848#endif
849  void *inPtr = inData->GetScalarPointerForExtent(outExt);
850  void *outPtr = outData->GetScalarPointerForExtent(outExt);
851
852  int x1;
853
854  x1 = inData->GetNumberOfScalarComponents();
855  if (x1 != 1)
856    {
857    vtkErrorMacro(<<"Input has "<<x1<<" instead of 1 scalar component.");
858    return;
859    }
860
861  /* Need short data */
862  s = inData->GetScalarType();
863  if (s != VTK_SHORT)
864    {
865    vtkErrorMacro("Warning: Input scalars are type "<<s
866      << " instead of "<<VTK_SHORT);
867    return;
868    }
869
870  vtkImageConnectivityExecute(this, inData, (short *)inPtr,
871    outData, (short *)(outPtr), outExt);
872}
873
874//----------------------------------------------------------------------------
875void vtkImageConnectivity::PrintSelf(ostream& os, vtkIndent indent)
876{
877  this->Superclass::PrintSelf(os,indent);
878
879  os << indent << "Background:        " << this->Background << "\n";
880  os << indent << "MinForeground:     " << this->MinForeground << "\n";
881  os << indent << "MaxForeground:     " << this->MaxForeground << "\n";
882  os << indent << "LargestIslandSize: " << this->LargestIslandSize << "\n";
883  os << indent << "IslandSize:        " << this->IslandSize << "\n";
884  os << indent << "MinSize:           " << this->MinSize << "\n";
885  os << indent << "OutputLabel:       " << this->OutputLabel << "\n";
886  os << indent << "Seed[0]:           " << this->Seed[0] << "\n";
887  os << indent << "Seed[1]:           " << this->Seed[1] << "\n";
888  os << indent << "Seed[2]:           " << this->Seed[2] << "\n";
889  os << indent << "Function:          " << this->Function << "\n";
890}