[Insight-developers] help rendering spatial objects

[Richard Beare]

Hi,

This is the first time I've attempted to use spatial objects and I
suspect I'm going about it slightly the wrong way.

My task is relatively simple - I have a 3D binary object and I want to
replace each slice with its best fitting ellipse. My current approach
is to extract each slice, feed it through the ImageMomentsCalculator
to compute the major and minor axes lengt and angle and create the
corresponding ellipse spatial object. I'm rendering this object to a
2D image and copying the region to the corresponding slice in the 3D
image. The transforms I'm applying seem to be producing the correct
size and angle ellipse (verified by saving the slices). However I run
into problems when I have a rotation that pushes the bounding box
corners out of the bounds of the image into which I want to copy the
result - it seems that the ComputeBoundingBox function returns the
bounds of the transformed bounding box, rather than computing the
bounding box of the transformed ellipse.

Currently I'm using the bounding box to define a region in each image
to loop over, copying values back to the 3D image. However this
becomes extremely messy if one of the regions is no longer inside the
corresponding image. These problems make me think I'm going about it
the wrong way, particularly in my use of spatial object transforms. Do
any errors leap out?

Thanks

Here's my ellipse manipulation code:

    ell->SetRadius(erad);
    if (u11 != 0)
      ell->GetObjectToParentTransform()->Rotate2D(theta, true);
    ell->GetObjectToParentTransform()->SetOffset(moment->GetCenterOfGravity());
    ell->ComputeObjectToWorldTransform();
    ell->ComputeBoundingBox();
    typename EllipseType::BoundingBoxType * ebox = ell->GetBoundingBox();

//   rendering code

    typename MomentType::VectorType boxCorner;

    typename SliceType::SpacingType spacing = slicer->GetOutput()->GetSpacing();
    typename SliceType::SizeType bbsize;
    for (unsigned k = 0; k < SliceType::ImageDimension; k++)
      {
      bbsize[k]= (long)((ebox->GetBounds()[2*k+1] -
ebox->GetBounds()[2*k])/ spacing[k]) + 2;
      boxCorner[k]=ebox->GetBounds()[2*k];
      }
    typename SliceType::PointType COGp;
    COGp.Fill(0);
    COGp += boxCorner;

    // render this into our slice
    renderer->SetInput(ell);
    //renderer->SetSize(slicer->GetOutput()->GetLargestPossibleRegion().GetSize());
    renderer->SetSize(bbsize);
    renderer->SetSpacing(slicer->GetOutput()->GetSpacing());
    renderer->SetOrigin(COGp);
    renderer->Update();

    // now copy the results into output marker
    typename ImType::RegionType WritingRegion3D;
    typename SliceType::RegionType SliceRegion =
renderer->GetOutput()->GetLargestPossibleRegion();

    typename ImType::RegionType::SizeType sz;
    typename ImType::RegionType::IndexType Orig;
    typename SliceType::IndexType SlIdx;
    slicer->GetOutput()->TransformPhysicalPointToIndex(COGp, SlIdx);
    std::cout << "SlIdx " << SlIdx << std::endl;
    writeIm<SliceType>(renderer->GetOutput(), "/tmp/rend.tif");
    Orig[0]=(long)round(SlIdx[0]);
    Orig[1]=(long)round(SlIdx[1]);
    Orig[2]=ExtractionRegion.GetIndex()[2];

    sz[0]=SliceRegion.GetSize()[0];
    sz[1]=SliceRegion.GetSize()[1];
    sz[2]=ExtractionRegion.GetSize()[2];

    WritingRegion3D.SetSize(sz);
    WritingRegion3D.SetIndex(Orig);
    std::cout << SliceRegion << std::endl;
    std::cout << WritingRegion3D << std::endl;

    It3DType it3(outputmarker, WritingRegion3D);
    It2DType it2(renderer->GetOutput(), SliceRegion);
    for (it3.GoToBegin(), it2.GoToBegin(); !it2.IsAtEnd(); ++it3, ++it2)
      {
      typename SliceType::PixelType V = it2.Get();
      if (V)
	{
	it3.Set(V);
	}
      }
    }