Added check for table height. TODO: first connect CT and CBCT table height then transform based on markers.

SeekTransformModule/SeekTransformModule.py

@@ -446,7 +446,102 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                 #print(f"Detected Centroid (RAS): {centroid}, Max Intensity: {intensity}")
 
             return unique_centroids
-    
+
+        def find_table_top_z(ct_volume_name, writefilecheck, yesCbct):
+            """
+            Najde višino zgornjega roba mize v CT/CBCT volumnu in doda markerje.
+            
+            :param ct_volume_name: Ime volumna v slicerju
+            :param writefilecheck: Če je True, zapiše rezultat v CSV
+            :param yesCbct: Če je True, uporabi CBCT thresholde
+            :return: Višina zgornjega roba mize v mm
+            """
+            # Pridobi volumen
+            ct_volume_node = slicer.util.getNode(ct_volume_name)
+            image_data = ct_volume_node.GetImageData()
+            spacing = ct_volume_node.GetSpacing()
+            dims = image_data.GetDimensions()
+            #origin = ct_volume_node.GetOrigin()
+
+            # Pretvori volumen v numpy array
+            np_array = slicer.util.arrayFromVolume(ct_volume_node)
+
+            # Izračunaj sredinske IJK koordinate
+            mid_ijk = [dims[0] // 2, dims[1] // 2, dims[2] // 2]
+
+            # Preveri, da so IJK indeksi v mejah volumna
+            mid_ijk = [max(0, min(dims[i] - 1, mid_ijk[i])) for i in range(3)]
+
+            # Pretvorba IJK → RAS
+            ijkToRasMatrix = vtk.vtkMatrix4x4()
+            ct_volume_node.GetIJKToRASMatrix(ijkToRasMatrix)
+
+            #mid_ras = np.array(ijkToRasMatrix.MultiplyPoint([*mid_ijk, 1]))[:3]
+
+            # Sredinski Z slice
+            mid_z_voxel = mid_ijk[2]
+            slice_data = np_array[mid_z_voxel, :, :]  # (Y, X)
+
+            # Sredinski stolpec
+            mid_x_voxel = mid_ijk[0] - 15  # 15 pikslov levo od sredine da ne merimo pri hrbtenici
+            column_values = slice_data[:, mid_x_voxel]  # Y smer
+
+            # Doda marker v RAS koordinatah
+            #mid_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"Sredina_{ct_volume_name}")
+            #mid_node.AddControlPoint(mid_ras)
+
+            # Določi threshold glede na CBCT ali CT
+            threshold = -300 #if yesCbct else 0
+
+            # Poišči rob mize (drugi dvig)
+            previous_value = -1000
+            edge_count = 0
+            table_top_y = None
+            min_jump = 100
+
+            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}")
+                        break
+                previous_value = column_values[y]
+
+            if table_top_y is None:
+                print("❌ Zgornji rob mize ni bil najden!")
+                return None
+
+            # Pretvorba Y IJK → RAS
+            table_ijk = [mid_x_voxel, table_top_y, mid_z_voxel]
+            table_ras = np.array(ijkToRasMatrix.MultiplyPoint([*table_ijk, 1]))[:3]
+
+            # Doda marker za višino mize
+            table_node = slicer.mrmlScene.AddNewNodeByClass("vtkMRMLMarkupsFiducialNode", f"VišinaMize_{ct_volume_name}")
+            table_node.AddControlPoint(table_ras)
+
+            # 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.")
+            # 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:
+                    writer = csv.writer(file)
+                    writer.writerow([f"Upper part of table detected at Z = {mm_offset:.2f} mm, {pixel_offset} pixels"])
+
+            return mm_offset
         
         
         # Globalni seznami za končno statistiko
@@ -474,10 +569,6 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
             # Iteracija čez vse volumne v posameznem studyju
             for j in range(volumeItems.GetNumberOfIds()):
                 intermediateItem = volumeItems.GetId(j)
-    
-                # Preveri, ali je to dejanska skupina volumnov (npr. "No study description")
-                intermediateName = shNode.GetItemName(intermediateItem)
-                #print(f"Checking intermediate item: {intermediateName}")
 
                 finalVolumeItems = vtk.vtkIdList()
                 shNode.GetItemChildren(intermediateItem, finalVolumeItems)  # Išči globlje!
@@ -496,30 +587,15 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                 
                     if not volumeNode or not volumeNode.IsA("vtkMRMLScalarVolumeNode"):
                         print("Can't find volumeNode")
-                        continue  # Preskoči, če ni veljaven volume
-
-                    # Preveri, če volume ima StorageNode (drugače `.GetFileName()` vrže napako)
-                    storageNode = volumeNode.GetStorageNode()
-                    
-                    
-                    if not storageNode:
-                        print("Can't find storageNode")
-                        continue  # Preskoči, če volume nima shranjenih DICOM podatkov
+                        #continue  # Preskoči, če ni veljaven volume
+                        
+                        
                     volumeName = volumeNode.GetName()
-                    #print(volumeName)
                     imageItem = shNode.GetItemByDataNode(volumeNode)
-                    #print(imageItem)
-                    #dicomUIDsList = volumeNode.GetAttribute("DICOM.instanceUIDs").split()
-                    # Preverimo modaliteto volumna (DICOM metapodatki)
-                    #modality = slicer.dicomDatabase.fileValue(storageNode.GetFileName(), "0008,0060")              #prazen   
-                    #modality = volumeNode.GetAttribute("DICOM.Modality")                                           #None
-                    #modality = slicer.dicomDatabase.fileValue(uid, "0008,0060")  # Modality                        #prazen
-                    #modality = slicer.dicomDatabase.fileValue(dicomUIDsList[0], "0008,0060")                       #prazen
                     modality = shNode.GetItemAttribute(imageItem, "DICOM.Modality")                                 #deluje!
-                    #print(modality)
 
-                    dimensions = volumeNode.GetImageData().GetDimensions()
-                    spacing = volumeNode.GetSpacing()
+                    #dimensions = volumeNode.GetImageData().GetDimensions()
+                    #spacing = volumeNode.GetSpacing()
                     #print(f"Volume {volumeNode.GetName()} - Dimenzije: {dimensions}, Spacing: {spacing}")
 
                     if modality != "CT":
@@ -542,12 +618,12 @@ class MyTransformModuleLogic(ScriptedLoadableModuleLogic):
                         cbct_list.append(volumeName)
                         scan_type = "CBCT"
                         yesCbct = True
-                        print("CBCT")
+                        #print("CBCT")
                     else:  # Siemens ali Philips
                         ct_list.append(volumeName)
                         scan_type = "CT"
                         yesCbct = False
-                        print("CT")
+                        #print("CT")
 
                     # Detekcija točk v volumnu
                     grouped_points = detect_points_region_growing(volumeName, yesCbct, intensity_threshold=3000)
@@ -558,15 +634,22 @@ 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)        
+                
                 ct_points = [centroid for centroid, _ in volume_points_dict[("CT", ct_volume_name)]]
 
                 if len(ct_points) < 3:
                     print(f"CT volume {ct_volume_name} doesn't have enough points for registration.")
                 else:
                     for cbct_volume_name in cbct_list:
+                        print(f"\nProcessing CBCT Volume: {cbct_volume_name}")
+                        yesCbct = True
                         cbct_points = [centroid for centroid, _ in volume_points_dict[("CBCT", cbct_volume_name)]]
+                                                
+                        find_table_top_z(cbct_volume_name, writefilecheck, yesCbct)
 
-                        print(f"\nProcessing CBCT Volume: {cbct_volume_name}")
                         if len(cbct_points) < 3:
                             print(f"CBCT Volume '{cbct_volume_name}' doesn't have enough points for registration.")
                             continue