topng.py

   1 import cv2
   2 import os
   3 from os.path import isdir, join
   4 from os import mkdir
   5 from timeit import timeit
   6 import pydicom
   7 from pprint import pprint
   8 import numpy as np
   9 import pathlib
  10 import json
  11
  12 from scipy import ndimage
  13 from scipy import misc
  14
  15
  16 from decimal import Decimal as d, ROUND_HALF_UP as rhu
  17
  18 # locals
  19 from helpers import drawcontours
  20
  21 np.set_printoptions(edgeitems=372, linewidth=1200)# hey numpy, take advantage of my large screen
  22
  23 PNG16_MAX = pow(2, 16) - 1 # here if thinks the heaviest weight bit is for transparency or something not in use with dicom imgs
  24 PNG8_MAX = pow(2, 8+1) - 1 # heaviest weight bit is 8 => 2**8, but dont forget the others: the reason of +1
  25
  26 INPUT_DIR = '../../Data/Images_anonymous/Case_0002/'
  27 OUT_DIR = './generated/'
  28
  29 CROP_SIZE = (45, 45) # (width, height)
  30 RED_COLOR = 100
  31
  32 def roundall(*l):
  33         return (int(round(e)) for e in l)
  34
  35 def ftrim(Mat):
  36         # ---------- FTRIM ----------
  37         # private func to trim the Matrix, but in one direction, vertically | horizontally
  38         # return the slice, don't affect the Matrix
  39         # y | : for (top to bottom), x -> : for (left to right)
  40         #   v
  41         # not my fault, numpy architecture
  42
  43         y1 = y2 = i = 0
  44
  45         while i < len(Mat):# search by top
  46                 if Mat[i].max() > 0:
  47                         y1 = i
  48                         break
  49                 i += 1
  50
  51         i = len(Mat) - 1
  52         while i >= 0:# search by bottom
  53                 if Mat[i].max() > 0:
  54                         y2 = i
  55                         break
  56                 i -= 1
  57         # print('y1, y2', y1, y2)
  58         return slice(y1, y2+1)# +1 to stop at y2
  59
  60 def trim(Mat):
  61         # most use who make implicit call of ftrim twice
  62         horizontal = ftrim(Mat)
  63         vertical = ftrim(Mat.T)
  64         # print('horizontal:vertical', horizontal, vertical)
  65         return Mat[horizontal, vertical]
  66
  67
  68 def getxy(file):
  69         """
  70                 {
  71                         'Image00001': [{
  72                                 'color': '#ff0000',
  73                                 'points': [
  74                                         {'x': 94.377, 'y': 137.39},
  75                                         {'x': 100.38, 'y': 139.55},
  76                                         {'x': 103.26, 'y': 142.67},
  77                                         {'x': 105.91, 'y': 147.95},
  78                                         {'x': 105.42, 'y': 152.76},
  79                                         {'x': 100.62, 'y': 156.84},
  80                                         {'x': 95.338, 'y': 159.96},
  81                                         {'x': 89.573, 'y': 158.52},
  82                                         {'x': 84.53, 'y': 153},
  83                                         {'x': 82.848, 'y': 149.15},
  84                                         {'x': 82.368, 'y': 142.91},
  85                                         {'x': 85.01, 'y': 138.11},
  86                                         {'x': 89.813, 'y': 137.39},
  87                                         {'x': 94.377, 'y': 137.39}
  88                                 ]
  89                         }]
  90                 }
  91                 return [(94.377, 137.39), ...]
  92         """
  93         with open(file) as jsonfile:
  94                 data = json.load(jsonfile)
  95                 # pprint(data)
  96
  97                 for imgXXX in data: pass  # get the value of key ~= "Image00001", cause it's dynamic
  98
  99                 for obj in data[imgXXX]: pass  # get the object that contains the points, cause it's a list
 100                 points = obj['points']
 101
 102                 # print("print, ", data)
 103                 # print("imgXXX, ", imgXXX)
 104                 # print("points, ", points)
 105                 # print("imgXXX, ", obj['points'])
 106                 tmp = [np.array( (round(pt['x']), round(pt['y'])) ).astype(np.int32) for pt in
 107                            points]  # extract x,y. {'x': 94.377, 'y': 137.39} => (94.377, 137.39)
 108                 return np.array(tmp, dtype=np.int32)
 109                 # return tmp
 110
 111 def minmax(file):
 112         r = getxy(file)
 113         if r is not None:
 114                 # print(r) # log
 115                 xmin, ymin = np.min(r, axis=0)
 116                 xmax, ymax = np.max(r, axis=0)
 117
 118                 # print('xmax, ymax', xmax, ymax)
 119                 # print('xmin, ymin', xmin, ymin)
 120
 121                 return roundall(xmin, ymin, xmax, ymax)
 122
 123 def crop(Mat, maskfile, size=None):
 124         """
 125                 size : if filled with a (45, 45), it will search the center of maskfile then crop by center Mat
 126         """
 127         xmin, ymin, xmax, ymax = minmax(maskfile)
 128         if size:
 129                 xcenter = (xmin + xmax) / 2
 130                 ycenter = (ymin + ymax) / 2
 131
 132                 # crop coords
 133                 ymin, xmin, ymax, xmax = roundall(xcenter - size[0], ycenter - size[0],
 134                                                                   xcenter + size[1], ycenter + size[1])
 135
 136         return Mat[xmin:xmax, ymin:ymax]
 137
 138 def contour(Mat, maskfiles, color, thickness=1):
 139         # ---------- CONTOUR POLYGON ----------
 140         # thickness = -1 => fill it
 141         #                       = 1 => just draw the contour with color
 142         #
 143         # return new Mat
 144         contours = [getxy(maskfile) for maskfile in maskfiles] # list of list of coords,
 145         # [
 146         #       [ (3,43), (3,4) ],
 147         #       [ (33,43) ]
 148         # ]
 149         # print("log: contours", contours)
 150         if contours is not None:
 151                 # print('type contours', type(contours))
 152                 # print('contours', contours)
 153                 msk = np.zeros(Mat.shape, dtype=np.int32) # init all the mask with True...
 154                 cv2.drawContours(msk, contours, -1, True, thickness) # affect Mat, fill the polygone with False
 155                 msk = msk.astype(np.bool)
 156                 # print('msk')
 157                 # print(msk)
 158                 # print("bool", timeit(lambda:msk, number=1000))
 159                 # print("normal", timeit(lambda:msk.astype(np.bool), number=1000))
 160                 # Mat = cv2.cvtColor(Mat, cv2.COLOR_RGB2GRAY) # apply gray
 161
 162                 # Mat = drawcontours(Mat, contours)
 163                 # pass
 164                 Mat[msk] = color # each True in msk means 0 in Mat. msk is like a selector
 165         return Mat
 166
 167 def mask(Mat, maskfile, color=0, thickness=-1):
 168         # ---------- MASK POLYGON ----------
 169         # thickness = -1 => fill it
 170         #                       = 1 => just draw the contour with color
 171         #
 172         # return new Mat
 173         contours = getxy(maskfile)
 174         # print("log: contours", contours)
 175         if contours is not None:
 176                 # print('type contours', type(contours))
 177                 # print('contours', contours)
 178                 msk = np.ones(Mat.shape, dtype=np.int32) # init all the mask with True...
 179                 cv2.drawContours(msk, [contours], -1, False, thickness) # affect msk, fill the polygone with False, => further, don't touch where is False
 180                 msk = msk.astype(np.bool)
 181                 # print('msk')
 182                 # print(msk)
 183                 # print("bool", timeit(lambda:msk, number=1000))
 184                 # print("normal", timeit(lambda:msk.astype(np.bool), number=1000))
 185                 # Mat = cv2.cvtColor(Mat, cv2.COLOR_RGB2GRAY) # apply gray
 186
 187                 # Mat = drawcontours(Mat, contours)
 188                 # pass
 189                 Mat[msk] = color # each True in msk means 0 in Mat. msk is like a selector
 190         return Mat
 191
 192
 193 def hollowmask(Mat, epifile, endofile, color=0, thickness=-1):
 194         # ---------- MASK POLYGON ----------
 195         # thickness = -1 => fill it
 196         #                       = 1 => just draw the contour with color
 197         #
 198         # return new Mat
 199         epicontours = getxy(epifile)
 200         endocontours = getxy(endofile)
 201         if epicontours is not None and endocontours is not None:
 202                 msk = np.ones(Mat.shape, dtype=np.int32) # init all the mask with True...
 203                 cv2.drawContours(msk, [epicontours], -1, False, thickness) # affect msk, fill the polygone with False, => further, don't touch where is False
 204                 cv2.drawContours(msk, [endocontours], -1, True, thickness) #                     fill the intern polygone with True, (the holow in the larger polygon), => further, color where is True with black for example
 205                 msk = msk.astype(np.bool)
 206
 207                 Mat[msk] = color # each True in msk means 0 in Mat. msk is like a selector
 208         return Mat
 209
 210
 211 def sqrmask1(Mat, maskfile):
 212         # ---------- SQUARE MASK 1st approach ----------
 213         # print( timeit( lambda:sqrmask1(img, epimask), number=1000 ) ) # 0.48110522600000005
 214         # return new Mat
 215         xmin, ymin, xmax, ymax = minmax(maskfile)
 216         # print("xmin, ymin, xmax, ymax", xmin, ymin, xmax, ymax)
 217
 218         i = 0
 219         while i < ymin:# search by top
 220                 Mat[i].fill(0)# paint the row in black
 221                 i += 1
 222
 223         i = len(Mat) - 1
 224         while i > ymax:# search by bottom
 225                 Mat[i].fill(0)# paint the row in black
 226                 i -= 1
 227
 228         i = 0
 229         while i < xmin:# search by top
 230                 Mat.T[i, ymin:ymax+1].fill(0) # paint the column (row of the Transpose) in black, ymin and ymax to optimize, cause, I previously painted a part
 231                 i += 1
 232
 233         i = len(Mat.T) - 1
 234         while i > xmax:# search by bottom
 235                 Mat.T[i, ymin:ymax+1].fill(0) # paint the column (row of the Transpose) in black, ymin and ymax to optimize, cause, I previously painted a part
 236                 i -= 1
 237
 238         return Mat
 239
 240 def sqrmask2(Mat, maskfile):
 241         # ---------- SQUARE MASK 2nd and best approach ----------
 242         # print( timeit( lambda:sqrmask2(img, epimask), number=1000 ) ) # 0.3097705270000001
 243         # return new Mat
 244         xmin, ymin, xmax, ymax = minmax(maskfile)
 245         # print("xmin, ymin, xmax, ymax", xmin, ymin, xmax, ymax)
 246
 247         msk = np.ones(Mat.shape, dtype=np.int32) # init all the mask with True...
 248         msk = msk.astype(np.bool)
 249
 250         msk[ymin:ymax, xmin:xmax] = False # for after, don't touch between min and max region
 251         Mat[msk] = 0
 252
 253         return Mat
 254
 255 def map16(array):# can be useful in future
 256         return array * 16
 257
 258 # def dround(n):
 259 #       return d(str(n)).quantize(d('1'), rounding=rhu)# safe round without
 260
 261 def affine(Mat, ab, cd):
 262         """
 263                 Affine transformation
 264                 ab: the 'from' interval (2, 384) or {begin:2, end:384} begin is 0, and end is 1
 265                 cd: the 'to' interval (0, 1024) {begin:0, end:1024}
 266         """
 267         a, b = ab[0], ab[1]
 268         c, d = cd[0], cd[1]
 269
 270         with np.nditer(Mat, op_flags=['readwrite']) as M:
 271                 for x in M:
 272                         x[...] = max( 0, round( (x-a) * (d-c) / (b-a) + c ) ) # could not be negative
 273
 274
 275 def getdir(filepath):
 276         return '/'.join(filepath.split('/')[:-1]) + '/'
 277
 278 def readpng(inputfile):# just a specific func to preview a "shape = (X,Y,3)" image
 279         image = misc.imread(inputfile)
 280
 281         print("image")
 282         print("image.shape", image.shape)
 283
 284         for tab in image:
 285                 for i, row in enumerate(tab):
 286                         print(row[0]*65536 + row[0]*256 + row[0], end=" " if i % image.shape[0] != 0 else "\n")
 287
 288 def topng(inputfile, outfile=None, overwrite=True, verbose=False, epimask='', endomask='', centercrop=None, blackmask=False, square=False, redcontour=False):
 289         """
 290                 (verbose) return (64, 64) : the width and height
 291                 (not verbose) return (img, dicimg) : the image and the dicimg objects
 292                 centercrop : it's a size (45, 45), to mention I want to crop by center of epimask and by size.
 293                 blackmask : draw the outside with black
 294
 295         """
 296         try:
 297                 dicimg = pydicom.read_file(inputfile) # read dicom image
 298         except pydicom.errors.InvalidDicomError as e:
 299                 # @TODO: log, i can't read this file
 300                 return
 301         # img = trim(dicimg.pixel_array)# get image array (12bits), Don't trim FOR THE MOMENT
 302         img = dicimg.pixel_array# get image array (12bits)
 303
 304         # test <<
 305         # return img.shape # $$ COMMENT OR REMOVE THIS LINE
 306         # print('img', img)
 307         # print('img', type(img))
 308         # print('min', img.min(), 'max', img.max())
 309         # dicimg.convert_pixel_data() # same as using dicimg.pixel_array
 310         # pixa = dicimg._pixel_array
 311         # print('dicimg._pixel_array', pixa)
 312         # print('dicimg.pixel_array==pixa', dicimg.pixel_array==pixa)
 313         # test >>
 314
 315         # affine transfo to png 16 bits, func affects img variable
 316         maxdepth = pow(2, dicimg.BitsStored) - 1
 317         # print('dicimg.BitsStored, (img.min(), img.max())', dicimg.BitsStored, (img.min(), img.max())) # testing..
 318         if int(dicimg.BitsStored) < 12:
 319                 print("\n\n\n-----{} Bits-----\n\n\n\n".format(dicimg.BitsStored))
 320         affine(img,
 321                 (0, maxdepth), # img.min() replace 0 may be, but not sure it would be good choice
 322                 (0, PNG16_MAX if maxdepth > PNG8_MAX else PNG8_MAX)
 323         )
 324
 325         # test <<
 326         # imgray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
 327         # ret,thresh = cv2.threshold(img,127,255,0)
 328         # im2, contours, hierarchy = cv2.findContours(thresh,cv2.RETR_TREE,cv2.CHAIN_APPROX_SIMPLE)
 329         # print("im2, contours, hierarchy", im2, contours, hierarchy)
 330         # test >>
 331
 332         # return img
 333
 334         # print("log: epimask", epimask)
 335
 336         if epimask:
 337                 if redcontour:
 338
 339                         contours = [epimask] # list of list of coords
 340                         if endomask:contours.append(endomask)
 341
 342                         img = contour(img, contours, color=RED_COLOR, thickness=1)
 343
 344                 if blackmask:# if is there a mask => apply
 345                         if square:
 346                                 img = sqrmask2(img, epimask)
 347                         else:
 348                                 if endomask:
 349                                         img = hollowmask(img, epimask, endomask)
 350                                 else:
 351                                         img = mask(img, epimask)
 352
 353                 if centercrop:
 354                         img = crop(img, epimask, centercrop)
 355
 356
 357         # return
 358         # test
 359         # if verbose:
 360         #       return img, dicimg, minmax(epimask)
 361         # else:
 362         #       return img.shape
 363
 364         savepath = (outfile or inputfile) + '.png'
 365         savedir = getdir(savepath)
 366         if overwrite and not isdir( savedir ):
 367                 pathlib.Path(savedir).mkdir(parents=True, exist_ok=True)
 368
 369
 370         # test <<
 371         # tmp = np.array(img) # to get eye on the numpy format of img
 372         # tmp = np.array(img) # to get eye on the numpy format of img
 373         # print("img[0,0]", img[0,0])
 374         # img[0,0] = 0
 375         # tmp.dtype = 'uint32'
 376         # np.savetxt(savepath + '.npy', img)
 377         # test >>
 378
 379         if np.count_nonzero(img) > 0: # matrix not full of zero
 380                 cv2.imwrite(savepath, img, [cv2.IMWRITE_PNG_COMPRESSION, 0]) # write png image
 381
 382         # print("ndimage.imread(savepath)", ndimage.imread(savepath).shape)
 383         # print( ndimage.imread(savepath) )
 384         # print("np.expand_dims(ndimage.imread(savepath), 0)", np.expand_dims(ndimage.imread(savepath), 0).shape)
 385         # print(np.expand_dims(ndimage.imread(savepath), 0))
 386
 387         # test
 388         if verbose:
 389                 return img, dicimg, minmax(epimask)
 390         else:
 391                 return img.shape
 392
 393 def topngs(inputdir, outdir):
 394         """
 395                 inputdir : directory which contains directly dicom files
 396         """
 397         files = [f for f in os.listdir(inputdir)]
 398
 399         for f in files:
 400                 topng( inputdir + f, join(outdir, f) )
 401
 402 if __name__ == '__main__':
 403         # topngs( INPUT_DIR, join(OUT_DIR, INPUT_DIR.split('/')[-2]) )
 404         readpng(OUT_DIR+'aug_circle.png')
 405         # topng(INPUT_DIR+'Image00003', OUT_DIR + INPUT_DIR.split('/')[-2] +'-Image00003', epimask="/Users/user/Desktop/Master/Stage/Data/json/json_GT/0_Case_0002/contours/1.2.3.4.5.6/31/-85.9968/Epicardic.json")