Merge pull request #39 from AllenInstitute/automate

Automate

Author: RussTorres

acpreprocessing/stitching_modules/stitch/graph_stitching_corr.py

@@ -0,0 +1,101 @@
+from cmath import phase
+from acpreprocessing.utils import io
+import argschema
+from argschema.fields import Str
+import os
+import matplotlib.pyplot as plt
+import numpy as np
+
+example_input = {
+    "stitchDir": "/ACdata/processed/iSPIM2/MN8_S10_220211_high_res/stitch-s3/iter0",
+    "filename": "pairwise-stitched.json",
+    "d_name": "MN8_S10_220211_high_res",
+    "save_dir": "/ACdata/processed/iSPIM2/MN8_S10_220211_high_res/stitch-s3/stitching_eval/"
+}
+
+def normalize(data):
+    norm_data = []
+    dmin, dmax = min(data), max(data)
+    for i, val in enumerate(data):
+        norm_data.append((val-dmin) / (dmax-dmin))
+    return norm_data
+
+class GraphStitchCorr(argschema.ArgSchema):
+    stitchDir = Str(required=True, description='stitch iter directory')
+    filename = Str(required=False, default="pairwise-stitched.json", description='which json file to grab stitch data from')
+    d_name = Str(required=True, description='dataset name')
+
+class GraphStitchCorr(argschema.ArgSchemaParser):
+    default_schema = GraphStitchCorr
+
+    def run(self):
+        pairwise = io.read_json(os.path.join(self.args['stitchDir'], self.args["filename"]))
+        n_pairs= len(pairwise)
+        labels = []
+        corrs = []
+        variances = []
+        displacement_z = []
+        # displacement_x = []
+        # displacement_y = []
+        # phase_corrs = []
+        for pair in pairwise:
+            pair_label="{}&{}".format(pair[0]["tilePair"]["tilePair"][0]["index"],pair[0]["tilePair"]["tilePair"][1]["index"])
+            labels.append(pair_label)
+            variances.append(pair[0]["variance"])
+            corrs.append(pair[0]["crossCorrelation"])
+            # displacement_x.append(abs(pair[0]["displacement"][0]))
+            # displacement_y.append(abs(pair[0]["displacement"][1]))
+            displacement_z.append(abs(pair[0]["displacement"][2]))
+            # phase_corrs.append(pair[0]["phaseCorrelation"])
+        
+
+        norm_var = normalize(variances)
+        # norm_displ_x = normalize(displacement_x)
+        # norm_displ_y = normalize(displacement_y)
+        norm_displ_z = normalize(displacement_z)
+        # norm_phase = normalize(phase_corrs)
+
+        plt.plot(labels, corrs, linestyle='--', marker='x', label="cross correlation")
+        plt.legend(loc="upper left")
+        plt.xticks(rotation = 90)
+        plt.xlabel("Tile Pair")
+        plt.ylabel("Cross Correlation")
+        plt.title("Cross correlation for {}".format(self.args['d_name']))
+        plt.savefig(self.args["save_dir"]+'crossCorr.png', bbox_inches = "tight")
+
+        plt.figure(2)
+        plt.plot(labels, corrs, linestyle='--', marker='x', label="cross correlation")
+        # plt.plot(labels, norm_displ_x, linestyle='--', marker='+', label="norm displacement x")
+        # plt.plot(labels, norm_displ_y, linestyle='--', marker='+', label="norm displacement y")
+        plt.plot(labels, norm_displ_z, linestyle='--', marker='+', label="norm displacement z")
+        plt.legend(loc="upper left")
+        plt.xticks(rotation = 90)
+        plt.xlabel("Tile Pair")
+        plt.ylabel("Values 0-1")
+        plt.title("Cross correlation and Z Normalized Displacements for {}".format(self.args['d_name']))
+        plt.savefig(self.args["save_dir"]+'with_z_displ.png', bbox_inches = "tight")
+
+        plt.figure(3)
+        plt.plot(labels, corrs, linestyle='--', marker='x', label="cross correlation")
+        plt.plot(labels, norm_var, linestyle='--', marker='.', label="norm variance")
+        plt.legend(loc="lower left")
+        plt.xticks(rotation = 90)
+        plt.xlabel("Tile Pair")
+        plt.ylabel("Values 0-1")
+        plt.title("Cross correlation and Normalized Variance for {}".format(self.args['d_name']))
+        plt.savefig(self.args["save_dir"]+'with_norm_vars.png', bbox_inches = "tight")
+        
+        # plt.figure(4)
+        # plt.plot(labels, corrs, linestyle='--', marker='x', label="cross correlation")
+        # plt.plot(labels, phase_corrs, linestyle='--', marker='.', label="phase correlation")
+        # plt.legend(loc="upper left")
+        # plt.xticks(rotation = 90)
+        # plt.xlabel("Tile Pair")
+        # plt.ylabel("Values 0-1")
+        # plt.title("Cross correlation and Phase Correlation for {}".format(self.args['d_name']))
+        # plt.savefig(self.args["save_dir"]+'with_norm_phase.png', bbox_inches = "tight")
+
+
+if __name__ == '__main__':
+    mod = GraphStitchCorr(example_input)
+    mod.run()

acpreprocessing/utils/8bit.py

@@ -0,0 +1,38 @@
+"""
+Basic script to convert existing n5 data to 8bit
+Used this for evaluation/experimentation of methods
+
+@author: shbhas
+"""
+
+import numpy as np
+import z5py
+
+output_n5 = "/ACdata/processed/MN7_RH_3_9_S33_220405_high_res/MN7_RH_3_9_S33_220405_high_res.n5/"
+out_8bit = "/ACdata/processed/MN7_RH_3_9_S33_220405_high_res8bit/MN7_RH_3_9_S33_220405_high_res8bit.n5"
+fout = z5py.File(out_8bit)
+
+for s in range(18,19):
+    with z5py.File(output_n5, mode='r') as f:
+        s0shape = f[f"setup{s}"]['timepoint0']["s5"].shape
+        t = 4
+        u = 9
+        xyf = s0shape[1]/t
+        zf = s0shape[0]/u
+        print(s0shape)
+        num_ints = np.iinfo(np.uint16).max + 1
+        lut = np.uint8(np.sqrt(np.arange(num_ints)))
+
+        fout.create_group(f"setup{s}")
+        fout[f"setup{s}"].create_group("timepoint0")
+        ds = fout[f"setup{s}"]["timepoint0"].create_dataset("s5", shape=s0shape, chunks = (int(zf),int(xyf),int(xyf)), dtype='uint8')
+        
+        for x in range(t):
+            for y in range(t):
+                for z in range(u):
+                    a = int(z*zf)
+                    bx = int(x*xyf)
+                    by = int(y*xyf)
+                    imvol = f[f"setup{s}"]['timepoint0']["s5"][a:int(a+zf), bx:int(bx+xyf), by:int(by+xyf)]
+                    fout[f"setup{s}"]["timepoint0"]["s5"][a:int(a+zf), bx:int(bx+xyf), by:int(by+xyf)] = lut[imvol]
+                    # print((a,int(a+zf),bx,int(bx+xyf),by,int(by+xyf)))

acpreprocessing/utils/metadata/parse_metadata.py

@@ -42,6 +42,17 @@ def get_number_of_positions(self):
         """Return number of position strips in section"""
         return len(self.md['positions'])
 
+    def get_number_of_channels(self):
+        """Return number of channels of data for section"""
+        return self.md['channels']
+
+    def get_direction(self):
+        """Get direction of y direction for position strips in section"""
+        if (self.md["positions"][0]["y_start_um"]) > (self.md["positions"][1]["y_start_um"]):
+            return True
+        else:
+            return False
+    
     def get_size(self):
         """Get tiff width, height, and number of frames in stack"""
         sz = [self.md['settings']['image_size_xy'][0],

acpreprocessing/utils/nglink/create_layer.py

@@ -6,7 +6,8 @@
     "outputDir": "/ACdata/processed/demoModules/output/",
     "position": 2,
     "rootDir": "/ACdata/processed/demoModules/raw/",
-    "md_filename": "acqinfo_metadata.json"
+    "md_filename": "acqinfo_metadata.json",
+    "resolution": "high_res"
     }
 
 
@@ -65,6 +66,7 @@ class CreateLayerSchema(argschema.ArgSchema):
     reverse = Boolean(required=False,default=False, description="Whether to reverse direction of stitching or not")
     deskew = Int(required=False,default=0, description="deskew factor (0 if want to leave undeskewed)")
     channel = Int(required=True, default=1, description="channel number")
+    resolution = Str(default="high_res", description="whether data is high_res or overview")
 
 class NgLayer(argschema.ArgSchemaParser):
     default_schema = CreateLayerSchema
@@ -79,8 +81,10 @@ def run(self, state=None):
         md = parse_metadata.ParseMetadata(input_data=md_input)
         pr = md.get_pixel_resolution()
         sz = md.get_size()
-        ypos = sz[1]-md.get_overlap()  # subtract height of image by pixel overlap to get yposition
-
+        if self.args["resolution"] == "high_res":
+            ypos = sz[1]-md.get_overlap()  # subtract height of image by pixel overlap to get yposition
+        elif self.args["resolution"] == "overview":
+            ypos = sz[1]
         if self.args["reverse"]:
             layer0 = create_layer(self.args['outputDir'], self.args['position'],
                                   -1*ypos, pr, self.args['deskew'], self.args['channel'])
@@ -100,7 +104,10 @@ def run_consolidate(self, state=None):
         pr = md.get_pixel_resolution()
         sz = md.get_size()
         n_pos = md.get_number_of_positions()
-        ypos = sz[1]-md.get_overlap()  # subtract height of image by pixel overlap to get yposition
+        if self.args["resolution"] == "high_res":
+            ypos = sz[1]-md.get_overlap()  # subtract height of image by pixel overlap to get yposition
+        elif self.args["resolution"] == "overview":
+            ypos = sz[1]
         # Create the layer
         if self.args["reverse"]:
             layer = create_layer(self.args['outputDir'], self.args['position'],

acpreprocessing/utils/nglink/create_overview.py

@@ -0,0 +1,68 @@
+from acpreprocessing.utils.nglink import create_layer, create_nglink
+from argschema.fields import Dict
+import argschema
+import os
+
+example_input = {
+    "run_input": {
+        "outputDir": "/ACdata/processed/demoModules/output/",
+        "rootDir": "/ACdata/processed/demoModules/raw/",
+        "ds_name": "MN6_2_S83_220531_high_res",
+        "mip_level": 3,
+        "md_filename": "/ACdata/processed/demoModules/raw/acqinfo_metadata.json",
+        "consolidate_pos": True,
+        "reverse_stitch": False,
+        "deskew": False
+        }
+    }
+
+
+class CreateOverviewSchema(argschema.ArgSchema):
+    run_input = argschema.fields.Dict(required=True, description='Input json for processing')
+
+class Overview(argschema.ArgSchemaParser):
+    default_schema = CreateOverviewSchema
+
+    def run(self, n_channels, n_pos, dirname, deskew, state=None):
+        if not state:
+            state = {"showDefaultAnnotations": False, "layers": []}
+        for channel in range(n_channels):
+            for pos in range(n_pos):
+                if self.args["run_input"]["consolidate_pos"]:
+                    layer_input = {
+                            "position": 0,
+                            "outputDir": self.args["run_input"]['outputDir']+dirname+".n5/channel"+str(channel),
+                            "rootDir": self.args["run_input"]['rootDir'],
+                            "reverse": self.args["run_input"]["reverse_stitch"],
+                            "deskew": deskew,
+                            "channel": channel,
+                            "resolution": dirname.split("_")[-1]
+                            }
+                    create_layer.NgLayer(input_data=layer_input).run_consolidate(state)
+                    break
+                else:
+                    layer_input = {
+                        "position": pos,
+                        "outputDir": self.args["run_input"]['outputDir']+dirname+".n5/channel"+str(channel),
+                        "rootDir": self.args["run_input"]['rootDir'],
+                        "reverse": self.args["run_input"]["reverse_stitch"],
+                        "deskew": deskew,
+                        "channel": channel,
+                        "resolution": dirname.split("_")[-1]
+                        }
+                    create_layer.NgLayer(input_data=layer_input).run(state)
+
+        # Create nglink from created state and estimated positions (overview link)
+        nglink_input = {
+                "outputDir": self.args["run_input"]['outputDir'],
+                "fname": self.args["run_input"]["nglink_name"],
+                "state_json": self.args["run_input"]["state_json"]
+                }
+        if not os.path.exists(os.path.join(nglink_input['outputDir'], nglink_input['fname'])):
+            create_nglink.Nglink(input_data=nglink_input).run(state)
+        else:
+            print("nglink txt file already exists!")
+
+
+if __name__ == '__main__':
+    mod = Overview(input_data=example_input).run()

acpreprocessing/utils/split_channels.py

@@ -0,0 +1,48 @@
+"""
+Basic script to split channels of acquisition (on tiff level)
+Used this before channel splitting was incorporated into n5 conversion
+
+@author: shbhas
+"""
+
+import os
+from natsort import natsorted, ns
+from acpreprocessing.utils import io
+from acpreprocessing.utils.metadata import parse_metadata
+
+def sort_files(filedir):
+    filelist = os.listdir(filedir)
+    return natsorted(filelist, alg=ns.IGNORECASE)
+
+outputdir = "/ACdata/processed/iSPIM2/MN7_RH_3_9_S19_220606_high_res/split_channels/"
+inputdir = "/ispim2_data/workflow_data/iSPIM2/MN7_RH_3_9_S19_220606_high_res/"
+md_input = {
+            "rootDir": "/ispim2_data/workflow_data/iSPIM2/MN7_RH_3_9_S19_220606_high_res/",
+            "fname": "acqinfo_metadata.json"
+            }
+
+n_pos = parse_metadata.ParseMetadata(input_data=md_input).get_number_of_positions()
+print(n_pos)
+if not os.path.isdir(outputdir):
+    os.makedirs(outputdir)
+
+for s in range(n_pos):
+    index = 0
+    posdir = inputdir+f"high_res_Pos{s}/"
+    for inputfile in sort_files(posdir):
+        print(inputfile)
+        if inputfile[0]=='.':
+            continue
+        I = io.get_tiff_image(posdir+inputfile)
+        if index%2==0:
+            chan0 = I[0::2, : , :]
+            chan1 = I[1::2, : , :]
+        else:
+            chan1 = I[0::2, : , :]
+            chan0 = I[1::2, : , :]
+        fname0 = outputdir + f"ch0/high_res_Pos{s}/"+"{0:05d}.tif".format(index)
+        fname1 = outputdir + f"ch1/high_res_Pos{s}/"+"{0:05d}.tif".format(index)
+        print(fname0)
+        io.save_tiff_image(chan0, fname0)
+        io.save_tiff_image(chan1, fname1)
+        index += 1