24.09., 9:00 - 11:00: Due to updates GitLab will be unavailable for some minutes between 09:00 and 11:00.

Commit 31fc08ee authored by Bianca-Maria Borza's avatar Bianca-Maria Borza

my modifications added

parent 19d733a0
Pipeline #289821 failed with stages
in 47 seconds
......@@ -59,8 +59,9 @@ namespace elsa
// enables small optimizations
bool hasRotation = (std::abs(phi) + std::abs(theta) + std::abs(psi)) > 0;
// convert to radians
if(hasRotation || !hasRotation)
EllipsoidGenerator<data_t>::drawFilledEllipsoid(dc,amplitude,center,sizes);
else{// convert to radians
auto phiRad = phi * pi<double> / 180.0;
auto thetaRad = theta * pi<double> / 180.0;
auto psiRad = psi * pi<double> / 180.0;
......@@ -145,6 +146,7 @@ namespace elsa
}
}
}
}
}
template <typename data_t>
......
......@@ -2,6 +2,7 @@
#include "elsaDefines.h"
#include "DataContainer.h"
#include "EllipsoidGenerator.h"
namespace elsa
{
......
#include "EllipsoidGenerator.h"
#include "Timer.h"
#include "Logger.h"
#include <stdexcept>
namespace elsa
{
template <typename data_t>
void EllipsoidGenerator<data_t>::drawFilledEllipsoid(DataContainer<data_t>& dc,
data_t amplitude, Vec3 center, Vec3 sizes)
{
auto twoSizeXSquared = 2*sizes[0]*sizes[0];
auto twoSizeYSquared = 2*sizes[1]*sizes[1];
auto twoSizeZSquared = 2*sizes[2]*sizes[2];
auto x = sizes[0];
auto y = sizes[1];
index_t z = 0;
auto xChange = sizes[2]*sizes[2]*(1-2*sizes[0]);
auto xzChange = sizes[0]*sizes[0];
auto yChange = sizes[2]*sizes[2]*(1-2*sizes[1]);
auto yzChange = sizes[1]*sizes[1];
index_t ellipseErrorX = 0;
index_t ellipseErrorY = 0;
auto xStop = twoSizeZSquared * sizes[0];
index_t xzStop = 0;
auto yStop = twoSizeZSquared * sizes[1];
index_t yzStop = 0;
while ( xStop >= xzStop && yStop >= yzStop )
{
drawFilledEllipse(dc,amplitude,center,x,y,center[2]+z);
if ( z != 0 ){
drawFilledEllipse(dc,amplitude,center,x,y,center[2]-z);
}
z += 1;
xzStop += twoSizeXSquared;
yzStop += twoSizeYSquared;
ellipseErrorX += xzChange;
ellipseErrorY += yzChange;
xzChange += twoSizeXSquared;
yzChange += twoSizeYSquared;
if ((2*ellipseErrorX + xChange) > 0){
x -= 1;
xStop -= twoSizeZSquared;
ellipseErrorX += xChange;
xChange += twoSizeZSquared;
}
if ((2*ellipseErrorY + yChange) > 0){
y -= 1;
yStop -= twoSizeZSquared;
ellipseErrorY += yChange;
yChange += twoSizeZSquared;
}
}
if ( xStop >= xzStop){
index_t listY[sizes[2]];
int i = 0;
y = -1;
yChange = sizes[2]*sizes[2];
yzChange = sizes[1]*sizes[1]*(1-2*sizes[2]);
ellipseErrorY = 0;
yStop = 0;
yzStop = twoSizeYSquared * sizes[2];
while (yStop <yzStop){
y += 1;
yStop += twoSizeZSquared;
ellipseErrorY += yChange;
yChange += twoSizeZSquared;
if ((2*ellipseErrorY + yzChange) > 0){
listY[i] = y;
i += 1;
yzStop -= twoSizeYSquared;
ellipseErrorY += yzChange;
yzChange += twoSizeYSquared;
}
}
i -= 1;
while (xStop >= xzStop){
drawfilledEllipse(dc,amplitude,center,x,listY[i],center[2]+z);
if(z != 0){
drawfilledEllipse(dc,amplitude,center,x,listY[i],center[2]-z);
}
i -= 1;
z += 1;
xzStop += twoSizeXSquared;
ellipseErrorX += xzChange;
xzChange += twoSizeXSquared;
if((2*ellipseErrorX + xChange) > 0 ){
x -= 1;
xStop -= twoSizeZSquared;
ellipseErrorX += xChange;
xChange += twoSizeZSquared;
}
}
i = 0;
x = -1;
z = sizes[2];
xChange = sizes[2]*sizes[2];
xzChange = sizes[0]*sizes[0]*(1-2*sizes[2]);
ellipseErrorX = 0;
xStop = 0;
xzStop = twoSizeXSquared*sizes[2];
while(xStop < xzStop){
while((2*ellipseErrorX+xzChange)<=0){
x += 1;
xStop += twoSizeZSquared;
ellipseErrorX += xChange;
xChange += twoSizeZSquared;
}
drawfilledEllipse(dc,amplitude,center,x,listY[i],center[2]+z);
if(z!=0){
drawfilledEllipse(dc,amplitude,center,x,listY[i],center[2]-z);
}
i += 1;
z -= 1;
xzStop -= twoSizeXSquared;
ellipseErrorX += xzChange;
xzChange += twoSizeXSquared;
}
}
else{
index_t listX [sizes[2]];
int i = 0;
x = -1;
xChange = sizes[2]*sizes[2];
xzChange = sizes[0]*sizes[0]*(1-2*sizes[2]);
ellipseErrorX = 0;
xStop = 0;
xzStop = twoSizeXSquared*sizes[2];
while(xStop < xzStop){
x += 1;
xStop += twoSizeZSquared;
ellipseErrorX += xChange;
xChange += twoSizeZSquared;
if ((2*ellipseErrorX+xzChange)>0){
listX[i] = x;
i += 1;
xzStop -= twoSizeXSquared;
ellipseErrorX += xzChange;
xzChange += twoSizeXSquared;
}
}
i -= 1;
while(yStop >= yzStop){
drawfilledEllipse(dc,amplitude,center,listX[i],y,center[2]+z);
if(z!=0){
drawfilledEllipse(dc,amplitude,center,listX[i],y,center[2]-z);
}
i -= 1;
z += 1;
yzStop += twoSizeYSquared;
ellipseErrorY += yzChange;
yzChange += twoSizeYSquared;
if((2*ellipseErrorY + yChange) > 0){
y -= 1;
yStop -= twoSizeZSquared;
ellipseErrorY += yChange;
yChange += twoSizeZSquared;
}
}
i = 0;
y = -1;
z = sizes[2];
yChange = sizes[2]*sizes[2];
yzChange = sizes[1]*sizes[1]*(1-2*sizes[2]);
ellipseErrorY = 0;
yStop = 0;
yzStop = twoSizeYSquared*sizes[2];
while(yStop <yzStop){
while((2*ellipseErrorY+yzChange) <= 0){
y += 1;
yStop += twoSizeZSquared;
ellipseErrorY += yChange;
yChange += twoSizeZSquared;
}
drawfilledEllipse(dc, amplitude, center,listX[i],y,center[2]+z);
if(z!=0){
drawfilledEllipse(dc,amplitude,center,listX[i],y,center[2]-z);
}
i += 1;
z -= 1;
yzStop -= twoSizeYSquared;
ellipseErrorY += yzChange;
yzChange += twoSizeYSquared;
}
}
}
template <typename data_t>
void EllipsoidGenerator<data_t>::drawFilledEllipse(DataContainer<data_t>& dc,
data_t amplitude, Vec3 center, data_t radiusX, data_t radiusY, index_t z){
auto twoSizeXSquared = 2 * radiusX * radiusX;
auto twoSizeYSquared = 2 * radiusY * radiusY;
auto x = radiusX;
index_t y = 0;
auto xChange = radiusY * radiusY * (1 - 2 * radiusX);
auto yChange = radiusX * radiusX;
index_t ellipseError = 0;
auto xStop = twoSizeYSquared * radiusX;
index_t yStop = 0;
while (xStop >= yStop) {
drawLine(dc, amplitude, center, x, y, z);
y += 1;
yStop += twoSizeXSquared;
ellipseError += yChange;
yChange += twoSizeXSquared;
if ((2 * ellipseError + xChange) > 0) {
x -= 1;
xStop -= twoSizeYSquared;
ellipseError += xChange;
xChange += twoSizeYSquared;
}
}
x = -1;
y = radiusY;
xChange = radiusY * radiusY;
yChange = radiusX * radiusX * (1 - 2 * radiusY);
ellipseError = 0;
xStop = 0;
yStop = twoSizeXSquared * radiusY;
while (xStop < yStop) {
x += 1;
xStop += twoSizeYSquared;
ellipseError += xChange;
xChange += twoSizeYSquared;
if ((2 * ellipseError + yChange) > 0) {
drawLine(dc, amplitude, center, x, y, z);
y -= 1;
yStop -= twoSizeXSquared;
ellipseError += yChange;
yChange += twoSizeXSquared;
}
}
}
template <typename data_t>
void EllipsoidGenerator<data_t>::drawLine(DataContainer<data_t>&dc,
data_t amplitude, Vec3 center, index_t xOffset, index_t yOffset, index_t z){
IndexVector_t coord(3);
coord[2] = z;
for (index_t x = center[0] - xOffset; x <= center[0] + xOffset; ++x) {
coord[0] = x;
coord[1] = center[1] + yOffset;
if (coord[0] < 0
|| coord[0] >= dc.getDataDescriptor().getNumberOfCoefficientsPerDimension()[0])
throw std::invalid_argument("EllipseGenerator::drawShearedLinePairs2d: drawing "
"coordinate (x) out of bounds");
if (coord[1] < 0
|| coord[1] >= dc.getDataDescriptor().getNumberOfCoefficientsPerDimension()[1])
throw std::invalid_argument("EllipseGenerator::drawShearedLinePairs2d: drawing "
"coordinate (y) out of bounds");
dc(coord) += amplitude;
if (yOffset != 0) {
coord[1] = center[1] - yOffset;
if (coord[1] < 0
|| coord[1] >= dc.getDataDescriptor().getNumberOfCoefficientsPerDimension()[1])
throw std::invalid_argument("EllipseGenerator::drawShearedLinePairs2d: drawing "
"coordinate (y) out of bounds");
dc(coord) += amplitude;
}
}
}
// ------------------------------------------
// explicit template instantiation
template class EllipsoidGenerator<float>;
template class EllipsoidGenerator<double>;
}
\ No newline at end of file
#pragma once
#include "elsaDefines.h"
#include "DataContainer.h"
namespace elsa
{
/**
* \brief This class draws 3d spheres into DataContainers
*
* \tparam data_t data type for the DataContainer, defaulting to real_t
*
* Drawing is done by adding a given constant to the appropriate places in the 3d image.
*/
template <typename data_t = real_t>
class EllipsoidGenerator
{
public:
/// short-hand for 3d vector
using Vec3 = Eigen::Matrix<index_t, 3, 1>;
/**
* \brief Draw a filled 3d sphere
*
* \param[in,out] dc the DataContainer where the ellipsoid should be drawn in
* \param[in] amplitude the "color" of the sphere and its filling
* \param[in] center the 3d index of where to place the center of the sphere in dc
* \param[in] sizes of the radius (in x/y/z) of the sphere
*
*/
static void drawFilledEllipsoid(DataContainer<data_t>& dc, data_t amplitude, Vec3 center, Vec3 sizes);
static void drawFilledEllipse(DataContainer<data_t>& dc, data_t amplitude, Vec3 center, data_t radiusX, data_t radiusY, index_t z);
private:
static void drawLine(DataContainer<data_t>&dc, data_t amplitude, Vec3 center, index_t x, index_t y, index_t z);
};
}
\ No newline at end of file
......@@ -9,13 +9,13 @@ using namespace elsa;
void example2d()
{
// generate 2d phantom
IndexVector_t size(2);
size << 128, 128;
IndexVector_t size(3);
size << 128, 128, 128;
auto phantom = PhantomGenerator<real_t>::createModifiedSheppLogan(size);
auto& volumeDescriptor = phantom.getDataDescriptor();
// write the phantom out
EDF::write(phantom, "2dphantom.edf");
MHD::write(phantom, "23dphantom.mhd", "23dphantom.raw");
// generate circular trajectory
index_t numAngles{180}, arc{360};
......@@ -34,7 +34,7 @@ void example2d()
auto sinogram = projector.apply(phantom);
// write the sinogram out
EDF::write(sinogram, "2dsinogram.edf");
MHD::write(sinogram, "23dsinogram.mhd", "23dsinogram.raw");
// setup reconstruction problem
WLSProblem problem(projector, sinogram);
......@@ -48,7 +48,7 @@ void example2d()
auto reconstruction = solver.solve(noIterations);
// write the reconstruction out
EDF::write(reconstruction, "2dreconstruction.edf");
MHD::write(reconstruction, "23dreconstruction.mhd", "23dreconstruction.raw");
}
int main()
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment