Mercurial > hg > orthanc-stone
changeset 126:c9e88e7935a4 wasm
rt-struct projection
| author | Sebastien Jodogne <s.jodogne@gmail.com> |
|---|---|
| date | Wed, 08 Nov 2017 17:27:38 +0100 |
| parents | 44fc253d4876 |
| children | b963f3a9a86c |
| files | Applications/Samples/SingleVolumeApplication.h Framework/Layers/DicomStructureSetRendererFactory.cpp Framework/Toolbox/DicomStructureSet.cpp Framework/Toolbox/DicomStructureSet.h |
| diffstat | 4 files changed, 140 insertions(+), 4 deletions(-) [+] |
line wrap: on
line diff
--- a/Applications/Samples/SingleVolumeApplication.h Wed Nov 08 16:13:55 2017 +0100 +++ b/Applications/Samples/SingleVolumeApplication.h Wed Nov 08 17:27:38 2017 +0100 @@ -203,7 +203,7 @@ //ct->ScheduleLoadSeries("dd069910-4f090474-7d2bba07-e5c10783-f9e4fb1d"); ct->ScheduleLoadSeries("a04ecf01-79b2fc33-58239f7e-ad9db983-28e81afa"); // IBA //ct->ScheduleLoadSeries("03677739-1d8bca40-db1daf59-d74ff548-7f6fc9c0"); // 0522c0001 TCIA - + std::auto_ptr<OrthancVolumeImage> pet(new OrthancVolumeImage(context.GetWebService(), true)); //pet->ScheduleLoadSeries("aabad2e7-80702b5d-e599d26c-4f13398e-38d58a9e"); pet->ScheduleLoadInstance("830a69ff-8e4b5ee3-b7f966c8-bccc20fb-d322dceb"); // IBA 1 @@ -220,6 +220,7 @@ context.AddInteractor(new Interactor(*pet, *widget, projection, 1)); //context.AddInteractor(new VolumeImageInteractor(*ct, *widget, projection)); + context.AddVolume(ct.release()); context.AddVolume(pet.release());
--- a/Framework/Layers/DicomStructureSetRendererFactory.cpp Wed Nov 08 16:13:55 2017 +0100 +++ b/Framework/Layers/DicomStructureSetRendererFactory.cpp Wed Nov 08 17:27:38 2017 +0100 @@ -54,6 +54,11 @@ for (size_t k = 0; k < structureSet_.GetStructureCount(); k++) { + /*if (structureSet_.GetStructureName(k) != "CORD") + { + continue; + }*/ + std::vector< std::vector<DicomStructureSet::PolygonPoint> > polygons; if (structureSet_.ProjectStructure(polygons, k, slice_))
--- a/Framework/Toolbox/DicomStructureSet.cpp Wed Nov 08 16:13:55 2017 +0100 +++ b/Framework/Toolbox/DicomStructureSet.cpp Wed Nov 08 17:27:38 2017 +0100 @@ -28,8 +28,61 @@ #include <Plugins/Samples/Common/FullOrthancDataset.h> #include <Plugins/Samples/Common/DicomDatasetReader.h> +#include <limits> #include <stdio.h> #include <boost/lexical_cast.hpp> +#include <boost/geometry.hpp> +#include <boost/geometry/geometries/point_xy.hpp> +#include <boost/geometry/geometries/polygon.hpp> + +typedef boost::geometry::model::d2::point_xy<double> BoostPoint; +typedef boost::geometry::model::polygon<BoostPoint> BoostPolygon; +typedef boost::geometry::model::multi_polygon<BoostPolygon> BoostMultiPolygon; + + +static void Union(BoostMultiPolygon& output, + std::vector<BoostPolygon>& input) +{ + for (size_t i = 0; i < input.size(); i++) + { + boost::geometry::correct(input[i]); + } + + if (input.size() == 0) + { + output.clear(); + } + else if (input.size() == 1) + { + output.resize(1); + output[0] = input[0]; + } + else + { + boost::geometry::union_(input[0], input[1], output); + + for (size_t i = 0; i < input.size(); i++) + { + BoostMultiPolygon tmp; + boost::geometry::union_(output, input[i], tmp); + output = tmp; + } + } +} + + +static BoostPolygon CreateRectangle(float x1, float y1, + float x2, float y2) +{ + BoostPolygon r; + boost::geometry::append(r, BoostPoint(x1, y1)); + boost::geometry::append(r, BoostPoint(x1, y2)); + boost::geometry::append(r, BoostPoint(x2, y2)); + boost::geometry::append(r, BoostPoint(x2, y1)); + return r; +} + + namespace OrthancStone { @@ -449,7 +502,7 @@ Vector DicomStructureSet::GetNormal() const { - if (!referencedSlices_.empty()) + if (referencedSlices_.empty()) { Vector v; GeometryToolbox::AssignVector(v, 0, 0, 1); @@ -540,9 +593,81 @@ else if (GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisX()) || GeometryToolbox::IsParallelOrOpposite(isOpposite, normal, slice.GetAxisY())) { - // Sagittal or coronal projections + // Sagittal or coronal projection + + std::vector<BoostPolygon> projected; + + for (Polygons::iterator polygon = structure.polygons_.begin(); + polygon != structure.polygons_.end(); ++polygon) + { + polygon->UpdateReferencedSlice(referencedSlices_); + + // First, check that the polygon intersects the input slice + double zmin = std::numeric_limits<double>::infinity(); + double zmax = -std::numeric_limits<double>::infinity(); + + double zref = slice.ProjectAlongNormal(slice.GetOrigin()); + + for (Points::const_iterator p = polygon->GetPoints().begin(); + p != polygon->GetPoints().end(); ++p) + { + double z = slice.ProjectAlongNormal(*p); + zmin = std::min(zmin, z); + zmax = std::max(zmax, z); + } + + if (zmin < zref && zref < zmax) + { + // The polygon intersect the input slice, project the polygon + double xmin = std::numeric_limits<double>::infinity(); + double xmax = -std::numeric_limits<double>::infinity(); + double ymin = std::numeric_limits<double>::infinity(); + double ymax = -std::numeric_limits<double>::infinity(); + + for (Points::const_iterator p = polygon->GetPoints().begin(); + p != polygon->GetPoints().end(); ++p) + { + double x, y; + slice.ProjectPoint(x, y, *p); + xmin = std::min(xmin, x); + xmax = std::max(xmax, x); + ymin = std::min(ymin, y); + ymax = std::max(ymax, y); + } + + if (GeometryToolbox::IsNear(xmin, xmax)) + { + projected.push_back(CreateRectangle(xmin - polygon->GetSliceThickness() / 2.0, ymin, + xmax + polygon->GetSliceThickness() / 2.0, ymax)); + } + else if (GeometryToolbox::IsNear(ymin, ymax)) + { + projected.push_back(CreateRectangle(xmin, ymin - polygon->GetSliceThickness() / 2.0, + xmax, ymax + polygon->GetSliceThickness() / 2.0)); + } + else + { + // Should not happen + } + } + } + + BoostMultiPolygon merged; + Union(merged, projected); + + polygons.resize(merged.size()); + for (size_t i = 0; i < merged.size(); i++) + { + const std::vector<BoostPoint>& outer = merged[i].outer(); + + polygons[i].resize(outer.size()); + for (size_t j = 0; j < outer.size(); j++) + { + polygons[i][j] = std::make_pair(outer[j].x(), outer[j].y()); + } + } - return false; + return true; } else {