diff --git a/Modules/PhotoacousticsLib/include/mitkPAVolumeManipulator.h b/Modules/PhotoacousticsLib/include/mitkPAVolumeManipulator.h index 4b1eada957..3a02d50974 100644 --- a/Modules/PhotoacousticsLib/include/mitkPAVolumeManipulator.h +++ b/Modules/PhotoacousticsLib/include/mitkPAVolumeManipulator.h @@ -1,58 +1,58 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITKPHOTOACOUSTIC3DVOLUMEMANIPULATOR_H #define MITKPHOTOACOUSTIC3DVOLUMEMANIPULATOR_H #include #include #include "mitkPAVolume.h" #include "mitkPAInSilicoTissueVolume.h" namespace mitk { namespace pa { class MITKPHOTOACOUSTICSLIB_EXPORT VolumeManipulator final { public: /** * @brief ThresholdImage applies a binary threshold filter to this image. * @param threshold */ static void ThresholdImage(Volume::Pointer image, double threshold); /** * @brief Multiplies the image with a given factor * @param factor */ static void MultiplyImage(Volume::Pointer image, double factor); static void GaussianBlur3D(Volume::Pointer paVolume, double sigma); static void Log10Image(Volume::Pointer image); static void RescaleImage(InSilicoTissueVolume::Pointer image, double ratio); - static Volume::Pointer RescaleImage(Volume::Pointer image, double ratio); + static Volume::Pointer RescaleImage(Volume::Pointer image, double ratio, double sigma); private: VolumeManipulator(); virtual ~VolumeManipulator(); }; } } #endif // MITKPHOTOACOUSTIC3DVOLUMEMANIPULATOR_H diff --git a/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp b/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp index 7e9090acb7..402a2e02ee 100644 --- a/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp +++ b/Modules/PhotoacousticsLib/src/Domain/Vessel/mitkPAVessel.cpp @@ -1,117 +1,117 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPAVessel.h" #include #define _USE_MATH_DEFINES #include #include mitk::pa::Vessel::Vessel(VesselProperties::Pointer initialProperties) : m_RangeDistribution(M_PI / 16, M_PI / 8), m_SignDistribution(-1, 1) { m_Finished = false; //Copy this so it may be reused for other vessels. m_VesselProperties = VesselProperties::New(initialProperties); m_RadiusRangeDistribution = std::uniform_real_distribution<>(NEW_RADIUS_MINIMUM_RELATIVE_SIZE, NEW_RADIUS_MAXIMUM_RELATIVE_SIZE); m_VesselMeanderStrategy = VesselMeanderStrategy::New(); m_WalkedDistance = 0; m_VesselDrawer = VesselDrawer::New(); } mitk::pa::Vessel::~Vessel() { m_VesselProperties = nullptr; m_VesselMeanderStrategy = nullptr; } void mitk::pa::Vessel::ExpandVessel(InSilicoTissueVolume::Pointer volume, CalculateNewVesselPositionCallback calculateNewPosition, double bendingFactor, std::mt19937* rng) { m_VesselDrawer->DrawVesselInVolume(m_VesselProperties, volume); (m_VesselMeanderStrategy->*calculateNewPosition)(m_VesselProperties->GetPositionVector(), m_VesselProperties->GetDirectionVector(), bendingFactor, rng); - m_WalkedDistance += m_VesselProperties->GetDirectionVector()->GetNorm(); + m_WalkedDistance += (m_VesselProperties->GetDirectionVector()->GetNorm() / volume->GetSpacing()); } bool mitk::pa::Vessel::CanBifurcate() { return m_VesselProperties->GetBifurcationFrequency() < m_WalkedDistance; } int mitk::pa::Vessel::GetSign(std::mt19937 *rng) { if (m_SignDistribution(*rng) < 0) return -1; return 1; } mitk::pa::Vessel::Pointer mitk::pa::Vessel::Bifurcate(std::mt19937* rng) { VesselProperties::Pointer vesselParams = VesselProperties::New(m_VesselProperties); double thetaChange = m_RangeDistribution(*rng) * GetSign(rng); double phiChange = m_RangeDistribution(*rng) * GetSign(rng); vesselParams->GetDirectionVector()->Rotate(thetaChange, phiChange); m_VesselProperties->GetDirectionVector()->Rotate(-thetaChange, -phiChange); double newRadius = m_RadiusRangeDistribution(*rng)*m_VesselProperties->GetRadiusInVoxel(); vesselParams->SetRadiusInVoxel(newRadius); m_VesselProperties->SetRadiusInVoxel( sqrt(m_VesselProperties->GetRadiusInVoxel()*m_VesselProperties->GetRadiusInVoxel() - newRadius*newRadius)); m_WalkedDistance = 0; return Vessel::New(vesselParams); } bool mitk::pa::Vessel::IsFinished() { return m_VesselProperties->GetRadiusInVoxel() < MINIMUM_VESSEL_RADIUS; } bool mitk::pa::Equal(const Vessel::Pointer leftHandSide, const Vessel::Pointer rightHandSide, double eps, bool verbose) { MITK_INFO(verbose) << "=== mitk::pa::Vessel Equal ==="; if (rightHandSide.IsNull() || leftHandSide.IsNull()) { MITK_INFO(verbose) << "Cannot compare nullpointers"; return false; } if (leftHandSide->IsFinished() != rightHandSide->IsFinished()) { MITK_INFO(verbose) << "Not same finished state."; return false; } if (leftHandSide->CanBifurcate() != rightHandSide->CanBifurcate()) { MITK_INFO(verbose) << "Not same bifurcation state."; return false; } if (!Equal(leftHandSide->GetVesselProperties(), rightHandSide->GetVesselProperties(), eps, verbose)) { MITK_INFO(verbose) << "Vesselproperties not equal"; return false; } return true; } diff --git a/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp b/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp index 41ce37435e..e618155824 100644 --- a/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp +++ b/Modules/PhotoacousticsLib/src/Generator/mitkPATissueGenerator.cpp @@ -1,176 +1,178 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPATissueGenerator.h" #include "mitkPAVector.h" #include "mitkPAVolumeManipulator.h" mitk::pa::InSilicoTissueVolume::Pointer mitk::pa::InSilicoTissueGenerator::GenerateInSilicoData( TissueGeneratorParameters::Pointer parameters) { MITK_DEBUG << "Initializing Empty Volume"; + const double RESAMPLING_FACTOR = 2; + if (parameters->GetDoPartialVolume()) { - parameters->SetXDim(parameters->GetXDim() * 10); - parameters->SetYDim(parameters->GetYDim() * 10); - parameters->SetZDim(parameters->GetZDim() * 10); - parameters->SetVesselBifurcationFrequency(parameters->GetVesselBifurcationFrequency() * 10); - parameters->SetVoxelSpacingInCentimeters(parameters->GetVoxelSpacingInCentimeters() / 10); + parameters->SetXDim(parameters->GetXDim() * RESAMPLING_FACTOR); + parameters->SetYDim(parameters->GetYDim() * RESAMPLING_FACTOR); + parameters->SetZDim(parameters->GetZDim() * RESAMPLING_FACTOR); + parameters->SetVesselBifurcationFrequency(parameters->GetVesselBifurcationFrequency() * RESAMPLING_FACTOR); + parameters->SetVoxelSpacingInCentimeters(parameters->GetVoxelSpacingInCentimeters() / RESAMPLING_FACTOR); } auto generatedVolume = mitk::pa::InSilicoTissueVolume::New(parameters); const double DIRECTION_VECTOR_INITIAL_VARIANCE = 0.2; std::mt19937 randomNumberGenerator; std::random_device randomDevice; if (parameters->GetUseRngSeed()) { randomNumberGenerator.seed(parameters->GetRngSeed()); } else { if (randomDevice.entropy() > 0.1) { randomNumberGenerator.seed(randomDevice()); } else { randomNumberGenerator.seed(std::chrono::duration_cast(std::chrono::high_resolution_clock::now().time_since_epoch()).count()); } } std::uniform_int_distribution randomNumVesselDistribution(parameters->GetMinNumberOfVessels(), parameters->GetMaxNumberOfVessels()); std::uniform_real_distribution randomBendingDistribution(parameters->GetMinVesselBending(), parameters->GetMaxVesselBending()); std::uniform_real_distribution randomAbsorptionDistribution(parameters->GetMinVesselAbsorption(), parameters->GetMaxVesselAbsorption()); std::uniform_real_distribution randomRadiusDistribution(parameters->GetMinVesselRadiusInMillimeters(), parameters->GetMaxVesselRadiusInMillimeters()); std::uniform_real_distribution randomScatteringDistribution(parameters->GetMinVesselScattering(), parameters->GetMaxVesselScattering()); std::uniform_real_distribution randomAnisotropyDistribution(parameters->GetMinVesselAnisotropy(), parameters->GetMaxVesselAnisotropy()); std::uniform_int_distribution borderTypeDistribution(0, 3); int numberOfBloodVessels = randomNumVesselDistribution(randomNumberGenerator); generatedVolume->AddIntProperty("numberOfVessels", numberOfBloodVessels); generatedVolume->AddIntProperty("bifurcationFrequency", parameters->GetVesselBifurcationFrequency()); MITK_INFO << "Simulating " << numberOfBloodVessels << " vessel structures"; for (int vesselNumber = 0; vesselNumber < numberOfBloodVessels; vesselNumber++) { Vector::Pointer initialPosition = Vector::New(); Vector::Pointer initialDirection = Vector::New(); double initialRadius = randomRadiusDistribution(randomNumberGenerator) / parameters->GetVoxelSpacingInCentimeters() / 10; std::stringstream radiusString; radiusString << "vessel_" << vesselNumber + 1 << "_radius"; generatedVolume->AddDoubleProperty(radiusString.str(), initialRadius); double absorptionCoefficient = randomAbsorptionDistribution(randomNumberGenerator); std::stringstream absorptionString; absorptionString << "vessel_" << vesselNumber + 1 << "_absorption"; generatedVolume->AddDoubleProperty(absorptionString.str(), absorptionCoefficient); double bendingFactor = randomBendingDistribution(randomNumberGenerator); std::stringstream bendingString; bendingString << "vessel_" << vesselNumber + 1 << "_bendingFactor"; generatedVolume->AddDoubleProperty(bendingString.str(), bendingFactor); double vesselScattering = randomScatteringDistribution(randomNumberGenerator); std::stringstream scatteringString; scatteringString << "vessel_" << vesselNumber + 1 << "_scattering"; generatedVolume->AddDoubleProperty(scatteringString.str(), vesselScattering); double vesselAnisotropy = randomAnisotropyDistribution(randomNumberGenerator); std::stringstream anisotropyString; anisotropyString << "vessel_" << vesselNumber + 1 << "_anisotropy"; generatedVolume->AddDoubleProperty(anisotropyString.str(), vesselAnisotropy); /*The vessel tree shall start at one of the 4 sides of the volume. * The vessels will always be completely contained in the volume * when starting to meander. * They will meander in a direction perpendicular to the * plane they started from (within the limits of the * DIRECTION_VECTOR_INITIAL_VARIANCE) */ int borderType = borderTypeDistribution(randomNumberGenerator); switch (borderType) { case 0: initialPosition->Randomize(0, 0, initialRadius, parameters->GetYDim() - initialRadius, parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &randomNumberGenerator); initialDirection->Randomize(1, 2, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator); break; case 1: initialPosition->Randomize(parameters->GetXDim(), parameters->GetXDim(), initialRadius, parameters->GetYDim() - initialRadius, parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &randomNumberGenerator); initialDirection->Randomize(-2, -1, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator); break; case 2: initialPosition->Randomize(initialRadius, parameters->GetXDim() - initialRadius, 0, 0, parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &randomNumberGenerator); initialDirection->Randomize(-DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, 1, 2, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator); break; case 3: initialPosition->Randomize(initialRadius, parameters->GetXDim() - initialRadius, parameters->GetYDim(), parameters->GetYDim(), parameters->GetMinVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), parameters->GetMaxVesselZOrigin() / parameters->GetVoxelSpacingInCentimeters(), &randomNumberGenerator); initialDirection->Randomize(-DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, -2, -1, -DIRECTION_VECTOR_INITIAL_VARIANCE, DIRECTION_VECTOR_INITIAL_VARIANCE, &randomNumberGenerator); break; } VesselProperties::Pointer vesselParams = VesselProperties::New(); vesselParams->SetDirectionVector(initialDirection); vesselParams->SetPositionVector(initialPosition); vesselParams->SetRadiusInVoxel(initialRadius); vesselParams->SetAbsorptionCoefficient(absorptionCoefficient); vesselParams->SetScatteringCoefficient(vesselScattering); vesselParams->SetAnisotopyCoefficient(vesselAnisotropy); vesselParams->SetBifurcationFrequency(parameters->GetVesselBifurcationFrequency()); vesselParams->SetDoPartialVolume(parameters->GetDoPartialVolume()); VesselTree::Pointer vesselTree = VesselTree::New(vesselParams); while (!vesselTree->IsFinished()) { vesselTree->Step(generatedVolume, parameters->GetCalculateNewVesselPositionCallback(), bendingFactor, &randomNumberGenerator); } } if (parameters->GetDoPartialVolume()) { - VolumeManipulator::RescaleImage(generatedVolume, 0.1); - parameters->SetXDim(parameters->GetXDim() / 10); - parameters->SetYDim(parameters->GetYDim() / 10); - parameters->SetZDim(parameters->GetZDim() / 10); - parameters->SetVesselBifurcationFrequency(parameters->GetVesselBifurcationFrequency() / 10); - parameters->SetVoxelSpacingInCentimeters(parameters->GetVoxelSpacingInCentimeters() * 10); + VolumeManipulator::RescaleImage(generatedVolume, (1.0 / RESAMPLING_FACTOR)); + parameters->SetXDim(parameters->GetXDim() / RESAMPLING_FACTOR); + parameters->SetYDim(parameters->GetYDim() / RESAMPLING_FACTOR); + parameters->SetZDim(parameters->GetZDim() / RESAMPLING_FACTOR); + parameters->SetVesselBifurcationFrequency(parameters->GetVesselBifurcationFrequency() / RESAMPLING_FACTOR); + parameters->SetVoxelSpacingInCentimeters(parameters->GetVoxelSpacingInCentimeters() * RESAMPLING_FACTOR); } generatedVolume->FinalizeVolume(); return generatedVolume; } mitk::pa::InSilicoTissueGenerator::InSilicoTissueGenerator() { } mitk::pa::InSilicoTissueGenerator::~InSilicoTissueGenerator() { } diff --git a/Modules/PhotoacousticsLib/src/Utils/mitkPAVolumeManipulator.cpp b/Modules/PhotoacousticsLib/src/Utils/mitkPAVolumeManipulator.cpp index 9a2abb9380..c5f728505e 100644 --- a/Modules/PhotoacousticsLib/src/Utils/mitkPAVolumeManipulator.cpp +++ b/Modules/PhotoacousticsLib/src/Utils/mitkPAVolumeManipulator.cpp @@ -1,225 +1,236 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #include "mitkPAVolumeManipulator.h" #include "mitkPAExceptions.h" #include "mitkPAInSilicoTissueVolume.h" #include "itkResampleImageFilter.h" #include "itkGaussianInterpolateImageFunction.h" // Includes for image casting between ITK and MITK #include #include #include mitk::pa::VolumeManipulator::VolumeManipulator() {} mitk::pa::VolumeManipulator::~VolumeManipulator() {} void mitk::pa::VolumeManipulator::ThresholdImage(Volume::Pointer image, double threshold) { for (unsigned int z = 0; z < image->GetZDim(); z++) for (unsigned int y = 0; y < image->GetYDim(); y++) for (unsigned int x = 0; x < image->GetXDim(); x++) if (image->GetData(x, y, z) > threshold) image->SetData(1, x, y, z); else image->SetData(0, x, y, z); } void mitk::pa::VolumeManipulator::MultiplyImage(Volume::Pointer image, double factor) { for (unsigned int z = 0; z < image->GetZDim(); z++) for (unsigned int y = 0; y < image->GetYDim(); y++) for (unsigned int x = 0; x < image->GetXDim(); x++) image->SetData(image->GetData(x, y, z)*factor, x, y, z); } void mitk::pa::VolumeManipulator::Log10Image(Volume::Pointer image) { for (unsigned int z = 0; z < image->GetZDim(); z++) for (unsigned int y = 0; y < image->GetYDim(); y++) for (unsigned int x = 0; x < image->GetXDim(); x++) { if (image->GetData(x, y, z) < 0) { MITK_ERROR << "Signal was smaller than 0. There is an error in the input image file."; throw InvalidValueException("Signal was smaller than 0. There is an error in the input image file."); } image->SetData(log10(image->GetData(x, y, z)), x, y, z); } } void mitk::pa::VolumeManipulator::RescaleImage(InSilicoTissueVolume::Pointer image, double ratio) { - image->SetScatteringVolume(RescaleImage(image->GetScatteringVolume(), ratio)); - image->SetAnisotropyVolume(RescaleImage(image->GetAnisotropyVolume(), ratio)); - image->SetSegmentationVolume(RescaleImage(image->GetSegmentationVolume(), ratio)); - image->SetAbsorptionVolume(RescaleImage(image->GetAbsorptionVolume(), ratio)); + MITK_INFO << "Rescaling images.."; + double sigma = image->GetSpacing() / (ratio * 2); + image->SetAbsorptionVolume(RescaleImage(image->GetAbsorptionVolume(), ratio, sigma)); + image->SetScatteringVolume(RescaleImage(image->GetScatteringVolume(), ratio, sigma)); + image->SetAnisotropyVolume(RescaleImage(image->GetAnisotropyVolume(), ratio, sigma)); + image->SetSegmentationVolume(RescaleImage(image->GetSegmentationVolume(), ratio, 0)); + MITK_INFO << "Rescaling images..[Done]"; } -mitk::pa::Volume::Pointer mitk::pa::VolumeManipulator::RescaleImage(Volume::Pointer image, double ratio) +mitk::pa::Volume::Pointer mitk::pa::VolumeManipulator::RescaleImage(Volume::Pointer image, double ratio, double sigma) { + MITK_INFO << "Rescaling image.."; typedef itk::Image ImageType; typedef itk::ResampleImageFilter FilterType; typedef itk::GaussianInterpolateImageFunction InterpolatorType; auto input = image->AsMitkImage(); ImageType::Pointer itkInput = ImageType::New(); mitk::CastToItkImage(input, itkInput); ImageType::SizeType outputSize; outputSize[0] = input->GetDimensions()[0] * ratio; outputSize[1] = input->GetDimensions()[1] * ratio; outputSize[2] = input->GetDimensions()[2] * ratio; FilterType::Pointer resampleImageFilter = FilterType::New(); resampleImageFilter->SetInput(itkInput); resampleImageFilter->SetSize(outputSize); - resampleImageFilter->SetInterpolator(InterpolatorType::New()); + if (sigma > mitk::eps) + { + auto interpolator = InterpolatorType::New(); + interpolator->SetSigma(sigma); + resampleImageFilter->SetInterpolator(interpolator); + } resampleImageFilter->SetOutputSpacing(input->GetGeometry()->GetSpacing()[0] / ratio); + MITK_INFO << "Update.."; resampleImageFilter->UpdateLargestPossibleRegion(); + MITK_INFO << "Update..[Done]"; ImageType::Pointer output = resampleImageFilter->GetOutput(); mitk::Image::Pointer mitkOutput = mitk::Image::New(); GrabItkImageMemory(output, mitkOutput); - + MITK_INFO << "Rescaling image..[Done]"; return Volume::New(mitkOutput); } /** * @brief Fast 3D Gaussian convolution IIR approximation * @param paVolume * @param sigma * @author Pascal Getreuer * * Copyright (c) 2011, Pascal Getreuer * All rights reserved. * * This program is free software: you can redistribute it and/or modify it * under the terms of the simplified BSD license. * * You should have received a copy of these licenses along with this program. * If not, see . */ void mitk::pa::VolumeManipulator::GaussianBlur3D(mitk::pa::Volume::Pointer paVolume, double sigma) { double* volume = paVolume->GetData(); long width = paVolume->GetYDim(); long height = paVolume->GetXDim(); long depth = paVolume->GetZDim(); const long plane = width*height; const long numel = plane*depth; double lambda, dnu; double nu, boundaryscale, postscale; double *ptr; long i, x, y, z; int step; if (sigma <= 0) return; lambda = (sigma*sigma) / (8.0); dnu = (1.0 + 2.0*lambda - sqrt(1.0 + 4.0*lambda)) / (2.0*lambda); nu = dnu; boundaryscale = 1.0 / (1.0 - dnu); postscale = pow(dnu / lambda, 12); /* Filter horizontally along each row */ for (z = 0; z < depth; z++) { for (y = 0; y < height; y++) { for (step = 0; step < 4; step++) { ptr = volume + width*(y + height*z); ptr[0] *= boundaryscale; /* Filter rightwards */ for (x = 1; x < width; x++) { ptr[x] += nu*ptr[x - 1]; } ptr[x = width - 1] *= boundaryscale; /* Filter leftwards */ for (; x > 0; x--) { ptr[x - 1] += nu*ptr[x]; } } } } /* Filter vertically along each column */ for (z = 0; z < depth; z++) { for (x = 0; x < width; x++) { for (step = 0; step < 4; step++) { ptr = volume + x + plane*z; ptr[0] *= boundaryscale; /* Filter downwards */ for (i = width; i < plane; i += width) { ptr[i] += nu*ptr[i - width]; } ptr[i = plane - width] *= boundaryscale; /* Filter upwards */ for (; i > 0; i -= width) { ptr[i - width] += nu*ptr[i]; } } } } /* Filter along z-dimension */ for (y = 0; y < height; y++) { for (x = 0; x < width; x++) { for (step = 0; step < 4; step++) { ptr = volume + x + width*y; ptr[0] *= boundaryscale; for (i = plane; i < numel; i += plane) { ptr[i] += nu*ptr[i - plane]; } ptr[i = numel - plane] *= boundaryscale; for (; i > 0; i -= plane) { ptr[i - plane] += nu*ptr[i]; } } } } for (i = 0; i < numel; i++) { volume[i] *= postscale; } }