User talk:Jasminium

Welcome to Openrtk! We hope you will contribute much and well. You will probably want to read the help pages. Again, welcome and have fun! Matthew.bowman (talk) 08:45, 25 November 2014 (EST)

Forward Projection (RayCast vs RayBox)

Hi I have a question regarding forward projections. I tried to reconcile RayBox and RayCast Forward Projections but got wrong results from RayCast as soon as I introduce volume rotation. I put the projection output and code below, since projection is analytic, they should coincide.

Has anybody got an idea what went wrong?

Thanks a lot Steve

Projection from RayCast

z: 0

              0:           1:           2:           3:
 0:     166.6788     176.3021     176.3021     166.6788 
 1:     166.0071     175.5859     175.5859     166.0071 
 2:     166.0071     175.5859     175.5859     166.0071 
 3:     166.6788     176.3021     176.3021     166.6788 

Projection from RayBox

z: 0

              0:           1:           2:           3:
 0:      178.787     240.9991     240.9991      178.787 
 1:     178.0665     240.0199     240.0199     178.0665 
 2:     178.0665     240.0199     240.0199     178.0665 
 3:      178.787     240.9991     240.9991      178.787 

void test_forward_projection() {

 const unsigned int Dimension = 3;
 typedef double                                    OutputPixelType;
 typedef itk::Image< OutputPixelType, Dimension > OutputImageType;
 typedef itk::Vector<double, 3>                   VectorType;
 int nProj = 1;

double angular_arc=360; double first_angle=45;

 // Constant image sources
 typedef rtk::ConstantImageSource< OutputImageType > ConstantImageSourceType;
 ConstantImageSourceType::PointType origin;
 ConstantImageSourceType::SizeType size;
 ConstantImageSourceType::SpacingType spacing;

 // The test projects a volume filled with ones. The forward projector should
 // then return the intersection of the ray with the box and it is compared
 // with the analytical intersection of a box with a ray.

 // Create volume input.
 const ConstantImageSourceType::Pointer volInput = ConstantImageSourceType::New();
 origin[0] = -96;
 origin[1] = -96;
 origin[2] = -96;

 size[0] = 4;
 size[1] = 4;
 size[2] = 4;
 spacing[0] = 64.;
 spacing[1] = 64.;
 spacing[2] = 64.;

 volInput->SetOrigin( origin );
 volInput->SetSpacing( spacing );
 volInput->SetSize( size );
 volInput->SetConstant( 1. );
 volInput->UpdateOutputInformation();

 // Initialization Volume, it is used in the forward projector and in the
 // Ray Box Intersection Filter in order to initialize the stack of projections.
 const ConstantImageSourceType::Pointer projInput = ConstantImageSourceType::New();
 size[2] = nProj;
 projInput->SetOrigin( origin );
 projInput->SetSpacing( spacing );
 projInput->SetSize( size );
 projInput->SetConstant( 0. );
 projInput->Update();

 // Forward Projection filter
 typedef rtk::RayCastInterpolatorForwardProjectionImageFilter<OutputImageType, OutputImageType> FPType;
 FPType::Pointer fp = FPType::New();
 fp->InPlaceOff();
 fp->SetInput( projInput->GetOutput() );
 fp->SetInput( 1, volInput->GetOutput() );

 // Ray Box Intersection filter (reference)
 typedef rtk::RayBoxIntersectionImageFilter<OutputImageType, OutputImageType> RBIType;
 RBIType::Pointer rbi = RBIType::New();
 rbi->InPlaceOff();
 rbi->SetInput( projInput->GetOutput() );
 VectorType boxMin, boxMax;
 boxMin[0] = -96.0;
 boxMin[1] = -96.0;
 boxMin[2] = -96.0;
 boxMax[0] =  96.0;
 boxMax[1] =  96.0;
 boxMax[2] =  96.0;
 rbi->SetBoxMin(boxMin);
 rbi->SetBoxMax(boxMax);

double dAngle = angular_arc / double(nProj);

 // Geometry
 typedef rtk::ThreeDCircularProjectionGeometry GeometryType;
 GeometryType::Pointer geometry = GeometryType::New();
 for( int i=0; i < nProj; i++)
{
    double angle = first_angle + double(i) * dAngle;
    // cout << angle << endl;
   geometry->AddProjection(500., 1000., -angle);
  }

 rbi->SetGeometry( geometry );
 rbi->Update();

 fp->SetGeometry( geometry );
 fp->Update();

 CheckImageQuality<OutputImageType>(rbi->GetOutput(), fp->GetOutput());

// std::cout << "\n\nTest PASSED! " << std::endl;

 WriteITK(fp->GetOutput());
 WriteITK(rbi->GetOutput());

}