diff --git a/nuc_morph_analysis/analyses/density/extra_checks/visually_validate_watershed_psuedo_cell_seg_workflow.py b/nuc_morph_analysis/analyses/density/extra_checks/visually_validate_watershed_psuedo_cell_seg_workflow.py
index a87ea3d9..716b17c3 100644
--- a/nuc_morph_analysis/analyses/density/extra_checks/visually_validate_watershed_psuedo_cell_seg_workflow.py
+++ b/nuc_morph_analysis/analyses/density/extra_checks/visually_validate_watershed_psuedo_cell_seg_workflow.py
@@ -20,7 +20,6 @@ def get_contours_from_pair_of_2d_seg_image(nuc_mip,cell_mip,dft=None):
                 #ask if the label_img is in the dataframe
                 if label_img not in dft['label_img'].values:
                     continue
-            # color = np.float64(rgb_array0_255[label_img % len(rgb_array0_255)])
             color = np.float64(cmapper(label_img))
 
             # get the nucleus boundary
@@ -96,8 +95,6 @@ def plot_colorized_image_with_contours(img_dict,dft,feature,cmapstr,colony='test
         cmap = cm.get_cmap(cmapstr)
         if categorical:
             cimg = np.round(cimg).astype('uint16')
-
-        # rgb = np.take(np.uint16(cmaparr*255),cimg.astype('uint16'),axis=0)
     
         # create the figure
         fig, axlist = plt.subplots(1, 1, figsize=(6, 4))
@@ -245,6 +242,4 @@ def run_validation_and_plot(TIMEPOINT=48,colony='medium',RESOLUTION_LEVEL=1,plot
 if __name__ == '__main__':
     # set the details
     dft_test = run_validation_and_plot(testing=True,plot_everything=False)
-    dft0 = run_validation_and_plot(testing=False,plot_everything=False)
-    # dft0 = run_validation_and_plot(247,colony='small',RESOLUTION_LEVEL=1,plot_everything=True)
-    # dft0 = run_validation_and_plot(110,colony='medium',RESOLUTION_LEVEL=1,plot_everything=True)
\ No newline at end of file
+    dft0 = run_validation_and_plot(testing=False,plot_everything=False)
\ No newline at end of file
diff --git a/nuc_morph_analysis/analyses/density/figure_watershed_based_density_schematic.py b/nuc_morph_analysis/analyses/density/figure_watershed_based_density_schematic.py
index 55ac9e02..6419976c 100644
--- a/nuc_morph_analysis/analyses/density/figure_watershed_based_density_schematic.py
+++ b/nuc_morph_analysis/analyses/density/figure_watershed_based_density_schematic.py
@@ -17,6 +17,8 @@
 from nuc_morph_analysis.analyses.density.extra_checks.visually_validate_watershed_psuedo_cell_seg_workflow import get_contours_from_pair_of_2d_seg_image, draw_contours_on_image
 from nuc_morph_analysis.analyses.dataset_images_for_figures.figure_helper import INTENSITIES_DICT
 
+from nuc_morph_analysis.lib.visualization.example_tracks import EXAMPLE_TRACKS
+
 import matplotlib
 matplotlib.rcParams['pdf.fonttype'] = 42
 matplotlib.rcParams['ps.fonttype'] = 42
@@ -61,7 +63,7 @@ def determine_colormaps(img,key,crop_exp):
     return cmap, vmin, vmax
 
 
-def run_validation_and_plot(TIMEPOINT=88,track=81463,colony='medium',RESOLUTION_LEVEL=1,plot_everything=False, testing=False):
+def run_validation_and_plot(TIMEPOINT=88,track=EXAMPLE_TRACKS["pseudocell_density_example"],colony='medium',RESOLUTION_LEVEL=1,plot_everything=False, testing=False):
     """
     run an image through the watershed based pseudo cell segmentation and examine the outputs (as full fov and crop within that fov)
     optionally, run a test image through the same pipeline
diff --git a/nuc_morph_analysis/analyses/height/figure_3_s4_workflow.py b/nuc_morph_analysis/analyses/height/figure_3_s4_workflow.py
index 7f7243b7..fd37c399 100644
--- a/nuc_morph_analysis/analyses/height/figure_3_s4_workflow.py
+++ b/nuc_morph_analysis/analyses/height/figure_3_s4_workflow.py
@@ -1,15 +1,16 @@
 # %%
 from nuc_morph_analysis.analyses.colony_context.colony_context_analysis import plot_radial_profile
 from nuc_morph_analysis.lib.preprocessing.global_dataset_filtering import load_dataset_with_features
-from nuc_morph_analysis.lib.preprocessing import load_data, filter_data, global_dataset_filtering
+from nuc_morph_analysis.lib.preprocessing import load_data, filter_data
 from nuc_morph_analysis.lib.preprocessing.load_data import get_dataset_pixel_size
 from nuc_morph_analysis.analyses.height import plot
-from nuc_morph_analysis.analyses.height.plot_crowding import plot_density_schematic
 from nuc_morph_analysis.analyses.height.toymodel import toymodel
 from nuc_morph_analysis.analyses.height.centroid import (
     get_centroid,
     get_neighbor_centroids,
 )
+from nuc_morph_analysis.analyses.density import figure_watershed_based_density_schematic
+from nuc_morph_analysis.analyses.neighbor_of_X import figure_mitotic_filtering_examples
 from nuc_morph_analysis.lib.preprocessing.load_data import get_dataset_pixel_size
 from pathlib import Path
 import matplotlib.pyplot as plt
@@ -67,13 +68,12 @@
 
 # plot density schematic
 pix_size = get_dataset_pixel_size("medium")
-plot_density_schematic(
-    df_timepoint, track_centroid, neighbor_centroids, frame_centroids, pix_size, figdir,
-)
 
 # plot colony-averaged density over aligned colony time
-plot.density_colony_time_alignment(df_all, pixel_size, interval, time_axis="colony_time",
-                                   use_old_density=False)
+plot.density_colony_time_alignment(df_all, pixel_size, interval, time_axis="colony_time")
+
+figure_watershed_based_density_schematic.run_validation_and_plot() #SuppFigS4 panel C
+figure_mitotic_filtering_examples.run_validation_and_plot() #SuppFigS4 panel D, this code takes ~6 min to run
 
 # %%
 # Run and plot toy model
diff --git a/nuc_morph_analysis/analyses/height/plot.py b/nuc_morph_analysis/analyses/height/plot.py
index f5864e0e..cac15359 100644
--- a/nuc_morph_analysis/analyses/height/plot.py
+++ b/nuc_morph_analysis/analyses/height/plot.py
@@ -96,7 +96,7 @@ def height_colony_time_alignment(
     )
 
 
-def calculate_mean_density(df, scale, use_old_density=False):
+def calculate_mean_density(df, scale):
     """
     Calculate the mean height for a given index_sequence (i.e. timepoint) and the standard deviation of the mean.
 
@@ -106,8 +106,6 @@ def calculate_mean_density(df, scale, use_old_density=False):
         DataFrame containing the data.
     pixel_size : float
         Pixel size in microns.
-    use_old_density : bool
-        Whether to use the old density calculation method ('density') or the new method ('2d_area_nuc_cell_ratio')
 
     Returns
     -------
@@ -118,7 +116,7 @@ def calculate_mean_density(df, scale, use_old_density=False):
     """
     mean = []
     standard_dev = []
-    feature_col = "density" if use_old_density else "2d_area_nuc_cell_ratio"
+    feature_col = "2d_area_nuc_cell_ratio"
     for _, df_frame in df.groupby("index_sequence"):
         density = df_frame[feature_col].values * scale
         mean.append(np.nanmean(density))
@@ -134,7 +132,6 @@ def density_colony_time_alignment(
     show_legend=False,
     error="percentile",
     figdir="height/figures",
-    use_old_density=False,
 ):
     """
     Plot the mean nuclear height across the colony over time for each colony. This is done in real time and colony time.
@@ -159,9 +156,6 @@ def density_colony_time_alignment(
     error: str
         "std" or percentile
 
-    use_old_density: bool
-        Whether to use the old density calculation method ('density') or the new method ('2d_area_nuc_cell_ratio')
-
     Returns
     -------
     Plot of mean nuclear height across the colony over time for each colony.
@@ -169,7 +163,7 @@ def density_colony_time_alignment(
     plt.close()
     fig, ax = plt.subplots(1, 1, figsize=(5, 4))
 
-    feature_col = "density" if use_old_density else "2d_area_nuc_cell_ratio"
+    feature_col = "2d_area_nuc_cell_ratio"
     scale, label, units, _ = get_plot_labels_for_metric(feature_col)
 
     for colony, df_colony in df.groupby("colony"):
@@ -212,7 +206,6 @@ def density_colony_time_alignment(
             time, mean_density, linewidth=1.2, color=color, label=COLONY_LABELS[colony], zorder=20
         )
 
-    # ax.set_ylim(0.0005, 0.0065)
     ax.set_ylabel(f"Average Density \n Across Colony {units}")
     ax.set_xlabel(x_label)
     if show_legend is True:
diff --git a/nuc_morph_analysis/analyses/height/plot_crowding.py b/nuc_morph_analysis/analyses/height/plot_crowding.py
deleted file mode 100644
index e4b7e986..00000000
--- a/nuc_morph_analysis/analyses/height/plot_crowding.py
+++ /dev/null
@@ -1,139 +0,0 @@
-import numpy as np
-import matplotlib.pyplot as plt
-from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
-
-from nuc_morph_analysis.lib.visualization.movie_tools import get_colorized_img
-from nuc_morph_analysis.lib.visualization.notebook_tools import save_and_show_plot
-
-
-def get_axis_limits_from_points(points, padding=250):
-    """
-    Get the axis limits for a plot based on a set of points.
-
-    Parameters
-    ----------
-    points : numpy.ndarray
-        An array of shape (N, 2) representing the x and y coordinates of the points.
-    padding : int, optional
-        The amount of padding to add to the axis limits (default is 250).
-
-    Returns
-    -------
-    xlim : list
-        A list containing the minimum and maximum x-axis limits.
-    ylim : list
-        A list containing the minimum and maximum y-axis limits.
-    """
-    min_x = np.min(points[:, 0])
-    max_x = np.max(points[:, 0])
-    min_y = np.min(points[:, 1])
-    max_y = np.max(points[:, 1])
-    xlim = [min_x - padding, max_x + padding]
-    ylim = [min_y - padding, max_y + padding]
-    return xlim, ylim
-
-
-def plot_density_schematic(
-    df_timepoint, track_centroid, neighbor_centroids, frame_centroids, pix_size, figdir
-):
-    """
-    Plot the schematic showing how density is calculated
-
-    Parameters
-    ----------
-    df_timepoint : DataFrame
-        The timepoint data containing information about the nuclei.
-    track_centroid : tuple
-        The centroid coordinates of the tracked nucleus.
-    neighbor_centroids : array-like
-        The centroid coordinates of the neighboring nuclei.
-    frame_centroids : array-like
-        The centroid coordinates of all nuclei in the frame.
-    pix_size : float
-        The pixel size in micrometers.
-    figdir : str
-        The directory to save the generated figure.
-
-    Returns
-    -------
-    None
-        This function does not return anything. The figure is saved in the specified directory.
-    """
-    img = get_colorized_img(df_timepoint, "index_sequence", 1, filter_vals=False)
-    fig, axs = plt.subplots(1, 2, figsize=(8, 5), dpi=300)
-    for ct, ax in enumerate(axs):
-        ax.imshow(img, cmap="Greys_r", origin="upper")
-        ax.tick_params(
-            axis="both",
-            which="both",
-            bottom=False,
-            top=False,
-            left=False,
-            right=False,
-            labelbottom=False,
-            labelleft=False,
-        )
-
-        # draw arrows to each neighbor
-        for centroid in neighbor_centroids:
-            ax.arrow(
-                track_centroid[0],
-                track_centroid[1],
-                centroid[0] - track_centroid[0],
-                centroid[1] - track_centroid[1],
-                head_width=10,
-                head_length=10,
-                fc="blue",
-                ec="blue",
-            )
-
-        if ct == 0:  # plot entire colony with all nuclei
-            xlim, ylim = get_axis_limits_from_points(frame_centroids, padding=250)
-            ax.set_xlim(xlim)
-            ax.set_ylim(ylim[1], ylim[0])
-            sz = 10
-            scale = 50
-            size_vertical = 30
-            scale_bar_y_pos = 0.02
-        else:  # plot zoomed in view of nucleus and neighbors
-            avg_neighbor_distance = (
-                np.mean(np.linalg.norm(neighbor_centroids - track_centroid, axis=1)) * pix_size
-            )
-            density = 1 / avg_neighbor_distance**2
-            xlim, ylim = get_axis_limits_from_points(neighbor_centroids, padding=50)
-            ax.set_xlim(xlim)
-            ax.set_ylim(ylim[1], ylim[0])
-            sz = 50
-            scale = 10
-            size_vertical = 6
-            title_str = (
-                f"Mean distance to neighbors $(r)$: {avg_neighbor_distance:.2f}$\\mu m$"
-                + "\n"
-                + f"Density $(1 / r ^2)$: {density:.2g}/$\\mu m^2$"
-            )
-            scale_bar_y_pos = -0.1
-
-            ax.set_title(title_str)
-
-        # create scalebar
-        scalebar = AnchoredSizeBar(
-            ax.transData,
-            scale / pix_size,
-            f"{scale} $\\mu m$",
-            "lower right",
-            bbox_to_anchor=(0.95, scale_bar_y_pos),
-            bbox_transform=ax.transAxes,
-            color="black",
-            pad=0,
-            borderpad=0,
-            frameon=False,
-            size_vertical=size_vertical,
-            label_top=False,
-        )
-        ax.add_artist(scalebar)
-
-        # draw nucleus centroid
-        ax.scatter(track_centroid[0], track_centroid[1], color="red", s=sz)
-
-    plt.tight_layout(pad=0)
-    save_and_show_plot(f"{figdir}/density_schematic", figure=fig, bbox_inches="tight")
diff --git a/nuc_morph_analysis/analyses/height/toymodel.py b/nuc_morph_analysis/analyses/height/toymodel.py
index d203d574..24d60abd 100644
--- a/nuc_morph_analysis/analyses/height/toymodel.py
+++ b/nuc_morph_analysis/analyses/height/toymodel.py
@@ -19,7 +19,7 @@
 plt.rcParams["font.family"] = "Arial"
 
 
-def toymodel(nc_ratio=0.29, cvol=0.5e6, cell_H_mod=0, num_workers=1, old_unfiltered_method=False):
+def toymodel(nc_ratio=0.29, cvol=0.5e6, cell_H_mod=0, num_workers=1):
     """
     main function to run
 
@@ -29,12 +29,11 @@ def toymodel(nc_ratio=0.29, cvol=0.5e6, cell_H_mod=0, num_workers=1, old_unfilte
     cvol: target volume of nuclei to fit toy model to, default is 0.5e6 pixels^3 ~ 630 um^3
     cell_H_mod: height difference between top of nucleus and top of cell to use for toy model. default is 0
     num_workers: how many workers to use for multiprocessing
-    old_unfiltered_method: whether to use the old (unfiltered) method that does not remove bad pseudo cells
     """
     save_path = Path(__file__).parent / "figures"
     save_path.mkdir(parents=True, exist_ok=True)
     pix_size = load_data.get_dataset_pixel_size("all_baseline")
-    data, cvol = get_data(cvol, save_path, num_workers, old_unfiltered_method)
+    data, cvol = get_data(cvol, save_path, num_workers)
     cvol_um = cvol * pix_size**3
     cvol_um_cell = cvol_um * 1 / nc_ratio
     # rescale features
@@ -43,8 +42,8 @@ def toymodel(nc_ratio=0.29, cvol=0.5e6, cell_H_mod=0, num_workers=1, old_unfilte
 
     stats = get_toy_model(data, cvol_um_cell, cell_H_mod)
 
-    plot_toy_model(data, stats, save_path, old_unfiltered_method)
-    plot_growth_rate(data, save_path, old_unfiltered_method)
+    plot_toy_model(data, stats, save_path)
+    plot_growth_rate(data, save_path)
 
 
 def get_toy_model(data, cvol_um_cell, cell_H_mod=0):
@@ -63,7 +62,6 @@ def get_toy_model(data, cvol_um_cell, cell_H_mod=0):
     cell_H_mod: height difference between top of nucleus and top of cell to use for toy model. default is 0
     """
     H = np.linspace(data["height"].min(), data["height"].max(), 100)
-    # cell_H = H + 2
     cell_H = H + cell_H_mod
     # Volume of cylinder is pi r^2 H
     # distance = np.sqrt(cvol/(np.pi * H))*2 for cylinder
@@ -89,7 +87,7 @@ def get_toy_model(data, cvol_um_cell, cell_H_mod=0):
     return stats
 
 
-def plot_toy_model(data, toy_stats, save_path=Path("./"), old_unfiltered_method=False):
+def plot_toy_model(data, toy_stats, save_path=Path("./")):
     """
     Plot real data and toy model curves
     For the real data, plot a gaussian weighted moving average
@@ -99,7 +97,6 @@ def plot_toy_model(data, toy_stats, save_path=Path("./"), old_unfiltered_method=
     data: dataframe with real data
     toy_stats: dataframe with toy model fits
     save_path: path to save pdf
-    old_unfiltered_method: whether to use the old (unfiltered) method that does not remove bad pseudo cells
     """
     data = data.sort_values(by="height")
     x = data["height"]
@@ -130,16 +127,14 @@ def plot_toy_model(data, toy_stats, save_path=Path("./"), old_unfiltered_method=
     )
     axes.set_ylim(10, 40)
     axes.legend()
-    suffix = "unfiltered.pdf" if old_unfiltered_method else ".pdf"
-    savename = f"toymodel{suffix}"
+    savename = f"toymodel.pdf"
     fig.savefig(save_path / savename, bbox_inches="tight")
 
 
-def plot_growth_rate(data, save_path=Path("./"), old_unfiltered_method=False):
+def plot_growth_rate(data, save_path=Path("./")):
     """
     Plot growth rate vs crowding
     crowding is calculated as a mean distance to neighbors
-    old_unfiltered_method: whether to use the old (unfiltered) method that does not remove bad pseudo cells
     """
     data = data.sort_values(by="growth_rate")
     y = data["growth_rate"]
@@ -157,8 +152,7 @@ def plot_growth_rate(data, save_path=Path("./"), old_unfiltered_method=False):
     axes.plot(bins, average, label=f"gaussian weighted moving average", c="tab:blue")
     axes.fill_between(bins, lower_bound, upper_bound, alpha=0.3, edgecolor="none")
     axes.legend()
-    savename = "growthrate_vs_crowding_unfiltered.pdf" if old_unfiltered_method else "growthrate_vs_crowding.pdf"
-    fig.savefig(save_path / savename, bbox_inches="tight")
+    fig.savefig(save_path / "growthrate_vs_crowding.pdf", bbox_inches="tight")
 
 
 def weighted_moving_average(x, y, step_size=0.05, width=1):
@@ -183,7 +177,7 @@ def gaussian(x, amp=1, mean=0, sigma=1):
     return (bin_centers, bin_avg, bin_std)
 
 
-def get_data(cvol, save_path=Path("./"), num_workers=1, old_unfiltered_method=False):
+def get_data(cvol, save_path=Path("./"), num_workers=1):
     """
     Load all datasets, remove outliers and edge cells, fitler to tracks > 120
     frames, select all nuclei that are close to a set volume (cvol), and compute
@@ -198,7 +192,6 @@ def get_data(cvol, save_path=Path("./"), num_workers=1, old_unfiltered_method=Fa
     
     df = global_dataset_filtering.load_dataset_with_features(remove_growth_outliers=False)
     df = filter_data.filter_all_outliers(df)
-    # df = filter_data.filter_out_cells_entering_or_exiting_mitosis(df)
 
     df_full = filter_data.all_timepoints_full_tracks(df)
     df_ft = df_full[df_full["colony"].isin(["small", "medium", "large"])].reset_index()
@@ -226,7 +219,7 @@ def get_data(cvol, save_path=Path("./"), num_workers=1, old_unfiltered_method=Fa
     df_all = df.loc[df.index.isin(neighbor_ids)].reset_index()
 
     # compute distance/density metric
-    neigh_stats = compute_density(df_cvol, df_all, num_workers, old_unfiltered_method)
+    neigh_stats = compute_density(df_cvol, df_all, num_workers)
     neigh_stats = neigh_stats.reset_index()
 
     return neigh_stats, cvol
@@ -247,6 +240,4 @@ def determine_volume_at_middle_of_cell_cycle():
 
 
 if __name__ == "__main__":
-    for cell_H_mod in [0, 1.5, 2]:
-        toymodel(cell_H_mod=cell_H_mod)
-        toymodel(cell_H_mod=cell_H_mod,old_unfiltered_method=True)
+    toymodel()
diff --git a/nuc_morph_analysis/analyses/neighbor_of_X/figure_mitotic_filtering_examples.py b/nuc_morph_analysis/analyses/neighbor_of_X/figure_mitotic_filtering_examples.py
index 52aee2f5..67080620 100644
--- a/nuc_morph_analysis/analyses/neighbor_of_X/figure_mitotic_filtering_examples.py
+++ b/nuc_morph_analysis/analyses/neighbor_of_X/figure_mitotic_filtering_examples.py
@@ -1,22 +1,20 @@
 #%%
 #SuppFigS4 panel D, this code takes ~6 min to run
-from nuc_morph_analysis.lib.preprocessing.twoD_zMIP_area import watershed_workflow, pseudo_cell_helper, pseudo_cell_testing_helper
+from nuc_morph_analysis.lib.preprocessing.twoD_zMIP_area import watershed_workflow
 from pathlib import Path
-import pandas as pd
 from nuc_morph_analysis.lib.preprocessing import global_dataset_filtering
-from nuc_morph_analysis.lib.preprocessing import filter_data, load_data
+from nuc_morph_analysis.lib.preprocessing import load_data
 from matplotlib.colors import ListedColormap
 
 import numpy as np
 import matplotlib.pyplot as plt
 import matplotlib.cm as cm
-from nuc_morph_analysis.lib.visualization.plotting_tools import colorize_image, get_plot_labels_for_metric
 from nuc_morph_analysis.lib.visualization.notebook_tools import save_and_show_plot
 
 from nuc_morph_analysis.analyses.dataset_images_for_figures.figure_helper import return_glasbey_on_dark
+from nuc_morph_analysis.lib.visualization.example_tracks import EXAMPLE_TRACKS
 
-
-from nuc_morph_analysis.analyses.density.visually_validate_watershed_psuedo_cell_seg_workflow import get_contours_from_pair_of_2d_seg_image, draw_contours_on_image
+from nuc_morph_analysis.analyses.density.extra_checks.visually_validate_watershed_psuedo_cell_seg_workflow import get_contours_from_pair_of_2d_seg_image, draw_contours_on_image
 from nuc_morph_analysis.analyses.dataset_images_for_figures.figure_helper import INTENSITIES_DICT
 from tqdm import tqdm
 
@@ -64,7 +62,7 @@ def determine_colormaps(img,key,crop_exp):
         vmax= INTENSITIES_DICT['egfp_max'][1]
     return cmap, vmin, vmax
 
-def run_validation_and_plot(track_id=87135,RESOLUTION_LEVEL=1,frames_before=0,frames_after=7,w=300,plot_everything=False):
+def run_validation_and_plot(track_id=EXAMPLE_TRACKS['pseudocell_mitoticfilter_example'],RESOLUTION_LEVEL=1,frames_before=0,frames_after=7,w=300,plot_everything=False):
     """
     run an image through the watershed based pseudo cell segmentation and examine the outputs
     optionally, run a test image through the same pipeline
@@ -92,7 +90,6 @@ def run_validation_and_plot(track_id=87135,RESOLUTION_LEVEL=1,frames_before=0,fr
 
    # load the tracking dataframe and apply appropriate filters
     df = global_dataset_filtering.load_dataset_with_features(dataset='all_baseline')
-    # df = filter_data.all_timepoints_minimal_filtering(df)
     
     dftrack = df.loc[df['track_id']==track_id].copy()
     timepoint = int(dftrack['predicted_breakdown'].values[0])
@@ -102,10 +99,6 @@ def run_validation_and_plot(track_id=87135,RESOLUTION_LEVEL=1,frames_before=0,fr
     dfc = df.loc[df['colony']==colony].copy()
     dfm = dfc.loc[(dfc['index_sequence'].isin(time_list))].copy()
 
-    # set figure directory
-    figdir = Path(__file__).parent / "figures" / "SuppFigS4_mitotic_filtering"
-    figdir.mkdir(exist_ok=True,parents=True)
-
     # artificially set all nuclei to have predicted_breakdown and predicted_formation to -1
     dft = dfm[dfm['index_sequence']==timepoint]
     x1,y1,w,h = determine_crop_size_and_location_from_track(track_id,dft,crop_w=w,crop_h=w,RESOLUTION_LEVEL=1)
@@ -178,9 +171,6 @@ def run_validation_and_plot(track_id=87135,RESOLUTION_LEVEL=1,frames_before=0,fr
                 colormap_dict.update({'frame_of_breakdown':('frame_of_breakdown',True,8,(1.0,1.0,0.0),f"breakdown event")})
                 colormap_dict.update({'frame_of_formation':('frame_of_formation',True,9,(0.0,1.0,1.0),f"formation event")})
 
-                # now update colors in contour_list based on colormap_dict
-                # contour_list.append((label_img,nuc_contours,cell_contours,color))
-
                 new_colors = np.zeros((np.max([x[2] for x in colormap_dict.values()])+1,3))
                 for col in colormap_dict.keys():
                     new_colors[colormap_dict[col][2]] = colormap_dict[col][3]
diff --git a/nuc_morph_analysis/lib/preprocessing/neighbor_analysis/neighbor_analysis.py b/nuc_morph_analysis/lib/preprocessing/neighbor_analysis/neighbor_analysis.py
index 867fff5b..07aaedd5 100644
--- a/nuc_morph_analysis/lib/preprocessing/neighbor_analysis/neighbor_analysis.py
+++ b/nuc_morph_analysis/lib/preprocessing/neighbor_analysis/neighbor_analysis.py
@@ -88,7 +88,7 @@ def compute_neigh(tup):
             return return_dict
 
 
-def compute_density(df, global_df, num_workers=1, old_unfiltered_method=False):
+def compute_density(df, global_df, num_workers=1):
     """
     Main function to run
     Given a dataframe, compute density for each row
@@ -97,8 +97,6 @@ def compute_density(df, global_df, num_workers=1, old_unfiltered_method=False):
     df: dataframe of cell ids with location and height information
     global_df: dataframe of all neighboring cell ids with location and height information
     num_workers: number of workers for multiprocessing
-    old_unfiltered_method: whether to use the old (unfiltered) method that does not remove bad pseudo cells
-    (e.g. cells next to mitotic cells where there are missing segmentations)
     """
     feature_keys = [
         "centroid_x",
@@ -130,15 +128,12 @@ def compute_density(df, global_df, num_workers=1, old_unfiltered_method=False):
     neigh_stats = pd.concat(neigh_stats, axis=0)
     neigh_stats = neigh_stats.groupby(["CellId"]).mean()
 
-
-    print('old_unfiltered_method', old_unfiltered_method)
-    if not old_unfiltered_method:
-        # now remove CellIds that have bad pseudo cell segmentation (i.e. bad_pseudo_cells_segmentation)
-        # 'uncaught_pseudo_cell_artifact','bad_pseudo_cells_segmentation'
-        print("Removing bad pseudo cells")
-        print('Before removing bad pseudo cells: ', neigh_stats.shape[0])
-        cell_ids_to_remove = df[df['bad_pseudo_cells_segmentation'] == True]['CellId'].values
-        neigh_stats = neigh_stats[~neigh_stats.index.isin(cell_ids_to_remove)]
-        print('After removing bad pseudo cells: ', neigh_stats.shape[0])
+    # now remove CellIds that have bad pseudo cell segmentation (i.e. bad_pseudo_cells_segmentation)
+    # 'uncaught_pseudo_cell_artifact','bad_pseudo_cells_segmentation'
+    print("Removing bad pseudo cells")
+    print('Before removing bad pseudo cells: ', neigh_stats.shape[0])
+    cell_ids_to_remove = df[df['bad_pseudo_cells_segmentation'] == True]['CellId'].values
+    neigh_stats = neigh_stats[~neigh_stats.index.isin(cell_ids_to_remove)]
+    print('After removing bad pseudo cells: ', neigh_stats.shape[0])
 
     return neigh_stats
diff --git a/nuc_morph_analysis/lib/visualization/example_tracks.py b/nuc_morph_analysis/lib/visualization/example_tracks.py
index 32ea6b82..52366a00 100644
--- a/nuc_morph_analysis/lib/visualization/example_tracks.py
+++ b/nuc_morph_analysis/lib/visualization/example_tracks.py
@@ -23,4 +23,6 @@
     "delta_v_BC_low": 86418,
     "transition_point_supplement": 82210,
     "sample_full_trajectories": [97942, 85296, 9808, 77656, 83322],
+    "pseudocell_density_example": 81463,
+    "pseudocell_mitoticfilter_example": 87135
 }
diff --git a/run_all_manuscript_workflows.py b/run_all_manuscript_workflows.py
index 45701047..f6217cc4 100644
--- a/run_all_manuscript_workflows.py
+++ b/run_all_manuscript_workflows.py
@@ -22,11 +22,6 @@ def figure_s2_s3_s18_segmentation_model_validation():
     def figure_3_s4_height_density():
         # figure 3 images generated using timelapse feauture explorer
         import nuc_morph_analysis.analyses.height.figure_3_s4_workflow
-        from nuc_morph_analysis.analyses.density import figure_watershed_based_density_schematic
-        from nuc_morph_analysis.analyses.neighbor_of_X import figure_mitotic_filtering_examples
-
-        figure_watershed_based_density_schematic.run_validation_and_plot() #SuppFigS4 panel C
-        figure_mitotic_filtering_examples.run_validation_and_plot() #SuppFigS4 panel D, this code takes ~6 min to run
         
     # figure s5 image data can be found on quilt
     
diff --git a/test/lib/preprocessing/test_voronoi.py b/test/lib/preprocessing/test_voronoi.py
index 23d19745..fe95855b 100644
--- a/test/lib/preprocessing/test_voronoi.py
+++ b/test/lib/preprocessing/test_voronoi.py
@@ -117,7 +117,6 @@ def test_voronoi_real_densities():
     df_colony_metrics["old_density"] = df_colony_metrics["density"].apply(lambda x: 4 * np.sqrt(x))
     assert np.allclose(df_colony_metrics["old_density"], df_colony_metrics.density_gt, rtol=0.2)
 
-
 def test_voronoi_synthetic_distance_density():
     """
     This test uses a set of cells laid out in the following pattern.