Mercurial > hg > orthanc-stone
diff OrthancStone/Sources/Toolbox/DicomInstanceParameters.cpp @ 1512:244ad1e4e76a
reorganization of folders
author | Sebastien Jodogne <s.jodogne@gmail.com> |
---|---|
date | Tue, 07 Jul 2020 16:21:02 +0200 |
parents | Framework/Toolbox/DicomInstanceParameters.cpp@970ee51fe01f |
children | 85e117739eca |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/OrthancStone/Sources/Toolbox/DicomInstanceParameters.cpp Tue Jul 07 16:21:02 2020 +0200 @@ -0,0 +1,530 @@ +/** + * Stone of Orthanc + * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics + * Department, University Hospital of Liege, Belgium + * Copyright (C) 2017-2020 Osimis S.A., Belgium + * + * This program is free software: you can redistribute it and/or + * modify it under the terms of the GNU Affero General Public License + * as published by the Free Software Foundation, either version 3 of + * the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Affero General Public License for more details. + * + * You should have received a copy of the GNU Affero General Public License + * along with this program. If not, see <http://www.gnu.org/licenses/>. + **/ + + +#include "DicomInstanceParameters.h" + +#include "../Scene2D/ColorTextureSceneLayer.h" +#include "../Scene2D/FloatTextureSceneLayer.h" +#include "../Toolbox/GeometryToolbox.h" +#include "../Toolbox/ImageToolbox.h" + +#include <Images/Image.h> +#include <Images/ImageProcessing.h> +#include <Logging.h> +#include <OrthancException.h> +#include <Toolbox.h> + + +namespace OrthancStone +{ + void DicomInstanceParameters::Data::ComputeDoseOffsets(const Orthanc::DicomMap& dicom) + { + // http://dicom.nema.org/medical/Dicom/2016a/output/chtml/part03/sect_C.8.8.3.2.html + + { + std::string increment; + + if (dicom.LookupStringValue(increment, Orthanc::DICOM_TAG_FRAME_INCREMENT_POINTER, false)) + { + Orthanc::Toolbox::ToUpperCase(increment); + if (increment != "3004,000C") // This is the "Grid Frame Offset Vector" tag (DICOM_TAG_GRID_FRAME_OFFSET_VECTOR) + { + LOG(ERROR) << "RT-DOSE: Bad value for the \"FrameIncrementPointer\" tag"; + return; + } + } + } + + if (!LinearAlgebra::ParseVector(frameOffsets_, dicom, Orthanc::DICOM_TAG_GRID_FRAME_OFFSET_VECTOR) || + frameOffsets_.size() < imageInformation_.GetNumberOfFrames()) + { + LOG(ERROR) << "RT-DOSE: No information about the 3D location of some slice(s)"; + frameOffsets_.clear(); + } + else + { + if (frameOffsets_.size() >= 2) + { + thickness_ = std::abs(frameOffsets_[1] - frameOffsets_[0]); + } + } + } + + + DicomInstanceParameters::Data::Data(const Orthanc::DicomMap& dicom) : + imageInformation_(dicom) + { + if (imageInformation_.GetNumberOfFrames() <= 0) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + + if (!dicom.LookupStringValue(studyInstanceUid_, Orthanc::DICOM_TAG_STUDY_INSTANCE_UID, false) || + !dicom.LookupStringValue(seriesInstanceUid_, Orthanc::DICOM_TAG_SERIES_INSTANCE_UID, false) || + !dicom.LookupStringValue(sopInstanceUid_, Orthanc::DICOM_TAG_SOP_INSTANCE_UID, false)) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + + std::string s; + if (!dicom.LookupStringValue(s, Orthanc::DICOM_TAG_SOP_CLASS_UID, false)) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadFileFormat); + } + else + { + sopClassUid_ = StringToSopClassUid(s); + } + + if (!dicom.ParseDouble(thickness_, Orthanc::DICOM_TAG_SLICE_THICKNESS)) + { + thickness_ = 100.0 * std::numeric_limits<double>::epsilon(); + } + + GeometryToolbox::GetPixelSpacing(pixelSpacingX_, pixelSpacingY_, dicom); + + std::string position, orientation; + if (dicom.LookupStringValue(position, Orthanc::DICOM_TAG_IMAGE_POSITION_PATIENT, false) && + dicom.LookupStringValue(orientation, Orthanc::DICOM_TAG_IMAGE_ORIENTATION_PATIENT, false)) + { + geometry_ = CoordinateSystem3D(position, orientation); + } + + if (sopClassUid_ == SopClassUid_RTDose) + { + ComputeDoseOffsets(dicom); + + static const Orthanc::DicomTag DICOM_TAG_DOSE_UNITS(0x3004, 0x0002); + + if (!dicom.LookupStringValue(doseUnits_, DICOM_TAG_DOSE_UNITS, false)) + { + LOG(ERROR) << "Tag DoseUnits (0x3004, 0x0002) is missing in " << sopInstanceUid_; + doseUnits_ = ""; + } + } + + isColor_ = (imageInformation_.GetPhotometricInterpretation() != Orthanc::PhotometricInterpretation_Monochrome1 && + imageInformation_.GetPhotometricInterpretation() != Orthanc::PhotometricInterpretation_Monochrome2); + + if (dicom.ParseDouble(rescaleIntercept_, Orthanc::DICOM_TAG_RESCALE_INTERCEPT) && + dicom.ParseDouble(rescaleSlope_, Orthanc::DICOM_TAG_RESCALE_SLOPE)) + { + if (sopClassUid_ == SopClassUid_RTDose) + { + LOG(INFO) << "DOSE HAS Rescale*: rescaleIntercept_ = " << rescaleIntercept_ << " rescaleSlope_ = " << rescaleSlope_; + // WE SHOULD NOT TAKE THE RESCALE VALUE INTO ACCOUNT IN THE CASE OF DOSES + hasRescale_ = false; + } + else + { + hasRescale_ = true; + } + + } + else + { + hasRescale_ = false; + } + + if (dicom.ParseDouble(doseGridScaling_, Orthanc::DICOM_TAG_DOSE_GRID_SCALING)) + { + if (sopClassUid_ == SopClassUid_RTDose) + { + LOG(INFO) << "DOSE HAS DoseGridScaling: doseGridScaling_ = " << doseGridScaling_; + } + } + else + { + doseGridScaling_ = 1.0; + if (sopClassUid_ == SopClassUid_RTDose) + { + LOG(ERROR) << "Tag DoseGridScaling (0x3004, 0x000e) is missing in " << sopInstanceUid_ << " doseGridScaling_ will be set to 1.0"; + } + } + + Vector c, w; + if (LinearAlgebra::ParseVector(c, dicom, Orthanc::DICOM_TAG_WINDOW_CENTER) && + LinearAlgebra::ParseVector(w, dicom, Orthanc::DICOM_TAG_WINDOW_WIDTH) && + c.size() > 0 && + w.size() > 0) + { + hasDefaultWindowing_ = true; + defaultWindowingCenter_ = static_cast<float>(c[0]); + defaultWindowingWidth_ = static_cast<float>(w[0]); + } + else + { + hasDefaultWindowing_ = false; + defaultWindowingCenter_ = 0; + defaultWindowingWidth_ = 0; + } + + if (sopClassUid_ == SopClassUid_RTDose) + { + switch (imageInformation_.GetBitsStored()) + { + case 16: + expectedPixelFormat_ = Orthanc::PixelFormat_Grayscale16; + break; + + case 32: + expectedPixelFormat_ = Orthanc::PixelFormat_Grayscale32; + break; + + default: + throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); + } + } + else if (isColor_) + { + expectedPixelFormat_ = Orthanc::PixelFormat_RGB24; + } + else if (imageInformation_.IsSigned()) + { + expectedPixelFormat_ = Orthanc::PixelFormat_SignedGrayscale16; + } + else + { + expectedPixelFormat_ = Orthanc::PixelFormat_Grayscale16; + } + + // This computes the "IndexInSeries" metadata from Orthanc (check + // out "Orthanc::ServerIndex::Store()") + hasIndexInSeries_ = ( + dicom.ParseUnsignedInteger32(indexInSeries_, Orthanc::DICOM_TAG_INSTANCE_NUMBER) || + dicom.ParseUnsignedInteger32(indexInSeries_, Orthanc::DICOM_TAG_IMAGE_INDEX)); + } + + + CoordinateSystem3D DicomInstanceParameters::Data::GetFrameGeometry(unsigned int frame) const + { + if (frame == 0) + { + return geometry_; + } + else if (frame >= imageInformation_.GetNumberOfFrames()) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + else if (sopClassUid_ == SopClassUid_RTDose) + { + if (frame >= frameOffsets_.size()) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_InternalError); + } + + return CoordinateSystem3D( + geometry_.GetOrigin() + frameOffsets_[frame] * geometry_.GetNormal(), + geometry_.GetAxisX(), + geometry_.GetAxisY()); + } + else + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_NotImplemented); + } + } + + + bool DicomInstanceParameters::Data::IsPlaneWithinSlice(unsigned int frame, + const CoordinateSystem3D& plane) const + { + if (frame >= imageInformation_.GetNumberOfFrames()) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_ParameterOutOfRange); + } + + CoordinateSystem3D tmp = geometry_; + + if (frame != 0) + { + tmp = GetFrameGeometry(frame); + } + + double distance; + + return (CoordinateSystem3D::ComputeDistance(distance, tmp, plane) && + distance <= thickness_ / 2.0); + } + + void DicomInstanceParameters::Data::ApplyRescaleAndDoseScaling(Orthanc::ImageAccessor& image, + bool useDouble) const + { + if (image.GetFormat() != Orthanc::PixelFormat_Float32) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageFormat); + } + + double factor = doseGridScaling_; + double offset = 0.0; + + if (hasRescale_) + { + factor *= rescaleSlope_; + offset = rescaleIntercept_; + } + + if ( (factor != 1.0) || (offset != 0.0) ) + { + const unsigned int width = image.GetWidth(); + const unsigned int height = image.GetHeight(); + + for (unsigned int y = 0; y < height; y++) + { + float* p = reinterpret_cast<float*>(image.GetRow(y)); + + if (useDouble) + { + // Slower, accurate implementation using double + for (unsigned int x = 0; x < width; x++, p++) + { + double value = static_cast<double>(*p); + *p = static_cast<float>(value * factor + offset); + } + } + else + { + // Fast, approximate implementation using float + for (unsigned int x = 0; x < width; x++, p++) + { + *p = (*p) * static_cast<float>(factor) + static_cast<float>(offset); + } + } + } + } + } + + double DicomInstanceParameters::GetRescaleIntercept() const + { + if (data_.hasRescale_) + { + return data_.rescaleIntercept_; + } + else + { + LOG(ERROR) << "DicomInstanceParameters::GetRescaleIntercept(): !data_.hasRescale_"; + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + } + + + double DicomInstanceParameters::GetRescaleSlope() const + { + if (data_.hasRescale_) + { + return data_.rescaleSlope_; + } + else + { + LOG(ERROR) << "DicomInstanceParameters::GetRescaleSlope(): !data_.hasRescale_"; + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + } + + + float DicomInstanceParameters::GetDefaultWindowingCenter() const + { + if (data_.hasDefaultWindowing_) + { + return data_.defaultWindowingCenter_; + } + else + { + LOG(ERROR) << "DicomInstanceParameters::GetDefaultWindowingCenter(): no default windowing"; + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + } + + + float DicomInstanceParameters::GetDefaultWindowingWidth() const + { + if (data_.hasDefaultWindowing_) + { + return data_.defaultWindowingWidth_; + } + else + { + LOG(ERROR) << "DicomInstanceParameters::GetDefaultWindowingWidth(): no default windowing"; + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + } + + + Orthanc::ImageAccessor* DicomInstanceParameters::ConvertToFloat(const Orthanc::ImageAccessor& pixelData) const + { + std::unique_ptr<Orthanc::Image> converted(new Orthanc::Image(Orthanc::PixelFormat_Float32, + pixelData.GetWidth(), + pixelData.GetHeight(), + false)); + Orthanc::ImageProcessing::Convert(*converted, pixelData); + + + // Correct rescale slope/intercept if need be + //data_.ApplyRescaleAndDoseScaling(*converted, (pixelData.GetFormat() == Orthanc::PixelFormat_Grayscale32)); + data_.ApplyRescaleAndDoseScaling(*converted, false); + + return converted.release(); + } + + + TextureBaseSceneLayer* DicomInstanceParameters::CreateTexture + (const Orthanc::ImageAccessor& pixelData) const + { + // { + // const Orthanc::ImageAccessor& source = pixelData; + // const void* sourceBuffer = source.GetConstBuffer(); + // intptr_t sourceBufferInt = reinterpret_cast<intptr_t>(sourceBuffer); + // int sourceWidth = source.GetWidth(); + // int sourceHeight = source.GetHeight(); + // int sourcePitch = source.GetPitch(); + + // // TODO: turn error into trace below + // LOG(ERROR) << "ConvertGrayscaleToFloat | source:" + // << " W = " << sourceWidth << " H = " << sourceHeight + // << " P = " << sourcePitch << " B = " << sourceBufferInt + // << " B % 4 == " << sourceBufferInt % 4; + // } + + assert(sizeof(float) == 4); + + Orthanc::PixelFormat sourceFormat = pixelData.GetFormat(); + + if (sourceFormat != GetExpectedPixelFormat()) + { + throw Orthanc::OrthancException(Orthanc::ErrorCode_IncompatibleImageFormat); + } + + if (sourceFormat == Orthanc::PixelFormat_RGB24) + { + // This is the case of a color image. No conversion has to be done. + return new ColorTextureSceneLayer(pixelData); + } + else + { + // This is the case of a grayscale frame. Convert it to Float32. + std::unique_ptr<FloatTextureSceneLayer> texture; + + if (pixelData.GetFormat() == Orthanc::PixelFormat_Float32) + { + texture.reset(new FloatTextureSceneLayer(pixelData)); + } + else + { + std::unique_ptr<Orthanc::ImageAccessor> converted(ConvertToFloat(pixelData)); + texture.reset(new FloatTextureSceneLayer(*converted)); + } + + if (data_.hasDefaultWindowing_) + { + texture->SetCustomWindowing(data_.defaultWindowingCenter_, + data_.defaultWindowingWidth_); + } + + + if (data_.imageInformation_.GetPhotometricInterpretation() + == Orthanc::PhotometricInterpretation_Monochrome1) + { + texture->SetInverted(true); + } + else if (data_.imageInformation_.GetPhotometricInterpretation() + == Orthanc::PhotometricInterpretation_Monochrome2) + { + texture->SetInverted(false); + } + + return texture.release(); + } + } + + + LookupTableTextureSceneLayer* DicomInstanceParameters::CreateLookupTableTexture + (const Orthanc::ImageAccessor& pixelData) const + { + std::unique_ptr<FloatTextureSceneLayer> texture; + + if (pixelData.GetFormat() == Orthanc::PixelFormat_Float32) + { + return new LookupTableTextureSceneLayer(pixelData); + } + else + { + std::unique_ptr<Orthanc::ImageAccessor> converted(ConvertToFloat(pixelData)); + return new LookupTableTextureSceneLayer(*converted); + } + } + + + unsigned int DicomInstanceParameters::GetIndexInSeries() const + { + if (data_.hasIndexInSeries_) + { + return data_.indexInSeries_; + } + else + { + LOG(ERROR) << "DicomInstanceParameters::GetIndexInSeries(): !data_.hasIndexInSeries_"; + throw Orthanc::OrthancException(Orthanc::ErrorCode_BadSequenceOfCalls); + } + } + + + double DicomInstanceParameters::Data::ApplyRescale(double value) const + { + double factor = doseGridScaling_; + double offset = 0.0; + + if (hasRescale_) + { + factor *= rescaleSlope_; + offset = rescaleIntercept_; + } + + return (value * factor + offset); + } + + + bool DicomInstanceParameters::Data::ComputeRegularSpacing(double& spacing) const + { + if (frameOffsets_.size() == 0) // Not a RT-DOSE + { + return false; + } + else if (frameOffsets_.size() == 1) + { + spacing = 1; // Edge case: RT-DOSE with one single frame + return true; + } + else + { + spacing = std::abs(frameOffsets_[1] - frameOffsets_[0]); + + for (size_t i = 1; i + 1 < frameOffsets_.size(); i++) + { + double s = frameOffsets_[i + 1] - frameOffsets_[i]; + if (!LinearAlgebra::IsNear(spacing, s, 0.001)) + { + return false; + } + } + + return true; + } + } +}