Fix #32 by computing grid size once for native aligned grid, then copying the…

Fix #32 by computing grid size once for native aligned grid, then copying the grid geometry into ITK container

src/plastimatch/base/bspline_header.cxx

@@ -61,7 +61,6 @@ Bspline_header::set (
     this->num_coeff = this->cdims[0] * this->cdims[1] * this->cdims[2] * 3;
 }
 
-
 void
 Bspline_header::set (
     const Plm_image_header *pih,
@@ -96,3 +95,70 @@ Bspline_header::set (
     this->set (img_origin, img_spacing, img_dim,
         roi_offset, roi_dim, vox_per_rgn, direction_cosines);
 }
+
+void
+Bspline_header::set_unaligned (
+    const float img_origin[3],
+    const float img_spacing[3],
+    const plm_long img_dim[3],
+    const plm_long roi_offset[3],
+    const plm_long roi_dim[3],
+    const float grid_spac[3],
+    const float direction_cosines[9]
+)
+{
+    this->dc.set (direction_cosines);
+    for (int d = 0; d < 3; d++) {
+        /* copy input parameters over */
+        this->img_origin[d] = img_origin[d];
+        this->img_spacing[d] = img_spacing[d];
+        this->img_dim[d] = img_dim[d];
+        this->roi_offset[d] = roi_offset[d];
+        this->roi_dim[d] = roi_dim[d];
+
+        /* vox_per_rgn is unused for unaligned grids */
+        this->vox_per_rgn[d] = 0;
+
+        /* rdims is the number of regions */
+        float img_ext = (img_dim[d] - 1) * fabs (img_spacing[d]);
+        this->rdims[d] = 1 + (int) floor (img_ext / grid_spac[d]);
+
+        /* cdims is the number of control points */
+        this->cdims[d] = 3 + this->rdims[d];
+        
+    }
+
+    /* total number of control points & coefficients */
+    this->num_knots = this->cdims[0] * this->cdims[1] * this->cdims[2];
+    this->num_coeff = this->cdims[0] * this->cdims[1] * this->cdims[2] * 3;
+}
+
+void
+Bspline_header::set_unaligned (
+    const Plm_image_header *pih,
+    const float grid_spac[3]
+)
+{
+    float img_origin[3];
+    float img_spacing[3];
+    plm_long img_dim[3];
+    plm_long roi_offset[3];
+    plm_long roi_dim[3];
+    plm_long vox_per_rgn[3];
+    float direction_cosines[9];
+
+    pih->get_origin (img_origin);
+    pih->get_dim (img_dim);
+    pih->get_spacing (img_spacing);
+    pih->get_direction_cosines (direction_cosines);
+
+    for (int d = 0; d < 3; d++) {
+        /* Old ROI was whole image */
+        roi_offset[d] = 0;
+        roi_dim[d] = img_dim[d];
+    }
+    
+    this->set_unaligned (
+        img_origin, img_spacing, img_dim, roi_offset, roi_dim,
+        grid_spac, direction_cosines);
+}

src/plastimatch/base/bspline_header.h

@@ -58,9 +58,10 @@ public:
         const plm_long vox_per_rgn[3],
         const float direction_cosines[9]
     );
+
     /** Set the internal geometry.  
-        This version of the function gets used when creating a B-Spline 
-        with a specified grid spacing. 
+        This version of the function gets used to calculate the 
+        geometry of an ALIGNED B-Spline with a specified grid spacing. 
         \param pih The image geometry associated with B-spline
         \param grid_spac The B-Spline grid spacing (in mm)
     */
@@ -68,6 +69,33 @@ public:
         const Plm_image_header *pih,
         const float grid_spac[3]
     );
+
+    /** Set the internal geometry.  
+        This version of the function gets used to calculate the 
+        geometry of an UNALIGNED B-Spline with a specified grid spacing. 
+        \param pih The image geometry associated with B-spline
+        \param grid_spac The B-Spline grid spacing (in mm)
+    */
+    void set_unaligned (
+        const float img_origin[3],
+        const float img_spacing[3],
+        const plm_long img_dim[3],
+        const plm_long roi_offset[3],
+        const plm_long roi_dim[3],
+        const float grid_spac[3],
+        const float direction_cosines[9]
+    );
+
+    /** Set the internal geometry.  
+        This version of the function gets used to calculate the 
+        geometry of an UNALIGNED B-Spline with a specified grid spacing. 
+        \param pih The image geometry associated with B-spline
+        \param grid_spac The B-Spline grid spacing (in mm)
+    */
+    void set_unaligned (
+        const Plm_image_header *pih,
+        const float grid_spac[3]
+    );
     void get_volume_header (Volume_header *vh);
 };
 

src/plastimatch/base/xform.cxx

@@ -836,10 +836,10 @@ xform_aff_to_itk_bsp_bulk (
 /* Convert xf to vector field to bspline */
 static void
 xform_any_to_itk_bsp_nobulk (
-    Xform *xf_out, Xform* xf_in,
-    const Plm_image_header* pih,
-    const float* grid_spac)
+    Xform *xf_out, Xform *xf_in,
+    const Bspline_header *bh)
 {
+#if defined (commentout)
     int d;
     Xform xf_tmp;
     Plm_image_header pih_bsp;
@@ -849,11 +849,6 @@ xform_any_to_itk_bsp_nobulk (
     itk_bsp_set_grid_img (xf_out, pih, grid_spac);
     BsplineTransformType::Pointer bsp_out = xf_out->get_itk_bsp();
 
-    /* Create temporary array for output coefficients */
-    const unsigned int num_parms = bsp_out->GetNumberOfParameters();
-    BsplineTransformType::ParametersType bsp_coeff;
-    bsp_coeff.SetSize (num_parms);
-
     /* Compute bspline grid specifications */
     BsplineTransformType::OriginType bsp_origin;
     BsplineTransformType::SpacingType bsp_spacing;
@@ -861,18 +856,47 @@ xform_any_to_itk_bsp_nobulk (
     BsplineTransformType::DirectionType bsp_direction;
     bsp_grid_from_img_grid (bsp_origin, bsp_spacing,
         bsp_region, bsp_direction, pih, grid_spac);
+#endif
 
+    /* Create ITK B-Spline container */
+    xform_itk_bsp_init_default (xf_out);
+    BsplineTransformType::Pointer bsp_out = xf_out->get_itk_bsp();
+
+    /* Copy geometry and allocate coefficients */
+    BsplineTransformType::OriginType bsp_origin;
+    BsplineTransformType::SpacingType bsp_spacing;
+    BsplineTransformType::RegionType bsp_region;
+    BsplineTransformType::DirectionType bsp_direction;
+    BsplineTransformType::RegionType::SizeType bsp_size;
+    for (int d=0; d<3; d++) {
+        bsp_origin[d] = bh->img_origin[d];
+        bsp_spacing[d] = bh->grid_spac[d];
+        bsp_size[d] = bh->cdims[d];
+    }
+    bsp_region.SetSize (bsp_size);
+    itk_direction_from_dc (&bsp_direction, bh->dc);
+    xf_out->itk_bsp_set_grid (bsp_origin, bsp_spacing, bsp_region,
+        bsp_direction);
+    
     /* Make a vector field at bspline grid spacing */
+    Plm_image_header pih_bsp;
     pih_bsp.set (bsp_region, bsp_origin, bsp_spacing, bsp_direction);
+    Xform xf_tmp;
     xform_to_itk_vf (&xf_tmp, xf_in, &pih_bsp);
 
-    /* Vector field is interleaved.  We need planar for decomposition. */
+    /* Vector field is interleaved.  Create temporary images to hold the 
+       non-interleaved images of coefficients */
     FloatImageType::Pointer img = itk_image_create<float> (pih_bsp);
 
+    /* Create temporary array for storing output coefficients */
+    const unsigned int num_parms = bsp_out->GetNumberOfParameters();
+    BsplineTransformType::ParametersType bsp_coeff;
+    bsp_coeff.SetSize (num_parms);
+
     /* Loop through planes */
     unsigned int counter = 0;
     DeformationFieldType::Pointer vf = xf_tmp.get_itk_vf();
-    for (d = 0; d < 3; d++) {
+    for (int d = 0; d < 3; d++) {
         /* Copy a single VF plane into img */
         typedef itk::ImageRegionIterator< FloatImageType > FloatIteratorType;
         typedef itk::ImageRegionIterator< DeformationFieldType > VFIteratorType;
@@ -910,6 +934,18 @@ xform_any_to_itk_bsp_nobulk (
     bsp_out->SetParametersByValue (bsp_coeff);
 }
 
+/* This function converts to itk_bsp using an unaligned grid */
+static void
+xform_any_to_itk_bsp_nobulk (
+    Xform *xf_out, Xform* xf_in,
+    const Plm_image_header* pih,
+    const float* grid_spac)
+{
+    Bspline_header bh;
+    bh.set_unaligned (pih, grid_spac);
+    xform_any_to_itk_bsp_nobulk (xf_out, xf_in, &bh);
+}
+
 /* This function extends the B-spline grid, padding with zeros,
    so that the grid contains the specified Region "roi."
    This is sometimes needed so that the B-spline can warp an
@@ -1426,7 +1462,7 @@ xform_any_to_gpuit_bsp (
 
         /* Create itk_bsp xf using image specifications */
         Xform xf_tmp;
-        xform_any_to_itk_bsp_nobulk (&xf_tmp, xf_in, pih, bxf_new->grid_spac);
+        xform_any_to_itk_bsp_nobulk (&xf_tmp, xf_in, bxf_new);
 
         /* Copy from ITK coefficient array to gpuit coefficient array */
         int k = 0;