Line data    Source code

       1             : #include "Phantoms.h"
       2             : #include "EllipseGenerator.h"
       3             : #include "Logger.h"
       4             : #include "VolumeDescriptor.h"
       5             : #include "CartesianIndices.h"
       6             : #include "ForbildPhantom.h"
       7             : #include "ForbildData.h"
       8             : #include "Ellipsoid.h"
       9             : #include "Blobs.h"
      10             : 
      11             : #include <cmath>
      12             : #include <stdexcept>
      13             : 
      14             : namespace elsa::phantoms
      15             : {
      16             : 
      17             :     // scale sizes from [0,1] to the (square) phantom size, producing indices (integers)
      18             :     template <typename data_t>
      19             :     index_t scale(const DataDescriptor& dd, data_t value)
      20         550 :     {
      21         550 :         return std::lround(
      22         550 :             value * static_cast<data_t>(dd.getNumberOfCoefficientsPerDimension()[0] - 1) / 2.0f);
      23         550 :     }
      24             : 
      25             :     // scale and shift center coordinates to the (square) phantom size, producing indices
      26             :     // (integers)
      27             :     template <typename data_t>
      28             :     index_t scaleShift(const DataDescriptor& dd, data_t value)
      29         442 :     {
      30         442 :         return std::lround(static_cast<double>(value)
      31         442 :                            * static_cast<double>(dd.getNumberOfCoefficientsPerDimension()[0] - 1)
      32         442 :                            / 2.0)
      33         442 :                + (dd.getNumberOfCoefficientsPerDimension()[0] / 2);
      34         442 :     }
      35             : 
      36             :     template <typename data_t>
      37             :     DataContainer<data_t> circular(IndexVector_t volumesize, data_t radius)
      38          32 :     {
      39          32 :         VolumeDescriptor dd(volumesize);
      40          32 :         DataContainer<data_t> dc(dd);
      41          32 :         dc = 0;
      42             : 
      43          32 :         const Vector_t<data_t> sizef = volumesize.template cast<data_t>();
      44          32 :         const auto center = (sizef.array() / 2).matrix();
      45             : 
      46       69632 :         for (auto pos : CartesianIndices(volumesize)) {
      47       69632 :             const Vector_t<data_t> p = pos.template cast<data_t>();
      48       69632 :             if ((p - center).norm() <= radius) {
      49       12038 :                 dc(pos) = 1;
      50       12038 :             }
      51       69632 :         }
      52             : 
      53          32 :         return dc;
      54          32 :     }
      55             : 
      56             :     template <typename data_t>
      57             :     DataContainer<data_t> rectangle(IndexVector_t volumesize, IndexVector_t lower,
      58             :                                     IndexVector_t upper)
      59           8 :     {
      60           8 :         VolumeDescriptor dd(volumesize);
      61           8 :         DataContainer<data_t> dc(dd);
      62           8 :         dc = 0;
      63             : 
      64        9856 :         for (auto pos : CartesianIndices(lower, upper)) {
      65        9856 :             dc(pos) = 1;
      66        9856 :         }
      67             : 
      68           8 :         return dc;
      69           8 :     }
      70             : 
      71             :     template <typename data_t>
      72             :     DataContainer<data_t> modifiedSheppLogan(IndexVector_t sizes)
      73          24 :     {
      74             :         // sanity check
      75          24 :         if (sizes.size() < 2 || sizes.size() > 3)
      76           0 :             throw InvalidArgumentError("phantom::modifiedSheppLogan: only 2d or 3d supported");
      77          24 :         if (sizes.size() == 2 && sizes[0] != sizes[1])
      78           0 :             throw InvalidArgumentError("phantom::modifiedSheppLogan: 2d size has to be square");
      79          24 :         if (sizes.size() == 3 && (sizes[0] != sizes[1] || sizes[0] != sizes[2]))
      80           0 :             throw InvalidArgumentError("phantom::modifiedSheppLogan: 3d size has to be cubed");
      81             : 
      82          24 :         Logger::get("phantom::modifiedSheppLogan")
      83          24 :             ->info("creating modified Shepp Logan phantom of size {}^{}", sizes[0], sizes.size());
      84             : 
      85          24 :         VolumeDescriptor dd(sizes);
      86          24 :         DataContainer<data_t> dc(dd);
      87          24 :         dc = 0;
      88             : 
      89          24 :         if (sizes.size() == 2) {
      90          17 :             EllipseGenerator<data_t>::drawFilledEllipse2d(dc, 1.0,
      91          17 :                                                           {scaleShift(dd, 0), scaleShift(dd, 0)},
      92          17 :                                                           {scale(dd, 0.69f), scale(dd, 0.92f)}, 0);
      93          17 :             EllipseGenerator<data_t>::drawFilledEllipse2d(
      94          17 :                 dc, -0.8f, {scaleShift(dd, 0), scaleShift(dd, -0.0184f)},
      95          17 :                 {scale(dd, 0.6624f), scale(dd, 0.8740f)}, 0);
      96          17 :             EllipseGenerator<data_t>::drawFilledEllipse2d(
      97          17 :                 dc, -0.2f, {scaleShift(dd, 0.22f), scaleShift(dd, 0)},
      98          17 :                 {scale(dd, 0.11f), scale(dd, 0.31f)}, -18);
      99          17 :             EllipseGenerator<data_t>::drawFilledEllipse2d(
     100          17 :                 dc, -0.2f, {scaleShift(dd, -0.22f), scaleShift(dd, 0)},
     101          17 :                 {scale(dd, 0.16f), scale(dd, 0.41f)}, 18);
     102          17 :             EllipseGenerator<data_t>::drawFilledEllipse2d(
     103          17 :                 dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, 0.35f)},
     104          17 :                 {scale(dd, 0.21f), scale(dd, 0.25)}, 0);
     105          17 :             EllipseGenerator<data_t>::drawFilledEllipse2d(
     106          17 :                 dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, 0.1f)},
     107          17 :                 {scale(dd, 0.046f), scale(dd, 0.046f)}, 0);
     108          17 :             EllipseGenerator<data_t>::drawFilledEllipse2d(
     109          17 :                 dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, -0.1f)},
     110          17 :                 {scale(dd, 0.046f), scale(dd, 0.046f)}, 0);
     111          17 :             EllipseGenerator<data_t>::drawFilledEllipse2d(
     112          17 :                 dc, 0.1f, {scaleShift(dd, -0.08f), scaleShift(dd, -0.605f)},
     113          17 :                 {scale(dd, 0.046f), scale(dd, 0.023f)}, 0);
     114          17 :             EllipseGenerator<data_t>::drawFilledEllipse2d(
     115          17 :                 dc, 0.1f, {scaleShift(dd, 0), scaleShift(dd, -0.606f)},
     116          17 :                 {scale(dd, 0.023f), scale(dd, 0.023f)}, 0);
     117          17 :             EllipseGenerator<data_t>::drawFilledEllipse2d(
     118          17 :                 dc, 0.1f, {scaleShift(dd, 0.06f), scaleShift(dd, -0.605f)},
     119          17 :                 {scale(dd, 0.023f), scale(dd, 0.046f)}, 0);
     120          17 :         }
     121             : 
     122          24 :         if (sizes.size() == 3) {
     123          70 :             for (std::array<data_t, 10> e : modifiedSheppLoganParameters<data_t>) {
     124             : 
     125          70 :                 Vec3X<data_t> halfAxis{data_t(scale(dd, e[1])), data_t(scale(dd, e[2])),
     126          70 :                                        data_t(scale(dd, e[3]))};
     127             : 
     128          70 :                 if (halfAxis[0] < 1 || halfAxis[1] < 1 || halfAxis[2] < 1 || e[0] == data_t(0)) {
     129          36 :                     Logger::get("phantom::modifiedSheppLogan")
     130          36 :                         ->warn(
     131          36 :                             "Ellipsoid will not be rendered, because of amplitude=0 or an invalid "
     132          36 :                             "half axis! amplitude {}, half axis ({},{},{}) ",
     133          36 :                             e[0], halfAxis[0], halfAxis[1], halfAxis[2]);
     134          36 :                     continue;
     135          36 :                 }
     136             : 
     137          34 :                 Ellipsoid<data_t> ellipsoid{
     138          34 :                     e[0],
     139          34 :                     {scaleShift(dd, e[4]), scaleShift(dd, e[5]), scaleShift(dd, e[6])},
     140          34 :                     halfAxis,
     141          34 :                     {e[7], e[8], e[9]}};
     142          34 :                 Logger::get("phantom::modifiedSheppLogan")->info("rasterize {}", ellipsoid);
     143          34 :                 rasterize<Blending::ADDITION>(ellipsoid, dd, dc);
     144          34 :             }
     145           7 :         }
     146             : 
     147          24 :         return dc;
     148          24 :     }
     149             : 
     150             :     template <typename data_t>
     151             :     DataContainer<data_t> smoothBlob(IndexVector_t sizes, double radius_manipulation)
     152           0 :     {
     153             :         // sanity check
     154           0 :         if (sizes.size() < 2 || sizes.size() > 3)
     155           0 :             throw InvalidArgumentError("phantom::smoothBlob: only 2d or 3d supported");
     156           0 :         if (sizes.size() == 2 && sizes[0] != sizes[1])
     157           0 :             throw InvalidArgumentError("phantom::smoothBlob: 2d size has to be square");
     158           0 :         if (sizes.size() == 3 && (sizes[0] != sizes[1] || sizes[0] != sizes[2]))
     159           0 :             throw InvalidArgumentError("phantom::smoothBlob: 3d size has to be cubed");
     160             : 
     161           0 :         Logger::get("phantom::smoothBlob")
     162           0 :             ->info("creating smooth Blob phantom of size {}^{}", sizes[0], sizes.size());
     163             : 
     164           0 :         VolumeDescriptor dd(sizes);
     165           0 :         DataContainer<data_t> dc(dd);
     166           0 :         dc = 0;
     167             : 
     168           0 :         auto radius = sizes[0] / 2;
     169             : 
     170           0 : #pragma omp parallel for
     171           0 :         for (index_t i = 0; i < dc.getSize(); i++) {
     172           0 :             const RealVector_t coord =
     173           0 :                 dd.getCoordinateFromIndex(i).template cast<real_t>().array() + 0.5;
     174           0 :             data_t distance_from_center = (coord - dd.getLocationOfOrigin()).norm();
     175           0 :             dc[i] =
     176           0 :                 blobs::blob_evaluate(distance_from_center, radius * radius_manipulation, 10.83f, 2);
     177           0 :         }
     178             : 
     179           0 :         return dc;
     180           0 :     }
     181             : 
     182             :     template <typename data_t>
     183             :     DataContainer<data_t> forbild(IndexVector_t sizes, ForbildPhantom<data_t>& fp)
     184           3 :     {
     185             :         // sanity checks
     186           3 :         if (sizes.size() != 3 || sizes[0] != sizes[1] || sizes[0] != sizes[2]) {
     187           0 :             throw InvalidArgumentError("phantom::forbild: 3d size has to be cubed");
     188           0 :         }
     189             : 
     190           3 :         VolumeDescriptor dd(sizes);
     191           3 :         DataContainer<data_t> dc(dd);
     192           3 :         dc = 0;
     193             : 
     194         905 :         for (int order = 0; order < fp.getMaxOrderIndex(); order++) {
     195         902 :             if (fp.getEllipsoids().count(order)) {
     196             :                 // Ellipsoids
     197          15 :                 auto e = fp.getEllipsoids().at(order);
     198          15 :                 Logger::get("phantoms::forbild")->info("rasterize {}", e);
     199          15 :                 rasterize<Blending::OVERWRITE>(e, dd, dc);
     200         887 :             } else if (fp.getEllipsoidsClippedX().count(order)) {
     201             :                 // Ellipsoids clipped at x axis
     202           1 :                 auto [el, clipping] = fp.getEllipsoidsClippedX().at(order);
     203           1 :                 Logger::get("phantoms::forbild")->info("rasterize {} with clipping", el);
     204           1 :                 rasterizeWithClipping<Blending::OVERWRITE>(el, dd, dc, clipping);
     205         886 :             } else if (fp.getSpheres().count(order)) {
     206             :                 // Spheres
     207          42 :                 auto s = fp.getSpheres().at(order);
     208          42 :                 Logger::get("phantoms::forbild")->info("rasterize {}", s);
     209          42 :                 rasterize<Blending::OVERWRITE>(s, dd, dc);
     210         844 :             } else if (fp.getEllipCylinders().count(order)) {
     211             :                 // EllipCylinders
     212           4 :                 auto e = fp.getEllipCylinders().at(order);
     213           4 :                 Logger::get("phantoms::forbild")->info("rasterize {}", e);
     214           4 :                 rasterize<Blending::OVERWRITE>(e, dd, dc);
     215         840 :             } else if (fp.getEllipCylindersFree().count(order)) {
     216             :                 // EllipCylinders free
     217           5 :                 auto e = fp.getEllipCylindersFree().at(order);
     218           5 :                 Logger::get("phantoms::forbild")->info("rasterize {}", e);
     219           5 :                 rasterize<Blending::OVERWRITE>(e, dd, dc);
     220         835 :             } else if (fp.getCylinders().count(order)) {
     221             :                 // Cylinders
     222          43 :                 auto c = fp.getCylinders().at(order);
     223          43 :                 Logger::get("phantoms::forbild")->info("rasterize {}", c);
     224          43 :                 rasterize<Blending::OVERWRITE>(c, dd, dc);
     225         792 :             } else if (fp.getCylindersFree().count(order)) {
     226             :                 // Cylinders free
     227         185 :                 auto cf = fp.getCylindersFree().at(order);
     228         185 :                 Logger::get("phantoms::forbild")->info("rasterize {}", cf);
     229         185 :                 rasterize<Blending::OVERWRITE>(cf, dd, dc);
     230         607 :             } else if (fp.getBoxes().count(order)) {
     231             :                 // Boxes
     232          40 :                 auto b = fp.getBoxes().at(order);
     233          40 :                 Logger::get("phantoms::forbild")->info("rasterize {}", b);
     234          40 :                 rasterize<Blending::OVERWRITE>(b, dd, dc);
     235         567 :             } else {
     236         567 :                 Logger::get("phantoms::forbild")->warn("No object found for order index {}", order);
     237         567 :             }
     238         902 :         }
     239             : 
     240           3 :         return dc;
     241           3 :     }
     242             : 
     243             :     template <typename data_t>
     244             :     DataContainer<data_t> forbildHead(IndexVector_t sizes)
     245           1 :     {
     246             : 
     247           1 :         Logger::get("phantom::forbildHead")
     248           1 :             ->info("creating forbild of size {}^{}", sizes[0], sizes.size());
     249             :         // field view from http://ftp.imp.uni-erlangen.de/phantoms/head/head.html
     250           1 :         ForbildPhantom<data_t> fp{sizes[0] - 1, data_t(25.6f), forbilddata::maxOrderIndex_head};
     251           1 :         fp.addEllipsoids(forbilddata::ellipsoids_head<data_t>);
     252           1 :         fp.addEllipsoidsClippedX(forbilddata::ellipsoidsClippingX_head<data_t>);
     253           1 :         fp.addSpheres(forbilddata::spheres_head<data_t>);
     254           1 :         fp.addEllipCylinders(forbilddata::ellipCyls_head<data_t>);
     255           1 :         fp.addEllipCylindersFree(forbilddata::ellipCylsFree_head<data_t>);
     256           1 :         fp.addCylinders(forbilddata::cyls_head<data_t>);
     257           1 :         fp.addCylindersFree(forbilddata::cylsFree_head<data_t>);
     258           1 :         fp.addBoxes(forbilddata::boxes_head<data_t>);
     259           1 :         return std::move(forbild<data_t>(sizes, fp));
     260           1 :     }
     261             : 
     262             :     template <typename data_t>
     263             :     DataContainer<data_t> forbildThorax(IndexVector_t sizes)
     264           1 :     {
     265             : 
     266           1 :         Logger::get("phantom::forbildThorax")
     267           1 :             ->info("creating forbild of size {}^{}", sizes[0], sizes.size());
     268             :         // field view from http://ftp.imp.uni-erlangen.de/phantoms/thorax/thorax.htm
     269           1 :         ForbildPhantom<data_t> fp{sizes[0] - 1, data_t(50.f), forbilddata::maxOrderIndex_thorax};
     270           1 :         fp.addEllipsoids(forbilddata::ellipsoids_thorax<data_t>);
     271           1 :         fp.addSpheres(forbilddata::spheres_thorax<data_t>);
     272           1 :         fp.addEllipCylinders(forbilddata::ellipCyls_thorax<data_t>);
     273           1 :         fp.addEllipCylindersFree(forbilddata::ellipCylsFree_thorax<data_t>);
     274           1 :         fp.addCylinders(forbilddata::cyls_thorax<data_t>);
     275           1 :         fp.addCylindersFree(forbilddata::cylsFree_thorax<data_t>);
     276           1 :         fp.addBoxes(forbilddata::boxes_thorax<data_t>);
     277             : 
     278           1 :         return std::move(forbild<data_t>(sizes, fp));
     279           1 :     }
     280             : 
     281             :     template <typename data_t>
     282             :     DataContainer<data_t> forbildAbdomen(IndexVector_t sizes)
     283           1 :     {
     284             : 
     285           1 :         Logger::get("phantom::forbildAbdomen")
     286           1 :             ->info("creating forbild of size {}^{}", sizes[0], sizes.size());
     287             :         // field view from http://ftp.imp.uni-erlangen.de/phantoms/abdomen/abdomen.html
     288           1 :         ForbildPhantom<data_t> fp{sizes[0] - 1, data_t(40.f), forbilddata::maxOrderIndex_abdomen};
     289           1 :         fp.addEllipsoids(forbilddata::ellipsoids_abdomen<data_t>);
     290           1 :         fp.addSpheres(forbilddata::spheres_abdomen<data_t>);
     291           1 :         fp.addEllipCylinders(forbilddata::ellipCyls_abdomen<data_t>);
     292           1 :         fp.addEllipCylindersFree(forbilddata::ellipCylsFree_abdomen<data_t>);
     293           1 :         fp.addCylinders(forbilddata::cyls_abdomen<data_t>);
     294           1 :         fp.addCylindersFree(forbilddata::cylsFree_abdomen<data_t>);
     295           1 :         fp.addBoxes(forbilddata::boxes_abdomen<data_t>);
     296           1 :         return std::move(forbild<data_t>(sizes, fp));
     297           1 :     }
     298             : 
     299             :     // ------------------------------------------
     300             :     // explicit template instantiation
     301             :     template DataContainer<float> circular<float>(IndexVector_t volumesize, float radius);
     302             :     template DataContainer<double> circular<double>(IndexVector_t volumesize, double radius);
     303             : 
     304             :     template DataContainer<float> modifiedSheppLogan<float>(IndexVector_t sizes);
     305             :     template DataContainer<double> modifiedSheppLogan<double>(IndexVector_t sizes);
     306             : 
     307             :     template DataContainer<float> forbildHead<float>(IndexVector_t sizes);
     308             :     template DataContainer<double> forbildHead<double>(IndexVector_t sizes);
     309             : 
     310             :     template DataContainer<float> forbildAbdomen<float>(IndexVector_t sizes);
     311             :     template DataContainer<double> forbildAbdomen<double>(IndexVector_t sizes);
     312             : 
     313             :     template DataContainer<float> forbildThorax<float>(IndexVector_t sizes);
     314             :     template DataContainer<double> forbildThorax<double>(IndexVector_t sizes);
     315             : 
     316             :     template DataContainer<float> smoothBlob<float>(IndexVector_t sizes,
     317             :                                                     double radius_manipulation);
     318             :     template DataContainer<double> smoothBlob<double>(IndexVector_t sizes,
     319             :                                                       double radius_manipulation);
     320             : 
     321             :     template DataContainer<float> rectangle<float>(IndexVector_t volumesize, IndexVector_t lower,
     322             :                                                    IndexVector_t upper);
     323             :     template DataContainer<double> rectangle<double>(IndexVector_t volumesize, IndexVector_t lower,
     324             :                                                      IndexVector_t upper);
     325             : } // namespace elsa::phantoms