<jupyter_start><jupyter_text>Surface Distances and Hausdorff Distances<jupyter_code>import vtk
import pandas as pd
import numpy as np
from pathlib import Path
import os<jupyter_output><empty_output><jupyter_text>Files to compare in pairThe 1st file listed in inputFiles_0 will compare with the 1st file in inputFiles_1.The 2nd file in inputFiles_0 will compare with the 2nd file in inputFiles_1 and so on.<jupyter_code>inputFiles_0 = ["phantom0.vtp","phantom0.vtp","phantom0.vtp","phantom0.vtp"]
inputFiles_1 = ["phantom0.vtp","phantom1.vtp",
                "phantom2.vtp","phantom3.vtp"]

## path for exporting vtp files with surface distance scalar mapping
pathToColorDistanceScalarVTP = "VTP-Result-Exported"

distanceMode = "point to cell"
#distanceMode = "point to point"
def surfaceDistance(points_A,points_B,mode):
    information = vtk.vtkInformation()
    inputA = points_A.GetPoints()
    inputB = points_B
    if mode == "point to point":
        locator_A = vtk.vtkKdTreePointLocator()
        locator_A.SetDataSet(points_A)
        locator_A.BuildLocator()
        locator_B = vtk.vtkKdTreePointLocator()
        locator_B.SetDataSet(points_B)
        locator_B.BuildLocator()
    elif mode == "point to cell":
        cellLocatorB = vtk.vtkCellLocator()
        cellLocatorB.SetDataSet(inputB)
        cellLocatorB.BuildLocator()
        cellId = vtk.mutable(0)
        subId = vtk.mutable(0)
        dist2 = vtk.mutable(0)
        cell = [0.0, 0.0, 0.0]
        closestPoint = np.zeros(3)
        
    distanceAToB = []
    for i in range (inputA.GetNumberOfPoints()):
        currentPoint = inputA.GetPoint(i)
        if mode == "point to point":
            closestPointId = locator_B.FindClosestPoint(currentPoint)
            closestPoint = inputB.GetPoint(closestPointId)
        elif mode == "point to cell":
            cellLocatorB.FindClosestPoint(currentPoint,closestPoint, cellId, subId, dist2)
    
        dist = np.sqrt(np.power(currentPoint[0] - closestPoint[0], 2) +
        np.power(currentPoint[1] - closestPoint[1], 2) +
        np.power(currentPoint[2] - closestPoint[2], 2))
        distanceAToB.append(dist)
    return(distanceAToB)
def distanceAnalysis(fileName0,fileName1,sf_0to1,sf_1to0):
        parameters = []
        sf_0to1_min = min(sf_0to1)
        sf_1to0_min = min(sf_1to0)
        sf_0to1_max = max(sf_0to1)
        sf_1to0_max = max(sf_1to0)
        sf_0to1_mean = np.mean(sf_0to1)
        sf_1to0_mean = np.mean(sf_1to0)
        #sf_0to1_25 = np.percentile(sf_0to1, 25)
        #sf_1to0_25 = np.percentile(sf_1to0, 25)
        sf_0to1_75 = np.percentile(sf_0to1, 75)
        sf_1to0_75 = np.percentile(sf_1to0, 75)
        sf_0to1_95 = np.percentile(sf_0to1, 95)
        sf_1to0_95 = np.percentile(sf_1to0, 95)
                
        parameters = [ 
                       sf_0to1_min,sf_1to0_min,
                       sf_0to1_max,sf_1to0_max,
                       sf_0to1_mean,sf_1to0_mean,
                       sf_0to1_95,sf_1to0_95
                     ]
        return parameters
def setScalar(polydata,distance):
    nPoints = polydata.GetNumberOfPoints()
    assert(nPoints == len(distance))
    # VTK has its own array format. Convert the input
    # array (z) to a vtkFloatArray.
    distances = vtk.vtkFloatArray()
    distances.SetName("Distance")
    distances.SetNumberOfValues(nPoints)
    for i in range(nPoints):
        distances.SetValue(i, distance[i])
    # Assign the scalar array.
    polydata.GetPointData().SetScalars(distances)
    
    return polydata<jupyter_output><empty_output><jupyter_text>Main Program<jupyter_code>inputData_0 =[]
inputData_1 =[]
# Read files into list
for x in range(len(inputFiles_0)):
    reader = vtk.vtkXMLPolyDataReader()
    reader.SetFileName(inputFiles_0[x])
    reader.Update()
    inputData_0.append([inputFiles_0[x],reader.GetOutput()])
for x in range(len(inputFiles_1)):
    reader = vtk.vtkXMLPolyDataReader()
    reader.SetFileName(inputFiles_1[x])
    reader.Update()
    inputData_1.append([inputFiles_1[x],reader.GetOutput()])

# Calculation Result
in_0 = inputData_0
in_1 = inputData_1
resultList = []
rawDistanceList_0 = []
rawDistanceList_1 = []
for x in range(len(in_0)):
    sf_0to1 = surfaceDistance(in_0[x][1],in_1[x][1],distanceMode)
    sf_1to0 = surfaceDistance(in_1[x][1],in_0[x][1],distanceMode)
    parameters = distanceAnalysis(in_0[x][1],in_1[x][1],sf_0to1,sf_1to0)
    rawDistanceList_0.append([in_0[x][0],in_1[x][0],sf_0to1])
    rawDistanceList_1.append([in_1[x][0],in_0[x][0],sf_1to0])
    resultList.append([in_0[x][0],in_1[x][0],
                   parameters[0],parameters[1],
                   parameters[2],parameters[3],
                   max([parameters[2],parameters[3]]),
                   parameters[4],parameters[5],
                   parameters[6],parameters[7]])
    print("Processed files - {} and {}".format(in_0[x][0],in_1[x][0]))
    
#export color scalar vtp bidirectional
Path(pathToColorDistanceScalarVTP).mkdir(parents=True, exist_ok=True)

in_0 = inputData_0
in_1 = inputData_1
exportedDataList = []
for x in range(len(in_0)):
    data = in_0[x][1]
    distance = rawDistanceList_0[x][2]
    basename_0 = os.path.basename(rawDistanceList_0[x][0]).split(".", 1)[0]
    basename_1 = os.path.basename(rawDistanceList_0[x][1]).split(".", 1)[0]
    filename =  '{}_to_{}'.format(basename_0,basename_1)
    exportdata = setScalar(data,distance)
    exportedDataList.append(exportdata)
    writer = vtk.vtkXMLPolyDataWriter()
    writer.SetInputData(exportdata)
    writer.SetFileName('{}/{}.vtp'.format(pathToColorDistanceScalarVTP,filename))
    writer.Update()
    
for x in range(len(in_1)):
    data = in_1[x][1]
    distance = rawDistanceList_1[x][2]
    basename_0 = os.path.basename(rawDistanceList_1[x][0]).split(".", 1)[0]
    basename_1 = os.path.basename(rawDistanceList_1[x][1]).split(".", 1)[0]
    filename =  '{}_to_{}'.format(basename_0,basename_1)
    exportdata = setScalar(data,distance)
    writer = vtk.vtkXMLPolyDataWriter()
    writer.SetInputData(exportdata)
    writer.SetFileName('{}/{}.vtp'.format(pathToColorDistanceScalarVTP,filename))
    writer.Update()

n = len(in_1) + len(in_1)
print("Total {} files have been exported with distance scalars in folder {}".format(n,pathToColorDistanceScalarVTP))<jupyter_output>Processed files - phantom0.vtp and phantom0.vtp
Processed files - phantom0.vtp and phantom1.vtp
Processed files - phantom0.vtp and phantom2.vtp
Processed files - phantom0.vtp and phantom3.vtp
Total 8 files have been exported with distance scalars in folder VTP-Result-Exported<jupyter_text>Result<jupyter_code>df = pd.DataFrame(resultList, columns =['File 1','File 2',
                                             'Min HF (File 1 to 2)','Min HF (File 2 to 1)',
                                             'Max HF (File 1 to 2)','Max HF (File 2 to 1)','Max HF (bi-directional)',
                                             'Mean HF (File 1 to 2)','Mean HF (File 2 to 1)',
                                             '95% HF (File 1 to 2)','95% HF (File 2 to 1)' ])
df.index += 1  
df_list = []
n = len(inputFiles_0) // 10
if n ==0:
    display(df)
else:
    for x in range(n):
        index0 = 10 * x
        index1 = index0 + 10
        df_list.append(df[index0:index1])
    if not df[index1:].empty:
        df_list.append(df[index1:])
    for i in df_list:
        display(i)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.12304582210242587, 'fraction_numerical': 0.03948787061994609, 'mean_word_length': 3.157422969187675, 'pattern_counts': {'":': 4, '<': 16, '<?xml version=': 0, '>': 16, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '4', 'dir': 'jupyter-structured-clean-dedup', 'id': '92147', 'license': 'MIT', 'n_tokens_mistral': 2714, 'n_tokens_neox': 2614, 'n_words': 474, 'path': '.ipynb_checkpoints/VTP_HausdorffDistanceReport-checkpoint.ipynb', 'repo_name': 'rolandshc/vtp-HausdorffDistance-metric', 'size': '7420'}