vtkVolumeRayCastSpaceLeapingImageFilter.cxx 37.1 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkVolumeRayCastSpaceLeapingImageFilter.cxx

  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
  All rights reserved.
  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.

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

=========================================================================*/
#include "vtkVolumeRayCastSpaceLeapingImageFilter.h"

#include "vtkImageData.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkObjectFactory.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkDataArray.h"
23
#include "vtkPointData.h"
24 25 26
#include <sstream>
#include <fstream>
#include <iostream>
27 28 29 30

#ifdef vtkVolumeRayCastSpaceLeapingImageFilter_DEBUG
#include "vtkMetaImageWriter.h"
#endif
31

32
#include <cmath>
33 34 35 36 37 38 39 40 41 42 43 44 45 46

// Space leaping block size
#define VTK_SL_BLK 4

//----------------------------------------------------------------------------
vtkStandardNewMacro(vtkVolumeRayCastSpaceLeapingImageFilter);
vtkCxxSetObjectMacro(vtkVolumeRayCastSpaceLeapingImageFilter,
                     CurrentScalars, vtkDataArray);

//----------------------------------------------------------------------------
vtkVolumeRayCastSpaceLeapingImageFilter::vtkVolumeRayCastSpaceLeapingImageFilter()
{
  this->ComputeMinMax = 0;
  this->ComputeGradientOpacity = 0;
47
  this->UpdateGradientOpacityFlags = 0;
48
  this->IndependentComponents  = 1;
49 50 51 52
  this->CurrentScalars = nullptr;
  this->MinNonZeroScalarIndex = nullptr;
  this->MinNonZeroGradientMagnitudeIndex = nullptr;
  this->GradientMagnitude = nullptr;
53
  for (int i = 0; i < 4; i++)
54
  {
55 56 57
    this->TableSize[i] = 0;
    this->TableShift[i] = 0;
    this->TableScale[i] = 1;
58 59
    this->ScalarOpacityTable[i] = nullptr;
    this->GradientOpacityTable[i] = nullptr;
60
  }
61
  this->Cache = nullptr;
62 63 64 65 66

  // Ensure that no splits occur along X or Y axes when multithreading,
  // there seems to be a bug with the increments that requires this.
  this->SplitPath[0] = 2;
  this->SplitPathLength = 1;
67 68 69 70 71
}

//----------------------------------------------------------------------------
vtkVolumeRayCastSpaceLeapingImageFilter::~vtkVolumeRayCastSpaceLeapingImageFilter()
{
72
  this->SetCurrentScalars(nullptr);
73 74
  delete [] this->MinNonZeroScalarIndex;
  delete [] this->MinNonZeroGradientMagnitudeIndex;
75 76 77 78 79 80 81 82 83
}

//----------------------------------------------------------------------------
void vtkVolumeRayCastSpaceLeapingImageFilter::PrintSelf(ostream& os, vtkIndent indent)
{
  this->Superclass::PrintSelf(os, indent);

  os << indent << "ComputeMinMax: " << this->ComputeMinMax << "\n";
  os << indent << "ComputeGradientOpacity: " << this->ComputeGradientOpacity << "\n";
84
  os << indent << "UpdateGradientOpacityFlags: " << this->UpdateGradientOpacityFlags << "\n";
85 86
  os << indent << "IndependentComponents: " << this->IndependentComponents << "\n";
  os << indent << "CurrentScalars: " << this->CurrentScalars << "\n";
87 88 89 90 91 92 93
  // this->TableShift
  // this->TableScale
  // this->TableSize
  // this->ScalarOpacityTable
  // this->GradientOpacityTable
  // this->MinNonZeroScalarIndex
  // this->MinNonZeroGradientMagnitudeIndex
94 95 96 97 98 99
}

//----------------------------------------------------------------------------
int vtkVolumeRayCastSpaceLeapingImageFilter::RequestUpdateExtent (
  vtkInformation * vtkNotUsed(request),
  vtkInformationVector **inputVector,
100
  vtkInformationVector *vtkNotUsed(outputVector))
101 102 103 104 105 106 107 108 109 110 111 112 113
{
  // get the info objects
  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);

  // Ask for the whole input

  int wholeExtent[6];
  inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), wholeExtent);
  inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), wholeExtent, 6);

  return 1;
}

114 115 116 117 118 119 120 121
//----------------------------------------------------------------------------
void vtkVolumeRayCastSpaceLeapingImageFilter
::SetCache( vtkImageData * cache )
{
  // Do not reference count it to avoid reference counting loops
  this->Cache = cache;
}

122 123 124 125 126 127 128 129 130
//----------------------------------------------------------------------------
void vtkVolumeRayCastSpaceLeapingImageFilter
::InternalRequestUpdateExtent( int *inExt,
                               int *wholeExtent)
{
  int dim[3];

  // We group four cells (which require 5 samples) into one element in the min/max tree
  for ( int i = 0; i < 3; i++ )
131
  {
132 133 134 135 136 137
    // size of the input image.
    dim[i] = wholeExtent[2*i+1] - wholeExtent[2*i] + 1;

    inExt[2*i] = 0; // The output extent is 0 based.
    inExt[2*i+1] = (dim[i] < 2) ? (0) :
        (static_cast<int>((dim[i] - 2)/VTK_SL_BLK));
138
  }
139 140 141
}

//----------------------------------------------------------------------------
142
static void
143
vtkVolumeRayCastSpaceLeapingImageFilterClearOutput(vtkImageData *outData,
144
                                                   int outExt[6],
145
                                                   int nComponents )
146 147 148 149 150 151 152
{
  unsigned short *tmpPtr = static_cast< unsigned short * >(
                outData->GetScalarPointerForExtent(outExt));

  // Get increments to march through the thread's output extents

  vtkIdType outInc0, outInc1, outInc2;
153
  outData->GetContinuousIncrements(outExt, outInc0, outInc1, outInc2);
154 155 156
  // A. Initialize the arrays with a blank flag.

  int i,j,k;
157
  int c;
158
  for (k = outExt[4]; k <= outExt[5]; ++k, tmpPtr += outInc2)
159
  {
160
    for (j = outExt[2]; j <= outExt[3]; ++j, tmpPtr += outInc1)
161
    {
162
      for (i = outExt[0]; i <= outExt[1]; ++i)
163
      {
164
        for ( c = 0; c < nComponents; ++c )
165
        {
166 167 168
          *(tmpPtr++) = 0xffff;  // Min Scalar
          *(tmpPtr++) = 0;       // Max Scalar
          *(tmpPtr++) = 0;       // Max Gradient Magnitude and
169
        }                      // Flag computed from transfer functions
170 171
      }
    }
172
  }
173 174 175 176
}

//----------------------------------------------------------------------------
void vtkVolumeRayCastSpaceLeapingImageFilter
177
::ComputeInputExtentsForOutput( int inExt[6], int inDim[3],
178 179 180 181 182
                                int outExt[6], vtkImageData *inData )
{
  int inWholeExt[6];
  inData->GetExtent(inWholeExt);

183
  for ( int i = 0; i < 3; i++ )
184
  {
185 186 187 188 189 190 191 192
    inExt[2*i] = outExt[2*i] * VTK_SL_BLK + inWholeExt[2*i];

    // Extra +1 needed here since we group four cells (which require 5
    // samples) into one element in the min/max tree
    inExt[2*i+1] = (outExt[2*i+1]+1) * VTK_SL_BLK + inWholeExt[2*i] + 1;

    // Clip the extents with the whole extent.
    if (inExt[2*i] < inWholeExt[2*i])
193
    {
194
      inExt[2*i] = inWholeExt[2*i];
195
    }
196
    if (inExt[2*i+1] > inWholeExt[2*i+1])
197
    {
198
      inExt[2*i+1] = inWholeExt[2*i+1];
199
    }
200 201

    inDim[i] = inExt[2*i+1] - inExt[2*i] + 1;
202
  }
203 204 205 206 207 208 209 210 211 212 213 214
}

//----------------------------------------------------------------------------
// Fill in the min-max space leaping information.
template <class T>
void
vtkVolumeRayCastSpaceLeapingImageFilterMinMaxExecute(
    vtkVolumeRayCastSpaceLeapingImageFilter *self,
    vtkImageData *inData,
    vtkImageData *outData, int outExt[6],
    T )
{
215

216 217 218 219 220
  // the number of independent components for which we need to keep track of
  // min/max
  vtkDataArray * scalars = self->GetCurrentScalars();
  const int components = scalars->GetNumberOfComponents();
  const int independent = self->GetIndependentComponents();
221
  const int nComponents = (independent) ? components : 1;
222 223 224 225 226 227

  // B. Now fill in the max-min-gradient volume structure

  // B.1 First compute the extents of the input that contribute to this structure

  int inExt[6], inWholeExt[6];
228
  int inDim[3];
229 230 231 232 233 234 235 236 237 238 239 240 241 242
  int outWholeDim[3];
  vtkVolumeRayCastSpaceLeapingImageFilter::ComputeInputExtentsForOutput(
    inExt, inDim, outExt, inData );
  inData->GetExtent(inWholeExt);
  outData->GetDimensions(outWholeDim);

  float shift[4], scale[4];
  self->GetTableShift(shift);
  self->GetTableScale(scale);


  // B.2 Get increments to march through the input extents

  vtkIdType inInc0, inInc1, inInc2;
243 244
  inData->GetContinuousIncrements(scalars,
                                  inExt, inInc0, inInc1, inInc2);
245 246 247

  // Get increments to march through the output extents

248 249 250
  const vtkIdType outInc0 = 3*nComponents;
  const vtkIdType outInc1 = outInc0*outWholeDim[0];
  const vtkIdType outInc2 = outInc1*outWholeDim[1];
251 252 253

  // B.3 Now fill in the min-max volume.

254 255 256 257
  int i, j, k;
  int c;
  int sx1, sx2, sy1, sy2, sz1, sz2;
  int x, y, z;
258 259 260 261 262 263 264 265 266 267 268 269 270

  T *dptr = static_cast< T * >(scalars->GetVoidPointer(0));
  unsigned short val;
  unsigned short *outBasePtr = static_cast< unsigned short * >(
                                outData->GetScalarPointer());

  // Initialize pointer to the starting extents given by inExt.
  dptr += self->ComputeOffset( inExt, inWholeExt, nComponents );

  // The pointer into the space-leaping output volume.
  unsigned short *tmpPtr, *tmpPtrK, *tmpPtrJ, *tmpPtrI;

  for ( k = 0; k < inDim[2]; k++, dptr += inInc2 )
271
  {
272 273
    sz1 = (k < 1)?(0):((k-1)/4);
    sz2 =             ((k  )/4);
274 275 276 277 278
    sz2 = ( k == inDim[2]-1 )?(sz1):(sz2);

    sz1 += outExt[4];
    sz2 += outExt[4];

279
    // Bounds check
280
    if (sz2 > outExt[5])
281
    {
282
      sz2 = outExt[5];
283
    }
284

285 286 287
    tmpPtrK = outBasePtr + sz1 * outInc2;

    for ( j = 0; j < inDim[1]; j++, dptr+= inInc1 )
288
    {
289 290
      sy1 = (j < 1)?(0):((j-1)/4);
      sy2 =             ((j  )/4);
291 292 293 294 295
      sy2 = ( j == inDim[1]-1 )?(sy1):(sy2);

      sy1 += outExt[2];
      sy2 += outExt[2];

296
      // Bounds check
297
      if (sy2 > outExt[3])
298
      {
299
        sy2 = outExt[3];
300
      }
301

302 303 304
      tmpPtrJ = tmpPtrK + sy1 * outInc1;

      for ( i = 0; i < inDim[0]; i++ )
305
      {
306 307
        sx1 = (i < 1)?(0):((i-1)/4);
        sx2 =             ((i  )/4);
308 309 310 311 312
        sx2 = ( i == inDim[0]-1 )?(sx1):(sx2);

        sx1 += outExt[0];
        sx2 += outExt[0];

313
        // Bounds check
314
        if (sx2 > outExt[1])
315
        {
316
          sx2 = outExt[1];
317
        }
318

319 320 321
        tmpPtrI = tmpPtrJ + sx1 * outInc0;

        for ( c = 0; c < nComponents; c++, tmpPtrI += 3 )
322
        {
323
          if ( independent )
324
          {
325 326
            val = static_cast<unsigned short>((*dptr + shift[c]) * scale[c]);
            ++dptr;
327
          }
328
          else
329
          {
330 331 332
            val = static_cast<unsigned short>((*(dptr+components-1) +
                   shift[components-1]) * scale[components-1]);
            dptr += components;
333
          }
334 335

          for ( z = sz1; z <= sz2; z++ )
336
          {
337
            for ( y = sy1; y <= sy2; y++ )
338
            {
339 340
              tmpPtr = tmpPtrI + (z-sz1)*outInc2 + (y-sy1)*outInc1;
              for ( x = sx1; x <= sx2; x++, tmpPtr += outInc0 )
341
              {
342
                if (val < tmpPtr[0])
343
                {
344
                  tmpPtr[0] = val;
345
                }
346
                if (val > tmpPtr[1])
347
                {
348 349 350 351 352 353 354 355
                  tmpPtr[1] = val;
                }
              }
            }
          }
        }
      }
    }
356
  }
357 358 359 360 361 362 363 364 365 366 367 368 369 370
}

//----------------------------------------------------------------------------
// Fill in the maximum gradient magnitude space leaping information.
template <class T>
void
vtkVolumeRayCastSpaceLeapingImageFilterMaxGradientMagnitudeExecute(
    vtkVolumeRayCastSpaceLeapingImageFilter *self,
    vtkImageData *inData,
    vtkImageData *outData, int outExt[6],
    T )
{
  // the number of independent components for which we need to keep track of
  // min/max
371
  const int nComponents = self->GetNumberOfIndependentComponents();
372 373 374 375 376 377 378


  // B. Now fill in the max-min-gradient volume structure

  // B.1 First compute the extents of the input that contribute to this structure

  int inExt[6], inWholeExt[6];
379
  int inDim[3];
380 381 382 383 384 385 386 387 388 389 390 391 392 393
  int outWholeDim[3];
  vtkVolumeRayCastSpaceLeapingImageFilter::ComputeInputExtentsForOutput(
    inExt, inDim, outExt, inData );
  inData->GetExtent(inWholeExt);
  outData->GetDimensions(outWholeDim);

  float shift[4], scale[4];
  self->GetTableShift(shift);
  self->GetTableScale(scale);


  // B.2 Get increments to march through the input extents

  vtkIdType inInc0, inInc1, inInc2;
394 395
  inData->GetContinuousIncrements(self->GetCurrentScalars(),
                                  inExt, inInc0, inInc1, inInc2);
396 397 398

  // Get increments to march through the output extents

399 400 401
  const vtkIdType outInc0 = 3*nComponents;
  const vtkIdType outInc1 = outInc0*outWholeDim[0];
  const vtkIdType outInc2 = outInc1*outWholeDim[1];
402 403 404

  // B.3 Now fill in the min-max volume.

405 406 407 408
  int i, j, k;
  int c;
  int sx1, sx2, sy1, sy2, sz1, sz2;
  int x, y, z;
409 410 411 412 413 414 415 416 417 418 419 420 421 422 423

  unsigned char val;
  unsigned short *outBasePtr = static_cast< unsigned short * >(
                                outData->GetScalarPointer());

  // The pointer into the space-leaping output volume.
  unsigned short *tmpPtr, *tmpPtrK, *tmpPtrJ, *tmpPtrI;

  // pointer to the slice of the gradient magnitude
  unsigned char **gsptr = self->GetGradientMagnitude();

  // Initialize pointer to the starting extents given by inExt.
  gsptr += (inExt[4]-inWholeExt[4]);

  for ( k = 0; k < inDim[2]; k++, ++gsptr )
424
  {
425 426
    sz1 = (k < 1)?(0):((k-1)/4);
    sz2 =             ((k  )/4);
427 428 429 430 431
    sz2 = ( k == inDim[2]-1 )?(sz1):(sz2);

    sz1 += outExt[4];
    sz2 += outExt[4];

432
    // Bounds check
433
    if (sz2 > outExt[5])
434
    {
435
      sz2 = outExt[5];
436
    }
437

438 439 440 441 442
    tmpPtrK = outBasePtr + sz1 * outInc2;

    unsigned char *gptr = *gsptr;

    for ( j = 0; j < inDim[1]; j++, gptr+= inInc1 )
443
    {
444 445
      sy1 = (j < 1)?(0):((j-1)/4);
      sy2 =             ((j  )/4);
446 447 448 449 450
      sy2 = ( j == inDim[1]-1 )?(sy1):(sy2);

      sy1 += outExt[2];
      sy2 += outExt[2];

451
      // Bounds check
452
      if (sy2 > outExt[3])
453
      {
454
        sy2 = outExt[3];
455
      }
456

457 458 459
      tmpPtrJ = tmpPtrK + sy1 * outInc1;

      for ( i = 0; i < inDim[0]; i++ )
460
      {
461 462
        sx1 = (i < 1)?(0):((i-1)/4);
        sx2 =             ((i  )/4);
463 464 465 466 467
        sx2 = ( i == inDim[0]-1 )?(sx1):(sx2);

        sx1 += outExt[0];
        sx2 += outExt[0];

468
        // Bounds check
469
        if (sx2 > outExt[1])
470
        {
471
          sx2 = outExt[1];
472
        }
473

474 475 476
        tmpPtrI = tmpPtrJ + sx1 * outInc0;

        for ( c = 0; c < nComponents; c++, tmpPtrI += 3 )
477
        {
478 479 480 481
          val = *gptr;
          ++gptr;

          for ( z = sz1; z <= sz2; z++ )
482
          {
483
            for ( y = sy1; y <= sy2; y++ )
484
            {
485 486
              tmpPtr = tmpPtrI + (z-sz1)*outInc2 + (y-sy1)*outInc1;
              for ( x = sx1; x <= sx2; x++, tmpPtr += outInc0 )
487
              {
488 489 490 491 492

                // Need to keep track of max gradient magnitude in upper
                // eight bits. No need to preserve lower eight (the flag)
                // since we will be recomputing this.
                if (val>(tmpPtr[2]>>8))
493
                {
494 495
                  tmpPtr[2] = (val<<8);
                }
496

497 498 499 500 501 502
              }
            }
          }
        }
      }
    }
503
  }
504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522
}

//----------------------------------------------------------------------------
// Optimized method that does both the following in one pass
// - Fill in the min-max space leaping information.
// - Fill in the maximum gradient magnitude space leaping information.
template <class T>
void
vtkVolumeRayCastSpaceLeapingImageFilterMinMaxAndMaxGradientMagnitudeExecute(
    vtkVolumeRayCastSpaceLeapingImageFilter *self,
    vtkImageData *inData,
    vtkImageData *outData, int outExt[6],
    T )
{
  // the number of independent components for which we need to keep track of
  // min/max
  vtkDataArray * scalars = self->GetCurrentScalars();
  const int components = scalars->GetNumberOfComponents();
  const int independent = self->GetIndependentComponents();
523
  const int nComponents = (independent) ? components : 1;
524 525 526 527 528


  // B.1 First compute the extents of the input that contribute to this structure

  int inExt[6], inWholeExt[6];
529
  int inDim[3];
530 531 532 533 534 535 536 537 538 539 540 541 542 543
  int outWholeDim[3];
  vtkVolumeRayCastSpaceLeapingImageFilter::ComputeInputExtentsForOutput(
    inExt, inDim, outExt, inData );
  inData->GetExtent(inWholeExt);
  outData->GetDimensions(outWholeDim);

  float shift[4], scale[4];
  self->GetTableShift(shift);
  self->GetTableScale(scale);


  // B.2 Get increments to march through the input extents

  vtkIdType inInc0, inInc1, inInc2;
544 545
  inData->GetContinuousIncrements(scalars,
                                  inExt, inInc0, inInc1, inInc2);
546 547 548

  // Get increments to march through the output extents

549 550 551
  const vtkIdType outInc0 = 3*nComponents;
  const vtkIdType outInc1 = outInc0*outWholeDim[0];
  const vtkIdType outInc2 = outInc1*outWholeDim[1];
552 553 554

  // B.3 Now fill in the min-max and gradient max structure

555 556 557 558
  int i, j, k;
  int c;
  int sx1, sx2, sy1, sy2, sz1, sz2;
  int x, y, z;
559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577


  T *dptr = static_cast< T * >(scalars->GetVoidPointer(0));
  unsigned char val;
  unsigned short minMaxVal;
  unsigned short *outBasePtr = static_cast< unsigned short * >(
                                outData->GetScalarPointer());

  // pointer to the slice of the gradient magnitude
  unsigned char **gsptr = self->GetGradientMagnitude();

  // Initialize pointers to the starting extents given by inExt.
  gsptr += (inExt[4]-inWholeExt[4]);  // pointer to slice gradient
  dptr += self->ComputeOffset( inExt, inWholeExt, nComponents );

  // The pointer into the space-leaping output volume.
  unsigned short *tmpPtr, *tmpPtrK, *tmpPtrJ, *tmpPtrI;

  for ( k = 0; k < inDim[2]; k++, dptr += inInc2, ++gsptr )
578
  {
579 580
    sz1 = (k < 1)?(0):((k-1)/4);
    sz2 =             ((k  )/4);
581 582 583 584 585
    sz2 = ( k == inDim[2]-1 )?(sz1):(sz2);

    sz1 += outExt[4];
    sz2 += outExt[4];

586
    // Bounds check
587
    if (sz2 > outExt[5])
588
    {
589
      sz2 = outExt[5];
590
    }
591

592 593 594 595 596
    tmpPtrK = outBasePtr + sz1 * outInc2;

    unsigned char *gptr = *gsptr;

    for ( j = 0; j < inDim[1]; j++, dptr+= inInc1, gptr+= inInc1 )
597
    {
598 599
      sy1 = (j < 1)?(0):((j-1)/4);
      sy2 =             ((j  )/4);
600 601 602 603 604
      sy2 = ( j == inDim[1]-1 )?(sy1):(sy2);

      sy1 += outExt[2];
      sy2 += outExt[2];

605
      // Bounds check
606
      if (sy2 > outExt[3])
607
      {
608
        sy2 = outExt[3];
609
      }
610

611 612 613
      tmpPtrJ = tmpPtrK + sy1 * outInc1;

      for ( i = 0; i < inDim[0]; i++ )
614
      {
615 616
        sx1 = (i < 1)?(0):((i-1)/4);
        sx2 =             ((i  )/4);
617 618 619 620 621
        sx2 = ( i == inDim[0]-1 )?(sx1):(sx2);

        sx1 += outExt[0];
        sx2 += outExt[0];

622
        // Bounds check
623
        if (sx2 > outExt[1])
624
        {
625
          sx2 = outExt[1];
626
        }
627

628 629 630
        tmpPtrI = tmpPtrJ + sx1 * outInc0;

        for ( c = 0; c < nComponents; c++, tmpPtrI += 3 )
631
        {
632 633 634 635
          val = *gptr;
          ++gptr;

          if ( independent )
636
          {
637 638
            minMaxVal = static_cast<unsigned short>((*dptr + shift[c]) * scale[c]);
            ++dptr;
639
          }
640
          else
641
          {
642 643 644
            minMaxVal = static_cast<unsigned short>((*(dptr+components-1) +
                   shift[components-1]) * scale[components-1]);
            dptr += components;
645
          }
646 647

          for ( z = sz1; z <= sz2; z++ )
648
          {
649
            for ( y = sy1; y <= sy2; y++ )
650
            {
651 652 653

              tmpPtr = tmpPtrI + (z-sz1)*outInc2 + (y-sy1)*outInc1;
              for ( x = sx1; x <= sx2; x++, tmpPtr += outInc0 )
654
              {
655 656

                if (minMaxVal<tmpPtr[0])
657
                {
658
                  tmpPtr[0] = minMaxVal;
659
                }
660
                if (minMaxVal>tmpPtr[1])
661
                {
662
                  tmpPtr[1] = minMaxVal;
663
                }
664
                if (val>(tmpPtr[2]>>8))
665
                {
666 667 668 669 670 671 672 673
                  tmpPtr[2] = (val<<8);
                }
              }
            }
          }
        }
      }
    }
674
  }
675 676 677 678 679 680 681 682 683
}

//----------------------------------------------------------------------------
void vtkVolumeRayCastSpaceLeapingImageFilter
::FillScalarAndGradientOpacityFlags( vtkImageData *outData, int outExt[6] )
{
  // Get increments to march through the output

  vtkIdType outInc0, outInc1, outInc2;
684
  outData->GetContinuousIncrements(outExt, outInc0, outInc1, outInc2);
685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703

  // Now process the flags

  unsigned short *tmpPtr = static_cast< unsigned short * >(
                outData->GetScalarPointerForExtent(outExt));
  unsigned short *minNonZeroScalarIndex
                     = this->GetMinNonZeroScalarIndex();
  unsigned char  *minNonZeroGradientMagnitudeIndex
                     = this->GetMinNonZeroGradientMagnitudeIndex();

  int i, j, k, c, loop;

  // the number of independent components for which we need to keep track of
  // min/max/gradient
  const int nComponents = this->GetNumberOfIndependentComponents();

  // Loop over the data with in the supplied extents

  for (k = outExt[4]; k <= outExt[5]; ++k, tmpPtr += outInc2)
704
  {
705
    for (j = outExt[2]; j <= outExt[3]; ++j, tmpPtr += outInc1)
706
    {
707
      for (i = outExt[0]; i <= outExt[1]; ++i)
708
      {
709
        for ( c = 0; c < nComponents; ++c, tmpPtr += 3 )
710
        {
711 712 713 714 715

          // We definite have 0 opacity because our maximum scalar value in
          // this region is below the minimum scalar value with non-zero opacity
          // for this component
          if ( tmpPtr[1] < minNonZeroScalarIndex[c] )
716
          {
717
            tmpPtr[2] &= 0xff00;
718
          }
719 720 721 722
          // We have 0 opacity because we are using gradient magnitudes and
          // the maximum gradient magnitude in this area is below the minimum
          // gradient magnitude with non-zero opacity for this component
          else if ( (tmpPtr[2]>>8) < minNonZeroGradientMagnitudeIndex[c] )
723
          {
724
            tmpPtr[2] &= 0xff00;
725
          }
726 727 728 729 730
          // We definitely have non-zero opacity because our minimum scalar
          // value is lower than our first scalar with non-zero opacity, and
          // the maximum scalar value is greater than this threshold - so
          // we must encounter scalars with opacity in between
          else if ( tmpPtr[0] < minNonZeroScalarIndex[c] )
731
          {
732 733
            tmpPtr[2] &= 0xff00;
            tmpPtr[2] |= 0x0001;
734
          }
735 736 737 738 739
          // We have to search between min scalar value and the
          // max scalar stored in the minmax volume to look for non-zero
          // opacity since both values must be above our first non-zero
          // threshold so we don't have information in this area
          else
740
          {
741
            for ( loop = tmpPtr[0]; loop <= tmpPtr[1]; ++loop )
742
            {
743
              if ( this->ScalarOpacityTable[c][loop] )
744
              {
745 746
                break;
              }
747
            }
748
            if ( loop <= tmpPtr[1] )
749
            {
750 751
              tmpPtr[2] &= 0xff00;
              tmpPtr[2] |= 0x0001;
752
            }
753
            else
754
            {
755 756 757 758 759 760
              tmpPtr[2] &= 0xff00;
            }
          }
        }
      }
    }
761
  }
762 763 764 765 766 767 768 769 770
}

//----------------------------------------------------------------------------
void vtkVolumeRayCastSpaceLeapingImageFilter
::FillScalarOpacityFlags( vtkImageData *outData, int outExt[6] )
{
  // Get increments to march through the output

  vtkIdType outInc0, outInc1, outInc2;
771
  outData->GetContinuousIncrements(outExt, outInc0, outInc1, outInc2);
772

773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788
  // Now process the flags

  unsigned short *tmpPtr = static_cast< unsigned short * >(
                outData->GetScalarPointerForExtent(outExt));
  unsigned short *minNonZeroScalarIndex
                     = this->GetMinNonZeroScalarIndex();

  int i, j, k, c, loop;

  // the number of independent components for which we need to keep track of
  // min/max/gradient
  const int nComponents = this->GetNumberOfIndependentComponents();

  // Loop over the data with in the supplied extents

  for (k = outExt[4]; k <= outExt[5]; ++k, tmpPtr += outInc2)
789
  {
790
    for (j = outExt[2]; j <= outExt[3]; ++j, tmpPtr += outInc1)
791
    {
792
      for (i = outExt[0]; i <= outExt[1]; ++i)
793
      {
794
        for ( c = 0; c < nComponents; ++c, tmpPtr += 3 )
795
        {
796 797 798 799 800

          // We definite have 0 opacity because our maximum scalar value in
          // this region is below the minimum scalar value with non-zero opacity
          // for this component
          if ( tmpPtr[1] < minNonZeroScalarIndex[c] )
801
          {
802
            tmpPtr[2] &= 0xff00;
803
          }
804 805 806 807 808
          // We definitely have non-zero opacity because our minimum scalar
          // value is lower than our first scalar with non-zero opacity, and
          // the maximum scalar value is greater than this threshold - so
          // we must encounter scalars with opacity in between
          else if ( tmpPtr[0] < minNonZeroScalarIndex[c] )
809
          {
810 811
            tmpPtr[2] &= 0xff00;
            tmpPtr[2] |= 0x0001;
812
          }
813 814 815 816 817
          // We have to search between min scalar value and the
          // max scalar stored in the minmax volume to look for non-zero
          // opacity since both values must be above our first non-zero
          // threshold so we don't have information in this area
          else
818
          {
819
            for ( loop = tmpPtr[0]; loop <= tmpPtr[1]; ++loop )
820
            {
821
              if ( this->ScalarOpacityTable[c][loop] )
822
              {
823 824
                break;
              }
825
            }
826
            if ( loop <= tmpPtr[1] )
827
            {
828 829
              tmpPtr[2] &= 0xff00;
              tmpPtr[2] |= 0x0001;
830
            }
831
            else
832
            {
833 834 835 836 837 838
              tmpPtr[2] &= 0xff00;
            }
          }
        }
      }
    }
839
  }
840 841 842 843 844
}

//----------------------------------------------------------------------------
void vtkVolumeRayCastSpaceLeapingImageFilter::ThreadedRequestData(
  vtkInformation *vtkNotUsed(request),
845 846
  vtkInformationVector **vtkNotUsed(inputVector),
  vtkInformationVector *vtkNotUsed(outputVector),
847 848
  vtkImageData ***inData,
  vtkImageData **outData,
849
  int outExt[6], int vtkNotUsed(id))
850
{
851 852 853
#ifdef vtkVolumeRayCastSpaceLeapingImageFilter_DEBUG
  std::cout << "Thread id = " << id << std::endl;
#endif
854 855 856 857 858

  // A. Initialize the data with a blank flag.

  // - Get the number of independent components for which we need to keep
  //   track of min/max
859
  if (!this->GetCurrentScalars())
860
  {
861
    return;
862
  }
863

864
  const int components = this->GetCurrentScalars()->GetNumberOfComponents();
865
  const int nComponents = (this->GetIndependentComponents()) ? components : 1;
866

867 868 869 870
  // Clear the output if we are computing the min-max. In other cases, we
  // will be re-using the cache. (See the method AllocateOutputData)

  if (this->ComputeMinMax)
871
  {
872 873
    vtkVolumeRayCastSpaceLeapingImageFilterClearOutput(
      outData[0], outExt, nComponents );
874
  }
875 876 877


  // If only scalar min-max need to be re-computed
878

879
  if (this->ComputeMinMax && !this->ComputeGradientOpacity)
880
  {
881 882
    int scalarType   = this->CurrentScalars->GetDataType();
    switch (scalarType)
883
    {
884 885 886 887 888 889 890 891
      vtkTemplateMacro(
        vtkVolumeRayCastSpaceLeapingImageFilterMinMaxExecute(
          this, inData[0][0], outData[0], outExt, static_cast<VTK_TT>(0))
        );
      default:
        vtkErrorMacro("Unknown scalar type");
        return;
    }
892
  }
893 894

  // If only gradient max needs to be recomputed
895

896
  else if (this->ComputeGradientOpacity && !this->ComputeMinMax)
897
  {
898 899
    int scalarType   = this->CurrentScalars->GetDataType();
    switch (scalarType)
900
    {
901 902 903 904 905 906 907 908
      vtkTemplateMacro(
        vtkVolumeRayCastSpaceLeapingImageFilterMaxGradientMagnitudeExecute(
          this, inData[0][0], outData[0], outExt, static_cast<VTK_TT>(0))
        );
      default:
        vtkErrorMacro("Unknown scalar type");
        return;
    }
909
  }
910 911 912

  // If both scalar min-max need to be re-computed and gradient max needs to
  // be re-computed
913

914
  else if (this->ComputeGradientOpacity && this->ComputeMinMax)
915
  {
916 917
    int scalarType   = this->CurrentScalars->GetDataType();
    switch (scalarType)
918
    {
919 920 921 922 923 924 925
      vtkTemplateMacro(
        vtkVolumeRayCastSpaceLeapingImageFilterMinMaxAndMaxGradientMagnitudeExecute(
          this, inData[0][0], outData[0], outExt, static_cast<VTK_TT>(0))
        );
      default:
        vtkErrorMacro("Unknown scalar type");
        return;
926
    }
927
  }
928

929

930 931 932 933 934
  // Update the flags now for this extent. There are two specialized methods
  // here, depending on what mode we are in, so that we may do the flag update
  // in one pass through the data.

  if (this->UpdateGradientOpacityFlags)
935
  {
936 937
    // Process the flags based on the computed min-max volume
    this->FillScalarAndGradientOpacityFlags(outData[0], outExt);
938
  }
939
  else
940
  {
941
    this->FillScalarOpacityFlags(outData[0], outExt);
942
  }
943 944 945 946 947 948 949 950 951
}

//----------------------------------------------------------------------------
// Override superclass method to maintain a last successful execution time
int vtkVolumeRayCastSpaceLeapingImageFilter::RequestData(
  vtkInformation* request,
  vtkInformationVector** inputVector,
  vtkInformationVector* outputVector)
{
952 953 954 955
#ifdef vtkVolumeRayCastSpaceLeapingImageFilter_DEBUG
  cout << "ComputingGradientOpacity: " << ComputeGradientOpacity;
  cout << " ComputingMinMax: " << ComputeMinMax << " UpdatingFlags: 1" << endl;
#endif
956 957 958 959 960 961

  // Find the first non-zero scalar opacity and gradient opacity points on
  // the respective transfer functions

  this->ComputeFirstNonZeroOpacityIndices();

962
  // The actual work is done in the line below.
963
  if (this->Superclass::RequestData(request, inputVector, outputVector))
964
  {
965 966 967

    // if we recomputed the first two shorts in the output, update this.
    if (this->ComputeGradientOpacity || this->ComputeMinMax)
968
    {
969
      this->LastMinMaxBuildTime.Modified();
970
    }
971 972 973 974 975

    // flags were rebuilt. update this.
    this->LastMinMaxFlagTime.Modified();

    return 1;
976
  }
977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014 1015 1016 1017 1018

  return 0;
}

//----------------------------------------------------------------------------
int vtkVolumeRayCastSpaceLeapingImageFilter::RequestInformation (
  vtkInformation       * request,
  vtkInformationVector** inputVector,
  vtkInformationVector * outputVector)
{
  this->vtkImageAlgorithm::RequestInformation(request, inputVector, outputVector);

  vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
  vtkInformation* outInfo = outputVector->GetInformationObject(0);

  // Output scalar type is unsigned short,
  // 3 unsigned short values are needed to represent the min, max and gradient,
  // flag values. This is to be done for each independent component.

  vtkDataObject::SetPointDataActiveScalarInfo(outInfo,
      VTK_UNSIGNED_SHORT,
      3 * this->GetNumberOfIndependentComponents() );

  // The whole extent of the output is the whole extent of the input divided
  // by the block size along each dimension

  int outWholeExt[6], inWholeExtent[6];
  inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), inWholeExtent);
  this->InternalRequestUpdateExtent(outWholeExt, inWholeExtent);

  outInfo->Set(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), outWholeExt, 6);
  outInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), outWholeExt, 6);

  return 1;
}

//----------------------------------------------------------------------------
int vtkVolumeRayCastSpaceLeapingImageFilter
::GetNumberOfIndependentComponents()
{
  // the number of independent components for which we need to keep track of
  // min/max
1019
  if(this->CurrentScalars)
1020
  {
1021 1022
    const int components = this->CurrentScalars->GetNumberOfComponents();
    return ((this->IndependentComponents) ? components : 1);
1023
  }
1024
  return 0;
1025 1026 1027 1028 1029 1030 1031 1032 1033
}

//----------------------------------------------------------------------------
void vtkVolumeRayCastSpaceLeapingImageFilter
::ComputeFirstNonZeroOpacityIndices()
{
  // Find the first non-zero scalar opacity and gradient opacity points on
  // the respective transfer functions

1034
  const int nComponents = this->GetNumberOfIndependentComponents();
1035 1036

  // Initialize these arrays.
1037
  delete [] this->MinNonZeroScalarIndex;
1038
  this->MinNonZeroScalarIndex = nullptr;
1039
  delete [] this->MinNonZeroGradientMagnitudeIndex;
1040
  this->MinNonZeroGradientMagnitudeIndex = nullptr;
1041 1042

  // Update the flags now
1043
  int i;
1044
  this->MinNonZeroScalarIndex = new unsigned short [nComponents];
1045
  for ( int c = 0; c < nComponents; c++ )
1046
  {
1047
    for ( i = 0; i < this->TableSize[c]; i++ )
1048
    {
1049
      if ( this->ScalarOpacityTable[c][i] )
1050
      {
1051 1052 1053
        break;
      }
    }
1054 1055
    this->MinNonZeroScalarIndex[c] = i;
  }
1056 1057

  this->MinNonZeroGradientMagnitudeIndex = new unsigned char [nComponents];
1058
  for ( int c = 0; c < nComponents; c++ )
1059
  {
1060
    for ( i = 0; i < 256; i++ )
1061
    {
1062
      if ( this->GradientOpacityTable[c][i] )
1063
      {
1064 1065 1066
        break;
      }
    }
1067 1068
    this->MinNonZeroGradientMagnitudeIndex[c] = i;
  }
1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110 1111 1112 1113 1114 1115 1116 1117
}

//----------------------------------------------------------------------------
unsigned short * vtkVolumeRayCastSpaceLeapingImageFilter
::GetMinNonZeroScalarIndex()
{
  return this->MinNonZeroScalarIndex;
}

//----------------------------------------------------------------------------
unsigned char * vtkVolumeRayCastSpaceLeapingImageFilter
::GetMinNonZeroGradientMagnitudeIndex()
{
  return this->MinNonZeroGradientMagnitudeIndex;
}

//----------------------------------------------------------------------------
void vtkVolumeRayCastSpaceLeapingImageFilter
::SetGradientMagnitude( unsigned char ** gradientMagnitude )
{
  this->GradientMagnitude = gradientMagnitude;
}

//----------------------------------------------------------------------------
unsigned char ** vtkVolumeRayCastSpaceLeapingImageFilter
::GetGradientMagnitude()
{
  return this->GradientMagnitude;
}

//----------------------------------------------------------------------------
void vtkVolumeRayCastSpaceLeapingImageFilter
::SetScalarOpacityTable( int c, unsigned short * t )
{
  this->ScalarOpacityTable[c] = t;
}

//----------------------------------------------------------------------------
void vtkVolumeRayCastSpaceLeapingImageFilter
::SetGradientOpacityTable( int c, unsigned short * t )
{
  this->GradientOpacityTable[c] = t;
}

//----------------------------------------------------------------------------
unsigned short * vtkVolumeRayCastSpaceLeapingImageFilter
::GetMinMaxVolume( int size[4] )
{
  if (vtkImageData *output = this->GetOutput())
1118
  {
1119 1120 1121 1122 1123 1124 1125 1126
    int dims[3];
    output->GetDimensions(dims);
    size[0] = dims[0];
    size[1] = dims[1];
    size[2] = dims[2];
    size[3] = this->GetNumberOfIndependentComponents();

    return static_cast< unsigned short * >(output->GetScalarPointer());
1127
  }
1128
  return nullptr;
1129 1130 1131 1132
}

//----------------------------------------------------------------------------
// Fill in the min-max space leaping information.
1133 1134
vtkIdType vtkVolumeRayCastSpaceLeapingImageFilter
::ComputeOffset( const int ext[6], const int wholeExt[6], int nComponents )
1135
{
1136 1137 1138 1139
  int wDim[3] = { wholeExt[1]-wholeExt[0]+1,
                  wholeExt[3]-wholeExt[2]+1,
                  wholeExt[5]-wholeExt[4]+1 };

1140 1141
  // computation is done in parts to avoid int overflow
  vtkIdType offset = ext[4]-wholeExt[4];
1142 1143 1144 1145 1146 1147
  offset *= wDim[1];
  offset += ext[2]-wholeExt[2];
  offset *= wDim[0];
  offset += ext[0]-wholeExt[0];
  offset *= nComponents;

1148 1149
  return offset;
}
1150 1151

//----------------------------------------------------------------------------
1152
// Allocate the output data, caching if necessary. Caching may result in
1153 1154 1155
// invalid outputs and should be turned on, only when this filter is used
// as an internal ivar of the vtkFixedPointVolumeRayCastMapper.
void vtkVolumeRayCastSpaceLeapingImageFilter
1156 1157 1158 1159
::AllocateOutputData(vtkImageData *output,
                     vtkInformation* outInfo,
                     int *uExtent)
{
1160 1161 1162 1163
  // set the extent to be the update extent
  output->SetExtent(uExtent);

  if (this->Cache)
1164
  {
1165

1166 1167 1168 1169 1170 1171 1172
    int extent[6];
    this->Cache->GetExtent(extent);
    if (extent[0] == uExtent[0] && extent[1] == uExtent[1] &&
        extent[2] == uExtent[2] && extent[3] == uExtent[3] &&
        extent[4] == uExtent[4] && extent[5] == uExtent[5] &&
        this->Cache->GetNumberOfScalarComponents() ==
                    output->GetNumberOfScalarComponents())
1173
    {
1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186
      // Reuse the cache since it has the same dimensions. We may not be
      // updating all flags

      // This is absolutely scary code, if used as a standard VTK imaging filter,
      // but since the filter will be used only as an iVar of
      // vtkFixedPointVolumeRayCastMapper, we need a caching mechanism to avoid
      // reallocation of memory and re-update of certain bits in the Min-max
      // structure. In the interest of speed, we resort to a wee bit of ugly
      // code.
      output->GetPointData()->SetScalars(
          this->Cache->GetPointData()->GetScalars() );
      return;
    }
1187
  }
1188 1189

  // Otherwise allocate output afresh
1190
  output->AllocateScalars(outInfo);
1191 1192 1193 1194
}

//----------------------------------------------------------------------------
vtkImageData *vtkVolumeRayCastSpaceLeapingImageFilter
1195
::AllocateOutputData(vtkDataObject *output, vtkInformation *outInfo)
1196
{
1197
  // Call the superclass method
1198
  return vtkImageAlgorithm::AllocateOutputData(output, outInfo);
1199 1200 1201 1202 1203
}

//----------------------------------------------------------------------------
#ifdef vtkVolumeRayCastSpaceLeapingImageFilter_DEBUG
void vtkVolumeRayCastSpaceLeapingImageFilter
1204
::WriteMinMaxVolume(
1205 1206 1207 1208 1209 1210 1211 1212 1213 1214
  int component, unsigned short *minMaxVolume,
  int minMaxVolumeSize[4],
  const char *filename )
{
  vtkImageData *image = vtkImageData::New();
  image->SetExtent(0, minMaxVolumeSize[0]-1,
                   0, minMaxVolumeSize[1]-1,
                   0, minMaxVolumeSize[2]-1);
  image->SetScalarTypeToUnsignedShort();
  image->AllocateScalars();
1215 1216 1217

  const int nComponents = minMaxVolumeSize[3];
  const int inc = nComponents * 3;
1218 1219 1220
  unsigned short *pSrc = minMaxVolume + component;
  unsigned short *pDst = static_cast< unsigned short * >(
                                image->GetScalarPointer());
1221 1222 1223 1224 1225 1226
  // Do computation in parts to avoid int overfloat
  vtkIdType nVoxels = minMaxVolumeSize[0];
  nVoxels *= minMaxVolumeSize[1]
  nVoxels *= minMaxVolumeSize[2];

  for (vtkIdType i = 0; i < nVoxels; ++i, pSrc += inc, ++pDst)
1227
  {
1228
    *pDst = *pSrc;
1229
  }
1230 1231 1232 1233 1234 1235 1236 1237 1238

  vtkMetaImageWriter *writer = vtkMetaImageWriter::New();
  writer->SetFileName(filename);
  writer->SetInput(image);
  writer->Write();

  writer->Delete();
  image->Delete();
}
1239

1240
#endif