Merge pull request flatironinstitute#610 from phildong/hv_vis

add visualizations based on holoviews

caiman/source_extraction/cnmf/estimates.py

@@ -381,6 +381,55 @@ def nb_view_components(self, Yr=None, img=None, idx=None,
                                                         thr=thr, denoised_color=denoised_color, cmap=cmap)
         return self
 
+    def hv_view_components(self, Yr=None, img=None, idx=None,
+                           denoised_color=None, cmap='viridis'):
+        """view spatial and temporal components interactively in a notebook
+
+        Args:
+            Yr :    np.ndarray
+                movie in format pixels (d) x frames (T)
+
+            img :   np.ndarray
+                background image for contour plotting. Default is the mean
+                image of all spatial components (d1 x d2)
+
+            idx :   list
+                list of components to be plotted
+
+            denoised_color: string or None
+                color name (e.g. 'red') or hex color code (e.g. '#F0027F')
+
+            cmap: string
+                name of colormap (e.g. 'viridis') used to plot image_neurons
+        """
+        if 'csc_matrix' not in str(type(self.A)):
+            self.A = scipy.sparse.csc_matrix(self.A)
+
+        plt.ion()
+        nr, T = self.C.shape
+        if self.R is None:
+            self.R = self.YrA
+        if self.R.shape != [nr, T]:
+            if self.YrA is None:
+                self.compute_residuals(Yr)
+            else:
+                self.R = self.YrA
+
+        if img is None:
+            img = np.reshape(np.array(self.A.mean(axis=1)), self.dims, order='F')
+
+        if idx is None:
+            hv_plot = caiman.utils.visualization.hv_view_patches(
+                Yr, self.A, self.C, self.b, self.f, self.dims[0], self.dims[1],
+                YrA=self.R, image_neurons=img, denoised_color=denoised_color,
+                cmap=cmap)
+        else:
+            hv_plot = caiman.utils.visualization.hv_view_patches(
+                Yr, self.A.tocsc()[:, idx], self.C[idx], self.b, self.f,
+                self.dims[0], self.dims[1], YrA=self.R[idx], image_neurons=img,
+                denoised_color=denoised_color, cmap=cmap)
+        return hv_plot
+
     def nb_view_components_3d(self, Yr=None, image_type='mean', dims=None,
                               max_projection=False, axis=0,
                               denoised_color=None, cmap='jet', thr=0.9):

caiman/utils/visualization.py

@@ -30,6 +30,8 @@
 from tempfile import NamedTemporaryFile
 from typing import Dict
 from warnings import warn
+import holoviews as hv
+import functools as fct
 
 from ..base.rois import com
 from ..summary_images import local_correlations
@@ -248,6 +250,76 @@ def nb_view_patches(Yr, A, C, b, f, d1, d2, YrA=None, image_neurons=None, thr=0.
     return Y_r
 
 
+def hv_view_patches(Yr, A, C, b, f, d1, d2, YrA=None, image_neurons=None, denoised_color=None, cmap='viridis'):
+    """
+    Interactive plotting utility for ipython notebook
+
+    Args:
+        Yr: np.ndarray
+            movie
+
+        A,C,b,f: np.ndarrays
+            outputs of matrix factorization algorithm
+
+        d1,d2: floats
+            dimensions of movie (x and y)
+
+        YrA:   np.ndarray
+            ROI filtered residual as it is given from update_temporal_components
+            If not given, then it is computed (K x T)
+
+        image_neurons: np.ndarray
+            image to be overlaid to neurons (for instance the average)
+
+        denoised_color: string or None
+            color name (e.g. 'red') or hex color code (e.g. '#F0027F')
+
+        cmap: string
+            name of colormap (e.g. 'viridis') used to plot image_neurons
+    """
+
+    nr, T = C.shape
+    nA2 = np.ravel(np.power(A, 2).sum(0)) if type(
+        A) == np.ndarray else np.ravel(A.power(2).sum(0))
+    b = np.squeeze(b)
+    f = np.squeeze(f)
+    if YrA is None:
+        Y_r = np.array(
+            spdiags(old_div(1, nA2), 0, nr, nr) *
+            (A.T * np.matrix(Yr) -
+             (A.T * np.matrix(b[:, np.newaxis])) * np.matrix(f[np.newaxis]) -
+             A.T.dot(A) * np.matrix(C)) + C)
+    else:
+        Y_r = C + YrA
+    if image_neurons is None:
+        image_neurons = A.mean(1).reshape((d1, d2), order='F')
+    smp = cm.ScalarMappable(cmap=cmap)
+    im_rgb = smp.to_rgba(image_neurons)[:, :, :3]
+    im_hsv = mpl.colors.rgb_to_hsv(im_rgb)
+
+    def norm(a, rg=(0, 1)):
+        a_norm = (a - a.min()) / (a.max() - a.min())
+        return a_norm * (rg[1] - rg[0]) + rg[0]
+
+    Ad = np.asarray(A.todense()).reshape((d1, d2, -1), order='F')
+
+    def plot_unit(uid, scl):
+        trace = (
+            hv.Curve(Y_r[uid, :], kdims='time').opts(framewise=True)
+            * (hv.Curve(C[uid, :], kdims='time')
+               .opts(color=denoised_color, framewise=True))
+        ).opts(aspect=2, frame_height=200)
+        A_scl = norm(Ad[:, :, uid], (scl, 1))
+        im_hsv_scl = im_hsv.copy()
+        im_hsv_scl[:, :, 2] = im_hsv[:, :, 2] * A_scl
+        im_u = (hv.HSV(im_hsv_scl, kdims=['height', 'width'])
+                .opts(aspect='equal', frame_height=200))
+        return im_u + trace
+
+    return (hv.DynamicMap(plot_unit, kdims=['unit_id', 'scale'])
+            .redim.range(unit_id=(0, nr), scale=(0.0, 1.0)))
+
+
 def get_contours(A, dims, thr=0.9, thr_method='nrg', swap_dim=False):
     """Gets contour of spatial components and returns their coordinates
 
@@ -992,6 +1064,36 @@ def GetBox(centers, R, dims):
         ax.axis('off')
     pl.subplots_adjust(0, 0, 1, 1, .06, .06)
 
+
+def nb_inspect_correlation_pnr(corr, pnr):
+    """
+    inspect correlation and pnr images to infer the min_corr, min_pnr
+
+    Args:
+        corr: ndarray
+            correlation image created with caiman.summary_images.correlation_pnr
+
+        pnr: ndarray
+            peak-to-noise image created with caiman.summary_images.correlation_pnr
+    """
+    hv_corr = hv.Image(corr, vdims='corr', label='correlation')
+    hv_pnr = hv.Image(pnr, vdims='pnr', label='pnr')
+
+    def hist(im, rx, ry):
+        obj = im.select(x=rx, y=ry) if rx and ry else im
+        return hv.operation.histogram(obj)
+
+    str_corr = (hv.streams.RangeXY(source=hv_corr)
+                .rename(x_range='rx', y_range='ry'))
+    str_pnr = (hv.streams.RangeXY(source=hv_pnr)
+               .rename(x_range='rx', y_range='ry'))
+    hist_corr = hv.DynamicMap(
+        fct.partial(hist, im=hv_corr), streams=[str_corr])
+    hist_pnr = hv.DynamicMap(
+        fct.partial(hist, im=hv_pnr), streams=[str_pnr])
+    return (hv_corr << hist_corr) + (hv_pnr << hist_pnr)
+
+
 def inspect_correlation_pnr(correlation_image_pnr, pnr_image):
     """
     inspect correlation and pnr images to infer the min_corr, min_pnr

demos/notebooks/demo_pipeline_cnmfE.ipynb

@@ -43,7 +43,7 @@
     "import caiman as cm\n",
     "from caiman.source_extraction import cnmf\n",
     "from caiman.utils.utils import download_demo\n",
-    "from caiman.utils.visualization import inspect_correlation_pnr\n",
+    "from caiman.utils.visualization import inspect_correlation_pnr, nb_inspect_correlation_pnr\n",
     "from caiman.motion_correction import MotionCorrect\n",
     "from caiman.source_extraction.cnmf import params as params\n",
     "from caiman.utils.visualization import plot_contours, nb_view_patches, nb_plot_contour\n",
@@ -54,7 +54,9 @@
     "except:\n",
     "    pass\n",
     "import bokeh.plotting as bpl\n",
-    "bpl.output_notebook()"
+    "import holoviews as hv\n",
+    "bpl.output_notebook()\n",
+    "hv.notebook_extension('bokeh')"
    ]
   },
   {
@@ -306,7 +308,7 @@
     "# compute some summary images (correlation and peak to noise)\n",
     "cn_filter, pnr = cm.summary_images.correlation_pnr(images[::1], gSig=gSig[0], swap_dim=False) # change swap dim if output looks weird, it is a problem with tiffile\n",
     "# inspect the summary images and set the parameters\n",
-    "inspect_correlation_pnr(cn_filter, pnr)"
+    "nb_inspect_correlation_pnr(cn_filter, pnr)"
    ]
   },
   {
@@ -440,7 +442,7 @@
    "outputs": [],
    "source": [
     "# accepted components\n",
-    "cnm.estimates.nb_view_components(img=cn_filter, idx=cnm.estimates.idx_components,\n",
+    "cnm.estimates.hv_view_components(img=cn_filter, idx=cnm.estimates.idx_components,\n",
     "                                denoised_color='red', cmap='gray')"
    ]
   },
@@ -453,7 +455,7 @@
    "outputs": [],
    "source": [
     "# rejected components\n",
-    "cnm.estimates.nb_view_components(img=cn_filter, idx=cnm.estimates.idx_components_bad,\n",
+    "cnm.estimates.hv_view_components(img=cn_filter, idx=cnm.estimates.idx_components_bad,\n",
     "                                denoised_color='red', cmap='gray')"
    ]
   },