Rap method implementation

pyproject.toml

@@ -1,3 +1,23 @@
 [build-system]
 requires = ["setuptools"]
 build-backend = "setuptools.build_meta"
+
+[project]
+name = "quactography"
+version = "0.1"
+description = "Quactography : quantum tracking"
+authors = [{ name = "SCIL Team" }]
+readme = "README.md"
+requires-python = ">=3.10, <3.13"
+license-files = ["LICENSE"]
+dynamic = ["dependencies", "scripts"]
+classifiers = [
+    "Development Status :: 3 - Alpha",
+    "Environment :: Console",
+    "Intended Audience :: Science/Research",
+    "Operating System :: OS Independent",
+    "Programming Language :: Python",
+    "Topic :: Scientific/Engineering"
+]
+
+

quactography/adj_matrix/filter.py

@@ -160,4 +160,4 @@ def remove_zero_columns_rows(mat: np.ndarray):
     zero_rows = np.all(mat == 0, axis=1)
     non_zero_cols = np.where(~zero_cols)[0]
     non_zero_rows = np.where(~zero_rows)[0]
-    return mat[np.ix_(non_zero_rows, non_zero_cols)]
+    return mat[np.ix_(non_zero_rows, non_zero_cols)]

quactography/adj_matrix/reconst.py

@@ -1,6 +1,8 @@
 import numpy as np
 from dipy.reconst.shm import sh_to_sf
 from dipy.core.sphere import Sphere
+from quactography.graph.utils import get_output_nodes
+
 
 
 def build_adjacency_matrix(nodes_mask):
@@ -16,6 +18,7 @@ def build_adjacency_matrix(nodes_mask):
     adj_matrix : np.ndarray
         Adjacency matrix of the nodes.
     """
+
     # 1st. Assign labels to non-zero voxels (nodes)
     node_indices = np.flatnonzero(nodes_mask)
 
@@ -47,7 +50,18 @@ def build_adjacency_matrix(nodes_mask):
                             x + x_offset, y + y_offset, z + z_offset, i, nodes_mask, labels_volume, adj_matrix
                         )
 
-    return adj_matrix, node_indices
+        # Adds possibility of an edge to the actual node closest neighbour in 26 directions
+        for x_offset in [-1, 0, 1]:
+            for y_offset in [-1, 0, 1]:
+                for z_offset in [-1, 0, 1]:
+                    if x_offset == 0 and y_offset == 0 and z_offset == 0:
+                        continue
+                    else:
+                        adj_matrix = _add_edge_perhaps(
+                            x + x_offset, y + y_offset, z + z_offset, i, nodes_mask, labels_volume, adj_matrix
+                        )
+
+    return adj_matrix, node_indices, labels_volume
 
 
 def build_weighted_graph(adj_matrix, node_indices, sh, sh_order=12):
@@ -135,6 +149,33 @@ def _add_edge_perhaps(
     return adj_matrix
 
 
+def add_end_point_edge(adj_matrix, end, labels):
+    """
+    Add a node and edges to the end points of the ROI in the adjacency matrix.
+
+    Parameters
+    ----------
+    adj_matrix : np.ndarray
+        Adjacency matrix of the graph.
+    end_points : list of tuples
+        List of end points as (x, y, z) coordinates.
+    labels_volume : np.ndarray
+        List of every column in the image:
+
+    Returns
+    -------
+    np.ndarray
+        Updated adjacency matrix with end point edges added.
+    """
+
+    adj_matrix = np.lib.pad(adj_matrix, (0, 1), 'constant', constant_values=(0))
+    for i in end:
+        start = labels[i[0], i[1], i[2]]
+        adj_matrix[start, -1] = 1
+        adj_matrix[-1, start] = 1  # Assuming undirected graph
+    return adj_matrix
+
+
 def _is_valid_pos(pos_x, pos_y, pos_z, nodes_mask):
     if pos_x < 0 or pos_x >= nodes_mask.shape[0]:
         return False
@@ -143,3 +184,16 @@ def _is_valid_pos(pos_x, pos_y, pos_z, nodes_mask):
     if pos_z < 0 or pos_z >= nodes_mask.shape[2]:
         return False
     return nodes_mask[pos_x, pos_y, pos_z]
+
+if __name__ == "__main__":
+    # little demo code here and proto test 
+    mask = np.zeros((5, 5, 5), dtype=bool)
+    mask[1:3, 1:3, 1:3] = True
+    mat, nodes, labels = build_adjacency_matrix(mask)
+    entry_node = np.array([1, 1, 1])
+    propagation_direction = np.array([0, 1, 1])
+    angle_rad = np.pi / 8
+
+    indices = get_output_nodes(mask, entry_node, propagation_direction, angle_rad)
+    new = add_end_point_edge(mat,indices,labels)
+    print(new)

quactography/graph/utils.py

@@ -48,3 +48,15 @@ def get_output_nodes(mask, entry_node, propagation_direction, angle_rad):
     is_valid_edge = directions_cos_angle > propagation_cos_angle
     valid_indices = edge_indices[is_valid_edge.reshape((-1,))]
     return valid_indices
+
+
+if __name__ == '__main__':
+    # little demo code here
+    mask = np.zeros((25, 25, 25), dtype=bool)
+    mask[5:10, 5:10, 5:10] = True
+    entry_node = np.array([5, 5, 5])
+    propagation_direction = np.array([0, 1, 1])
+    angle_rad = np.pi / 8
+
+    indices = get_output_nodes(mask, entry_node, propagation_direction, angle_rad)
+    print(indices)

quactography/hamiltonian/hamiltonian_qubit_edge.py

@@ -328,22 +328,20 @@ def get_exact_sol(self):
 
 # # # # TEST:--------------------------------------------------------------------------
 
+# # from quactography.hamiltonian.hamiltonian_qubit_node import Hamiltonian_qubit_node
+# import numpy as np
 
-# # # mat = np.array([[0, 1, 1, 0], [1, 0, 0, 5], [1, 0, 0, 6], [0, 5, 6, 0]])
+# from quactography.adj_matrix.io import save_graph
+# mat = np.array([[0, 1, 1, 0], [1, 0, 0, 5], [1, 0, 0, 6], [0, 5, 6, 0]])
 
 # # # # # This is the given format you should use to save the graph, for mat:
-# # # # save_graph(mat, np.array([0, 1, 2, 3]), np.array([4, 4]), "rand_graph.npz")
+# save_graph(mat, np.array([0, 1, 2, 3]), np.array([4, 4]), "rand_graph.npz")
 # import sys
 
-# sys.path.append(r"C:\Users\harsh\quactography")
-
 # from quactography.graph.undirected_graph import Graph
 # from quactography.adj_matrix.io import load_graph
 
-# # from quactography.hamiltonian.hamiltonian_qubit_node import Hamiltonian_qubit_node
-# import numpy as np
 
-# from quactography.adj_matrix.io import save_graph
 
 # my_graph_class = Graph(
 #     np.array(
@@ -407,14 +405,14 @@ def get_exact_sol(self):
 # print_hamiltonian_circuit(h.hint_c, ["10101"])
 
 
-# # # print("total2")
-# # # print_hamiltonian_circuit(h.total_hamiltonian, ["11111"])
-# # # print("mandatory2")
-# # # print_hamiltonian_circuit(h.mandatory_c, ["11111"])
-# # # print("start2")
-# # # print_hamiltonian_circuit(h.starting_node_c, ["11111"])
-# # # print("finish2")
-# # # print_hamiltonian_circuit(h.ending_node_c, ["11111"])
-# # # print("int2")
-# # # print_hamiltonian_circuit(h.hint_c, ["11111"])
+# print("total2")
+# print_hamiltonian_circuit(h.total_hamiltonian, ["11111"])
+# print("mandatory2")
+# print_hamiltonian_circuit(h.mandatory_c, ["11111"])
+# print("start2")
+# print_hamiltonian_circuit(h.starting_node_c, ["11111"])
+# print("finish2")
+# print_hamiltonian_circuit(h.ending_node_c, ["11111"])
+# print("int2")
+# print_hamiltonian_circuit(h.hint_c, ["11111"])
 # # # ------------------------------------------------------------------------------------------------------------------------------------------

quactography/solver/Dijkstra.py

@@ -0,0 +1,74 @@
+import heapq
+import networkx as nx
+
+
+def dijkstra_stepwise(Graph, start, target):
+    """
+    Perform Dijkstra's algorithm in a stepwise manner on a given graph.
+    Parameters
+    ----------
+    G : networkx.Graph
+        The input graph where each node is typically a tuple (e.g., (x, y)).
+        Edges may have a 'weight' attribute indicating traversal cost
+        (default is 1).
+    start : hashable
+        The starting node for the pathfinding algorithm.
+    target : hashable
+        The target node to reach.
+    diagonal_mode : str, optional
+        Neighbor retrieval mode. Must be either:
+        - "nondiagonal": only 4-directional neighbors (up, down, left, right)
+        - "diagonal": include diagonal neighbors (8 directions total)
+        Default is "nondiagonal".
+    Returns
+    -------
+    evaluated_nodes : list
+        The list of nodes in the order they were evaluated by the algorithm.
+    path_to_current : list of list
+        A list containing the reconstructed path from the start node to each
+        node evaluated so far, in evaluation order.
+    current_distance : float
+        The final distance to the target node if found, or to the last
+        evaluated node otherwise. Returns None if no path is found.
+    """
+
+    G = nx.from_numpy_array(Graph)
+    distances = {node: float('inf') for node in G.nodes()}
+    distances[start] = 0
+    previous_nodes = {node: None for node in G.nodes()}
+    evaluated_nodes = []
+    path_to_current = []
+    priority_queue = [(0, start)]
+    heapq.heapify(priority_queue)
+
+    while priority_queue:
+        current_distance, current_node = heapq.heappop(priority_queue)
+
+        if current_node not in evaluated_nodes:
+            evaluated_nodes.append(current_node)
+
+        temp_path = []
+        node = current_node
+        while node is not None:
+            temp_path.append(node)
+            if node in previous_nodes:
+                node = previous_nodes[node]
+            else:
+                break
+        temp_path.reverse()
+        path_to_current.append(temp_path)
+
+        if current_node == target:
+            break
+
+        neighbors = list(G.neighbors(current_node))
+
+        for neighbor in neighbors:
+            if neighbor not in evaluated_nodes:
+                edge_weight = G[current_node][neighbor].get("weight", 1)
+                new_distance = current_distance + edge_weight
+                if new_distance < distances[neighbor]:
+                    distances[neighbor] = new_distance
+                    previous_nodes[neighbor] = current_node
+                    heapq.heappush(priority_queue, (new_distance, neighbor))
+    return path_to_current[-1]