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