update vicinity_connect and add cleaning of disconnected cilia

R/detectCilia.R

@@ -333,7 +333,8 @@ detectCilia <- function(
   # Use half of square root of a small cilium area as max distance to combine
   # two cilium points
   if(is.null(vicinity_connect)){
-    vicinity_connect <- floor(0.5 + sqrt(min_cilium_area_in_pixels))
+    # vicinity_connect <- floor(0.5 + sqrt(min_cilium_area_in_pixels))
+    vicinity_connect <- floor(0.5 * sqrt(min_cilium_area_in_pixels))
   }
 
   # Determine the nuclei mask area
@@ -615,7 +616,7 @@ detectCilia <- function(
     # threshold_connect <- quantile(Image_cilia_layer_connect,
     #                               (1-ratio_of_cilia_pixels), na.rm = TRUE)
     # threshold_connect <- as.numeric(threshold_connect)
-    threshold_connect <- 0.1*threshold_find #0.5 0.2
+    threshold_connect <- 0.1*threshold_find #0.5 0.2 # 0.1
 
     print(paste("The new threshold values are: threshold_find = ",
                 threshold_find, " and threshold_connect = ",
@@ -1608,6 +1609,187 @@ detectCilia <- function(
       !df_cilium_information$ciliumNumber %in% cilium_numbers_to_be_removed,]
   }
 
+  ###new from here
+  # Delete all cilia that are disconnected
+  # In every found cilium (cluster), delete separated cilium parts that
+  # are not directly connected (being vicinity_connect apart) to the
+  # brightest part of that cilium
+
+  df_cilium_all$disconnectedCilia <- FALSE
+  df_cilium_all$ClusterNumber <- 0
+
+  for(i in unique(df_cilium_all$ciliumNumber)){
+    # print(i)
+    df_dummy <- df_cilium_all[df_cilium_all$ciliumNumber == i,]
+
+    # # Create a data frame with (x, y) coordinates
+    # coordinates <- data.frame(x = df_dummy$pos_x, y = df_dummy$pos_y)
+    # 
+    # # Calculate pairwise distances
+    # distances <- as.matrix(dist(coordinates, method = "manhattan"))
+    # distances[!lower.tri(distances)] <- NA
+    # 
+    # # Set the threshold distance (e.g., 2 units)
+    # threshold_distance <- 1
+    # 
+    # # Initialize labels
+    # labels <- rep(NA, length(df_dummy$pos_x))
+    # current_label <- 1
+    # 
+    # # Assign labels based on connectivity
+    # for (k in 1:(length(x_coords) - 1)) {
+    #   if(is.na(labels[k])) {
+    #     neighbors <- which(distances[,k] <= threshold_distance)
+    #     labels[neighbors] <- current_label
+    #     current_label <- current_label + 1
+    #   }
+    # }
+    # 
+    # 
+    # 
+    # 
+    # 
+    # 
+    # distances <- as.matrix(dist(cbind(df_dummy$pos_x, df_dummy$pos_y), upper = FALSE, method = "manhattan"))
+    # distances[!lower.tri(distances)] <- NA
+    # 
+    # threshold_distance <- 1
+    # 
+    # connected_pairs <- which(distances <= threshold_distance, arr.ind = TRUE)
+    # 
+    # for (k in 1:nrow(connected_pairs)) {
+    #   point1 <- connected_pairs[k, "row"]
+    #   point2 <- connected_pairs[k, "col"]
+    #   
+    #   df_dummy
+    # }
+
+    j <- 1
+
+    rows_gone_through <- j
+    rows_to_go_through <- 1:length(df_dummy$pos_x)
+    rows_to_go_through <- rows_to_go_through[rows_to_go_through != j]
+
+    threshold_distance <- 1
+    df_dummy$ClusterNumber[j] <- 1
+
+    if(length(df_dummy$ClusterNumber) > 1){
+
+      next_point <- rows_to_go_through[
+        rows_to_go_through %in%
+          which(df_dummy$pos_x %in% (df_dummy$pos_x[j] - threshold_distance) : (df_dummy$pos_x[j] + threshold_distance) &
+                  df_dummy$pos_y %in% (df_dummy$pos_y[j] - threshold_distance) : (df_dummy$pos_y[j] + threshold_distance))][1]
+      if(is.na(next_point)){
+        j <- rows_to_go_through[1]
+      }else{
+        j <- next_point
+      }
+
+      while(0 %in% df_dummy$ClusterNumber){
+
+        .pos_x_distance <- df_dummy$pos_x[j] -
+          df_dummy$pos_x
+
+        .pos_y_distance <- df_dummy$pos_y[j] -
+          df_dummy$pos_y
+
+        .pos_x_distance[abs(.pos_x_distance) <= threshold_distance] <- 0
+        .pos_y_distance[abs(.pos_y_distance) <= threshold_distance] <- 0
+
+        .distance <- abs(.pos_x_distance) + abs(.pos_y_distance)
+
+        # Get the cluster number (close cilium that has been detected)
+        ClusterNumber_dummy <-
+          unique(df_dummy$ClusterNumber[.distance == 0])
+
+        if(length(ClusterNumber_dummy) == 1 && ClusterNumber_dummy == 0){
+          # Advance cluster number number because there is no cluster close by
+          ClusterNumber <- max(df_dummy$ClusterNumber) + 1
+        }else{
+          # Points belong to already existing cluster
+          ClusterNumber <- ClusterNumber_dummy[!(ClusterNumber_dummy == 0)]
+        }
+
+        # if(length(ClusterNumber) > 1){
+        #   print(paste0("The following clusters of cilium ", i,
+        #                " are now one: ", paste0(ClusterNumber, collapse = ", ")))
+        #   for(k in 2:length(ClusterNumber)){
+        #     df_dummy$ClusterNumber[
+        #       df_dummy$ClusterNumber == ClusterNumber[k]] <-
+        #       ClusterNumber[1]
+        #   }
+        #   rm(k)
+        #   
+        # }else{
+          df_dummy$ClusterNumber[
+            df_dummy$ClusterNumber == 0 & .distance == 0] <- ClusterNumber
+        # }
+        # 
+          rows_to_go_through <- rows_to_go_through[rows_to_go_through != j]
+          next_point <- rows_to_go_through[
+            rows_to_go_through %in%
+              which(df_dummy$pos_x %in% (df_dummy$pos_x[j] - threshold_distance) : (df_dummy$pos_x[j] + threshold_distance) &
+                      df_dummy$pos_y %in% (df_dummy$pos_y[j] - threshold_distance) : (df_dummy$pos_y[j] + threshold_distance))][1]
+          if(is.na(next_point)){
+            j <- rows_to_go_through[1]
+          }else{
+            j <- next_point
+          }
+
+      }
+      rm(list = c("j", "ClusterNumber"))
+
+      # df_test <- df_dummy %>% 
+      #   dplyr::group_by(ClusterNumber) %>%
+      #   dplyr::summarise(meanIntensity = mean(.data$fluorescence_intensity),
+      #                    clusterSize = length(.data$fluorescence_intensity))
+      # 
+      # # Remove all clusters that are too small
+      # df_test <- df_test[!(df_test$clusterSize < min_cilium_area_in_pixels),]
+      # 
+      # cluster_number_with_heighest_mean_intensity <- df_test$ClusterNumber[
+      #   df_test$meanIntensity == max(df_test$meanIntensity)]
+      # 
+      # if(length(cluster_number_with_heighest_mean_intensity)>1){
+      #   print("Please check the determination of the cluster number.")
+      # }
+      # 
+      # df_dummy$disconnectedPart[
+      #   df_dummy$ClusterNumber != cluster_number_with_heighest_mean_intensity] <- TRUE
+      # 
+      # df_dummy <- df_dummy[!df_dummy$disconnectedPart,]
+      # df_dummy <- df_dummy[,!names(df_dummy)=="ClusterNumber"]
+
+      # # Keep only rows that have been found
+      # if(i == unique(df_cilium_points$ciliumNumber)[1]){
+      #   df_cilium_points2 <- df_dummy
+      # }else{
+      #   df_cilium_points2 <- rbind(df_cilium_points2, df_dummy)
+      # }
+
+      if(max(df_dummy$ClusterNumber) > 1){
+        df_cilium_all$disconnectedCilia[df_cilium_all$ciliumNumber == i] <- TRUE
+      }
+
+    }
+
+  }
+  rm(list = c("i", "df_dummy"))
+
+  # Remove cilia with disconnected cilia
+  disconnectedCilia <- unique(df_cilium_all$ciliumNumber[df_cilium_all$disconnectedCilia])
+
+  df_cilium_information <- df_cilium_information[!df_cilium_information$ciliumNumber %in% disconnectedCilia, ]
+  df_cilium_all <- df_cilium_all[!df_cilium_all$ciliumNumber %in% disconnectedCilia,]
+
+  # Drop extra columns
+  df_cilium_all <-
+    df_cilium_all[,!names(df_cilium_all)=="disconnectedCilia"]
+  df_cilium_all <-
+    df_cilium_all[,!names(df_cilium_all)=="ClusterNumber"]
+
+  ###new up to here
+
   # Renumber cilia
   .number <- 1
   df_cilium_information$ciliumNumber_NEW <- NA