prepare.py

#coding=utf-8
import os
import shutil
import numpy as np
import sys
# sys.path.append('./')
# import config_training
from config_training import config 
# config = config_training.config
from scipy.io import loadmat
import numpy as np
import matplotlib.pyplot as plt
import h5py
import pandas
import scipy
from scipy.ndimage.interpolation import zoom
from skimage import measure
import SimpleITK as sitk
from scipy.ndimage.morphology import binary_dilation,generate_binary_structure
from skimage.morphology import convex_hull_image
import pandas
from multiprocessing import Pool
from functools import partial
import torch
import warnings
#######configuration of GPU 
# import os
# import tensorflow as tf
# from keras.backend.tensorflow_backend import set_session
# os.environ["CUDA_VISIBLE_DEVICES"]="3"
# config_GPU = tf.ConfigProto()
# config_GPU.gpu_options.per_process_gpu_memory_fraction = 0.5
# config_GPU.gpu_options.allow_growth = True
# set_session(tf.Session(config=config_GPU))
##########
def resample(imgs, spacing, new_spacing,order=2):
    # print('imgs.shape',imgs.shape)
    if len(imgs.shape)==3:
        new_shape = np.round(imgs.shape * spacing / new_spacing)
        true_spacing = spacing * imgs.shape / new_shape
        # print('000',true_spacing,spacing * imgs.shape,new_shape)#float
        resize_factor = new_shape / imgs.shape
        imgs = zoom(imgs, resize_factor, mode = 'nearest',order=order)
        return imgs, true_spacing
    elif len(imgs.shape)==4:
        n = imgs.shape[-1]
        newimg = []
        for i in range(n):
            slice = imgs[:,:,:,i]
            newslice,true_spacing = resample(slice,spacing,new_spacing)
            newimg.append(newslice)
        newimg=np.transpose(np.array(newimg),[1,2,3,0])
        return newimg,true_spacing
    else:
        raise ValueError('wrong shape')
def worldToVoxelCoord(worldCoord, origin, spacing):
     
    stretchedVoxelCoord = np.absolute(worldCoord - origin)
    voxelCoord = stretchedVoxelCoord / spacing
    return voxelCoord

def load_itk_image(filename):
    with open(filename) as f:
        contents = f.readlines()
        # print('contents',filename)#/public/share/jiezhao/Minerva/Lung/data/luna16/subset_data/subset0/1.3.6.1.4.1.14519.5.2.1.6279.6001.1
        # print('contents',contents)
        '''contents ['ObjectType = Image\n', 'NDims = 3\n', 'BinaryData = True\n', 'BinaryDataByteOrderMSB = False\n', 'CompressedData = False\n', 'TransformMatrix = 1 0 0 0 1 0 0 0 1\n', 'Offset = -185 -205 -380.80999800000001\n', 'CenterOfRotation = 0 0 0\n', 'AnatomicalOrientation = RAI\n', 'ElementSpacing = 0.7226560115814209 0.7226560115814209 2.5\n', 'DimSize = 512 512 140\n', 'ElementType = MET_SHORT\n', 'ElementDataFile = 1.3.6.1.4.1.14519.5.2.1.6279.6001.139713436241461669335487719526.raw\n']'''

        line = [k for k in contents if k.startswith('TransformMatrix')][0]#TransformMatrix = 1 0 0 0 1 0 0 0 1
        
        transformM = np.array(line.split(' = ')[1].split(' ')).astype('float')
        '''transformM [1. 0. 0. 0. 1. 0. 0. 0. 1.],图像矩阵是否翻转的标志.(在现实中CT扫描的时候有的人是正卧,有的人是仰卧的,所以会导致图像会出现翻转的情况.)'''
        
        transformM = np.round(transformM)
        if np.any( transformM!=np.array([1,0,0, 0, 1, 0, 0, 0, 1])):
            isflip = True
        else:
            isflip = False

    itkimage = sitk.ReadImage(filename)
    numpyImage = sitk.GetArrayFromImage(itkimage)
     
    numpyOrigin = np.array(list(reversed(itkimage.GetOrigin())))
    numpySpacing = np.array(list(reversed(itkimage.GetSpacing())))#numpyOrigin [-298.774994 -143.100006 -180.600006] [0.625      0.74218798 0.74218798]Z-X-Y
    
    return numpyImage, numpyOrigin, numpySpacing,isflip

def process_mask(mask):
    convex_mask = np.copy(mask)
    for i_layer in range(convex_mask.shape[0]):
        mask1  = np.ascontiguousarray(mask[i_layer])
        if np.sum(mask1)>0:
            mask2 = convex_hull_image(mask1)#将图片中的所有目标看作一个整体，因此计算出来只有一个最小凸多边形
            if np.sum(mask2)>1.5*np.sum(mask1):
                mask2 = mask1
        else:
            mask2 = mask1
        convex_mask[i_layer] = mask2
    struct = generate_binary_structure(3,1)  
    dilatedMask = binary_dilation(convex_mask,structure=struct,iterations=10) 
    return dilatedMask


def lumTrans(img):
    lungwin = np.array([-1200.,600.])
    # print('np.max(img), np.min(img)',np.max(img), np.min(img))
    newimg = (img-lungwin[0])/(lungwin[1]-lungwin[0])
    # print('np.max(newimg), np.min(newimg)',np.max(newimg), np.min(newimg))
    newimg[newimg<0]=0
    newimg[newimg>1]=1
    newimg = (newimg*255).astype('uint8')
    return newimg

def binarize_per_slice(image, spacing, intensity_th=-600, sigma=1, area_th=30, eccen_th=0.99, bg_patch_size=10):
    bw = np.zeros(image.shape, dtype=bool)
    
    # prepare a mask, with all corner values set to nan
    image_size = image.shape[1]
    print('image_size',image_size)
    grid_axis = np.linspace(-image_size/2+0.5, image_size/2-0.5, image_size)
    x, y = np.meshgrid(grid_axis, grid_axis)
    d = (x**2+y**2)**0.5
    nan_mask = (d<image_size/2).astype(float)
    nan_mask[nan_mask == 0] = np.nan
    for i in range(image.shape[0]):
        # Check if corner pixels are identical, if so the slice  before Gaussian filtering
        if len(np.unique(image[i, 0:bg_patch_size, 0:bg_patch_size])) == 1:
            current_bw = scipy.ndimage.filters.gaussian_filter(np.multiply(image[i].astype('float32'), nan_mask), sigma, truncate=2.0) < intensity_th
        else:
            current_bw = scipy.ndimage.filters.gaussian_filter(image[i].astype('float32'), sigma, truncate=2.0) < intensity_th
        
        # select proper components
        label = measure.label(current_bw)
        properties = measure.regionprops(label)
        valid_label = set()
        for prop in properties:
            if prop.area * spacing[1] * spacing[2] > area_th and prop.eccentricity < eccen_th:
                valid_label.add(prop.label)
        current_bw = np.in1d(label, list(valid_label)).reshape(label.shape)
        bw[i] = current_bw
        
    return bw

def all_slice_analysis(bw, spacing, cut_num=0, vol_limit=[0.68, 8.2], area_th=6e3, dist_th=62):
    # in some cases, several top layers need to be removed first
    if cut_num > 0:
        bw0 = np.copy(bw)
        bw[-cut_num:] = False
    label = measure.label(bw, connectivity=1)
    # remove components access to corners
    mid = int(label.shape[2] / 2)
    bg_label = set([label[0, 0, 0], label[0, 0, -1], label[0, -1, 0], label[0, -1, -1], \
                    label[-1-cut_num, 0, 0], label[-1-cut_num, 0, -1], label[-1-cut_num, -1, 0], label[-1-cut_num, -1, -1], \
                    label[0, 0, mid], label[0, -1, mid], label[-1-cut_num, 0, mid], label[-1-cut_num, -1, mid]])
    for l in bg_label:
        label[label == l] = 0
        
    # select components based on volume
    properties = measure.regionprops(label)
    for prop in properties:
        if prop.area * spacing.prod() < vol_limit[0] * 1e6 or prop.area * spacing.prod() > vol_limit[1] * 1e6:
            label[label == prop.label] = 0
            
    # prepare a distance map for further analysis
    x_axis = np.linspace(-label.shape[1]/2+0.5, label.shape[1]/2-0.5, label.shape[1]) * spacing[1]
    y_axis = np.linspace(-label.shape[2]/2+0.5, label.shape[2]/2-0.5, label.shape[2]) * spacing[2]
    x, y = np.meshgrid(x_axis, y_axis)
    d = (x**2+y**2)**0.5
    vols = measure.regionprops(label)
    valid_label = set()
    # select components based on their area and distance to center axis on all slices
    for vol in vols:
        single_vol = label == vol.label
        slice_area = np.zeros(label.shape[0])
        min_distance = np.zeros(label.shape[0])
        for i in range(label.shape[0]):
            slice_area[i] = np.sum(single_vol[i]) * np.prod(spacing[1:3])
            min_distance[i] = np.min(single_vol[i] * d + (1 - single_vol[i]) * np.max(d))
        
        if np.average([min_distance[i] for i in range(label.shape[0]) if slice_area[i] > area_th]) < dist_th:
            valid_label.add(vol.label)
            
    bw = np.in1d(label, list(valid_label)).reshape(label.shape)
    
    # fill back the parts removed earlier
    if cut_num > 0:
        # bw1 is bw with removed slices, bw2 is a dilated version of bw, part of their intersection is returned as final mask
        bw1 = np.copy(bw)
        bw1[-cut_num:] = bw0[-cut_num:]
        bw2 = np.copy(bw)
        bw2 = scipy.ndimage.binary_dilation(bw2, iterations=cut_num)
        bw3 = bw1 & bw2
        label = measure.label(bw, connectivity=1)
        label3 = measure.label(bw3, connectivity=1)
        l_list = list(set(np.unique(label)) - {0})
        valid_l3 = set()
        for l in l_list:
            indices = np.nonzero(label==l)
            l3 = label3[indices[0][0], indices[1][0], indices[2][0]]
            if l3 > 0:
                valid_l3.add(l3)
        bw = np.in1d(label3, list(valid_l3)).reshape(label3.shape)
    
    return bw, len(valid_label)

def fill_hole(bw):
    # fill 3d holes
    label = measure.label(~bw)
    # idendify corner components
    bg_label = set([label[0, 0, 0], label[0, 0, -1], label[0, -1, 0], label[0, -1, -1], \
                    label[-1, 0, 0], label[-1, 0, -1], label[-1, -1, 0], label[-1, -1, -1]])
    bw = ~np.in1d(label, list(bg_label)).reshape(label.shape)
    
    return bw

def two_lung_only(bw, spacing, max_iter=22, max_ratio=4.8):    
    def extract_main(bw, cover=0.95):
        for i in range(bw.shape[0]):
            current_slice = bw[i]
            label = measure.label(current_slice)
            properties = measure.regionprops(label)
            properties.sort(key=lambda x: x.area, reverse=True)
            area = [prop.area for prop in properties]
            count = 0
            sum = 0
            while sum < np.sum(area)*cover:
                sum = sum+area[count]
                count = count+1
            filter = np.zeros(current_slice.shape, dtype=bool)
            for j in range(count):
                bb = properties[j].bbox
                filter[bb[0]:bb[2], bb[1]:bb[3]] = filter[bb[0]:bb[2], bb[1]:bb[3]] | properties[j].convex_image
            bw[i] = bw[i] & filter
           
        label = measure.label(bw)
        properties = measure.regionprops(label)
        properties.sort(key=lambda x: x.area, reverse=True)
        bw = label==properties[0].label

        return bw
    
    def fill_2d_hole(bw):
        for i in range(bw.shape[0]):
            current_slice = bw[i]
            label = measure.label(current_slice)
            properties = measure.regionprops(label)
            for prop in properties:
                bb = prop.bbox
                current_slice[bb[0]:bb[2], bb[1]:bb[3]] = current_slice[bb[0]:bb[2], bb[1]:bb[3]] | prop.filled_image
            bw[i] = current_slice

        return bw
    
    found_flag = False
    iter_count = 0
    bw0 = np.copy(bw)
    while not found_flag and iter_count < max_iter:
        label = measure.label(bw, connectivity=2)
        properties = measure.regionprops(label)
        properties.sort(key=lambda x: x.area, reverse=True)
        if len(properties) > 1 and properties[0].area/properties[1].area < max_ratio:
            found_flag = True
            bw1 = label == properties[0].label
            bw2 = label == properties[1].label
        else:
            bw = scipy.ndimage.binary_erosion(bw)
            iter_count = iter_count + 1
    
    if found_flag:
        d1 = scipy.ndimage.morphology.distance_transform_edt(bw1 == False, sampling=spacing)
        d2 = scipy.ndimage.morphology.distance_transform_edt(bw2 == False, sampling=spacing)
        bw1 = bw0 & (d1 < d2)
        bw2 = bw0 & (d1 > d2)
                
        bw1 = extract_main(bw1)
        bw2 = extract_main(bw2)
        
    else:
        bw1 = bw0
        bw2 = np.zeros(bw.shape).astype('bool')
        
    bw1 = fill_2d_hole(bw1)
    bw2 = fill_2d_hole(bw2)
    bw = bw1 | bw2

    return bw1, bw2, bw

def step1_python_tianchi(case_path):
    # case = load_scan(case_path)
    # case_pixels, spacing = get_pixels_hu(case)
    ''' For the mhd file reader '''
    resolution = np.array([1,1,1])
    sliceim,origin,spacing,isflip = load_itk_image(case_path+'.mhd')
    if isflip:
        sliceim = sliceim[:,::-1,::-1]
        print('flip!')
    # sliceim = lumTrans(sliceim)
    # sliceim1,_ = resample(sliceim,spacing,resolution,order=1)
    case_pixels = np.array(sliceim)

    bw = binarize_per_slice(case_pixels, spacing)
    flag = 0
    cut_num = 0
    cut_step = 2
    bw0 = np.copy(bw)
    while flag == 0 and cut_num < bw.shape[0]:
        bw = np.copy(bw0)
        bw, flag = all_slice_analysis(bw, spacing, cut_num=cut_num, vol_limit=[0.68,7.5])
        cut_num = cut_num + cut_step

    bw = fill_hole(bw)
    bw1, bw2, bw = two_lung_only(bw, spacing)
    return case_pixels, bw1, bw2, spacing, origin, isflip

def savenpy(id, annos, filelist, data_path, prep_folder):    
    resolution = np.array([1, 1, 1])
    name = filelist[id]
    im, m1, m2, spacing, origin, isflip = step1_python_tianchi(os.path.join(data_path, name))
    missingmask = False
    if os.path.exists(os.path.join(prep_folder,name+'_clean.npy')) and \
        os.path.exists(os.path.join(prep_folder,name+'_originbox.npy')) and \
        os.path.exists(os.path.join(prep_folder,name+'_spacing.npy')) and \
        os.path.exists(os.path.join(prep_folder,name+'_origin.npy')) and \
        os.path.exists(os.path.join(prep_folder,name+'_label.npy')):
        if not isflip:
            print('skip', name)
            return
        else:
            missingmask = True

    print('process', name)
    label = annos[annos[:,0]==name]
    # label = label.astype('float')
    label = label[:, [3,1,2,4]].astype('float') # z, y, x, d
    
    
    Mask = m1 + m2
    
    newshape = np.round(np.array(Mask.shape) * spacing / resolution)
    xx,yy,zz = np.where(Mask)
    if xx.size == 0 or yy.size == 0 or zz.size == 0:
        print(name)
        assert 1 == 0

    box = np.array([[np.min(xx), np.max(xx)], [np.min(yy), np.max(yy)], [np.min(zz), np.max(zz)]])
    box = box * np.expand_dims(spacing, 1) / np.expand_dims(resolution, 1)
    box = np.floor(box).astype('int')
    margin = 5
    extendbox = np.vstack([np.max([[0,0,0],box[:,0]-margin],0),np.min([newshape,box[:,1]+2*margin],axis=0).T]).T
    extendbox = extendbox.astype('int')
    if extendbox[0,0] == extendbox[0,1] or extendbox[1,0] == extendbox[1,1] or extendbox[2,0] == extendbox[2,1]:
        print(name)
        assert 1==0

    convex_mask = m1
    dm1 = process_mask(m1)
    dm2 = process_mask(m2)
    dilatedMask = dm1+dm2
    Mask = m1+m2
    if missingmask:
        np.save(os.path.join(prep_folder,name+'_mask.npy'), Mask)
        print('skip', name)
        return
    extramask = dilatedMask - Mask
    bone_thresh = 210
    pad_value = 170
    im[np.isnan(im)]=-2000
    sliceim = lumTrans(im)
    sliceim = sliceim*dilatedMask+pad_value*(1-dilatedMask).astype('uint8')
    bones = sliceim*extramask>bone_thresh
    sliceim[bones] = pad_value
    sliceim1,_ = resample(sliceim,spacing,resolution,order=1)
    sliceim2 = sliceim1[extendbox[0,0]:extendbox[0,1],
                extendbox[1,0]:extendbox[1,1],
                extendbox[2,0]:extendbox[2,1]]
    sliceim = sliceim2[np.newaxis,...]
    print('-------------',sliceim.shape,extendbox.shape,spacing.shape,origin.shape)
    # np.save(os.path.join(prep_folder,name+'_clean.npy'), sliceim)#经过
    # np.save(os.path.join(prep_folder,name+'_originbox.npy'), extendbox)
    # np.save(os.path.join(prep_folder,name+'_spacing.npy'), spacing)
    # np.save(os.path.join(prep_folder,name+'_origin.npy'), origin)
    print(im.shape, '_clean', sliceim.shape, '_originbox', extendbox.shape, '_space', spacing, '_origin', origin)

    this_annos = np.copy(annos[annos[:,0]==name])
    label = []
    print('label', this_annos.shape, name)
    if len(this_annos)>0:
        
        for c in this_annos:
            pos = worldToVoxelCoord(c[1:4][::-1],origin=origin,spacing=spacing)
            if isflip:
                pos[1:] = Mask.shape[1:3]-pos[1:]
            label.append(np.concatenate([pos,[c[4]/spacing[1]]]))
        
    label = np.array(label)
    if len(label)==0:
        label2 = np.array([[0,0,0,0]])
    else:
        label2 = np.copy(label).T
        label2[:3] = label2[:3]*np.expand_dims(spacing,1)/np.expand_dims(resolution,1)
        label2[3] = label2[3]*spacing[1]/resolution[1]
        label2[:3] = label2[:3]-np.expand_dims(extendbox[:,0],1)
        label2 = label2[:4].T
    print('-------------',label2.shape)
    # np.save(os.path.join(prep_folder,name+'_label.npy'),label2)
    print(name)

def full_prep(train=True, val=True, test=True):
    warnings.filterwarnings("ignore")
    #preprocess_result_path = './prep_result'
    train_prep_folder = config['train_preprocess_result_path']
    val_prep_folder = config['val_preprocess_result_path']
    test_prep_folder = config['test_preprocess_result_path']
    
    train_data_path = config['train_data_path']
    val_data_path = config['val_data_path']
    test_data_path = config['test_data_path']
    
    finished_flag = '.flag_preptianchi'
    
    if not os.path.exists(finished_flag):
        trainlabelfiles = config['train_annos_path']
        vallabelfiles = config['val_annos_path']
        testlabelfiles = config['test_annos_path']

        traincontent = np.array(pandas.read_csv(trainlabelfiles))
        traincontent = traincontent[traincontent[:, 0] != np.nan]
        trainalllabel = traincontent[1:, :] # filename, x, y, z, d
        trainfilelist = []
        for f in os.listdir(config['train_data_path']):
            if f.endswith('.mhd'):
                if f[:-4] in config['black_list']:
                    continue
                trainfilelist.append(f[:-4])

        valcontent = np.array(pandas.read_csv(vallabelfiles))
        valcontent = valcontent[valcontent[:, 0] != np.nan]
        valalllabel = valcontent[1:, :] # filename, x, y, z, d
        valfilelist = []
        for f in os.listdir(config['val_data_path']):
            if f.endswith('.mhd'):
                if f[:-4] in config['black_list']:
                    continue
                valfilelist.append(f[:-4])

        testcontent = np.array(pandas.read_csv(testlabelfiles))
        testcontent = testcontent[testcontent[:, 0] != np.nan]
        testalllabel = testcontent[1:, :] # filename, x, y, z, d
        testfilelist = []
        for f in os.listdir(config['test_data_path']):
            if f.endswith('.mhd'):
                if f[:-4] in config['black_list']:
                    continue
                testfilelist.append(f[:-4])

        if not os.path.exists(train_prep_folder):
            os.mkdir(train_prep_folder)
        if not os.path.exists(val_prep_folder):
            os.mkdir(val_prep_folder)
        if not os.path.exists(test_prep_folder):
            os.mkdir(test_prep_folder)
        #eng.addpath('preprocessing/',nargout=0)
        if train:
            print('starting train preprocessing')
            pool = Pool(10)
            partial_savenpy = partial(savenpy, annos=trainalllabel, filelist=trainfilelist, data_path=train_data_path, prep_folder=train_prep_folder)
            N = len(trainfilelist)
            savenpy(1)
            _=pool.map(partial_savenpy, range(N))
            print('end train preprocessing')
        if val:
            print('starting val preprocessing')
            partial_savenpy = partial(savenpy, annos=valalllabel, filelist=valfilelist, data_path=val_data_path, prep_folder=val_prep_folder)
            N = len(valfilelist)
            savenpy(1)
            _=pool.map(partial_savenpy, range(N))
            print('end val preprocessing')
        if test:
            print('starting test preprocessing')
            partial_savenpy = partial(savenpy, annos=testalllabel, filelist=testfilelist, data_path=test_data_path, prep_folder=test_prep_folder)
            N = len(testfilelist)
            savenpy(1)
            _=pool.map(partial_savenpy, range(N))
            pool.close()
            pool.join()
            print('end test preprocessing')
    f= open(finished_flag,"w+")        
    
def splitvaltestcsv():
    testfiles = []
    for f in os.listdir(config['test_data_path']):
        if f.endswith('.mhd'):
            testfiles.append(f[:-4])
    valcsvlines = []
    testcsvlines = []
    import csv 
    valf = open(config['val_annos_path'], 'r')
    valfcsv = csv.reader(valf)
    for line in valfcsv:
        if line[0] in testfiles:
            testcsvlines.append(line)
        else:
            valcsvlines.append(line)
    valf.close()
    testf = open(config['test_annos_path']+'annotations.csv', 'w')
    testfcsv = csv.writer(testf)
    for line in testcsvlines:
        testfcsv.writerow(line)
    testf.close()
    valf = open(config['val_annos_path'], 'w')
    valfcsv = csv.writer(valf)
    for line in valcsvlines:
        valfcsv.writerow(line)
    valf.close()

def savenpy_luna(id, annos, filelist, luna_segment, luna_data,savepath):
    #annos=annos, filelist=filelist, luna_segment=luna_segment, luna_data=luna_data+'subset'+str(setidx)+'/', savepath=savepath+'subset'+str(setidx)+'/'
    islabel = True
    isClean = True
    resolution = np.array([1,1,1])
#     resolution = np.array([2,2,2])
    name = filelist[id]
    # print('name',id, name)#name 0 1.3.6.1.4.1.14519.5.2.1.6279.6001.303421828981831854739626597495
    #name 1 1.3.6.1.4.1.14519.5.2.1.6279.6001.219087313261026510628926082729
    
    sliceim,origin,spacing,isflip = load_itk_image(os.path.join(luna_data,name+'.mhd'))
    # plt.imshow(sliceim[100])
    Mask,origin,spacing,isflip = load_itk_image(os.path.join(luna_segment,name+'.mhd'))
    if isflip:
        Mask = Mask[:,::-1,::-1]
    newshape = np.round(np.array(Mask.shape)*spacing/resolution).astype('int')#Uniform Thickness
    m1 = Mask==3#左肺
    m2 = Mask==4#右肺
    Mask = m1+m2
    
    xx,yy,zz= np.where(Mask)
    box = np.array([[np.min(xx),np.max(xx)],[np.min(yy),np.max(yy)],[np.min(zz),np.max(zz)]])#[[ 38 640][ 88 425][ 16 495]]
    
    box = box*np.expand_dims(spacing,1)/np.expand_dims(resolution,1)
    #[[ 19. 320.][ 53.453125 258.15429688][9.71875 300.67382812]]
    box = np.floor(box).astype('int')#[[ 19 320][53 258][9 300]]
    margin = 5
    extendbox = np.vstack([np.max([[0,0,0],box[:,0]-margin],0),np.min([newshape,box[:,1]+2*margin],axis=0).T]).T#[[ 14 330][ 48 268][4 310]]
    this_annos = np.copy(annos[annos[:,0]==(name)])#[['1.3.6.1.4.1.14519.5.2.1.6279.6001.317087518531899043292346860596' -87.92743528 -49.5343381 -117.00081509999998 9.443896197]]     
    
    if isClean:
        convex_mask = m1
        dm1 = process_mask(m1)#对掩码采取膨胀操作
        dm2 = process_mask(m2)
        dilatedMask = dm1+dm2
        Mask = m1+m2#dilated

        extramask = dilatedMask ^ Mask#异或操作
        bone_thresh = 210
        pad_value = 170

        if isflip:
            sliceim = sliceim[:,::-1,::-1]
            print('flip!')
        #预处理Three automated preprocessing steps are employed for the input CT images. First, we clip the raw data into [−1200, 600]. Second, we transform the range linearlyinto [0, 1]. Finally, we use LUNA16’s given segmentation ground truth and remove the background.
        sliceim = lumTrans(sliceim)
        
        sliceim = sliceim*dilatedMask+pad_value*(1-dilatedMask).astype('uint8')#膨胀肺区，其他地方=170
        bones = (sliceim*extramask)>bone_thresh
        sliceim[bones] = pad_value#肺边缘大于bone_thresh=170
        
        sliceim1,_ = resample(sliceim,spacing,resolution,order=1)
        sliceim2 = sliceim1[extendbox[0,0]:extendbox[0,1],
                    extendbox[1,0]:extendbox[1,1],
                    extendbox[2,0]:extendbox[2,1]]#(290, 254, 334) 
        # plt.imshow(sliceim2[100])
        ###################################
        # new_ct = sitk.GetImageFromArray(sliceim2)
        # new_ct.SetDirection(ct.GetDirection())
        # new_ct.SetOrigin(origin)
        # new_ct.SetSpacing((1, 1, 1))
        # new_ct_name = 'volume-' + str(random.randint) + '.nii'
        # sitk.WriteImage(new_ct, os.path.join('./', new_ct_name))
        ##########################################
        sliceim = sliceim2[np.newaxis,...]#(1, 290, 254, 334)
        # print('-------------1',sliceim)#(1, 262, 243, 346) (133, 512, 512) (3, 2) (3,) (3,)
        # print('-------------1',sliceim.shape,Mask.shape,extendbox.shape,spacing.shape,origin.shape)#(1, 262, 243, 346) (133, 512, 512) (3, 2) (3,) (3,)
        print('clean,spacing,extendbox,origin,mask,label to',savepath)#dtype=uint8)
        np.save(os.path.join(savepath, name+'_clean.npy'), sliceim)#dtype=uint8)
        np.save(os.path.join(savepath, name+'_spacing.npy'), spacing)#float
        np.save(os.path.join(savepath, name+'_extendbox.npy'), extendbox)#int
        np.save(os.path.join(savepath, name+'_origin.npy'), origin)#float
        np.save(os.path.join(savepath, name+'_mask.npy'), Mask)#False True
        # print(stop)
    if islabel:
        this_annos = np.copy(annos[annos[:,0]==(name)])
        label = []
        if len(this_annos)>0:
            
            for c in this_annos:#['1.3.6.1.4.1.14519.5.2.1.6279.6001.313334055029671473836954456733' 68.05574564 43.703103399999996 -266.0215516 12.68580657]
                #[-77.98638617 32.80911691 88.23145221] to [222.786891864 270.481017322138 389.3338931958254]
                pos = worldToVoxelCoord(c[1:4][::-1],origin=origin,spacing=spacing)
                # print('c',c)
                # print('pos',pos)
                if isflip:
                    print('isflip',pos[1:],Mask.shape[1:3],pos[1:])
                    pos[1:] = Mask.shape[1:3]-pos[1:]
                
                # print('this_annos',c[4])
                # print('this_annos',spacing[1])
                
                label.append(np.concatenate([pos,[c[4]/spacing[1]]]))
        
        label = np.array(label)#some of [77.63902124 249.74407667712 208.87253984000003 6.534554632959999]
        if len(label)==0:
            label2 = np.array([[0,0,0,0]])
        else:
            #lable
            '''[[68.33709408000001 285.52261784497364 426.1599867758768 24.773647110830687]
               [56.542471040000024 337.2125031836933  408.310142823606277.913589855427735]]'''
            label2 = np.copy(label).T#变成列向量
            '''[[68.33709408000001 56.542471040000024]
               [285.52261784497364 337.2125031836933]
               [426.1599867758768 408.31014282360627]
               [24.773647110830687 7.913589855427735]]'''

            label2[:3] = label2[:3]*np.expand_dims(spacing,1)/np.expand_dims(resolution,1)
            label2[3] = label2[3]*spacing[1]/resolution[1]
            # print('0',label2[:3]*spacing[1], resolution[1])#float
            label2[:3] = label2[:3]-np.expand_dims(extendbox[:,0],1)
            # print('1',label2[:3])
            # print(stop)
            label2 = label2[:4].T
        np.save(os.path.join(savepath,name+'_label.npy'), label2)
        # print('-------------2',label2)#float
        
    print(name)

def preprocess_luna():
    luna_segment = config['luna_segment']#'/public/share/jiezhao/Minerva/Lung/data/luna16/seg-lungs-LUNA16/'
    savepath = config['preprocess_result_path']#'/public/share/jiezhao/Minerva/Lung/data/luna16/LUNA16PROPOCESSPATH/'
    luna_data = config['luna_data']#'/public/share/jiezhao/Minerva/Lung/data/luna16/subset_data/'
    luna_label = config['luna_label']#'/public/share/jiezhao/Minerva/Lung/data/luna16/CSVFILES/annotations.csv'
    finished_flag = '.1flag_preprocessluna'
    print('starting preprocessing luna')
    if not os.path.exists(finished_flag):
        annos = np.array(pandas.read_csv(luna_label))
        pool = Pool()
        if not os.path.exists(savepath):
            os.mkdir(savepath)
        for setidx in range(10):
            print('process subset', setidx)
            filelist = [f.split('.mhd')[0] for f in os.listdir(luna_data+'subset'+str(setidx)) if f.endswith('.mhd') ]#1.3.6.1.4.1.14519.5.2.1.6279.6001.105756658031515062000744821260
            if not os.path.exists(savepath+'subset'+str(setidx)):
                os.mkdir(savepath+'subset'+str(setidx))
            partial_savenpy_luna = partial(savenpy_luna, annos=annos, filelist=filelist, luna_segment=luna_segment, luna_data=luna_data+'subset'+str(setidx)+'/', savepath=savepath+'subset'+str(setidx)+'/')
            N = len(filelist)
            #savenpy(1)
            _=pool.map(partial_savenpy_luna,range(N))
        pool.close()
        pool.join()
    print('end preprocessing luna')
    # f= open(finished_flag,"w+")


if __name__=='__main__':
    preprocess_luna()