Temp check in

Author: cindygr

Image_based/HandleFileNames.py

@@ -193,6 +193,22 @@ def get_image_name(self, path, index, b_add_tag=True):
 
         return im_name
 
+    def get_edge_image_name(self, path, index, b_add_tag=True):
+        """ Get the image name corresponding to the index given by (subdirectory index, image index, -)
+        @param path should be one of self.path, self.path_calculated, or path_debug
+        @param index (tuple, either 2 dim or 3 dim, index into sorted lists)
+        @param b_add_tag - add the image tag, y/n
+        @return full image name with path"""
+
+        im_name = path
+        if len(self.sub_dirs[index[0]]) > 0:
+            im_name = im_name + self.sub_dirs[index[0]] + "/"
+        im_name = im_name + self.image_names[index[0]][index[1]] + "_edge"
+        if b_add_tag:
+            im_name = im_name + self.image_tag
+
+        return im_name
+
     def get_mask_name(self, path, index, b_add_tag=True):
         """ Get the mask name corresponding to the index given by (subdirectory index, image index, mask name, mask id)
         @param path should be one of self.path, self.path_calculated, or path_debug

Image_based/bezier_cyl_2d.py

@@ -47,12 +47,13 @@ def __init__(self, start_pt=None, end_pt=None, radius=1, mid_pt=None):
         else:
             self.p1 = np.array(mid_pt)
         self.start_radius = radius
-        self.end_radius = 100
-    #want to pass every single point from start to end and return radius along that path
+        self.end_radius = radius
+
     def radius(self, t):
-        return (1-t)*(self.start_radius) + (t*self.end_radius)
-        
-    
+        """ Radius is a linear interpolation of two end radii
+        @param t - t between 0 and 1"""
+        return (1 - t) * self.start_radius + t * self.end_radius
+
     @staticmethod
     def _orientation(start_pt, end_pt):
         """Set the orientation and ensure left-right or down-up
@@ -103,22 +104,22 @@ def edge_pts(self, t):
         right_pt = [pt[0] + vec_step[1], pt[1] - vec_step[0]]
         return left_pt, right_pt
 
-    def edge_offset_pt(self, t, perc_in_out, dir):
+    def edge_offset_pt(self, t, perc_in_out, direction):
         """ Go in/out of the edge point a given percentage
         @param t - t value along the curve (in range 0, 1)
         @param perc_in_out - if 1, get point on edge. If 0.5, get halfway to centerline. If 2.0 get 2 width
-        @param dir - 'Left' is the left direction, 'Right' is the right direction
+        @param direction - 'Left' is the left direction, 'Right' is the right direction
         @return numpy array x,y """
         pt_edge = self.pt_axis(t)
         vec_tang = self.tangent_axis(t)
         vec_step = perc_in_out * self.radius(t) * vec_tang / np.sqrt(vec_tang[0] * vec_tang[0] + vec_tang[1] * vec_tang[1])
 
-        if dir == "Left":
+        if direction == "Left":
             return np.array([pt_edge[0] + vec_step[1], pt_edge[1] - vec_step[0]])
         return np.array([pt_edge[0] - vec_step[1], pt_edge[1] + vec_step[0]])
 
     @staticmethod
-    def _rect_in_image(im, r, pad=2):
+    def rect_in_image(im, r, pad=2):
         """ See if the rectangle is within the image boundaries
         @im - image (for width and height)
         @r - the rectangle
@@ -233,7 +234,7 @@ def interior_rects_mask(self, image_shape, step_size=40, perc_width=0.3):
         return ret_im_mask
 
     @staticmethod
-    def _image_cutout(im, rect, step_size, height):
+    def image_cutout(im, rect, step_size, height):
         """Cutout a warped bit of the image and return it
         @param im - the image rect is in
         @param rect - four corners of the rectangle to cut out
@@ -430,9 +431,9 @@ def make_mask_image(self, im_mask, step_size=20, perc_fuzzy=0.2):
         @param perc_fuzzy How much of the boundary to make fuzzy
         """
         self.draw_interior_rects_filled(im_mask, b_solid=True,
-                                                 col_solid=(255, 255, 255),
-                                                 step_size=step_size,
-                                                 perc_width=1.0)
+                                        col_solid=(255, 255, 255),
+                                        step_size=step_size,
+                                        perc_width=1.0)
         self.draw_boundary_rects_filled(im_mask, b_solid=True,
                                         col_solid=(128, 128, 128),
                                         step_size=step_size,

Image_based/extract_curves.py

@@ -8,15 +8,12 @@
 #     - Store as t, percentage in/out from radius (so we can re-use at different scales)
 
 import numpy as np
-from glob import glob
-import csv
 import cv2
 import json
 from os.path import exists
-from bezier_cyl_2d import BezierCyl2D
 from line_seg_2d import LineSeg2D
-from fit_bezier_cyl_2d_edge import FitBezierCyl2DEdge
 from HandleFileNames import HandleFileNames
+from fit_bezier_cyl_2d_edge import FitBezierCyl2DEdge
 
 
 class ExtractCurves:
@@ -53,7 +50,7 @@ def __init__(self, fname_rgb_image, fname_edge_image, fname_mask_image, params=N
         if b_recalc or not fname_calculated or not exists(self.fname_full_edge_stats):
             # Recalculate and write
             self.edge_stats = ExtractCurves.full_edge_stats(self.bezier_edge.image_edge,
-                                                            self.bezier_edge.fit_bezier_edge,
+                                                            self.bezier_edge.bezier_crv_fit_to_edge,
                                                             self.params)
             # Write out the bezier curve
             if fname_calculated:
@@ -83,30 +80,34 @@ def __init__(self, fname_rgb_image, fname_edge_image, fname_mask_image, params=N
 
         if fname_debug:
             # Draw the mask with the initial and fitted curve
-            im_covert_back = cv2.cvtColor(self.bezier_edge.edge_image, cv2.COLOR_GRAY2RGB)
-            im_rgb = cv2.copyTo(self.bezier_edge.rgb_image)
+            im_covert_back = cv2.cvtColor(self.bezier_edge.image_edge, cv2.COLOR_GRAY2RGB)
+            im_rgb = np.copy(self.bezier_edge.image_rgb)
             for pix in self.edge_stats["pixs_edge"]:
                 im_covert_back[pix[0], pix[1], :] = (255, 0, 0)
                 im_rgb[pix[0], pix[1], :] = (255, 255, 255)
 
             for do_both_crv, do_both_name in [(self.left_curve, "Left"), (self.right_curve, "Right")]:
                 for pt_e1, pt_e2 in zip(do_both_crv[0:-1], do_both_crv[1:]):
-                    pt1 = self.bezier_edge.bezier_fit.edge_offset_pt(pt_e1[0], pt_e1[1], do_both_name)
-                    pt2 = self.bezier_edge.bezier_fit.edge_offset_pt(pt_e2[0], pt_e2[1], do_both_name)
+                    pt1 = self.bezier_edge.bezier_crv_fit_to_edge.edge_offset_pt(pt_e1[0], pt_e1[1], do_both_name)
+                    pt2 = self.bezier_edge.bezier_crv_fit_to_edge.edge_offset_pt(pt_e2[0], pt_e2[1], do_both_name)
                     LineSeg2D.draw_line(im_covert_back, pt1, pt2, color=(255, 255, 0))
                     LineSeg2D.draw_line(im_rgb, pt1, pt2, color=(255, 255, 0))
             im_both = np.hstack([im_covert_back, im_rgb])
             cv2.imwrite(fname_debug, im_both)
 
-    @staticmethod def full_edge_stats(image_edge, bezier_edge, params):
+    @staticmethod
+    def full_edge_stats(image_edge, bezier_edge, params):
         """ Get the best pixel offset (if any) for each point along the edge
         @param image_edge - the edge image
         @param bezier_edge - the Bezier curve
         @param params - parameters for extraction"""
         bdry_rects1, ts1 = bezier_edge.boundary_rects(step_size=params["step_size"], perc_width=params["perc_width"])
         bdry_rects2, ts2 = bezier_edge.boundary_rects(step_size=params["step_size"], perc_width=params["perc_width"], offset=True)
         n_bdry1 = len(bdry_rects1)
-        t_step = (ts1[-1] - ts1[0]) / (len(bdry_rects1) // 2)
+        try:
+            t_step = ts1[2] - ts1[0]
+        except IndexError:
+            t_step = 1.0
 
         bdry_rects1.extend(bdry_rects2)
         ts1.extend(ts2)
@@ -120,22 +121,30 @@ def __init__(self, fname_rgb_image, fname_edge_image, fname_mask_image, params=N
         for i_rect, r in enumerate(bdry_rects1):
             # b_rect_inside = BezierCyl2D._rect_in_image(image_edge, r, pad=2)
 
-            im_warp, tform3_back = bezier_edge._image_cutout(image_edge, r, step_size=width, height=height)
+            im_warp, tform3_back = bezier_edge.image_cutout(image_edge, r, step_size=width, height=height)
             # Actual hough transform on the cut-out image
             lines = cv2.HoughLines(im_warp, 1, np.pi / 180.0, 10)
 
+            # Check for any lines in the cutout image
+            if lines is None:
+                continue
+            # .. and check if any of those are horizontal
+            ret_pts = FitBezierCyl2DEdge.get_horizontal_lines_from_hough(lines, tform3_back, width, height)
+            if ret_pts is []:
+                continue
             i_side = i_rect % 2
 
             if i_rect == n_bdry1:
-                t_step = (ts2[-1] - ts2[0]) / (len(bdry_rects2) // 2)
+                try:
+                    t_step = ts2[2] - ts2[0]
+                except IndexError:
+                    t_step = 1.0
 
-            if not lines:
-                continue
             max_y = im_warp.max(axis=0)
 
-            ts_seg = np.linspace(ts1[i_rect] - t_step * 0.5, ts1[i_rect] - t_step * 0.5, len(max_y))
+            ts_seg = np.linspace(ts1[i_rect] - t_step * 0.5, ts1[i_rect] + t_step * 0.5, len(max_y))
             for i_col, y in enumerate(max_y):
-                if max_y > params["edge_max"]:
+                if y > params["edge_max"]:
                     ids = np.where(im_warp[:, i_col] == y)
                     if i_side == 0:
                         tag = "left"
@@ -146,13 +155,17 @@ def __init__(self, fname_rgb_image, fname_edge_image, fname_mask_image, params=N
                     h_max = (1 + params["perc_width"]) * bezier_edge.radius(ts1[i_col])
                     ret_stats[tag + "_perc"].append(h_min + h_max + ids[0] / width)
 
-                    p1_in = np.transpose(np.array([ids[0], i_col, 1.0]))
+                    p1_in = np.transpose(np.array([ids[0][0], i_col, 1.0]))
                     p1_back = tform3_back @ p1_in
 
                     ret_stats["pixs_edge"].append([p1_back[0], p1_back[1]])
 
-        sort(zip(ret_stats["ts_left"], ret_stats["left_perc"]), key=0)
-        sort(zip(ret_stats["ts_right"], ret_stats["right"]), key=0)
+        np.array([ret_stats["ts_left"], ret_stats["left_perc"]])
+        np.sort(axis=0)
+        ret_stats["ts_left"] = list(np.array[0, :])
+        ret_stats["left_perc"] = list(np.array[1, :])
+        # sort(zip(ret_stats["ts_left"], ret_stats["left_perc"]), key=0)
+        # sort(zip(ret_stats["ts_right"], ret_stats["right"]), key=0)
         return ret_stats
 
     @staticmethod
@@ -168,35 +181,37 @@ def curves_from_stats(stats_edge, params):
         for ts, ps in [(stats_edge["ts_left"], stats_edge["Left_perc"]), (stats_edge["ts_right"], stats_edge["Right_perc"])]:
             ps_filter = np.array(ps)
             for i_filter in range(0, params["n_filter"]):
-                np.filter(ps_filter)
+                # np.filter(ps_filter)
                 ts_crvs = np.linspace(0, 1, params["n_samples"])
                 ps_crvs = np.interp(ts_crvs, ts, ps_filter)
                 crvs.append([(t, p) for t, p in zip(ts_crvs, ps_crvs)])
         return crvs[0], crvs[1]
 
 
 if __name__ == '__main__':
-    #path_bpd = "./data/trunk_segmentation_names.json"
+    # path_bpd = "./data/trunk_segmentation_names.json"
     path_bpd = "./data/forcindy_fnames.json"
     all_files = HandleFileNames.read_filenames(path_bpd)
 
     b_do_debug = True
-    b_do_recalc = False
+    b_do_recalc = True
     for ind in all_files.loop_masks():
-        mask_fname = all_files.get_mask_name(path=all_files.path, index=ind, b_add_tag=True)
         rgb_fname = all_files.get_image_name(path=all_files.path, index=ind, b_add_tag=True)
-        edge_fname = all_files.get_image_name(path=all_files.path_calculated, index=ind, b_add_tag=False) + "_edge.png"
-        mask_fname_debug = all_files.get_mask_name(path=all_files.path_debug, index=ind, b_add_tag=False)
+        edge_fname = all_files.get_edge_image_name(path=all_files.path_calculated, index=ind, b_add_tag=True)
+        mask_fname = all_files.get_mask_name(path=all_files.path, index=ind, b_add_tag=True)
+        ec_fname_debug = all_files.get_mask_name(path=all_files.path_debug, index=ind, b_add_tag=False)
         if not b_do_debug:
-            mask_fname_debug = ""
+            ec_fname_debug = ""
         else:
-            mask_fname_debug = mask_fname_debug + "_crv.png"
+            ec_fname_debug = ec_fname_debug + "_extract_profile.png"
 
-        mask_fname_calculate = all_files.get_mask_name(path=all_files.path_calculated, index=ind, b_add_tag=False)
+        ec_fname_calculate = all_files.get_mask_name(path=all_files.path_calculated, index=ind, b_add_tag=False)
 
         if not exists(mask_fname):
             raise ValueError(f"Error, file {mask_fname} does not exist")
         if not exists(rgb_fname):
             raise ValueError(f"Error, file {rgb_fname} does not exist")
 
-        b_stats = ExtractCurves(rgb_fname, edge_fname, mask_fname, mask_fname_calculate, mask_fname_debug, b_recalc=b_do_recalc)
+        profile_crvs = ExtractCurves(rgb_fname, edge_fname, mask_fname,
+                                     fname_calculated=ec_fname_calculate,
+                                     fname_debug=ec_fname_debug, b_recalc=b_do_recalc)

Image_based/fit_bezier_cyl_2d_edge.py

@@ -41,7 +41,7 @@ def __init__(self, fname_mask_image, fname_calculated=None, fname_debug=None, b_
         #   This is the curve that will be fit to the mask
         self.bezier_crv_fit_to_mask = FitBezierCyl2D(self.bezier_crv_initial)
 
-        # Fit a quad to the mask, using the end points of the base image as a starting point
+        # Create the calculated file names
         print(f"Fitting bezier curve to mask image {fname_mask_image}")
         self.fname_bezier_cyl_initial = None    # The actual quadratic bezier
         self.fname_bezier_cyl_fit_to_mask = None    # The actual quadratic bezier
@@ -98,7 +98,7 @@ def __init__(self, fname_mask_image, fname_calculated=None, fname_debug=None, b_
             cv2.imwrite(fname_debug, im_both)
 
         self.score = self.score_mask_fit(self.stats_dict.mask_image)
-        print(f"Mask {mask_fname}, score {self.score}")
+        print(f"Mask {fname_mask_image}, score {self.score}")
 
     @staticmethod
     def create_bezier_crv_from_eigen_vectors(stats):
@@ -131,17 +131,17 @@ def _adjust_bezier_crv_by_mask(fit_bezier_crv, im_mask, step_size=40, perc_width
         @param step_size - how many pixels to step along
         @param perc_width - how much wider than the radius to look in the mask
         @returns how much the points moved"""
-        height = int(fit_bezier_crv.radius_2d)
+        height = int(fit_bezier_crv.radius(0.5))
         rects, ts = fit_bezier_crv.interior_rects(step_size=step_size, perc_width=perc_width)
 
         # Set up the matrix - include the 3 current points plus the centers of the mask
         a_constraints, b_rhs = fit_bezier_crv.setup_least_squares(ts)
 
         x_grid, y_grid = np.meshgrid(range(0, step_size), range(0, height))
         for i, r in enumerate(rects):
-            b_rect_inside = BezierCyl2D._rect_in_image(im_mask, r, pad=2)
+            b_rect_inside = BezierCyl2D.rect_in_image(im_mask, r, pad=2)
 
-            im_warp, tform_inv = fit_bezier_crv._image_cutout(im_mask, r, step_size=step_size, height=height)
+            im_warp, tform_inv = fit_bezier_crv.image_cutout(im_mask, r, step_size=step_size, height=height)
             if b_rect_inside and np.sum(im_warp > 0) > 0:
                 x_mean = np.mean(x_grid[im_warp > 0])
                 y_mean = np.mean(y_grid[im_warp > 0])
@@ -182,17 +182,12 @@ def score_mask_fit(self, im_mask):
         # edges
         im_bezier_mask = np.zeros((self.stats_dict.mask_image.shape[0], self.stats_dict.mask_image.shape[1]), np.uint8)
 
-        self.bezier_crv_fit_to_mask.make_mask_image(im_bezier_mask,
-                                                     step_size=10,
-                                                     perc_fuzzy=0.25)
+        self.bezier_crv_fit_to_mask.make_mask_image(im_bezier_mask, step_size=10, perc_fuzzy=0.25)
 
         pixs_in_mask_not_bezier = np.logical_and(im_mask > 0, im_bezier_mask == 0)
         pixs_in_bezier_not_mask = np.logical_and(im_mask == 0, im_bezier_mask == 255)
         pixs_in_both_masks = np.logical_and(im_mask > 0, im_bezier_mask == 255)
         pixs_in_union = np.logical_or(np.logical_or(pixs_in_bezier_not_mask, pixs_in_mask_not_bezier), pixs_in_both_masks)
-        n_in_mask_not_bezier = np.count_nonzero(pixs_in_mask_not_bezier)
-        n_in_bezier_not_mask = np.count_nonzero(pixs_in_bezier_not_mask)
-        n_in_both_masks = np.count_nonzero(pixs_in_both_masks)
         n_in_union = np.count_nonzero(pixs_in_union)
         if n_in_union == 0:
             return 0
@@ -201,7 +196,7 @@ def score_mask_fit(self, im_mask):
 
 
 if __name__ == '__main__':
-    #path_bpd = "./data/trunk_segmentation_names.json"
+    # path_bpd = "./data/trunk_segmentation_names.json"
     path_bpd = "./data/forcindy_fnames.json"
     all_files = HandleFileNames.read_filenames(path_bpd)
 

Image_based/fit_bezier_cyl_2d_mask.py

@@ -24,9 +24,9 @@ def __init__(self, p1=(0.0, 0.0, 0.0), p2=(0.5, 0.75, 0.5), p3=(1.0, 1.0, 1.0),
         @param end_radius - ending radius"""
 
         # Information about the current branch/trunk
-        self.pt1 = p1
-        self.pt2 = p2
-        self.pt3 = p3
+        self.pt1 = np.array(p1)
+        self.pt2 = np.array(p2)
+        self.pt3 = np.array(p3)
         self.start_radii = start_radius
         self.end_radii = end_radius
 
@@ -37,6 +37,24 @@ def __init__(self, p1=(0.0, 0.0, 0.0), p2=(0.5, 0.75, 0.5), p3=(1.0, 1.0, 1.0),
         self.vertex_locs = np.zeros((self.n_along, self.n_around, 3))
         self.make_mesh()
 
+    def copy(self, bezier_crv=None, b_compute_mesh=False):
+        """Return a copy of self - mostly just copy, rather than =, the points/vertices
+        @param bezier_crv - use this bezier curve versus making a new one
+        @param b_compute_mesh - compute/copy the mesh, y/n
+        @return: New curve"""
+        if bezier_crv is None:
+            bezier_crv = BezierCyl3D(np.copy(self.pt1), np.copy(self.pt2), np.copy(self.pt3), self.start_radii, self.end_radii)
+        for k, v in self.__dict__:
+            try:
+                if v.size > 1:
+                    pass
+            except TypeError:
+                setattr(bezier_crv, k, v)
+
+        if b_compute_mesh:
+            bezier_crv.vertex_locs = np.copy(self.vertex_locs)
+        return bezier_crv
+
     def n_vertices(self):
         return self.n_along * self.n_around
 
@@ -51,9 +69,9 @@ def set_pts(self, pt1, pt2, pt3):
         @param pt2 Mid point
         @param pt3 End point
         """
-        self.pt1 = pt1
-        self.pt2 = pt2
-        self.pt3 = pt3
+        self.pt1 = np.array(pt1)
+        self.pt2 = np.array(pt2)
+        self.pt3 = np.array(pt3)
 
     def set_pts_from_pt_tangent(self, pt1, vec1, pt3):
         """Set the points from a starting point/tangent
@@ -132,7 +150,7 @@ def _calc_radii(self):
 
     def _calc_cyl_vertices(self):
         """Calculate the cylinder vertices"""
-        pt = np.ones(shape=(4))
+        pt = np.ones(shape=(4,))
         radii = self._calc_radii()
 
         for it, t in enumerate(np.linspace(0, 1.0, self.n_along)):
@@ -175,6 +193,8 @@ def write_json(self, fname):
         @param fname - file name"""
         fix_nparray = []
         for k, v in self.__dict__.items():
+            if k == "vertex_locs":
+                continue
             try:
                 if v.size == 3:
                     fix_nparray.append([k, v])
@@ -207,7 +227,7 @@ def read_json(fname, bezier_crv=None, compute_mesh=True):
                         setattr(bezier_crv, k, v)
                 except TypeError:
                     setattr(bezier_crv, k, v)
-
+        bezier_crv.set_dims(bezier_crv.n_along, bezier_crv.n_around)
         if compute_mesh:
             bezier_crv.make_mesh()
         return bezier_crv

bezier_cyl_3d.py

@@ -0,0 +1,171 @@
+#!/usr/bin/env python3
+
+# a 3D Bezier branch/trunk
+#  Inherits from 3D bezier curve
+#     - The 3D bezier curve
+#     - Start and end radii
+#  Adds
+#    - placing buds/bud geometry
+#    - profile curve - making the outside "bumpy"
+#
+# Can:
+#   - Evaluate points along the curve
+#   - Make a cylinderical mesh
+#   - Project self into an image
+
+import numpy as np
+from json import load
+from bezier_cyl_3d import BezierCyl3D
+from scipy.spatial.transform import Rotation as R
+
+
+class BezierCyl3DWithDetail(BezierCyl3D):
+    _profiles = None
+    _type_names = {"trunk", "sidebranch", "branch"}
+    _bud_shape = None
+    _data_dir = "./data/"
+
+    def __init__(self, bezier_crv=None):
+        """ Initialize a 3D curve, built from a quadratic Bezier with radii
+        @param bezier_crv - 3D bezier curve to start with"""
+
+        if bezier_crv:
+            super(BezierCyl3DWithDetail, self).__init__()
+            bezier_crv.copy(bezier_crv=self, b_compute_mesh=True)
+        else:
+            super(BezierCyl3DWithDetail, self).__init__()
+
+        self.start_is_junction = False
+        self.end_is_bud = False
+        self.start_bud = 0.7
+        self.bud_angle = 0.8 * np.pi / 2
+        self.bud_length = 0.1
+
+        # Read in global data/profiles
+        BezierCyl3DWithDetail._read_profiles(self._data_dir)
+        BezierCyl3DWithDetail._make_bud_shape()
+
+    @staticmethod
+    def _read_profiles(data_dir):
+        """ Read in all the profile curves for the various branch types"""
+        if BezierCyl3DWithDetail._profiles is not None:
+            return
+
+        BezierCyl3DWithDetail._profiles = {}
+        for t in BezierCyl3DWithDetail._type_names:
+            try:
+                fname = data_dir + "/" + t + "_profiles.json"
+                with open(fname, "r") as fp:
+                    BezierCyl3DWithDetail._profiles[t] = load(fp)
+            except FileNotFoundError:
+                pass
+
+    @staticmethod
+    def _make_bud_shape():
+        if BezierCyl3DWithDetail._bud_shape is None:
+            n_pts = 10
+            BezierCyl3DWithDetail._bud_shape = np.zeros((2, n_pts))
+            BezierCyl3DWithDetail._bud_shape[0, :] = np.linspace(0, 1.0, n_pts)
+            BezierCyl3DWithDetail._bud_shape[0, -2] = 0.5 * BezierCyl3DWithDetail._bud_shape[0, -2] + 0.5 * BezierCyl3DWithDetail._bud_shape[0, -1]
+            BezierCyl3DWithDetail._bud_shape[1, 0] = 1.0
+            BezierCyl3DWithDetail._bud_shape[1, 1] = 0.95
+            BezierCyl3DWithDetail._bud_shape[1, 2] = 1.05
+            BezierCyl3DWithDetail._bud_shape[1, 3] = 1.1
+            BezierCyl3DWithDetail._bud_shape[1, 4] = 1.05
+            BezierCyl3DWithDetail._bud_shape[1, 5] = 0.8
+            BezierCyl3DWithDetail._bud_shape[1, 6] = 0.7
+            BezierCyl3DWithDetail._bud_shape[1, 7] = 0.5
+            BezierCyl3DWithDetail._bud_shape[1, 8] = 0.3
+            BezierCyl3DWithDetail._bud_shape[1, 9] = 0.0
+
+    def _make_cyl(self, profiles):
+        """ Make a 3D generalized cylinder
+        @param profiles - variations to the radii """
+        if profiles:
+            self._calc_cyl_vertices()
+        else:
+            self._calc_cyl_vertices()
+
+    def set_radii_and_junction(self, start_radius=1.0, end_radius=1.0, b_start_is_junction=False, b_end_is_bud=False):
+        """ Set the radius of the branch
+        @param start_radius - radius at pt1
+        @param end_radius - radius at pt3
+        @param b_start_is_junction - is the start of the curve a junction?
+        @param b_end_is_bud - is the end a bud? """
+        self.set_radii(start_radius, end_radius)
+        self.start_is_junction = b_start_is_junction
+        self.end_is_bud = b_end_is_bud
+
+    def make_branch_segment(self, pt1, pt2, pt3, radius_start, radius_end, start_is_junction, end_is_bud):
+        """ Output a 3D generalized cylinder"""
+        self.set_pts(pt1, pt2, pt3)
+        self.set_radii_and_junction(start_radius=radius_start, end_radius=radius_end, b_start_is_junction=start_is_junction, b_end_is_bud=end_is_bud)
+        try:
+            self._make_cyl(BezierCyl3DWithDetail._profiles["sidebranches"])
+        except KeyError:
+            self._make_cyl(None)
+
+    def place_buds(self, locs):
+        """ Position and orientation of buds,
+        @param locs - t along, radius loc tuples in a list
+        @
+        @return [(pt1, pt2, pt3) """
+
+        ts = np.linspace(0, 1, self.n_along)
+        radii = self._calc_radii()
+
+        pt = np.ones((4,))
+        zero_pt = np.ones((4,))
+        zero_pt[0:3] = 0.0
+        vec = np.zeros((4,))
+        ret_list = []
+        for loc in locs:
+            mat = self.frenet_frame(loc[0])
+            r = np.interp(loc[0], ts, radii)
+            pt_on_crv = mat @ zero_pt
+            pt[0] = np.cos(loc[1]) * r
+            pt[1] = np.sin(loc[1]) * r
+            pt[2] = 0
+            pt_on_srf = mat @ pt
+            vec[0] = np.cos(loc[1])
+            vec[1] = np.sin(loc[1])
+            vec[2] = 0
+            vec_rotate = np.cross(vec[0:3], np.array([0, 0, 1]))
+            vec_rotate = vec_rotate / np.linalg.norm(vec_rotate)
+            mat_rotate_bud = R.from_rotvec(self.bud_angle * vec_rotate)
+            # Note - newer versions use as_matrix
+            mat_rot = mat_rotate_bud.as_dcm()
+            vec[0:3] = mat_rot @ vec[0:3]
+            vec_on_crv = mat @ vec
+            vec_on_crv = vec_on_crv * (self.bud_length / np.linalg.norm(vec_on_crv))
+            pt_end_bud = pt_on_srf + vec_on_crv
+            pt_mid = 0.7 * pt_on_srf + 0.3 * pt_on_crv
+            ret_list.append((pt_mid[0:3], pt_on_srf[0:3], pt_end_bud[0:3]))
+
+        return ret_list
+
+
+if __name__ == '__main__':
+
+    branch = BezierCyl3DWithDetail()
+
+    branch.set_pts([506.5, 156.0, 0.0], [457.49999996771703, 478.9999900052037, 0.0], [521.5, 318.0, 0.0])
+    branch.set_radii_and_junction(start_radius=10.5, end_radius=8.25, b_start_is_junction=True, b_end_is_bud=False)
+    branch.write_mesh("data/jos.obj")
+
+    branch.set_pts([-0.5, 0.0, 0.0], [0.0, 0.1, 0.05], [0.5, 0.0, 0.0])
+    branch.set_radii_and_junction(start_radius=0.5, end_radius=0.25, b_start_is_junction=True, b_end_is_bud=False)
+    branch.write_mesh("data/cyl.obj")
+
+    branch.set_dims(n_along=30, n_radial=32)
+    branch.set_pts([-0.5, 0.0, 0.0], [0.0, 0.1, 0.05], [0.5, 0.0, 0.0])
+    branch.set_radii_and_junction(start_radius=0.1, end_radius=0.075, b_start_is_junction=False, b_end_is_bud=True)
+    branch.write_mesh("data/cyl_bud.obj")
+
+    bud_loc = branch.place_buds(((0.2, 0), (0.3, np.pi/4), (0.4, 3.0 * np.pi/4)))
+    bud = BezierCyl3DWithDetail()
+    bud.start_bud = 0.2
+    for i, b in enumerate(bud_loc):
+        bud.set_pts(b[0], b[1], b[2])
+        bud.make_branch_segment(b[0], b[1], b[2], radius_start=0.025, radius_end=0.03, start_is_junction=False, end_is_bud=True)
+        bud.write_mesh(f"data/bud_{i}.obj")