C++ Pointer::GetVolumeData方法代码示例

mitk::DataNode::Pointer mitk::Tool::CreateEmptySegmentationNode( Image* original, const std::string& organName, const mitk::Color& color )
{
  // we NEED a reference image for size etc.
  if (!original) return NULL;

  // actually create a new empty segmentation
  PixelType pixelType(mitk::MakeScalarPixelType<DefaultSegmentationDataType>() );
  Image::Pointer segmentation = Image::New();

  if (original->GetDimension() == 2)
  {
    const unsigned int dimensions[] = { original->GetDimension(0), original->GetDimension(1), 1 };
    segmentation->Initialize(pixelType, 3, dimensions);
  }
  else
  {
    segmentation->Initialize(pixelType, original->GetDimension(), original->GetDimensions());
  }

  unsigned int byteSize = sizeof(DefaultSegmentationDataType);

  if(segmentation->GetDimension() < 4)
  {
    for (unsigned int dim = 0; dim < segmentation->GetDimension(); ++dim)
    {
      byteSize *= segmentation->GetDimension(dim);
    }

    mitk::ImageWriteAccessor writeAccess(segmentation, segmentation->GetVolumeData(0));

    memset( writeAccess.GetData(), 0, byteSize );
  }
  else
  {//if we have a time-resolved image we need to set memory to 0 for each time step
    for (unsigned int dim = 0; dim < 3; ++dim)
    {
      byteSize *= segmentation->GetDimension(dim);
    }

    for( unsigned int volumeNumber = 0; volumeNumber < segmentation->GetDimension(3); volumeNumber++)
    {
      mitk::ImageWriteAccessor writeAccess(segmentation, segmentation->GetVolumeData(volumeNumber));

      memset( writeAccess.GetData(), 0, byteSize );
    }
  }

  if (original->GetTimeGeometry() )
  {
    TimeGeometry::Pointer originalGeometry = original->GetTimeGeometry()->Clone();
    segmentation->SetTimeGeometry( originalGeometry );
  }
  else
  {
    Tool::ErrorMessage("Original image does not have a 'Time sliced geometry'! Cannot create a segmentation.");
    return NULL;
  }

  return CreateSegmentationNode( segmentation, organName, color );
}

mitk::Image::Pointer mitk::CompressedImageContainer::GetImage()
{
  if (m_ByteBuffers.empty())
    return nullptr;

  // uncompress image data, create an Image
  Image::Pointer image = Image::New();
  unsigned int dims[20]; // more than 20 dimensions and bang
  for (unsigned int dim = 0; dim < m_ImageDimension; ++dim)
    dims[dim] = m_ImageDimensions[dim];

  image->Initialize(*m_PixelType, m_ImageDimension, dims); // this IS needed, right ?? But it does allocate memory ->
                                                           // does create one big lump of memory (also in windows)

  unsigned int timeStep(0);
  for (auto iter = m_ByteBuffers.begin(); iter != m_ByteBuffers.end(); ++iter, ++timeStep)
  {
    ImageReadAccessor imgAcc(image, image->GetVolumeData(timeStep));
    auto *dest((unsigned char *)imgAcc.GetData());
    ::uLongf destLen(m_OneTimeStepImageSizeInBytes);
    ::Bytef *source(iter->first);
    ::uLongf sourceLen(iter->second);
    int zlibRetVal = ::uncompress(dest, &destLen, source, sourceLen);
    if (itk::Object::GetDebug())
    {
      if (zlibRetVal == Z_OK)
      {
        MITK_INFO << "Success, destLen now " << destLen << " bytes" << std::endl;
      }
      else
      {
        switch (zlibRetVal)
        {
          case Z_DATA_ERROR:
            MITK_ERROR << "compressed data corrupted" << std::endl;
            break;
          case Z_MEM_ERROR:
            MITK_ERROR << "not enough memory" << std::endl;
            break;
          case Z_BUF_ERROR:
            MITK_ERROR << "output buffer too small" << std::endl;
            break;
          default:
            MITK_ERROR << "other, unspecified error" << std::endl;
            break;
        }
      }
    }
  }

  image->SetGeometry(m_ImageGeometry);
  image->Modified();

  return image;
}

static void ReadPixel(const PixelType&, Image::Pointer image, const itk::Index<3>& index, ScalarType* returnValue)
{
  switch (image->GetDimension())
  {
  case 2:
    {
      ImagePixelReadAccessor<T, 2> readAccess(image, image->GetSliceData(0));
      *returnValue = readAccess.GetPixelByIndex(reinterpret_cast<const itk::Index<2>&>(index));
      break;
    }

  case 3:
    {
      ImagePixelReadAccessor<T, 3> readAccess(image, image->GetVolumeData(0));
      *returnValue = readAccess.GetPixelByIndex(index);
      break;
    }

  default:
    *returnValue = 0;
    break;
  }
}

void
mitk::TimeFramesRegistrationHelper::Generate()
{
  CheckValidInputs();

  //prepare processing
  mitk::Image::Pointer targetFrame = GetFrameImage(this->m_4DImage, 0);

  this->m_Registered4DImage = this->m_4DImage->Clone();

  Image::ConstPointer mask;

  if (m_TargetMask.IsNotNull())
  {
    if (m_TargetMask->GetTimeSteps() > 1)
    {
      mask = GetFrameImage(m_TargetMask, 0);
    }
    else
    {
      mask = m_TargetMask;
    }
  }

  double progressDelta = 1.0 / ((this->m_4DImage->GetTimeSteps() - 1) * 3.0);
  m_Progress = 0.0;

  //process the frames
  for (unsigned int i = 1; i < this->m_4DImage->GetTimeSteps(); ++i)
  {
    Image::Pointer movingFrame = GetFrameImage(this->m_4DImage, i);
    Image::Pointer mappedFrame;

    IgnoreListType::iterator finding = std::find(m_IgnoreList.begin(), m_IgnoreList.end(), i);


    if (finding == m_IgnoreList.end())
    {
      //frame should be processed
      RegistrationPointer reg = DoFrameRegistration(movingFrame, targetFrame, mask);

      m_Progress += progressDelta;
      this->InvokeEvent(::mitk::FrameRegistrationEvent(0,
                        "Registred frame #" +::map::core::convert::toStr(i)));

      mappedFrame = DoFrameMapping(movingFrame, reg, targetFrame);

      m_Progress += progressDelta;
      this->InvokeEvent(::mitk::FrameMappingEvent(0,
                        "Mapped frame #" + ::map::core::convert::toStr(i)));

      mitk::ImageReadAccessor accessor(mappedFrame, mappedFrame->GetVolumeData(0, 0, nullptr,
                                       mitk::Image::ReferenceMemory));


      this->m_Registered4DImage->SetVolume(accessor.GetData(), i);
      this->m_Registered4DImage->GetTimeGeometry()->SetTimeStepGeometry(mappedFrame->GetGeometry(), i);

      m_Progress += progressDelta;
    }
    else
    {
      m_Progress += 3 * progressDelta;
    }

    this->InvokeEvent(::itk::ProgressEvent());

  }

};

mitk::GIFVolumetricStatistics::FeatureListType mitk::GIFVolumetricStatistics::CalculateFeatures(const Image::Pointer & image, const Image::Pointer &mask)
{
  FeatureListType featureList;
  if (image->GetDimension() < 3)
  {
    return featureList;
  }


  AccessByItk_3(image, CalculateVolumeStatistic, mask, featureList, FeatureDescriptionPrefix());
  AccessByItk_3(mask, CalculateLargestDiameter, image, featureList, FeatureDescriptionPrefix());

  vtkSmartPointer<vtkImageMarchingCubes> mesher = vtkSmartPointer<vtkImageMarchingCubes>::New();
  vtkSmartPointer<vtkMassProperties> stats = vtkSmartPointer<vtkMassProperties>::New();
  mesher->SetInputData(mask->GetVtkImageData());
  mesher->SetValue(0, 0.5);
  stats->SetInputConnection(mesher->GetOutputPort());
  stats->Update();

  double pi = vnl_math::pi;

  double meshVolume = stats->GetVolume();
  double meshSurf = stats->GetSurfaceArea();
  double pixelVolume = featureList[1].second;
  double pixelSurface = featureList[3].second;

  MITK_INFO << "Surface: " << pixelSurface << " Volume: " << pixelVolume;

  double compactness1 = pixelVolume / (std::sqrt(pi) * std::pow(meshSurf, 2.0 / 3.0));
  double compactness1Pixel = pixelVolume / (std::sqrt(pi) * std::pow(pixelSurface, 2.0 / 3.0));
  //This is the definition used by Aertz. However, due to 2/3 this feature is not demensionless. Use compactness3 instead.

  double compactness2 = 36 * pi*pixelVolume*pixelVolume / meshSurf / meshSurf / meshSurf;
  double compactness2MeshMesh = 36 * pi*meshVolume*meshVolume / meshSurf / meshSurf / meshSurf;
  double compactness2Pixel = 36 * pi*pixelVolume*pixelVolume / pixelSurface / pixelSurface / pixelSurface;
  double compactness3 = pixelVolume / (std::sqrt(pi) * std::pow(meshSurf, 3.0 / 2.0));
  double compactness3MeshMesh = meshVolume / (std::sqrt(pi) * std::pow(meshSurf, 3.0 / 2.0));
  double compactness3Pixel = pixelVolume / (std::sqrt(pi) * std::pow(pixelSurface, 3.0 / 2.0));

  double sphericity = std::pow(pi, 1 / 3.0) *std::pow(6 * pixelVolume, 2.0 / 3.0) / meshSurf;
  double sphericityMesh = std::pow(pi, 1 / 3.0) *std::pow(6 * meshVolume, 2.0 / 3.0) / meshSurf;
  double sphericityPixel = std::pow(pi, 1 / 3.0) *std::pow(6 * pixelVolume, 2.0 / 3.0) / pixelSurface;
  double surfaceToVolume = meshSurf / meshVolume;
  double surfaceToVolumePixel = pixelSurface / pixelVolume;
  double sphericalDisproportion = meshSurf / 4 / pi / std::pow(3.0 / 4.0 / pi * pixelVolume, 2.0 / 3.0);
  double sphericalDisproportionMesh = meshSurf / 4 / pi / std::pow(3.0 / 4.0 / pi * meshVolume, 2.0 / 3.0);
  double sphericalDisproportionPixel = pixelSurface / 4 / pi / std::pow(3.0 / 4.0 / pi * pixelVolume, 2.0 / 3.0);
  double asphericity = std::pow(1.0/compactness2, (1.0 / 3.0)) - 1;
  double asphericityMesh = std::pow(1.0 / compactness2MeshMesh, (1.0 / 3.0)) - 1;
  double asphericityPixel = std::pow(1.0/compactness2Pixel, (1.0 / 3.0)) - 1;

  //Calculate center of mass shift
  int xx = mask->GetDimensions()[0];
  int yy = mask->GetDimensions()[1];
  int zz = mask->GetDimensions()[2];

  double xd = mask->GetGeometry()->GetSpacing()[0];
  double yd = mask->GetGeometry()->GetSpacing()[1];
  double zd = mask->GetGeometry()->GetSpacing()[2];

  vtkSmartPointer<vtkDoubleArray> dataset1Arr = vtkSmartPointer<vtkDoubleArray>::New();
  vtkSmartPointer<vtkDoubleArray> dataset2Arr = vtkSmartPointer<vtkDoubleArray>::New();
  vtkSmartPointer<vtkDoubleArray> dataset3Arr = vtkSmartPointer<vtkDoubleArray>::New();
  dataset1Arr->SetNumberOfComponents(1);
  dataset2Arr->SetNumberOfComponents(1);
  dataset3Arr->SetNumberOfComponents(1);
  dataset1Arr->SetName("M1");
  dataset2Arr->SetName("M2");
  dataset3Arr->SetName("M3");

  vtkSmartPointer<vtkDoubleArray> dataset1ArrU = vtkSmartPointer<vtkDoubleArray>::New();
  vtkSmartPointer<vtkDoubleArray> dataset2ArrU = vtkSmartPointer<vtkDoubleArray>::New();
  vtkSmartPointer<vtkDoubleArray> dataset3ArrU = vtkSmartPointer<vtkDoubleArray>::New();
  dataset1ArrU->SetNumberOfComponents(1);
  dataset2ArrU->SetNumberOfComponents(1);
  dataset3ArrU->SetNumberOfComponents(1);
  dataset1ArrU->SetName("M1");
  dataset2ArrU->SetName("M2");
  dataset3ArrU->SetName("M3");

  for (int x = 0; x < xx; x++)
  {
    for (int y = 0; y < yy; y++)
    {
      for (int z = 0; z < zz; z++)
      {
        itk::Image<int,3>::IndexType index;

        index[0] = x;
        index[1] = y;
        index[2] = z;

        mitk::ScalarType pxImage;
        mitk::ScalarType pxMask;

        mitkPixelTypeMultiplex5(
              mitk::FastSinglePixelAccess,
              image->GetChannelDescriptor().GetPixelType(),
              image,
              image->GetVolumeData(),
//.........这里部分代码省略.........