vedo.volume
Work with volumetric datasets (voxel data).
1import glob 2import os 3from weakref import ref as weak_ref_to 4 5import numpy as np 6 7import vedo.vtkclasses as vtk 8 9import vedo 10from vedo import transformations 11from vedo import utils 12from vedo.mesh import Mesh 13from vedo.core import VolumeAlgorithms 14from vedo.visual import VolumeVisual 15 16 17__docformat__ = "google" 18 19__doc__ = """ 20Work with volumetric datasets (voxel data). 21 22 23""" 24 25__all__ = ["Volume"] 26 27 28########################################################################## 29class Volume(VolumeVisual, VolumeAlgorithms): 30 """ 31 Class to describe dataset that are defined on "voxels", 32 the 3D equivalent of 2D pixels. 33 """ 34 def __init__( 35 self, 36 inputobj=None, 37 dims=None, 38 origin=None, 39 spacing=None, 40 ): 41 """ 42 This class can be initialized with a numpy object, 43 a `vtkImageData` or a list of 2D bmp files. 44 45 Arguments: 46 origin : (list) 47 set volume origin coordinates 48 spacing : (list) 49 voxel dimensions in x, y and z. 50 dims : (list) 51 specify the dimensions of the volume. 52 53 Example: 54 ```python 55 from vedo import Volume 56 vol = Volume("path/to/mydata/rec*.bmp") 57 vol.show() 58 ``` 59 60 Examples: 61 - [numpy2volume1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/numpy2volume1.py) 62 63  64 65 - [read_volume2.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/read_volume2.py) 66 67  68 69 .. note:: 70 if a `list` of values is used for `alphas` this is interpreted 71 as a transfer function along the range of the scalar. 72 """ 73 self.name = "Volume" 74 self.filename = "" 75 self.info = {} 76 self.time = 0 77 self.rendered_at = set() 78 79 self.actor = vtk.vtkVolume() 80 self.actor.retrieve_object = weak_ref_to(self) 81 self.properties = self.actor.GetProperty() 82 83 ################### 84 if isinstance(inputobj, str): 85 if "https://" in inputobj: 86 inputobj = vedo.file_io.download(inputobj, verbose=False) # fpath 87 elif os.path.isfile(inputobj): 88 self.filename = inputobj 89 else: 90 inputobj = sorted(glob.glob(inputobj)) 91 92 ################### 93 inputtype = str(type(inputobj)) 94 95 # print('Volume inputtype', inputtype, c='b') 96 97 if inputobj is None: 98 img = vtk.vtkImageData() 99 100 elif utils.is_sequence(inputobj): 101 102 if isinstance(inputobj[0], str) and ".bmp" in inputobj[0].lower(): 103 # scan sequence of BMP files 104 ima = vtk.new("ImageAppend") 105 ima.SetAppendAxis(2) 106 pb = utils.ProgressBar(0, len(inputobj)) 107 for i in pb.range(): 108 f = inputobj[i] 109 if "_rec_spr" in f: # OPT specific 110 continue 111 picr = vtk.new("BMPReader") 112 picr.SetFileName(f) 113 picr.Update() 114 mgf = vtk.new("ImageMagnitude") 115 mgf.SetInputData(picr.GetOutput()) 116 mgf.Update() 117 ima.AddInputData(mgf.GetOutput()) 118 pb.print("loading...") 119 ima.Update() 120 img = ima.GetOutput() 121 122 else: 123 124 if len(inputobj.shape) == 1: 125 varr = utils.numpy2vtk(inputobj) 126 else: 127 varr = utils.numpy2vtk(inputobj.ravel(order="F")) 128 varr.SetName("input_scalars") 129 130 img = vtk.vtkImageData() 131 if dims is not None: 132 img.SetDimensions(dims[2], dims[1], dims[0]) 133 else: 134 if len(inputobj.shape) == 1: 135 vedo.logger.error("must set dimensions (dims keyword) in Volume") 136 raise RuntimeError() 137 img.SetDimensions(inputobj.shape) 138 img.GetPointData().AddArray(varr) 139 img.GetPointData().SetActiveScalars(varr.GetName()) 140 141 elif isinstance(inputobj, vtk.vtkImageData): 142 img = inputobj 143 144 elif isinstance(inputobj, str): 145 if "https://" in inputobj: 146 inputobj = vedo.file_io.download(inputobj, verbose=False) 147 img = vedo.file_io.loadImageData(inputobj) 148 self.filename = inputobj 149 150 else: 151 vedo.logger.error(f"cannot understand input type {inputtype}") 152 return 153 154 if dims is not None: 155 img.SetDimensions(dims) 156 157 if origin is not None: 158 img.SetOrigin(origin) 159 160 if spacing is not None: 161 img.SetSpacing(spacing) 162 163 self.dataset = img 164 self.transform = None 165 166 ##################################### 167 mapper = vtk.new("SmartVolumeMapper") 168 mapper.SetInputData(img) 169 self.actor.SetMapper(mapper) 170 171 if img.GetPointData().GetScalars(): 172 if img.GetPointData().GetScalars().GetNumberOfComponents() == 1: 173 self.properties.SetShade(True) 174 self.properties.SetInterpolationType(1) 175 self.cmap("RdBu_r") 176 self.alpha([0.0, 0.0, 0.2, 0.4, 0.8, 1.0]) 177 self.alpha_gradient(None) 178 self.properties.SetScalarOpacityUnitDistance(1.0) 179 180 self.pipeline = utils.OperationNode( 181 "Volume", comment=f"dims={tuple(self.dimensions())}", c="#4cc9f0" 182 ) 183 ####################################################################### 184 185 @property 186 def mapper(self): 187 """Return the underlying `vtkMapper` object.""" 188 return self.actor.GetMapper() 189 190 @mapper.setter 191 def mapper(self, mapper): 192 """ 193 Set the underlying `vtkMapper` object. 194 195 Arguments: 196 mapper : (str, vtkMapper) 197 either 'gpu', 'opengl_gpu', 'fixed' or 'smart' 198 """ 199 if isinstance(mapper, 200 (vtk.get_class("Mapper"), 201 vtk.get_class("ImageResliceMapper", 202 ) )): 203 pass 204 elif "gpu" in mapper: 205 mapper = vtk.new("GPUVolumeRayCastMapper") 206 elif "opengl_gpu" in mapper: 207 mapper = vtk.new("OpenGLGPUVolumeRayCastMapper") 208 elif "smart" in mapper: 209 mapper = vtk.new("SmartVolumeMapper") 210 elif "fixed" in mapper: 211 mapper = vtk.new("FixedPointVolumeRayCastMapper") 212 else: 213 print("Error unknown mapper type", [mapper]) 214 raise RuntimeError() 215 self.actor.SetMapper(mapper) 216 217 def c(self, *args, **kwargs): 218 """Deprecated. Use `Volume.cmap()` instead.""" 219 vedo.logger.warning("Volume.c() is deprecated, use Volume.cmap() instead") 220 return self.cmap(*args, **kwargs) 221 222 def _update(self, data): 223 self.dataset = data 224 self.mapper.SetInputData(data) 225 self.dataset.GetPointData().Modified() 226 self.mapper.Modified() 227 self.mapper.Update() 228 return self 229 230 def __str__(self): 231 """Print a summary for the Volume object.""" 232 module = self.__class__.__module__ 233 name = self.__class__.__name__ 234 out = vedo.printc( 235 f"{module}.{name} at ({hex(self.memory_address())})".ljust(75), 236 c="c", bold=True, invert=True, return_string=True, 237 ) 238 out += "\x1b[0m\x1b[36;1m" 239 240 out+= "name".ljust(14) + ": " + str(self.name) + "\n" 241 if self.filename: 242 out+= "filename".ljust(14) + ": " + str(self.filename) + "\n" 243 244 out+= "dimensions".ljust(14) + ": " + str(self.shape) + "\n" 245 246 out+= "origin".ljust(14) + ": " 247 out+= utils.precision(self.origin(), 6) + "\n" 248 249 out+= "center".ljust(14) + ": " 250 out+= utils.precision(self.center(), 6) + "\n" 251 252 out+= "spacing".ljust(14) + ": " 253 out+= utils.precision(self.spacing(), 6) + "\n" 254 255 bnds = self.bounds() 256 bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3) 257 by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3) 258 bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3) 259 out+= "bounds".ljust(14) + ":" 260 out+= " x=(" + bx1 + ", " + bx2 + ")," 261 out+= " y=(" + by1 + ", " + by2 + ")," 262 out+= " z=(" + bz1 + ", " + bz2 + ")\n" 263 264 out+= "memory size".ljust(14) + ": " 265 out+= str(int(self.dataset.GetActualMemorySize()/1024+0.5))+" MB\n" 266 267 st = self.dataset.GetScalarTypeAsString() 268 out+= "scalar size".ljust(14) + ": " 269 out+= str(self.dataset.GetScalarSize()) + f" bytes ({st})\n" 270 out+= "scalar range".ljust(14) + ": " 271 out+= str(self.dataset.GetScalarRange()) + "\n" 272 273 #utils.print_histogram(self, logscale=True, bins=8, height=15, c="b", bold=True) 274 return out.rstrip() + "\x1b[0m" 275 276 def _repr_html_(self): 277 """ 278 HTML representation of the Volume object for Jupyter Notebooks. 279 280 Returns: 281 HTML text with the image and some properties. 282 """ 283 import io 284 import base64 285 from PIL import Image 286 287 library_name = "vedo.volume.Volume" 288 help_url = "https://vedo.embl.es/docs/vedo/volume.html" 289 290 arr = self.thumbnail(azimuth=0, elevation=-60, zoom=1.4, axes=True) 291 292 im = Image.fromarray(arr) 293 buffered = io.BytesIO() 294 im.save(buffered, format="PNG", quality=100) 295 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 296 url = "data:image/png;base64," + encoded 297 image = f"<img src='{url}'></img>" 298 299 # statisitics 300 bounds = "<br/>".join( 301 [ 302 utils.precision(min_x, 4) + " ... " + utils.precision(max_x, 4) 303 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 304 ] 305 ) 306 307 help_text = "" 308 if self.name: 309 help_text += f"<b> {self.name}:   </b>" 310 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 311 if self.filename: 312 dots = "" 313 if len(self.filename) > 30: 314 dots = "..." 315 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 316 317 pdata = "" 318 if self.dataset.GetPointData().GetScalars(): 319 if self.dataset.GetPointData().GetScalars().GetName(): 320 name = self.dataset.GetPointData().GetScalars().GetName() 321 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 322 323 cdata = "" 324 if self.dataset.GetCellData().GetScalars(): 325 if self.dataset.GetCellData().GetScalars().GetName(): 326 name = self.dataset.GetCellData().GetScalars().GetName() 327 cdata = "<tr><td><b> voxel data array </b></td><td>" + name + "</td></tr>" 328 329 img = self.dataset 330 331 allt = [ 332 "<table>", 333 "<tr>", 334 "<td>", 335 image, 336 "</td>", 337 "<td style='text-align: center; vertical-align: center;'><br/>", 338 help_text, 339 "<table>", 340 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 341 "<tr><td><b> dimensions </b></td><td>" + str(img.GetDimensions()) + "</td></tr>", 342 "<tr><td><b> voxel spacing </b></td><td>" 343 + utils.precision(img.GetSpacing(), 3) 344 + "</td></tr>", 345 "<tr><td><b> in memory size </b></td><td>" 346 + str(int(img.GetActualMemorySize() / 1024)) 347 + "MB</td></tr>", 348 pdata, 349 cdata, 350 "<tr><td><b> scalar range </b></td><td>" 351 + utils.precision(img.GetScalarRange(), 4) 352 + "</td></tr>", 353 "</table>", 354 "</table>", 355 ] 356 return "\n".join(allt) 357 358 def copy(self, deep=True): 359 """Return a copy of the Volume. Alias of `clone()`.""" 360 return self.clone(deep=deep) 361 362 def clone(self, deep=True): 363 """Return a clone copy of the Volume. Alias of `copy()`.""" 364 if deep: 365 newimg = vtk.vtkImageData() 366 newimg.CopyStructure(self.dataset) 367 newimg.CopyAttributes(self.dataset) 368 newvol = Volume(newimg) 369 else: 370 newvol = Volume(self.dataset) 371 372 prop = vtk.vtkVolumeProperty() 373 prop.DeepCopy(self.properties) 374 newvol.actor.SetProperty(prop) 375 newvol.properties = prop 376 377 newvol.pipeline = utils.OperationNode("clone", parents=[self], c="#bbd0ff", shape="diamond") 378 return newvol 379 380 def component_weight(self, i, weight): 381 """Set the scalar component weight in range [0,1].""" 382 self.properties.SetComponentWeight(i, weight) 383 return self 384 385 def xslice(self, i): 386 """Extract the slice at index `i` of volume along x-axis.""" 387 vslice = vtk.new("ImageDataGeometryFilter") 388 vslice.SetInputData(self.dataset) 389 nx, ny, nz = self.dataset.GetDimensions() 390 if i > nx - 1: 391 i = nx - 1 392 vslice.SetExtent(i, i, 0, ny, 0, nz) 393 vslice.Update() 394 m = Mesh(vslice.GetOutput()) 395 m.pipeline = utils.OperationNode(f"xslice {i}", parents=[self], c="#4cc9f0:#e9c46a") 396 return m 397 398 def yslice(self, j): 399 """Extract the slice at index `j` of volume along y-axis.""" 400 vslice = vtk.new("ImageDataGeometryFilter") 401 vslice.SetInputData(self.dataset) 402 nx, ny, nz = self.dataset.GetDimensions() 403 if j > ny - 1: 404 j = ny - 1 405 vslice.SetExtent(0, nx, j, j, 0, nz) 406 vslice.Update() 407 m = Mesh(vslice.GetOutput()) 408 m.pipeline = utils.OperationNode(f"yslice {j}", parents=[self], c="#4cc9f0:#e9c46a") 409 return m 410 411 def zslice(self, k): 412 """Extract the slice at index `i` of volume along z-axis.""" 413 vslice = vtk.new("ImageDataGeometryFilter") 414 vslice.SetInputData(self.dataset) 415 nx, ny, nz = self.dataset.GetDimensions() 416 if k > nz - 1: 417 k = nz - 1 418 vslice.SetExtent(0, nx, 0, ny, k, k) 419 vslice.Update() 420 m = Mesh(vslice.GetOutput()) 421 m.pipeline = utils.OperationNode(f"zslice {k}", parents=[self], c="#4cc9f0:#e9c46a") 422 return m 423 424 def slice_plane(self, origin=(0, 0, 0), normal=(1, 1, 1), autocrop=False): 425 """ 426 Extract the slice along a given plane position and normal. 427 428 Example: 429 - [slice_plane1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/slice_plane1.py) 430 431  432 """ 433 reslice = vtk.new("ImageReslice") 434 reslice.SetInputData(self.dataset) 435 reslice.SetOutputDimensionality(2) 436 newaxis = utils.versor(normal) 437 pos = np.array(origin) 438 initaxis = (0, 0, 1) 439 crossvec = np.cross(initaxis, newaxis) 440 angle = np.arccos(np.dot(initaxis, newaxis)) 441 T = vtk.vtkTransform() 442 T.PostMultiply() 443 T.RotateWXYZ(np.rad2deg(angle), crossvec) 444 T.Translate(pos) 445 M = T.GetMatrix() 446 reslice.SetResliceAxes(M) 447 reslice.SetInterpolationModeToLinear() 448 reslice.SetAutoCropOutput(not autocrop) 449 reslice.Update() 450 vslice = vtk.new("ImageDataGeometryFilter") 451 vslice.SetInputData(reslice.GetOutput()) 452 vslice.Update() 453 msh = Mesh(vslice.GetOutput()) 454 msh.apply_transform(T) 455 msh.pipeline = utils.OperationNode( 456 "slice_plane", parents=[self], c="#4cc9f0:#e9c46a") 457 return msh 458 459 def slab(self, slice_range=(), axis='z', operation="mean"): 460 """ 461 Extract a slab from a `Volume` by combining 462 all of the slices of an image to create a single slice. 463 464 Returns a `Mesh` containing metadata which 465 can be accessed with e.g. `mesh.metadata["slab_range"]`. 466 467 Metadata: 468 slab_range : (list) 469 contains the range of slices extracted 470 slab_axis : (str) 471 contains the axis along which the slab was extracted 472 slab_operation : (str) 473 contains the operation performed on the slab 474 slab_bounding_box : (list) 475 contains the bounding box of the slab 476 477 Arguments: 478 slice_range : (list) 479 range of slices to extract 480 axis : (str) 481 axis along which to extract the slab 482 operation : (str) 483 operation to perform on the slab, 484 allowed values are: "sum", "min", "max", "mean". 485 """ 486 if len(slice_range) != 2: 487 vedo.logger.error("in slab(): slice_range is empty or invalid") 488 raise ValueError() 489 490 islab = vtk.new("ImageSlab") 491 islab.SetInputData(self.dataset) 492 493 if operation in ["+", "add", "sum"]: 494 islab.SetOperationToSum() 495 elif "min" in operation: 496 islab.SetOperationToMin() 497 elif "max" in operation: 498 islab.SetOperationToMax() 499 elif "mean" in operation: 500 islab.SetOperationToMean() 501 else: 502 vedo.logger.error(f"in slab(): unknown operation {operation}") 503 raise ValueError() 504 505 dims = self.dimensions() 506 if axis == 'x': 507 islab.SetOrientationToX() 508 if slice_range[0] > dims[0]-1: 509 slice_range[0] = dims[0]-1 510 if slice_range[1] > dims[0]-1: 511 slice_range[1] = dims[0]-1 512 elif axis == 'y': 513 islab.SetOrientationToY() 514 if slice_range[0] > dims[1]-1: 515 slice_range[0] = dims[1]-1 516 if slice_range[1] > dims[1]-1: 517 slice_range[1] = dims[1]-1 518 elif axis == 'z': 519 islab.SetOrientationToZ() 520 if slice_range[0] > dims[2]-1: 521 slice_range[0] = dims[2]-1 522 if slice_range[1] > dims[2]-1: 523 slice_range[1] = dims[2]-1 524 else: 525 vedo.logger.error(f"Error in slab(): unknown axis {axis}") 526 raise RuntimeError() 527 528 islab.SetSliceRange(slice_range) 529 islab.Update() 530 531 msh = Mesh(islab.GetOutput()).lighting('off') 532 msh.mapper.SetLookupTable(utils.ctf2lut(self, msh)) 533 msh.mapper.SetScalarRange(self.scalar_range()) 534 535 msh.metadata["slab_range"] = slice_range 536 msh.metadata["slab_axis"] = axis 537 msh.metadata["slab_operation"] = operation 538 539 # compute bounds of slab 540 origin = self.origin() 541 spacing = self.spacing() 542 if axis == 'x': 543 msh.metadata["slab_bounding_box"] = [ 544 origin[0] + slice_range[0]*spacing[0], 545 origin[0] + slice_range[1]*spacing[0], 546 origin[1], 547 origin[1] + dims[1]*spacing[1], 548 origin[2], 549 origin[2] + dims[2]*spacing[2], 550 ] 551 elif axis == 'y': 552 msh.metadata["slab_bounding_box"] = [ 553 origin[0], 554 origin[0] + dims[0]*spacing[0], 555 origin[1] + slice_range[0]*spacing[1], 556 origin[1] + slice_range[1]*spacing[1], 557 origin[2], 558 origin[2] + dims[2]*spacing[2], 559 ] 560 elif axis == 'z': 561 msh.metadata["slab_bounding_box"] = [ 562 origin[0], 563 origin[0] + dims[0]*spacing[0], 564 origin[1], 565 origin[1] + dims[1]*spacing[1], 566 origin[2] + slice_range[0]*spacing[2], 567 origin[2] + slice_range[1]*spacing[2], 568 ] 569 570 msh.pipeline = utils.OperationNode( 571 f"slab{slice_range}", 572 comment=f"axis={axis}, operation={operation}", 573 parents=[self], 574 c="#4cc9f0:#e9c46a", 575 ) 576 msh.name = "SlabMesh" 577 return msh 578 579 580 def warp(self, source, target, sigma=1, mode="3d", fit=True): 581 """ 582 Warp volume scalars within a Volume by specifying 583 source and target sets of points. 584 585 Arguments: 586 source : (Points, list) 587 the list of source points 588 target : (Points, list) 589 the list of target points 590 fit : (bool) 591 fit/adapt the old bounding box to the warped geometry 592 """ 593 if isinstance(source, vedo.Points): 594 source = source.vertices 595 if isinstance(target, vedo.Points): 596 target = target.vertices 597 598 NLT = transformations.NonLinearTransform() 599 NLT.source_points = source 600 NLT.target_points = target 601 NLT.sigma = sigma 602 NLT.mode = mode 603 604 self.apply_transform(NLT, fit=fit) 605 self.pipeline = utils.OperationNode("warp", parents=[self], c="#4cc9f0") 606 return self 607 608 def apply_transform(self, T, fit=True, interpolation="cubic"): 609 """ 610 Apply a transform to the scalars in the volume. 611 612 Arguments: 613 T : (LinearTransform, NonLinearTransform) 614 The transformation to be applied 615 fit : (bool) 616 fit/adapt the old bounding box to the modified geometry 617 interpolation : (str) 618 one of the following: "linear", "nearest", "cubic" 619 """ 620 TI = T.compute_inverse() 621 reslice = vtk.new("ImageReslice") 622 reslice.SetInputData(self.dataset) 623 reslice.SetResliceTransform(TI.T) 624 reslice.SetOutputDimensionality(3) 625 if "lin" in interpolation.lower(): 626 reslice.SetInterpolationModeToLinear() 627 elif "near" in interpolation.lower(): 628 reslice.SetInterpolationModeToNearestNeighbor() 629 elif "cubic" in interpolation.lower(): 630 reslice.SetInterpolationModeToCubic() 631 else: 632 vedo.logger.error( 633 f"in apply_transform: unknown interpolation mode {interpolation}") 634 raise ValueError() 635 reslice.SetAutoCropOutput(fit) 636 reslice.Update() 637 self._update(reslice.GetOutput()) 638 self.transform = T 639 self.pipeline = utils.OperationNode( 640 "apply_transform", parents=[self], c="#4cc9f0") 641 return self 642 643 def imagedata(self): 644 """ 645 DEPRECATED: 646 Use `Volume.dataset` instead. 647 648 Return the underlying `vtkImagaData` object. 649 """ 650 print("Volume.imagedata() is deprecated, use Volume.dataset instead") 651 return self.dataset 652 653 def modified(self): 654 """ 655 Mark the object as modified. 656 657 Example: 658 659 - [numpy2volume0.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/numpy2volume0.py) 660 """ 661 self.dataset.GetPointData().GetScalars().Modified() 662 return self 663 664 def tonumpy(self): 665 """ 666 Get read-write access to voxels of a Volume object as a numpy array. 667 668 When you set values in the output image, you don't want numpy to reallocate the array 669 but instead set values in the existing array, so use the [:] operator. 670 671 Example: 672 `arr[:] = arr*2 + 15` 673 674 If the array is modified add a call to: 675 `volume.modified()` 676 when all your modifications are completed. 677 """ 678 narray_shape = tuple(reversed(self.dataset.GetDimensions())) 679 680 scals = self.dataset.GetPointData().GetScalars() 681 comps = scals.GetNumberOfComponents() 682 if comps == 1: 683 narray = utils.vtk2numpy(scals).reshape(narray_shape) 684 narray = np.transpose(narray, axes=[2, 1, 0]) 685 else: 686 narray = utils.vtk2numpy(scals).reshape(*narray_shape, comps) 687 narray = np.transpose(narray, axes=[2, 1, 0, 3]) 688 689 # narray = utils.vtk2numpy(self.dataset.GetPointData().GetScalars()).reshape(narray_shape) 690 # narray = np.transpose(narray, axes=[2, 1, 0]) 691 692 return narray 693 694 @property 695 def shape(self): 696 """Return the nr. of voxels in the 3 dimensions.""" 697 return np.array(self.dataset.GetDimensions()) 698 699 def dimensions(self): 700 """Return the nr. of voxels in the 3 dimensions.""" 701 return np.array(self.dataset.GetDimensions()) 702 703 def scalar_range(self): 704 """Return the range of the scalar values.""" 705 return np.array(self.dataset.GetScalarRange()) 706 707 def spacing(self, s=None): 708 """Set/get the voxels size in the 3 dimensions.""" 709 if s is not None: 710 self.dataset.SetSpacing(s) 711 return self 712 return np.array(self.dataset.GetSpacing()) 713 714 def origin(self, s=None): 715 """ 716 Set/get the origin of the volumetric dataset. 717 718 The origin is the position in world coordinates of the point index (0,0,0). 719 This point does not have to be part of the dataset, in other words, 720 the dataset extent does not have to start at (0,0,0) and the origin 721 can be outside of the dataset bounding box. 722 The origin plus spacing determine the position in space of the points. 723 """ 724 if s is not None: 725 self.dataset.SetOrigin(s) 726 return self 727 return np.array(self.dataset.GetOrigin()) 728 729 def center(self): 730 """Get the center of the volumetric dataset.""" 731 # note that this does not have the set method like origin and spacing 732 return np.array(self.dataset.GetCenter()) 733 734 def get_cell_from_ijk(self, ijk): 735 """ 736 Get the voxel id number at the given ijk coordinates. 737 738 Arguments: 739 ijk : (list) 740 the ijk coordinates of the voxel 741 """ 742 return self.ComputeCellId(ijk) 743 744 def get_point_from_ijk(self, ijk): 745 """ 746 Get the point id number at the given ijk coordinates. 747 748 Arguments: 749 ijk : (list) 750 the ijk coordinates of the voxel 751 """ 752 return self.ComputePointId(ijk) 753 754 def permute_axes(self, x, y, z): 755 """ 756 Reorder the axes of the Volume by specifying 757 the input axes which are supposed to become the new X, Y, and Z. 758 """ 759 imp = vtk.new("ImagePermute") 760 imp.SetFilteredAxes(x, y, z) 761 imp.SetInputData(self.dataset) 762 imp.Update() 763 self._update(imp.GetOutput()) 764 self.pipeline = utils.OperationNode( 765 f"permute_axes({(x,y,z)})", parents=[self], c="#4cc9f0" 766 ) 767 return self 768 769 def resample(self, new_spacing, interpolation=1): 770 """ 771 Resamples a `Volume` to be larger or smaller. 772 773 This method modifies the spacing of the input. 774 Linear interpolation is used to resample the data. 775 776 Arguments: 777 new_spacing : (list) 778 a list of 3 new spacings for the 3 axes 779 interpolation : (int) 780 0=nearest_neighbor, 1=linear, 2=cubic 781 """ 782 rsp = vtk.new("ImageResample") 783 oldsp = self.spacing() 784 for i in range(3): 785 if oldsp[i] != new_spacing[i]: 786 rsp.SetAxisOutputSpacing(i, new_spacing[i]) 787 rsp.InterpolateOn() 788 rsp.SetInterpolationMode(interpolation) 789 rsp.OptimizationOn() 790 rsp.Update() 791 self._update(rsp.GetOutput()) 792 self.pipeline = utils.OperationNode( 793 "resample", comment=f"spacing: {tuple(new_spacing)}", parents=[self], c="#4cc9f0" 794 ) 795 return self 796 797 798 def threshold(self, above=None, below=None, replace=None, replace_value=None): 799 """ 800 Binary or continuous volume thresholding. 801 Find the voxels that contain a value above/below the input values 802 and replace them with a new value (default is 0). 803 """ 804 th = vtk.new("ImageThreshold") 805 th.SetInputData(self.dataset) 806 807 # sanity checks 808 if above is not None and below is not None: 809 if above == below: 810 return self 811 if above > below: 812 vedo.logger.warning("in volume.threshold(), above > below, skip.") 813 return self 814 815 ## cases 816 if below is not None and above is not None: 817 th.ThresholdBetween(above, below) 818 819 elif above is not None: 820 th.ThresholdByUpper(above) 821 822 elif below is not None: 823 th.ThresholdByLower(below) 824 825 ## 826 if replace is not None: 827 th.SetReplaceIn(True) 828 th.SetInValue(replace) 829 else: 830 th.SetReplaceIn(False) 831 832 if replace_value is not None: 833 th.SetReplaceOut(True) 834 th.SetOutValue(replace_value) 835 else: 836 th.SetReplaceOut(False) 837 838 th.Update() 839 self._update(th.GetOutput()) 840 self.pipeline = utils.OperationNode("threshold", parents=[self], c="#4cc9f0") 841 return self 842 843 def crop(self, left=None, right=None, back=None, front=None, bottom=None, top=None, VOI=()): 844 """ 845 Crop a `Volume` object. 846 847 Arguments: 848 left : (float) 849 fraction to crop from the left plane (negative x) 850 right : (float) 851 fraction to crop from the right plane (positive x) 852 back : (float) 853 fraction to crop from the back plane (negative y) 854 front : (float) 855 fraction to crop from the front plane (positive y) 856 bottom : (float) 857 fraction to crop from the bottom plane (negative z) 858 top : (float) 859 fraction to crop from the top plane (positive z) 860 VOI : (list) 861 extract Volume Of Interest expressed in voxel numbers 862 863 Example: 864 `vol.crop(VOI=(xmin, xmax, ymin, ymax, zmin, zmax)) # all integers nrs` 865 """ 866 extractVOI = vtk.new("ExtractVOI") 867 extractVOI.SetInputData(self.dataset) 868 869 if VOI: 870 extractVOI.SetVOI(VOI) 871 else: 872 d = self.dataset.GetDimensions() 873 bx0, bx1, by0, by1, bz0, bz1 = 0, d[0]-1, 0, d[1]-1, 0, d[2]-1 874 if left is not None: bx0 = int((d[0]-1)*left) 875 if right is not None: bx1 = int((d[0]-1)*(1-right)) 876 if back is not None: by0 = int((d[1]-1)*back) 877 if front is not None: by1 = int((d[1]-1)*(1-front)) 878 if bottom is not None: bz0 = int((d[2]-1)*bottom) 879 if top is not None: bz1 = int((d[2]-1)*(1-top)) 880 extractVOI.SetVOI(bx0, bx1, by0, by1, bz0, bz1) 881 extractVOI.Update() 882 self._update(extractVOI.GetOutput()) 883 884 self.pipeline = utils.OperationNode( 885 "crop", parents=[self], c="#4cc9f0", comment=f"dims={tuple(self.dimensions())}" 886 ) 887 return self 888 889 def append(self, volumes, axis="z", preserve_extents=False): 890 """ 891 Take the components from multiple inputs and merges them into one output. 892 Except for the append axis, all inputs must have the same extent. 893 All inputs must have the same number of scalar components. 894 The output has the same origin and spacing as the first input. 895 The origin and spacing of all other inputs are ignored. 896 All inputs must have the same scalar type. 897 898 Arguments: 899 axis : (int, str) 900 axis expanded to hold the multiple images 901 preserve_extents : (bool) 902 if True, the extent of the inputs is used to place 903 the image in the output. The whole extent of the output is the union of the input 904 whole extents. Any portion of the output not covered by the inputs is set to zero. 905 The origin and spacing is taken from the first input. 906 907 Example: 908 ```python 909 from vedo import Volume, dataurl 910 vol = Volume(dataurl+'embryo.tif') 911 vol.append(vol, axis='x').show().close() 912 ``` 913  914 """ 915 ima = vtk.new("ImageAppend") 916 ima.SetInputData(self.dataset) 917 if not utils.is_sequence(volumes): 918 volumes = [volumes] 919 for volume in volumes: 920 if isinstance(volume, vtk.vtkImageData): 921 ima.AddInputData(volume) 922 else: 923 ima.AddInputData(volume.dataset) 924 ima.SetPreserveExtents(preserve_extents) 925 if axis == "x": 926 axis = 0 927 elif axis == "y": 928 axis = 1 929 elif axis == "z": 930 axis = 2 931 ima.SetAppendAxis(axis) 932 ima.Update() 933 self._update(ima.GetOutput()) 934 935 self.pipeline = utils.OperationNode( 936 "append", 937 parents=[self, *volumes], 938 c="#4cc9f0", 939 comment=f"dims={tuple(self.dimensions())}", 940 ) 941 return self 942 943 def pad(self, voxels=10, value=0): 944 """ 945 Add the specified number of voxels at the `Volume` borders. 946 Voxels can be a list formatted as `[nx0, nx1, ny0, ny1, nz0, nz1]`. 947 948 Arguments: 949 voxels : (int, list) 950 number of voxels to be added (or a list of length 4) 951 value : (int) 952 intensity value (gray-scale color) of the padding 953 954 Example: 955 ```python 956 from vedo import Volume, dataurl, show 957 iso = Volume(dataurl+'embryo.tif').isosurface() 958 vol = iso.binarize(spacing=(100, 100, 100)).pad(10) 959 vol.dilate([15,15,15]) 960 show(iso, vol.isosurface(), N=2, axes=1) 961 ``` 962  963 """ 964 x0, x1, y0, y1, z0, z1 = self.dataset.GetExtent() 965 pf = vtk.new("ImageConstantPad") 966 pf.SetInputData(self.dataset) 967 pf.SetConstant(value) 968 if utils.is_sequence(voxels): 969 pf.SetOutputWholeExtent( 970 x0 - voxels[0], x1 + voxels[1], 971 y0 - voxels[2], y1 + voxels[3], 972 z0 - voxels[4], z1 + voxels[5], 973 ) 974 else: 975 pf.SetOutputWholeExtent( 976 x0 - voxels, x1 + voxels, 977 y0 - voxels, y1 + voxels, 978 z0 - voxels, z1 + voxels, 979 ) 980 pf.Update() 981 self._update(pf.GetOutput()) 982 self.pipeline = utils.OperationNode( 983 "pad", comment=f"{voxels} voxels", parents=[self], c="#f28482" 984 ) 985 return self 986 987 def resize(self, newdims): 988 """Increase or reduce the number of voxels of a Volume with interpolation.""" 989 rsz = vtk.new("ImageResize") 990 rsz.SetResizeMethodToOutputDimensions() 991 rsz.SetInputData(self.dataset) 992 rsz.SetOutputDimensions(newdims) 993 rsz.Update() 994 self.dataset = rsz.GetOutput() 995 self._update(self.dataset) 996 self.pipeline = utils.OperationNode( 997 "resize", parents=[self], c="#4cc9f0", comment=f"dims={tuple(self.dimensions())}" 998 ) 999 return self 1000 1001 def normalize(self): 1002 """Normalize that scalar components for each point.""" 1003 norm = vtk.new("ImageNormalize") 1004 norm.SetInputData(self.dataset) 1005 norm.Update() 1006 self._update(norm.GetOutput()) 1007 self.pipeline = utils.OperationNode("normalize", parents=[self], c="#4cc9f0") 1008 return self 1009 1010 def mirror(self, axis="x"): 1011 """ 1012 Mirror flip along one of the cartesian axes. 1013 """ 1014 img = self.dataset 1015 1016 ff = vtk.new("ImageFlip") 1017 ff.SetInputData(img) 1018 if axis.lower() == "x": 1019 ff.SetFilteredAxis(0) 1020 elif axis.lower() == "y": 1021 ff.SetFilteredAxis(1) 1022 elif axis.lower() == "z": 1023 ff.SetFilteredAxis(2) 1024 else: 1025 vedo.logger.error("mirror must be set to either x, y, z or n") 1026 raise RuntimeError() 1027 ff.Update() 1028 self._update(ff.GetOutput()) 1029 self.pipeline = utils.OperationNode(f"mirror {axis}", parents=[self], c="#4cc9f0") 1030 return self 1031 1032 def operation(self, operation, volume2=None): 1033 """ 1034 Perform operations with `Volume` objects. 1035 Keyword `volume2` can be a constant float. 1036 1037 Possible operations are: 1038 ``` 1039 +, -, /, 1/x, sin, cos, exp, log, 1040 abs, **2, sqrt, min, max, atan, atan2, median, 1041 mag, dot, gradient, divergence, laplacian. 1042 ``` 1043 1044 Examples: 1045 - [volume_operations.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/volume_operations.py) 1046 """ 1047 op = operation.lower() 1048 image1 = self.dataset 1049 1050 mf = None 1051 if op in ["median"]: 1052 mf = vtk.new("ImageMedian3D") 1053 mf.SetInputData(image1) 1054 elif op in ["mag"]: 1055 mf = vtk.new("ImageMagnitude") 1056 mf.SetInputData(image1) 1057 elif op in ["dot", "dotproduct"]: 1058 mf = vtk.new("ImageDotProduct") 1059 mf.SetInput1Data(image1) 1060 mf.SetInput2Data(volume2.dataset) 1061 elif op in ["grad", "gradient"]: 1062 mf = vtk.new("ImageGradient") 1063 mf.SetDimensionality(3) 1064 mf.SetInputData(image1) 1065 elif op in ["div", "divergence"]: 1066 mf = vtk.new("ImageDivergence") 1067 mf.SetInputData(image1) 1068 elif op in ["laplacian"]: 1069 mf = vtk.new("ImageLaplacian") 1070 mf.SetDimensionality(3) 1071 mf.SetInputData(image1) 1072 1073 if mf is not None: 1074 mf.Update() 1075 vol = Volume(mf.GetOutput()) 1076 vol.pipeline = utils.OperationNode( 1077 "operation", comment=f"{op}", parents=[self], c="#4cc9f0", shape="cylinder" 1078 ) 1079 return vol ########################### 1080 1081 mat = vtk.new("ImageMathematics") 1082 mat.SetInput1Data(image1) 1083 1084 K = None 1085 1086 if utils.is_number(volume2): 1087 K = volume2 1088 mat.SetConstantK(K) 1089 mat.SetConstantC(K) 1090 1091 elif volume2 is not None: # assume image2 is a constant value 1092 mat.SetInput2Data(volume2.dataset) 1093 1094 # ########################### 1095 if op in ["+", "add", "plus"]: 1096 if K: 1097 mat.SetOperationToAddConstant() 1098 else: 1099 mat.SetOperationToAdd() 1100 1101 elif op in ["-", "subtract", "minus"]: 1102 if K: 1103 mat.SetConstantC(-float(K)) 1104 mat.SetOperationToAddConstant() 1105 else: 1106 mat.SetOperationToSubtract() 1107 1108 elif op in ["*", "multiply", "times"]: 1109 if K: 1110 mat.SetOperationToMultiplyByK() 1111 else: 1112 mat.SetOperationToMultiply() 1113 1114 elif op in ["/", "divide"]: 1115 if K: 1116 mat.SetConstantK(1.0 / K) 1117 mat.SetOperationToMultiplyByK() 1118 else: 1119 mat.SetOperationToDivide() 1120 1121 elif op in ["1/x", "invert"]: 1122 mat.SetOperationToInvert() 1123 elif op in ["sin"]: 1124 mat.SetOperationToSin() 1125 elif op in ["cos"]: 1126 mat.SetOperationToCos() 1127 elif op in ["exp"]: 1128 mat.SetOperationToExp() 1129 elif op in ["log"]: 1130 mat.SetOperationToLog() 1131 elif op in ["abs"]: 1132 mat.SetOperationToAbsoluteValue() 1133 elif op in ["**2", "square"]: 1134 mat.SetOperationToSquare() 1135 elif op in ["sqrt", "sqr"]: 1136 mat.SetOperationToSquareRoot() 1137 elif op in ["min"]: 1138 mat.SetOperationToMin() 1139 elif op in ["max"]: 1140 mat.SetOperationToMax() 1141 elif op in ["atan"]: 1142 mat.SetOperationToATAN() 1143 elif op in ["atan2"]: 1144 mat.SetOperationToATAN2() 1145 else: 1146 vedo.logger.error(f"unknown operation {operation}") 1147 raise RuntimeError() 1148 mat.Update() 1149 1150 self._update(mat.GetOutput()) 1151 1152 self.pipeline = utils.OperationNode( 1153 "operation", comment=f"{op}", parents=[self, volume2], shape="cylinder", c="#4cc9f0" 1154 ) 1155 return self 1156 1157 def frequency_pass_filter(self, low_cutoff=None, high_cutoff=None, order=1): 1158 """ 1159 Low-pass and high-pass filtering become trivial in the frequency domain. 1160 A portion of the pixels/voxels are simply masked or attenuated. 1161 This function applies a high pass Butterworth filter that attenuates the 1162 frequency domain image. 1163 1164 The gradual attenuation of the filter is important. 1165 A simple high-pass filter would simply mask a set of pixels in the frequency domain, 1166 but the abrupt transition would cause a ringing effect in the spatial domain. 1167 1168 Arguments: 1169 low_cutoff : (list) 1170 the cutoff frequencies for x, y and z 1171 high_cutoff : (list) 1172 the cutoff frequencies for x, y and z 1173 order : (int) 1174 order determines sharpness of the cutoff curve 1175 """ 1176 # https://lorensen.github.io/VTKExamples/site/Cxx/ImageProcessing/IdealHighPass 1177 fft = vtk.new("ImageFFT") 1178 fft.SetInputData(self.dataset) 1179 fft.Update() 1180 out = fft.GetOutput() 1181 1182 if high_cutoff: 1183 blp = vtk.new("ImageButterworthLowPass") 1184 blp.SetInputData(out) 1185 blp.SetCutOff(high_cutoff) 1186 blp.SetOrder(order) 1187 blp.Update() 1188 out = blp.GetOutput() 1189 1190 if low_cutoff: 1191 bhp = vtk.new("ImageButterworthHighPass") 1192 bhp.SetInputData(out) 1193 bhp.SetCutOff(low_cutoff) 1194 bhp.SetOrder(order) 1195 bhp.Update() 1196 out = bhp.GetOutput() 1197 1198 rfft = vtk.new("ImageRFFT") 1199 rfft.SetInputData(out) 1200 rfft.Update() 1201 1202 ecomp = vtk.new("ImageExtractComponents") 1203 ecomp.SetInputData(rfft.GetOutput()) 1204 ecomp.SetComponents(0) 1205 ecomp.Update() 1206 self._update(ecomp.GetOutput()) 1207 self.pipeline = utils.OperationNode("frequency_pass_filter", parents=[self], c="#4cc9f0") 1208 return self 1209 1210 def smooth_gaussian(self, sigma=(2, 2, 2), radius=None): 1211 """ 1212 Performs a convolution of the input Volume with a gaussian. 1213 1214 Arguments: 1215 sigma : (float, list) 1216 standard deviation(s) in voxel units. 1217 A list can be given to smooth in the three direction differently. 1218 radius : (float, list) 1219 radius factor(s) determine how far out the gaussian 1220 kernel will go before being clamped to zero. A list can be given too. 1221 """ 1222 gsf = vtk.new("ImageGaussianSmooth") 1223 gsf.SetDimensionality(3) 1224 gsf.SetInputData(self.dataset) 1225 if utils.is_sequence(sigma): 1226 gsf.SetStandardDeviations(sigma) 1227 else: 1228 gsf.SetStandardDeviation(sigma) 1229 if radius is not None: 1230 if utils.is_sequence(radius): 1231 gsf.SetRadiusFactors(radius) 1232 else: 1233 gsf.SetRadiusFactor(radius) 1234 gsf.Update() 1235 self._update(gsf.GetOutput()) 1236 self.pipeline = utils.OperationNode("smooth_gaussian", parents=[self], c="#4cc9f0") 1237 return self 1238 1239 def smooth_median(self, neighbours=(2, 2, 2)): 1240 """ 1241 Median filter that replaces each pixel with the median value 1242 from a rectangular neighborhood around that pixel. 1243 """ 1244 imgm = vtk.new("ImageMedian3D") 1245 imgm.SetInputData(self.dataset) 1246 if utils.is_sequence(neighbours): 1247 imgm.SetKernelSize(neighbours[0], neighbours[1], neighbours[2]) 1248 else: 1249 imgm.SetKernelSize(neighbours, neighbours, neighbours) 1250 imgm.Update() 1251 self._update(imgm.GetOutput()) 1252 self.pipeline = utils.OperationNode("smooth_median", parents=[self], c="#4cc9f0") 1253 return self 1254 1255 def erode(self, neighbours=(2, 2, 2)): 1256 """ 1257 Replace a voxel with the minimum over an ellipsoidal neighborhood of voxels. 1258 If `neighbours` of an axis is 1, no processing is done on that axis. 1259 1260 Examples: 1261 - [erode_dilate.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/erode_dilate.py) 1262 1263  1264 """ 1265 ver = vtk.new("ImageContinuousErode3D") 1266 ver.SetInputData(self.dataset) 1267 ver.SetKernelSize(neighbours[0], neighbours[1], neighbours[2]) 1268 ver.Update() 1269 self._update(ver.GetOutput()) 1270 self.pipeline = utils.OperationNode("erode", parents=[self], c="#4cc9f0") 1271 return self 1272 1273 def dilate(self, neighbours=(2, 2, 2)): 1274 """ 1275 Replace a voxel with the maximum over an ellipsoidal neighborhood of voxels. 1276 If `neighbours` of an axis is 1, no processing is done on that axis. 1277 1278 Check also `erode()` and `pad()`. 1279 1280 Examples: 1281 - [erode_dilate.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/erode_dilate.py) 1282 """ 1283 ver = vtk.new("ImageContinuousDilate3D") 1284 ver.SetInputData(self.dataset) 1285 ver.SetKernelSize(neighbours[0], neighbours[1], neighbours[2]) 1286 ver.Update() 1287 self._update(ver.GetOutput()) 1288 self.pipeline = utils.OperationNode("dilate", parents=[self], c="#4cc9f0") 1289 return self 1290 1291 def magnitude(self): 1292 """Colapses components with magnitude function.""" 1293 imgm = vtk.new("ImageMagnitude") 1294 imgm.SetInputData(self.dataset) 1295 imgm.Update() 1296 self._update(imgm.GetOutput()) 1297 self.pipeline = utils.OperationNode("magnitude", parents=[self], c="#4cc9f0") 1298 return self 1299 1300 def topoints(self): 1301 """ 1302 Extract all image voxels as points. 1303 This function takes an input `Volume` and creates an `Mesh` 1304 that contains the points and the point attributes. 1305 1306 Examples: 1307 - [vol2points.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/vol2points.py) 1308 """ 1309 v2p = vtk.new("ImageToPoints") 1310 v2p.SetInputData(self.dataset) 1311 v2p.Update() 1312 mpts = vedo.Points(v2p.GetOutput()) 1313 mpts.pipeline = utils.OperationNode("topoints", parents=[self], c="#4cc9f0:#e9c46a") 1314 return mpts 1315 1316 def euclidean_distance(self, anisotropy=False, max_distance=None): 1317 """ 1318 Implementation of the Euclidean DT (Distance Transform) using Saito's algorithm. 1319 The distance map produced contains the square of the Euclidean distance values. 1320 The algorithm has a O(n^(D+1)) complexity over n x n x...x n images in D dimensions. 1321 1322 Check out also: https://en.wikipedia.org/wiki/Distance_transform 1323 1324 Arguments: 1325 anisotropy : bool 1326 used to define whether Spacing should be used in the 1327 computation of the distances. 1328 max_distance : bool 1329 any distance bigger than max_distance will not be 1330 computed but set to this specified value instead. 1331 1332 Examples: 1333 - [euclidian_dist.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/euclidian_dist.py) 1334 """ 1335 euv = vtk.new("ImageEuclideanDistance") 1336 euv.SetInputData(self.dataset) 1337 euv.SetConsiderAnisotropy(anisotropy) 1338 if max_distance is not None: 1339 euv.InitializeOn() 1340 euv.SetMaximumDistance(max_distance) 1341 euv.SetAlgorithmToSaito() 1342 euv.Update() 1343 vol = Volume(euv.GetOutput()) 1344 vol.pipeline = utils.OperationNode("euclidean_distance", parents=[self], c="#4cc9f0") 1345 return vol 1346 1347 def correlation_with(self, vol2, dim=2): 1348 """ 1349 Find the correlation between two volumetric data sets. 1350 Keyword `dim` determines whether the correlation will be 3D, 2D or 1D. 1351 The default is a 2D Correlation. 1352 1353 The output size will match the size of the first input. 1354 The second input is considered the correlation kernel. 1355 """ 1356 imc = vtk.new("ImageCorrelation") 1357 imc.SetInput1Data(self.dataset) 1358 imc.SetInput2Data(vol2.dataset) 1359 imc.SetDimensionality(dim) 1360 imc.Update() 1361 vol = Volume(imc.GetOutput()) 1362 1363 vol.pipeline = utils.OperationNode("correlation_with", parents=[self, vol2], c="#4cc9f0") 1364 return vol 1365 1366 def scale_voxels(self, scale=1): 1367 """Scale the voxel content by factor `scale`.""" 1368 rsl = vtk.new("ImageReslice") 1369 rsl.SetInputData(self.dataset) 1370 rsl.SetScalarScale(scale) 1371 rsl.Update() 1372 self._update(rsl.GetOutput()) 1373 self.pipeline = utils.OperationNode( 1374 "scale_voxels", comment=f"scale={scale}", parents=[self], c="#4cc9f0" 1375 ) 1376 return self
30class Volume(VolumeVisual, VolumeAlgorithms): 31 """ 32 Class to describe dataset that are defined on "voxels", 33 the 3D equivalent of 2D pixels. 34 """ 35 def __init__( 36 self, 37 inputobj=None, 38 dims=None, 39 origin=None, 40 spacing=None, 41 ): 42 """ 43 This class can be initialized with a numpy object, 44 a `vtkImageData` or a list of 2D bmp files. 45 46 Arguments: 47 origin : (list) 48 set volume origin coordinates 49 spacing : (list) 50 voxel dimensions in x, y and z. 51 dims : (list) 52 specify the dimensions of the volume. 53 54 Example: 55 ```python 56 from vedo import Volume 57 vol = Volume("path/to/mydata/rec*.bmp") 58 vol.show() 59 ``` 60 61 Examples: 62 - [numpy2volume1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/numpy2volume1.py) 63 64  65 66 - [read_volume2.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/read_volume2.py) 67 68  69 70 .. note:: 71 if a `list` of values is used for `alphas` this is interpreted 72 as a transfer function along the range of the scalar. 73 """ 74 self.name = "Volume" 75 self.filename = "" 76 self.info = {} 77 self.time = 0 78 self.rendered_at = set() 79 80 self.actor = vtk.vtkVolume() 81 self.actor.retrieve_object = weak_ref_to(self) 82 self.properties = self.actor.GetProperty() 83 84 ################### 85 if isinstance(inputobj, str): 86 if "https://" in inputobj: 87 inputobj = vedo.file_io.download(inputobj, verbose=False) # fpath 88 elif os.path.isfile(inputobj): 89 self.filename = inputobj 90 else: 91 inputobj = sorted(glob.glob(inputobj)) 92 93 ################### 94 inputtype = str(type(inputobj)) 95 96 # print('Volume inputtype', inputtype, c='b') 97 98 if inputobj is None: 99 img = vtk.vtkImageData() 100 101 elif utils.is_sequence(inputobj): 102 103 if isinstance(inputobj[0], str) and ".bmp" in inputobj[0].lower(): 104 # scan sequence of BMP files 105 ima = vtk.new("ImageAppend") 106 ima.SetAppendAxis(2) 107 pb = utils.ProgressBar(0, len(inputobj)) 108 for i in pb.range(): 109 f = inputobj[i] 110 if "_rec_spr" in f: # OPT specific 111 continue 112 picr = vtk.new("BMPReader") 113 picr.SetFileName(f) 114 picr.Update() 115 mgf = vtk.new("ImageMagnitude") 116 mgf.SetInputData(picr.GetOutput()) 117 mgf.Update() 118 ima.AddInputData(mgf.GetOutput()) 119 pb.print("loading...") 120 ima.Update() 121 img = ima.GetOutput() 122 123 else: 124 125 if len(inputobj.shape) == 1: 126 varr = utils.numpy2vtk(inputobj) 127 else: 128 varr = utils.numpy2vtk(inputobj.ravel(order="F")) 129 varr.SetName("input_scalars") 130 131 img = vtk.vtkImageData() 132 if dims is not None: 133 img.SetDimensions(dims[2], dims[1], dims[0]) 134 else: 135 if len(inputobj.shape) == 1: 136 vedo.logger.error("must set dimensions (dims keyword) in Volume") 137 raise RuntimeError() 138 img.SetDimensions(inputobj.shape) 139 img.GetPointData().AddArray(varr) 140 img.GetPointData().SetActiveScalars(varr.GetName()) 141 142 elif isinstance(inputobj, vtk.vtkImageData): 143 img = inputobj 144 145 elif isinstance(inputobj, str): 146 if "https://" in inputobj: 147 inputobj = vedo.file_io.download(inputobj, verbose=False) 148 img = vedo.file_io.loadImageData(inputobj) 149 self.filename = inputobj 150 151 else: 152 vedo.logger.error(f"cannot understand input type {inputtype}") 153 return 154 155 if dims is not None: 156 img.SetDimensions(dims) 157 158 if origin is not None: 159 img.SetOrigin(origin) 160 161 if spacing is not None: 162 img.SetSpacing(spacing) 163 164 self.dataset = img 165 self.transform = None 166 167 ##################################### 168 mapper = vtk.new("SmartVolumeMapper") 169 mapper.SetInputData(img) 170 self.actor.SetMapper(mapper) 171 172 if img.GetPointData().GetScalars(): 173 if img.GetPointData().GetScalars().GetNumberOfComponents() == 1: 174 self.properties.SetShade(True) 175 self.properties.SetInterpolationType(1) 176 self.cmap("RdBu_r") 177 self.alpha([0.0, 0.0, 0.2, 0.4, 0.8, 1.0]) 178 self.alpha_gradient(None) 179 self.properties.SetScalarOpacityUnitDistance(1.0) 180 181 self.pipeline = utils.OperationNode( 182 "Volume", comment=f"dims={tuple(self.dimensions())}", c="#4cc9f0" 183 ) 184 ####################################################################### 185 186 @property 187 def mapper(self): 188 """Return the underlying `vtkMapper` object.""" 189 return self.actor.GetMapper() 190 191 @mapper.setter 192 def mapper(self, mapper): 193 """ 194 Set the underlying `vtkMapper` object. 195 196 Arguments: 197 mapper : (str, vtkMapper) 198 either 'gpu', 'opengl_gpu', 'fixed' or 'smart' 199 """ 200 if isinstance(mapper, 201 (vtk.get_class("Mapper"), 202 vtk.get_class("ImageResliceMapper", 203 ) )): 204 pass 205 elif "gpu" in mapper: 206 mapper = vtk.new("GPUVolumeRayCastMapper") 207 elif "opengl_gpu" in mapper: 208 mapper = vtk.new("OpenGLGPUVolumeRayCastMapper") 209 elif "smart" in mapper: 210 mapper = vtk.new("SmartVolumeMapper") 211 elif "fixed" in mapper: 212 mapper = vtk.new("FixedPointVolumeRayCastMapper") 213 else: 214 print("Error unknown mapper type", [mapper]) 215 raise RuntimeError() 216 self.actor.SetMapper(mapper) 217 218 def c(self, *args, **kwargs): 219 """Deprecated. Use `Volume.cmap()` instead.""" 220 vedo.logger.warning("Volume.c() is deprecated, use Volume.cmap() instead") 221 return self.cmap(*args, **kwargs) 222 223 def _update(self, data): 224 self.dataset = data 225 self.mapper.SetInputData(data) 226 self.dataset.GetPointData().Modified() 227 self.mapper.Modified() 228 self.mapper.Update() 229 return self 230 231 def __str__(self): 232 """Print a summary for the Volume object.""" 233 module = self.__class__.__module__ 234 name = self.__class__.__name__ 235 out = vedo.printc( 236 f"{module}.{name} at ({hex(self.memory_address())})".ljust(75), 237 c="c", bold=True, invert=True, return_string=True, 238 ) 239 out += "\x1b[0m\x1b[36;1m" 240 241 out+= "name".ljust(14) + ": " + str(self.name) + "\n" 242 if self.filename: 243 out+= "filename".ljust(14) + ": " + str(self.filename) + "\n" 244 245 out+= "dimensions".ljust(14) + ": " + str(self.shape) + "\n" 246 247 out+= "origin".ljust(14) + ": " 248 out+= utils.precision(self.origin(), 6) + "\n" 249 250 out+= "center".ljust(14) + ": " 251 out+= utils.precision(self.center(), 6) + "\n" 252 253 out+= "spacing".ljust(14) + ": " 254 out+= utils.precision(self.spacing(), 6) + "\n" 255 256 bnds = self.bounds() 257 bx1, bx2 = utils.precision(bnds[0], 3), utils.precision(bnds[1], 3) 258 by1, by2 = utils.precision(bnds[2], 3), utils.precision(bnds[3], 3) 259 bz1, bz2 = utils.precision(bnds[4], 3), utils.precision(bnds[5], 3) 260 out+= "bounds".ljust(14) + ":" 261 out+= " x=(" + bx1 + ", " + bx2 + ")," 262 out+= " y=(" + by1 + ", " + by2 + ")," 263 out+= " z=(" + bz1 + ", " + bz2 + ")\n" 264 265 out+= "memory size".ljust(14) + ": " 266 out+= str(int(self.dataset.GetActualMemorySize()/1024+0.5))+" MB\n" 267 268 st = self.dataset.GetScalarTypeAsString() 269 out+= "scalar size".ljust(14) + ": " 270 out+= str(self.dataset.GetScalarSize()) + f" bytes ({st})\n" 271 out+= "scalar range".ljust(14) + ": " 272 out+= str(self.dataset.GetScalarRange()) + "\n" 273 274 #utils.print_histogram(self, logscale=True, bins=8, height=15, c="b", bold=True) 275 return out.rstrip() + "\x1b[0m" 276 277 def _repr_html_(self): 278 """ 279 HTML representation of the Volume object for Jupyter Notebooks. 280 281 Returns: 282 HTML text with the image and some properties. 283 """ 284 import io 285 import base64 286 from PIL import Image 287 288 library_name = "vedo.volume.Volume" 289 help_url = "https://vedo.embl.es/docs/vedo/volume.html" 290 291 arr = self.thumbnail(azimuth=0, elevation=-60, zoom=1.4, axes=True) 292 293 im = Image.fromarray(arr) 294 buffered = io.BytesIO() 295 im.save(buffered, format="PNG", quality=100) 296 encoded = base64.b64encode(buffered.getvalue()).decode("utf-8") 297 url = "data:image/png;base64," + encoded 298 image = f"<img src='{url}'></img>" 299 300 # statisitics 301 bounds = "<br/>".join( 302 [ 303 utils.precision(min_x, 4) + " ... " + utils.precision(max_x, 4) 304 for min_x, max_x in zip(self.bounds()[::2], self.bounds()[1::2]) 305 ] 306 ) 307 308 help_text = "" 309 if self.name: 310 help_text += f"<b> {self.name}:   </b>" 311 help_text += '<b><a href="' + help_url + '" target="_blank">' + library_name + "</a></b>" 312 if self.filename: 313 dots = "" 314 if len(self.filename) > 30: 315 dots = "..." 316 help_text += f"<br/><code><i>({dots}{self.filename[-30:]})</i></code>" 317 318 pdata = "" 319 if self.dataset.GetPointData().GetScalars(): 320 if self.dataset.GetPointData().GetScalars().GetName(): 321 name = self.dataset.GetPointData().GetScalars().GetName() 322 pdata = "<tr><td><b> point data array </b></td><td>" + name + "</td></tr>" 323 324 cdata = "" 325 if self.dataset.GetCellData().GetScalars(): 326 if self.dataset.GetCellData().GetScalars().GetName(): 327 name = self.dataset.GetCellData().GetScalars().GetName() 328 cdata = "<tr><td><b> voxel data array </b></td><td>" + name + "</td></tr>" 329 330 img = self.dataset 331 332 allt = [ 333 "<table>", 334 "<tr>", 335 "<td>", 336 image, 337 "</td>", 338 "<td style='text-align: center; vertical-align: center;'><br/>", 339 help_text, 340 "<table>", 341 "<tr><td><b> bounds </b> <br/> (x/y/z) </td><td>" + str(bounds) + "</td></tr>", 342 "<tr><td><b> dimensions </b></td><td>" + str(img.GetDimensions()) + "</td></tr>", 343 "<tr><td><b> voxel spacing </b></td><td>" 344 + utils.precision(img.GetSpacing(), 3) 345 + "</td></tr>", 346 "<tr><td><b> in memory size </b></td><td>" 347 + str(int(img.GetActualMemorySize() / 1024)) 348 + "MB</td></tr>", 349 pdata, 350 cdata, 351 "<tr><td><b> scalar range </b></td><td>" 352 + utils.precision(img.GetScalarRange(), 4) 353 + "</td></tr>", 354 "</table>", 355 "</table>", 356 ] 357 return "\n".join(allt) 358 359 def copy(self, deep=True): 360 """Return a copy of the Volume. Alias of `clone()`.""" 361 return self.clone(deep=deep) 362 363 def clone(self, deep=True): 364 """Return a clone copy of the Volume. Alias of `copy()`.""" 365 if deep: 366 newimg = vtk.vtkImageData() 367 newimg.CopyStructure(self.dataset) 368 newimg.CopyAttributes(self.dataset) 369 newvol = Volume(newimg) 370 else: 371 newvol = Volume(self.dataset) 372 373 prop = vtk.vtkVolumeProperty() 374 prop.DeepCopy(self.properties) 375 newvol.actor.SetProperty(prop) 376 newvol.properties = prop 377 378 newvol.pipeline = utils.OperationNode("clone", parents=[self], c="#bbd0ff", shape="diamond") 379 return newvol 380 381 def component_weight(self, i, weight): 382 """Set the scalar component weight in range [0,1].""" 383 self.properties.SetComponentWeight(i, weight) 384 return self 385 386 def xslice(self, i): 387 """Extract the slice at index `i` of volume along x-axis.""" 388 vslice = vtk.new("ImageDataGeometryFilter") 389 vslice.SetInputData(self.dataset) 390 nx, ny, nz = self.dataset.GetDimensions() 391 if i > nx - 1: 392 i = nx - 1 393 vslice.SetExtent(i, i, 0, ny, 0, nz) 394 vslice.Update() 395 m = Mesh(vslice.GetOutput()) 396 m.pipeline = utils.OperationNode(f"xslice {i}", parents=[self], c="#4cc9f0:#e9c46a") 397 return m 398 399 def yslice(self, j): 400 """Extract the slice at index `j` of volume along y-axis.""" 401 vslice = vtk.new("ImageDataGeometryFilter") 402 vslice.SetInputData(self.dataset) 403 nx, ny, nz = self.dataset.GetDimensions() 404 if j > ny - 1: 405 j = ny - 1 406 vslice.SetExtent(0, nx, j, j, 0, nz) 407 vslice.Update() 408 m = Mesh(vslice.GetOutput()) 409 m.pipeline = utils.OperationNode(f"yslice {j}", parents=[self], c="#4cc9f0:#e9c46a") 410 return m 411 412 def zslice(self, k): 413 """Extract the slice at index `i` of volume along z-axis.""" 414 vslice = vtk.new("ImageDataGeometryFilter") 415 vslice.SetInputData(self.dataset) 416 nx, ny, nz = self.dataset.GetDimensions() 417 if k > nz - 1: 418 k = nz - 1 419 vslice.SetExtent(0, nx, 0, ny, k, k) 420 vslice.Update() 421 m = Mesh(vslice.GetOutput()) 422 m.pipeline = utils.OperationNode(f"zslice {k}", parents=[self], c="#4cc9f0:#e9c46a") 423 return m 424 425 def slice_plane(self, origin=(0, 0, 0), normal=(1, 1, 1), autocrop=False): 426 """ 427 Extract the slice along a given plane position and normal. 428 429 Example: 430 - [slice_plane1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/slice_plane1.py) 431 432  433 """ 434 reslice = vtk.new("ImageReslice") 435 reslice.SetInputData(self.dataset) 436 reslice.SetOutputDimensionality(2) 437 newaxis = utils.versor(normal) 438 pos = np.array(origin) 439 initaxis = (0, 0, 1) 440 crossvec = np.cross(initaxis, newaxis) 441 angle = np.arccos(np.dot(initaxis, newaxis)) 442 T = vtk.vtkTransform() 443 T.PostMultiply() 444 T.RotateWXYZ(np.rad2deg(angle), crossvec) 445 T.Translate(pos) 446 M = T.GetMatrix() 447 reslice.SetResliceAxes(M) 448 reslice.SetInterpolationModeToLinear() 449 reslice.SetAutoCropOutput(not autocrop) 450 reslice.Update() 451 vslice = vtk.new("ImageDataGeometryFilter") 452 vslice.SetInputData(reslice.GetOutput()) 453 vslice.Update() 454 msh = Mesh(vslice.GetOutput()) 455 msh.apply_transform(T) 456 msh.pipeline = utils.OperationNode( 457 "slice_plane", parents=[self], c="#4cc9f0:#e9c46a") 458 return msh 459 460 def slab(self, slice_range=(), axis='z', operation="mean"): 461 """ 462 Extract a slab from a `Volume` by combining 463 all of the slices of an image to create a single slice. 464 465 Returns a `Mesh` containing metadata which 466 can be accessed with e.g. `mesh.metadata["slab_range"]`. 467 468 Metadata: 469 slab_range : (list) 470 contains the range of slices extracted 471 slab_axis : (str) 472 contains the axis along which the slab was extracted 473 slab_operation : (str) 474 contains the operation performed on the slab 475 slab_bounding_box : (list) 476 contains the bounding box of the slab 477 478 Arguments: 479 slice_range : (list) 480 range of slices to extract 481 axis : (str) 482 axis along which to extract the slab 483 operation : (str) 484 operation to perform on the slab, 485 allowed values are: "sum", "min", "max", "mean". 486 """ 487 if len(slice_range) != 2: 488 vedo.logger.error("in slab(): slice_range is empty or invalid") 489 raise ValueError() 490 491 islab = vtk.new("ImageSlab") 492 islab.SetInputData(self.dataset) 493 494 if operation in ["+", "add", "sum"]: 495 islab.SetOperationToSum() 496 elif "min" in operation: 497 islab.SetOperationToMin() 498 elif "max" in operation: 499 islab.SetOperationToMax() 500 elif "mean" in operation: 501 islab.SetOperationToMean() 502 else: 503 vedo.logger.error(f"in slab(): unknown operation {operation}") 504 raise ValueError() 505 506 dims = self.dimensions() 507 if axis == 'x': 508 islab.SetOrientationToX() 509 if slice_range[0] > dims[0]-1: 510 slice_range[0] = dims[0]-1 511 if slice_range[1] > dims[0]-1: 512 slice_range[1] = dims[0]-1 513 elif axis == 'y': 514 islab.SetOrientationToY() 515 if slice_range[0] > dims[1]-1: 516 slice_range[0] = dims[1]-1 517 if slice_range[1] > dims[1]-1: 518 slice_range[1] = dims[1]-1 519 elif axis == 'z': 520 islab.SetOrientationToZ() 521 if slice_range[0] > dims[2]-1: 522 slice_range[0] = dims[2]-1 523 if slice_range[1] > dims[2]-1: 524 slice_range[1] = dims[2]-1 525 else: 526 vedo.logger.error(f"Error in slab(): unknown axis {axis}") 527 raise RuntimeError() 528 529 islab.SetSliceRange(slice_range) 530 islab.Update() 531 532 msh = Mesh(islab.GetOutput()).lighting('off') 533 msh.mapper.SetLookupTable(utils.ctf2lut(self, msh)) 534 msh.mapper.SetScalarRange(self.scalar_range()) 535 536 msh.metadata["slab_range"] = slice_range 537 msh.metadata["slab_axis"] = axis 538 msh.metadata["slab_operation"] = operation 539 540 # compute bounds of slab 541 origin = self.origin() 542 spacing = self.spacing() 543 if axis == 'x': 544 msh.metadata["slab_bounding_box"] = [ 545 origin[0] + slice_range[0]*spacing[0], 546 origin[0] + slice_range[1]*spacing[0], 547 origin[1], 548 origin[1] + dims[1]*spacing[1], 549 origin[2], 550 origin[2] + dims[2]*spacing[2], 551 ] 552 elif axis == 'y': 553 msh.metadata["slab_bounding_box"] = [ 554 origin[0], 555 origin[0] + dims[0]*spacing[0], 556 origin[1] + slice_range[0]*spacing[1], 557 origin[1] + slice_range[1]*spacing[1], 558 origin[2], 559 origin[2] + dims[2]*spacing[2], 560 ] 561 elif axis == 'z': 562 msh.metadata["slab_bounding_box"] = [ 563 origin[0], 564 origin[0] + dims[0]*spacing[0], 565 origin[1], 566 origin[1] + dims[1]*spacing[1], 567 origin[2] + slice_range[0]*spacing[2], 568 origin[2] + slice_range[1]*spacing[2], 569 ] 570 571 msh.pipeline = utils.OperationNode( 572 f"slab{slice_range}", 573 comment=f"axis={axis}, operation={operation}", 574 parents=[self], 575 c="#4cc9f0:#e9c46a", 576 ) 577 msh.name = "SlabMesh" 578 return msh 579 580 581 def warp(self, source, target, sigma=1, mode="3d", fit=True): 582 """ 583 Warp volume scalars within a Volume by specifying 584 source and target sets of points. 585 586 Arguments: 587 source : (Points, list) 588 the list of source points 589 target : (Points, list) 590 the list of target points 591 fit : (bool) 592 fit/adapt the old bounding box to the warped geometry 593 """ 594 if isinstance(source, vedo.Points): 595 source = source.vertices 596 if isinstance(target, vedo.Points): 597 target = target.vertices 598 599 NLT = transformations.NonLinearTransform() 600 NLT.source_points = source 601 NLT.target_points = target 602 NLT.sigma = sigma 603 NLT.mode = mode 604 605 self.apply_transform(NLT, fit=fit) 606 self.pipeline = utils.OperationNode("warp", parents=[self], c="#4cc9f0") 607 return self 608 609 def apply_transform(self, T, fit=True, interpolation="cubic"): 610 """ 611 Apply a transform to the scalars in the volume. 612 613 Arguments: 614 T : (LinearTransform, NonLinearTransform) 615 The transformation to be applied 616 fit : (bool) 617 fit/adapt the old bounding box to the modified geometry 618 interpolation : (str) 619 one of the following: "linear", "nearest", "cubic" 620 """ 621 TI = T.compute_inverse() 622 reslice = vtk.new("ImageReslice") 623 reslice.SetInputData(self.dataset) 624 reslice.SetResliceTransform(TI.T) 625 reslice.SetOutputDimensionality(3) 626 if "lin" in interpolation.lower(): 627 reslice.SetInterpolationModeToLinear() 628 elif "near" in interpolation.lower(): 629 reslice.SetInterpolationModeToNearestNeighbor() 630 elif "cubic" in interpolation.lower(): 631 reslice.SetInterpolationModeToCubic() 632 else: 633 vedo.logger.error( 634 f"in apply_transform: unknown interpolation mode {interpolation}") 635 raise ValueError() 636 reslice.SetAutoCropOutput(fit) 637 reslice.Update() 638 self._update(reslice.GetOutput()) 639 self.transform = T 640 self.pipeline = utils.OperationNode( 641 "apply_transform", parents=[self], c="#4cc9f0") 642 return self 643 644 def imagedata(self): 645 """ 646 DEPRECATED: 647 Use `Volume.dataset` instead. 648 649 Return the underlying `vtkImagaData` object. 650 """ 651 print("Volume.imagedata() is deprecated, use Volume.dataset instead") 652 return self.dataset 653 654 def modified(self): 655 """ 656 Mark the object as modified. 657 658 Example: 659 660 - [numpy2volume0.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/numpy2volume0.py) 661 """ 662 self.dataset.GetPointData().GetScalars().Modified() 663 return self 664 665 def tonumpy(self): 666 """ 667 Get read-write access to voxels of a Volume object as a numpy array. 668 669 When you set values in the output image, you don't want numpy to reallocate the array 670 but instead set values in the existing array, so use the [:] operator. 671 672 Example: 673 `arr[:] = arr*2 + 15` 674 675 If the array is modified add a call to: 676 `volume.modified()` 677 when all your modifications are completed. 678 """ 679 narray_shape = tuple(reversed(self.dataset.GetDimensions())) 680 681 scals = self.dataset.GetPointData().GetScalars() 682 comps = scals.GetNumberOfComponents() 683 if comps == 1: 684 narray = utils.vtk2numpy(scals).reshape(narray_shape) 685 narray = np.transpose(narray, axes=[2, 1, 0]) 686 else: 687 narray = utils.vtk2numpy(scals).reshape(*narray_shape, comps) 688 narray = np.transpose(narray, axes=[2, 1, 0, 3]) 689 690 # narray = utils.vtk2numpy(self.dataset.GetPointData().GetScalars()).reshape(narray_shape) 691 # narray = np.transpose(narray, axes=[2, 1, 0]) 692 693 return narray 694 695 @property 696 def shape(self): 697 """Return the nr. of voxels in the 3 dimensions.""" 698 return np.array(self.dataset.GetDimensions()) 699 700 def dimensions(self): 701 """Return the nr. of voxels in the 3 dimensions.""" 702 return np.array(self.dataset.GetDimensions()) 703 704 def scalar_range(self): 705 """Return the range of the scalar values.""" 706 return np.array(self.dataset.GetScalarRange()) 707 708 def spacing(self, s=None): 709 """Set/get the voxels size in the 3 dimensions.""" 710 if s is not None: 711 self.dataset.SetSpacing(s) 712 return self 713 return np.array(self.dataset.GetSpacing()) 714 715 def origin(self, s=None): 716 """ 717 Set/get the origin of the volumetric dataset. 718 719 The origin is the position in world coordinates of the point index (0,0,0). 720 This point does not have to be part of the dataset, in other words, 721 the dataset extent does not have to start at (0,0,0) and the origin 722 can be outside of the dataset bounding box. 723 The origin plus spacing determine the position in space of the points. 724 """ 725 if s is not None: 726 self.dataset.SetOrigin(s) 727 return self 728 return np.array(self.dataset.GetOrigin()) 729 730 def center(self): 731 """Get the center of the volumetric dataset.""" 732 # note that this does not have the set method like origin and spacing 733 return np.array(self.dataset.GetCenter()) 734 735 def get_cell_from_ijk(self, ijk): 736 """ 737 Get the voxel id number at the given ijk coordinates. 738 739 Arguments: 740 ijk : (list) 741 the ijk coordinates of the voxel 742 """ 743 return self.ComputeCellId(ijk) 744 745 def get_point_from_ijk(self, ijk): 746 """ 747 Get the point id number at the given ijk coordinates. 748 749 Arguments: 750 ijk : (list) 751 the ijk coordinates of the voxel 752 """ 753 return self.ComputePointId(ijk) 754 755 def permute_axes(self, x, y, z): 756 """ 757 Reorder the axes of the Volume by specifying 758 the input axes which are supposed to become the new X, Y, and Z. 759 """ 760 imp = vtk.new("ImagePermute") 761 imp.SetFilteredAxes(x, y, z) 762 imp.SetInputData(self.dataset) 763 imp.Update() 764 self._update(imp.GetOutput()) 765 self.pipeline = utils.OperationNode( 766 f"permute_axes({(x,y,z)})", parents=[self], c="#4cc9f0" 767 ) 768 return self 769 770 def resample(self, new_spacing, interpolation=1): 771 """ 772 Resamples a `Volume` to be larger or smaller. 773 774 This method modifies the spacing of the input. 775 Linear interpolation is used to resample the data. 776 777 Arguments: 778 new_spacing : (list) 779 a list of 3 new spacings for the 3 axes 780 interpolation : (int) 781 0=nearest_neighbor, 1=linear, 2=cubic 782 """ 783 rsp = vtk.new("ImageResample") 784 oldsp = self.spacing() 785 for i in range(3): 786 if oldsp[i] != new_spacing[i]: 787 rsp.SetAxisOutputSpacing(i, new_spacing[i]) 788 rsp.InterpolateOn() 789 rsp.SetInterpolationMode(interpolation) 790 rsp.OptimizationOn() 791 rsp.Update() 792 self._update(rsp.GetOutput()) 793 self.pipeline = utils.OperationNode( 794 "resample", comment=f"spacing: {tuple(new_spacing)}", parents=[self], c="#4cc9f0" 795 ) 796 return self 797 798 799 def threshold(self, above=None, below=None, replace=None, replace_value=None): 800 """ 801 Binary or continuous volume thresholding. 802 Find the voxels that contain a value above/below the input values 803 and replace them with a new value (default is 0). 804 """ 805 th = vtk.new("ImageThreshold") 806 th.SetInputData(self.dataset) 807 808 # sanity checks 809 if above is not None and below is not None: 810 if above == below: 811 return self 812 if above > below: 813 vedo.logger.warning("in volume.threshold(), above > below, skip.") 814 return self 815 816 ## cases 817 if below is not None and above is not None: 818 th.ThresholdBetween(above, below) 819 820 elif above is not None: 821 th.ThresholdByUpper(above) 822 823 elif below is not None: 824 th.ThresholdByLower(below) 825 826 ## 827 if replace is not None: 828 th.SetReplaceIn(True) 829 th.SetInValue(replace) 830 else: 831 th.SetReplaceIn(False) 832 833 if replace_value is not None: 834 th.SetReplaceOut(True) 835 th.SetOutValue(replace_value) 836 else: 837 th.SetReplaceOut(False) 838 839 th.Update() 840 self._update(th.GetOutput()) 841 self.pipeline = utils.OperationNode("threshold", parents=[self], c="#4cc9f0") 842 return self 843 844 def crop(self, left=None, right=None, back=None, front=None, bottom=None, top=None, VOI=()): 845 """ 846 Crop a `Volume` object. 847 848 Arguments: 849 left : (float) 850 fraction to crop from the left plane (negative x) 851 right : (float) 852 fraction to crop from the right plane (positive x) 853 back : (float) 854 fraction to crop from the back plane (negative y) 855 front : (float) 856 fraction to crop from the front plane (positive y) 857 bottom : (float) 858 fraction to crop from the bottom plane (negative z) 859 top : (float) 860 fraction to crop from the top plane (positive z) 861 VOI : (list) 862 extract Volume Of Interest expressed in voxel numbers 863 864 Example: 865 `vol.crop(VOI=(xmin, xmax, ymin, ymax, zmin, zmax)) # all integers nrs` 866 """ 867 extractVOI = vtk.new("ExtractVOI") 868 extractVOI.SetInputData(self.dataset) 869 870 if VOI: 871 extractVOI.SetVOI(VOI) 872 else: 873 d = self.dataset.GetDimensions() 874 bx0, bx1, by0, by1, bz0, bz1 = 0, d[0]-1, 0, d[1]-1, 0, d[2]-1 875 if left is not None: bx0 = int((d[0]-1)*left) 876 if right is not None: bx1 = int((d[0]-1)*(1-right)) 877 if back is not None: by0 = int((d[1]-1)*back) 878 if front is not None: by1 = int((d[1]-1)*(1-front)) 879 if bottom is not None: bz0 = int((d[2]-1)*bottom) 880 if top is not None: bz1 = int((d[2]-1)*(1-top)) 881 extractVOI.SetVOI(bx0, bx1, by0, by1, bz0, bz1) 882 extractVOI.Update() 883 self._update(extractVOI.GetOutput()) 884 885 self.pipeline = utils.OperationNode( 886 "crop", parents=[self], c="#4cc9f0", comment=f"dims={tuple(self.dimensions())}" 887 ) 888 return self 889 890 def append(self, volumes, axis="z", preserve_extents=False): 891 """ 892 Take the components from multiple inputs and merges them into one output. 893 Except for the append axis, all inputs must have the same extent. 894 All inputs must have the same number of scalar components. 895 The output has the same origin and spacing as the first input. 896 The origin and spacing of all other inputs are ignored. 897 All inputs must have the same scalar type. 898 899 Arguments: 900 axis : (int, str) 901 axis expanded to hold the multiple images 902 preserve_extents : (bool) 903 if True, the extent of the inputs is used to place 904 the image in the output. The whole extent of the output is the union of the input 905 whole extents. Any portion of the output not covered by the inputs is set to zero. 906 The origin and spacing is taken from the first input. 907 908 Example: 909 ```python 910 from vedo import Volume, dataurl 911 vol = Volume(dataurl+'embryo.tif') 912 vol.append(vol, axis='x').show().close() 913 ``` 914  915 """ 916 ima = vtk.new("ImageAppend") 917 ima.SetInputData(self.dataset) 918 if not utils.is_sequence(volumes): 919 volumes = [volumes] 920 for volume in volumes: 921 if isinstance(volume, vtk.vtkImageData): 922 ima.AddInputData(volume) 923 else: 924 ima.AddInputData(volume.dataset) 925 ima.SetPreserveExtents(preserve_extents) 926 if axis == "x": 927 axis = 0 928 elif axis == "y": 929 axis = 1 930 elif axis == "z": 931 axis = 2 932 ima.SetAppendAxis(axis) 933 ima.Update() 934 self._update(ima.GetOutput()) 935 936 self.pipeline = utils.OperationNode( 937 "append", 938 parents=[self, *volumes], 939 c="#4cc9f0", 940 comment=f"dims={tuple(self.dimensions())}", 941 ) 942 return self 943 944 def pad(self, voxels=10, value=0): 945 """ 946 Add the specified number of voxels at the `Volume` borders. 947 Voxels can be a list formatted as `[nx0, nx1, ny0, ny1, nz0, nz1]`. 948 949 Arguments: 950 voxels : (int, list) 951 number of voxels to be added (or a list of length 4) 952 value : (int) 953 intensity value (gray-scale color) of the padding 954 955 Example: 956 ```python 957 from vedo import Volume, dataurl, show 958 iso = Volume(dataurl+'embryo.tif').isosurface() 959 vol = iso.binarize(spacing=(100, 100, 100)).pad(10) 960 vol.dilate([15,15,15]) 961 show(iso, vol.isosurface(), N=2, axes=1) 962 ``` 963  964 """ 965 x0, x1, y0, y1, z0, z1 = self.dataset.GetExtent() 966 pf = vtk.new("ImageConstantPad") 967 pf.SetInputData(self.dataset) 968 pf.SetConstant(value) 969 if utils.is_sequence(voxels): 970 pf.SetOutputWholeExtent( 971 x0 - voxels[0], x1 + voxels[1], 972 y0 - voxels[2], y1 + voxels[3], 973 z0 - voxels[4], z1 + voxels[5], 974 ) 975 else: 976 pf.SetOutputWholeExtent( 977 x0 - voxels, x1 + voxels, 978 y0 - voxels, y1 + voxels, 979 z0 - voxels, z1 + voxels, 980 ) 981 pf.Update() 982 self._update(pf.GetOutput()) 983 self.pipeline = utils.OperationNode( 984 "pad", comment=f"{voxels} voxels", parents=[self], c="#f28482" 985 ) 986 return self 987 988 def resize(self, newdims): 989 """Increase or reduce the number of voxels of a Volume with interpolation.""" 990 rsz = vtk.new("ImageResize") 991 rsz.SetResizeMethodToOutputDimensions() 992 rsz.SetInputData(self.dataset) 993 rsz.SetOutputDimensions(newdims) 994 rsz.Update() 995 self.dataset = rsz.GetOutput() 996 self._update(self.dataset) 997 self.pipeline = utils.OperationNode( 998 "resize", parents=[self], c="#4cc9f0", comment=f"dims={tuple(self.dimensions())}" 999 ) 1000 return self 1001 1002 def normalize(self): 1003 """Normalize that scalar components for each point.""" 1004 norm = vtk.new("ImageNormalize") 1005 norm.SetInputData(self.dataset) 1006 norm.Update() 1007 self._update(norm.GetOutput()) 1008 self.pipeline = utils.OperationNode("normalize", parents=[self], c="#4cc9f0") 1009 return self 1010 1011 def mirror(self, axis="x"): 1012 """ 1013 Mirror flip along one of the cartesian axes. 1014 """ 1015 img = self.dataset 1016 1017 ff = vtk.new("ImageFlip") 1018 ff.SetInputData(img) 1019 if axis.lower() == "x": 1020 ff.SetFilteredAxis(0) 1021 elif axis.lower() == "y": 1022 ff.SetFilteredAxis(1) 1023 elif axis.lower() == "z": 1024 ff.SetFilteredAxis(2) 1025 else: 1026 vedo.logger.error("mirror must be set to either x, y, z or n") 1027 raise RuntimeError() 1028 ff.Update() 1029 self._update(ff.GetOutput()) 1030 self.pipeline = utils.OperationNode(f"mirror {axis}", parents=[self], c="#4cc9f0") 1031 return self 1032 1033 def operation(self, operation, volume2=None): 1034 """ 1035 Perform operations with `Volume` objects. 1036 Keyword `volume2` can be a constant float. 1037 1038 Possible operations are: 1039 ``` 1040 +, -, /, 1/x, sin, cos, exp, log, 1041 abs, **2, sqrt, min, max, atan, atan2, median, 1042 mag, dot, gradient, divergence, laplacian. 1043 ``` 1044 1045 Examples: 1046 - [volume_operations.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/volume_operations.py) 1047 """ 1048 op = operation.lower() 1049 image1 = self.dataset 1050 1051 mf = None 1052 if op in ["median"]: 1053 mf = vtk.new("ImageMedian3D") 1054 mf.SetInputData(image1) 1055 elif op in ["mag"]: 1056 mf = vtk.new("ImageMagnitude") 1057 mf.SetInputData(image1) 1058 elif op in ["dot", "dotproduct"]: 1059 mf = vtk.new("ImageDotProduct") 1060 mf.SetInput1Data(image1) 1061 mf.SetInput2Data(volume2.dataset) 1062 elif op in ["grad", "gradient"]: 1063 mf = vtk.new("ImageGradient") 1064 mf.SetDimensionality(3) 1065 mf.SetInputData(image1) 1066 elif op in ["div", "divergence"]: 1067 mf = vtk.new("ImageDivergence") 1068 mf.SetInputData(image1) 1069 elif op in ["laplacian"]: 1070 mf = vtk.new("ImageLaplacian") 1071 mf.SetDimensionality(3) 1072 mf.SetInputData(image1) 1073 1074 if mf is not None: 1075 mf.Update() 1076 vol = Volume(mf.GetOutput()) 1077 vol.pipeline = utils.OperationNode( 1078 "operation", comment=f"{op}", parents=[self], c="#4cc9f0", shape="cylinder" 1079 ) 1080 return vol ########################### 1081 1082 mat = vtk.new("ImageMathematics") 1083 mat.SetInput1Data(image1) 1084 1085 K = None 1086 1087 if utils.is_number(volume2): 1088 K = volume2 1089 mat.SetConstantK(K) 1090 mat.SetConstantC(K) 1091 1092 elif volume2 is not None: # assume image2 is a constant value 1093 mat.SetInput2Data(volume2.dataset) 1094 1095 # ########################### 1096 if op in ["+", "add", "plus"]: 1097 if K: 1098 mat.SetOperationToAddConstant() 1099 else: 1100 mat.SetOperationToAdd() 1101 1102 elif op in ["-", "subtract", "minus"]: 1103 if K: 1104 mat.SetConstantC(-float(K)) 1105 mat.SetOperationToAddConstant() 1106 else: 1107 mat.SetOperationToSubtract() 1108 1109 elif op in ["*", "multiply", "times"]: 1110 if K: 1111 mat.SetOperationToMultiplyByK() 1112 else: 1113 mat.SetOperationToMultiply() 1114 1115 elif op in ["/", "divide"]: 1116 if K: 1117 mat.SetConstantK(1.0 / K) 1118 mat.SetOperationToMultiplyByK() 1119 else: 1120 mat.SetOperationToDivide() 1121 1122 elif op in ["1/x", "invert"]: 1123 mat.SetOperationToInvert() 1124 elif op in ["sin"]: 1125 mat.SetOperationToSin() 1126 elif op in ["cos"]: 1127 mat.SetOperationToCos() 1128 elif op in ["exp"]: 1129 mat.SetOperationToExp() 1130 elif op in ["log"]: 1131 mat.SetOperationToLog() 1132 elif op in ["abs"]: 1133 mat.SetOperationToAbsoluteValue() 1134 elif op in ["**2", "square"]: 1135 mat.SetOperationToSquare() 1136 elif op in ["sqrt", "sqr"]: 1137 mat.SetOperationToSquareRoot() 1138 elif op in ["min"]: 1139 mat.SetOperationToMin() 1140 elif op in ["max"]: 1141 mat.SetOperationToMax() 1142 elif op in ["atan"]: 1143 mat.SetOperationToATAN() 1144 elif op in ["atan2"]: 1145 mat.SetOperationToATAN2() 1146 else: 1147 vedo.logger.error(f"unknown operation {operation}") 1148 raise RuntimeError() 1149 mat.Update() 1150 1151 self._update(mat.GetOutput()) 1152 1153 self.pipeline = utils.OperationNode( 1154 "operation", comment=f"{op}", parents=[self, volume2], shape="cylinder", c="#4cc9f0" 1155 ) 1156 return self 1157 1158 def frequency_pass_filter(self, low_cutoff=None, high_cutoff=None, order=1): 1159 """ 1160 Low-pass and high-pass filtering become trivial in the frequency domain. 1161 A portion of the pixels/voxels are simply masked or attenuated. 1162 This function applies a high pass Butterworth filter that attenuates the 1163 frequency domain image. 1164 1165 The gradual attenuation of the filter is important. 1166 A simple high-pass filter would simply mask a set of pixels in the frequency domain, 1167 but the abrupt transition would cause a ringing effect in the spatial domain. 1168 1169 Arguments: 1170 low_cutoff : (list) 1171 the cutoff frequencies for x, y and z 1172 high_cutoff : (list) 1173 the cutoff frequencies for x, y and z 1174 order : (int) 1175 order determines sharpness of the cutoff curve 1176 """ 1177 # https://lorensen.github.io/VTKExamples/site/Cxx/ImageProcessing/IdealHighPass 1178 fft = vtk.new("ImageFFT") 1179 fft.SetInputData(self.dataset) 1180 fft.Update() 1181 out = fft.GetOutput() 1182 1183 if high_cutoff: 1184 blp = vtk.new("ImageButterworthLowPass") 1185 blp.SetInputData(out) 1186 blp.SetCutOff(high_cutoff) 1187 blp.SetOrder(order) 1188 blp.Update() 1189 out = blp.GetOutput() 1190 1191 if low_cutoff: 1192 bhp = vtk.new("ImageButterworthHighPass") 1193 bhp.SetInputData(out) 1194 bhp.SetCutOff(low_cutoff) 1195 bhp.SetOrder(order) 1196 bhp.Update() 1197 out = bhp.GetOutput() 1198 1199 rfft = vtk.new("ImageRFFT") 1200 rfft.SetInputData(out) 1201 rfft.Update() 1202 1203 ecomp = vtk.new("ImageExtractComponents") 1204 ecomp.SetInputData(rfft.GetOutput()) 1205 ecomp.SetComponents(0) 1206 ecomp.Update() 1207 self._update(ecomp.GetOutput()) 1208 self.pipeline = utils.OperationNode("frequency_pass_filter", parents=[self], c="#4cc9f0") 1209 return self 1210 1211 def smooth_gaussian(self, sigma=(2, 2, 2), radius=None): 1212 """ 1213 Performs a convolution of the input Volume with a gaussian. 1214 1215 Arguments: 1216 sigma : (float, list) 1217 standard deviation(s) in voxel units. 1218 A list can be given to smooth in the three direction differently. 1219 radius : (float, list) 1220 radius factor(s) determine how far out the gaussian 1221 kernel will go before being clamped to zero. A list can be given too. 1222 """ 1223 gsf = vtk.new("ImageGaussianSmooth") 1224 gsf.SetDimensionality(3) 1225 gsf.SetInputData(self.dataset) 1226 if utils.is_sequence(sigma): 1227 gsf.SetStandardDeviations(sigma) 1228 else: 1229 gsf.SetStandardDeviation(sigma) 1230 if radius is not None: 1231 if utils.is_sequence(radius): 1232 gsf.SetRadiusFactors(radius) 1233 else: 1234 gsf.SetRadiusFactor(radius) 1235 gsf.Update() 1236 self._update(gsf.GetOutput()) 1237 self.pipeline = utils.OperationNode("smooth_gaussian", parents=[self], c="#4cc9f0") 1238 return self 1239 1240 def smooth_median(self, neighbours=(2, 2, 2)): 1241 """ 1242 Median filter that replaces each pixel with the median value 1243 from a rectangular neighborhood around that pixel. 1244 """ 1245 imgm = vtk.new("ImageMedian3D") 1246 imgm.SetInputData(self.dataset) 1247 if utils.is_sequence(neighbours): 1248 imgm.SetKernelSize(neighbours[0], neighbours[1], neighbours[2]) 1249 else: 1250 imgm.SetKernelSize(neighbours, neighbours, neighbours) 1251 imgm.Update() 1252 self._update(imgm.GetOutput()) 1253 self.pipeline = utils.OperationNode("smooth_median", parents=[self], c="#4cc9f0") 1254 return self 1255 1256 def erode(self, neighbours=(2, 2, 2)): 1257 """ 1258 Replace a voxel with the minimum over an ellipsoidal neighborhood of voxels. 1259 If `neighbours` of an axis is 1, no processing is done on that axis. 1260 1261 Examples: 1262 - [erode_dilate.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/erode_dilate.py) 1263 1264  1265 """ 1266 ver = vtk.new("ImageContinuousErode3D") 1267 ver.SetInputData(self.dataset) 1268 ver.SetKernelSize(neighbours[0], neighbours[1], neighbours[2]) 1269 ver.Update() 1270 self._update(ver.GetOutput()) 1271 self.pipeline = utils.OperationNode("erode", parents=[self], c="#4cc9f0") 1272 return self 1273 1274 def dilate(self, neighbours=(2, 2, 2)): 1275 """ 1276 Replace a voxel with the maximum over an ellipsoidal neighborhood of voxels. 1277 If `neighbours` of an axis is 1, no processing is done on that axis. 1278 1279 Check also `erode()` and `pad()`. 1280 1281 Examples: 1282 - [erode_dilate.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/erode_dilate.py) 1283 """ 1284 ver = vtk.new("ImageContinuousDilate3D") 1285 ver.SetInputData(self.dataset) 1286 ver.SetKernelSize(neighbours[0], neighbours[1], neighbours[2]) 1287 ver.Update() 1288 self._update(ver.GetOutput()) 1289 self.pipeline = utils.OperationNode("dilate", parents=[self], c="#4cc9f0") 1290 return self 1291 1292 def magnitude(self): 1293 """Colapses components with magnitude function.""" 1294 imgm = vtk.new("ImageMagnitude") 1295 imgm.SetInputData(self.dataset) 1296 imgm.Update() 1297 self._update(imgm.GetOutput()) 1298 self.pipeline = utils.OperationNode("magnitude", parents=[self], c="#4cc9f0") 1299 return self 1300 1301 def topoints(self): 1302 """ 1303 Extract all image voxels as points. 1304 This function takes an input `Volume` and creates an `Mesh` 1305 that contains the points and the point attributes. 1306 1307 Examples: 1308 - [vol2points.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/vol2points.py) 1309 """ 1310 v2p = vtk.new("ImageToPoints") 1311 v2p.SetInputData(self.dataset) 1312 v2p.Update() 1313 mpts = vedo.Points(v2p.GetOutput()) 1314 mpts.pipeline = utils.OperationNode("topoints", parents=[self], c="#4cc9f0:#e9c46a") 1315 return mpts 1316 1317 def euclidean_distance(self, anisotropy=False, max_distance=None): 1318 """ 1319 Implementation of the Euclidean DT (Distance Transform) using Saito's algorithm. 1320 The distance map produced contains the square of the Euclidean distance values. 1321 The algorithm has a O(n^(D+1)) complexity over n x n x...x n images in D dimensions. 1322 1323 Check out also: https://en.wikipedia.org/wiki/Distance_transform 1324 1325 Arguments: 1326 anisotropy : bool 1327 used to define whether Spacing should be used in the 1328 computation of the distances. 1329 max_distance : bool 1330 any distance bigger than max_distance will not be 1331 computed but set to this specified value instead. 1332 1333 Examples: 1334 - [euclidian_dist.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/euclidian_dist.py) 1335 """ 1336 euv = vtk.new("ImageEuclideanDistance") 1337 euv.SetInputData(self.dataset) 1338 euv.SetConsiderAnisotropy(anisotropy) 1339 if max_distance is not None: 1340 euv.InitializeOn() 1341 euv.SetMaximumDistance(max_distance) 1342 euv.SetAlgorithmToSaito() 1343 euv.Update() 1344 vol = Volume(euv.GetOutput()) 1345 vol.pipeline = utils.OperationNode("euclidean_distance", parents=[self], c="#4cc9f0") 1346 return vol 1347 1348 def correlation_with(self, vol2, dim=2): 1349 """ 1350 Find the correlation between two volumetric data sets. 1351 Keyword `dim` determines whether the correlation will be 3D, 2D or 1D. 1352 The default is a 2D Correlation. 1353 1354 The output size will match the size of the first input. 1355 The second input is considered the correlation kernel. 1356 """ 1357 imc = vtk.new("ImageCorrelation") 1358 imc.SetInput1Data(self.dataset) 1359 imc.SetInput2Data(vol2.dataset) 1360 imc.SetDimensionality(dim) 1361 imc.Update() 1362 vol = Volume(imc.GetOutput()) 1363 1364 vol.pipeline = utils.OperationNode("correlation_with", parents=[self, vol2], c="#4cc9f0") 1365 return vol 1366 1367 def scale_voxels(self, scale=1): 1368 """Scale the voxel content by factor `scale`.""" 1369 rsl = vtk.new("ImageReslice") 1370 rsl.SetInputData(self.dataset) 1371 rsl.SetScalarScale(scale) 1372 rsl.Update() 1373 self._update(rsl.GetOutput()) 1374 self.pipeline = utils.OperationNode( 1375 "scale_voxels", comment=f"scale={scale}", parents=[self], c="#4cc9f0" 1376 ) 1377 return self
Class to describe dataset that are defined on "voxels", the 3D equivalent of 2D pixels.
35 def __init__( 36 self, 37 inputobj=None, 38 dims=None, 39 origin=None, 40 spacing=None, 41 ): 42 """ 43 This class can be initialized with a numpy object, 44 a `vtkImageData` or a list of 2D bmp files. 45 46 Arguments: 47 origin : (list) 48 set volume origin coordinates 49 spacing : (list) 50 voxel dimensions in x, y and z. 51 dims : (list) 52 specify the dimensions of the volume. 53 54 Example: 55 ```python 56 from vedo import Volume 57 vol = Volume("path/to/mydata/rec*.bmp") 58 vol.show() 59 ``` 60 61 Examples: 62 - [numpy2volume1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/numpy2volume1.py) 63 64  65 66 - [read_volume2.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/read_volume2.py) 67 68  69 70 .. note:: 71 if a `list` of values is used for `alphas` this is interpreted 72 as a transfer function along the range of the scalar. 73 """ 74 self.name = "Volume" 75 self.filename = "" 76 self.info = {} 77 self.time = 0 78 self.rendered_at = set() 79 80 self.actor = vtk.vtkVolume() 81 self.actor.retrieve_object = weak_ref_to(self) 82 self.properties = self.actor.GetProperty() 83 84 ################### 85 if isinstance(inputobj, str): 86 if "https://" in inputobj: 87 inputobj = vedo.file_io.download(inputobj, verbose=False) # fpath 88 elif os.path.isfile(inputobj): 89 self.filename = inputobj 90 else: 91 inputobj = sorted(glob.glob(inputobj)) 92 93 ################### 94 inputtype = str(type(inputobj)) 95 96 # print('Volume inputtype', inputtype, c='b') 97 98 if inputobj is None: 99 img = vtk.vtkImageData() 100 101 elif utils.is_sequence(inputobj): 102 103 if isinstance(inputobj[0], str) and ".bmp" in inputobj[0].lower(): 104 # scan sequence of BMP files 105 ima = vtk.new("ImageAppend") 106 ima.SetAppendAxis(2) 107 pb = utils.ProgressBar(0, len(inputobj)) 108 for i in pb.range(): 109 f = inputobj[i] 110 if "_rec_spr" in f: # OPT specific 111 continue 112 picr = vtk.new("BMPReader") 113 picr.SetFileName(f) 114 picr.Update() 115 mgf = vtk.new("ImageMagnitude") 116 mgf.SetInputData(picr.GetOutput()) 117 mgf.Update() 118 ima.AddInputData(mgf.GetOutput()) 119 pb.print("loading...") 120 ima.Update() 121 img = ima.GetOutput() 122 123 else: 124 125 if len(inputobj.shape) == 1: 126 varr = utils.numpy2vtk(inputobj) 127 else: 128 varr = utils.numpy2vtk(inputobj.ravel(order="F")) 129 varr.SetName("input_scalars") 130 131 img = vtk.vtkImageData() 132 if dims is not None: 133 img.SetDimensions(dims[2], dims[1], dims[0]) 134 else: 135 if len(inputobj.shape) == 1: 136 vedo.logger.error("must set dimensions (dims keyword) in Volume") 137 raise RuntimeError() 138 img.SetDimensions(inputobj.shape) 139 img.GetPointData().AddArray(varr) 140 img.GetPointData().SetActiveScalars(varr.GetName()) 141 142 elif isinstance(inputobj, vtk.vtkImageData): 143 img = inputobj 144 145 elif isinstance(inputobj, str): 146 if "https://" in inputobj: 147 inputobj = vedo.file_io.download(inputobj, verbose=False) 148 img = vedo.file_io.loadImageData(inputobj) 149 self.filename = inputobj 150 151 else: 152 vedo.logger.error(f"cannot understand input type {inputtype}") 153 return 154 155 if dims is not None: 156 img.SetDimensions(dims) 157 158 if origin is not None: 159 img.SetOrigin(origin) 160 161 if spacing is not None: 162 img.SetSpacing(spacing) 163 164 self.dataset = img 165 self.transform = None 166 167 ##################################### 168 mapper = vtk.new("SmartVolumeMapper") 169 mapper.SetInputData(img) 170 self.actor.SetMapper(mapper) 171 172 if img.GetPointData().GetScalars(): 173 if img.GetPointData().GetScalars().GetNumberOfComponents() == 1: 174 self.properties.SetShade(True) 175 self.properties.SetInterpolationType(1) 176 self.cmap("RdBu_r") 177 self.alpha([0.0, 0.0, 0.2, 0.4, 0.8, 1.0]) 178 self.alpha_gradient(None) 179 self.properties.SetScalarOpacityUnitDistance(1.0) 180 181 self.pipeline = utils.OperationNode( 182 "Volume", comment=f"dims={tuple(self.dimensions())}", c="#4cc9f0" 183 ) 184 #######################################################################
This class can be initialized with a numpy object,
a vtkImageData
or a list of 2D bmp files.
Arguments:
- origin : (list) set volume origin coordinates
- spacing : (list) voxel dimensions in x, y and z.
- dims : (list) specify the dimensions of the volume.
Example:
from vedo import Volume vol = Volume("path/to/mydata/rec*.bmp") vol.show()
Examples:
if a list
of values is used for alphas
this is interpreted
as a transfer function along the range of the scalar.
218 def c(self, *args, **kwargs): 219 """Deprecated. Use `Volume.cmap()` instead.""" 220 vedo.logger.warning("Volume.c() is deprecated, use Volume.cmap() instead") 221 return self.cmap(*args, **kwargs)
Deprecated. Use Volume.cmap()
instead.
359 def copy(self, deep=True): 360 """Return a copy of the Volume. Alias of `clone()`.""" 361 return self.clone(deep=deep)
Return a copy of the Volume. Alias of clone()
.
363 def clone(self, deep=True): 364 """Return a clone copy of the Volume. Alias of `copy()`.""" 365 if deep: 366 newimg = vtk.vtkImageData() 367 newimg.CopyStructure(self.dataset) 368 newimg.CopyAttributes(self.dataset) 369 newvol = Volume(newimg) 370 else: 371 newvol = Volume(self.dataset) 372 373 prop = vtk.vtkVolumeProperty() 374 prop.DeepCopy(self.properties) 375 newvol.actor.SetProperty(prop) 376 newvol.properties = prop 377 378 newvol.pipeline = utils.OperationNode("clone", parents=[self], c="#bbd0ff", shape="diamond") 379 return newvol
Return a clone copy of the Volume. Alias of copy()
.
381 def component_weight(self, i, weight): 382 """Set the scalar component weight in range [0,1].""" 383 self.properties.SetComponentWeight(i, weight) 384 return self
Set the scalar component weight in range [0,1].
386 def xslice(self, i): 387 """Extract the slice at index `i` of volume along x-axis.""" 388 vslice = vtk.new("ImageDataGeometryFilter") 389 vslice.SetInputData(self.dataset) 390 nx, ny, nz = self.dataset.GetDimensions() 391 if i > nx - 1: 392 i = nx - 1 393 vslice.SetExtent(i, i, 0, ny, 0, nz) 394 vslice.Update() 395 m = Mesh(vslice.GetOutput()) 396 m.pipeline = utils.OperationNode(f"xslice {i}", parents=[self], c="#4cc9f0:#e9c46a") 397 return m
Extract the slice at index i
of volume along x-axis.
399 def yslice(self, j): 400 """Extract the slice at index `j` of volume along y-axis.""" 401 vslice = vtk.new("ImageDataGeometryFilter") 402 vslice.SetInputData(self.dataset) 403 nx, ny, nz = self.dataset.GetDimensions() 404 if j > ny - 1: 405 j = ny - 1 406 vslice.SetExtent(0, nx, j, j, 0, nz) 407 vslice.Update() 408 m = Mesh(vslice.GetOutput()) 409 m.pipeline = utils.OperationNode(f"yslice {j}", parents=[self], c="#4cc9f0:#e9c46a") 410 return m
Extract the slice at index j
of volume along y-axis.
412 def zslice(self, k): 413 """Extract the slice at index `i` of volume along z-axis.""" 414 vslice = vtk.new("ImageDataGeometryFilter") 415 vslice.SetInputData(self.dataset) 416 nx, ny, nz = self.dataset.GetDimensions() 417 if k > nz - 1: 418 k = nz - 1 419 vslice.SetExtent(0, nx, 0, ny, k, k) 420 vslice.Update() 421 m = Mesh(vslice.GetOutput()) 422 m.pipeline = utils.OperationNode(f"zslice {k}", parents=[self], c="#4cc9f0:#e9c46a") 423 return m
Extract the slice at index i
of volume along z-axis.
425 def slice_plane(self, origin=(0, 0, 0), normal=(1, 1, 1), autocrop=False): 426 """ 427 Extract the slice along a given plane position and normal. 428 429 Example: 430 - [slice_plane1.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/slice_plane1.py) 431 432  433 """ 434 reslice = vtk.new("ImageReslice") 435 reslice.SetInputData(self.dataset) 436 reslice.SetOutputDimensionality(2) 437 newaxis = utils.versor(normal) 438 pos = np.array(origin) 439 initaxis = (0, 0, 1) 440 crossvec = np.cross(initaxis, newaxis) 441 angle = np.arccos(np.dot(initaxis, newaxis)) 442 T = vtk.vtkTransform() 443 T.PostMultiply() 444 T.RotateWXYZ(np.rad2deg(angle), crossvec) 445 T.Translate(pos) 446 M = T.GetMatrix() 447 reslice.SetResliceAxes(M) 448 reslice.SetInterpolationModeToLinear() 449 reslice.SetAutoCropOutput(not autocrop) 450 reslice.Update() 451 vslice = vtk.new("ImageDataGeometryFilter") 452 vslice.SetInputData(reslice.GetOutput()) 453 vslice.Update() 454 msh = Mesh(vslice.GetOutput()) 455 msh.apply_transform(T) 456 msh.pipeline = utils.OperationNode( 457 "slice_plane", parents=[self], c="#4cc9f0:#e9c46a") 458 return msh
460 def slab(self, slice_range=(), axis='z', operation="mean"): 461 """ 462 Extract a slab from a `Volume` by combining 463 all of the slices of an image to create a single slice. 464 465 Returns a `Mesh` containing metadata which 466 can be accessed with e.g. `mesh.metadata["slab_range"]`. 467 468 Metadata: 469 slab_range : (list) 470 contains the range of slices extracted 471 slab_axis : (str) 472 contains the axis along which the slab was extracted 473 slab_operation : (str) 474 contains the operation performed on the slab 475 slab_bounding_box : (list) 476 contains the bounding box of the slab 477 478 Arguments: 479 slice_range : (list) 480 range of slices to extract 481 axis : (str) 482 axis along which to extract the slab 483 operation : (str) 484 operation to perform on the slab, 485 allowed values are: "sum", "min", "max", "mean". 486 """ 487 if len(slice_range) != 2: 488 vedo.logger.error("in slab(): slice_range is empty or invalid") 489 raise ValueError() 490 491 islab = vtk.new("ImageSlab") 492 islab.SetInputData(self.dataset) 493 494 if operation in ["+", "add", "sum"]: 495 islab.SetOperationToSum() 496 elif "min" in operation: 497 islab.SetOperationToMin() 498 elif "max" in operation: 499 islab.SetOperationToMax() 500 elif "mean" in operation: 501 islab.SetOperationToMean() 502 else: 503 vedo.logger.error(f"in slab(): unknown operation {operation}") 504 raise ValueError() 505 506 dims = self.dimensions() 507 if axis == 'x': 508 islab.SetOrientationToX() 509 if slice_range[0] > dims[0]-1: 510 slice_range[0] = dims[0]-1 511 if slice_range[1] > dims[0]-1: 512 slice_range[1] = dims[0]-1 513 elif axis == 'y': 514 islab.SetOrientationToY() 515 if slice_range[0] > dims[1]-1: 516 slice_range[0] = dims[1]-1 517 if slice_range[1] > dims[1]-1: 518 slice_range[1] = dims[1]-1 519 elif axis == 'z': 520 islab.SetOrientationToZ() 521 if slice_range[0] > dims[2]-1: 522 slice_range[0] = dims[2]-1 523 if slice_range[1] > dims[2]-1: 524 slice_range[1] = dims[2]-1 525 else: 526 vedo.logger.error(f"Error in slab(): unknown axis {axis}") 527 raise RuntimeError() 528 529 islab.SetSliceRange(slice_range) 530 islab.Update() 531 532 msh = Mesh(islab.GetOutput()).lighting('off') 533 msh.mapper.SetLookupTable(utils.ctf2lut(self, msh)) 534 msh.mapper.SetScalarRange(self.scalar_range()) 535 536 msh.metadata["slab_range"] = slice_range 537 msh.metadata["slab_axis"] = axis 538 msh.metadata["slab_operation"] = operation 539 540 # compute bounds of slab 541 origin = self.origin() 542 spacing = self.spacing() 543 if axis == 'x': 544 msh.metadata["slab_bounding_box"] = [ 545 origin[0] + slice_range[0]*spacing[0], 546 origin[0] + slice_range[1]*spacing[0], 547 origin[1], 548 origin[1] + dims[1]*spacing[1], 549 origin[2], 550 origin[2] + dims[2]*spacing[2], 551 ] 552 elif axis == 'y': 553 msh.metadata["slab_bounding_box"] = [ 554 origin[0], 555 origin[0] + dims[0]*spacing[0], 556 origin[1] + slice_range[0]*spacing[1], 557 origin[1] + slice_range[1]*spacing[1], 558 origin[2], 559 origin[2] + dims[2]*spacing[2], 560 ] 561 elif axis == 'z': 562 msh.metadata["slab_bounding_box"] = [ 563 origin[0], 564 origin[0] + dims[0]*spacing[0], 565 origin[1], 566 origin[1] + dims[1]*spacing[1], 567 origin[2] + slice_range[0]*spacing[2], 568 origin[2] + slice_range[1]*spacing[2], 569 ] 570 571 msh.pipeline = utils.OperationNode( 572 f"slab{slice_range}", 573 comment=f"axis={axis}, operation={operation}", 574 parents=[self], 575 c="#4cc9f0:#e9c46a", 576 ) 577 msh.name = "SlabMesh" 578 return msh
Extract a slab from a Volume
by combining
all of the slices of an image to create a single slice.
Returns a Mesh
containing metadata which
can be accessed with e.g. mesh.metadata["slab_range"]
.
Metadata:
slab_range : (list) contains the range of slices extracted slab_axis : (str) contains the axis along which the slab was extracted slab_operation : (str) contains the operation performed on the slab slab_bounding_box : (list) contains the bounding box of the slab
Arguments:
- slice_range : (list) range of slices to extract
- axis : (str) axis along which to extract the slab
- operation : (str) operation to perform on the slab, allowed values are: "sum", "min", "max", "mean".
581 def warp(self, source, target, sigma=1, mode="3d", fit=True): 582 """ 583 Warp volume scalars within a Volume by specifying 584 source and target sets of points. 585 586 Arguments: 587 source : (Points, list) 588 the list of source points 589 target : (Points, list) 590 the list of target points 591 fit : (bool) 592 fit/adapt the old bounding box to the warped geometry 593 """ 594 if isinstance(source, vedo.Points): 595 source = source.vertices 596 if isinstance(target, vedo.Points): 597 target = target.vertices 598 599 NLT = transformations.NonLinearTransform() 600 NLT.source_points = source 601 NLT.target_points = target 602 NLT.sigma = sigma 603 NLT.mode = mode 604 605 self.apply_transform(NLT, fit=fit) 606 self.pipeline = utils.OperationNode("warp", parents=[self], c="#4cc9f0") 607 return self
Warp volume scalars within a Volume by specifying source and target sets of points.
Arguments:
- source : (Points, list) the list of source points
- target : (Points, list) the list of target points
- fit : (bool) fit/adapt the old bounding box to the warped geometry
609 def apply_transform(self, T, fit=True, interpolation="cubic"): 610 """ 611 Apply a transform to the scalars in the volume. 612 613 Arguments: 614 T : (LinearTransform, NonLinearTransform) 615 The transformation to be applied 616 fit : (bool) 617 fit/adapt the old bounding box to the modified geometry 618 interpolation : (str) 619 one of the following: "linear", "nearest", "cubic" 620 """ 621 TI = T.compute_inverse() 622 reslice = vtk.new("ImageReslice") 623 reslice.SetInputData(self.dataset) 624 reslice.SetResliceTransform(TI.T) 625 reslice.SetOutputDimensionality(3) 626 if "lin" in interpolation.lower(): 627 reslice.SetInterpolationModeToLinear() 628 elif "near" in interpolation.lower(): 629 reslice.SetInterpolationModeToNearestNeighbor() 630 elif "cubic" in interpolation.lower(): 631 reslice.SetInterpolationModeToCubic() 632 else: 633 vedo.logger.error( 634 f"in apply_transform: unknown interpolation mode {interpolation}") 635 raise ValueError() 636 reslice.SetAutoCropOutput(fit) 637 reslice.Update() 638 self._update(reslice.GetOutput()) 639 self.transform = T 640 self.pipeline = utils.OperationNode( 641 "apply_transform", parents=[self], c="#4cc9f0") 642 return self
Apply a transform to the scalars in the volume.
Arguments:
- T : (LinearTransform, NonLinearTransform) The transformation to be applied
- fit : (bool) fit/adapt the old bounding box to the modified geometry
- interpolation : (str) one of the following: "linear", "nearest", "cubic"
644 def imagedata(self): 645 """ 646 DEPRECATED: 647 Use `Volume.dataset` instead. 648 649 Return the underlying `vtkImagaData` object. 650 """ 651 print("Volume.imagedata() is deprecated, use Volume.dataset instead") 652 return self.dataset
DEPRECATED:
Use Volume.dataset
instead.
Return the underlying vtkImagaData
object.
665 def tonumpy(self): 666 """ 667 Get read-write access to voxels of a Volume object as a numpy array. 668 669 When you set values in the output image, you don't want numpy to reallocate the array 670 but instead set values in the existing array, so use the [:] operator. 671 672 Example: 673 `arr[:] = arr*2 + 15` 674 675 If the array is modified add a call to: 676 `volume.modified()` 677 when all your modifications are completed. 678 """ 679 narray_shape = tuple(reversed(self.dataset.GetDimensions())) 680 681 scals = self.dataset.GetPointData().GetScalars() 682 comps = scals.GetNumberOfComponents() 683 if comps == 1: 684 narray = utils.vtk2numpy(scals).reshape(narray_shape) 685 narray = np.transpose(narray, axes=[2, 1, 0]) 686 else: 687 narray = utils.vtk2numpy(scals).reshape(*narray_shape, comps) 688 narray = np.transpose(narray, axes=[2, 1, 0, 3]) 689 690 # narray = utils.vtk2numpy(self.dataset.GetPointData().GetScalars()).reshape(narray_shape) 691 # narray = np.transpose(narray, axes=[2, 1, 0]) 692 693 return narray
Get read-write access to voxels of a Volume object as a numpy array.
When you set values in the output image, you don't want numpy to reallocate the array but instead set values in the existing array, so use the [:] operator.
Example:
arr[:] = arr*2 + 15
If the array is modified add a call to:
volume.modified()
when all your modifications are completed.
700 def dimensions(self): 701 """Return the nr. of voxels in the 3 dimensions.""" 702 return np.array(self.dataset.GetDimensions())
Return the nr. of voxels in the 3 dimensions.
704 def scalar_range(self): 705 """Return the range of the scalar values.""" 706 return np.array(self.dataset.GetScalarRange())
Return the range of the scalar values.
708 def spacing(self, s=None): 709 """Set/get the voxels size in the 3 dimensions.""" 710 if s is not None: 711 self.dataset.SetSpacing(s) 712 return self 713 return np.array(self.dataset.GetSpacing())
Set/get the voxels size in the 3 dimensions.
715 def origin(self, s=None): 716 """ 717 Set/get the origin of the volumetric dataset. 718 719 The origin is the position in world coordinates of the point index (0,0,0). 720 This point does not have to be part of the dataset, in other words, 721 the dataset extent does not have to start at (0,0,0) and the origin 722 can be outside of the dataset bounding box. 723 The origin plus spacing determine the position in space of the points. 724 """ 725 if s is not None: 726 self.dataset.SetOrigin(s) 727 return self 728 return np.array(self.dataset.GetOrigin())
Set/get the origin of the volumetric dataset.
The origin is the position in world coordinates of the point index (0,0,0). This point does not have to be part of the dataset, in other words, the dataset extent does not have to start at (0,0,0) and the origin can be outside of the dataset bounding box. The origin plus spacing determine the position in space of the points.
730 def center(self): 731 """Get the center of the volumetric dataset.""" 732 # note that this does not have the set method like origin and spacing 733 return np.array(self.dataset.GetCenter())
Get the center of the volumetric dataset.
735 def get_cell_from_ijk(self, ijk): 736 """ 737 Get the voxel id number at the given ijk coordinates. 738 739 Arguments: 740 ijk : (list) 741 the ijk coordinates of the voxel 742 """ 743 return self.ComputeCellId(ijk)
Get the voxel id number at the given ijk coordinates.
Arguments:
- ijk : (list) the ijk coordinates of the voxel
745 def get_point_from_ijk(self, ijk): 746 """ 747 Get the point id number at the given ijk coordinates. 748 749 Arguments: 750 ijk : (list) 751 the ijk coordinates of the voxel 752 """ 753 return self.ComputePointId(ijk)
Get the point id number at the given ijk coordinates.
Arguments:
- ijk : (list) the ijk coordinates of the voxel
755 def permute_axes(self, x, y, z): 756 """ 757 Reorder the axes of the Volume by specifying 758 the input axes which are supposed to become the new X, Y, and Z. 759 """ 760 imp = vtk.new("ImagePermute") 761 imp.SetFilteredAxes(x, y, z) 762 imp.SetInputData(self.dataset) 763 imp.Update() 764 self._update(imp.GetOutput()) 765 self.pipeline = utils.OperationNode( 766 f"permute_axes({(x,y,z)})", parents=[self], c="#4cc9f0" 767 ) 768 return self
Reorder the axes of the Volume by specifying the input axes which are supposed to become the new X, Y, and Z.
770 def resample(self, new_spacing, interpolation=1): 771 """ 772 Resamples a `Volume` to be larger or smaller. 773 774 This method modifies the spacing of the input. 775 Linear interpolation is used to resample the data. 776 777 Arguments: 778 new_spacing : (list) 779 a list of 3 new spacings for the 3 axes 780 interpolation : (int) 781 0=nearest_neighbor, 1=linear, 2=cubic 782 """ 783 rsp = vtk.new("ImageResample") 784 oldsp = self.spacing() 785 for i in range(3): 786 if oldsp[i] != new_spacing[i]: 787 rsp.SetAxisOutputSpacing(i, new_spacing[i]) 788 rsp.InterpolateOn() 789 rsp.SetInterpolationMode(interpolation) 790 rsp.OptimizationOn() 791 rsp.Update() 792 self._update(rsp.GetOutput()) 793 self.pipeline = utils.OperationNode( 794 "resample", comment=f"spacing: {tuple(new_spacing)}", parents=[self], c="#4cc9f0" 795 ) 796 return self
Resamples a Volume
to be larger or smaller.
This method modifies the spacing of the input. Linear interpolation is used to resample the data.
Arguments:
- new_spacing : (list) a list of 3 new spacings for the 3 axes
- interpolation : (int) 0=nearest_neighbor, 1=linear, 2=cubic
799 def threshold(self, above=None, below=None, replace=None, replace_value=None): 800 """ 801 Binary or continuous volume thresholding. 802 Find the voxels that contain a value above/below the input values 803 and replace them with a new value (default is 0). 804 """ 805 th = vtk.new("ImageThreshold") 806 th.SetInputData(self.dataset) 807 808 # sanity checks 809 if above is not None and below is not None: 810 if above == below: 811 return self 812 if above > below: 813 vedo.logger.warning("in volume.threshold(), above > below, skip.") 814 return self 815 816 ## cases 817 if below is not None and above is not None: 818 th.ThresholdBetween(above, below) 819 820 elif above is not None: 821 th.ThresholdByUpper(above) 822 823 elif below is not None: 824 th.ThresholdByLower(below) 825 826 ## 827 if replace is not None: 828 th.SetReplaceIn(True) 829 th.SetInValue(replace) 830 else: 831 th.SetReplaceIn(False) 832 833 if replace_value is not None: 834 th.SetReplaceOut(True) 835 th.SetOutValue(replace_value) 836 else: 837 th.SetReplaceOut(False) 838 839 th.Update() 840 self._update(th.GetOutput()) 841 self.pipeline = utils.OperationNode("threshold", parents=[self], c="#4cc9f0") 842 return self
Binary or continuous volume thresholding. Find the voxels that contain a value above/below the input values and replace them with a new value (default is 0).
844 def crop(self, left=None, right=None, back=None, front=None, bottom=None, top=None, VOI=()): 845 """ 846 Crop a `Volume` object. 847 848 Arguments: 849 left : (float) 850 fraction to crop from the left plane (negative x) 851 right : (float) 852 fraction to crop from the right plane (positive x) 853 back : (float) 854 fraction to crop from the back plane (negative y) 855 front : (float) 856 fraction to crop from the front plane (positive y) 857 bottom : (float) 858 fraction to crop from the bottom plane (negative z) 859 top : (float) 860 fraction to crop from the top plane (positive z) 861 VOI : (list) 862 extract Volume Of Interest expressed in voxel numbers 863 864 Example: 865 `vol.crop(VOI=(xmin, xmax, ymin, ymax, zmin, zmax)) # all integers nrs` 866 """ 867 extractVOI = vtk.new("ExtractVOI") 868 extractVOI.SetInputData(self.dataset) 869 870 if VOI: 871 extractVOI.SetVOI(VOI) 872 else: 873 d = self.dataset.GetDimensions() 874 bx0, bx1, by0, by1, bz0, bz1 = 0, d[0]-1, 0, d[1]-1, 0, d[2]-1 875 if left is not None: bx0 = int((d[0]-1)*left) 876 if right is not None: bx1 = int((d[0]-1)*(1-right)) 877 if back is not None: by0 = int((d[1]-1)*back) 878 if front is not None: by1 = int((d[1]-1)*(1-front)) 879 if bottom is not None: bz0 = int((d[2]-1)*bottom) 880 if top is not None: bz1 = int((d[2]-1)*(1-top)) 881 extractVOI.SetVOI(bx0, bx1, by0, by1, bz0, bz1) 882 extractVOI.Update() 883 self._update(extractVOI.GetOutput()) 884 885 self.pipeline = utils.OperationNode( 886 "crop", parents=[self], c="#4cc9f0", comment=f"dims={tuple(self.dimensions())}" 887 ) 888 return self
Crop a Volume
object.
Arguments:
- left : (float) fraction to crop from the left plane (negative x)
- right : (float) fraction to crop from the right plane (positive x)
- back : (float) fraction to crop from the back plane (negative y)
- front : (float) fraction to crop from the front plane (positive y)
- bottom : (float) fraction to crop from the bottom plane (negative z)
- top : (float) fraction to crop from the top plane (positive z)
- VOI : (list) extract Volume Of Interest expressed in voxel numbers
Example:
vol.crop(VOI=(xmin, xmax, ymin, ymax, zmin, zmax)) # all integers nrs
890 def append(self, volumes, axis="z", preserve_extents=False): 891 """ 892 Take the components from multiple inputs and merges them into one output. 893 Except for the append axis, all inputs must have the same extent. 894 All inputs must have the same number of scalar components. 895 The output has the same origin and spacing as the first input. 896 The origin and spacing of all other inputs are ignored. 897 All inputs must have the same scalar type. 898 899 Arguments: 900 axis : (int, str) 901 axis expanded to hold the multiple images 902 preserve_extents : (bool) 903 if True, the extent of the inputs is used to place 904 the image in the output. The whole extent of the output is the union of the input 905 whole extents. Any portion of the output not covered by the inputs is set to zero. 906 The origin and spacing is taken from the first input. 907 908 Example: 909 ```python 910 from vedo import Volume, dataurl 911 vol = Volume(dataurl+'embryo.tif') 912 vol.append(vol, axis='x').show().close() 913 ``` 914  915 """ 916 ima = vtk.new("ImageAppend") 917 ima.SetInputData(self.dataset) 918 if not utils.is_sequence(volumes): 919 volumes = [volumes] 920 for volume in volumes: 921 if isinstance(volume, vtk.vtkImageData): 922 ima.AddInputData(volume) 923 else: 924 ima.AddInputData(volume.dataset) 925 ima.SetPreserveExtents(preserve_extents) 926 if axis == "x": 927 axis = 0 928 elif axis == "y": 929 axis = 1 930 elif axis == "z": 931 axis = 2 932 ima.SetAppendAxis(axis) 933 ima.Update() 934 self._update(ima.GetOutput()) 935 936 self.pipeline = utils.OperationNode( 937 "append", 938 parents=[self, *volumes], 939 c="#4cc9f0", 940 comment=f"dims={tuple(self.dimensions())}", 941 ) 942 return self
Take the components from multiple inputs and merges them into one output. Except for the append axis, all inputs must have the same extent. All inputs must have the same number of scalar components. The output has the same origin and spacing as the first input. The origin and spacing of all other inputs are ignored. All inputs must have the same scalar type.
Arguments:
- axis : (int, str) axis expanded to hold the multiple images
- preserve_extents : (bool) if True, the extent of the inputs is used to place the image in the output. The whole extent of the output is the union of the input whole extents. Any portion of the output not covered by the inputs is set to zero. The origin and spacing is taken from the first input.
Example:
from vedo import Volume, dataurl vol = Volume(dataurl+'embryo.tif') vol.append(vol, axis='x').show().close()
944 def pad(self, voxels=10, value=0): 945 """ 946 Add the specified number of voxels at the `Volume` borders. 947 Voxels can be a list formatted as `[nx0, nx1, ny0, ny1, nz0, nz1]`. 948 949 Arguments: 950 voxels : (int, list) 951 number of voxels to be added (or a list of length 4) 952 value : (int) 953 intensity value (gray-scale color) of the padding 954 955 Example: 956 ```python 957 from vedo import Volume, dataurl, show 958 iso = Volume(dataurl+'embryo.tif').isosurface() 959 vol = iso.binarize(spacing=(100, 100, 100)).pad(10) 960 vol.dilate([15,15,15]) 961 show(iso, vol.isosurface(), N=2, axes=1) 962 ``` 963  964 """ 965 x0, x1, y0, y1, z0, z1 = self.dataset.GetExtent() 966 pf = vtk.new("ImageConstantPad") 967 pf.SetInputData(self.dataset) 968 pf.SetConstant(value) 969 if utils.is_sequence(voxels): 970 pf.SetOutputWholeExtent( 971 x0 - voxels[0], x1 + voxels[1], 972 y0 - voxels[2], y1 + voxels[3], 973 z0 - voxels[4], z1 + voxels[5], 974 ) 975 else: 976 pf.SetOutputWholeExtent( 977 x0 - voxels, x1 + voxels, 978 y0 - voxels, y1 + voxels, 979 z0 - voxels, z1 + voxels, 980 ) 981 pf.Update() 982 self._update(pf.GetOutput()) 983 self.pipeline = utils.OperationNode( 984 "pad", comment=f"{voxels} voxels", parents=[self], c="#f28482" 985 ) 986 return self
Add the specified number of voxels at the Volume
borders.
Voxels can be a list formatted as [nx0, nx1, ny0, ny1, nz0, nz1]
.
Arguments:
- voxels : (int, list) number of voxels to be added (or a list of length 4)
- value : (int) intensity value (gray-scale color) of the padding
Example:
from vedo import Volume, dataurl, show iso = Volume(dataurl+'embryo.tif').isosurface() vol = iso.binarize(spacing=(100, 100, 100)).pad(10) vol.dilate([15,15,15]) show(iso, vol.isosurface(), N=2, axes=1)
988 def resize(self, newdims): 989 """Increase or reduce the number of voxels of a Volume with interpolation.""" 990 rsz = vtk.new("ImageResize") 991 rsz.SetResizeMethodToOutputDimensions() 992 rsz.SetInputData(self.dataset) 993 rsz.SetOutputDimensions(newdims) 994 rsz.Update() 995 self.dataset = rsz.GetOutput() 996 self._update(self.dataset) 997 self.pipeline = utils.OperationNode( 998 "resize", parents=[self], c="#4cc9f0", comment=f"dims={tuple(self.dimensions())}" 999 ) 1000 return self
Increase or reduce the number of voxels of a Volume with interpolation.
1002 def normalize(self): 1003 """Normalize that scalar components for each point.""" 1004 norm = vtk.new("ImageNormalize") 1005 norm.SetInputData(self.dataset) 1006 norm.Update() 1007 self._update(norm.GetOutput()) 1008 self.pipeline = utils.OperationNode("normalize", parents=[self], c="#4cc9f0") 1009 return self
Normalize that scalar components for each point.
1011 def mirror(self, axis="x"): 1012 """ 1013 Mirror flip along one of the cartesian axes. 1014 """ 1015 img = self.dataset 1016 1017 ff = vtk.new("ImageFlip") 1018 ff.SetInputData(img) 1019 if axis.lower() == "x": 1020 ff.SetFilteredAxis(0) 1021 elif axis.lower() == "y": 1022 ff.SetFilteredAxis(1) 1023 elif axis.lower() == "z": 1024 ff.SetFilteredAxis(2) 1025 else: 1026 vedo.logger.error("mirror must be set to either x, y, z or n") 1027 raise RuntimeError() 1028 ff.Update() 1029 self._update(ff.GetOutput()) 1030 self.pipeline = utils.OperationNode(f"mirror {axis}", parents=[self], c="#4cc9f0") 1031 return self
Mirror flip along one of the cartesian axes.
1033 def operation(self, operation, volume2=None): 1034 """ 1035 Perform operations with `Volume` objects. 1036 Keyword `volume2` can be a constant float. 1037 1038 Possible operations are: 1039 ``` 1040 +, -, /, 1/x, sin, cos, exp, log, 1041 abs, **2, sqrt, min, max, atan, atan2, median, 1042 mag, dot, gradient, divergence, laplacian. 1043 ``` 1044 1045 Examples: 1046 - [volume_operations.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/volume_operations.py) 1047 """ 1048 op = operation.lower() 1049 image1 = self.dataset 1050 1051 mf = None 1052 if op in ["median"]: 1053 mf = vtk.new("ImageMedian3D") 1054 mf.SetInputData(image1) 1055 elif op in ["mag"]: 1056 mf = vtk.new("ImageMagnitude") 1057 mf.SetInputData(image1) 1058 elif op in ["dot", "dotproduct"]: 1059 mf = vtk.new("ImageDotProduct") 1060 mf.SetInput1Data(image1) 1061 mf.SetInput2Data(volume2.dataset) 1062 elif op in ["grad", "gradient"]: 1063 mf = vtk.new("ImageGradient") 1064 mf.SetDimensionality(3) 1065 mf.SetInputData(image1) 1066 elif op in ["div", "divergence"]: 1067 mf = vtk.new("ImageDivergence") 1068 mf.SetInputData(image1) 1069 elif op in ["laplacian"]: 1070 mf = vtk.new("ImageLaplacian") 1071 mf.SetDimensionality(3) 1072 mf.SetInputData(image1) 1073 1074 if mf is not None: 1075 mf.Update() 1076 vol = Volume(mf.GetOutput()) 1077 vol.pipeline = utils.OperationNode( 1078 "operation", comment=f"{op}", parents=[self], c="#4cc9f0", shape="cylinder" 1079 ) 1080 return vol ########################### 1081 1082 mat = vtk.new("ImageMathematics") 1083 mat.SetInput1Data(image1) 1084 1085 K = None 1086 1087 if utils.is_number(volume2): 1088 K = volume2 1089 mat.SetConstantK(K) 1090 mat.SetConstantC(K) 1091 1092 elif volume2 is not None: # assume image2 is a constant value 1093 mat.SetInput2Data(volume2.dataset) 1094 1095 # ########################### 1096 if op in ["+", "add", "plus"]: 1097 if K: 1098 mat.SetOperationToAddConstant() 1099 else: 1100 mat.SetOperationToAdd() 1101 1102 elif op in ["-", "subtract", "minus"]: 1103 if K: 1104 mat.SetConstantC(-float(K)) 1105 mat.SetOperationToAddConstant() 1106 else: 1107 mat.SetOperationToSubtract() 1108 1109 elif op in ["*", "multiply", "times"]: 1110 if K: 1111 mat.SetOperationToMultiplyByK() 1112 else: 1113 mat.SetOperationToMultiply() 1114 1115 elif op in ["/", "divide"]: 1116 if K: 1117 mat.SetConstantK(1.0 / K) 1118 mat.SetOperationToMultiplyByK() 1119 else: 1120 mat.SetOperationToDivide() 1121 1122 elif op in ["1/x", "invert"]: 1123 mat.SetOperationToInvert() 1124 elif op in ["sin"]: 1125 mat.SetOperationToSin() 1126 elif op in ["cos"]: 1127 mat.SetOperationToCos() 1128 elif op in ["exp"]: 1129 mat.SetOperationToExp() 1130 elif op in ["log"]: 1131 mat.SetOperationToLog() 1132 elif op in ["abs"]: 1133 mat.SetOperationToAbsoluteValue() 1134 elif op in ["**2", "square"]: 1135 mat.SetOperationToSquare() 1136 elif op in ["sqrt", "sqr"]: 1137 mat.SetOperationToSquareRoot() 1138 elif op in ["min"]: 1139 mat.SetOperationToMin() 1140 elif op in ["max"]: 1141 mat.SetOperationToMax() 1142 elif op in ["atan"]: 1143 mat.SetOperationToATAN() 1144 elif op in ["atan2"]: 1145 mat.SetOperationToATAN2() 1146 else: 1147 vedo.logger.error(f"unknown operation {operation}") 1148 raise RuntimeError() 1149 mat.Update() 1150 1151 self._update(mat.GetOutput()) 1152 1153 self.pipeline = utils.OperationNode( 1154 "operation", comment=f"{op}", parents=[self, volume2], shape="cylinder", c="#4cc9f0" 1155 ) 1156 return self
Perform operations with Volume
objects.
Keyword volume2
can be a constant float.
Possible operations are:
+, -, /, 1/x, sin, cos, exp, log,
abs, **2, sqrt, min, max, atan, atan2, median,
mag, dot, gradient, divergence, laplacian.
Examples:
1158 def frequency_pass_filter(self, low_cutoff=None, high_cutoff=None, order=1): 1159 """ 1160 Low-pass and high-pass filtering become trivial in the frequency domain. 1161 A portion of the pixels/voxels are simply masked or attenuated. 1162 This function applies a high pass Butterworth filter that attenuates the 1163 frequency domain image. 1164 1165 The gradual attenuation of the filter is important. 1166 A simple high-pass filter would simply mask a set of pixels in the frequency domain, 1167 but the abrupt transition would cause a ringing effect in the spatial domain. 1168 1169 Arguments: 1170 low_cutoff : (list) 1171 the cutoff frequencies for x, y and z 1172 high_cutoff : (list) 1173 the cutoff frequencies for x, y and z 1174 order : (int) 1175 order determines sharpness of the cutoff curve 1176 """ 1177 # https://lorensen.github.io/VTKExamples/site/Cxx/ImageProcessing/IdealHighPass 1178 fft = vtk.new("ImageFFT") 1179 fft.SetInputData(self.dataset) 1180 fft.Update() 1181 out = fft.GetOutput() 1182 1183 if high_cutoff: 1184 blp = vtk.new("ImageButterworthLowPass") 1185 blp.SetInputData(out) 1186 blp.SetCutOff(high_cutoff) 1187 blp.SetOrder(order) 1188 blp.Update() 1189 out = blp.GetOutput() 1190 1191 if low_cutoff: 1192 bhp = vtk.new("ImageButterworthHighPass") 1193 bhp.SetInputData(out) 1194 bhp.SetCutOff(low_cutoff) 1195 bhp.SetOrder(order) 1196 bhp.Update() 1197 out = bhp.GetOutput() 1198 1199 rfft = vtk.new("ImageRFFT") 1200 rfft.SetInputData(out) 1201 rfft.Update() 1202 1203 ecomp = vtk.new("ImageExtractComponents") 1204 ecomp.SetInputData(rfft.GetOutput()) 1205 ecomp.SetComponents(0) 1206 ecomp.Update() 1207 self._update(ecomp.GetOutput()) 1208 self.pipeline = utils.OperationNode("frequency_pass_filter", parents=[self], c="#4cc9f0") 1209 return self
Low-pass and high-pass filtering become trivial in the frequency domain. A portion of the pixels/voxels are simply masked or attenuated. This function applies a high pass Butterworth filter that attenuates the frequency domain image.
The gradual attenuation of the filter is important. A simple high-pass filter would simply mask a set of pixels in the frequency domain, but the abrupt transition would cause a ringing effect in the spatial domain.
Arguments:
- low_cutoff : (list) the cutoff frequencies for x, y and z
- high_cutoff : (list) the cutoff frequencies for x, y and z
- order : (int) order determines sharpness of the cutoff curve
1211 def smooth_gaussian(self, sigma=(2, 2, 2), radius=None): 1212 """ 1213 Performs a convolution of the input Volume with a gaussian. 1214 1215 Arguments: 1216 sigma : (float, list) 1217 standard deviation(s) in voxel units. 1218 A list can be given to smooth in the three direction differently. 1219 radius : (float, list) 1220 radius factor(s) determine how far out the gaussian 1221 kernel will go before being clamped to zero. A list can be given too. 1222 """ 1223 gsf = vtk.new("ImageGaussianSmooth") 1224 gsf.SetDimensionality(3) 1225 gsf.SetInputData(self.dataset) 1226 if utils.is_sequence(sigma): 1227 gsf.SetStandardDeviations(sigma) 1228 else: 1229 gsf.SetStandardDeviation(sigma) 1230 if radius is not None: 1231 if utils.is_sequence(radius): 1232 gsf.SetRadiusFactors(radius) 1233 else: 1234 gsf.SetRadiusFactor(radius) 1235 gsf.Update() 1236 self._update(gsf.GetOutput()) 1237 self.pipeline = utils.OperationNode("smooth_gaussian", parents=[self], c="#4cc9f0") 1238 return self
Performs a convolution of the input Volume with a gaussian.
Arguments:
- sigma : (float, list) standard deviation(s) in voxel units. A list can be given to smooth in the three direction differently.
- radius : (float, list) radius factor(s) determine how far out the gaussian kernel will go before being clamped to zero. A list can be given too.
1240 def smooth_median(self, neighbours=(2, 2, 2)): 1241 """ 1242 Median filter that replaces each pixel with the median value 1243 from a rectangular neighborhood around that pixel. 1244 """ 1245 imgm = vtk.new("ImageMedian3D") 1246 imgm.SetInputData(self.dataset) 1247 if utils.is_sequence(neighbours): 1248 imgm.SetKernelSize(neighbours[0], neighbours[1], neighbours[2]) 1249 else: 1250 imgm.SetKernelSize(neighbours, neighbours, neighbours) 1251 imgm.Update() 1252 self._update(imgm.GetOutput()) 1253 self.pipeline = utils.OperationNode("smooth_median", parents=[self], c="#4cc9f0") 1254 return self
Median filter that replaces each pixel with the median value from a rectangular neighborhood around that pixel.
1256 def erode(self, neighbours=(2, 2, 2)): 1257 """ 1258 Replace a voxel with the minimum over an ellipsoidal neighborhood of voxels. 1259 If `neighbours` of an axis is 1, no processing is done on that axis. 1260 1261 Examples: 1262 - [erode_dilate.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/erode_dilate.py) 1263 1264  1265 """ 1266 ver = vtk.new("ImageContinuousErode3D") 1267 ver.SetInputData(self.dataset) 1268 ver.SetKernelSize(neighbours[0], neighbours[1], neighbours[2]) 1269 ver.Update() 1270 self._update(ver.GetOutput()) 1271 self.pipeline = utils.OperationNode("erode", parents=[self], c="#4cc9f0") 1272 return self
Replace a voxel with the minimum over an ellipsoidal neighborhood of voxels.
If neighbours
of an axis is 1, no processing is done on that axis.
Examples:
1274 def dilate(self, neighbours=(2, 2, 2)): 1275 """ 1276 Replace a voxel with the maximum over an ellipsoidal neighborhood of voxels. 1277 If `neighbours` of an axis is 1, no processing is done on that axis. 1278 1279 Check also `erode()` and `pad()`. 1280 1281 Examples: 1282 - [erode_dilate.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/erode_dilate.py) 1283 """ 1284 ver = vtk.new("ImageContinuousDilate3D") 1285 ver.SetInputData(self.dataset) 1286 ver.SetKernelSize(neighbours[0], neighbours[1], neighbours[2]) 1287 ver.Update() 1288 self._update(ver.GetOutput()) 1289 self.pipeline = utils.OperationNode("dilate", parents=[self], c="#4cc9f0") 1290 return self
Replace a voxel with the maximum over an ellipsoidal neighborhood of voxels.
If neighbours
of an axis is 1, no processing is done on that axis.
Examples:
1292 def magnitude(self): 1293 """Colapses components with magnitude function.""" 1294 imgm = vtk.new("ImageMagnitude") 1295 imgm.SetInputData(self.dataset) 1296 imgm.Update() 1297 self._update(imgm.GetOutput()) 1298 self.pipeline = utils.OperationNode("magnitude", parents=[self], c="#4cc9f0") 1299 return self
Colapses components with magnitude function.
1301 def topoints(self): 1302 """ 1303 Extract all image voxels as points. 1304 This function takes an input `Volume` and creates an `Mesh` 1305 that contains the points and the point attributes. 1306 1307 Examples: 1308 - [vol2points.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/vol2points.py) 1309 """ 1310 v2p = vtk.new("ImageToPoints") 1311 v2p.SetInputData(self.dataset) 1312 v2p.Update() 1313 mpts = vedo.Points(v2p.GetOutput()) 1314 mpts.pipeline = utils.OperationNode("topoints", parents=[self], c="#4cc9f0:#e9c46a") 1315 return mpts
Extract all image voxels as points.
This function takes an input Volume
and creates an Mesh
that contains the points and the point attributes.
Examples:
1317 def euclidean_distance(self, anisotropy=False, max_distance=None): 1318 """ 1319 Implementation of the Euclidean DT (Distance Transform) using Saito's algorithm. 1320 The distance map produced contains the square of the Euclidean distance values. 1321 The algorithm has a O(n^(D+1)) complexity over n x n x...x n images in D dimensions. 1322 1323 Check out also: https://en.wikipedia.org/wiki/Distance_transform 1324 1325 Arguments: 1326 anisotropy : bool 1327 used to define whether Spacing should be used in the 1328 computation of the distances. 1329 max_distance : bool 1330 any distance bigger than max_distance will not be 1331 computed but set to this specified value instead. 1332 1333 Examples: 1334 - [euclidian_dist.py](https://github.com/marcomusy/vedo/tree/master/examples/volumetric/euclidian_dist.py) 1335 """ 1336 euv = vtk.new("ImageEuclideanDistance") 1337 euv.SetInputData(self.dataset) 1338 euv.SetConsiderAnisotropy(anisotropy) 1339 if max_distance is not None: 1340 euv.InitializeOn() 1341 euv.SetMaximumDistance(max_distance) 1342 euv.SetAlgorithmToSaito() 1343 euv.Update() 1344 vol = Volume(euv.GetOutput()) 1345 vol.pipeline = utils.OperationNode("euclidean_distance", parents=[self], c="#4cc9f0") 1346 return vol
Implementation of the Euclidean DT (Distance Transform) using Saito's algorithm. The distance map produced contains the square of the Euclidean distance values. The algorithm has a O(n^(D+1)) complexity over n x n x...x n images in D dimensions.
Check out also: https://en.wikipedia.org/wiki/Distance_transform
Arguments:
- anisotropy : bool used to define whether Spacing should be used in the computation of the distances.
- max_distance : bool any distance bigger than max_distance will not be computed but set to this specified value instead.
Examples:
1348 def correlation_with(self, vol2, dim=2): 1349 """ 1350 Find the correlation between two volumetric data sets. 1351 Keyword `dim` determines whether the correlation will be 3D, 2D or 1D. 1352 The default is a 2D Correlation. 1353 1354 The output size will match the size of the first input. 1355 The second input is considered the correlation kernel. 1356 """ 1357 imc = vtk.new("ImageCorrelation") 1358 imc.SetInput1Data(self.dataset) 1359 imc.SetInput2Data(vol2.dataset) 1360 imc.SetDimensionality(dim) 1361 imc.Update() 1362 vol = Volume(imc.GetOutput()) 1363 1364 vol.pipeline = utils.OperationNode("correlation_with", parents=[self, vol2], c="#4cc9f0") 1365 return vol
Find the correlation between two volumetric data sets.
Keyword dim
determines whether the correlation will be 3D, 2D or 1D.
The default is a 2D Correlation.
The output size will match the size of the first input. The second input is considered the correlation kernel.
1367 def scale_voxels(self, scale=1): 1368 """Scale the voxel content by factor `scale`.""" 1369 rsl = vtk.new("ImageReslice") 1370 rsl.SetInputData(self.dataset) 1371 rsl.SetScalarScale(scale) 1372 rsl.Update() 1373 self._update(rsl.GetOutput()) 1374 self.pipeline = utils.OperationNode( 1375 "scale_voxels", comment=f"scale={scale}", parents=[self], c="#4cc9f0" 1376 ) 1377 return self
Scale the voxel content by factor scale
.