Sfoglia il codice sorgente

Added table height correction and more statistics. Bug fixes and optimization as well.

Luka 1 settimana fa
parent
commit
e9dd6a032e
1 ha cambiato i file con 62 aggiunte e 34 eliminazioni
  1. 62 34
      SeekTransformModule/SeekTransformModule.py

+ 62 - 34
SeekTransformModule/SeekTransformModule.py

@@ -443,6 +443,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             markups_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"Markers_{volume_name}")
             for centroid, intensity in unique_centroids:
                 markups_node.AddControlPoint(*centroid)
+                markups_node.SetDisplayVisibility(False)
                 #print(f"Detected Centroid (RAS): {centroid}, Max Intensity: {intensity}")
 
             return unique_centroids
@@ -491,28 +492,29 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             #mid_node.AddControlPoint(mid_ras)
 
             # Določi threshold glede na CBCT ali CT
-            threshold = -300 #if yesCbct else 0
+            threshold = -300 if yesCbct else -100 # če je cbct iščemo vrednost -300, sicer pri CT 0
 
             # Poišči rob mize (drugi dvig)
             previous_value = -1000
             edge_count = 0
             table_top_y = None
-            min_jump = 100
+            min_jump = 100 if yesCbct else 50  # Minimalni skok za CBCT in CT
 
             for y in range(len(column_values) - 1, -1, -1):  # Od spodaj navzgor
                 intensity = column_values[y]
                 #if column_values[y] > threshold and previous_value <= threshold:
                 if (intensity - previous_value) > min_jump and intensity > threshold:  # Namesto primerjave s threshold
                     if yesCbct:
-                        table_top_y = y
-                        print(f"Zgornji rob mize najden pri Y = {y}")
-                        break
-                    edge_count += 1
-                    print(f"Zaznan rob mize pri Y = {y}")
-                    if edge_count == 2:  # Drugi dvig = zgornji rob mize
-                        table_top_y = y
-                        print(f"Zgornji rob mize najden pri Y = {y}")
+                        table_top_y = y + 1 #Da nismo že v trupu
+                        #print(f"Zgornji rob mize najden pri Y = {table_top_y}") # CBCT
                         break
+                    if edge_count == 0 or (edge_count == 1 and previous_value < -200):  # Check if intensity is back lower than -400
+                        edge_count += 1
+                        #print(f"Zaznan rob mize pri X, Y, Z = {mid_x_voxel},{y}, {mid_z_voxel}")
+                        if edge_count == 2:  # Drugi dvig = zgornji rob mize
+                            table_top_y = y + 1
+                            #print(f"Zgornji rob mize najden pri Y = {table_top_y}")
+                            break
                 previous_value = column_values[y]
 
             if table_top_y is None:
@@ -526,22 +528,24 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             # Doda marker za višino mize
             table_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"VišinaMize_{ct_volume_name}")
             table_node.AddControlPoint(table_ras)
+            table_node.SetDisplayVisibility(False)
 
             # Izračun višine v mm
             image_center_y = dims[1] // 2
             pixel_offset = table_top_y - image_center_y
             mm_offset = pixel_offset * spacing[1]
 
-            print(f"📏 Miza je {abs(mm_offset):.2f} mm {'nižja' if mm_offset > 0 else 'višja'} od središča.")
-            print(f"📏 Miza je {abs(pixel_offset)} pixel {'nižja' if pixel_offset > 0 else 'višja'} od središča.")
+            #print(f"📏 Miza je {abs(mm_offset):.2f} mm {'nižja' if mm_offset > 0 else 'višja'} od središča.")
+            #print(f"📏 Miza je {abs(pixel_offset)} pixel {'nižja' if pixel_offset > 0 else 'višja'} od središča.")
             # Shrani v CSV
             if writefilecheck:
                 file_path = os.path.join(os.path.dirname(__file__), "heightdata.csv")
-                with open(file_path, mode='w', newline='') as file:
+                with open(file_path, mode='a', newline='') as file:
                     writer = csv.writer(file)
-                    writer.writerow([f"Upper part of table detected at Z = {mm_offset:.2f} mm, {pixel_offset} pixels"])
+                    modality = "CBCT " if yesCbct else "CT "
+                    writer.writerow([modality, ct_volume_name, f" Upper part of table detected at Z = {mm_offset:.2f} mm, {pixel_offset} pixels"])
 
-            return mm_offset
+            return mm_offset, pixel_offset
         
         
         # Globalni seznami za končno statistiko
@@ -561,11 +565,12 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             cbct_list = []
             ct_list = []
             volume_points_dict = {}
+            CT_offset = 0
 
             # Get child items of the study item
             volumeItems = vtk.vtkIdList()
             shNode.GetItemChildren(studyItem, volumeItems)
-
+            
             # Iteracija čez vse volumne v posameznem studyju
             for j in range(volumeItems.GetNumberOfIds()):
                 intermediateItem = volumeItems.GetId(j)
@@ -582,12 +587,7 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                         print("❌ This is an NRRD volume!")
                         continue  # Preskoči, če ni DICOM volume
                     
-                    if volumeNode and volumeNode.IsA("vtkMRMLScalarVolumeNode"):
-                        print(f"✔️ Found volume: {volumeNode.GetName()} (ID: {volumeItem})")
-                
-                    if not volumeNode or not volumeNode.IsA("vtkMRMLScalarVolumeNode"):
-                        print("Can't find volumeNode")
-                        #continue  # Preskoči, če ni veljaven volume
+                    
                         
                         
                     volumeName = volumeNode.GetName()
@@ -613,18 +613,46 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                     #manufacturer = volumeNode.GetAttribute("DICOM.Manufacturer")
                     #manufacturer = slicer.dicomDatabase.fileValue(uid, "0008,0070")
                     #print(manufacturer)
+                    
                     # Določimo, ali gre za CBCT ali CT
                     if "varian" in manufacturer.lower() or "elekta" in manufacturer.lower():
                         cbct_list.append(volumeName)
                         scan_type = "CBCT"
                         yesCbct = True
-                        #print("CBCT")
                     else:  # Siemens ali Philips
                         ct_list.append(volumeName)
                         scan_type = "CT"
                         yesCbct = False
-                        #print("CT")
 
+                    if volumeNode and volumeNode.IsA("vtkMRMLScalarVolumeNode"):
+                        print(f"✔️ {scan_type} {volumeNode.GetName()} (ID: {volumeItem})")
+                
+                    if not volumeNode or not volumeNode.IsA("vtkMRMLScalarVolumeNode"):
+                        print("Can't find volumeNode")
+                        #continue  # Preskoči, če ni veljaven volume
+                    
+                    mm_offset, pixel_offset = find_table_top_z(volumeName, writefilecheck, yesCbct)
+                        
+                    if scan_type == "CT":
+                        CT_offset = pixel_offset #Določi koliko je višina CT slike od središča
+                    else:
+                        # Poravnava CBCT z višino CT
+                        CBCT_offset = pixel_offset
+                        alignment_offset = CT_offset - CBCT_offset
+                        print(f"Poravnavam CBCT z CT. Offset: {alignment_offset}")
+
+                        # Uporabi alignment_offset za premik CBCT
+                        transform = vtk.vtkTransform()
+                        transform.Translate(0, alignment_offset, 0)  # Premik samo po y-osi (višina)
+                        
+                        #transformacija volumna oziroma poravnava po višini
+                        transformNode = slicer.vtkMRMLTransformNode()
+                        slicer.mrmlScene.AddNode(transformNode)
+                        transformNode.SetAndObserveTransformToParent(transform)
+                        volumeNode.SetAndObserveTransformNodeID(transformNode.GetID())  # Prilepi transformacijo na CBCT
+                        slicer.mrmlScene.RemoveNode(transformNode)  # Odstrani transformacijo iz scene
+                        
+                        
                     # Detekcija točk v volumnu
                     grouped_points = detect_points_region_growing(volumeName, yesCbct, intensity_threshold=3000)
                     #print(f"Populating volume_points_dict with key ('{scan_type}', '{volumeName}')")
@@ -634,9 +662,9 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             # Če imamo oba tipa volumna (CBCT in CT) **znotraj istega studyja**
             if cbct_list and ct_list:
                 ct_volume_name = ct_list[0]  # Uporabi prvi CT kot referenco
-                print("\nProcessing CT")
-                yesCbct = False
-                find_table_top_z(ct_volume_name, writefilecheck, yesCbct)        
+                print(f"\nProcessing CT: {ct_volume_name}")
+                #yesCbct = False
+                    
                 
                 ct_points = [centroid for centroid, _ in volume_points_dict[("CT", ct_volume_name)]]
 
@@ -645,10 +673,10 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                 else:
                     for cbct_volume_name in cbct_list:
                         print(f"\nProcessing CBCT Volume: {cbct_volume_name}")
-                        yesCbct = True
+                        #yesCbct = True
                         cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", cbct_volume_name)]]
                                                 
-                        find_table_top_z(cbct_volume_name, writefilecheck, yesCbct)
+                        #find_table_top_z(cbct_volume_name, writefilecheck, yesCbct)
 
                         if len(cbct_points) < 3:
                             print(f"CBCT Volume '{cbct_volume_name}' doesn't have enough points for registration.")
@@ -662,12 +690,12 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
 
                         ct_volume_node = slicer.util.getNode(ct_volume_name)
                         cbct_volume_node = slicer.util.getNode(cbct_volume_name)
-                        ct_spacing = ct_volume_node.GetSpacing()  # (x_spacing, y_spacing, z_spacing)
-                        cbct_spacing = cbct_volume_node.GetSpacing()  # (x_spacing, y_spacing, z_spacing)
+                        #ct_spacing = ct_volume_node.GetSpacing()  # (x_spacing, y_spacing, z_spacing)
+                        #cbct_spacing = cbct_volume_node.GetSpacing()  # (x_spacing, y_spacing, z_spacing)
                         
-                        ct_scale_factor = np.array(ct_spacing)  # Spacing za CT (x, y, z)                        
-                        cbct_scale_factor = np.array(cbct_spacing)  # Spacing za CBCT (x, y, z)
-                        print(ct_scale_factor, cbct_scale_factor)
+                        #ct_scale_factor = np.array(ct_spacing)  # Spacing za CT (x, y, z)                        
+                        #cbct_scale_factor = np.array(cbct_spacing)  # Spacing za CBCT (x, y, z)
+                        #print(ct_scale_factor, cbct_scale_factor)
 
                         
                         # Sortiramo točke po Z-koordinati (ali X/Y, če raje uporabljaš drugo os)