Mercurial > hg > orthanc-stone
view OrthancStone/Sources/Toolbox/DicomStructureSetUtils.cpp @ 1757:28979b77ce90
Suppressed one pass on the dose distribution computation
author | bgo@SHARKNADO.localdomain |
---|---|
date | Mon, 26 Apr 2021 17:37:54 +0200 |
parents | 9ac2a65d4172 |
children | 3889ae96d2e9 |
line wrap: on
line source
/** * Stone of Orthanc * Copyright (C) 2012-2016 Sebastien Jodogne, Medical Physics * Department, University Hospital of Liege, Belgium * Copyright (C) 2017-2021 Osimis S.A., Belgium * * This program is free software: you can redistribute it and/or * modify it under the terms of the GNU Lesser 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 * Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with this program. If not, see * <http://www.gnu.org/licenses/>. **/ #include "DicomStructureSetUtils.h" namespace OrthancStone { #if 0 void DicomStructure2::PartitionRectangleList(std::vector< std::vector<size_t> > & sets, const std::vector<RtStructRectanglesInSlab> slabCuts) { // map position ( <slabIndex,rectIndex> )--> disjoint set index std::map<std::pair<size_t, size_t>, size_t> posToIndex; // disjoint set index --> position std::map<size_t, std::pair<size_t, size_t> > indexToPos; size_t nextIndex = 0; for (size_t i = 0; i < slabCuts.size(); ++i) { for (size_t j = 0; j < slabCuts[i].size(); ++j) { std::pair<size_t, size_t> pos(i, j); posToIndex<pos> = nextIndex; indexToPos<nextIndex> = pos; } } // nextIndex is now the total rectangle count DisjointDataSet ds(nextIndex); // we loop on all slabs (except the last one) and we connect all rectangles if (slabCuts.size() < 2) { #error write special case } else { for (size_t i = 0; i < slabCuts.size() - 1; ++i) { for (size_t j = 0; j < slabCuts[i].size(); ++j) { const RtStructRectangleInSlab& r1 = slabCuts[i][j]; const size_t r1i = posToIndex(std::pair<size_t, size_t>(i, j)); for (size_t k = 0; k < slabCuts[i + 1].size(); ++k) { const RtStructRectangleInSlab& r2 = slabCuts[i + 1][k]; const size_t r2i = posToIndex(std::pair<size_t, size_t>(i, j)); // rect.xmin <= rectBottom.xmax && rectBottom.xmin <= rect.xmax if ((r1.xmin <= r2.xmax) && (r2.xmin <= r1.xmax)) { #error now go! } } } } } #endif /* compute list of segments : numberOfRectsFromHereOn = 0 possibleNext = {in_k,in_kplus1} for all boundaries: - we create a vertical segment and we push it - if boundary is a start, numberOfRectsFromHereOn += 1. - if we switch from 0 to 1, we start a segment - if we switch from 1 to 2, we end the current segment and we record it - if boundary is an end, numberOfRectsFromHereOn -= 1. - if we switch from 1 to 0, we end the current segment and we record it - if we switch from 2 to 1, we start a segment */ // static void AddSlabBoundaries( std::vector<std::pair<double, RectangleBoundaryKind> > & boundaries, const std::vector<RtStructRectanglesInSlab> & slabCuts, size_t iSlab) { if (iSlab < slabCuts.size()) { const RtStructRectanglesInSlab& slab = slabCuts[iSlab]; for (size_t iRect = 0; iRect < slab.size(); ++iRect) { const RtStructRectangleInSlab& rect = slab[iRect]; { std::pair<double, RectangleBoundaryKind> boundary(rect.xmin, RectangleBoundaryKind_Start); boundaries.insert(std::lower_bound(boundaries.begin(), boundaries.end(), boundary), boundary); } { std::pair<double, RectangleBoundaryKind> boundary(rect.xmax, RectangleBoundaryKind_End); boundaries.insert(std::lower_bound(boundaries.begin(), boundaries.end(), boundary), boundary); } } } } // static void ProcessBoundaryList( std::vector< std::pair<Point2D, Point2D> > & segments, const std::vector<std::pair<double, RectangleBoundaryKind> > & boundaries, double y) { Point2D start; Point2D end; int curNumberOfSegments = 0; // we count the number of segments. we only draw if it is 1 (not 0 or 2) for (size_t i = 0; i < boundaries.size(); ++i) { switch (boundaries[i].second) { case RectangleBoundaryKind_Start: curNumberOfSegments += 1; switch (curNumberOfSegments) { case 0: assert(false); break; case 1: // a new segment has begun! start.x = boundaries[i].first; start.y = y; break; case 2: // an extra segment has begun : stop the current one (we don't draw overlaps) end.x = boundaries[i].first; end.y = y; segments.push_back(std::pair<Point2D, Point2D>(start, end)); break; default: //assert(false); // seen IRL ! break; } break; case RectangleBoundaryKind_End: curNumberOfSegments -= 1; switch (curNumberOfSegments) { case 0: // a lone (thus active) segment has ended. end.x = boundaries[i].first; end.y = y; segments.push_back(std::pair<Point2D, Point2D>(start, end)); break; case 1: // an extra segment has ended : start a new one one start.x = boundaries[i].first; start.y = y; break; default: // this should not happen! //assert(false); break; } break; default: assert(false); break; } } } #if 0 void ConvertListOfSlabsToSegments( std::vector< std::pair<Point2D, Point2D> >& segments, const std::vector<RtStructRectanglesInSlab>& slabCuts, const size_t totalRectCount) { #error to delete } #else // See https://www.dropbox.com/s/bllco6q8aazxk44/2019-09-18-rtstruct-cut-algorithm-rect-merge.png void ConvertListOfSlabsToSegments( std::vector< std::pair<Point2D, Point2D> > & segments, const std::vector<RtStructRectanglesInSlab> & slabCuts, const size_t totalRectCount) { if (slabCuts.size() == 0) return; if (totalRectCount > 0) segments.reserve(4 * totalRectCount); // worst case, but common. /* VERTICAL */ for (size_t iSlab = 0; iSlab < slabCuts.size(); ++iSlab) { for (size_t iRect = 0; iRect < slabCuts[iSlab].size(); ++iRect) { const RtStructRectangleInSlab& rect = slabCuts[iSlab][iRect]; { Point2D p1(rect.xmin, rect.ymin); Point2D p2(rect.xmin, rect.ymax); segments.push_back(std::pair<Point2D, Point2D>(p1, p2)); } { Point2D p1(rect.xmax, rect.ymin); Point2D p2(rect.xmax, rect.ymax); segments.push_back(std::pair<Point2D, Point2D>(p1, p2)); } } } /* HORIZONTAL */ // if we have N slabs, we have N+1 potential vertical positions for horizontal segments // - one for top of slab 0 // - N-1 for all positions between two slabs // - one for bottom of slab N-1 // this adds all the horizontal segments for the tops of 3the rectangles // in row 0 if (slabCuts[0].size() > 0) { std::vector<std::pair<double, RectangleBoundaryKind> > boundaries; AddSlabBoundaries(boundaries, slabCuts, 0); ProcessBoundaryList(segments, boundaries, slabCuts[0][0].ymin); } // this adds all the horizontal segments belonging to two slabs for (size_t iSlab = 0; iSlab < slabCuts.size() - 1; ++iSlab) { std::vector<std::pair<double, RectangleBoundaryKind> > boundaries; AddSlabBoundaries(boundaries, slabCuts, iSlab); AddSlabBoundaries(boundaries, slabCuts, iSlab + 1); double curY = 0; if (slabCuts[iSlab].size() > 0) { curY = slabCuts[iSlab][0].ymax; ProcessBoundaryList(segments, boundaries, curY); } else if (slabCuts[iSlab + 1].size() > 0) { curY = slabCuts[iSlab + 1][0].ymin; ProcessBoundaryList(segments, boundaries, curY); } else { // nothing to do!! : both slab lists are empty! } } // this adds all the horizontal segments for the BOTTOM of the rectangles // on last row if (slabCuts[slabCuts.size() - 1].size() > 0) { std::vector<std::pair<double, RectangleBoundaryKind> > boundaries; AddSlabBoundaries(boundaries, slabCuts, slabCuts.size() - 1); ProcessBoundaryList(segments, boundaries, slabCuts[slabCuts.size() - 1][0].ymax); } } #endif }