Chen, Elvis C.S.

Nisar, Hareem

Peters, Terry M

English

Scholarship@Western

Western University Scholarship@Western Robarts Imaging Publications Robarts Research Institute 2021 CHARACTERIZATION & CALIBRATION OF FORESIGHT ICE Hareem Nisar Western University, hnisar3@uwo.ca Terry M Peters Imaging Research Laboratories, Robarts Research Institute, London, Ontario, Canada & Biomedical Engineering Graduate Program, Western University, London, Ontario, Canada & Department of Medical Biophysics, Western University, London, Ontario, Canada Elvis C.S. Chen Western University Follow this and additional works at: https://ir.lib.uwo.ca/robartspub  Part of the Bioimaging and Biomedical Optics Commons Citation of this paper: Nisar, Hareem; Peters, Terry M; and Chen, Elvis C.S., "CHARACTERIZATION & CALIBRATION OF FORESIGHT ICE" (2021). Robarts Imaging Publications. 27. https://ir.lib.uwo.ca/robartspub/27  CHARACTERIZATION & CALIBRATION OF FORESIGHTTM ICE Written by: Hareem Nisar Supervisor: Dr. Terry Peters Advisor: Dr. Elvis Chen   JULY 8, 2021 ROBARTS RESEARCH INSTITUTE Western University   Table of Contents                               Image geometry  θ θφ φImage Orientation Image Data Acquisition Image Display Size  𝑅𝑎𝑑𝑖𝑢𝑠 (𝑖𝑛 𝑝𝑖𝑥𝑒𝑙𝑠) = −0.06φ2  +  10.88φ −  56.92,φ𝑅𝑎𝑑𝑖𝑢𝑠 (𝑖𝑛 𝑝𝑖𝑥𝑒𝑙𝑠) = (−0.14φ2  +  25.31φ −  132.38) ∗ 103 ∗ (𝑅𝑚𝑎𝑥),φ• • • • First Scan-Line Consistency  • • • • • Imaging Angle Consistency ±Image Artifacts Volume Reconstruction & Rendering  Surface Texture mapping DICOM reconstruction and custom interpolation       Temporal calibration  Method:Results:   Spatial calibration Methods: φθ φValidation  Calibration Precision φφPoint Source Localization Spherical Phantom Centroid Results:  φφDiscussion Challenges Mitigation Characterization of ICE Calibration Repeatability of Spatial Calibration Intra-procedural positioning  Calibration Uniformity wrt Imaging angle TRE  (mm) Radial distance Tracking accuracy wrt Imaging Angle Appendix A  Code for Conical Volume Reconstruction From a 2D Image #                                                             # #                                                             # #------  Slicer function for Conical Image Reconstruction ----# #                                                             # #                                                             # import slicer.util import numpy as np import math    def reconstruction(self, inputVolume, outputVolume, anglePhi):   #both volumes are scalar.         phi = float(anglePhi)          imgArr = slicer.util.arrayFromVolume(inputVolume)     imgArr = np.squeeze(imgArr)     Row, Col = imgArr.shape          #Defining the origin/apex of the cone as the middle point of the image.      origin = ((int(Row/2), int(Col/2)))     edge   = ((Row, Col))          #Calculating the maximum value of z based on the max displacement of the image from the center      # i.e from the middle of the image (cone apex) to the corner of the image.     max_d = math.sqrt(sum( (edge - origin)**2 for edge, origin in zip(edge, origin)))     max_z = max_d*math.tan(math.radians(90-phi))          #making an empty volume (choose an appropriate value of z)     vol = np.zeros((int(max_z+1), Row, Col), dtype=np.uint16)          #reconstruction: keeping x-axis and y-axis as same,      #calculation a height for each pixel value and populating volume      #origin(at apex) stays where it is, peripherals moved along z-axis.     result = np.where(imgArr) #get index of all voxels, for label map replace imgArr==1     X = result[0] #row     Y = result[1] #col     D = np.sqrt((X-origin[0])**2 + (Y-origin[1])**2) # 6.6 us for X= np.zeros(200)      Z = D*math.tan(math.radians(90-phi))     Z = Z.astype(int)     vol[Z,X, Y] = imgArr[X,Y]           #Update the volume(data) in the volume node with new array      slicer.util.updateVolumeFromArray(outputVolume, vol)          # -- Optional Volume Rendering --- #     volumeNode = outputVolume     slicer.util.updateVolumeFromArray(volumeNode, vol)     volRenLogic = slicer.modules.volumerendering.logic()     volRenWidget = slicer.modules.volumerendering.widgetRepresentation()     displayNode = volRenLogic.CreateDefaultVolumeRenderingNodes(volumeNode)     volumePropertyNode = displayNode.GetVolumePropertyNode()     volumePropertyNodeWidget = slicer.util.findChild(volRenWidget, 'VolumePropertyNodeWidget')     volumePropertyNodeWidget.setMRMLVolumePropertyNode(volumePropertyNode)          return outputVolume Code for Texture Mapping an Image to a Cone Model #                                                             # #                                                             # #------  VTK-based, Slicer function for Texture Mapping ------# #                                                             # #                                                             # import vtk, math    def textureMap(self, inputVolumeNode, anglePhi):     dim = inputVolumeNode.GetImageData().GetDimensions()     radius = int(dim[0]/2)     height = radius*(math.tan(math.radians(90-anglePhi)))          print("Texture Mapping function")     print("Cone Angle is ", anglePhi, " Radius is ", radius, " and Height is ", height)          ##create a VTK coneSource / conical plane     self.coneSource = vtk.vtkConeSource()     self.coneSource.SetHeight(height)     self.coneSource.SetRadius(radius)     self.coneSource.SetResolution(100)          #Center of the vtkConeSource is halfway up the height axis.      #Actual origin should be at the apex of the cone/middle of circle     self.coneSource.SetCenter(radius,radius, height/2)      self.coneSource.CappingOff()     self.coneSource.SetAngle(anglePhi)          #This part defines the axis: the line drawn from the center of the base up towards the apex/vertex.     self.coneSource.SetDirection(0,0,-1) #-z axis     self.coneSource.Update()          #Set the input image as teh texture     self.texture = vtk.vtkTexture()     self.texture.SetInputConnection(inputVolumeNode.GetImageDataConnection())          #Set the cone source as the plane to be texture mapped     self.texturePlane = vtk.vtkTextureMapToPlane()     self.texturePlane.SetInputConnection(self.coneSource.GetOutputPort())      #By default, the plane is centered at the origin and perpendicular to the z-axis,      #with width and height of length 1 and resolutions set to 1.     #Match the origin to the corner off of the cone and set two points at the other corners     self.texturePlane.SetOrigin( 0, 0, 0)     self.texturePlane.SetPoint1(  2*radius, 0, 0)     self.texturePlane.SetPoint2( 0,  2*radius, 0)          self.model = None     if self.model is None:         ##creating a new model in Slicer and set display properties         self.model = slicer.modules.models.logic().AddModel(self.texturePlane.GetOutputPort())         self.model.SetName("Cone Model")         self.modelDisp = self.model.GetDisplayNode()         self.modelDisp.SetAmbient(1)         self.modelDisp.SetDiffuse(1)         self.modelDisp.SetSpecular(1)         self.modelDisp.SetPower(1)         self.modelDisp.SetFrontfaceCulling(0)         self.modelDisp.SetBackfaceCulling(0)          #Apply texture image to the plane (cone)     self.modelDisp.SetTextureImageDataConnection(inputVolumeNode.GetImageDataConnection())          return self.model Code for Volume Reconstruction from a DICOM file  #                                                             # #                                                             # #------ Slicer function for Image Reconstruction from DICOM --# #                                                             # #                                                             # import scipy from scipy import interpolate import pydicom import numpy as np import math    def dicomReconstruction(self, filePath, outputVolumeNode, interpFactor, interpEnabled):     #fname = 'C:/Users/haree/Desktop/2.25.110163955407302945943070993686554932643.570.dcm'     fname = filePath     ds = pydicom.dcmread(fname)     img = ds.pixel_array #numpy array with pixel data     [row, col] = img.shape  # [TH, R]     R   = img.shape[1] # or R = cols, (representing each echo along radius R)     TH  = img.shape[0] # or TH = rows (representing each theta)          # Reading DICOM tags and extracting data     phi       = ds[0x15, 0x1000].value #float      theta_en  = ds[0x15, 0x1004].value #list #encoded [3,8... 1021]     numFrames = ds[0x28, 0x0008].value # pydicom.valuerep.IS     rows      = ds[0x28, 0x0010].value #int (THETA)     cols      = ds[0x28, 0x0011].value #int (along RADIUS)     spacing   = ds[0x28, 0x0030].value #pydicom.multival.MultiValue     TH_spacing= float(ds.PixelSpacing[0]) #float     R_spacing = float(ds.PixelSpacing[1]) #float     theta = theta_en/1024*360 # decoded theta values in degrees      #--------------- interpolating ----------------------------------------------------     if interpEnabled:       print("DCM shape before interpolation:  ",np.shape(img), np.shape(theta), R, TH)       rows2add = int(interpFactor) - 1       [img, theta] = self.interpImage(img, theta, rows2add) #function to interpolate img data before reconstruction       R   = img.shape[1] # or R = cols       TH  = img.shape[0] # or TH = rows       print("DCM shape after interpolation:   ",np.shape(img), np.shape(theta), R, TH)     #-----------------------------------------------------------------------------------          #needed for spacing later     maxC, minC = 0, 2000          # Min size requirement of 3D volume is based on simple geometry (see diagram)     [volx, voly, volz] = [2*R*math.sin(math.radians(phi)), 2*R*math.sin(math.radians(phi)), R*math.cos(math.radians(phi))]     vol = np.zeros((int(volx), int(voly), int(volz)), dtype=np.uint16)          [ox, oy, oz] = [round(volx/2), round(voly/2), 0]          vector = np.arange(0, TH)     sample = np.arange(0, R)      for v in vector:       #getting polar angle theta for each vector       th = theta[v]       for s in sample:         r = s*R_spacing #radial position for each sample          #convert to cartesian          a = r*math.sin(math.radians(phi))*math.cos(math.radians(th))         b = r*math.sin(math.radians(phi))*math.sin(math.radians(th))         c = r*math.cos(math.radians(phi))                  #scaling by 80(optimum number) to avoid overalpping or maybe FOV?         a=int(round(a*80) + ox)         b=int(round(b*80) + oy)         c=int(round(c*80) + oz)                  # dicom image to volume pixel mapping         vol[a,b,c] = img[v, s]                  if maxC<c:           maxC = c         if minC>c:           minC = c              sp = R*R_spacing*10*math.cos(math.radians(phi))/(maxC-minC)          # update the volume(data) in the volume node with new array      slicer.util.updateVolumeFromArray(outputVolumeNode, vol)     outputVolumeNode.SetSpacing((sp,sp,sp)) # setting spacing using cone geometry to ensure the radial length remains the same     outputVolumeNode.SetOrigin((0,0,0))     # now perform volume rendering through Slicer module          return outputVolumeNode      def interpImage(self, img, th, numRows):     #     #-------------------INTERPOLATING THETA VALUES --------------------     #      #Getting theta information from the original data, duplicating the last row, and then upsampling theta array          TH = img.shape[0]     # th = decoded theta values in degrees      #2. Adding last row = first row (so that interpolation is on a closed loop)     th = np.append(th, [th[0]+360])     TH += 1  #adding one to the total number of vectors TH=th.shape[0]          # 3. Upsampling theta values     intrp_func = interpolate.interp1d(np.array(range(TH)), th, kind='linear', axis=0)     R2add = numRows #e.g. 15 ROWS or VECTORS TO BE ADDED IN BETWEEN TWO VECTORS     new_row = np.linspace(0, TH-1, R2add*(TH-1)+TH)     intrp_th = np.array(intrp_func(new_row))     TH = intrp_th.shape[0]      # making sure that all theta values are under 360 (interpolation at the end)     for i in range(TH-1-R2add , TH):#intrp_th[-(R2add+1):]:       if intrp_th[i]>=360:         intrp_th[i] -= 360       #end if     #end for            #     #--------------------- INTERPOLATING IMAGE ---------------------------     #     # Before using interpolation on the rectangle image, copying the first vector as an additional last vector     # because of cone geometry or loop - first and the last vector will coincide now.      # Copy last row to end     img = np.concatenate((img, np.expand_dims(img[1,:], axis=0)), axis=0)     TH = img.shape[0]          # Upsampling - there is a duplicated last row which we will ignore in reconstruction     intrp_func = interpolate.interp1d(np.array(range(TH)), img, kind='linear', axis=0)          # ROWS or VECTORS TO BE ADDED IN BETWEEN TWO VECTORS (must be same as before)     new_row = np.linspace(0, TH-1, R2add*(TH-1)+TH)     intrp_img = np.array(intrp_func(new_row))     TH = intrp_img.shape[0]          return intrp_img, intrp_th   

CHARACTERIZATION & CALIBRATION OF FORESIGHT ICE

https://ir.lib.uwo.ca/cgi/viewcontent.cgi?article=1027&context=robartspub

CHARACTERIZATION & CALIBRATION OF FORESIGHT ICE

Abstract

Similar works

Full text

Available Versions

Scholarship@Western