vtkImageResliceMapper.cxx 48.2 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11
/*=========================================================================

  Program:   Visualization Toolkit
  Module:    vtkImageResliceMapper.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
12
     PURPOSE.  See the side copyright notice for more information.
13 14 15 16

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

17
#include "vtkImageSliceMapper.h"
18 19 20 21 22 23 24
#include "vtkRenderer.h"
#include "vtkCamera.h"
#include "vtkImageSlice.h"
#include "vtkImageData.h"
#include "vtkImageProperty.h"
#include "vtkLookupTable.h"
#include "vtkMath.h"
25
#include "vtkPoints.h"
26
#include "vtkMatrix4x4.h"
27
#include "vtkLinearTransform.h"
28 29 30 31 32 33
#include "vtkPlane.h"
#include "vtkStreamingDemandDrivenPipeline.h"
#include "vtkGarbageCollector.h"
#include "vtkInformation.h"
#include "vtkInformationVector.h"
#include "vtkImageResliceToColors.h"
34
#include "vtkAbstractImageInterpolator.h"
35
#include "vtkObjectFactory.h"
36

37 38 39
// A tolerance to compensate for roundoff errors
#define VTK_RESLICE_MAPPER_VOXEL_TOL 7.62939453125e-06

40
vtkStandardNewMacro(vtkImageResliceMapper);
41 42 43 44

//----------------------------------------------------------------------------
vtkImageResliceMapper::vtkImageResliceMapper()
{
45
  this->SliceMapper = vtkImageSliceMapper::New();
46 47 48 49
  this->ImageReslice = vtkImageResliceToColors::New();
  this->ResliceMatrix = vtkMatrix4x4::New();
  this->WorldToDataMatrix = vtkMatrix4x4::New();
  this->SliceToWorldMatrix = vtkMatrix4x4::New();
50

51
  this->JumpToNearestSlice = 0;
52
  this->AutoAdjustImageQuality = 1;
53
  this->SeparateWindowLevelOperation = 1;
54 55 56 57
  this->SlabType = VTK_IMAGE_SLAB_MEAN;
  this->SlabThickness = 0.0;
  this->SlabSampleFactor = 2;
  this->ImageSampleFactor = 1;
58
  this->ResampleToScreenPixels = 1;
59
  this->InternalResampleToScreenPixels = 0;
60
  this->ResliceNeedUpdate = 0;
61 62 63

  // streaming requires an output port
  this->SetNumberOfOutputPorts(1);
64 65 66 67 68
}

//----------------------------------------------------------------------------
vtkImageResliceMapper::~vtkImageResliceMapper()
{
69 70 71 72
  if (this->SliceMapper)
    {
    this->SliceMapper->Delete();
    }
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114
  if (this->ImageReslice)
    {
    this->ImageReslice->Delete();
    }
  if (this->ResliceMatrix)
    {
    this->ResliceMatrix->Delete();
    }
  if (this->WorldToDataMatrix)
    {
    this->WorldToDataMatrix->Delete();
    }
  if (this->SliceToWorldMatrix)
    {
    this->SliceToWorldMatrix->Delete();
    }
}

//----------------------------------------------------------------------------
void vtkImageResliceMapper::SetSlicePlane(vtkPlane *plane)
{
  if (this->SlicePlane == plane)
    {
    return;
    }
  if (this->SlicePlane)
    {
    this->SlicePlane->Delete();
    }
  if (!plane)
    {
    this->SlicePlane = vtkPlane::New();
    }
  else
    {
    this->SlicePlane = plane;
    plane->Register(this);
    }

  this->Modified();
}

115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
//----------------------------------------------------------------------------
void vtkImageResliceMapper::SetInterpolator(
  vtkAbstractImageInterpolator *interpolator)
{
  unsigned long mtime = this->ImageReslice->GetMTime();

  this->ImageReslice->SetInterpolator(interpolator);

  if (this->ImageReslice->GetMTime() > mtime)
    {
    this->Modified();
    }
}

//----------------------------------------------------------------------------
vtkAbstractImageInterpolator *vtkImageResliceMapper::GetInterpolator()
{
  return this->ImageReslice->GetInterpolator();
}

135
//----------------------------------------------------------------------------
136
void vtkImageResliceMapper::ReleaseGraphicsResources(vtkWindow *win)
137
{
138
  this->SliceMapper->ReleaseGraphicsResources(win);
139 140 141
}

//----------------------------------------------------------------------------
142
void vtkImageResliceMapper::Render(vtkRenderer *ren, vtkImageSlice *prop)
143
{
144 145
  if (this->ResliceNeedUpdate)
    {
146 147
    this->ImageReslice->SetInputConnection(
      this->GetInputConnection(0, 0));
148
    this->ImageReslice->UpdateWholeExtent();
149
    this->ResliceNeedUpdate = 0;
150 151
    }

152
  // apply checkerboard pattern (should have timestamps)
153
  vtkImageProperty *property = prop->GetProperty();
154 155
  if (property && property->GetCheckerboard() &&
      this->InternalResampleToScreenPixels &&
156
      !this->SeparateWindowLevelOperation &&
157
      this->SliceFacesCamera)
158 159 160 161
    {
    this->CheckerboardImage(this->ImageReslice->GetOutput(),
      ren->GetActiveCamera(), property);
    }
162

163
  // delegate to vtkImageSliceMapper
164 165
  this->SliceMapper->SetInputConnection(
    this->ImageReslice->GetOutputPort());
166 167
  this->SliceMapper->GetDataToWorldMatrix()->DeepCopy(
    this->SliceToWorldMatrix);
168 169 170 171
  // the mapper uses SliceFacesCamera to decide whether to use a polygon
  // for the texture versus using a quad the size of the window
  this->SliceMapper->SetSliceFacesCamera(
    (this->SliceFacesCamera && !this->SeparateWindowLevelOperation));
172 173 174
  this->SliceMapper->SetExactPixelMatch(this->InternalResampleToScreenPixels);
  this->SliceMapper->SetBorder( (this->Border ||
                                 this->InternalResampleToScreenPixels) );
175 176 177
  this->SliceMapper->SetBackground( (this->Background &&
    !(this->SliceFacesCamera && this->InternalResampleToScreenPixels &&
      !this->SeparateWindowLevelOperation) ) );
178
  this->SliceMapper->SetPassColorData(!this->SeparateWindowLevelOperation);
179 180 181 182 183 184 185 186
  this->SliceMapper->SetDisplayExtent(this->ImageReslice->GetOutputExtent());

  // render pass info for members of vtkImageStack
  this->SliceMapper->MatteEnable = this->MatteEnable;
  this->SliceMapper->ColorEnable = this->ColorEnable;
  this->SliceMapper->DepthEnable = this->DepthEnable;

  // let vtkImageSliceMapper do the rest of the work
187
  this->SliceMapper->SetNumberOfThreads(this->NumberOfThreads);
188
  this->SliceMapper->SetClippingPlanes(this->ClippingPlanes);
189
  this->SliceMapper->Render(ren, prop);
190 191
}

192
//----------------------------------------------------------------------------
193
void vtkImageResliceMapper::Update(int port)
194 195 196
{
  // I don't like to override Update, or call Modified() in Update,
  // but this allows updates to be forced where MTimes can't be used
197
  bool resampleToScreenPixels = (this->ResampleToScreenPixels != 0);
198
  vtkRenderer *ren = 0;
199 200 201 202 203 204

  if (this->AutoAdjustImageQuality && resampleToScreenPixels)
    {
    // only use image-size texture if image is smaller than render window,
    // since otherwise there is far less advantage in doing so
    vtkImageSlice *prop = this->GetCurrentProp();
205
    ren = this->GetCurrentRenderer();
206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246
    if (ren && prop)
      {
      int *rsize = ren->GetSize();
      int maxrsize = (rsize[0] > rsize[1] ? rsize[0] : rsize[1]);
      int *isize = this->GetInput()->GetDimensions();
      int maxisize = (isize[0] > isize[1] ? isize[0] : isize[1]);
      maxisize = (isize[2] > maxisize ? isize[2] : maxisize);
      if (maxisize <= maxrsize && maxisize <= 1024)
        {
        resampleToScreenPixels = (prop->GetAllocatedRenderTime() >= 1.0);
        }
      }
    }

  if (resampleToScreenPixels)
    {
    // force update if quality has increased to "ResampleToScreenPixels"
    if (!this->InternalResampleToScreenPixels)
      {
      this->Modified();
      }
    else
      {
      // force update if renderer size has changes, since the texture
      // size is equal to the renderer size for "ResampleToScreenPixels"
      if (!ren)
        {
        ren = this->GetCurrentRenderer();
        }
      if (ren)
        {
        int *extent = this->ImageReslice->GetOutputExtent();
        int *size = ren->GetSize();
        if (size[0] != (extent[1] - extent[0] + 1) ||
            size[1] != (extent[3] - extent[2] + 1))
          {
          this->Modified();
          }
        }
      }
    }
247 248
  else if (this->InternalResampleToScreenPixels)
    {
249
    // if execution reaches this point in the code, then the
250 251 252 253 254 255
    // rendering has just switched to interactive quality, and it is
    // necessary to force update if modified since the last update
    if (this->GetMTime() > this->UpdateTime.GetMTime())
      {
      this->Modified();
      }
256 257 258 259 260 261
    else
      {
      // don't switch yet: wait until the camera changes position,
      // which will cause the MTime to change
      resampleToScreenPixels = true;
      }
262 263
    }

264 265
  this->InternalResampleToScreenPixels = resampleToScreenPixels;

266 267 268 269 270 271 272
  // Always update if something else caused the input to update
  vtkImageData *input = this->GetInput();
  if (input && input->GetUpdateTime() > this->UpdateTime.GetMTime())
    {
    this->Modified();
    }

273
  this->Superclass::Update(port);
274
  this->UpdateTime.Modified();
275 276
}

277 278 279 280 281
void vtkImageResliceMapper::Update()
{
  this->Superclass::Update();
}

282 283 284 285 286 287 288
//----------------------------------------------------------------------------
int vtkImageResliceMapper::ProcessRequest(
  vtkInformation* request, vtkInformationVector** inputVector,
  vtkInformationVector* outputVector)
{
  if (request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_INFORMATION()))
    {
289
    // use superclass method to update some info
290
    this->Superclass::ProcessRequest(request, inputVector, outputVector);
291

292 293 294 295 296 297
    // need the prop and renderer
    vtkImageSlice *prop = this->GetCurrentProp();
    vtkRenderer *ren = this->GetCurrentRenderer();

    if (ren && prop)
      {
298 299
      vtkImageProperty *property = prop->GetProperty();

300 301
      // Get point/normal from camera
      if (this->SliceFacesCamera || this->SliceAtFocalPoint)
302
        {
303 304
        vtkCamera *camera = ren->GetActiveCamera();

305 306 307 308 309 310 311 312 313
        if (this->SliceFacesCamera)
          {
          double normal[3];
          camera->GetDirectionOfProjection(normal);
          normal[0] = -normal[0];
          normal[1] = -normal[1];
          normal[2] = -normal[2];
          this->SlicePlane->SetNormal(normal);
          }
314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371

        if (this->SliceAtFocalPoint)
          {
          double point[4];
          camera->GetFocalPoint(point);

          if (this->JumpToNearestSlice)
            {
            double normal[4];
            this->SlicePlane->GetNormal(normal);
            normal[3] = -vtkMath::Dot(point, normal);
            point[3] = 1.0;

            // convert normal to data coordinates
            double worldToData[16];
            vtkMatrix4x4 *dataToWorld = this->GetDataToWorldMatrix();
            vtkMatrix4x4::Transpose(*dataToWorld->Element, worldToData);
            vtkMatrix4x4::MultiplyPoint(worldToData, normal, normal);

            // find the slice orientation from the normal
            int k = 0;
            double maxsq = 0;
            double sumsq = 0;
            for (int i = 0; i < 3; i++)
              {
              double tmpsq = normal[i]*normal[i];
              sumsq += tmpsq;
              if (tmpsq > maxsq)
                {
                maxsq = tmpsq;
                k = i;
                }
              }

            // if the slice is not oblique
            if ((1.0 - maxsq/sumsq) < 1e-12)
              {
              // get the point in data coordinates
              vtkMatrix4x4::Invert(*dataToWorld->Element, worldToData);
              vtkMatrix4x4::MultiplyPoint(worldToData, point, point);

              // set the point to lie exactly on a slice
              double z = (point[k] - this->DataOrigin[k])/this->DataSpacing[k];
              if (z > VTK_INT_MIN && z < VTK_INT_MAX)
                {
                int j = vtkMath::Floor(z + 0.5);
                point[k] = j*this->DataSpacing[k] + this->DataOrigin[k];
                }

              // convert back to world coordinates
              dataToWorld->MultiplyPoint(point, point);
              }
            }

          this->SlicePlane->SetOrigin(point);
          }

        } // end of "Get point/normal from camera"
372 373 374 375 376 377 378 379 380 381 382

      // set the matrices
      this->UpdateResliceMatrix(ren, prop);

      // update the coords for the polygon to be textured
      this->UpdatePolygonCoords(ren);

      // set the reslice spacing/origin/extent/axes
      this->UpdateResliceInformation(ren);

      // set the reslice bits related to the property
383 384 385 386
      this->UpdateResliceInterpolation(property);

      // update anything related to the image coloring
      this->UpdateColorInformation(property);
387 388
      }

389 390
    // set the number of threads to use when executing
    this->ImageReslice->SetNumberOfThreads(this->NumberOfThreads);
391

392 393 394 395 396 397 398 399
    // delegate request to vtkImageReslice (generally not a good thing to
    // do, but I'm familiar with the vtkImageReslice code that gets called).
    return this->ImageReslice->ProcessRequest(
      request, inputVector, outputVector);
    }

  if(request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_UPDATE_EXTENT()))
    {
400 401 402 403 404 405 406 407 408 409 410 411 412 413 414
    if (this->Streaming)
      {
      // delegate request to vtkImageReslice (generally not a good thing to
      // do, but I'm familiar with the vtkImageReslice code that gets called).
      return this->ImageReslice->ProcessRequest(
        request, inputVector, outputVector);
      }
    else
      {
      vtkInformation *inInfo = inputVector[0]->GetInformationObject(0);
      int ext[6];
      inInfo->Get(vtkStreamingDemandDrivenPipeline::WHOLE_EXTENT(), ext);
      inInfo->Set(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT(), ext, 6);
      }
    return 1;
415 416
    }

417 418
  if(request->Has(vtkStreamingDemandDrivenPipeline::REQUEST_DATA()))
    {
419 420 421 422
    vtkInformation *outInfo = outputVector->GetInformationObject(0);
    vtkImageData *output = vtkImageData::SafeDownCast(
      outInfo->Get(vtkDataObject::DATA_OBJECT()));

423 424 425 426 427 428 429
    // set output extent to avoid re-execution
    output->GetInformation()->Set(vtkDataObject::DATA_EXTENT(),
      outInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_EXTENT()), 6);

    // do an update of Reslice on the next render
    this->ResliceNeedUpdate = 1;

430 431 432
    return 1;
    }

433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462
  return this->Superclass::ProcessRequest(request, inputVector, outputVector);
}

//----------------------------------------------------------------------------
// Update the WorldToData transformation matrix, which is just the
// inverse of the vtkProp3D matrix.
void vtkImageResliceMapper::UpdateWorldToDataMatrix(vtkImageSlice *prop)
{
  // copy the matrix, but only if it has changed (we do this to
  // preserve the modified time of the matrix)
  double tmpmat[16] = { 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1 };
  if (!prop->GetIsIdentity())
    {
    vtkMatrix4x4::Invert(*prop->GetMatrix()->Element, tmpmat);
    }
  double *mat = *this->WorldToDataMatrix->Element;
  for (int i = 0; i < 16; i++)
    {
    if (mat[i] != tmpmat[i])
      {
      this->WorldToDataMatrix->DeepCopy(tmpmat);
      break;
      }
    }
}

//----------------------------------------------------------------------------
// Update the SliceToWorld transformation matrix
void vtkImageResliceMapper::UpdateSliceToWorldMatrix(vtkCamera *camera)
{
463
  // Get slice plane in world coords by passing null as the prop matrix
464 465
  double plane[4];
  this->GetSlicePlaneInDataCoords(0, plane);
466 467

  // Make sure normal is facing towards camera
468
  vtkMatrix4x4 *viewMatrix = camera->GetViewTransformMatrix();
469
  double *ndop = viewMatrix->Element[2];
470
  if (vtkMath::Dot(ndop, plane) < 0)
471
    {
472 473 474 475
    plane[0] = -plane[0];
    plane[1] = -plane[1];
    plane[2] = -plane[2];
    plane[3] = -plane[3];
476 477
    }

478 479 480 481 482 483
  // The normal is the first three elements
  double *normal = plane;

  // The last element is -dot(normal, origin)
  double dp = -plane[3];

484 485 486 487 488 489 490 491 492 493 494 495 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 523 524 525 526 527 528 529 530 531 532 533 534 535 536
  // Compute rotation angle between camera axis and slice normal
  double vec[3];
  vtkMath::Cross(ndop, normal, vec);
  double costheta = vtkMath::Dot(ndop, normal);
  double sintheta = vtkMath::Norm(vec);
  double theta = atan2(sintheta, costheta);
  if (sintheta != 0)
    {
    vec[0] /= sintheta;
    vec[1] /= sintheta;
    vec[2] /= sintheta;
    }
  // convert to quaternion
  costheta = cos(0.5*theta);
  sintheta = sin(0.5*theta);
  double quat[4];
  quat[0] = costheta;
  quat[1] = vec[0]*sintheta;
  quat[2] = vec[1]*sintheta;
  quat[3] = vec[2]*sintheta;
  // convert to matrix
  double mat[3][3];
  vtkMath::QuaternionToMatrix3x3(quat, mat);

  // Create a slice-to-world transform matrix
  // The columns are v1, v2, normal
  vtkMatrix4x4 *sliceToWorld = this->SliceToWorldMatrix;

  double v1[3], v2[3];
  vtkMath::Multiply3x3(mat, viewMatrix->Element[0], v1);
  vtkMath::Multiply3x3(mat, viewMatrix->Element[1], v2);

  sliceToWorld->Element[0][0] = v1[0];
  sliceToWorld->Element[1][0] = v1[1];
  sliceToWorld->Element[2][0] = v1[2];
  sliceToWorld->Element[3][0] = 0.0;

  sliceToWorld->Element[0][1] = v2[0];
  sliceToWorld->Element[1][1] = v2[1];
  sliceToWorld->Element[2][1] = v2[2];
  sliceToWorld->Element[3][1] = 0.0;

  sliceToWorld->Element[0][2] = normal[0];
  sliceToWorld->Element[1][2] = normal[1];
  sliceToWorld->Element[2][2] = normal[2];
  sliceToWorld->Element[3][2] = 0.0;

  sliceToWorld->Element[0][3] = -dp*normal[0];
  sliceToWorld->Element[1][3] = -dp*normal[1];
  sliceToWorld->Element[2][3] = dp-dp*normal[2];
  sliceToWorld->Element[3][3] = 1.0;
}

537 538 539 540 541
//----------------------------------------------------------------------------
// Update the reslice matrix, which is the slice-to-data matrix.
void vtkImageResliceMapper::UpdateResliceMatrix(
  vtkRenderer *ren, vtkImageSlice *prop)
{
542 543 544 545 546
  // Save the old matrix
  double *matrixElements = *this->ResliceMatrix->Element;
  double oldMatrixElements[16];
  vtkMatrix4x4::DeepCopy(oldMatrixElements, matrixElements);

547 548 549
  // Get world-to-data matrix from the prop matrix
  this->UpdateWorldToDataMatrix(prop);

550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572
  // Check if prop matrix is orthonormal
  bool propMatrixIsOrthonormal = false;
  vtkMatrix4x4 *propMatrix = 0;
  if (!this->InternalResampleToScreenPixels)
    {
    static double tol = 1e-12;
    propMatrix = prop->GetMatrix();
    double *row0 = propMatrix->Element[0];
    double *row1 = propMatrix->Element[1];
    double *row2 = propMatrix->Element[2];
    propMatrixIsOrthonormal = (
      fabs(vtkMath::Dot(row0, row0) - 1.0) < tol &&
      fabs(vtkMath::Dot(row1, row1) - 1.0) < tol &&
      fabs(vtkMath::Dot(row2, row2) - 1.0) < tol &&
      fabs(vtkMath::Dot(row0, row1)) < tol &&
      fabs(vtkMath::Dot(row0, row2)) < tol &&
      fabs(vtkMath::Dot(row1, row2)) < tol);
    }

  // Compute SliceToWorld matrix from camera if prop matrix is not
  // orthonormal or if InternalResampleToScreenPixels is set
  if (this->InternalResampleToScreenPixels ||
      !propMatrixIsOrthonormal)
573 574 575 576 577
    {
    this->UpdateSliceToWorldMatrix(ren->GetActiveCamera());
    vtkMatrix4x4::Multiply4x4(
      this->WorldToDataMatrix, this->SliceToWorldMatrix, this->ResliceMatrix);
    }
578
  else
579
    {
580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609
    // Get the matrices used to compute the reslice matrix
    vtkMatrix4x4 *resliceMatrix = this->ResliceMatrix;
    vtkMatrix4x4 *viewMatrix =
      ren->GetActiveCamera()->GetViewTransformMatrix();

    // Get slice plane in world coords by passing null as the matrix
    double wplane[4];
    this->GetSlicePlaneInDataCoords(0, wplane);

    // Check whether normal is facing towards camera, the "ndop" is
    // the negative of the direction of projection for the camera
    double *ndop = viewMatrix->Element[2];
    double dotprod = vtkMath::Dot(ndop, wplane);

    // Get slice plane in data coords by passing the prop matrix, flip
    // normal to face the camera
    double plane[4];
    this->GetSlicePlaneInDataCoords(propMatrix, plane);
    if (dotprod < 0)
      {
      plane[0] = -plane[0];
      plane[1] = -plane[1];
      plane[2] = -plane[2];
      plane[3] = -plane[3];

      wplane[0] = -wplane[0];
      wplane[1] = -wplane[1];
      wplane[2] = -wplane[2];
      wplane[3] = -wplane[3];
      }
610

611 612 613 614 615 616 617 618 619 620 621 622
    // Find the largest component of the normal
    int maxi = 0;
    double maxv = 0.0;
    for (int i = 0; i < 3; i++)
      {
      double tmp = plane[i]*plane[i];
      if (tmp > maxv)
        {
        maxi = i;
        maxv = tmp;
        }
      }
623

624 625 626 627 628 629 630 631 632 633 634 635 636
    // Create the corresponding axis
    double axis[3];
    axis[0] = 0.0;
    axis[1] = 0.0;
    axis[2] = 0.0;
    axis[maxi] = ((plane[maxi] < 0.0) ? -1.0 : 1.0);

    // Create two orthogonal axes
    double saxis[3], taxis[3];
    taxis[0] = 0.0;
    taxis[1] = 1.0;
    taxis[2] = 0.0;
    if (maxi == 1)
637
      {
638 639
      taxis[1] = 0.0;
      taxis[2] = 1.0;
640
      }
641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658
    vtkMath::Cross(taxis, axis, saxis);

    // The normal is the first three elements
    double *normal = plane;

    // The last element is -dot(normal, origin)
    double dp = (-plane[3] +
                 wplane[0]*propMatrix->Element[0][3] +
                 wplane[1]*propMatrix->Element[1][3] +
                 wplane[2]*propMatrix->Element[2][3]);

    // Compute the rotation angle between the axis and the normal
    double vec[3];
    vtkMath::Cross(axis, normal, vec);
    double costheta = vtkMath::Dot(axis, normal);
    double sintheta = vtkMath::Norm(vec);
    double theta = atan2(sintheta, costheta);
    if (sintheta != 0)
659
      {
660 661 662
      vec[0] /= sintheta;
      vec[1] /= sintheta;
      vec[2] /= sintheta;
663
      }
664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713
    // convert to quaternion
    costheta = cos(0.5*theta);
    sintheta = sin(0.5*theta);
    double quat[4];
    quat[0] = costheta;
    quat[1] = vec[0]*sintheta;
    quat[2] = vec[1]*sintheta;
    quat[3] = vec[2]*sintheta;
    // convert to matrix
    double mat[3][3];
    vtkMath::QuaternionToMatrix3x3(quat, mat);

    // Create a slice-to-data transform matrix
    // The columns are v1, v2, normal
    double v1[3], v2[3];
    vtkMath::Multiply3x3(mat, saxis, v1);
    vtkMath::Multiply3x3(mat, taxis, v2);

    resliceMatrix->Element[0][0] = v1[0];
    resliceMatrix->Element[1][0] = v1[1];
    resliceMatrix->Element[2][0] = v1[2];
    resliceMatrix->Element[3][0] = 0.0;

    resliceMatrix->Element[0][1] = v2[0];
    resliceMatrix->Element[1][1] = v2[1];
    resliceMatrix->Element[2][1] = v2[2];
    resliceMatrix->Element[3][1] = 0.0;

    resliceMatrix->Element[0][2] = normal[0];
    resliceMatrix->Element[1][2] = normal[1];
    resliceMatrix->Element[2][2] = normal[2];
    resliceMatrix->Element[3][2] = 0.0;

    resliceMatrix->Element[0][3] = dp*(propMatrix->Element[2][0] - normal[0]) -
      (propMatrix->Element[0][3]*propMatrix->Element[0][0] +
       propMatrix->Element[1][3]*propMatrix->Element[1][0] +
       propMatrix->Element[2][3]*propMatrix->Element[2][0]);
    resliceMatrix->Element[1][3] = dp*(propMatrix->Element[2][1] - normal[1]) -
      (propMatrix->Element[0][3]*propMatrix->Element[0][1] +
       propMatrix->Element[1][3]*propMatrix->Element[1][1] +
       propMatrix->Element[2][3]*propMatrix->Element[2][1]);
    resliceMatrix->Element[2][3] = dp*(propMatrix->Element[2][2] - normal[2]) -
      (propMatrix->Element[0][3]*propMatrix->Element[0][2] +
       propMatrix->Element[1][3]*propMatrix->Element[1][2] +
       propMatrix->Element[2][3]*propMatrix->Element[2][2]);
    resliceMatrix->Element[3][3] = 1.0;

    // Compute the SliceToWorldMatrix
    vtkMatrix4x4::Multiply4x4(propMatrix, resliceMatrix,
      this->SliceToWorldMatrix);
714 715
    }

716 717 718
  // If matrix changed, mark as modified so that Reslice will update
  int matrixChanged = 0;
  for (int j = 0; j < 16; j++)
719
    {
720
    matrixChanged |= (matrixElements[j] != oldMatrixElements[j]);
721
    }
722
  if (matrixChanged)
723
    {
724
    this->ResliceMatrix->Modified();
725 726 727
    }
}

728 729 730 731 732 733 734
//----------------------------------------------------------------------------
// Do all the fancy math to set up the reslicing
void vtkImageResliceMapper::UpdateResliceInformation(vtkRenderer *ren)
{
  vtkMatrix4x4 *resliceMatrix = this->ResliceMatrix;
  vtkImageResliceToColors *reslice = this->ImageReslice;

735 736 737
  int extent[6];
  double spacing[3];
  double origin[3];
738

739 740 741 742 743
  // Get current spacing and origin
  reslice->GetOutputSpacing(spacing);
  reslice->GetOutputOrigin(origin);
  reslice->GetOutputExtent(extent);

744
  // Get the view matrix
745 746 747
  vtkCamera *camera = ren->GetActiveCamera();
  vtkMatrix4x4 *viewMatrix = camera->GetViewTransformMatrix();

748 749 750
  // Get slice plane in world coords by passing null as the matrix
  double plane[4];
  this->GetSlicePlaneInDataCoords(0, plane);
751

752 753
  // Check whether normal is facing towards camera, the "ndop" is
  // the negative of the direction of projection for the camera
754
  double *ndop = viewMatrix->Element[2];
755
  if (vtkMath::Dot(ndop, plane) < 0)
756
    {
757 758 759 760
    plane[0] = -plane[0];
    plane[1] = -plane[1];
    plane[2] = -plane[2];
    plane[3] = -plane[3];
761 762
    }

763 764 765
  // Get the z position of the slice in slice coords
  // (requires plane to be normalized by GetSlicePlaneInDataCoords)
  double z = (plane[2] - 2.0)*plane[3];
766

767
  if (this->InternalResampleToScreenPixels)
768
    {
769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798
    // Get the projection matrix
    double aspect = ren->GetTiledAspectRatio();
    vtkMatrix4x4 *projMatrix = camera->GetProjectionTransformMatrix(
                                 aspect, 0, 1);

    // Compute other useful matrices
    double worldToView[16];
    double viewToWorld[16];
    double planeWorldToView[16];
    vtkMatrix4x4::Multiply4x4(
      *projMatrix->Element, *viewMatrix->Element, worldToView);
    vtkMatrix4x4::Invert(worldToView, viewToWorld);
    vtkMatrix4x4::Transpose(viewToWorld, planeWorldToView);

    double worldToSlice[16];
    double viewToSlice[16];
    vtkMatrix4x4::Invert(*this->SliceToWorldMatrix->Element, worldToSlice);
    vtkMatrix4x4::Multiply4x4(worldToSlice, viewToWorld, viewToSlice);

    // Transform the plane into view coordinates, using the transpose
    // of the inverse of the world-to-view matrix
    vtkMatrix4x4::MultiplyPoint(planeWorldToView, plane, plane);

    // Compute the bounds in slice coords
    double xmin = VTK_DOUBLE_MAX;
    double xmax = -VTK_DOUBLE_MAX;
    double ymin = VTK_DOUBLE_MAX;
    double ymax = -VTK_DOUBLE_MAX;

    for (int i = 0; i < 4; i++)
799
      {
800 801 802 803 804 805 806 807 808
      // The four corners of the view
      double x = (((i & 1) == 0) ? -1.0 : 1.0);
      double y = (((i & 2) == 0) ? -1.0 : 1.0);

      double hpoint[4];
      hpoint[0] = x;
      hpoint[1] = y;
      hpoint[2] = 0.0;
      hpoint[3] = 1.0;
809

810
      if (fabs(plane[2]) < 1e-6)
811
        {
812 813 814
        // Looking at plane edge-on, just put some
        // points at front clipping plane, others at back plane
        hpoint[2] = (((i & 1) == 0) ? 0.0 : 1.0);
815
        }
816
      else
817
        {
818 819 820 821 822 823 824 825 826 827 828 829
        // Intersect with the slice plane
        hpoint[2] = - (x*plane[0] + y*plane[1] + plane[3])/plane[2];

        // Clip to the front and back clipping planes
        if (hpoint[2] < 0)
          {
          hpoint[2] = 0.0;
          }
        else if (hpoint[2] > 1)
          {
          hpoint[2] = 1.0;
          }
830 831
        }

832 833
      // Transform into slice coords
      vtkMatrix4x4::MultiplyPoint(viewToSlice, hpoint, hpoint);
834

835 836
      x = hpoint[0]/hpoint[3];
      y = hpoint[1]/hpoint[3];
837

838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867
      // Find min/max in slice coords
      if (x < xmin) { xmin = x; }
      if (x > xmax) { xmax = x; }
      if (y < ymin) { ymin = y; }
      if (y > ymax) { ymax = y; }
      }

    // The ResliceExtent is always set to the renderer size,
    // this is the maximum size ever required and sticking to
    // this size avoids any memory reallocation on GPU or CPU
    int *size = ren->GetSize();
    int xsize = ((size[0] <= 0) ? 1 : size[0]);
    int ysize = ((size[1] <= 0) ? 1 : size[1]);

    extent[0] = 0;
    extent[1] = xsize - 1;
    extent[2] = 0;
    extent[3] = ysize - 1;
    extent[4] = 0;
    extent[5] = 0;

    // Find the spacing
    spacing[0] = (xmax - xmin)/xsize;
    spacing[1] = (ymax - ymin)/ysize;

    // Corner of resliced plane, including half-pixel offset to
    // exactly match texels to pixels in the final rendering
    origin[0] = xmin + 0.5*spacing[0];
    origin[1] = ymin + 0.5*spacing[1];
    origin[2] = z;
868
    }
869 870 871 872 873 874 875 876 877 878 879 880 881
  else
    {
    // Compute texel spacing from image spacing
    double inputSpacing[3];
    this->GetInput()->GetSpacing(inputSpacing);
    inputSpacing[0] = fabs(inputSpacing[0]);
    inputSpacing[1] = fabs(inputSpacing[1]);
    inputSpacing[2] = fabs(inputSpacing[2]);
    for (int j = 0; j < 2; j++)
      {
      double xc = this->ResliceMatrix->Element[j][0];
      double yc = this->ResliceMatrix->Element[j][1];
      double zc = this->ResliceMatrix->Element[j][2];
882 883 884 885 886
      double s = (xc*xc*inputSpacing[0] +
                  yc*yc*inputSpacing[1] +
                  zc*zc*inputSpacing[2])/sqrt(xc*xc + yc*yc + zc*zc);
      s /= this->ImageSampleFactor;
      // only modify if difference is greater than roundoff tolerance
887
      if (fabs((s - spacing[j])/s) > 1e-12)
888 889 890
        {
        spacing[j] = s;
        }
891
      }
892

893 894 895 896 897
    // Find the bounds for the texture
    double xmin = VTK_DOUBLE_MAX;
    double xmax = -VTK_DOUBLE_MAX;
    double ymin = VTK_DOUBLE_MAX;
    double ymax = -VTK_DOUBLE_MAX;
898

899 900 901
    vtkPoints *points = this->SliceMapper->GetPoints();
    vtkIdType n = points->GetNumberOfPoints();
    if (n == 0)
902
      {
903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919
      double inputOrigin[3];
      this->GetInput()->GetOrigin(inputOrigin);
      xmin = inputOrigin[0];
      xmax = inputOrigin[0];
      ymin = inputOrigin[1];
      ymax = inputOrigin[1];
      }

    for (vtkIdType k = 0; k < n; k++)
      {
      double point[3];
      points->GetPoint(k, point);

      xmin = ((xmin < point[0]) ? xmin : point[0]);
      xmax = ((xmax > point[0]) ? xmax : point[0]);
      ymin = ((ymin < point[1]) ? ymin : point[1]);
      ymax = ((ymax > point[1]) ? ymax : point[1]);
920 921
      }

922
    double tol = VTK_RESLICE_MAPPER_VOXEL_TOL;
923 924 925 926 927 928 929 930 931 932
    int xsize = vtkMath::Floor((xmax - xmin)/spacing[0] + tol);
    int ysize = vtkMath::Floor((ymax - ymin)/spacing[1] + tol);
    if (this->Border == 0)
      {
      xsize += 1;
      ysize += 1;
      }
    if (xsize < 1) { xsize = 1; }
    if (ysize < 1) { ysize = 1; }

933 934 935
    // Keep old size if possible, to avoid memory reallocation
    if ((xsize - 1) > extent[1] || (ysize - 1) > extent[3] ||
        (0.9*extent[1]/xsize) > 1.0 || (0.9*extent[3]/ysize) > 1.0)
936 937 938 939
      {
      extent[1] = xsize - 1;
      extent[3] = ysize - 1;
      }
940 941 942 943
    extent[0] = 0;
    extent[2] = 0;
    extent[4] = 0;
    extent[5] = 0;
944

945 946
    double x0 = xmin + 0.5*spacing[0]*(this->Border != 0);
    double y0 = ymin + 0.5*spacing[1]*(this->Border != 0);
947

948 949 950 951 952 953 954 955 956 957 958 959
    double dx = x0 - origin[0];
    double dy = y0 - origin[1];
    double dz = z - origin[2];

    // only modify origin if it has changed by tolerance
    if (dx*dx + dy*dy + dz*dz > tol*tol*spacing[0]*spacing[1])
      {
      origin[0] = x0;
      origin[1] = y0;
      origin[2] = z;
      }
    }
960

961 962 963 964 965 966
  // Prepare for reslicing
  reslice->SetResliceAxes(resliceMatrix);
  reslice->SetOutputExtent(extent);
  reslice->SetOutputSpacing(spacing);
  reslice->SetOutputOrigin(origin);

967 968
  if ((this->SliceFacesCamera && this->InternalResampleToScreenPixels) ||
      this->SlabThickness > 0)
969 970 971 972 973 974 975 976 977 978
    {
    // if slice follows camera, use reslice to set the border
    reslice->SetBorder(this->Border);
    }
  else
    {
    // tell reslice to use a double-thickness border,
    // since the polygon geometry will dictate the actual size
    reslice->SetBorder(2);
    }
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
//----------------------------------------------------------------------------
// Do all the fancy math to set up the reslicing
void vtkImageResliceMapper::UpdateColorInformation(vtkImageProperty *property)
{
  vtkScalarsToColors *lookupTable = this->DefaultLookupTable;

  if (property)
    {
    double colorWindow = property->GetColorWindow();
    double colorLevel = property->GetColorLevel();
    if (property->GetLookupTable())
      {
      lookupTable = property->GetLookupTable();
      if (!property->GetUseLookupTableScalarRange())
        {
        lookupTable->SetRange(colorLevel - 0.5*colorWindow,
                              colorLevel + 0.5*colorWindow);
        }
      }
    else
      {
      lookupTable->SetRange(colorLevel - 0.5*colorWindow,
                            colorLevel + 0.5*colorWindow);
      }
    }
  else
    {
    lookupTable->SetRange(0, 255);
    }
1010
  this->ImageReslice->SetBypass(this->SeparateWindowLevelOperation != 0);
1011
  this->ImageReslice->SetLookupTable(lookupTable);
1012 1013 1014 1015 1016 1017 1018 1019 1020 1021
  double backgroundColor[4] = { 0.0, 0.0, 0.0, 0.0 };
  if (this->Background)
    {
    this->GetBackgroundColor(property, backgroundColor);
    backgroundColor[0] *= 255;
    backgroundColor[1] *= 255;
    backgroundColor[2] *= 255;
    backgroundColor[3] *= 255;
    }
  this->ImageReslice->SetBackgroundColor(backgroundColor);
1022 1023
}

1024
//----------------------------------------------------------------------------
1025
// Set the reslice interpolation and slab thickness
1026 1027 1028 1029 1030
void vtkImageResliceMapper::UpdateResliceInterpolation(
  vtkImageProperty *property)
{
  // set the interpolation mode and border
  int interpMode = VTK_RESLICE_NEAREST;
1031
  int slabSlices = 1;
1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046

  if (property)
    {
    switch(property->GetInterpolationType())
      {
      case VTK_NEAREST_INTERPOLATION:
        interpMode = VTK_RESLICE_NEAREST;
        break;
      case VTK_LINEAR_INTERPOLATION:
        interpMode = VTK_RESLICE_LINEAR;
        break;
      case VTK_CUBIC_INTERPOLATION:
        interpMode = VTK_RESLICE_CUBIC;
        break;
      }
1047
    }
1048

1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068
  // set up the slice spacing for slab views
  double spacing[3], inputSpacing[3];
  this->ImageReslice->GetOutputSpacing(spacing);
  this->GetInput()->GetSpacing(inputSpacing);
  inputSpacing[0] = fabs(inputSpacing[0]);
  inputSpacing[1] = fabs(inputSpacing[1]);
  inputSpacing[2] = fabs(inputSpacing[2]);
  double xc = this->ResliceMatrix->Element[2][0];
  double yc = this->ResliceMatrix->Element[2][1];
  double zc = this->ResliceMatrix->Element[2][2];
  spacing[2] = (xc*xc*inputSpacing[0] +
                yc*yc*inputSpacing[1] +
                zc*zc*inputSpacing[2])/sqrt(xc*xc + yc*yc + zc*zc);

  // slab slice spacing is half the input slice spacing
  int n = vtkMath::Ceil(this->SlabThickness/spacing[2]);
  slabSlices = 1 + this->SlabSampleFactor*n;
  if (slabSlices > 1)
    {
    spacing[2] = this->SlabThickness/(slabSlices - 1);
1069
    }
1070 1071
  this->ImageReslice->SetOutputSpacing(spacing);
  int slabMode = this->SlabType;
1072 1073 1074 1075 1076 1077 1078
  double scalarScale = 1.0;
  if (slabMode == VTK_IMAGE_SLAB_SUM)
    {
    // "sum" means integrating over the path length of each ray through
    // the volume, so we need to include the sample spacing as a factor
    scalarScale = spacing[2];
    }
1079 1080

  this->ImageReslice->SetInterpolationMode(interpMode);
1081 1082
  this->ImageReslice->SetSlabMode(slabMode);
  this->ImageReslice->SetSlabNumberOfSlices(slabSlices);
1083
  this->ImageReslice->SetScalarScale(scalarScale);
1084
  this->ImageReslice->SlabTrapezoidIntegrationOn();
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 1118 1119 1120 1121 1122 1123 1124 1125 1126 1127 1128 1129 1130 1131 1132 1133 1134 1135
//----------------------------------------------------------------------------
void vtkImageResliceMapper::CheckerboardImage(
  vtkImageData *input, vtkCamera *camera, vtkImageProperty *property)
{
  // Use focal point as center of checkerboard pattern.  This guarantees
  // exactly the same checkerboard for all images in the scene, which is
  // useful when doing multiple overlays.
  double focalPoint[4];
  camera->GetFocalPoint(focalPoint);
  focalPoint[3] = 1.0;

  double worldToSlice[16];
  vtkMatrix4x4::Invert(*this->SliceToWorldMatrix->Element, worldToSlice);

  vtkMatrix4x4::MultiplyPoint(worldToSlice, focalPoint, focalPoint);
  if (focalPoint[3] != 0.0)
    {
    focalPoint[0] /= focalPoint[3];
    focalPoint[1] /= focalPoint[3];
    focalPoint[2] /= focalPoint[3];
    }

  // Get the checkerboard spacing and apply the offset fraction
  double checkSpacing[2], checkOffset[2];
  property->GetCheckerboardSpacing(checkSpacing);
  property->GetCheckerboardOffset(checkOffset);
  checkOffset[0] = checkOffset[0]*checkSpacing[0] + focalPoint[0];
  checkOffset[1] = checkOffset[1]*checkSpacing[1] + focalPoint[1];

  // Adjust according to the origin and spacing of the slice data
  double origin[3], spacing[3];
  input->GetSpacing(spacing);
  input->GetOrigin(origin);
  checkOffset[0] = (checkOffset[0] - origin[0])/spacing[0];
  checkOffset[1] = (checkOffset[1] - origin[1])/spacing[1];
  checkSpacing[0] /= spacing[0],
  checkSpacing[1] /= spacing[1];

  // Apply the checkerboard to the data
  int extent[6];
  input->GetExtent(extent);
  unsigned char *data = static_cast<unsigned char *>(
    input->GetScalarPointerForExtent(extent));

  vtkImageMapper3D::CheckerboardRGBA(
    data, extent[1] - extent[0] + 1, extent[3] - extent[2] + 1,
    checkOffset[0], checkOffset[1], checkSpacing[0], checkSpacing[1]);
}

1136
//----------------------------------------------------------------------------
1137
// Compute the vertices of the polygon in the slice coordinate system
1138 1139
#define VTK_IRM_MAX_VERTS 32
#define VTK_IRM_MAX_COORDS 96
1140
void vtkImageResliceMapper::UpdatePolygonCoords(vtkRenderer *ren)
1141
{
1142 1143 1144 1145 1146 1147
  // Get the projection matrix
  double aspect = ren->GetTiledAspectRatio();
  vtkCamera *camera = ren->GetActiveCamera();
  vtkMatrix4x4 *viewMatrix = camera->GetViewTransformMatrix();
  vtkMatrix4x4 *projMatrix = camera->GetProjectionTransformMatrix(
                               aspect, 0, 1);
1148

1149 1150 1151 1152 1153 1154
  // Compute other useful matrices
  double worldToView[16];
  double viewToWorld[16];
  vtkMatrix4x4::Multiply4x4(
    *projMatrix->Element, *viewMatrix->Element, worldToView);
  vtkMatrix4x4::Invert(worldToView, viewToWorld);
1155

1156 1157 1158 1159
  double worldToSlice[16];
  double viewToSlice[16];
  vtkMatrix4x4::Invert(*this->SliceToWorldMatrix->Element, worldToSlice);
  vtkMatrix4x4::Multiply4x4(worldToSlice, viewToWorld, viewToSlice);
1160

1161 1162 1163
  // Get slice plane in world coords by passing null as the matrix
  double plane[4];
  this->GetSlicePlaneInDataCoords(0, plane);
1164

1165 1166 1167 1168 1169 1170 1171 1172 1173 1174 1175 1176 1177 1178 1179 1180 1181 1182 1183 1184 1185 1186 1187 1188 1189 1190 1191 1192 1193 1194 1195 1196 1197 1198 1199 1200 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212
  // Check whether normal is facing towards camera, the "ndop" is
  // the negative of the direction of projection for the camera
  double *ndop = viewMatrix->Element[2];
  if (vtkMath::Dot(ndop, plane) < 0)
    {
    plane[0] = -plane[0];
    plane[1] = -plane[1];
    plane[2] = -plane[2];
    plane[3] = -plane[3];
    }

  // Get the z position of the slice in slice coords
  // (requires plane to be normalized by GetSlicePlaneInDataCoords)
  double z = (plane[2] - 2.0)*plane[3];

  // Generate a tolerance based on the screen pixel size
  double fpoint[4];
  camera->GetFocalPoint(fpoint);
  fpoint[3] = 1.0;
  vtkMatrix4x4::MultiplyPoint(worldToView, fpoint, fpoint);
  fpoint[0] /= fpoint[3];
  fpoint[1] /= fpoint[3];
  fpoint[2] /= fpoint[3];
  fpoint[3] = 1.0;

  double topOfScreen[4], botOfScreen[4];
  fpoint[1] -= 1.0;
  vtkMatrix4x4::MultiplyPoint(viewToWorld, fpoint, topOfScreen);
  fpoint[1] += 2.0;
  vtkMatrix4x4::MultiplyPoint(viewToWorld, fpoint, botOfScreen);

  topOfScreen[0] /= topOfScreen[3];
  topOfScreen[1] /= topOfScreen[3];
  topOfScreen[2] /= topOfScreen[3];
  topOfScreen[3] = 1.0;

  botOfScreen[0] /= botOfScreen[3];
  botOfScreen[1] /= botOfScreen[3];
  botOfScreen[2] /= botOfScreen[3];
  botOfScreen[3] = 1.0;

  // height of view in world coords at focal point
  double viewHeight =
    sqrt(vtkMath::Distance2BetweenPoints(topOfScreen, botOfScreen));

  // height of view in pixels
  int height = ren->GetSize()[1];

1213
  double tol = (height == 0 ? 0.5 : viewHeight*0.5/height);
1214 1215

  // make the data bounding box (with or without border)
1216
  double b = (this->Border ? 0.5 : VTK_RESLICE_MAPPER_VOXEL_TOL);
1217 1218 1219
  double bounds[6];
  for (int ii = 0; ii < 3; ii++)
    {
1220
    double c = b*this->DataSpacing[ii];
1221 1222
    int lo = this->DataWholeExtent[2*ii];
    int hi = this->DataWholeExtent[2*ii+1];
1223
    if (lo == hi && tol > c)
1224
      { // apply tolerance to avoid degeneracy
1225
      c = tol;
1226
      }
1227 1228
    bounds[2*ii]   = lo*this->DataSpacing[ii] + this->DataOrigin[ii] - c;
    bounds[2*ii+1] = hi*this->DataSpacing[ii] + this->DataOrigin[ii] + c;
1229 1230 1231
    }

  // transform the vertices to the slice coord system
1232 1233 1234
  double xpoints[8], ypoints[8];
  double weights1[8], weights2[8];
  bool above[8], below[8];
1235
  double mat[16];
1236 1237 1238 1239
  vtkMatrix4x4::Multiply4x4(*this->WorldToDataMatrix->Element,
                            *this->SliceToWorldMatrix->Element, mat);
  vtkMatrix4x4::Invert(mat, mat);

1240 1241
  // arrays for the list of polygon points
  int n = 0;
1242 1243
  double newxpoints[VTK_IRM_MAX_VERTS];
  double newypoints[VTK_IRM_MAX_VERTS];
1244 1245 1246
  double cx = 0.0;
  double cy = 0.0;

1247 1248 1249 1250 1251 1252 1253 1254 1255 1256
  for (int i = 0; i < 8; i++)
    {
    double point[4];
    point[0] = bounds[0 + ((i>>0)&1)];
    point[1] = bounds[2 + ((i>>1)&1)];
    point[2] = bounds[4 + ((i>>2)&1)];
    point[3] = 1.0;
    vtkMatrix4x4::MultiplyPoint(mat, point, point);
    xpoints[i] = point[0]/point[3];
    ypoints[i] = point[1]/point[3];
1257 1258 1259 1260 1261 1262 1263 1264 1265 1266 1267 1268 1269
    weights1[i] = point[2]/point[3] - z - 0.5*this->SlabThickness;
    weights2[i] = weights1[i] + this->SlabThickness;
    below[i] = (weights1[i] < 0);
    above[i] = (weights2[i] >= 0);

    if (this->SlabThickness > 0 && above[i] && below[i])
      {
      newxpoints[n] = xpoints[i];
      newypoints[n] = ypoints[i];
      cx += xpoints[i];
      cy += ypoints[i];
      n++;
      }
1270 1271
    }

1272
  // go through the edges and find the new points
1273 1274 1275 1276 1277 1278
  for (int j = 0; j < 12; j++)
    {
    // verts from edges (sorry about this..)
    int i1 = (j & 3) | (((j<<1) ^ (j<<2)) & 4);
    int i2 = (i1 ^ (1 << (j>>2)));

1279
    double *weights = weights2;
1280
    bool *side = above;
1281 1282
    int m = 1 + (this->SlabThickness > 0);
    for (int k = 0; k < m; k++)
1283
      {
1284 1285 1286 1287 1288 1289 1290 1291 1292 1293 1294 1295 1296 1297
      if (side[i1] ^ side[i2])
        {
        double w1 = weights[i2];
        double w2 = -weights[i1];
        double x = (w1*xpoints[i1] + w2*xpoints[i2])/(w1 + w2);
        double y = (w1*ypoints[i1] + w2*ypoints[i2])/(w1 + w2);
        newxpoints[n] = x;
        newypoints[n] = y;
        cx += x;
        cy += y;
        n++;
        }
      weights = weights1;
      side = below;
1298 1299 1300
      }
    }

1301
  double coords[VTK_IRM_MAX_COORDS];
1302

1303 1304
  if (n > 0)
    {
1305
    // centroid
1306 1307 1308 1309
    cx /= n;
    cy /= n;

    // sort the points to make a convex polygon
1310
    double angles[VTK_IRM_MAX_VERTS];
1311 1312 1313 1314 1315 1316 1317 1318 1319 1320 1321 1322 1323 1324
    for (int k = 0; k < n; k++)
      {
      double x = newxpoints[k];
      double y = newypoints[k];
      double t = atan2(y - cy, x - cx);
      int kk;
      for (kk = 0; kk < k; kk++)
        {
        if (t < angles[kk]) { break; }
        }
      for (int jj = k; jj > kk; --jj)
        {
        int jj3 = jj*3;
        angles[jj] = angles[jj-1];
1325 1326 1327
        coords[jj3] = coords[jj3-3];
        coords[jj3+1] = coords[jj3-2];
        coords[jj3+2] = coords[jj3-1];
1328 1329 1330 1331 1332 1333 1334 1335
        }
      int kk3 = kk*3;
      angles[kk] = t;
      coords[kk3] = x;
      coords[kk3+1] = y;
      coords[kk3+2] = z;
      }
    }
1336

1337 1338 1339 1340 1341 1342 1343 1344 1345 1346 1347 1348 1349 1350 1351 1352 1353 1354 1355 1356 1357 1358 1359 1360 1361 1362 1363 1364 1365 1366 1367 1368 1369 1370 1371 1372 1373 1374 1375 1376 1377 1378 1379 1380 1381 1382 1383 1384 1385 1386 1387 1388 1389 1390 1391 1392 1393 1394 1395 1396 1397 1398 1399
  // remove degenerate points
  if (n > 0)
    {
    bool found = true;
    int m = 0;
    do
      {
      m = 0;
      double xl = coords[3*(n-1)+0];
      double yl = coords[3*(n-1)+1];
      for (int k = 0; k < n; k++)
        {
        double x = coords[3*k+0];
        double y = coords[3*k+1];

        if (((x - xl)*(x - xl) + (y - yl)*(y - yl)) > tol*tol)
          {
          coords[3*m+0] = x;
          coords[3*m+1] = y;
          xl = x;
          yl = y;
          m++;
          }
        }
      found = (m < n);
      n = m;
      }
    while (found && n > 0);
    }

  // find convex hull
  if (this->SlabThickness > 0 && n > 0)
    {
    bool found = true;
    int m = 0;
    do
      {
      m = 0;
      double xl = coords[3*(n-1)+0];
      double yl = coords[3*(n-1)+1];
      for (int k = 0; k < n; k++)
        {
        double x = coords[3*k+0];
        double y = coords[3*k+1];
        int k1 = (k + 1) % n;
        double xn = coords[3*k1+0];
        double yn = coords[3*k1+1];

        if ((xn-xl)*(y-yl) - (yn-yl)*(x-xl) < tol*tol)
          {
          coords[3*m+0] = x;
          coords[3*m+1] = y;
          xl = x;
          yl = y;
          m++;
          }
        }
      found = (m < n);
      n = m;
      }
    while (found && n > 0);
    }

1400 1401 1402 1403 1404 1405 1406 1407
  vtkPoints *points = this->SliceMapper->GetPoints();
  if (!points)
    {
    points = vtkPoints::New();
    points->SetDataTypeToDouble();
    this->SliceMapper->SetPoints(points);
    points->Delete();
    }
1408

1409 1410
  points->SetNumberOfPoints(n);
  for (int k = 0; k < n; k++)
1411
    {
1412
    points->SetPoint(k, &coords[3*k]);
1413 1414
    }
}
1415 1416 1417 1418 1419

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

1421 1422
  os << indent << "JumpToNearestSlice: "
     << (this->JumpToNearestSlice ? "On\n" : "Off\n");
1423 1424
  os << indent << "AutoAdjustImageQuality: "
     << (this->AutoAdjustImageQuality ? "On\n" : "Off\n");
1425 1426
  os << indent << "SeparateWindowLevelOperation: "
     << (this->SeparateWindowLevelOperation ? "On\n" : "Off\n");
1427 1428
  os << indent << "ResampleToScreenPixels: "
     << (this->ResampleToScreenPixels ? "On\n" : "Off\n");
1429 1430 1431 1432
  os << indent << "SlabThickness: " << this->SlabThickness << "\n";
  os << indent << "SlabType: " << this->GetSlabTypeAsString() << "\n";
  os << indent << "SlabSampleFactor: " << this->SlabSampleFactor << "\n";
  os << indent << "ImageSampleFactor: " << this->ImageSampleFactor << "\n";
1433
  os << indent << "Interpolator: " << this->GetInterpolator() << "\n";
1434 1435 1436 1437 1438 1439 1440 1441 1442 1443 1444
}

//----------------------------------------------------------------------------
const char *vtkImageResliceMapper::GetSlabTypeAsString()
{
  switch (this->SlabType)
    {
    case VTK_IMAGE_SLAB_MIN:
      return "Min";
    case VTK_IMAGE_SLAB_MAX:
      return "Max";
1445 1446
    case VTK_IMAGE_SLAB_MEAN:
      return "Mean";
1447 1448
    case VTK_IMAGE_SLAB_SUM:
      return "Sum";
1449 1450
    }
  return "";
1451 1452 1453 1454 1455 1456 1457
}

//----------------------------------------------------------------------------
unsigned long vtkImageResliceMapper::GetMTime()
{
  unsigned long mTime = this->Superclass::GetMTime();

1458 1459 1460 1461 1462 1463 1464 1465 1466 1467 1468 1469
  // Check whether interpolator has changed
  vtkAbstractImageInterpolator *interpolator =
    this->ImageReslice->GetInterpolator();
  if (interpolator)
    {
    unsigned long mTime2 = interpolator->GetMTime();
    if (mTime2 > mTime)
      {
      mTime = mTime2;
      }
    }

1470 1471
  // Include camera in MTime so that REQUEST_INFORMATION
  // will be called if the camera changes
1472 1473
  if (this->SliceFacesCamera || this->SliceAtFocalPoint ||
      this->InternalResampleToScreenPixels)
1474
    {
1475 1476
    vtkRenderer *ren = this->GetCurrentRenderer();
    if (ren)
1477
      {
1478
      vtkCamera *camera = ren->GetActiveCamera();
1479
      unsigned long mTime2 = camera->GetMTime();
1480
      mTime = (mTime2 > mTime ? mTime2 : mTime);
1481 1482 1483 1484 1485 1486
      }
    }

  if (!this->SliceFacesCamera || !this->SliceAtFocalPoint)
    {
    unsigned long sTime = this->SlicePlane->GetMTime();
1487 1488 1489 1490 1491 1492
    mTime = (sTime > mTime ? sTime : mTime);
    }

  vtkImageSlice *prop = this->GetCurrentProp();
  if (prop != NULL)
    {
1493 1494
    unsigned long mTime2 = prop->GetUserTransformMatrixMTime();
    mTime = (mTime2 > mTime ? mTime2 : mTime);
1495 1496 1497 1498

    vtkImageProperty *property = prop->GetProperty();
    if (property != NULL)
      {
1499 1500 1501 1502 1503 1504 1505 1506 1507 1508 1509
      bool useMTime = true;
      if (this->SeparateWindowLevelOperation)
        {
        // only care about property if interpolation mode has changed,
        // since interpolation is the only property-related operation
        // done by vtkImageReslice if SeparateWindowLevelOperation is on
        int imode = this->ImageReslice->GetInterpolationMode();
        this->UpdateResliceInterpolation(property);
        useMTime = (imode != this->ImageReslice->GetInterpolationMode());
        }
      if (useMTime)
1510
        {
1511
        mTime2 = property->GetMTime();
1512
        mTime = (mTime2 > mTime ? mTime2 : mTime);
1513 1514 1515 1516 1517 1518 1519 1520

        vtkScalarsToColors *lookupTable = property->GetLookupTable();
        if (lookupTable != NULL)
          {
          // check the lookup table mtime
          mTime2 = lookupTable->GetMTime();
          mTime = (mTime2 > mTime ? mTime2 : mTime);
          }
1521
        }
1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536 1537 1538 1539 1540 1541 1542 1543
      }
    }

  return mTime;
}

//----------------------------------------------------------------------------
double *vtkImageResliceMapper::GetBounds()
{
  // Modify to give just the slice bounds
  if (!this->GetInput())
    {
    vtkMath::UninitializeBounds(this->Bounds);
    return this->Bounds;
    }
  else
    {
    this->UpdateInformation();
    double *spacing = this->DataSpacing;
    double *origin = this->DataOrigin;
    int *extent = this->DataWholeExtent;

1544 1545
    // expand by half a pixel if border is on
    double border = 0.5*(this->Border != 0);