1. parse xml data

2. stack tiles
3. 8bpp and 32bpp
4. extract values from specific layer named 'other'
5. test 8bpp against corresponding 32bpp

imagej_tiff.py

 #!/usr/bin/env python3
 
-'''
-/**
- * @file imagej_tiff.py
- * @brief open multi layer tiff files, display layers and parse meta data
- * @par <b>License</b>:
- *  This program is free software: you can redistribute it and/or modify
- *  it under the terms of the GNU General Public License as published by
- *  the Free Software Foundation, either version 3 of the License, or
- *  (at your option) any later version.
- *
- *  This program is distributed in the hope that it will be useful,
- *  but WITHOUT ANY WARRANTY; without even the implied warranty of
- *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- *  GNU General Public License for more details.
- *
- *  You should have received a copy of the GNU General Public License
- *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
-*/
-'''
-
-__copyright__ = "Copyright 2018, Elphel, Inc."
-__license__   = "GPL-3.0+"
-__email__     = "oleg@elphel.com"
-
 from PIL import Image
 import xml.etree.ElementTree as ET
 import numpy as np
 import matplotlib.pyplot as plt
 
+import sys
+import xml.dom.minidom as minidom
+
+import time
+
+#http://stackoverflow.com/questions/287871/print-in-terminal-with-colors-using-python
+class bcolors:
+    HEADER = '\033[95m'
+    OKBLUE = '\033[94m'
+    OKGREEN = '\033[92m'
+    WARNING = '\033[38;5;214m'
+    FAIL = '\033[91m'
+    ENDC = '\033[0m'
+    BOLD = '\033[1m'
+    BOLDWHITE = '\033[1;37m'
+    UNDERLINE = '\033[4m'
+
+# reshape to tiles
+def get_tile_images(image, width=8, height=8):
+  _nrows, _ncols, depth = image.shape
+  _size = image.size
+  _strides = image.strides
+
+  nrows, _m = divmod(_nrows, height)
+  ncols, _n = divmod(_ncols, width)
+  if _m != 0 or _n != 0:
+    return None
+
+  return np.lib.stride_tricks.as_strided(
+    np.ravel(image),
+    shape=(nrows, ncols, height, width, depth),
+    strides=(height * _strides[0], width * _strides[1], *_strides),
+    writeable=False
+  )
+
 # TiffFile has no len exception
 #import imageio
 
@@ -67,6 +78,12 @@ class imagej_tiff:
     self.labels = []
     # infos will contain xml data Elphel stores in some of tiff files
     self.infos = []
+    # dictionary from decoded infos[0] xml data
+    self.props = {}
+
+    # bits per sample, type int
+    self.bpp = tif.tag[258][0]
+
     self.__split_labels(tif.n_frames,tif.tag)
     self.__parse_info()
     # image layers stacked along depth - (think RGB)
@@ -77,7 +94,21 @@ class imagej_tiff:
       tif.seek(i)
       a = np.array(tif)
       a = np.reshape(a,(a.shape[0],a.shape[1],1))
+
       #a = a[:,:,np.newaxis]
+
+      # scale for 8-bits
+      # exclude layer named 'other'
+      if self.bpp==8:
+        _min = self.data_min
+        _max = self.data_max
+        _MIN = 1
+        _MAX = 255
+        a = a.astype(float)
+        if self.labels[i]!='other':
+          a[a==0]=np.nan
+          a = (_max-_min)*(a-_MIN)/(_MAX-_MIN)+_min
+
       # init
       if i==0:
         self.image = a
@@ -88,12 +119,37 @@ class imagej_tiff:
     # init done, close the image
     tif.close()
 
+  # label == tiff layer name
+  def getvalues(self,label=""):
+    l = self.getstack([label],shape_as_tiles=True)
+    res = np.empty((l.shape[0],l.shape[1],3))
+
+    for i in range(res.shape[0]):
+      for j in range(res.shape[1]):
+        # 9x9 -> 81x1
+        m = np.ravel(l[i,j])
+        if self.bpp==32:
+          res[i,j,0] = m[0]
+          res[i,j,1] = m[2]
+          res[i,j,2] = m[4]
+        elif self.bpp==8:
+          res[i,j,0] = ((m[0]-128)*256+m[1])/128
+          res[i,j,1] = ((m[2]-128)*256+m[3])/128
+          res[i,j,2] = (m[4]*256+m[5])/65536.0
+        else:
+          res[i,j,0] = np.nan
+          res[i,j,1] = np.nan
+          res[i,j,2] = np.nan
+
+    return res
+
+
   # get ordered stack of images by provided items
   # by index or label name
-  def getstack(self,items=[]):
+  def getstack(self,items=[],shape_as_tiles=False):
     a = ()
     if len(items)==0:
-      return self.image
+      b = self.image
     else:
       for i in items:
         if type(i)==int:
@@ -102,8 +158,12 @@ class imagej_tiff:
           j = self.labels.index(i)
           a += (self.image[:,:,j],)
       # stack along depth
-      return np.stack(a,axis=2)
+      b = np.stack(a,axis=2)
+
+    if shape_as_tiles:
+      b = get_tile_images(b,self.tileW,self.tileH)
 
+    return b
 
   # get np.array of a channel
   # * do not handle out of bounds
@@ -152,7 +212,7 @@ class imagej_tiff:
     #tmp_im.show()
 
 
-  # puts etrees in infos
+  # puts etrees in infoss
   def __parse_info(self):
 
     infos = []
@@ -161,6 +221,21 @@ class imagej_tiff:
 
     self.infos = infos
 
+    # specifics
+    # properties dictionary
+    pd = {}
+
+    for child in infos[0]:
+      #print(child.tag+"::::::"+child.text)
+      pd[child.tag] = child.text
+
+    self.props = pd
+
+    # tiles are squares
+    self.tileW = int(self.props['tileWidth'])
+    self.tileH = int(self.props['tileWidth'])
+    self.data_min = float(self.props['data_min'])
+    self.data_max = float(self.props['data_max'])
 
   # makes arrays of labels (strings) and unparsed xml infos
   def __split_labels(self,n,tag):
@@ -188,24 +263,98 @@ class imagej_tiff:
         skip -= 1
       tag_labels = tag_labels[l:]
 
+
 #MAIN
 if __name__ == "__main__":
+
+  try:
+    fname = sys.argv[1]
+  except IndexError:
+    fname = "1521849031_093189-ML_DATA-32B-O-OFFS1.0.tiff"
+    fname = "1521849031_093189-ML_DATA-08B-O-OFFS1.0.tiff"
+
   #fname = "1521849031_093189-DISP_MAP-D0.0-46.tif"
-  fname = "test.tiff"
+  #fname = "1526905735_662795-ML_DATA-08B-AIOTD-OFFS2.0.tiff"
+  #fname = "test.tiff"
+
+  print(bcolors.BOLDWHITE+"time: "+str(time.time())+bcolors.ENDC)
 
   ijt = imagej_tiff(fname)
 
+  print(bcolors.BOLDWHITE+"time: "+str(time.time())+bcolors.ENDC)
+
   print(ijt.labels)
   print(ijt.infos)
+
+  rough_string = ET.tostring(ijt.infos[0], "utf-8")
+  reparsed = minidom.parseString(rough_string)
+  print(reparsed.toprettyxml(indent="\t"))
+
+  print(ijt.props)
+
+  # needed properties:
+  print("Tiles shape: "+str(ijt.tileW)+"x"+str(ijt.tileH))
+  print("Data min: "+str(ijt.data_min))
+  print("Data max: "+str(ijt.data_max))
+
   print(ijt.image.shape)
 
-  #ijt.show_images()
+  # layer order: ['diagm-pair', 'diago-pair', 'hor-pairs', 'vert-pairs', 'other']
+  # now split this into tiles:
+
+  #tiles = get_tile_images(ijt.image,ijt.tileW,ijt.tileH)
+  #print(tiles.shape)
+
+  tiles = ijt.getstack(['diagm-pair','diago-pair','hor-pairs','vert-pairs'],shape_as_tiles=True)
+  print("Stack of images shape: "+str(tiles.shape))
+
+  print(bcolors.BOLDWHITE+"time: "+str(time.time())+bcolors.ENDC)
+  # provide layer name
+  values = ijt.getvalues(label='other')
+  print("Stack of values shape: "+str(values.shape))
+
+  print(bcolors.BOLDWHITE+"time: "+str(time.time())+bcolors.ENDC)
+  #print(values)
+
+  #print(value_tiles[131,162].flatten())
+  #print(np.ravel(value_tiles[131,162]))
+
+  #values = np.empty((vt.shape[0],vt.shape[1],3))
+
+  #for i in range(values.shape[0]):
+  #  for j in range(values.shape[1]):
+  #    values[i,j,0] = get_v1()
+
+
+  #print(tiles[121,160,:,:,0].shape)
+  #_nrows = int(ijt.image.shape[0] / ijt.tileH)
+  #_ncols = int(ijt.image.shape[1] / ijt.tileW)
+  #_nrows = 32
+  #_ncols = 32
+  #print(str(_nrows)+" "+str(_ncols))
+  #fig, ax = plt.subplots(nrows=_nrows, ncols=_ncols)
+  #for i in range(_nrows):
+  #  for j in range(_ncols):
+  #    ax[i,j].imshow(tiles[i+100,j,:,:,0])
+  #    ax[i,j].set_axis_off()
+
+  #for i in range(5):
+  #  fig = plt.figure()
+  #  plt.imshow(tiles[121,160,:,:,i])
+  #  plt.colorbar()
+
+  #ijt.show_images(['other'])
 
   #ijt.show_images([0,3])
-  ijt.show_images(['X-corr','Y-corr'])
+  #ijt.show_images(['X-corr','Y-corr'])
   #ijt.show_images(['R-vign',3])
 
-  plt.show()
+  #ijt.show_images()
+  #plt.show()
+
+
+
+
 
   # Examples
 

test_8bpp_vs_32bpp.py

+#!/usr/bin/env python3
+
+from PIL import Image
+import xml.etree.ElementTree as ET
+import numpy as np
+import matplotlib.pyplot as plt
+
+import sys
+import xml.dom.minidom as minidom
+
+import time
+
+from imagej_tiff import imagej_tiff
+
+fname08bpp = "1521849031_093189-ML_DATA-08B-O-OFFS1.0.tiff"
+fname32bpp = "1521849031_093189-ML_DATA-32B-O-OFFS1.0.tiff"
+
+im08 = imagej_tiff(fname08bpp)
+im32 = imagej_tiff(fname32bpp)
+
+print(im08.props)
+
+THRESHOLD = 0.01
+THRESHOLD2 = 0.1
+
+values08 = im08.getvalues(label='other')
+values32 = im32.getvalues(label='other')
+
+# compare values
+
+print("test1: Compare values layers")
+
+values32[np.isnan(values32)] = 0
+values08[np.isnan(values08)] = 0
+values08[values08==-256.0] = 0
+
+delta = abs(values32 - values08)
+#print(delta)
+#print(delta.shape)
+
+error_counter = 0
+for i in range(delta.shape[0]):
+  for j in range(delta.shape[1]):
+    for k in range(delta.shape[2]):
+      if delta[i,j,k]>THRESHOLD:
+        error_counter += 1
+        print("Fail ("+str(error_counter)+"): "+str((i,j,k))+" "+str(delta[i,j,k]))
+        print(str(values32[i,j,k])+" vs "+str(values08[i,j,k]))
+
+if error_counter==0:
+  print("Values are OK")
+
+#v08 = im08.getstack(['other'],shape_as_tiles=True)
+#v32 = im32.getstack(['other'],shape_as_tiles=True)
+
+#print(v08.shape)
+
+#_nrows = 1
+#_ncols = 1
+
+#fig = plt.figure()
+#plt.imshow(v08[241,321,:,:,0])
+#plt.colorbar()
+
+#fig = plt.figure()
+#plt.imshow(v32[241,321,:,:,0])
+#plt.colorbar()
+
+# compare layers
+
+print("test2: Compare layers")
+
+tiles08 = im08.getstack(['diagm-pair','diago-pair','hor-pairs','vert-pairs'],shape_as_tiles=True)
+tiles32 = im32.getstack(['diagm-pair','diago-pair','hor-pairs','vert-pairs'],shape_as_tiles=True)
+
+delta2 = abs(tiles32-tiles08)
+error_counter = 0
+for i in range(delta2.shape[0]):
+  for j in range(delta2.shape[1]):
+    for k in range(delta2.shape[2]):
+      for l in range(delta2.shape[3]):
+        for m in range(delta2.shape[4]):
+          if delta2[i,j,k,l,m]>THRESHOLD2:
+            error_counter += 1
+            print("Fail ("+str(error_counter)+"): "+str((i,j,k,l,m))+" "+str(delta2[i,j,k,l,m]))
+            print(str(tiles32[i,j,k,l,m])+" vs "+str(tiles08[i,j,k,l,m]))
+
+if error_counter==0:
+  print("Layers are OK")
+
+#fig = plt.figure()
+#plt.imshow(tiles08[121,187,:,:,1])
+#plt.colorbar()
+
+#fig = plt.figure()
+#plt.imshow(tiles32[121,187,:,:,1])
+#plt.colorbar()
+
+plt.show()
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+