{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# **Point Couds with zarr**\n\nGenerate Point Clouds while leveraging efficient image loading with zarr. Since we will be able to use all surfaces let's see if we can denoise the point cloud. \n\nCredit: https://www.kaggle.com/code/brettolsen/efficient-image-loading-with-zarr/notebookimport os","metadata":{}},{"cell_type":"code","source":"import os\nimport shutil\nfrom tifffile import tifffile\nimport time\nimport numpy as np\nimport PIL.Image as Image\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as patches\nfrom tqdm import tqdm\nfrom ipywidgets import interact, fixed\nfrom IPython.display import HTML, display","metadata":{"execution":{"iopub.status.busy":"2023-04-15T18:41:06.644256Z","iopub.execute_input":"2023-04-15T18:41:06.644744Z","iopub.status.idle":"2023-04-15T18:41:06.816788Z","shell.execute_reply.started":"2023-04-15T18:41:06.644705Z","shell.execute_reply":"2023-04-15T18:41:06.815547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install open3d\n!pip install zarr","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import zarr\nimport open3d as o3","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:47:58.706567Z","iopub.execute_input":"2023-04-15T17:47:58.707237Z","iopub.status.idle":"2023-04-15T17:48:04.031051Z","shell.execute_reply.started":"2023-04-15T17:47:58.707171Z","shell.execute_reply":"2023-04-15T17:48:04.029957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"INPUT_FOLDER = \"/kaggle/input/vesuvius-challenge-ink-detection\"\nWORKING_FOLDER = \"/kaggle/working/\"\nTEMP_FOLDER = \"kaggle/temp/\"","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:48:08.542915Z","iopub.execute_input":"2023-04-15T17:48:08.543349Z","iopub.status.idle":"2023-04-15T17:48:08.550138Z","shell.execute_reply.started":"2023-04-15T17:48:08.543313Z","shell.execute_reply":"2023-04-15T17:48:08.548788Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class TimerError(Exception):\n    pass\n\nclass Timer():\n    def __init__(self, text=None):\n        if text is not None:\n            self.text = text + \": {:0.4f} seconds\"\n        else:\n            self.text = \"Elapsed time: {:0.4f} seconds\"\n        def logfunc(x):\n            print(x)\n        self.logger = logfunc\n        self._start_time = None\n\n    def start(self):\n        if self._start_time is not None:\n            raise TimerError(\"Timer is already running.  Use .stop() to stop it.\")\n        self._start_time = time.time()\n\n    def stop(self):\n        if self._start_time is None:\n            raise TimerError(\"Timer is not running.  Use .start() to start it.\")\n        elapsed_time = time.time() - self._start_time\n        self._start_time = None\n\n        if self.logger is not None:\n            self.logger(self.text.format(elapsed_time))\n\n        return elapsed_time\n\n    def __enter__(self):\n        self.start()\n        return self\n\n    def __exit__(self, exc_type, exc_value, exc_traceback):\n        self.stop()","metadata":{"execution":{"iopub.status.busy":"2023-04-15T18:41:12.861046Z","iopub.execute_input":"2023-04-15T18:41:12.861568Z","iopub.status.idle":"2023-04-15T18:41:12.873551Z","shell.execute_reply.started":"2023-04-15T18:41:12.861521Z","shell.execute_reply":"2023-04-15T18:41:12.872175Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class FragmentImageException(Exception):\n    pass\n\nclass FragmentImageData:\n    \"\"\"A general class that uses persistent zarr objects to store the surface volume data,\n    binary data mask, and for training sets, the truth data and infrared image of a papyrus\n    fragment, in a compressed and efficient way.\n    \"\"\"\n    def __init__(self, sample_type: str, sample_index: str, working: bool = True):\n        if sample_type not in (\"test, train\"):\n            raise FragmentImageException(\n                f\"Invalid sample type f{sample_type}, must be one of 'test' or 'train'\"\n            )\n        zarrpath = self._zarr_path(sample_type, sample_index, working)\n        if os.path.exists(zarrpath):\n            self.zarr = self.load_from_zarr(zarrpath)\n        else:\n            dirpath = os.path.join(INPUT_FOLDER, sample_type, sample_index)\n            if not os.path.exists(dirpath):\n                raise FragmentImageException(\n                    f\"No input data found at f{zarrpath} or f{dirpath}\"\n                )\n            self.zarr = self.load_from_directory(dirpath, zarrpath)\n    \n    @property\n    def surface_volume(self):\n        return self.zarr.surface_volume\n    \n    @property\n    def mask(self):\n        return self.zarr.mask\n    \n    @property\n    def truth(self):\n        return self.zarr.truth\n    \n    @property\n    def infrared(self):\n        return self.zarr.infrared\n    \n    @staticmethod\n    def _zarr_path(sample_type: str, sample_index: str, working: bool = True):\n        filename = f\"{sample_type}-{sample_index}.zarr\"\n        if working:\n            return os.path.join(WORKING_FOLDER, filename)\n        else:\n            return os.path.join(TEMP_FOLDER, filename)\n    \n    @staticmethod\n    def clean_zarr(sample_type: str, sample_index: str, working: bool = True):\n        zarrpath = FragmentImageData._zarr_path(sample_type, sample_index, working)\n        if os.path.exists(zarrpath):\n            shutil.rmtree(zarrpath)\n    \n    @staticmethod\n    def load_from_zarr(filepath):\n        with Timer(\"Loading from existing zarr\"):\n            return zarr.open(filepath, mode=\"r\")\n    \n    @staticmethod\n    def load_from_directory(dirpath, zarrpath):\n        if os.path.exists(zarrpath):\n            raise FragmentImageException(\n                f\"Trying to overwrite existing zarr at f{zarrpath}\"\n            )\n        # Initialize the root zarr group and write the file\n        root = zarr.open_group(zarrpath, mode=\"w\")\n        # Load in the surface volume tif files\n        with Timer(\"Surface volume loading\"):\n            init = True\n            imgfiles = sorted([\n                imgfile for imgfile in\n                os.listdir(os.path.join(dirpath, \"surface_volume\"))\n            ])\n            for imgfile in imgfiles:\n                print(f\"Loading file {imgfile}\", end=\"\\r\")\n#                 img_data = np.array(\n#                     Image.open(os.path.join(dirpath, \"surface_volume\", imgfile))\n#                 )\n                img_data = tifffile.imread(\n                    os.path.join(dirpath, \"surface_volume\", imgfile)\n                )\n                if init:\n                    surface_volume = root.zeros(\n                        name=\"surface_volume\",\n                        shape=(img_data.shape[0], img_data.shape[1], len(imgfiles)),\n                        chunks=(1000, 1000, 4),\n                        dtype=img_data.dtype,\n                        write_empty_chunks=False,\n                    )\n                    init = False\n                z_index = int(imgfile.split(\".\")[0])\n                surface_volume[:,:,z_index] = img_data\n        # Load in the mask\n        with Timer(\"Mask loading\"):\n            img_data = np.array(Image.open(os.path.join(dirpath, \"mask.png\")), dtype=bool)\n            mask = root.array(\n                name=\"mask\",\n                data=img_data,\n                shape=img_data.shape,\n                chunks=(1000, 1000),\n                dtype=img_data.dtype,\n                write_empty_chunks=False,\n            )\n        # Load in the truth set (if it exists)\n        with Timer(\"Truth set loading\"):\n            truthfile = os.path.join(dirpath, \"inklabels.png\")\n            if os.path.exists(truthfile):\n                img_data = np.array(Image.open(truthfile), dtype=bool)\n                truth = root.array(\n                    name=\"truth\",\n                    data=img_data,\n                    shape=img_data.shape,\n                    chunks=(1000, 1000),\n                    dtype=img_data.dtype,\n                    write_empty_chunks=False,\n                )\n        # Load in the infrared image (if it exists)\n        with Timer(\"Infrared image loading\"):\n            irfile = os.path.join(dirpath, \"ir.png\")\n            if os.path.exists(irfile):\n                img_data = np.array(Image.open(irfile))\n                infrared = root.array(\n                    name = \"infrared\",\n                    data = img_data,\n                    shape = img_data.shape,\n                    chunks = (1000, 1000),\n                    dtype=img_data.dtype,\n                    write_empty_chunks=False,\n                )\n        return root        ","metadata":{"execution":{"iopub.status.busy":"2023-04-15T18:41:28.323507Z","iopub.execute_input":"2023-04-15T18:41:28.325342Z","iopub.status.idle":"2023-04-15T18:41:28.353608Z","shell.execute_reply.started":"2023-04-15T18:41:28.325266Z","shell.execute_reply":"2023-04-15T18:41:28.352097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Load data","metadata":{}},{"cell_type":"code","source":"FragmentImageData.clean_zarr(\"train\", 1)\ndata = FragmentImageData(\"train\", '1')","metadata":{"execution":{"iopub.status.busy":"2023-04-15T18:41:35.317780Z","iopub.execute_input":"2023-04-15T18:41:35.318449Z","iopub.status.idle":"2023-04-15T18:44:15.005081Z","shell.execute_reply.started":"2023-04-15T18:41:35.318391Z","shell.execute_reply":"2023-04-15T18:44:15.003797Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(data.surface_volume.info)\nprint(data.mask.info)\nprint(data.truth.info)\nprint(data.infrared.info)","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:51:26.676824Z","iopub.execute_input":"2023-04-15T17:51:26.677294Z","iopub.status.idle":"2023-04-15T17:51:26.736602Z","shell.execute_reply.started":"2023-04-15T17:51:26.677243Z","shell.execute_reply":"2023-04-15T17:51:26.732850Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with Timer():\n    plt.imshow(data.mask, cmap=\"gray\")","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:51:26.739327Z","iopub.execute_input":"2023-04-15T17:51:26.739679Z","iopub.status.idle":"2023-04-15T17:51:29.423188Z","shell.execute_reply.started":"2023-04-15T17:51:26.739646Z","shell.execute_reply":"2023-04-15T17:51:29.421413Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with Timer():\n    plt.imshow(data.surface_volume[:,:,20], cmap=\"gray\")","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:52:40.171643Z","iopub.execute_input":"2023-04-15T17:52:40.172132Z","iopub.status.idle":"2023-04-15T17:52:43.380206Z","shell.execute_reply.started":"2023-04-15T17:52:40.172092Z","shell.execute_reply":"2023-04-15T17:52:43.378863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Plot vertical slices of the surface volumes","metadata":{}},{"cell_type":"code","source":"with Timer():\n    plt.figure(figsize=(10, 1))\n    plt.imshow(data.surface_volume[2000,:,:].T, cmap=\"gray\", aspect=\"auto\")","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:52:43.382650Z","iopub.execute_input":"2023-04-15T17:52:43.383120Z","iopub.status.idle":"2023-04-15T17:52:44.552863Z","shell.execute_reply.started":"2023-04-15T17:52:43.383069Z","shell.execute_reply":"2023-04-15T17:52:44.551600Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with Timer():\n    plt.figure(figsize=(10, 1))\n    plt.imshow(data.surface_volume[:,2000,:].T, cmap=\"gray\", aspect=\"auto\")","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:54:07.751869Z","iopub.execute_input":"2023-04-15T17:54:07.752353Z","iopub.status.idle":"2023-04-15T17:54:08.738941Z","shell.execute_reply.started":"2023-04-15T17:54:07.752309Z","shell.execute_reply":"2023-04-15T17:54:08.737398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create Point Cloud\n\n## Sample from Surface Volumes","metadata":{}},{"cell_type":"code","source":"ROWS = data.surface_volume.shape[0]\nCOLS = data.surface_volume.shape[1]\nZ_DIM = data.surface_volume.shape[2] # number of volume slices\n\nN_SAMPLES = 10000","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:59:14.069202Z","iopub.execute_input":"2023-04-15T17:59:14.069750Z","iopub.status.idle":"2023-04-15T17:59:14.080390Z","shell.execute_reply.started":"2023-04-15T17:59:14.069699Z","shell.execute_reply":"2023-04-15T17:59:14.078407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with Timer():\n    # sample from valid regions of surface volume\n    c = np.ravel(data.mask).cumsum()\n    samples = np.random.uniform(low=0, \n                                high=c[-1], \n                                size=(N_SAMPLES, Z_DIM)).astype(int)\n    # get valid indexes\n    x, y = np.unravel_index(c.searchsorted(samples), data.mask.shape)\n    x, y = x[np.newaxis, ...], y[np.newaxis, ...]\n\n    # get z dimensions from surface volume locations\n    z = np.arange(0, Z_DIM)\n    z = np.tile(z, N_SAMPLES).reshape(N_SAMPLES, -1)[np.newaxis, ...]","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:59:14.347744Z","iopub.execute_input":"2023-04-15T17:59:14.349347Z","iopub.status.idle":"2023-04-15T17:59:15.945646Z","shell.execute_reply.started":"2023-04-15T17:59:14.349268Z","shell.execute_reply":"2023-04-15T17:59:15.944065Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get point cloud\nxyz = np.vstack((x, y, z))\nxyz.shape","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:59:15.948145Z","iopub.execute_input":"2023-04-15T17:59:15.948813Z","iopub.status.idle":"2023-04-15T17:59:15.964402Z","shell.execute_reply.started":"2023-04-15T17:59:15.948758Z","shell.execute_reply":"2023-04-15T17:59:15.962904Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Get Normalized Intensities","metadata":{}},{"cell_type":"code","source":"intensities = np.zeros((N_SAMPLES, Z_DIM))\n\nwith Timer():\n    for i in range(Z_DIM):\n        img = data.surface_volume[:, :, i]\n        intensities[:, i] = img[xyz[0, :, i], xyz[1, :, i]] / 65535.0\n        \n    intensities = intensities.astype(np.float32)","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:59:16.489432Z","iopub.execute_input":"2023-04-15T17:59:16.490553Z","iopub.status.idle":"2023-04-15T17:59:31.654684Z","shell.execute_reply.started":"2023-04-15T17:59:16.490481Z","shell.execute_reply":"2023-04-15T17:59:31.653296Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Sanity Check","metadata":{}},{"cell_type":"code","source":"print(xyz[:, 20, 1], intensities[20, 1])\nprint(xyz.T.reshape((-1, 3))[20 + N_SAMPLES, :], intensities.T.reshape((-1))[20 + N_SAMPLES])","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:59:37.872681Z","iopub.execute_input":"2023-04-15T17:59:37.873901Z","iopub.status.idle":"2023-04-15T17:59:37.909011Z","shell.execute_reply.started":"2023-04-15T17:59:37.873847Z","shell.execute_reply":"2023-04-15T17:59:37.907381Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Reshape and Normalize","metadata":{}},{"cell_type":"code","source":"xyz = xyz.T.reshape((-1, 3))\nxyz = xyz/xyz.max(axis=0)\n\nintensities = intensities.T.reshape((-1)).repeat((3)).reshape((-1, 3))","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:59:40.140408Z","iopub.execute_input":"2023-04-15T17:59:40.140981Z","iopub.status.idle":"2023-04-15T17:59:40.205644Z","shell.execute_reply.started":"2023-04-15T17:59:40.140936Z","shell.execute_reply":"2023-04-15T17:59:40.204509Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Get Colormap and Convert to Point Cloud","metadata":{}},{"cell_type":"code","source":"colors = plt.get_cmap('bone') # also use 'cool', 'bone'\ncolors","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:59:42.260150Z","iopub.execute_input":"2023-04-15T17:59:42.260598Z","iopub.status.idle":"2023-04-15T17:59:42.286928Z","shell.execute_reply.started":"2023-04-15T17:59:42.260558Z","shell.execute_reply":"2023-04-15T17:59:42.285561Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pcd = o3.geometry.PointCloud()\npcd.points = o3.utility.Vector3dVector(xyz) \npcd.colors = o3.utility.Vector3dVector(colors(intensities)[:, 0, :3])\n\npcd","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:59:42.817622Z","iopub.execute_input":"2023-04-15T17:59:42.818030Z","iopub.status.idle":"2023-04-15T17:59:43.862426Z","shell.execute_reply.started":"2023-04-15T17:59:42.817993Z","shell.execute_reply":"2023-04-15T17:59:43.861148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Display Point Cloud","metadata":{}},{"cell_type":"code","source":"o3.visualization.draw_plotly([pcd])","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:59:43.917393Z","iopub.execute_input":"2023-04-15T17:59:43.918191Z","iopub.status.idle":"2023-04-15T17:59:46.843859Z","shell.execute_reply.started":"2023-04-15T17:59:43.918150Z","shell.execute_reply":"2023-04-15T17:59:46.841258Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Inspect distribution at each layer \n\n### Animate intensity Histograms for each surface layer\n\nAnimation code resused from: https://www.kaggle.com/code/leonidkulyk/eda-vc-id-volume-layers-animation","metadata":{}},{"cell_type":"code","source":"!pip install celluloid -q","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:54:22.370816Z","iopub.execute_input":"2023-04-15T17:54:22.371346Z","iopub.status.idle":"2023-04-15T17:54:34.707094Z","shell.execute_reply.started":"2023-04-15T17:54:22.371298Z","shell.execute_reply":"2023-04-15T17:54:34.705454Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from celluloid import Camera","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:54:34.709934Z","iopub.execute_input":"2023-04-15T17:54:34.710476Z","iopub.status.idle":"2023-04-15T17:54:34.724999Z","shell.execute_reply.started":"2023-04-15T17:54:34.710423Z","shell.execute_reply":"2023-04-15T17:54:34.723773Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax  = plt.subplots(1, 1)\n\ncamera = Camera(fig) # define the camera that gets the fig we'll plot\nfor i in range(Z_DIM):\n    cnts, bins, _ = plt.hist(np.ravel(data.surface_volume[:, :, i][data.mask])/65535.0, bins=100);\n    ax.set_title(f\"Surfacer Layer: {i}\");\n    ax.text(\n            0.5, 1.08, f\"Surfacer Layer: {i}\", fontweight='bold', fontsize=18,\n            transform=ax.transAxes, horizontalalignment='center'\n        )\n    camera.snap() # the camera takes a snapshot of the plot\n    \n    \nplt.close(fig) # close figure\n\nanimation = camera.animate() # get plt animation\n\nfix_video_adjust = '<style> video {margin: 0px; padding: 0px; width:100%; height:auto;} </style>'\ndisplay(\n    HTML(fix_video_adjust + animation.to_html5_video())\n) # displaying the animation","metadata":{"execution":{"iopub.status.busy":"2023-04-15T17:54:34.726063Z","iopub.execute_input":"2023-04-15T17:54:34.726379Z","iopub.status.idle":"2023-04-15T17:55:54.829250Z","shell.execute_reply.started":"2023-04-15T17:54:34.726349Z","shell.execute_reply":"2023-04-15T17:55:54.827606Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There seems to be a mix of modes, especially in the lower surface volume layers. There seems to be a single mode around 0.35 for all surface volumes, while the lower volumes contain a mode with a wider spread with a center that seems to move across each layer. \n\nThe 0.35 centered mode is the dominate mode for the upper surface volumes which seem to contain less papyrus according to the volume slice cuts above. Also ccording to these histograms, the majority of the variance comes from the lower surface cuts. This can also be seen in the point cloud and the surface volume cuts.\n\nCould this dominate mode be noisy data? Let's take a closer look at with a new point cloud.","metadata":{}},{"cell_type":"markdown","source":"# Denoise Point Cloud \n\nHypothesis/Idea: The top layers do not contain useful data, they are just nosie.\n\nIn  the surface volume slices and point cloud the top layers look like they may not contain much useful data. The histograms may also support this hypothesis, since the dominate mode at 0.35 is the only mode at the top layers. \n\nLet's use the top surface volume to estimate to estimate the distribution of the noisy data. All we need is the mean and standard deviation.","metadata":{}},{"cell_type":"code","source":"mu = np.mean(np.ravel(data.surface_volume[:, :, -1][data.mask])/65535.0)\nsig = np.std(np.ravel(data.surface_volume[:, :, -1][data.mask])/65535.0)\nmu, sig","metadata":{"execution":{"iopub.status.busy":"2023-04-15T18:25:15.104688Z","iopub.execute_input":"2023-04-15T18:25:15.105179Z","iopub.status.idle":"2023-04-15T18:25:16.331550Z","shell.execute_reply.started":"2023-04-15T18:25:15.105138Z","shell.execute_reply":"2023-04-15T18:25:16.330673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's get upper and lower bounds of data to remove. We will just go up to 4 sigma for now to get the lower and upper bounds to remove. \n\nWe should also be weary of just naively truncating data. We will need to investigate this further, but that will be for another time.","metadata":{}},{"cell_type":"code","source":"lower, upper = mu - 4*sig, mu + 4*sig\nlower, upper","metadata":{"execution":{"iopub.status.busy":"2023-04-15T18:50:17.138009Z","iopub.execute_input":"2023-04-15T18:50:17.139361Z","iopub.status.idle":"2023-04-15T18:50:17.148891Z","shell.execute_reply.started":"2023-04-15T18:50:17.139292Z","shell.execute_reply":"2023-04-15T18:50:17.147591Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get intensity mask\ni_mask = (intensities[:, 0] > lower) & (intensities[:, 0] < upper)\n\n# remove 0 intensities?\ni_mask = ~i_mask & (intensities[:, 0] != 0)","metadata":{"execution":{"iopub.status.busy":"2023-04-15T18:50:17.510153Z","iopub.execute_input":"2023-04-15T18:50:17.510678Z","iopub.status.idle":"2023-04-15T18:50:17.526448Z","shell.execute_reply.started":"2023-04-15T18:50:17.510635Z","shell.execute_reply":"2023-04-15T18:50:17.524980Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pcd_2 = o3.geometry.PointCloud()\npcd_2.points = o3.utility.Vector3dVector(xyz[i_mask]) \npcd_2.colors = o3.utility.Vector3dVector(colors(intensities[i_mask])[:, 0, :3])\n\npcd_2","metadata":{"execution":{"iopub.status.busy":"2023-04-15T18:50:19.830709Z","iopub.execute_input":"2023-04-15T18:50:19.831195Z","iopub.status.idle":"2023-04-15T18:50:20.014737Z","shell.execute_reply.started":"2023-04-15T18:50:19.831154Z","shell.execute_reply":"2023-04-15T18:50:20.013272Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"o3.visualization.draw_plotly([pcd_2])","metadata":{"execution":{"iopub.status.busy":"2023-04-15T18:50:21.140246Z","iopub.execute_input":"2023-04-15T18:50:21.140756Z","iopub.status.idle":"2023-04-15T18:50:21.560048Z","shell.execute_reply.started":"2023-04-15T18:50:21.140713Z","shell.execute_reply":"2023-04-15T18:50:21.558856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's go ahead and plot plot the mean and standard dev for each surface layer. \n\nCode from: https://www.kaggle.com/code/brettolsen/fast-efficient-image-storing-and-manipulation\n","metadata":{}},{"cell_type":"code","source":"with Timer():\n    zindices = np.arange(Z_DIM)\n    means = np.zeros_like(zindices)\n    stdevs = np.zeros_like(zindices)\n    # mask = data.mask[:,:,0]\n    for z in zindices:\n        print(z, end=\"\\r\")\n        array = data.surface_volume[:,:,z][data.mask]\n        means[z] = array.mean()\n        stdevs[z] = np.std(array)\n","metadata":{"execution":{"iopub.status.busy":"2023-04-15T18:29:45.452861Z","iopub.execute_input":"2023-04-15T18:29:45.453246Z","iopub.status.idle":"2023-04-15T18:30:49.117253Z","shell.execute_reply.started":"2023-04-15T18:29:45.453212Z","shell.execute_reply":"2023-04-15T18:30:49.115871Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure()\nplt.grid()\nplt.errorbar(zindices, means, stdevs, marker='o')\nplt.plot(zindices, stdevs)","metadata":{"execution":{"iopub.status.busy":"2023-04-15T18:30:54.618999Z","iopub.execute_input":"2023-04-15T18:30:54.619474Z","iopub.status.idle":"2023-04-15T18:30:54.920963Z","shell.execute_reply.started":"2023-04-15T18:30:54.619436Z","shell.execute_reply":"2023-04-15T18:30:54.919572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It's nice to be able to see everything at once. The variance clearly reduces around layer 30, but the mean seems to level off around layer 45. In ","metadata":{}}]}