Mercurial > repos > perssond > unmicst
changeset 0:6bec4fef6b2e draft
"planemo upload for repository https://github.com/ohsu-comp-bio/unmicst commit 73e4cae15f2d7cdc86719e77470eb00af4b6ebb7-dirty"
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/UnMicst.py Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,674 @@ +import numpy as np +from scipy import misc +import tensorflow.compat.v1 as tf +import shutil +import scipy.io as sio +import os, fnmatch, glob +import skimage.exposure as sk +import skimage.io +import argparse +import czifile +from nd2reader import ND2Reader +import tifffile +import sys +tf.disable_v2_behavior() +#sys.path.insert(0, 'C:\\Users\\Public\\Documents\\ImageScience') + +from toolbox.imtools import * +from toolbox.ftools import * +from toolbox.PartitionOfImage import PI2D +from toolbox import GPUselect + +def concat3(lst): + return tf.concat(lst, 3) + + +class UNet2D: + hp = None # hyper-parameters + nn = None # network + tfTraining = None # if training or not (to handle batch norm) + tfData = None # data placeholder + Session = None + DatasetMean = 0 + DatasetStDev = 0 + + def setupWithHP(hp): + UNet2D.setup(hp['imSize'], + hp['nChannels'], + hp['nClasses'], + hp['nOut0'], + hp['featMapsFact'], + hp['downSampFact'], + hp['ks'], + hp['nExtraConvs'], + hp['stdDev0'], + hp['nLayers'], + hp['batchSize']) + + def setup(imSize, nChannels, nClasses, nOut0, featMapsFact, downSampFact, kernelSize, nExtraConvs, stdDev0, + nDownSampLayers, batchSize): + UNet2D.hp = {'imSize': imSize, + 'nClasses': nClasses, + 'nChannels': nChannels, + 'nExtraConvs': nExtraConvs, + 'nLayers': nDownSampLayers, + 'featMapsFact': featMapsFact, + 'downSampFact': downSampFact, + 'ks': kernelSize, + 'nOut0': nOut0, + 'stdDev0': stdDev0, + 'batchSize': batchSize} + + nOutX = [UNet2D.hp['nChannels'], UNet2D.hp['nOut0']] + dsfX = [] + for i in range(UNet2D.hp['nLayers']): + nOutX.append(nOutX[-1] * UNet2D.hp['featMapsFact']) + dsfX.append(UNet2D.hp['downSampFact']) + + # -------------------------------------------------- + # downsampling layer + # -------------------------------------------------- + + with tf.name_scope('placeholders'): + UNet2D.tfTraining = tf.placeholder(tf.bool, name='training') + UNet2D.tfData = tf.placeholder("float", shape=[None, UNet2D.hp['imSize'], UNet2D.hp['imSize'], + UNet2D.hp['nChannels']], name='data') + + def down_samp_layer(data, index): + with tf.name_scope('ld%d' % index): + ldXWeights1 = tf.Variable( + tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index], nOutX[index + 1]], + stddev=stdDev0), name='kernel1') + ldXWeightsExtra = [] + for i in range(nExtraConvs): + ldXWeightsExtra.append(tf.Variable( + tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index + 1], nOutX[index + 1]], + stddev=stdDev0), name='kernelExtra%d' % i)) + + c00 = tf.nn.conv2d(data, ldXWeights1, strides=[1, 1, 1, 1], padding='SAME') + for i in range(nExtraConvs): + c00 = tf.nn.conv2d(tf.nn.relu(c00), ldXWeightsExtra[i], strides=[1, 1, 1, 1], padding='SAME') + + ldXWeightsShortcut = tf.Variable( + tf.truncated_normal([1, 1, nOutX[index], nOutX[index + 1]], stddev=stdDev0), name='shortcutWeights') + shortcut = tf.nn.conv2d(data, ldXWeightsShortcut, strides=[1, 1, 1, 1], padding='SAME') + + bn = tf.layers.batch_normalization(tf.nn.relu(c00 + shortcut), training=UNet2D.tfTraining) + + return tf.nn.max_pool(bn, ksize=[1, dsfX[index], dsfX[index], 1], + strides=[1, dsfX[index], dsfX[index], 1], padding='SAME', name='maxpool') + + # -------------------------------------------------- + # bottom layer + # -------------------------------------------------- + + with tf.name_scope('lb'): + lbWeights1 = tf.Variable(tf.truncated_normal( + [UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[UNet2D.hp['nLayers']], nOutX[UNet2D.hp['nLayers'] + 1]], + stddev=stdDev0), name='kernel1') + + def lb(hidden): + return tf.nn.relu(tf.nn.conv2d(hidden, lbWeights1, strides=[1, 1, 1, 1], padding='SAME'), name='conv') + + # -------------------------------------------------- + # downsampling + # -------------------------------------------------- + + with tf.name_scope('downsampling'): + dsX = [] + dsX.append(UNet2D.tfData) + + for i in range(UNet2D.hp['nLayers']): + dsX.append(down_samp_layer(dsX[i], i)) + + b = lb(dsX[UNet2D.hp['nLayers']]) + + # -------------------------------------------------- + # upsampling layer + # -------------------------------------------------- + + def up_samp_layer(data, index): + with tf.name_scope('lu%d' % index): + luXWeights1 = tf.Variable( + tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index + 1], nOutX[index + 2]], + stddev=stdDev0), name='kernel1') + luXWeights2 = tf.Variable(tf.truncated_normal( + [UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index] + nOutX[index + 1], nOutX[index + 1]], + stddev=stdDev0), name='kernel2') + luXWeightsExtra = [] + for i in range(nExtraConvs): + luXWeightsExtra.append(tf.Variable( + tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index + 1], nOutX[index + 1]], + stddev=stdDev0), name='kernel2Extra%d' % i)) + + outSize = UNet2D.hp['imSize'] + for i in range(index): + outSize /= dsfX[i] + outSize = int(outSize) + + outputShape = [UNet2D.hp['batchSize'], outSize, outSize, nOutX[index + 1]] + us = tf.nn.relu( + tf.nn.conv2d_transpose(data, luXWeights1, outputShape, strides=[1, dsfX[index], dsfX[index], 1], + padding='SAME'), name='conv1') + cc = concat3([dsX[index], us]) + cv = tf.nn.relu(tf.nn.conv2d(cc, luXWeights2, strides=[1, 1, 1, 1], padding='SAME'), name='conv2') + for i in range(nExtraConvs): + cv = tf.nn.relu(tf.nn.conv2d(cv, luXWeightsExtra[i], strides=[1, 1, 1, 1], padding='SAME'), + name='conv2Extra%d' % i) + return cv + + # -------------------------------------------------- + # final (top) layer + # -------------------------------------------------- + + with tf.name_scope('lt'): + ltWeights1 = tf.Variable(tf.truncated_normal([1, 1, nOutX[1], nClasses], stddev=stdDev0), name='kernel') + + def lt(hidden): + return tf.nn.conv2d(hidden, ltWeights1, strides=[1, 1, 1, 1], padding='SAME', name='conv') + + # -------------------------------------------------- + # upsampling + # -------------------------------------------------- + + with tf.name_scope('upsampling'): + usX = [] + usX.append(b) + + for i in range(UNet2D.hp['nLayers']): + usX.append(up_samp_layer(usX[i], UNet2D.hp['nLayers'] - 1 - i)) + + t = lt(usX[UNet2D.hp['nLayers']]) + + sm = tf.nn.softmax(t, -1) + UNet2D.nn = sm + + def train(imPath, logPath, modelPath, pmPath, nTrain, nValid, nTest, restoreVariables, nSteps, gpuIndex, + testPMIndex): + os.environ['CUDA_VISIBLE_DEVICES'] = '%d' % gpuIndex + + outLogPath = logPath + trainWriterPath = pathjoin(logPath, 'Train') + validWriterPath = pathjoin(logPath, 'Valid') + outModelPath = pathjoin(modelPath, 'model.ckpt') + outPMPath = pmPath + + batchSize = UNet2D.hp['batchSize'] + imSize = UNet2D.hp['imSize'] + nChannels = UNet2D.hp['nChannels'] + nClasses = UNet2D.hp['nClasses'] + + # -------------------------------------------------- + # data + # -------------------------------------------------- + + Train = np.zeros((nTrain, imSize, imSize, nChannels)) + Valid = np.zeros((nValid, imSize, imSize, nChannels)) + Test = np.zeros((nTest, imSize, imSize, nChannels)) + LTrain = np.zeros((nTrain, imSize, imSize, nClasses)) + LValid = np.zeros((nValid, imSize, imSize, nClasses)) + LTest = np.zeros((nTest, imSize, imSize, nClasses)) + + print('loading data, computing mean / st dev') + if not os.path.exists(modelPath): + os.makedirs(modelPath) + if restoreVariables: + datasetMean = loadData(pathjoin(modelPath, 'datasetMean.data')) + datasetStDev = loadData(pathjoin(modelPath, 'datasetStDev.data')) + else: + datasetMean = 0 + datasetStDev = 0 + for iSample in range(nTrain + nValid + nTest): + I = im2double(tifread('%s/I%05d_Img.tif' % (imPath, iSample))) + datasetMean += np.mean(I) + datasetStDev += np.std(I) + datasetMean /= (nTrain + nValid + nTest) + datasetStDev /= (nTrain + nValid + nTest) + saveData(datasetMean, pathjoin(modelPath, 'datasetMean.data')) + saveData(datasetStDev, pathjoin(modelPath, 'datasetStDev.data')) + + perm = np.arange(nTrain + nValid + nTest) + np.random.shuffle(perm) + + for iSample in range(0, nTrain): + path = '%s/I%05d_Img.tif' % (imPath, perm[iSample]) + im = im2double(tifread(path)) + Train[iSample, :, :, 0] = (im - datasetMean) / datasetStDev + path = '%s/I%05d_Ant.tif' % (imPath, perm[iSample]) + im = tifread(path) + for i in range(nClasses): + LTrain[iSample, :, :, i] = (im == i + 1) + + for iSample in range(0, nValid): + path = '%s/I%05d_Img.tif' % (imPath, perm[nTrain + iSample]) + im = im2double(tifread(path)) + Valid[iSample, :, :, 0] = (im - datasetMean) / datasetStDev + path = '%s/I%05d_Ant.tif' % (imPath, perm[nTrain + iSample]) + im = tifread(path) + for i in range(nClasses): + LValid[iSample, :, :, i] = (im == i + 1) + + for iSample in range(0, nTest): + path = '%s/I%05d_Img.tif' % (imPath, perm[nTrain + nValid + iSample]) + im = im2double(tifread(path)) + Test[iSample, :, :, 0] = (im - datasetMean) / datasetStDev + path = '%s/I%05d_Ant.tif' % (imPath, perm[nTrain + nValid + iSample]) + im = tifread(path) + for i in range(nClasses): + LTest[iSample, :, :, i] = (im == i + 1) + + # -------------------------------------------------- + # optimization + # -------------------------------------------------- + + tfLabels = tf.placeholder("float", shape=[None, imSize, imSize, nClasses], name='labels') + + globalStep = tf.Variable(0, trainable=False) + learningRate0 = 0.01 + decaySteps = 1000 + decayRate = 0.95 + learningRate = tf.train.exponential_decay(learningRate0, globalStep, decaySteps, decayRate, staircase=True) + + with tf.name_scope('optim'): + loss = tf.reduce_mean(-tf.reduce_sum(tf.multiply(tfLabels, tf.log(UNet2D.nn)), 3)) + updateOps = tf.get_collection(tf.GraphKeys.UPDATE_OPS) + # optimizer = tf.train.MomentumOptimizer(1e-3,0.9) + optimizer = tf.train.MomentumOptimizer(learningRate, 0.9) + # optimizer = tf.train.GradientDescentOptimizer(learningRate) + with tf.control_dependencies(updateOps): + optOp = optimizer.minimize(loss, global_step=globalStep) + + with tf.name_scope('eval'): + error = [] + for iClass in range(nClasses): + labels0 = tf.reshape(tf.to_int32(tf.slice(tfLabels, [0, 0, 0, iClass], [-1, -1, -1, 1])), + [batchSize, imSize, imSize]) + predict0 = tf.reshape(tf.to_int32(tf.equal(tf.argmax(UNet2D.nn, 3), iClass)), + [batchSize, imSize, imSize]) + correct = tf.multiply(labels0, predict0) + nCorrect0 = tf.reduce_sum(correct) + nLabels0 = tf.reduce_sum(labels0) + error.append(1 - tf.to_float(nCorrect0) / tf.to_float(nLabels0)) + errors = tf.tuple(error) + + # -------------------------------------------------- + # inspection + # -------------------------------------------------- + + with tf.name_scope('scalars'): + tf.summary.scalar('avg_cross_entropy', loss) + for iClass in range(nClasses): + tf.summary.scalar('avg_pixel_error_%d' % iClass, error[iClass]) + tf.summary.scalar('learning_rate', learningRate) + with tf.name_scope('images'): + split0 = tf.slice(UNet2D.nn, [0, 0, 0, 0], [-1, -1, -1, 1]) + split1 = tf.slice(UNet2D.nn, [0, 0, 0, 1], [-1, -1, -1, 1]) + if nClasses > 2: + split2 = tf.slice(UNet2D.nn, [0, 0, 0, 2], [-1, -1, -1, 1]) + tf.summary.image('pm0', split0) + tf.summary.image('pm1', split1) + if nClasses > 2: + tf.summary.image('pm2', split2) + merged = tf.summary.merge_all() + + # -------------------------------------------------- + # session + # -------------------------------------------------- + + saver = tf.train.Saver() + sess = tf.Session(config=tf.ConfigProto( + allow_soft_placement=True)) # config parameter needed to save variables when using GPU + + if os.path.exists(outLogPath): + shutil.rmtree(outLogPath) + trainWriter = tf.summary.FileWriter(trainWriterPath, sess.graph) + validWriter = tf.summary.FileWriter(validWriterPath, sess.graph) + + if restoreVariables: + saver.restore(sess, outModelPath) + print("Model restored.") + else: + sess.run(tf.global_variables_initializer()) + + # -------------------------------------------------- + # train + # -------------------------------------------------- + + batchData = np.zeros((batchSize, imSize, imSize, nChannels)) + batchLabels = np.zeros((batchSize, imSize, imSize, nClasses)) + for i in range(nSteps): + # train + + perm = np.arange(nTrain) + np.random.shuffle(perm) + + for j in range(batchSize): + batchData[j, :, :, :] = Train[perm[j], :, :, :] + batchLabels[j, :, :, :] = LTrain[perm[j], :, :, :] + + summary, _ = sess.run([merged, optOp], + feed_dict={UNet2D.tfData: batchData, tfLabels: batchLabels, UNet2D.tfTraining: 1}) + trainWriter.add_summary(summary, i) + + # validation + + perm = np.arange(nValid) + np.random.shuffle(perm) + + for j in range(batchSize): + batchData[j, :, :, :] = Valid[perm[j], :, :, :] + batchLabels[j, :, :, :] = LValid[perm[j], :, :, :] + + summary, es = sess.run([merged, errors], + feed_dict={UNet2D.tfData: batchData, tfLabels: batchLabels, UNet2D.tfTraining: 0}) + validWriter.add_summary(summary, i) + + e = np.mean(es) + print('step %05d, e: %f' % (i, e)) + + if i == 0: + if restoreVariables: + lowestError = e + else: + lowestError = np.inf + + if np.mod(i, 100) == 0 and e < lowestError: + lowestError = e + print("Model saved in file: %s" % saver.save(sess, outModelPath)) + + # -------------------------------------------------- + # test + # -------------------------------------------------- + + if not os.path.exists(outPMPath): + os.makedirs(outPMPath) + + for i in range(nTest): + j = np.mod(i, batchSize) + + batchData[j, :, :, :] = Test[i, :, :, :] + batchLabels[j, :, :, :] = LTest[i, :, :, :] + + if j == batchSize - 1 or i == nTest - 1: + + output = sess.run(UNet2D.nn, + feed_dict={UNet2D.tfData: batchData, tfLabels: batchLabels, UNet2D.tfTraining: 0}) + + for k in range(j + 1): + pm = output[k, :, :, testPMIndex] + gt = batchLabels[k, :, :, testPMIndex] + im = np.sqrt(normalize(batchData[k, :, :, 0])) + imwrite(np.uint8(255 * np.concatenate((im, np.concatenate((pm, gt), axis=1)), axis=1)), + '%s/I%05d.png' % (outPMPath, i - j + k + 1)) + + # -------------------------------------------------- + # save hyper-parameters, clean-up + # -------------------------------------------------- + + saveData(UNet2D.hp, pathjoin(modelPath, 'hp.data')) + + trainWriter.close() + validWriter.close() + sess.close() + + def deploy(imPath, nImages, modelPath, pmPath, gpuIndex, pmIndex): + os.environ['CUDA_VISIBLE_DEVICES'] = '%d' % gpuIndex + + variablesPath = pathjoin(modelPath, 'model.ckpt') + outPMPath = pmPath + + hp = loadData(pathjoin(modelPath, 'hp.data')) + UNet2D.setupWithHP(hp) + + batchSize = UNet2D.hp['batchSize'] + imSize = UNet2D.hp['imSize'] + nChannels = UNet2D.hp['nChannels'] + nClasses = UNet2D.hp['nClasses'] + + # -------------------------------------------------- + # data + # -------------------------------------------------- + + Data = np.zeros((nImages, imSize, imSize, nChannels)) + + datasetMean = loadData(pathjoin(modelPath, 'datasetMean.data')) + datasetStDev = loadData(pathjoin(modelPath, 'datasetStDev.data')) + + for iSample in range(0, nImages): + path = '%s/I%05d_Img.tif' % (imPath, iSample) + im = im2double(tifread(path)) + Data[iSample, :, :, 0] = (im - datasetMean) / datasetStDev + + # -------------------------------------------------- + # session + # -------------------------------------------------- + + saver = tf.train.Saver() + sess = tf.Session(config=tf.ConfigProto( + allow_soft_placement=True)) # config parameter needed to save variables when using GPU + + saver.restore(sess, variablesPath) + print("Model restored.") + + # -------------------------------------------------- + # deploy + # -------------------------------------------------- + + batchData = np.zeros((batchSize, imSize, imSize, nChannels)) + + if not os.path.exists(outPMPath): + os.makedirs(outPMPath) + + for i in range(nImages): + print(i, nImages) + + j = np.mod(i, batchSize) + + batchData[j, :, :, :] = Data[i, :, :, :] + + if j == batchSize - 1 or i == nImages - 1: + + output = sess.run(UNet2D.nn, feed_dict={UNet2D.tfData: batchData, UNet2D.tfTraining: 0}) + + for k in range(j + 1): + pm = output[k, :, :, pmIndex] + im = np.sqrt(normalize(batchData[k, :, :, 0])) + # imwrite(np.uint8(255*np.concatenate((im,pm),axis=1)),'%s/I%05d.png' % (outPMPath,i-j+k+1)) + imwrite(np.uint8(255 * im), '%s/I%05d_Im.png' % (outPMPath, i - j + k + 1)) + imwrite(np.uint8(255 * pm), '%s/I%05d_PM.png' % (outPMPath, i - j + k + 1)) + + # -------------------------------------------------- + # clean-up + # -------------------------------------------------- + + sess.close() + + def singleImageInferenceSetup(modelPath, gpuIndex,mean,std): + variablesPath = pathjoin(modelPath, 'model.ckpt') + + hp = loadData(pathjoin(modelPath, 'hp.data')) + UNet2D.setupWithHP(hp) + if mean ==-1: + UNet2D.DatasetMean = loadData(pathjoin(modelPath, 'datasetMean.data')) + else: + UNet2D.DatasetMean = mean + + if std == -1: + UNet2D.DatasetStDev = loadData(pathjoin(modelPath, 'datasetStDev.data')) + else: + UNet2D.DatasetStDev = std + print(UNet2D.DatasetMean) + print(UNet2D.DatasetStDev) + + # -------------------------------------------------- + # session + # -------------------------------------------------- + + saver = tf.train.Saver() + UNet2D.Session = tf.Session(config=tf.ConfigProto()) + # allow_soft_placement=True)) # config parameter needed to save variables when using GPU + + saver.restore(UNet2D.Session, variablesPath) + print("Model restored.") + + def singleImageInferenceCleanup(): + UNet2D.Session.close() + + def singleImageInference(image, mode, pmIndex): + print('Inference...') + + batchSize = UNet2D.hp['batchSize'] + imSize = UNet2D.hp['imSize'] + nChannels = UNet2D.hp['nChannels'] + + PI2D.setup(image, imSize, int(imSize / 8), mode) + PI2D.createOutput(nChannels) + + batchData = np.zeros((batchSize, imSize, imSize, nChannels)) + for i in range(PI2D.NumPatches): + j = np.mod(i, batchSize) + batchData[j, :, :, 0] = (PI2D.getPatch(i) - UNet2D.DatasetMean) / UNet2D.DatasetStDev + if j == batchSize - 1 or i == PI2D.NumPatches - 1: + output = UNet2D.Session.run(UNet2D.nn, feed_dict={UNet2D.tfData: batchData, UNet2D.tfTraining: 0}) + for k in range(j + 1): + pm = output[k, :, :, pmIndex] + PI2D.patchOutput(i - j + k, pm) + # PI2D.patchOutput(i-j+k,normalize(imgradmag(PI2D.getPatch(i-j+k),1))) + + return PI2D.getValidOutput() + + +if __name__ == '__main__': + parser = argparse.ArgumentParser() + parser.add_argument("imagePath", help="path to the .tif file") + parser.add_argument("--model", help="type of model. For example, nuclei vs cytoplasm",default = 'nucleiDAPI') + parser.add_argument("--outputPath", help="output path of probability map") + parser.add_argument("--channel", help="channel to perform inference on", type=int, default=0) + parser.add_argument("--classOrder", help="background, contours, foreground", type = int, nargs = '+', default=-1) + parser.add_argument("--mean", help="mean intensity of input image. Use -1 to use model", type=float, default=-1) + parser.add_argument("--std", help="mean standard deviation of input image. Use -1 to use model", type=float, default=-1) + parser.add_argument("--scalingFactor", help="factor by which to increase/decrease image size by", type=float, + default=1) + parser.add_argument("--stackOutput", help="save probability maps as separate files", action='store_true') + parser.add_argument("--GPU", help="explicitly select GPU", type=int, default = -1) + parser.add_argument("--outlier", + help="map percentile intensity to max when rescaling intensity values. Max intensity as default", + type=float, default=-1) + args = parser.parse_args() + + logPath = '' + scriptPath = os.path.dirname(os.path.realpath(__file__)) + modelPath = os.path.join(scriptPath, 'models', args.model) + # modelPath = os.path.join(scriptPath, 'models/cytoplasmINcell') + # modelPath = os.path.join(scriptPath, 'cytoplasmZeissNikon') + pmPath = '' + + if os.system('nvidia-smi') == 0: + if args.GPU == -1: + print("automatically choosing GPU") + GPU = GPUselect.pick_gpu_lowest_memory() + else: + GPU = args.GPU + print('Using GPU ' + str(GPU)) + + else: + if sys.platform == 'win32': # only 1 gpu on windows + if args.GPU==-1: + GPU = 0 + else: + GPU = args.GPU + print('Using GPU ' + str(GPU)) + else: + GPU=0 + print('Using CPU') + os.environ['CUDA_VISIBLE_DEVICES'] = '%d' % GPU + UNet2D.singleImageInferenceSetup(modelPath, GPU,args.mean,args.std) + nClass = UNet2D.hp['nClasses'] + imagePath = args.imagePath + dapiChannel = args.channel + dsFactor = args.scalingFactor + parentFolder = os.path.dirname(os.path.dirname(imagePath)) + fileName = os.path.basename(imagePath) + fileNamePrefix = fileName.split(os.extsep, 1) + print(fileName) + fileType = fileNamePrefix[1] + + if fileType=='ome.tif' or fileType == 'btf' : + I = skio.imread(imagePath, img_num=dapiChannel,plugin='tifffile') + elif fileType == 'tif' : + I = tifffile.imread(imagePath, key=dapiChannel) + elif fileType == 'czi': + with czifile.CziFile(imagePath) as czi: + image = czi.asarray() + I = image[0, 0, dapiChannel, 0, 0, :, :, 0] + elif fileType == 'nd2': + with ND2Reader(imagePath) as fullStack: + I = fullStack[dapiChannel] + + if args.classOrder == -1: + args.classOrder = range(nClass) + + rawI = I + print(type(I)) + hsize = int((float(I.shape[0]) * float(dsFactor))) + vsize = int((float(I.shape[1]) * float(dsFactor))) + I = resize(I, (hsize, vsize)) + if args.outlier == -1: + maxLimit = np.max(I) + else: + maxLimit = np.percentile(I, args.outlier) + I = im2double(sk.rescale_intensity(I, in_range=(np.min(I), maxLimit), out_range=(0, 0.983))) + rawI = im2double(rawI) / np.max(im2double(rawI)) + if not args.outputPath: + args.outputPath = parentFolder + '//probability_maps' + + if not os.path.exists(args.outputPath): + os.makedirs(args.outputPath) + + append_kwargs = { + 'bigtiff': True, + 'metadata': None, + 'append': True, + } + save_kwargs = { + 'bigtiff': True, + 'metadata': None, + 'append': False, + } + if args.stackOutput: + slice=0 + for iClass in args.classOrder[::-1]: + PM = np.uint8(255*UNet2D.singleImageInference(I, 'accumulate', iClass)) # backwards in order to align with ilastik... + PM = resize(PM, (rawI.shape[0], rawI.shape[1])) + if slice==0: + skimage.io.imsave(args.outputPath + '//' + fileNamePrefix[0] + '_Probabilities_' + str(dapiChannel) + '.tif', np.uint8(255 * PM),**save_kwargs) + else: + skimage.io.imsave(args.outputPath + '//' + fileNamePrefix[0] + '_Probabilities_' + str(dapiChannel) + '.tif',np.uint8(255 * PM),**append_kwargs) + if slice==1: + save_kwargs['append'] = False + skimage.io.imsave(args.outputPath + '//' + fileNamePrefix[0] + '_Preview_' + str(dapiChannel) + '.tif', np.uint8(255 * PM), **save_kwargs) + skimage.io.imsave(args.outputPath + '//' + fileNamePrefix[0] + '_Preview_' + str(dapiChannel) + '.tif', np.uint8(255 * rawI), **append_kwargs) + slice = slice + 1 + + else: + contours = np.uint8(255*UNet2D.singleImageInference(I, 'accumulate', args.classOrder[1])) + hsize = int((float(I.shape[0]) * float(1 / dsFactor))) + vsize = int((float(I.shape[1]) * float(1 / dsFactor))) + contours = resize(contours, (rawI.shape[0], rawI.shape[1])) + skimage.io.imsave(args.outputPath + '//' + fileNamePrefix[0] + '_ContoursPM_' + str(dapiChannel) + '.tif',np.uint8(255 * contours),**save_kwargs) + skimage.io.imsave(args.outputPath + '//' + fileNamePrefix[0] + '_ContoursPM_' + str(dapiChannel) + '.tif',np.uint8(255 * rawI), **append_kwargs) + del contours + nuclei = np.uint8(255*UNet2D.singleImageInference(I, 'accumulate', args.classOrder[2])) + nuclei = resize(nuclei, (rawI.shape[0], rawI.shape[1])) + skimage.io.imsave(args.outputPath + '//' + fileNamePrefix[0] + '_NucleiPM_' + str(dapiChannel) + '.tif',np.uint8(255 * nuclei), **save_kwargs) + del nuclei + UNet2D.singleImageInferenceCleanup() + +#aligned output files to reflect ilastik +#outputting all classes as single file +#handles multiple formats including tif, ome.tif, nd2, czi +#selectable models (human nuclei, mouse nuclei, cytoplasm) + +#added legacy function to save output files +#append save function to reduce memory footprint +#added --classOrder parameter to specify which class is background, contours, and nuclei respectively \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/batchUNet2DTMACycif.py Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,594 @@ +import numpy as np +from scipy import misc +import tensorflow as tf +import shutil +import scipy.io as sio +import os,fnmatch,PIL,glob +import skimage.exposure as sk + +import sys +sys.path.insert(0, 'C:\\Users\\Public\\Documents\\ImageScience') +from toolbox.imtools import * +from toolbox.ftools import * +from toolbox.PartitionOfImage import PI2D + + +def concat3(lst): + return tf.concat(lst,3) + +class UNet2D: + hp = None # hyper-parameters + nn = None # network + tfTraining = None # if training or not (to handle batch norm) + tfData = None # data placeholder + Session = None + DatasetMean = 0 + DatasetStDev = 0 + + def setupWithHP(hp): + UNet2D.setup(hp['imSize'], + hp['nChannels'], + hp['nClasses'], + hp['nOut0'], + hp['featMapsFact'], + hp['downSampFact'], + hp['ks'], + hp['nExtraConvs'], + hp['stdDev0'], + hp['nLayers'], + hp['batchSize']) + + def setup(imSize,nChannels,nClasses,nOut0,featMapsFact,downSampFact,kernelSize,nExtraConvs,stdDev0,nDownSampLayers,batchSize): + UNet2D.hp = {'imSize':imSize, + 'nClasses':nClasses, + 'nChannels':nChannels, + 'nExtraConvs':nExtraConvs, + 'nLayers':nDownSampLayers, + 'featMapsFact':featMapsFact, + 'downSampFact':downSampFact, + 'ks':kernelSize, + 'nOut0':nOut0, + 'stdDev0':stdDev0, + 'batchSize':batchSize} + + nOutX = [UNet2D.hp['nChannels'],UNet2D.hp['nOut0']] + dsfX = [] + for i in range(UNet2D.hp['nLayers']): + nOutX.append(nOutX[-1]*UNet2D.hp['featMapsFact']) + dsfX.append(UNet2D.hp['downSampFact']) + + + # -------------------------------------------------- + # downsampling layer + # -------------------------------------------------- + + with tf.name_scope('placeholders'): + UNet2D.tfTraining = tf.placeholder(tf.bool, name='training') + UNet2D.tfData = tf.placeholder("float", shape=[None,UNet2D.hp['imSize'],UNet2D.hp['imSize'],UNet2D.hp['nChannels']],name='data') + + def down_samp_layer(data,index): + with tf.name_scope('ld%d' % index): + ldXWeights1 = tf.Variable(tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index], nOutX[index+1]], stddev=stdDev0),name='kernel1') + ldXWeightsExtra = [] + for i in range(nExtraConvs): + ldXWeightsExtra.append(tf.Variable(tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index+1], nOutX[index+1]], stddev=stdDev0),name='kernelExtra%d' % i)) + + c00 = tf.nn.conv2d(data, ldXWeights1, strides=[1, 1, 1, 1], padding='SAME') + for i in range(nExtraConvs): + c00 = tf.nn.conv2d(tf.nn.relu(c00), ldXWeightsExtra[i], strides=[1, 1, 1, 1], padding='SAME') + + ldXWeightsShortcut = tf.Variable(tf.truncated_normal([1, 1, nOutX[index], nOutX[index+1]], stddev=stdDev0),name='shortcutWeights') + shortcut = tf.nn.conv2d(data, ldXWeightsShortcut, strides=[1, 1, 1, 1], padding='SAME') + + bn = tf.layers.batch_normalization(tf.nn.relu(c00+shortcut), training=UNet2D.tfTraining) + + return tf.nn.max_pool(bn, ksize=[1, dsfX[index], dsfX[index], 1], strides=[1, dsfX[index], dsfX[index], 1], padding='SAME',name='maxpool') + + # -------------------------------------------------- + # bottom layer + # -------------------------------------------------- + + with tf.name_scope('lb'): + lbWeights1 = tf.Variable(tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[UNet2D.hp['nLayers']], nOutX[UNet2D.hp['nLayers']+1]], stddev=stdDev0),name='kernel1') + def lb(hidden): + return tf.nn.relu(tf.nn.conv2d(hidden, lbWeights1, strides=[1, 1, 1, 1], padding='SAME'),name='conv') + + # -------------------------------------------------- + # downsampling + # -------------------------------------------------- + + with tf.name_scope('downsampling'): + dsX = [] + dsX.append(UNet2D.tfData) + + for i in range(UNet2D.hp['nLayers']): + dsX.append(down_samp_layer(dsX[i],i)) + + b = lb(dsX[UNet2D.hp['nLayers']]) + + # -------------------------------------------------- + # upsampling layer + # -------------------------------------------------- + + def up_samp_layer(data,index): + with tf.name_scope('lu%d' % index): + luXWeights1 = tf.Variable(tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index+1], nOutX[index+2]], stddev=stdDev0),name='kernel1') + luXWeights2 = tf.Variable(tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index]+nOutX[index+1], nOutX[index+1]], stddev=stdDev0),name='kernel2') + luXWeightsExtra = [] + for i in range(nExtraConvs): + luXWeightsExtra.append(tf.Variable(tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index+1], nOutX[index+1]], stddev=stdDev0),name='kernel2Extra%d' % i)) + + outSize = UNet2D.hp['imSize'] + for i in range(index): + outSize /= dsfX[i] + outSize = int(outSize) + + outputShape = [UNet2D.hp['batchSize'],outSize,outSize,nOutX[index+1]] + us = tf.nn.relu(tf.nn.conv2d_transpose(data, luXWeights1, outputShape, strides=[1, dsfX[index], dsfX[index], 1], padding='SAME'),name='conv1') + cc = concat3([dsX[index],us]) + cv = tf.nn.relu(tf.nn.conv2d(cc, luXWeights2, strides=[1, 1, 1, 1], padding='SAME'),name='conv2') + for i in range(nExtraConvs): + cv = tf.nn.relu(tf.nn.conv2d(cv, luXWeightsExtra[i], strides=[1, 1, 1, 1], padding='SAME'),name='conv2Extra%d' % i) + return cv + + # -------------------------------------------------- + # final (top) layer + # -------------------------------------------------- + + with tf.name_scope('lt'): + ltWeights1 = tf.Variable(tf.truncated_normal([1, 1, nOutX[1], nClasses], stddev=stdDev0),name='kernel') + def lt(hidden): + return tf.nn.conv2d(hidden, ltWeights1, strides=[1, 1, 1, 1], padding='SAME',name='conv') + + + # -------------------------------------------------- + # upsampling + # -------------------------------------------------- + + with tf.name_scope('upsampling'): + usX = [] + usX.append(b) + + for i in range(UNet2D.hp['nLayers']): + usX.append(up_samp_layer(usX[i],UNet2D.hp['nLayers']-1-i)) + + t = lt(usX[UNet2D.hp['nLayers']]) + + + sm = tf.nn.softmax(t,-1) + UNet2D.nn = sm + + + def train(imPath,logPath,modelPath,pmPath,nTrain,nValid,nTest,restoreVariables,nSteps,gpuIndex,testPMIndex): + os.environ['CUDA_VISIBLE_DEVICES']= '%d' % gpuIndex + + outLogPath = logPath + trainWriterPath = pathjoin(logPath,'Train') + validWriterPath = pathjoin(logPath,'Valid') + outModelPath = pathjoin(modelPath,'model.ckpt') + outPMPath = pmPath + + batchSize = UNet2D.hp['batchSize'] + imSize = UNet2D.hp['imSize'] + nChannels = UNet2D.hp['nChannels'] + nClasses = UNet2D.hp['nClasses'] + + # -------------------------------------------------- + # data + # -------------------------------------------------- + + Train = np.zeros((nTrain,imSize,imSize,nChannels)) + Valid = np.zeros((nValid,imSize,imSize,nChannels)) + Test = np.zeros((nTest,imSize,imSize,nChannels)) + LTrain = np.zeros((nTrain,imSize,imSize,nClasses)) + LValid = np.zeros((nValid,imSize,imSize,nClasses)) + LTest = np.zeros((nTest,imSize,imSize,nClasses)) + + print('loading data, computing mean / st dev') + if not os.path.exists(modelPath): + os.makedirs(modelPath) + if restoreVariables: + datasetMean = loadData(pathjoin(modelPath,'datasetMean.data')) + datasetStDev = loadData(pathjoin(modelPath,'datasetStDev.data')) + else: + datasetMean = 0 + datasetStDev = 0 + for iSample in range(nTrain+nValid+nTest): + I = im2double(tifread('%s/I%05d_Img.tif' % (imPath,iSample))) + datasetMean += np.mean(I) + datasetStDev += np.std(I) + datasetMean /= (nTrain+nValid+nTest) + datasetStDev /= (nTrain+nValid+nTest) + saveData(datasetMean, pathjoin(modelPath,'datasetMean.data')) + saveData(datasetStDev, pathjoin(modelPath,'datasetStDev.data')) + + perm = np.arange(nTrain+nValid+nTest) + np.random.shuffle(perm) + + for iSample in range(0, nTrain): + path = '%s/I%05d_Img.tif' % (imPath,perm[iSample]) + im = im2double(tifread(path)) + Train[iSample,:,:,0] = (im-datasetMean)/datasetStDev + path = '%s/I%05d_Ant.tif' % (imPath,perm[iSample]) + im = tifread(path) + for i in range(nClasses): + LTrain[iSample,:,:,i] = (im == i+1) + + for iSample in range(0, nValid): + path = '%s/I%05d_Img.tif' % (imPath,perm[nTrain+iSample]) + im = im2double(tifread(path)) + Valid[iSample,:,:,0] = (im-datasetMean)/datasetStDev + path = '%s/I%05d_Ant.tif' % (imPath,perm[nTrain+iSample]) + im = tifread(path) + for i in range(nClasses): + LValid[iSample,:,:,i] = (im == i+1) + + for iSample in range(0, nTest): + path = '%s/I%05d_Img.tif' % (imPath,perm[nTrain+nValid+iSample]) + im = im2double(tifread(path)) + Test[iSample,:,:,0] = (im-datasetMean)/datasetStDev + path = '%s/I%05d_Ant.tif' % (imPath,perm[nTrain+nValid+iSample]) + im = tifread(path) + for i in range(nClasses): + LTest[iSample,:,:,i] = (im == i+1) + + # -------------------------------------------------- + # optimization + # -------------------------------------------------- + + tfLabels = tf.placeholder("float", shape=[None,imSize,imSize,nClasses],name='labels') + + globalStep = tf.Variable(0,trainable=False) + learningRate0 = 0.01 + decaySteps = 1000 + decayRate = 0.95 + learningRate = tf.train.exponential_decay(learningRate0,globalStep,decaySteps,decayRate,staircase=True) + + with tf.name_scope('optim'): + loss = tf.reduce_mean(-tf.reduce_sum(tf.multiply(tfLabels,tf.log(UNet2D.nn)),3)) + updateOps = tf.get_collection(tf.GraphKeys.UPDATE_OPS) + # optimizer = tf.train.MomentumOptimizer(1e-3,0.9) + optimizer = tf.train.MomentumOptimizer(learningRate,0.9) + # optimizer = tf.train.GradientDescentOptimizer(learningRate) + with tf.control_dependencies(updateOps): + optOp = optimizer.minimize(loss,global_step=globalStep) + + with tf.name_scope('eval'): + error = [] + for iClass in range(nClasses): + labels0 = tf.reshape(tf.to_int32(tf.slice(tfLabels,[0,0,0,iClass],[-1,-1,-1,1])),[batchSize,imSize,imSize]) + predict0 = tf.reshape(tf.to_int32(tf.equal(tf.argmax(UNet2D.nn,3),iClass)),[batchSize,imSize,imSize]) + correct = tf.multiply(labels0,predict0) + nCorrect0 = tf.reduce_sum(correct) + nLabels0 = tf.reduce_sum(labels0) + error.append(1-tf.to_float(nCorrect0)/tf.to_float(nLabels0)) + errors = tf.tuple(error) + + # -------------------------------------------------- + # inspection + # -------------------------------------------------- + + with tf.name_scope('scalars'): + tf.summary.scalar('avg_cross_entropy', loss) + for iClass in range(nClasses): + tf.summary.scalar('avg_pixel_error_%d' % iClass, error[iClass]) + tf.summary.scalar('learning_rate', learningRate) + with tf.name_scope('images'): + split0 = tf.slice(UNet2D.nn,[0,0,0,0],[-1,-1,-1,1]) + split1 = tf.slice(UNet2D.nn,[0,0,0,1],[-1,-1,-1,1]) + if nClasses > 2: + split2 = tf.slice(UNet2D.nn,[0,0,0,2],[-1,-1,-1,1]) + tf.summary.image('pm0',split0) + tf.summary.image('pm1',split1) + if nClasses > 2: + tf.summary.image('pm2',split2) + merged = tf.summary.merge_all() + + + # -------------------------------------------------- + # session + # -------------------------------------------------- + + saver = tf.train.Saver() + sess = tf.Session(config=tf.ConfigProto(allow_soft_placement=True)) # config parameter needed to save variables when using GPU + + if os.path.exists(outLogPath): + shutil.rmtree(outLogPath) + trainWriter = tf.summary.FileWriter(trainWriterPath, sess.graph) + validWriter = tf.summary.FileWriter(validWriterPath, sess.graph) + + if restoreVariables: + saver.restore(sess, outModelPath) + print("Model restored.") + else: + sess.run(tf.global_variables_initializer()) + + # -------------------------------------------------- + # train + # -------------------------------------------------- + + batchData = np.zeros((batchSize,imSize,imSize,nChannels)) + batchLabels = np.zeros((batchSize,imSize,imSize,nClasses)) + for i in range(nSteps): + # train + + perm = np.arange(nTrain) + np.random.shuffle(perm) + + for j in range(batchSize): + batchData[j,:,:,:] = Train[perm[j],:,:,:] + batchLabels[j,:,:,:] = LTrain[perm[j],:,:,:] + + summary,_ = sess.run([merged,optOp],feed_dict={UNet2D.tfData: batchData, tfLabels: batchLabels, UNet2D.tfTraining: 1}) + trainWriter.add_summary(summary, i) + + # validation + + perm = np.arange(nValid) + np.random.shuffle(perm) + + for j in range(batchSize): + batchData[j,:,:,:] = Valid[perm[j],:,:,:] + batchLabels[j,:,:,:] = LValid[perm[j],:,:,:] + + summary, es = sess.run([merged, errors],feed_dict={UNet2D.tfData: batchData, tfLabels: batchLabels, UNet2D.tfTraining: 0}) + validWriter.add_summary(summary, i) + + e = np.mean(es) + print('step %05d, e: %f' % (i,e)) + + if i == 0: + if restoreVariables: + lowestError = e + else: + lowestError = np.inf + + if np.mod(i,100) == 0 and e < lowestError: + lowestError = e + print("Model saved in file: %s" % saver.save(sess, outModelPath)) + + + # -------------------------------------------------- + # test + # -------------------------------------------------- + + if not os.path.exists(outPMPath): + os.makedirs(outPMPath) + + for i in range(nTest): + j = np.mod(i,batchSize) + + batchData[j,:,:,:] = Test[i,:,:,:] + batchLabels[j,:,:,:] = LTest[i,:,:,:] + + if j == batchSize-1 or i == nTest-1: + + output = sess.run(UNet2D.nn,feed_dict={UNet2D.tfData: batchData, tfLabels: batchLabels, UNet2D.tfTraining: 0}) + + for k in range(j+1): + pm = output[k,:,:,testPMIndex] + gt = batchLabels[k,:,:,testPMIndex] + im = np.sqrt(normalize(batchData[k,:,:,0])) + imwrite(np.uint8(255*np.concatenate((im,np.concatenate((pm,gt),axis=1)),axis=1)),'%s/I%05d.png' % (outPMPath,i-j+k+1)) + + + # -------------------------------------------------- + # save hyper-parameters, clean-up + # -------------------------------------------------- + + saveData(UNet2D.hp,pathjoin(modelPath,'hp.data')) + + trainWriter.close() + validWriter.close() + sess.close() + + def deploy(imPath,nImages,modelPath,pmPath,gpuIndex,pmIndex): + os.environ['CUDA_VISIBLE_DEVICES']= '%d' % gpuIndex + + variablesPath = pathjoin(modelPath,'model.ckpt') + outPMPath = pmPath + + hp = loadData(pathjoin(modelPath,'hp.data')) + UNet2D.setupWithHP(hp) + + batchSize = UNet2D.hp['batchSize'] + imSize = UNet2D.hp['imSize'] + nChannels = UNet2D.hp['nChannels'] + nClasses = UNet2D.hp['nClasses'] + + # -------------------------------------------------- + # data + # -------------------------------------------------- + + Data = np.zeros((nImages,imSize,imSize,nChannels)) + + datasetMean = loadData(pathjoin(modelPath,'datasetMean.data')) + datasetStDev = loadData(pathjoin(modelPath,'datasetStDev.data')) + + for iSample in range(0, nImages): + path = '%s/I%05d_Img.tif' % (imPath,iSample) + im = im2double(tifread(path)) + Data[iSample,:,:,0] = (im-datasetMean)/datasetStDev + + # -------------------------------------------------- + # session + # -------------------------------------------------- + + saver = tf.train.Saver() + sess = tf.Session(config=tf.ConfigProto(allow_soft_placement=True)) # config parameter needed to save variables when using GPU + + saver.restore(sess, variablesPath) + print("Model restored.") + + # -------------------------------------------------- + # deploy + # -------------------------------------------------- + + batchData = np.zeros((batchSize,imSize,imSize,nChannels)) + + if not os.path.exists(outPMPath): + os.makedirs(outPMPath) + + for i in range(nImages): + print(i,nImages) + + j = np.mod(i,batchSize) + + batchData[j,:,:,:] = Data[i,:,:,:] + + if j == batchSize-1 or i == nImages-1: + + output = sess.run(UNet2D.nn,feed_dict={UNet2D.tfData: batchData, UNet2D.tfTraining: 0}) + + for k in range(j+1): + pm = output[k,:,:,pmIndex] + im = np.sqrt(normalize(batchData[k,:,:,0])) + # imwrite(np.uint8(255*np.concatenate((im,pm),axis=1)),'%s/I%05d.png' % (outPMPath,i-j+k+1)) + imwrite(np.uint8(255*im),'%s/I%05d_Im.png' % (outPMPath,i-j+k+1)) + imwrite(np.uint8(255*pm),'%s/I%05d_PM.png' % (outPMPath,i-j+k+1)) + + + # -------------------------------------------------- + # clean-up + # -------------------------------------------------- + + sess.close() + + def singleImageInferenceSetup(modelPath,gpuIndex): + os.environ['CUDA_VISIBLE_DEVICES']= '%d' % gpuIndex + + variablesPath = pathjoin(modelPath,'model.ckpt') + + hp = loadData(pathjoin(modelPath,'hp.data')) + UNet2D.setupWithHP(hp) + + UNet2D.DatasetMean = loadData(pathjoin(modelPath,'datasetMean.data')) + UNet2D.DatasetStDev = loadData(pathjoin(modelPath,'datasetStDev.data')) + print(UNet2D.DatasetMean) + print(UNet2D.DatasetStDev) + + # -------------------------------------------------- + # session + # -------------------------------------------------- + + saver = tf.train.Saver() + UNet2D.Session = tf.Session(config=tf.ConfigProto(allow_soft_placement=True)) # config parameter needed to save variables when using GPU + + saver.restore(UNet2D.Session, variablesPath) + print("Model restored.") + + def singleImageInferenceCleanup(): + UNet2D.Session.close() + + def singleImageInference(image,mode,pmIndex): + print('Inference...') + + batchSize = UNet2D.hp['batchSize'] + imSize = UNet2D.hp['imSize'] + nChannels = UNet2D.hp['nChannels'] + + PI2D.setup(image,imSize,int(imSize/8),mode) + PI2D.createOutput(nChannels) + + batchData = np.zeros((batchSize,imSize,imSize,nChannels)) + for i in range(PI2D.NumPatches): + j = np.mod(i,batchSize) + batchData[j,:,:,0] = (PI2D.getPatch(i)-UNet2D.DatasetMean)/UNet2D.DatasetStDev + if j == batchSize-1 or i == PI2D.NumPatches-1: + output = UNet2D.Session.run(UNet2D.nn,feed_dict={UNet2D.tfData: batchData, UNet2D.tfTraining: 0}) + for k in range(j+1): + pm = output[k,:,:,pmIndex] + PI2D.patchOutput(i-j+k,pm) + # PI2D.patchOutput(i-j+k,normalize(imgradmag(PI2D.getPatch(i-j+k),1))) + + return PI2D.getValidOutput() + + +if __name__ == '__main__': + logPath = 'D:\\LSP\\Sinem\\fromOlympus\\TFLogsssssssss' + modelPath = 'D:\\LSP\\UNet\\tonsil20x1bin1chan\\TFModel - 3class 16 kernels 5ks 2 layers' + pmPath = 'D:\\LSP\\Sinem\\fromOlympus\\TFProbMaps' + + + # ----- test 1 ----- + + # imPath = 'D:\\LSP\\Sinem\\trainingSetContours\\trainingSetSmallLarge' + # # UNet2D.setup(128,1,2,8,2,2,3,1,0.1,2,8) + # # UNet2D.train(imPath,logPath,modelPath,pmPath,500,100,40,True,20000,1,0) + # UNet2D.setup(128, 1, 2, 12, 2, 2, 3, 4, 0.1, 4, 8) + # UNet2D.train(imPath, logPath, modelPath, pmPath, 1600, 400, 500, False, 150000, 1, 0) + # UNet2D.deploy(imPath,100,modelPath,pmPath,1,0) + + # I = im2double(tifread('/home/mc457/files/CellBiology/IDAC/Marcelo/Etc/UNetTestSets/SinemSaka_NucleiSegmentation_SingleImageInferenceTest3.tif')) + # UNet2D.singleImageInferenceSetup(modelPath,0) + # J = UNet2D.singleImageInference(I,'accumulate',0) + # UNet2D.singleImageInferenceCleanup() + # # imshowlist([I,J]) + # # sys.exit(0) + # # tifwrite(np.uint8(255*I),'/home/mc457/Workspace/I1.tif') + # # tifwrite(np.uint8(255*J),'/home/mc457/Workspace/I2.tif') + # K = np.zeros((2,I.shape[0],I.shape[1])) + # K[0,:,:] = I + # K[1,:,:] = J + # tifwrite(np.uint8(255*K),'/home/mc457/Workspace/Sinem_NucSeg.tif') + + UNet2D.singleImageInferenceSetup(modelPath, 1) + imagePath ='Y:/sorger/data/RareCyte/Clarence/NKI_TMA' + sampleList = glob.glob(imagePath + '/ZTMA_18_810*') + dapiChannel = 0 + for iSample in sampleList: + # fileList = glob.glob(iSample + '//dearray//*.tif') + fileList = [x for x in glob.glob(iSample + '/dearray/*.tif') if x != (iSample+'/dearray\\TMA_MAP.tif')] + print(fileList) + for iFile in fileList: + fileName = os.path.basename(iFile) + fileNamePrefix = fileName.split(os.extsep, 1) + I = tifffile.imread(iFile, key=dapiChannel) + I = im2double(sk.rescale_intensity(I, in_range=(np.min(I), np.max(I)), out_range=(0, 64424))) + # I=np.moveaxis(I,0,-1) + # I=I[:,:,0] + hsize = int((float(I.shape[0])*float(1))) + vsize = int((float(I.shape[1])*float(1))) + I = resize(I,(hsize,vsize)) + #I = im2double(tifread('D:\\LSP\\cycif\\Unet\\Caitlin\\E - 04(fld 8 wv UV - DAPI)downsampled.tif')) + outputPath = iSample + '//prob_maps' + if not os.path.exists(outputPath): + os.makedirs(outputPath) + K = np.zeros((2,I.shape[0],I.shape[1])) + contours = UNet2D.singleImageInference(I,'accumulate',1) + K[1,:,:] = I + K[0,:,:] = contours + tifwrite(np.uint8(255 * K), + outputPath + '//' + fileNamePrefix[0] + '_ContoursPM_' + str(dapiChannel + 1) + '.tif') + del K + K = np.zeros((1, I.shape[0], I.shape[1])) + nuclei = UNet2D.singleImageInference(I,'accumulate',2) + K[0, :, :] = nuclei + tifwrite(np.uint8(255 * K), + outputPath + '//' + fileNamePrefix[0] + '_NucleiPM_' + str(dapiChannel + 1) + '.tif') + del K + UNet2D.singleImageInferenceCleanup() + + + # ----- test 2 ----- + + # imPath = '/home/mc457/files/CellBiology/IDAC/Marcelo/Etc/UNetTestSets/ClarenceYapp_NucleiSegmentation' + # UNet2D.setup(128,1,2,8,2,2,3,1,0.1,3,4) + # UNet2D.train(imPath,logPath,modelPath,pmPath,800,100,100,False,10,1) + # UNet2D.deploy(imPath,100,modelPath,pmPath,1) + + + # ----- test 3 ----- + + # imPath = '/home/mc457/files/CellBiology/IDAC/Marcelo/Etc/UNetTestSets/CarmanLi_CellTypeSegmentation' + # # UNet2D.setup(256,1,2,8,2,2,3,1,0.1,3,4) + # # UNet2D.train(imPath,logPath,modelPath,pmPath,1400,100,164,False,10000,1) + # UNet2D.deploy(imPath,164,modelPath,pmPath,1) + + + # ----- test 4 ----- + + # imPath = '/home/cicconet/Downloads/TrainSet1' + # UNet2D.setup(64,1,2,8,2,2,3,1,0.1,3,4) + # UNet2D.train(imPath,logPath,modelPath,pmPath,200,8,8,False,2000,1,0) + # # UNet2D.deploy(imPath,164,modelPath,pmPath,1) \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/batchUNet2DtCycif.py Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,553 @@ +import numpy as np +from scipy import misc +import tensorflow as tf +import shutil +import scipy.io as sio +import os,fnmatch,glob +import skimage.exposure as sk + +import sys +sys.path.insert(0, 'C:\\Users\\Clarence\\Documents\\UNet code\\ImageScience') +from toolbox.imtools import * +from toolbox.ftools import * +from toolbox.PartitionOfImage import PI2D + + +def concat3(lst): + return tf.concat(lst,3) + +class UNet2D: + hp = None # hyper-parameters + nn = None # network + tfTraining = None # if training or not (to handle batch norm) + tfData = None # data placeholder + Session = None + DatasetMean = 0 + DatasetStDev = 0 + + def setupWithHP(hp): + UNet2D.setup(hp['imSize'], + hp['nChannels'], + hp['nClasses'], + hp['nOut0'], + hp['featMapsFact'], + hp['downSampFact'], + hp['ks'], + hp['nExtraConvs'], + hp['stdDev0'], + hp['nLayers'], + hp['batchSize']) + + def setup(imSize,nChannels,nClasses,nOut0,featMapsFact,downSampFact,kernelSize,nExtraConvs,stdDev0,nDownSampLayers,batchSize): + UNet2D.hp = {'imSize':imSize, + 'nClasses':nClasses, + 'nChannels':nChannels, + 'nExtraConvs':nExtraConvs, + 'nLayers':nDownSampLayers, + 'featMapsFact':featMapsFact, + 'downSampFact':downSampFact, + 'ks':kernelSize, + 'nOut0':nOut0, + 'stdDev0':stdDev0, + 'batchSize':batchSize} + + nOutX = [UNet2D.hp['nChannels'],UNet2D.hp['nOut0']] + dsfX = [] + for i in range(UNet2D.hp['nLayers']): + nOutX.append(nOutX[-1]*UNet2D.hp['featMapsFact']) + dsfX.append(UNet2D.hp['downSampFact']) + + + # -------------------------------------------------- + # downsampling layer + # -------------------------------------------------- + + with tf.name_scope('placeholders'): + UNet2D.tfTraining = tf.placeholder(tf.bool, name='training') + UNet2D.tfData = tf.placeholder("float", shape=[None,UNet2D.hp['imSize'],UNet2D.hp['imSize'],UNet2D.hp['nChannels']],name='data') + + def down_samp_layer(data,index): + with tf.name_scope('ld%d' % index): + ldXWeights1 = tf.Variable(tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index], nOutX[index+1]], stddev=stdDev0),name='kernel1') + ldXWeightsExtra = [] + for i in range(nExtraConvs): + ldXWeightsExtra.append(tf.Variable(tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index+1], nOutX[index+1]], stddev=stdDev0),name='kernelExtra%d' % i)) + + c00 = tf.nn.conv2d(data, ldXWeights1, strides=[1, 1, 1, 1], padding='SAME') + for i in range(nExtraConvs): + c00 = tf.nn.conv2d(tf.nn.relu(c00), ldXWeightsExtra[i], strides=[1, 1, 1, 1], padding='SAME') + + ldXWeightsShortcut = tf.Variable(tf.truncated_normal([1, 1, nOutX[index], nOutX[index+1]], stddev=stdDev0),name='shortcutWeights') + shortcut = tf.nn.conv2d(data, ldXWeightsShortcut, strides=[1, 1, 1, 1], padding='SAME') + + bn = tf.layers.batch_normalization(tf.nn.relu(c00+shortcut), training=UNet2D.tfTraining) + + return tf.nn.max_pool(bn, ksize=[1, dsfX[index], dsfX[index], 1], strides=[1, dsfX[index], dsfX[index], 1], padding='SAME',name='maxpool') + + # -------------------------------------------------- + # bottom layer + # -------------------------------------------------- + + with tf.name_scope('lb'): + lbWeights1 = tf.Variable(tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[UNet2D.hp['nLayers']], nOutX[UNet2D.hp['nLayers']+1]], stddev=stdDev0),name='kernel1') + def lb(hidden): + return tf.nn.relu(tf.nn.conv2d(hidden, lbWeights1, strides=[1, 1, 1, 1], padding='SAME'),name='conv') + + # -------------------------------------------------- + # downsampling + # -------------------------------------------------- + + with tf.name_scope('downsampling'): + dsX = [] + dsX.append(UNet2D.tfData) + + for i in range(UNet2D.hp['nLayers']): + dsX.append(down_samp_layer(dsX[i],i)) + + b = lb(dsX[UNet2D.hp['nLayers']]) + + # -------------------------------------------------- + # upsampling layer + # -------------------------------------------------- + + def up_samp_layer(data,index): + with tf.name_scope('lu%d' % index): + luXWeights1 = tf.Variable(tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index+1], nOutX[index+2]], stddev=stdDev0),name='kernel1') + luXWeights2 = tf.Variable(tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index]+nOutX[index+1], nOutX[index+1]], stddev=stdDev0),name='kernel2') + luXWeightsExtra = [] + for i in range(nExtraConvs): + luXWeightsExtra.append(tf.Variable(tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index+1], nOutX[index+1]], stddev=stdDev0),name='kernel2Extra%d' % i)) + + outSize = UNet2D.hp['imSize'] + for i in range(index): + outSize /= dsfX[i] + outSize = int(outSize) + + outputShape = [UNet2D.hp['batchSize'],outSize,outSize,nOutX[index+1]] + us = tf.nn.relu(tf.nn.conv2d_transpose(data, luXWeights1, outputShape, strides=[1, dsfX[index], dsfX[index], 1], padding='SAME'),name='conv1') + cc = concat3([dsX[index],us]) + cv = tf.nn.relu(tf.nn.conv2d(cc, luXWeights2, strides=[1, 1, 1, 1], padding='SAME'),name='conv2') + for i in range(nExtraConvs): + cv = tf.nn.relu(tf.nn.conv2d(cv, luXWeightsExtra[i], strides=[1, 1, 1, 1], padding='SAME'),name='conv2Extra%d' % i) + return cv + + # -------------------------------------------------- + # final (top) layer + # -------------------------------------------------- + + with tf.name_scope('lt'): + ltWeights1 = tf.Variable(tf.truncated_normal([1, 1, nOutX[1], nClasses], stddev=stdDev0),name='kernel') + def lt(hidden): + return tf.nn.conv2d(hidden, ltWeights1, strides=[1, 1, 1, 1], padding='SAME',name='conv') + + + # -------------------------------------------------- + # upsampling + # -------------------------------------------------- + + with tf.name_scope('upsampling'): + usX = [] + usX.append(b) + + for i in range(UNet2D.hp['nLayers']): + usX.append(up_samp_layer(usX[i],UNet2D.hp['nLayers']-1-i)) + + t = lt(usX[UNet2D.hp['nLayers']]) + + + sm = tf.nn.softmax(t,-1) + UNet2D.nn = sm + + + def train(imPath,logPath,modelPath,pmPath,nTrain,nValid,nTest,restoreVariables,nSteps,gpuIndex,testPMIndex): + os.environ['CUDA_VISIBLE_DEVICES']= '%d' % gpuIndex + + outLogPath = logPath + trainWriterPath = pathjoin(logPath,'Train') + validWriterPath = pathjoin(logPath,'Valid') + outModelPath = pathjoin(modelPath,'model.ckpt') + outPMPath = pmPath + + batchSize = UNet2D.hp['batchSize'] + imSize = UNet2D.hp['imSize'] + nChannels = UNet2D.hp['nChannels'] + nClasses = UNet2D.hp['nClasses'] + + # -------------------------------------------------- + # data + # -------------------------------------------------- + + Train = np.zeros((nTrain,imSize,imSize,nChannels)) + Valid = np.zeros((nValid,imSize,imSize,nChannels)) + Test = np.zeros((nTest,imSize,imSize,nChannels)) + LTrain = np.zeros((nTrain,imSize,imSize,nClasses)) + LValid = np.zeros((nValid,imSize,imSize,nClasses)) + LTest = np.zeros((nTest,imSize,imSize,nClasses)) + + print('loading data, computing mean / st dev') + if not os.path.exists(modelPath): + os.makedirs(modelPath) + if restoreVariables: + datasetMean = loadData(pathjoin(modelPath,'datasetMean.data')) + datasetStDev = loadData(pathjoin(modelPath,'datasetStDev.data')) + else: + datasetMean = 0 + datasetStDev = 0 + for iSample in range(nTrain+nValid+nTest): + I = im2double(tifread('%s/I%05d_Img.tif' % (imPath,iSample))) + datasetMean += np.mean(I) + datasetStDev += np.std(I) + datasetMean /= (nTrain+nValid+nTest) + datasetStDev /= (nTrain+nValid+nTest) + saveData(datasetMean, pathjoin(modelPath,'datasetMean.data')) + saveData(datasetStDev, pathjoin(modelPath,'datasetStDev.data')) + + perm = np.arange(nTrain+nValid+nTest) + np.random.shuffle(perm) + + for iSample in range(0, nTrain): + path = '%s/I%05d_Img.tif' % (imPath,perm[iSample]) + im = im2double(tifread(path)) + Train[iSample,:,:,0] = (im-datasetMean)/datasetStDev + path = '%s/I%05d_Ant.tif' % (imPath,perm[iSample]) + im = tifread(path) + for i in range(nClasses): + LTrain[iSample,:,:,i] = (im == i+1) + + for iSample in range(0, nValid): + path = '%s/I%05d_Img.tif' % (imPath,perm[nTrain+iSample]) + im = im2double(tifread(path)) + Valid[iSample,:,:,0] = (im-datasetMean)/datasetStDev + path = '%s/I%05d_Ant.tif' % (imPath,perm[nTrain+iSample]) + im = tifread(path) + for i in range(nClasses): + LValid[iSample,:,:,i] = (im == i+1) + + for iSample in range(0, nTest): + path = '%s/I%05d_Img.tif' % (imPath,perm[nTrain+nValid+iSample]) + im = im2double(tifread(path)) + Test[iSample,:,:,0] = (im-datasetMean)/datasetStDev + path = '%s/I%05d_Ant.tif' % (imPath,perm[nTrain+nValid+iSample]) + im = tifread(path) + for i in range(nClasses): + LTest[iSample,:,:,i] = (im == i+1) + + # -------------------------------------------------- + # optimization + # -------------------------------------------------- + + tfLabels = tf.placeholder("float", shape=[None,imSize,imSize,nClasses],name='labels') + + globalStep = tf.Variable(0,trainable=False) + learningRate0 = 0.01 + decaySteps = 1000 + decayRate = 0.95 + learningRate = tf.train.exponential_decay(learningRate0,globalStep,decaySteps,decayRate,staircase=True) + + with tf.name_scope('optim'): + loss = tf.reduce_mean(-tf.reduce_sum(tf.multiply(tfLabels,tf.log(UNet2D.nn)),3)) + updateOps = tf.get_collection(tf.GraphKeys.UPDATE_OPS) + # optimizer = tf.train.MomentumOptimizer(1e-3,0.9) + optimizer = tf.train.MomentumOptimizer(learningRate,0.9) + # optimizer = tf.train.GradientDescentOptimizer(learningRate) + with tf.control_dependencies(updateOps): + optOp = optimizer.minimize(loss,global_step=globalStep) + + with tf.name_scope('eval'): + error = [] + for iClass in range(nClasses): + labels0 = tf.reshape(tf.to_int32(tf.slice(tfLabels,[0,0,0,iClass],[-1,-1,-1,1])),[batchSize,imSize,imSize]) + predict0 = tf.reshape(tf.to_int32(tf.equal(tf.argmax(UNet2D.nn,3),iClass)),[batchSize,imSize,imSize]) + correct = tf.multiply(labels0,predict0) + nCorrect0 = tf.reduce_sum(correct) + nLabels0 = tf.reduce_sum(labels0) + error.append(1-tf.to_float(nCorrect0)/tf.to_float(nLabels0)) + errors = tf.tuple(error) + + # -------------------------------------------------- + # inspection + # -------------------------------------------------- + + with tf.name_scope('scalars'): + tf.summary.scalar('avg_cross_entropy', loss) + for iClass in range(nClasses): + tf.summary.scalar('avg_pixel_error_%d' % iClass, error[iClass]) + tf.summary.scalar('learning_rate', learningRate) + with tf.name_scope('images'): + split0 = tf.slice(UNet2D.nn,[0,0,0,0],[-1,-1,-1,1]) + split1 = tf.slice(UNet2D.nn,[0,0,0,1],[-1,-1,-1,1]) + if nClasses > 2: + split2 = tf.slice(UNet2D.nn,[0,0,0,2],[-1,-1,-1,1]) + tf.summary.image('pm0',split0) + tf.summary.image('pm1',split1) + if nClasses > 2: + tf.summary.image('pm2',split2) + merged = tf.summary.merge_all() + + + # -------------------------------------------------- + # session + # -------------------------------------------------- + + saver = tf.train.Saver() + sess = tf.Session(config=tf.ConfigProto(allow_soft_placement=True)) # config parameter needed to save variables when using GPU + + if os.path.exists(outLogPath): + shutil.rmtree(outLogPath) + trainWriter = tf.summary.FileWriter(trainWriterPath, sess.graph) + validWriter = tf.summary.FileWriter(validWriterPath, sess.graph) + + if restoreVariables: + saver.restore(sess, outModelPath) + print("Model restored.") + else: + sess.run(tf.global_variables_initializer()) + + # -------------------------------------------------- + # train + # -------------------------------------------------- + + batchData = np.zeros((batchSize,imSize,imSize,nChannels)) + batchLabels = np.zeros((batchSize,imSize,imSize,nClasses)) + for i in range(nSteps): + # train + + perm = np.arange(nTrain) + np.random.shuffle(perm) + + for j in range(batchSize): + batchData[j,:,:,:] = Train[perm[j],:,:,:] + batchLabels[j,:,:,:] = LTrain[perm[j],:,:,:] + + summary,_ = sess.run([merged,optOp],feed_dict={UNet2D.tfData: batchData, tfLabels: batchLabels, UNet2D.tfTraining: 1}) + trainWriter.add_summary(summary, i) + + # validation + + perm = np.arange(nValid) + np.random.shuffle(perm) + + for j in range(batchSize): + batchData[j,:,:,:] = Valid[perm[j],:,:,:] + batchLabels[j,:,:,:] = LValid[perm[j],:,:,:] + + summary, es = sess.run([merged, errors],feed_dict={UNet2D.tfData: batchData, tfLabels: batchLabels, UNet2D.tfTraining: 0}) + validWriter.add_summary(summary, i) + + e = np.mean(es) + print('step %05d, e: %f' % (i,e)) + + if i == 0: + if restoreVariables: + lowestError = e + else: + lowestError = np.inf + + if np.mod(i,100) == 0 and e < lowestError: + lowestError = e + print("Model saved in file: %s" % saver.save(sess, outModelPath)) + + + # -------------------------------------------------- + # test + # -------------------------------------------------- + + if not os.path.exists(outPMPath): + os.makedirs(outPMPath) + + for i in range(nTest): + j = np.mod(i,batchSize) + + batchData[j,:,:,:] = Test[i,:,:,:] + batchLabels[j,:,:,:] = LTest[i,:,:,:] + + if j == batchSize-1 or i == nTest-1: + + output = sess.run(UNet2D.nn,feed_dict={UNet2D.tfData: batchData, tfLabels: batchLabels, UNet2D.tfTraining: 0}) + + for k in range(j+1): + pm = output[k,:,:,testPMIndex] + gt = batchLabels[k,:,:,testPMIndex] + im = np.sqrt(normalize(batchData[k,:,:,0])) + imwrite(np.uint8(255*np.concatenate((im,np.concatenate((pm,gt),axis=1)),axis=1)),'%s/I%05d.png' % (outPMPath,i-j+k+1)) + + + # -------------------------------------------------- + # save hyper-parameters, clean-up + # -------------------------------------------------- + + saveData(UNet2D.hp,pathjoin(modelPath,'hp.data')) + + trainWriter.close() + validWriter.close() + sess.close() + + def deploy(imPath,nImages,modelPath,pmPath,gpuIndex,pmIndex): + os.environ['CUDA_VISIBLE_DEVICES']= '%d' % gpuIndex + + variablesPath = pathjoin(modelPath,'model.ckpt') + outPMPath = pmPath + + hp = loadData(pathjoin(modelPath,'hp.data')) + UNet2D.setupWithHP(hp) + + batchSize = UNet2D.hp['batchSize'] + imSize = UNet2D.hp['imSize'] + nChannels = UNet2D.hp['nChannels'] + nClasses = UNet2D.hp['nClasses'] + + # -------------------------------------------------- + # data + # -------------------------------------------------- + + Data = np.zeros((nImages,imSize,imSize,nChannels)) + + datasetMean = loadData(pathjoin(modelPath,'datasetMean.data')) + datasetStDev = loadData(pathjoin(modelPath,'datasetStDev.data')) + + for iSample in range(0, nImages): + path = '%s/I%05d_Img.tif' % (imPath,iSample) + im = im2double(tifread(path)) + Data[iSample,:,:,0] = (im-datasetMean)/datasetStDev + + # -------------------------------------------------- + # session + # -------------------------------------------------- + + saver = tf.train.Saver() + sess = tf.Session(config=tf.ConfigProto(allow_soft_placement=True)) # config parameter needed to save variables when using GPU + + saver.restore(sess, variablesPath) + print("Model restored.") + + # -------------------------------------------------- + # deploy + # -------------------------------------------------- + + batchData = np.zeros((batchSize,imSize,imSize,nChannels)) + + if not os.path.exists(outPMPath): + os.makedirs(outPMPath) + + for i in range(nImages): + print(i,nImages) + + j = np.mod(i,batchSize) + + batchData[j,:,:,:] = Data[i,:,:,:] + + if j == batchSize-1 or i == nImages-1: + + output = sess.run(UNet2D.nn,feed_dict={UNet2D.tfData: batchData, UNet2D.tfTraining: 0}) + + for k in range(j+1): + pm = output[k,:,:,pmIndex] + im = np.sqrt(normalize(batchData[k,:,:,0])) + # imwrite(np.uint8(255*np.concatenate((im,pm),axis=1)),'%s/I%05d.png' % (outPMPath,i-j+k+1)) + imwrite(np.uint8(255*im),'%s/I%05d_Im.png' % (outPMPath,i-j+k+1)) + imwrite(np.uint8(255*pm),'%s/I%05d_PM.png' % (outPMPath,i-j+k+1)) + + + # -------------------------------------------------- + # clean-up + # -------------------------------------------------- + + sess.close() + + def singleImageInferenceSetup(modelPath,gpuIndex): + #os.environ['CUDA_VISIBLE_DEVICES']= '%d' % gpuIndex + + variablesPath = pathjoin(modelPath,'model.ckpt') + + hp = loadData(pathjoin(modelPath,'hp.data')) + UNet2D.setupWithHP(hp) + + UNet2D.DatasetMean = loadData(pathjoin(modelPath,'datasetMean.data')) + UNet2D.DatasetStDev = loadData(pathjoin(modelPath,'datasetStDev.data')) + print(UNet2D.DatasetMean) + print(UNet2D.DatasetStDev) + + # -------------------------------------------------- + # session + # -------------------------------------------------- + + saver = tf.train.Saver() + UNet2D.Session = tf.Session(config=tf.ConfigProto(allow_soft_placement=True)) # config parameter needed to save variables when using GPU + + saver.restore(UNet2D.Session, variablesPath) + print("Model restored.") + + def singleImageInferenceCleanup(): + UNet2D.Session.close() + + def singleImageInference(image,mode,pmIndex): + print('Inference...') + + batchSize = UNet2D.hp['batchSize'] + imSize = UNet2D.hp['imSize'] + nChannels = UNet2D.hp['nChannels'] + + PI2D.setup(image,imSize,int(imSize/8),mode) + PI2D.createOutput(nChannels) + + batchData = np.zeros((batchSize,imSize,imSize,nChannels)) + for i in range(PI2D.NumPatches): + j = np.mod(i,batchSize) + batchData[j,:,:,0] = (PI2D.getPatch(i)-UNet2D.DatasetMean)/UNet2D.DatasetStDev + if j == batchSize-1 or i == PI2D.NumPatches-1: + output = UNet2D.Session.run(UNet2D.nn,feed_dict={UNet2D.tfData: batchData, UNet2D.tfTraining: 0}) + for k in range(j+1): + pm = output[k,:,:,pmIndex] + PI2D.patchOutput(i-j+k,pm) + # PI2D.patchOutput(i-j+k,normalize(imgradmag(PI2D.getPatch(i-j+k),1))) + + return PI2D.getValidOutput() + + +if __name__ == '__main__': + logPath = 'C://Users//Clarence//Documents//UNet code//TFLogs' + modelPath = 'D:\\LSP\\UNet\\tonsil20x1bin1chan\\TFModel - 3class 16 kernels 5ks 2 layers' + pmPath = 'C://Users//Clarence//Documents//UNet code//TFProbMaps' + + + + UNet2D.singleImageInferenceSetup(modelPath, 0) + imagePath = 'D:\\LSP\\cycif\\testsets' + sampleList = glob.glob(imagePath + '//exemplar-001*') + dapiChannel = 0 + dsFactor = 1 + for iSample in sampleList: + fileList = glob.glob(iSample + '//registration//*.tif') + print(fileList) + for iFile in fileList: + fileName = os.path.basename(iFile) + fileNamePrefix = fileName.split(os.extsep, 1) + I = tifffile.imread(iFile, key=dapiChannel) + rawI = I + hsize = int((float(I.shape[0])*float(dsFactor))) + vsize = int((float(I.shape[1])*float(dsFactor))) + I = resize(I,(hsize,vsize)) + I = im2double(sk.rescale_intensity(I, in_range=(np.min(I), np.max(I)), out_range=(0, 0.983))) + rawI = im2double(rawI)/np.max(im2double(rawI)) + outputPath = iSample + '//prob_maps' + if not os.path.exists(outputPath): + os.makedirs(outputPath) + K = np.zeros((2,rawI.shape[0],rawI.shape[1])) + contours = UNet2D.singleImageInference(I,'accumulate',1) + hsize = int((float(I.shape[0]) * float(1/dsFactor))) + vsize = int((float(I.shape[1]) * float(1/dsFactor))) + contours = resize(contours, (rawI.shape[0], rawI.shape[1])) + K[1,:,:] = rawI + K[0,:,:] = contours + tifwrite(np.uint8(255 * K), + outputPath + '//' + fileNamePrefix[0] + '_ContoursPM_' + str(dapiChannel + 1) + '.tif') + del K + K = np.zeros((1, rawI.shape[0], rawI.shape[1])) + nuclei = UNet2D.singleImageInference(I,'accumulate',2) + nuclei = resize(nuclei, (rawI.shape[0], rawI.shape[1])) + K[0, :, :] = nuclei + tifwrite(np.uint8(255 * K), + outputPath + '//' + fileNamePrefix[0] + '_NucleiPM_' + str(dapiChannel + 1) + '.tif') + del K + UNet2D.singleImageInferenceCleanup() +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/batchUnMicst.py Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,588 @@ +import numpy as np +from scipy import misc +import tensorflow as tf +import shutil +import scipy.io as sio +import os, fnmatch, PIL, glob +import skimage.exposure as sk +import argparse + +import sys + +sys.path.insert(0, 'C:\\Users\\Public\\Documents\\ImageScience') +from toolbox.imtools import * +from toolbox.ftools import * +from toolbox.PartitionOfImage import PI2D + + +def concat3(lst): + return tf.concat(lst, 3) + + +class UNet2D: + hp = None # hyper-parameters + nn = None # network + tfTraining = None # if training or not (to handle batch norm) + tfData = None # data placeholder + Session = None + DatasetMean = 0 + DatasetStDev = 0 + + def setupWithHP(hp): + UNet2D.setup(hp['imSize'], + hp['nChannels'], + hp['nClasses'], + hp['nOut0'], + hp['featMapsFact'], + hp['downSampFact'], + hp['ks'], + hp['nExtraConvs'], + hp['stdDev0'], + hp['nLayers'], + hp['batchSize']) + + def setup(imSize, nChannels, nClasses, nOut0, featMapsFact, downSampFact, kernelSize, nExtraConvs, stdDev0, + nDownSampLayers, batchSize): + UNet2D.hp = {'imSize': imSize, + 'nClasses': nClasses, + 'nChannels': nChannels, + 'nExtraConvs': nExtraConvs, + 'nLayers': nDownSampLayers, + 'featMapsFact': featMapsFact, + 'downSampFact': downSampFact, + 'ks': kernelSize, + 'nOut0': nOut0, + 'stdDev0': stdDev0, + 'batchSize': batchSize} + + nOutX = [UNet2D.hp['nChannels'], UNet2D.hp['nOut0']] + dsfX = [] + for i in range(UNet2D.hp['nLayers']): + nOutX.append(nOutX[-1] * UNet2D.hp['featMapsFact']) + dsfX.append(UNet2D.hp['downSampFact']) + + # -------------------------------------------------- + # downsampling layer + # -------------------------------------------------- + + with tf.name_scope('placeholders'): + UNet2D.tfTraining = tf.placeholder(tf.bool, name='training') + UNet2D.tfData = tf.placeholder("float", shape=[None, UNet2D.hp['imSize'], UNet2D.hp['imSize'], + UNet2D.hp['nChannels']], name='data') + + def down_samp_layer(data, index): + with tf.name_scope('ld%d' % index): + ldXWeights1 = tf.Variable( + tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index], nOutX[index + 1]], + stddev=stdDev0), name='kernel1') + ldXWeightsExtra = [] + for i in range(nExtraConvs): + ldXWeightsExtra.append(tf.Variable( + tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index + 1], nOutX[index + 1]], + stddev=stdDev0), name='kernelExtra%d' % i)) + + c00 = tf.nn.conv2d(data, ldXWeights1, strides=[1, 1, 1, 1], padding='SAME') + for i in range(nExtraConvs): + c00 = tf.nn.conv2d(tf.nn.relu(c00), ldXWeightsExtra[i], strides=[1, 1, 1, 1], padding='SAME') + + ldXWeightsShortcut = tf.Variable( + tf.truncated_normal([1, 1, nOutX[index], nOutX[index + 1]], stddev=stdDev0), name='shortcutWeights') + shortcut = tf.nn.conv2d(data, ldXWeightsShortcut, strides=[1, 1, 1, 1], padding='SAME') + + bn = tf.layers.batch_normalization(tf.nn.relu(c00 + shortcut), training=UNet2D.tfTraining) + + return tf.nn.max_pool(bn, ksize=[1, dsfX[index], dsfX[index], 1], + strides=[1, dsfX[index], dsfX[index], 1], padding='SAME', name='maxpool') + + # -------------------------------------------------- + # bottom layer + # -------------------------------------------------- + + with tf.name_scope('lb'): + lbWeights1 = tf.Variable(tf.truncated_normal( + [UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[UNet2D.hp['nLayers']], nOutX[UNet2D.hp['nLayers'] + 1]], + stddev=stdDev0), name='kernel1') + + def lb(hidden): + return tf.nn.relu(tf.nn.conv2d(hidden, lbWeights1, strides=[1, 1, 1, 1], padding='SAME'), name='conv') + + # -------------------------------------------------- + # downsampling + # -------------------------------------------------- + + with tf.name_scope('downsampling'): + dsX = [] + dsX.append(UNet2D.tfData) + + for i in range(UNet2D.hp['nLayers']): + dsX.append(down_samp_layer(dsX[i], i)) + + b = lb(dsX[UNet2D.hp['nLayers']]) + + # -------------------------------------------------- + # upsampling layer + # -------------------------------------------------- + + def up_samp_layer(data, index): + with tf.name_scope('lu%d' % index): + luXWeights1 = tf.Variable( + tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index + 1], nOutX[index + 2]], + stddev=stdDev0), name='kernel1') + luXWeights2 = tf.Variable(tf.truncated_normal( + [UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index] + nOutX[index + 1], nOutX[index + 1]], + stddev=stdDev0), name='kernel2') + luXWeightsExtra = [] + for i in range(nExtraConvs): + luXWeightsExtra.append(tf.Variable( + tf.truncated_normal([UNet2D.hp['ks'], UNet2D.hp['ks'], nOutX[index + 1], nOutX[index + 1]], + stddev=stdDev0), name='kernel2Extra%d' % i)) + + outSize = UNet2D.hp['imSize'] + for i in range(index): + outSize /= dsfX[i] + outSize = int(outSize) + + outputShape = [UNet2D.hp['batchSize'], outSize, outSize, nOutX[index + 1]] + us = tf.nn.relu( + tf.nn.conv2d_transpose(data, luXWeights1, outputShape, strides=[1, dsfX[index], dsfX[index], 1], + padding='SAME'), name='conv1') + cc = concat3([dsX[index], us]) + cv = tf.nn.relu(tf.nn.conv2d(cc, luXWeights2, strides=[1, 1, 1, 1], padding='SAME'), name='conv2') + for i in range(nExtraConvs): + cv = tf.nn.relu(tf.nn.conv2d(cv, luXWeightsExtra[i], strides=[1, 1, 1, 1], padding='SAME'), + name='conv2Extra%d' % i) + return cv + + # -------------------------------------------------- + # final (top) layer + # -------------------------------------------------- + + with tf.name_scope('lt'): + ltWeights1 = tf.Variable(tf.truncated_normal([1, 1, nOutX[1], nClasses], stddev=stdDev0), name='kernel') + + def lt(hidden): + return tf.nn.conv2d(hidden, ltWeights1, strides=[1, 1, 1, 1], padding='SAME', name='conv') + + # -------------------------------------------------- + # upsampling + # -------------------------------------------------- + + with tf.name_scope('upsampling'): + usX = [] + usX.append(b) + + for i in range(UNet2D.hp['nLayers']): + usX.append(up_samp_layer(usX[i], UNet2D.hp['nLayers'] - 1 - i)) + + t = lt(usX[UNet2D.hp['nLayers']]) + + sm = tf.nn.softmax(t, -1) + UNet2D.nn = sm + + def train(imPath, logPath, modelPath, pmPath, nTrain, nValid, nTest, restoreVariables, nSteps, gpuIndex, + testPMIndex): + os.environ['CUDA_VISIBLE_DEVICES'] = '%d' % gpuIndex + + outLogPath = logPath + trainWriterPath = pathjoin(logPath, 'Train') + validWriterPath = pathjoin(logPath, 'Valid') + outModelPath = pathjoin(modelPath, 'model.ckpt') + outPMPath = pmPath + + batchSize = UNet2D.hp['batchSize'] + imSize = UNet2D.hp['imSize'] + nChannels = UNet2D.hp['nChannels'] + nClasses = UNet2D.hp['nClasses'] + + # -------------------------------------------------- + # data + # -------------------------------------------------- + + Train = np.zeros((nTrain, imSize, imSize, nChannels)) + Valid = np.zeros((nValid, imSize, imSize, nChannels)) + Test = np.zeros((nTest, imSize, imSize, nChannels)) + LTrain = np.zeros((nTrain, imSize, imSize, nClasses)) + LValid = np.zeros((nValid, imSize, imSize, nClasses)) + LTest = np.zeros((nTest, imSize, imSize, nClasses)) + + print('loading data, computing mean / st dev') + if not os.path.exists(modelPath): + os.makedirs(modelPath) + if restoreVariables: + datasetMean = loadData(pathjoin(modelPath, 'datasetMean.data')) + datasetStDev = loadData(pathjoin(modelPath, 'datasetStDev.data')) + else: + datasetMean = 0 + datasetStDev = 0 + for iSample in range(nTrain + nValid + nTest): + I = im2double(tifread('%s/I%05d_Img.tif' % (imPath, iSample))) + datasetMean += np.mean(I) + datasetStDev += np.std(I) + datasetMean /= (nTrain + nValid + nTest) + datasetStDev /= (nTrain + nValid + nTest) + saveData(datasetMean, pathjoin(modelPath, 'datasetMean.data')) + saveData(datasetStDev, pathjoin(modelPath, 'datasetStDev.data')) + + perm = np.arange(nTrain + nValid + nTest) + np.random.shuffle(perm) + + for iSample in range(0, nTrain): + path = '%s/I%05d_Img.tif' % (imPath, perm[iSample]) + im = im2double(tifread(path)) + Train[iSample, :, :, 0] = (im - datasetMean) / datasetStDev + path = '%s/I%05d_Ant.tif' % (imPath, perm[iSample]) + im = tifread(path) + for i in range(nClasses): + LTrain[iSample, :, :, i] = (im == i + 1) + + for iSample in range(0, nValid): + path = '%s/I%05d_Img.tif' % (imPath, perm[nTrain + iSample]) + im = im2double(tifread(path)) + Valid[iSample, :, :, 0] = (im - datasetMean) / datasetStDev + path = '%s/I%05d_Ant.tif' % (imPath, perm[nTrain + iSample]) + im = tifread(path) + for i in range(nClasses): + LValid[iSample, :, :, i] = (im == i + 1) + + for iSample in range(0, nTest): + path = '%s/I%05d_Img.tif' % (imPath, perm[nTrain + nValid + iSample]) + im = im2double(tifread(path)) + Test[iSample, :, :, 0] = (im - datasetMean) / datasetStDev + path = '%s/I%05d_Ant.tif' % (imPath, perm[nTrain + nValid + iSample]) + im = tifread(path) + for i in range(nClasses): + LTest[iSample, :, :, i] = (im == i + 1) + + # -------------------------------------------------- + # optimization + # -------------------------------------------------- + + tfLabels = tf.placeholder("float", shape=[None, imSize, imSize, nClasses], name='labels') + + globalStep = tf.Variable(0, trainable=False) + learningRate0 = 0.01 + decaySteps = 1000 + decayRate = 0.95 + learningRate = tf.train.exponential_decay(learningRate0, globalStep, decaySteps, decayRate, staircase=True) + + with tf.name_scope('optim'): + loss = tf.reduce_mean(-tf.reduce_sum(tf.multiply(tfLabels, tf.log(UNet2D.nn)), 3)) + updateOps = tf.get_collection(tf.GraphKeys.UPDATE_OPS) + # optimizer = tf.train.MomentumOptimizer(1e-3,0.9) + optimizer = tf.train.MomentumOptimizer(learningRate, 0.9) + # optimizer = tf.train.GradientDescentOptimizer(learningRate) + with tf.control_dependencies(updateOps): + optOp = optimizer.minimize(loss, global_step=globalStep) + + with tf.name_scope('eval'): + error = [] + for iClass in range(nClasses): + labels0 = tf.reshape(tf.to_int32(tf.slice(tfLabels, [0, 0, 0, iClass], [-1, -1, -1, 1])), + [batchSize, imSize, imSize]) + predict0 = tf.reshape(tf.to_int32(tf.equal(tf.argmax(UNet2D.nn, 3), iClass)), + [batchSize, imSize, imSize]) + correct = tf.multiply(labels0, predict0) + nCorrect0 = tf.reduce_sum(correct) + nLabels0 = tf.reduce_sum(labels0) + error.append(1 - tf.to_float(nCorrect0) / tf.to_float(nLabels0)) + errors = tf.tuple(error) + + # -------------------------------------------------- + # inspection + # -------------------------------------------------- + + with tf.name_scope('scalars'): + tf.summary.scalar('avg_cross_entropy', loss) + for iClass in range(nClasses): + tf.summary.scalar('avg_pixel_error_%d' % iClass, error[iClass]) + tf.summary.scalar('learning_rate', learningRate) + with tf.name_scope('images'): + split0 = tf.slice(UNet2D.nn, [0, 0, 0, 0], [-1, -1, -1, 1]) + split1 = tf.slice(UNet2D.nn, [0, 0, 0, 1], [-1, -1, -1, 1]) + if nClasses > 2: + split2 = tf.slice(UNet2D.nn, [0, 0, 0, 2], [-1, -1, -1, 1]) + tf.summary.image('pm0', split0) + tf.summary.image('pm1', split1) + if nClasses > 2: + tf.summary.image('pm2', split2) + merged = tf.summary.merge_all() + + # -------------------------------------------------- + # session + # -------------------------------------------------- + + saver = tf.train.Saver() + sess = tf.Session(config=tf.ConfigProto( + allow_soft_placement=True)) # config parameter needed to save variables when using GPU + + if os.path.exists(outLogPath): + shutil.rmtree(outLogPath) + trainWriter = tf.summary.FileWriter(trainWriterPath, sess.graph) + validWriter = tf.summary.FileWriter(validWriterPath, sess.graph) + + if restoreVariables: + saver.restore(sess, outModelPath) + print("Model restored.") + else: + sess.run(tf.global_variables_initializer()) + + # -------------------------------------------------- + # train + # -------------------------------------------------- + + batchData = np.zeros((batchSize, imSize, imSize, nChannels)) + batchLabels = np.zeros((batchSize, imSize, imSize, nClasses)) + for i in range(nSteps): + # train + + perm = np.arange(nTrain) + np.random.shuffle(perm) + + for j in range(batchSize): + batchData[j, :, :, :] = Train[perm[j], :, :, :] + batchLabels[j, :, :, :] = LTrain[perm[j], :, :, :] + + summary, _ = sess.run([merged, optOp], + feed_dict={UNet2D.tfData: batchData, tfLabels: batchLabels, UNet2D.tfTraining: 1}) + trainWriter.add_summary(summary, i) + + # validation + + perm = np.arange(nValid) + np.random.shuffle(perm) + + for j in range(batchSize): + batchData[j, :, :, :] = Valid[perm[j], :, :, :] + batchLabels[j, :, :, :] = LValid[perm[j], :, :, :] + + summary, es = sess.run([merged, errors], + feed_dict={UNet2D.tfData: batchData, tfLabels: batchLabels, UNet2D.tfTraining: 0}) + validWriter.add_summary(summary, i) + + e = np.mean(es) + print('step %05d, e: %f' % (i, e)) + + if i == 0: + if restoreVariables: + lowestError = e + else: + lowestError = np.inf + + if np.mod(i, 100) == 0 and e < lowestError: + lowestError = e + print("Model saved in file: %s" % saver.save(sess, outModelPath)) + + # -------------------------------------------------- + # test + # -------------------------------------------------- + + if not os.path.exists(outPMPath): + os.makedirs(outPMPath) + + for i in range(nTest): + j = np.mod(i, batchSize) + + batchData[j, :, :, :] = Test[i, :, :, :] + batchLabels[j, :, :, :] = LTest[i, :, :, :] + + if j == batchSize - 1 or i == nTest - 1: + + output = sess.run(UNet2D.nn, + feed_dict={UNet2D.tfData: batchData, tfLabels: batchLabels, UNet2D.tfTraining: 0}) + + for k in range(j + 1): + pm = output[k, :, :, testPMIndex] + gt = batchLabels[k, :, :, testPMIndex] + im = np.sqrt(normalize(batchData[k, :, :, 0])) + imwrite(np.uint8(255 * np.concatenate((im, np.concatenate((pm, gt), axis=1)), axis=1)), + '%s/I%05d.png' % (outPMPath, i - j + k + 1)) + + # -------------------------------------------------- + # save hyper-parameters, clean-up + # -------------------------------------------------- + + saveData(UNet2D.hp, pathjoin(modelPath, 'hp.data')) + + trainWriter.close() + validWriter.close() + sess.close() + + def deploy(imPath, nImages, modelPath, pmPath, gpuIndex, pmIndex): + os.environ['CUDA_VISIBLE_DEVICES'] = '%d' % gpuIndex + + variablesPath = pathjoin(modelPath, 'model.ckpt') + outPMPath = pmPath + + hp = loadData(pathjoin(modelPath, 'hp.data')) + UNet2D.setupWithHP(hp) + + batchSize = UNet2D.hp['batchSize'] + imSize = UNet2D.hp['imSize'] + nChannels = UNet2D.hp['nChannels'] + nClasses = UNet2D.hp['nClasses'] + + # -------------------------------------------------- + # data + # -------------------------------------------------- + + Data = np.zeros((nImages, imSize, imSize, nChannels)) + + datasetMean = loadData(pathjoin(modelPath, 'datasetMean.data')) + datasetStDev = loadData(pathjoin(modelPath, 'datasetStDev.data')) + + for iSample in range(0, nImages): + path = '%s/I%05d_Img.tif' % (imPath, iSample) + im = im2double(tifread(path)) + Data[iSample, :, :, 0] = (im - datasetMean) / datasetStDev + + # -------------------------------------------------- + # session + # -------------------------------------------------- + + saver = tf.train.Saver() + sess = tf.Session(config=tf.ConfigProto( + allow_soft_placement=True)) # config parameter needed to save variables when using GPU + + saver.restore(sess, variablesPath) + print("Model restored.") + + # -------------------------------------------------- + # deploy + # -------------------------------------------------- + + batchData = np.zeros((batchSize, imSize, imSize, nChannels)) + + if not os.path.exists(outPMPath): + os.makedirs(outPMPath) + + for i in range(nImages): + print(i, nImages) + + j = np.mod(i, batchSize) + + batchData[j, :, :, :] = Data[i, :, :, :] + + if j == batchSize - 1 or i == nImages - 1: + + output = sess.run(UNet2D.nn, feed_dict={UNet2D.tfData: batchData, UNet2D.tfTraining: 0}) + + for k in range(j + 1): + pm = output[k, :, :, pmIndex] + im = np.sqrt(normalize(batchData[k, :, :, 0])) + # imwrite(np.uint8(255*np.concatenate((im,pm),axis=1)),'%s/I%05d.png' % (outPMPath,i-j+k+1)) + imwrite(np.uint8(255 * im), '%s/I%05d_Im.png' % (outPMPath, i - j + k + 1)) + imwrite(np.uint8(255 * pm), '%s/I%05d_PM.png' % (outPMPath, i - j + k + 1)) + + # -------------------------------------------------- + # clean-up + # -------------------------------------------------- + + sess.close() + + def singleImageInferenceSetup(modelPath, gpuIndex): + os.environ['CUDA_VISIBLE_DEVICES'] = '%d' % gpuIndex + + variablesPath = pathjoin(modelPath, 'model.ckpt') + + hp = loadData(pathjoin(modelPath, 'hp.data')) + UNet2D.setupWithHP(hp) + + UNet2D.DatasetMean = loadData(pathjoin(modelPath, 'datasetMean.data')) + UNet2D.DatasetStDev = loadData(pathjoin(modelPath, 'datasetStDev.data')) + print(UNet2D.DatasetMean) + print(UNet2D.DatasetStDev) + + # -------------------------------------------------- + # session + # -------------------------------------------------- + + saver = tf.train.Saver() + UNet2D.Session = tf.Session(config=tf.ConfigProto( + allow_soft_placement=True)) # config parameter needed to save variables when using GPU + + saver.restore(UNet2D.Session, variablesPath) + print("Model restored.") + + def singleImageInferenceCleanup(): + UNet2D.Session.close() + + def singleImageInference(image, mode, pmIndex): + print('Inference...') + + batchSize = UNet2D.hp['batchSize'] + imSize = UNet2D.hp['imSize'] + nChannels = UNet2D.hp['nChannels'] + + PI2D.setup(image, imSize, int(imSize / 8), mode) + PI2D.createOutput(nChannels) + + batchData = np.zeros((batchSize, imSize, imSize, nChannels)) + for i in range(PI2D.NumPatches): + j = np.mod(i, batchSize) + batchData[j, :, :, 0] = (PI2D.getPatch(i) - UNet2D.DatasetMean) / UNet2D.DatasetStDev + if j == batchSize - 1 or i == PI2D.NumPatches - 1: + output = UNet2D.Session.run(UNet2D.nn, feed_dict={UNet2D.tfData: batchData, UNet2D.tfTraining: 0}) + for k in range(j + 1): + pm = output[k, :, :, pmIndex] + PI2D.patchOutput(i - j + k, pm) + # PI2D.patchOutput(i-j+k,normalize(imgradmag(PI2D.getPatch(i-j+k),1))) + + return PI2D.getValidOutput() + + +if __name__ == '__main__': + parser = argparse.ArgumentParser() + parser.add_argument("imagePath", help="path to the .tif file") + parser.add_argument("--channel", help="channel to perform inference on", type=int, default=0) + parser.add_argument("--TMA", help="specify if TMA", action="store_true") + parser.add_argument("--scalingFactor", help="factor by which to increase/decrease image size by", type=float, + default=1) + args = parser.parse_args() + + logPath = '' + modelPath = 'D:\\LSP\\UNet\\tonsil20x1bin1chan\\TFModel - 3class 16 kernels 5ks 2 layers' + pmPath = '' + + UNet2D.singleImageInferenceSetup(modelPath, 1) + imagePath = args.imagePath + sampleList = glob.glob(imagePath + '/exemplar*') + dapiChannel = args.channel + dsFactor = args.scalingFactor + for iSample in sampleList: + if args.TMA: + fileList = [x for x in glob.glob(iSample + '\\dearray\\*.tif') if x != (iSample + '\\dearray\\TMA_MAP.tif')] + print(iSample) + else: + fileList = glob.glob(iSample + '//registration//*ome.tif') + print(fileList) + for iFile in fileList: + fileName = os.path.basename(iFile) + fileNamePrefix = fileName.split(os.extsep, 1) + I = tifffile.imread(iFile, key=dapiChannel) + rawI = I + hsize = int((float(I.shape[0]) * float(dsFactor))) + vsize = int((float(I.shape[1]) * float(dsFactor))) + I = resize(I, (hsize, vsize)) + I = im2double(sk.rescale_intensity(I, in_range=(np.min(I), np.max(I)), out_range=(0, 0.983))) + rawI = im2double(rawI) / np.max(im2double(rawI)) + outputPath = iSample + '//prob_maps' + if not os.path.exists(outputPath): + os.makedirs(outputPath) + K = np.zeros((2, rawI.shape[0], rawI.shape[1])) + contours = UNet2D.singleImageInference(I, 'accumulate', 1) + hsize = int((float(I.shape[0]) * float(1 / dsFactor))) + vsize = int((float(I.shape[1]) * float(1 / dsFactor))) + contours = resize(contours, (rawI.shape[0], rawI.shape[1])) + K[1, :, :] = rawI + K[0, :, :] = contours + tifwrite(np.uint8(255 * K), + outputPath + '//' + fileNamePrefix[0] + '_ContoursPM_' + str(dapiChannel + 1) + '.tif') + del K + K = np.zeros((1, rawI.shape[0], rawI.shape[1])) + nuclei = UNet2D.singleImageInference(I, 'accumulate', 2) + nuclei = resize(nuclei, (rawI.shape[0], rawI.shape[1])) + K[0, :, :] = nuclei + tifwrite(np.uint8(255 * K), + outputPath + '//' + fileNamePrefix[0] + '_NucleiPM_' + str(dapiChannel + 1) + '.tif') + del K + UNet2D.singleImageInferenceCleanup()
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/macros.xml Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,28 @@ +<?xml version="1.0"?> +<macros> + <xml name="requirements"> + <requirements> + <requirement type="package" version="3.7">python</requirement> + <requirement type="package" version="1.15.0">tensorflow</requirement> + <requirement type="package" version="1.15.1">tensorflow-estimator</requirement> + <requirement type="package">cudnn</requirement> + <requirement type="package" version="10.0">cudatoolkit</requirement> + <requirement type="package" version="0.17.2">scikit-image</requirement> + <requirement type="package" version="1.4.1">scipy</requirement> + <requirement type="package" version="2020.7.24">tifffile</requirement> + <requirement type="package" version="2019.7.2">czifile</requirement> + <requirement type="package" version="3.2.3">nd2reader</requirement> + </requirements> + </xml> + + <xml name="version_cmd"> + <version_command>echo @VERSION@</version_command> + </xml> + <xml name="citations"> + <citations> + </citations> + </xml> + + <token name="@VERSION@">3.1.1</token> + <token name="@CMD_BEGIN@">python ${__tool_directory__}/UnMicst.py</token> +</macros>
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/models/CytoplasmIncell/checkpoint Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,2 @@ +model_checkpoint_path: "D:\\Dan\\CytoplasmIncell\\model.ckpt" +all_model_checkpoint_paths: "D:\\Dan\\CytoplasmIncell\\model.ckpt"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/models/CytoplasmIncell2/datasetMean.data Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,1 @@ +€G?±ë… ¸Qì. \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/models/CytoplasmIncell2/datasetStDev.data Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,1 @@ +€G?±ë… ¸Qì. \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/models/CytoplasmZeissNikon/checkpoint Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,2 @@ +model_checkpoint_path: "D:\\Dan\\CytoplasmZeissNikon\\model.ckpt" +all_model_checkpoint_paths: "D:\\Dan\\CytoplasmZeissNikon\\model.ckpt"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/models/mousenucleiDAPI/checkpoint Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,2 @@ +model_checkpoint_path: "D:\\Olesja\\UNet\\nuclei20x2bin1chan 3layers ks3 bs16 20input\\model.ckpt" +all_model_checkpoint_paths: "D:\\Olesja\\UNet\\nuclei20x2bin1chan 3layers ks3 bs16 20input\\model.ckpt"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/models/nucleiDAPI/checkpoint Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,2 @@ +model_checkpoint_path: "D:\\LSP\\UNet\\tonsil20x1bin1chan\\TFModel\\model.ckpt" +all_model_checkpoint_paths: "D:\\LSP\\UNet\\tonsil20x1bin1chan\\TFModel\\model.ckpt"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/models/nucleiDAPI1-5/checkpoint Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,2 @@ +model_checkpoint_path: "D:\\LSP\\UNet\\TuuliaLPTBdapiTFv2\\model.ckpt" +all_model_checkpoint_paths: "D:\\LSP\\UNet\\TuuliaLPTBdapiTFv2\\model.ckpt"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/models/nucleiDAPI1-5/datasetMean.data Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,1 @@ +€G?ÕÂ\(õÃ. \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/models/nucleiDAPILAMIN/checkpoint Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,2 @@ +model_checkpoint_path: "/home/cy101/files/CellBiology/IDAC/Clarence/LSP/UNet models/LPTCdapilamin5-36/model.ckpt" +all_model_checkpoint_paths: "/home/cy101/files/CellBiology/IDAC/Clarence/LSP/UNet models/LPTCdapilamin5-36/model.ckpt"
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/models/nucleiDAPILAMIN/datasetMean.data Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,3 @@ +€G?Ç +=p£× +. \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/models/nucleiDAPILAMIN/datasetStDev.data Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,1 @@ +€G?ÅÂ\(õÃ. \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolbox/GPUselect.py Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,19 @@ +import subprocess, re +import numpy as np + +def pick_gpu_lowest_memory(): + output = subprocess.Popen("nvidia-smi", stdout=subprocess.PIPE, shell=True).communicate()[0] + output=output.decode("ascii") + gpu_output = output[output.find("Memory-Usage"):] + # lines of the form + # | 0 8734 C python 11705MiB | + memory_regex = re.compile(r"[|]\s+?\D+?.+[ ](?P<gpu_memory>\d+)MiB /") + rows = gpu_output.split("\n") + result=[] + for row in gpu_output.split("\n"): + m = memory_regex.search(row) + if not m: + continue + gpu_memory = int(m.group("gpu_memory")) + result.append(gpu_memory) + return np.argsort(result)[0] \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolbox/PartitionOfImage.py Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,305 @@ +import numpy as np +from toolbox.imtools import * +# from toolbox.ftools import * +# import sys + +class PI2D: + Image = None + PaddedImage = None + PatchSize = 128 + Margin = 14 + SubPatchSize = 100 + PC = None # patch coordinates + NumPatches = 0 + Output = None + Count = None + NR = None + NC = None + NRPI = None + NCPI = None + Mode = None + W = None + + def setup(image,patchSize,margin,mode): + PI2D.Image = image + PI2D.PatchSize = patchSize + PI2D.Margin = margin + subPatchSize = patchSize-2*margin + PI2D.SubPatchSize = subPatchSize + + W = np.ones((patchSize,patchSize)) + W[[0,-1],:] = 0 + W[:,[0,-1]] = 0 + for i in range(1,2*margin): + v = i/(2*margin) + W[i,i:-i] = v + W[-i-1,i:-i] = v + W[i:-i,i] = v + W[i:-i,-i-1] = v + PI2D.W = W + + if len(image.shape) == 2: + nr,nc = image.shape + elif len(image.shape) == 3: # multi-channel image + nz,nr,nc = image.shape + + PI2D.NR = nr + PI2D.NC = nc + + npr = int(np.ceil(nr/subPatchSize)) # number of patch rows + npc = int(np.ceil(nc/subPatchSize)) # number of patch cols + + nrpi = npr*subPatchSize+2*margin # number of rows in padded image + ncpi = npc*subPatchSize+2*margin # number of cols in padded image + + PI2D.NRPI = nrpi + PI2D.NCPI = ncpi + + if len(image.shape) == 2: + PI2D.PaddedImage = np.zeros((nrpi,ncpi)) + PI2D.PaddedImage[margin:margin+nr,margin:margin+nc] = image + elif len(image.shape) == 3: + PI2D.PaddedImage = np.zeros((nz,nrpi,ncpi)) + PI2D.PaddedImage[:,margin:margin+nr,margin:margin+nc] = image + + PI2D.PC = [] # patch coordinates [r0,r1,c0,c1] + for i in range(npr): + r0 = i*subPatchSize + r1 = r0+patchSize + for j in range(npc): + c0 = j*subPatchSize + c1 = c0+patchSize + PI2D.PC.append([r0,r1,c0,c1]) + + PI2D.NumPatches = len(PI2D.PC) + PI2D.Mode = mode # 'replace' or 'accumulate' + + def getPatch(i): + r0,r1,c0,c1 = PI2D.PC[i] + if len(PI2D.PaddedImage.shape) == 2: + return PI2D.PaddedImage[r0:r1,c0:c1] + if len(PI2D.PaddedImage.shape) == 3: + return PI2D.PaddedImage[:,r0:r1,c0:c1] + + def createOutput(nChannels): + if nChannels == 1: + PI2D.Output = np.zeros((PI2D.NRPI,PI2D.NCPI),np.float16) + else: + PI2D.Output = np.zeros((nChannels,PI2D.NRPI,PI2D.NCPI),np.float16) + if PI2D.Mode == 'accumulate': + PI2D.Count = np.zeros((PI2D.NRPI,PI2D.NCPI),np.float16) + + def patchOutput(i,P): + r0,r1,c0,c1 = PI2D.PC[i] + if PI2D.Mode == 'accumulate': + PI2D.Count[r0:r1,c0:c1] += PI2D.W + if len(P.shape) == 2: + if PI2D.Mode == 'accumulate': + PI2D.Output[r0:r1,c0:c1] += np.multiply(P,PI2D.W) + elif PI2D.Mode == 'replace': + PI2D.Output[r0:r1,c0:c1] = P + elif len(P.shape) == 3: + if PI2D.Mode == 'accumulate': + for i in range(P.shape[0]): + PI2D.Output[i,r0:r1,c0:c1] += np.multiply(P[i,:,:],PI2D.W) + elif PI2D.Mode == 'replace': + PI2D.Output[:,r0:r1,c0:c1] = P + + def getValidOutput(): + margin = PI2D.Margin + nr, nc = PI2D.NR, PI2D.NC + if PI2D.Mode == 'accumulate': + C = PI2D.Count[margin:margin+nr,margin:margin+nc] + if len(PI2D.Output.shape) == 2: + if PI2D.Mode == 'accumulate': + return np.divide(PI2D.Output[margin:margin+nr,margin:margin+nc],C) + if PI2D.Mode == 'replace': + return PI2D.Output[margin:margin+nr,margin:margin+nc] + if len(PI2D.Output.shape) == 3: + if PI2D.Mode == 'accumulate': + for i in range(PI2D.Output.shape[0]): + PI2D.Output[i,margin:margin+nr,margin:margin+nc] = np.divide(PI2D.Output[i,margin:margin+nr,margin:margin+nc],C) + return PI2D.Output[:,margin:margin+nr,margin:margin+nc] + + + def demo(): + I = np.random.rand(128,128) + # PI2D.setup(I,128,14) + PI2D.setup(I,64,4,'replace') + + nChannels = 2 + PI2D.createOutput(nChannels) + + for i in range(PI2D.NumPatches): + P = PI2D.getPatch(i) + Q = np.zeros((nChannels,P.shape[0],P.shape[1])) + for j in range(nChannels): + Q[j,:,:] = P + PI2D.patchOutput(i,Q) + + J = PI2D.getValidOutput() + J = J[0,:,:] + + D = np.abs(I-J) + print(np.max(D)) + + K = cat(1,cat(1,I,J),D) + imshow(K) + + +class PI3D: + Image = None + PaddedImage = None + PatchSize = 128 + Margin = 14 + SubPatchSize = 100 + PC = None # patch coordinates + NumPatches = 0 + Output = None + Count = None + NR = None # rows + NC = None # cols + NZ = None # planes + NRPI = None + NCPI = None + NZPI = None + Mode = None + W = None + + def setup(image,patchSize,margin,mode): + PI3D.Image = image + PI3D.PatchSize = patchSize + PI3D.Margin = margin + subPatchSize = patchSize-2*margin + PI3D.SubPatchSize = subPatchSize + + W = np.ones((patchSize,patchSize,patchSize)) + W[[0,-1],:,:] = 0 + W[:,[0,-1],:] = 0 + W[:,:,[0,-1]] = 0 + for i in range(1,2*margin): + v = i/(2*margin) + W[[i,-i-1],i:-i,i:-i] = v + W[i:-i,[i,-i-1],i:-i] = v + W[i:-i,i:-i,[i,-i-1]] = v + + PI3D.W = W + + if len(image.shape) == 3: + nz,nr,nc = image.shape + elif len(image.shape) == 4: # multi-channel image + nz,nw,nr,nc = image.shape + + PI3D.NR = nr + PI3D.NC = nc + PI3D.NZ = nz + + npr = int(np.ceil(nr/subPatchSize)) # number of patch rows + npc = int(np.ceil(nc/subPatchSize)) # number of patch cols + npz = int(np.ceil(nz/subPatchSize)) # number of patch planes + + nrpi = npr*subPatchSize+2*margin # number of rows in padded image + ncpi = npc*subPatchSize+2*margin # number of cols in padded image + nzpi = npz*subPatchSize+2*margin # number of plns in padded image + + PI3D.NRPI = nrpi + PI3D.NCPI = ncpi + PI3D.NZPI = nzpi + + if len(image.shape) == 3: + PI3D.PaddedImage = np.zeros((nzpi,nrpi,ncpi)) + PI3D.PaddedImage[margin:margin+nz,margin:margin+nr,margin:margin+nc] = image + elif len(image.shape) == 4: + PI3D.PaddedImage = np.zeros((nzpi,nw,nrpi,ncpi)) + PI3D.PaddedImage[margin:margin+nz,:,margin:margin+nr,margin:margin+nc] = image + + PI3D.PC = [] # patch coordinates [z0,z1,r0,r1,c0,c1] + for iZ in range(npz): + z0 = iZ*subPatchSize + z1 = z0+patchSize + for i in range(npr): + r0 = i*subPatchSize + r1 = r0+patchSize + for j in range(npc): + c0 = j*subPatchSize + c1 = c0+patchSize + PI3D.PC.append([z0,z1,r0,r1,c0,c1]) + + PI3D.NumPatches = len(PI3D.PC) + PI3D.Mode = mode # 'replace' or 'accumulate' + + def getPatch(i): + z0,z1,r0,r1,c0,c1 = PI3D.PC[i] + if len(PI3D.PaddedImage.shape) == 3: + return PI3D.PaddedImage[z0:z1,r0:r1,c0:c1] + if len(PI3D.PaddedImage.shape) == 4: + return PI3D.PaddedImage[z0:z1,:,r0:r1,c0:c1] + + def createOutput(nChannels): + if nChannels == 1: + PI3D.Output = np.zeros((PI3D.NZPI,PI3D.NRPI,PI3D.NCPI)) + else: + PI3D.Output = np.zeros((PI3D.NZPI,nChannels,PI3D.NRPI,PI3D.NCPI)) + if PI3D.Mode == 'accumulate': + PI3D.Count = np.zeros((PI3D.NZPI,PI3D.NRPI,PI3D.NCPI)) + + def patchOutput(i,P): + z0,z1,r0,r1,c0,c1 = PI3D.PC[i] + if PI3D.Mode == 'accumulate': + PI3D.Count[z0:z1,r0:r1,c0:c1] += PI3D.W + if len(P.shape) == 3: + if PI3D.Mode == 'accumulate': + PI3D.Output[z0:z1,r0:r1,c0:c1] += np.multiply(P,PI3D.W) + elif PI3D.Mode == 'replace': + PI3D.Output[z0:z1,r0:r1,c0:c1] = P + elif len(P.shape) == 4: + if PI3D.Mode == 'accumulate': + for i in range(P.shape[1]): + PI3D.Output[z0:z1,i,r0:r1,c0:c1] += np.multiply(P[:,i,:,:],PI3D.W) + elif PI3D.Mode == 'replace': + PI3D.Output[z0:z1,:,r0:r1,c0:c1] = P + + def getValidOutput(): + margin = PI3D.Margin + nz, nr, nc = PI3D.NZ, PI3D.NR, PI3D.NC + if PI3D.Mode == 'accumulate': + C = PI3D.Count[margin:margin+nz,margin:margin+nr,margin:margin+nc] + if len(PI3D.Output.shape) == 3: + if PI3D.Mode == 'accumulate': + return np.divide(PI3D.Output[margin:margin+nz,margin:margin+nr,margin:margin+nc],C) + if PI3D.Mode == 'replace': + return PI3D.Output[margin:margin+nz,margin:margin+nr,margin:margin+nc] + if len(PI3D.Output.shape) == 4: + if PI3D.Mode == 'accumulate': + for i in range(PI3D.Output.shape[1]): + PI3D.Output[margin:margin+nz,i,margin:margin+nr,margin:margin+nc] = np.divide(PI3D.Output[margin:margin+nz,i,margin:margin+nr,margin:margin+nc],C) + return PI3D.Output[margin:margin+nz,:,margin:margin+nr,margin:margin+nc] + + + def demo(): + I = np.random.rand(128,128,128) + PI3D.setup(I,64,4,'accumulate') + + nChannels = 2 + PI3D.createOutput(nChannels) + + for i in range(PI3D.NumPatches): + P = PI3D.getPatch(i) + Q = np.zeros((P.shape[0],nChannels,P.shape[1],P.shape[2])) + for j in range(nChannels): + Q[:,j,:,:] = P + PI3D.patchOutput(i,Q) + + J = PI3D.getValidOutput() + J = J[:,0,:,:] + + D = np.abs(I-J) + print(np.max(D)) + + pI = I[64,:,:] + pJ = J[64,:,:] + pD = D[64,:,:] + + K = cat(1,cat(1,pI,pJ),pD) + imshow(K) +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolbox/ftools.py Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,55 @@ +from os.path import * +from os import listdir, makedirs, remove +import pickle +import shutil + +def fileparts(path): # path = file path + [p,f] = split(path) + [n,e] = splitext(f) + return [p,n,e] + +def listfiles(path,token): # path = folder path + l = [] + for f in listdir(path): + fullPath = join(path,f) + if isfile(fullPath) and token in f: + l.append(fullPath) + l.sort() + return l + +def listsubdirs(path): # path = folder path + l = [] + for f in listdir(path): + fullPath = join(path,f) + if isdir(fullPath): + l.append(fullPath) + l.sort() + return l + +def pathjoin(p,ne): # '/path/to/folder', 'name.extension' (or a subfolder) + return join(p,ne) + +def saveData(data,path): + print('saving data') + dataFile = open(path, 'wb') + pickle.dump(data, dataFile) + +def loadData(path): + print('loading data') + dataFile = open(path, 'rb') + return pickle.load(dataFile) + +def createFolderIfNonExistent(path): + if not exists(path): # from os.path + makedirs(path) + +def moveFile(fullPathSource,folderPathDestination): + [p,n,e] = fileparts(fullPathSource) + shutil.move(fullPathSource,pathjoin(folderPathDestination,n+e)) + +def copyFile(fullPathSource,folderPathDestination): + [p,n,e] = fileparts(fullPathSource) + shutil.copy(fullPathSource,pathjoin(folderPathDestination,n+e)) + +def removeFile(path): + remove(path) \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/toolbox/imtools.py Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,312 @@ +import matplotlib.pyplot as plt +import tifffile +import os +import numpy as np +from skimage import io as skio +from scipy.ndimage import * +from skimage.morphology import * +from skimage.transform import resize + +def tifread(path): + return tifffile.imread(path) + +def tifwrite(I,path): + tifffile.imsave(path, I) + +def imshow(I,**kwargs): + if not kwargs: + plt.imshow(I,cmap='gray') + else: + plt.imshow(I,**kwargs) + + plt.axis('off') + plt.show() + +def imshowlist(L,**kwargs): + n = len(L) + for i in range(n): + plt.subplot(1, n, i+1) + if not kwargs: + plt.imshow(L[i],cmap='gray') + else: + plt.imshow(L[i],**kwargs) + plt.axis('off') + plt.show() + +def imread(path): + return skio.imread(path) + +def imwrite(I,path): + skio.imsave(path,I) + +def im2double(I): + if I.dtype == 'uint16': + return I.astype('float64')/65535 + elif I.dtype == 'uint8': + return I.astype('float64')/255 + elif I.dtype == 'float32': + return I.astype('float64') + elif I.dtype == 'float64': + return I + else: + print('returned original image type: ', I.dtype) + return I + +def size(I): + return list(I.shape) + +def imresizeDouble(I,sizeOut): # input and output are double + return resize(I,(sizeOut[0],sizeOut[1]),mode='reflect') + +def imresize3Double(I,sizeOut): # input and output are double + return resize(I,(sizeOut[0],sizeOut[1],sizeOut[2]),mode='reflect') + +def imresizeUInt8(I,sizeOut): # input and output are UInt8 + return np.uint8(resize(I.astype(float),(sizeOut[0],sizeOut[1]),mode='reflect',order=0)) + +def imresize3UInt8(I,sizeOut): # input and output are UInt8 + return np.uint8(resize(I.astype(float),(sizeOut[0],sizeOut[1],sizeOut[2]),mode='reflect',order=0)) + +def normalize(I): + m = np.min(I) + M = np.max(I) + if M > m: + return (I-m)/(M-m) + else: + return I + +def snormalize(I): + m = np.mean(I) + s = np.std(I) + if s > 0: + return (I-m)/s + else: + return I + +def cat(a,I,J): + return np.concatenate((I,J),axis=a) + +def imerode(I,r): + return binary_erosion(I, disk(r)) + +def imdilate(I,r): + return binary_dilation(I, disk(r)) + +def imerode3(I,r): + return morphology.binary_erosion(I, ball(r)) + +def imdilate3(I,r): + return morphology.binary_dilation(I, ball(r)) + +def sphericalStructuralElement(imShape,fRadius): + if len(imShape) == 2: + return disk(fRadius,dtype=float) + if len(imShape) == 3: + return ball(fRadius,dtype=float) + +def medfilt(I,filterRadius): + return median_filter(I,footprint=sphericalStructuralElement(I.shape,filterRadius)) + +def maxfilt(I,filterRadius): + return maximum_filter(I,footprint=sphericalStructuralElement(I.shape,filterRadius)) + +def minfilt(I,filterRadius): + return minimum_filter(I,footprint=sphericalStructuralElement(I.shape,filterRadius)) + +def ptlfilt(I,percentile,filterRadius): + return percentile_filter(I,percentile,footprint=sphericalStructuralElement(I.shape,filterRadius)) + +def imgaussfilt(I,sigma,**kwargs): + return gaussian_filter(I,sigma,**kwargs) + +def imlogfilt(I,sigma,**kwargs): + return -gaussian_laplace(I,sigma,**kwargs) + +def imgradmag(I,sigma): + if len(I.shape) == 2: + dx = imgaussfilt(I,sigma,order=[0,1]) + dy = imgaussfilt(I,sigma,order=[1,0]) + return np.sqrt(dx**2+dy**2) + if len(I.shape) == 3: + dx = imgaussfilt(I,sigma,order=[0,0,1]) + dy = imgaussfilt(I,sigma,order=[0,1,0]) + dz = imgaussfilt(I,sigma,order=[1,0,0]) + return np.sqrt(dx**2+dy**2+dz**2) + +def localstats(I,radius,justfeatnames=False): + ptls = [10,30,50,70,90] + featNames = [] + for i in range(len(ptls)): + featNames.append('locPtl%d' % ptls[i]) + if justfeatnames == True: + return featNames + sI = size(I) + nFeats = len(ptls) + F = np.zeros((sI[0],sI[1],nFeats)) + for i in range(nFeats): + F[:,:,i] = ptlfilt(I,ptls[i],radius) + return F + +def localstats3(I,radius,justfeatnames=False): + ptls = [10,30,50,70,90] + featNames = [] + for i in range(len(ptls)): + featNames.append('locPtl%d' % ptls[i]) + if justfeatnames == True: + return featNames + sI = size(I) + nFeats = len(ptls) + F = np.zeros((sI[0],sI[1],sI[2],nFeats)) + for i in range(nFeats): + F[:,:,:,i] = ptlfilt(I,ptls[i],radius) + return F + +def imderivatives(I,sigmas,justfeatnames=False): + if type(sigmas) is not list: + sigmas = [sigmas] + derivPerSigmaFeatNames = ['d0','dx','dy','dxx','dxy','dyy','normGrad','normHessDiag'] + if justfeatnames == True: + featNames = []; + for i in range(len(sigmas)): + for j in range(len(derivPerSigmaFeatNames)): + featNames.append('derivSigma%d%s' % (sigmas[i],derivPerSigmaFeatNames[j])) + return featNames + nDerivativesPerSigma = len(derivPerSigmaFeatNames) + nDerivatives = len(sigmas)*nDerivativesPerSigma + sI = size(I) + D = np.zeros((sI[0],sI[1],nDerivatives)) + for i in range(len(sigmas)): + sigma = sigmas[i] + dx = imgaussfilt(I,sigma,order=[0,1]) + dy = imgaussfilt(I,sigma,order=[1,0]) + dxx = imgaussfilt(I,sigma,order=[0,2]) + dyy = imgaussfilt(I,sigma,order=[2,0]) + D[:,:,nDerivativesPerSigma*i ] = imgaussfilt(I,sigma) + D[:,:,nDerivativesPerSigma*i+1] = dx + D[:,:,nDerivativesPerSigma*i+2] = dy + D[:,:,nDerivativesPerSigma*i+3] = dxx + D[:,:,nDerivativesPerSigma*i+4] = imgaussfilt(I,sigma,order=[1,1]) + D[:,:,nDerivativesPerSigma*i+5] = dyy + D[:,:,nDerivativesPerSigma*i+6] = np.sqrt(dx**2+dy**2) + D[:,:,nDerivativesPerSigma*i+7] = np.sqrt(dxx**2+dyy**2) + return D + # derivatives are indexed by the last dimension, which is good for ML features but not for visualization, + # in which case the expected dimensions are [plane,channel,y(row),x(col)]; to obtain that ordering, do + # D = np.moveaxis(D,[0,3,1,2],[0,1,2,3]) + +def imderivatives3(I,sigmas,justfeatnames=False): + if type(sigmas) is not list: + sigmas = [sigmas] + + derivPerSigmaFeatNames = ['d0','dx','dy','dz','dxx','dxy','dxz','dyy','dyz','dzz','normGrad','normHessDiag'] + + # derivPerSigmaFeatNames = ['d0','normGrad','normHessDiag'] + + if justfeatnames == True: + featNames = []; + for i in range(len(sigmas)): + for j in range(len(derivPerSigmaFeatNames)): + featNames.append('derivSigma%d%s' % (sigmas[i],derivPerSigmaFeatNames[j])) + return featNames + nDerivativesPerSigma = len(derivPerSigmaFeatNames) + nDerivatives = len(sigmas)*nDerivativesPerSigma + sI = size(I) + D = np.zeros((sI[0],sI[1],sI[2],nDerivatives)) # plane, channel, y, x + for i in range(len(sigmas)): + sigma = sigmas[i] + dx = imgaussfilt(I,sigma,order=[0,0,1]) # z, y, x + dy = imgaussfilt(I,sigma,order=[0,1,0]) + dz = imgaussfilt(I,sigma,order=[1,0,0]) + dxx = imgaussfilt(I,sigma,order=[0,0,2]) + dyy = imgaussfilt(I,sigma,order=[0,2,0]) + dzz = imgaussfilt(I,sigma,order=[2,0,0]) + + D[:,:,:,nDerivativesPerSigma*i ] = imgaussfilt(I,sigma) + D[:,:,:,nDerivativesPerSigma*i+1 ] = dx + D[:,:,:,nDerivativesPerSigma*i+2 ] = dy + D[:,:,:,nDerivativesPerSigma*i+3 ] = dz + D[:,:,:,nDerivativesPerSigma*i+4 ] = dxx + D[:,:,:,nDerivativesPerSigma*i+5 ] = imgaussfilt(I,sigma,order=[0,1,1]) + D[:,:,:,nDerivativesPerSigma*i+6 ] = imgaussfilt(I,sigma,order=[1,0,1]) + D[:,:,:,nDerivativesPerSigma*i+7 ] = dyy + D[:,:,:,nDerivativesPerSigma*i+8 ] = imgaussfilt(I,sigma,order=[1,1,0]) + D[:,:,:,nDerivativesPerSigma*i+9 ] = dzz + D[:,:,:,nDerivativesPerSigma*i+10] = np.sqrt(dx**2+dy**2+dz**2) + D[:,:,:,nDerivativesPerSigma*i+11] = np.sqrt(dxx**2+dyy**2+dzz**2) + + # D[:,:,:,nDerivativesPerSigma*i ] = imgaussfilt(I,sigma) + # D[:,:,:,nDerivativesPerSigma*i+1 ] = np.sqrt(dx**2+dy**2+dz**2) + # D[:,:,:,nDerivativesPerSigma*i+2 ] = np.sqrt(dxx**2+dyy**2+dzz**2) + return D + # derivatives are indexed by the last dimension, which is good for ML features but not for visualization, + # in which case the expected dimensions are [plane,y(row),x(col)]; to obtain that ordering, do + # D = np.moveaxis(D,[2,0,1],[0,1,2]) + +def imfeatures(I=[],sigmaDeriv=1,sigmaLoG=1,locStatsRad=0,justfeatnames=False): + if type(sigmaDeriv) is not list: + sigmaDeriv = [sigmaDeriv] + if type(sigmaLoG) is not list: + sigmaLoG = [sigmaLoG] + derivFeatNames = imderivatives([],sigmaDeriv,justfeatnames=True) + nLoGFeats = len(sigmaLoG) + locStatsFeatNames = [] + if locStatsRad > 1: + locStatsFeatNames = localstats([],locStatsRad,justfeatnames=True) + nLocStatsFeats = len(locStatsFeatNames) + if justfeatnames == True: + featNames = derivFeatNames + for i in range(nLoGFeats): + featNames.append('logSigma%d' % sigmaLoG[i]) + for i in range(nLocStatsFeats): + featNames.append(locStatsFeatNames[i]) + return featNames + nDerivFeats = len(derivFeatNames) + nFeatures = nDerivFeats+nLoGFeats+nLocStatsFeats + sI = size(I) + F = np.zeros((sI[0],sI[1],nFeatures)) + F[:,:,:nDerivFeats] = imderivatives(I,sigmaDeriv) + for i in range(nLoGFeats): + F[:,:,nDerivFeats+i] = imlogfilt(I,sigmaLoG[i]) + if locStatsRad > 1: + F[:,:,nDerivFeats+nLoGFeats:] = localstats(I,locStatsRad) + return F + +def imfeatures3(I=[],sigmaDeriv=2,sigmaLoG=2,locStatsRad=0,justfeatnames=False): + if type(sigmaDeriv) is not list: + sigmaDeriv = [sigmaDeriv] + if type(sigmaLoG) is not list: + sigmaLoG = [sigmaLoG] + derivFeatNames = imderivatives3([],sigmaDeriv,justfeatnames=True) + nLoGFeats = len(sigmaLoG) + locStatsFeatNames = [] + if locStatsRad > 1: + locStatsFeatNames = localstats3([],locStatsRad,justfeatnames=True) + nLocStatsFeats = len(locStatsFeatNames) + if justfeatnames == True: + featNames = derivFeatNames + for i in range(nLoGFeats): + featNames.append('logSigma%d' % sigmaLoG[i]) + for i in range(nLocStatsFeats): + featNames.append(locStatsFeatNames[i]) + return featNames + nDerivFeats = len(derivFeatNames) + nFeatures = nDerivFeats+nLoGFeats+nLocStatsFeats + sI = size(I) + F = np.zeros((sI[0],sI[1],sI[2],nFeatures)) + F[:,:,:,:nDerivFeats] = imderivatives3(I,sigmaDeriv) + for i in range(nLoGFeats): + F[:,:,:,nDerivFeats+i] = imlogfilt(I,sigmaLoG[i]) + if locStatsRad > 1: + F[:,:,:,nDerivFeats+nLoGFeats:] = localstats3(I,locStatsRad) + return F + +def stack2list(S): + L = [] + for i in range(size(S)[2]): + L.append(S[:,:,i]) + return L + +def thrsegment(I,wsBlr,wsThr): # basic threshold segmentation + G = imgaussfilt(I,sigma=(1-wsBlr)+wsBlr*5) # min 1, max 5 + M = G > wsThr + return M \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/unmicst.xml Fri Mar 12 00:17:29 2021 +0000 @@ -0,0 +1,104 @@ +<tool id="unmicst" name="UnMicst" version="@VERSION@.1" profile="17.09"> + <description>UNet Model for Identifying Cells and Segmenting Tissue</description> + <macros> + <import>macros.xml</import> + </macros> + + <expand macro="requirements"/> + @VERSION_CMD@ + + <command detect_errors="exit_code"><![CDATA[ + #set $typeCorrected = str($image.name).replace('.ome.tiff','').replace('.ome.tif','').replace('.tiff','').replace('.tif','')+'.ome.tif' + + ln -s $image '$typeCorrected'; + + @CMD_BEGIN@ '$typeCorrected' + + #if $stackoutput + --stackOutput + #end if + + --outputPath `pwd` + --channel $channel + --model $model + --mean $mean + --std $stdev + --scalingFactor $scalingfactor; + + ## Move files to different files for from_work_dir differentiation + #if $stackoutput + mv *Probabilities*.tif Probabilities.tif; + mv *Preview*.tif Preview.tif + #else + mv *ContoursPM*.tif ContoursPM.tif; + mv *NucleiPM*.tif NucleiPM.tif + #end if + ]]></command> + + <inputs> + <param name="image" type="data" format="tiff" label="Registered TIFF"/> + <param name="model" type="select" label="Model"> + <option value="nucleiDAPI">nucleiDAPI</option> + <option value="mousenucleiDAPI">mousenucleiDAPI</option> + <option value="CytoplasmIncell">CytoplasmIncell</option> + <option value="CytoplasmZeissNikon">CytoplasmZeissNikon</option> + </param> + <param name="mean" type="float" value="-1" label="Mean (-1 for model default)"/> + <param name="stdev" type="float" value="-1" label="Standard Deviation (-1 for model default)"/> + <param name="channel" type="integer" value="0" label="Channel to perform inference on"/> + <param name="stackoutput" type="boolean" label="Stack probability map outputs"/> + <param name="scalingfactor" type="float" value="1.0" label="Factor to scale by"/> + </inputs> + + <outputs> + <data format="tiff" name="previews" from_work_dir="Preview.tif" label="${tool.name} on ${on_string}: Preview"> + <filter>stackoutput</filter> + </data> + <data format="tiff" name="probabilities" from_work_dir="Probabilities.tif" label="${tool.name} on ${on_string}: Probabilities"> + <filter>stackoutput</filter> + </data> + <data format="tiff" name="contours" from_work_dir="ContoursPM.tif" label="${tool.name} on ${on_string}: ContoursPM"> + <filter>not stackoutput</filter> + </data> + <data format="tiff" name="nuclei" from_work_dir="NucleiPM.tif" label="${tool.name} on ${on_string}: NucleiPM"> + <filter>not stackoutput</filter> + </data> + </outputs> + <help><![CDATA[ +UnMicst - UNet Model for Identifying Cells and Segmenting Tissue +Image Preprocessing +Images can be preprocessed by inferring nuclei contours via a pretrained UNet model. The model is trained on 3 classes : background, nuclei contours and nuclei centers. The resulting probability maps can then be loaded into any modular segmentation pipeline that may use (but not limited to) a marker controlled watershed algorithm. + +The only input file is: an .ome.tif or .tif (preferably flat field corrected, minimal saturated pixels, and in focus. The model is trained on images acquired at 20x with binning 2x2 or a pixel size of 0.65 microns/px. If your settings differ, you can upsample/downsample to some extent. + +Running as a Docker container + +The docker image is distributed through Dockerhub and includes UnMicst with all of its dependencies. Parallel images with and without gpu support are available. + +docker pull labsyspharm/unmicst:latest +docker pull labsyspharm/unmicst:latest-gpu +Instatiate a container and mount the input directory containing your image. + +docker run -it --runtime=nvidia -v /path/to/data:/data labsyspharm/unmicst:latest-gpu bash +When using the CPU-only image, --runtime=nvidia can be omitted: + +docker run -it -v /path/to/data:/data labsyspharm/unmicst:latest bash +UnMicst resides in the /app directory inside the container: + +root@0ea0cdc46c8f:/# python app/UnMicst.py /data/input/my.tif --outputPath /data/results +Running in a Conda environment + +If Docker is not available on your system, you can run the tool locally by creating a Conda environment. Ensure conda is installed on your system, then clone the repo and use conda.yml to create the environment. + +git clone https://github.com/HMS-IDAC/UnMicst.git +cd UnMicst +conda env create -f conda.yml +conda activate unmicst +python UnMicst.py /path/to/input.tif --outputPath /path/to/results/directory +References: +S Saka, Y Wang, J Kishi, A Zhu, Y Zeng, W Xie, K Kirli, C Yapp, M Cicconet, BJ Beliveau, SW Lapan, S Yin, M Lin, E Boyde, PS Kaeser, G Pihan, GM Church, P Yin, Highly multiplexed in situ protein imaging with signal amplification by Immuno-SABER, Nat Biotechnology (accepted) + +OHSU Wrapper Repo: https://github.com/ohsu-comp-bio/UnMicst + ]]></help> + <expand macro="citations" /> +</tool>