add classification, bad label exclusion and analyze deluxe

src/guanine_crystal_analysis/_widget.py

@@ -47,7 +47,10 @@ def analyze_image(image: "napari.types.ImageData",
     df_skimage["aspect_ratio"] = df_skimage["major_axis_length"]/df_skimage["minor_axis_length"]
     df_simpleitk = pd.DataFrame(label_statistics(image, label, size = True, intensity = True, perimeter = True, shape = True))
     results = pd.merge(df_skimage, df_simpleitk, on = "label", suffixes = ('_skimage', '_simpleitk'))
+    to_drop = ["extent", "local_centroid-0", "local_centroid-1","orientation"]
+    results = results.drop(to_drop, axis=1)
 
+    #visualize table in napari
     if napari_viewer is not None:
         from napari_workflows._workflow import _get_layer_from_data
         labels_layer = _get_layer_from_data(napari_viewer, label)
@@ -61,6 +64,141 @@ def analyze_image(image: "napari.types.ImageData",
         import pandas
         return pandas.DataFrame(results)
 
+def classify_objects(label: "napari.types.LabelsData",
+                    intensity: bool = False,
+                     shape: bool = True, 
+                     size: bool = True,
+                    classifier_filename:str = "object_classifier_20220523_",
+                    napari_viewer: "napari.Viewer" = None)-> "napari.types.LabelsData":
+    from napari_workflows._workflow import _get_layer_from_data
+    import pandas as pd
+    labels_layer = _get_layer_from_data(napari_viewer, label)
+    table = pd.DataFrame(labels_layer.properties)
+    print("alle keys:", table.keys())
+    print(type(table))
+    keep = table[['label', 'aspect_ratio', 'max_intensity', 'min_intensity', 'perimeter_skimage', 'area', 'mean_intensity', 'major_axis_length', 'minor_axis_length', 'circularity', 'solidity', 'eccentricity', 'roundness_skimage', 'median', 'sum', 'variance', 'perimeter_on_border','perimeter_on_border_ratio']]
+
+    #intensity table
+    df_intensity = keep[['label','max_intensity', 'mean_intensity', 'min_intensity','median', 'sum', 'variance']]
+
+    #size table
+    df_size = keep[['label','area']]
+
+    #shape table
+    df_shape = keep[['label', 'aspect_ratio','perimeter_skimage', 'major_axis_length',
+       'minor_axis_length', 'circularity', 'solidity', 'eccentricity',
+       'roundness_skimage', 'perimeter_on_border',
+       'perimeter_on_border_ratio']]
+
+    selected_table = None 
+
+    if size:
+        classifier_filename = classifier_filename + "size_"
+        if selected_table is not None:
+            selected_table = pd.merge(df_size, selected_table, on=('label'))
+        else:
+            selected_table = df_size
+
+    if shape:
+        classifier_filename = classifier_filename + "shape_"
+        if selected_table is not None:
+            selected_table = pd.merge(df_shape, selected_table, on=('label'))
+        else:
+            selected_table = df_shape
+
+    if intensity:
+        classifier_filename = classifier_filename + "intensity_"
+        if selected_table is not None:
+            selected_table = pd.merge(df_intensity, selected_table, on=('label'))
+        else:
+            selected_table = df_intensity
+
+    classifier_filename = classifier_filename[:-1]+".cl"
+    print("selected keys:", selected_table.keys())
+    print(classifier_filename)
+    return classifier(selected_table, label, classifier_filename)
+
+def classifier(table, label, tabrowcl_filename):
+    '''
+    this function applies the classifier to the testing data.
+    it creates:
+    - an output image
+    - writes the good crystal labels to a new table
+    '''
+
+    #import statements
+    from pyclesperanto_prototype import imshow, replace_intensities
+    import numpy as np
+    import pandas as pd
+    import apoc
+    import os
+    import pyclesperanto_prototype as cle
+    from skimage.segmentation import relabel_sequential
+    from pathlib import Path
+
+    #prediction
+    classifier_filename = Path(__file__).parent/"cl_classifier"/tabrowcl_filename
+    print(classifier_filename)
+
+    classifier = apoc.TableRowClassifier(opencl_filename= classifier_filename)
+    prediction = classifier.predict(table)
+
+    #append predicted classes as columns to the table
+    table['predicted_class'] = prediction
+
+    #add background to the prediction
+    predicted_class_with_background = [0] + prediction.tolist()
+
+    # connect prediction to label image
+    class_image = replace_intensities(label, predicted_class_with_background).astype(int)
+    return class_image
+
+
+def bad_label_exclusion(label: "napari.types.LabelsData",
+                        class_image: "napari.types.LabelsData")-> "napari.types.LabelsData":
+    import os
+    import numpy as np
+    from pyclesperanto_prototype import relabel_sequential
+
+    #change class of bad crystals (1) to background (0)
+    class_image = np.asarray(class_image)
+    class_image[class_image == 1]=0
+    class_image_mask = class_image.astype(bool)
+
+    #exclude bad labels from label image
+    label_image_filtered=np.copy(label)
+    label_image_filtered[class_image_mask==False]=0
+    return relabel_sequential(label_image_filtered)
+
+def analyze_deluxe(input_image: "napari.types.ImageData",
+                   normalize: bool = True,
+                   segmenter: str="object_segmenter_20220523_sorted.cl",
+                   minimum_pixel_count:int=50,
+                   intensity: bool = False,
+                    shape: bool = True, 
+                    size: bool = True,
+                    classifier_filename:str = "object_classifier_20220523_",
+                   napari_viewer: "napari.Viewer" = None):
+    if normalize:
+        image_maybe_normalized = normalization(input_image)
+        napari_viewer.add_image(image_maybe_normalized,name="normalized")
+    else:
+        image_maybe_normalized = input_image
+
+    segmentation_result = segmentation(image_maybe_normalized, segmenter, minimum_pixel_count)
+    napari_viewer.add_labels(segmentation_result)
+
+    analyze_image(image_maybe_normalized, segmentation_result, napari_viewer)
+    classification_result = classify_objects(segmentation_result, intensity, shape, size, classifier_filename, napari_viewer)
+    napari_viewer.add_labels(classification_result)
+
+    exclusion_result = bad_label_exclusion(segmentation_result, classification_result)
+    napari_viewer.add_labels(exclusion_result)
+
+    analyze_image(image_maybe_normalized, exclusion_result, napari_viewer)
+
+
+
 
 
 

src/guanine_crystal_analysis/cl_classifier/object_classifier_20220523_intensity.cl

@@ -0,0 +1,164 @@
+/*
+OpenCL RandomForestClassifier
+classifier_class_name = TableRowClassifier
+feature_specification = max_intensity mean_intensity min_intensity median sum variance
+num_ground_truth_dimensions = 1
+num_classes = 2
+num_features = 6
+max_depth = 2
+num_trees = 10
+apoc_version = 0.8.0
+*/
+__kernel void predict (IMAGE_in0_TYPE in0, IMAGE_in1_TYPE in1, IMAGE_in2_TYPE in2, IMAGE_in3_TYPE in3, IMAGE_in4_TYPE in4, IMAGE_in5_TYPE in5, IMAGE_out_TYPE out) {
+ sampler_t sampler = CLK_NORMALIZED_COORDS_FALSE | CLK_ADDRESS_CLAMP_TO_EDGE | CLK_FILTER_NEAREST;
+ const int x = get_global_id(0);
+ const int y = get_global_id(1);
+ const int z = get_global_id(2);
+ float i0 = READ_IMAGE(in0, sampler, POS_in0_INSTANCE(x,y,z,0)).x;
+ float i1 = READ_IMAGE(in1, sampler, POS_in1_INSTANCE(x,y,z,0)).x;
+ float i2 = READ_IMAGE(in2, sampler, POS_in2_INSTANCE(x,y,z,0)).x;
+ float i3 = READ_IMAGE(in3, sampler, POS_in3_INSTANCE(x,y,z,0)).x;
+ float i4 = READ_IMAGE(in4, sampler, POS_in4_INSTANCE(x,y,z,0)).x;
+ float i5 = READ_IMAGE(in5, sampler, POS_in5_INSTANCE(x,y,z,0)).x;
+ float s0=0;
+ float s1=0;
+if(i0<1.827868938446045){
+ s0+=25.0;
+} else {
+ if(i4<3159.510009765625){
+  s0+=24.0;
+  s1+=28.0;
+ } else {
+  s0+=22.0;
+  s1+=2.0;
+ }
+}
+if(i3<0.959505021572113){
+ if(i2<0.41803279519081116){
+  s0+=25.0;
+ } else {
+  s0+=3.0;
+  s1+=2.0;
+ }
+} else {
+ if(i4<2760.884765625){
+  s0+=13.0;
+  s1+=28.0;
+ } else {
+  s0+=21.0;
+  s1+=9.0;
+ }
+}
+if(i3<2.205306053161621){
+ if(i1<1.728459358215332){
+  s0+=38.0;
+  s1+=21.0;
+ } else {
+  s0+=23.0;
+  s1+=2.0;
+ }
+} else {
+ if(i1<4.788457870483398){
+  s1+=16.0;
+ } else {
+  s0+=1.0;
+ }
+}
+if(i3<2.673387050628662){
+ if(i4<3262.62939453125){
+  s0+=50.0;
+  s1+=26.0;
+ } else {
+  s0+=18.0;
+ }
+} else {
+ s1+=7.0;
+}
+if(i3<2.673387050628662){
+ if(i4<3258.214599609375){
+  s0+=46.0;
+  s1+=26.0;
+ } else {
+  s0+=20.0;
+  s1+=1.0;
+ }
+} else {
+ s1+=8.0;
+}
+if(i1<0.9806220531463623){
+ if(i4<271.11474609375){
+  s0+=23.0;
+ } else {
+  s0+=3.0;
+  s1+=1.0;
+ }
+} else {
+ if(i3<2.673387050628662){
+  s0+=43.0;
+  s1+=21.0;
+ } else {
+  s1+=10.0;
+ }
+}
+if(i5<0.11518747359514236){
+ if(i4<185.0370635986328){
+  s0+=21.0;
+ } else {
+  s0+=7.0;
+  s1+=4.0;
+ }
+} else {
+ if(i0<3.4861111640930176){
+  s0+=8.0;
+  s1+=11.0;
+ } else {
+  s0+=37.0;
+  s1+=13.0;
+ }
+}
+if(i3<2.2412109375){
+ if(i3<0.9633152484893799){
+  s0+=21.0;
+  s1+=1.0;
+ } else {
+  s0+=45.0;
+  s1+=26.0;
+ }
+} else {
+ s1+=8.0;
+}
+if(i1<1.0034186840057373){
+ if(i1<0.9806220531463623){
+  s0+=22.0;
+ } else {
+  s0+=3.0;
+  s1+=1.0;
+ }
+} else {
+ if(i5<0.42960822582244873){
+  s0+=9.0;
+  s1+=23.0;
+ } else {
+  s0+=29.0;
+  s1+=14.0;
+ }
+}
+if(i0<1.9183892011642456){
+ s0+=23.0;
+} else {
+ if(i2<0.18333333730697632){
+  s0+=27.0;
+  s1+=16.0;
+ } else {
+  s0+=6.0;
+  s1+=29.0;
+ }
+}
+ float max_s=s0;
+ int cls=1;
+ if (max_s < s1) {
+  max_s = s1;
+  cls=2;
+ }
+ WRITE_IMAGE (out, POS_out_INSTANCE(x,y,z,0), cls);
+}