User talk:Jasminium: Difference between revisions

From Openrtk
Jump to navigation Jump to search
(Replaced content with "'''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! [[Us...")
 
Line 3: Line 3:
You will probably want to read the [[Help:Contents|help pages]].
You will probably want to read the [[Help:Contents|help pages]].
Again, welcome and have fun! [[User:Matthew.bowman|Matthew.bowman]] ([[User talk:Matthew.bowman|talk]]) 08:45, 25 November 2014 (EST)
Again, welcome and have fun! [[User:Matthew.bowman|Matthew.bowman]] ([[User talk: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());
}

Latest revision as of 03:36, 26 November 2014

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)

Retrieved from "http://wiki.openrtk.org/index.php?title=User_talk:Jasminium&oldid=276"