Skip to content

Commit

Permalink
Merge pull request #6 from pnlbwh/wip
Browse files Browse the repository at this point in the history 
atlas outPrefix instead of outDir
  • Loading branch information
@tashrifbillah
tashrifbillah authored Jul 3, 2019
2 parents 87a408d + 36a036c commit c2e7395
Show file tree
Hide file tree
Showing 3 changed files with 150 additions and 124 deletions.
88 changes: 43 additions & 45 deletions TUTORIAL.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -173,76 +173,74 @@ space of the target T1/T2 image. The set of candidate masks (or labelmaps) are f


> ./atlas.py --help-all

Makes atlas image/labelmap pairs for a target image. Option to merge labelmaps via averaging
or AntsJointFusion.

Usage:
atlas.py [SWITCHES] [SUBCOMMAND [SWITCHES]] args...

Sub-commands:
args Specify training images and labelmaps via command line arguments.; see 'atlas.py args
--help' for more info
csv Specify training images and labelmaps via a csv file. Put the images with any header in the
first column, and labelmaps with proper headers in the consecutive columns. The headers in
the labelmap columns will be used to name the generated atlas labelmaps.; see 'atlas.py csv
--help' for more info
args Specify training images and labelmaps via command line arguments.; see 'atlas.py args --help' for more info
csv Specify training images and labelmaps via a csv file. Put the images with any header in the first column, and
labelmaps with proper headers in the consecutive columns. The headers in the labelmap columns will be used to
name the generated atlas labelmaps.; see 'atlas.py csv --help' for more info


======================================================================
======================================================================

Specify training images and labelmaps via command line arguments.

Usage:
atlas.py args [SWITCHES]
atlas.py args [SWITCHES]

Switches:
--fusion VALUE:{'avg', 'wavg', 'antsJointFusion'} Also create predicted labelmap(s) by combining the
atlas labelmaps: avg is naive mathematical average,
wavg is weighted average where weights are computed
from MI between the warped atlases and target image,
antsJointFusion is local weighted averaging; the
default is wavg
-i, --images VALUE:str list of images in quotations, e.g. "img1.nrrd
img2.nrrd"; required
-l, --labels VALUE:str list of labelmap images in quotations, e.g. "mask1.nrrd
mask2.nrrd cingr1.nrrd cingr2.nrrd"; required
-n, --nproc VALUE:str number of processes/threads to use (-1 for all
available); the default is 4
--names VALUE:str list of names for generated labelmaps, e.g. "atlasmask
atlascingr"; required
-o, --out VALUE:str output directory; required
-d Debug mode, saves intermediate labelmaps to atlas-debug-<pid> in output
directory
--fusion VALUE:{'avg', 'wavg', 'antsJointFusion'} Also create predicted labelmap(s) by combining the atlas labelmaps: avg
is naive mathematical average, wavg is weighted average where weights are
computed from MI between the warped atlases and target image,
antsJointFusion is local weighted averaging; the default is wavg
-i, --images VALUE:str list of images in quotations, e.g. "img1.nrrd img2.nrrd"; required
-l, --labels VALUE:str list of labelmap images in quotations, e.g. "mask1.nrrd mask2.nrrd
cingr1.nrrd cingr2.nrrd"; required
-n, --nproc VALUE:str number of processes/threads to use (-1 for all available); the default is
4
--names VALUE:str list of names for generated labelmaps, e.g. "atlasmask atlascingr";
required
-o, --outPrefix VALUE:str output prefix, output labelmaps are saved as outPrefix-mask.nii.gz,
outPrefix-cingr.nii.gz, ...; required
-t, --target VALUE:ExistingFile target image; required



======================================================================
======================================================================

Specify training images and labelmaps via a csv file.
Put the images with any header in the first column,
and labelmaps with proper headers in the consecutive columns.
Put the images with any header in the first column,
and labelmaps with proper headers in the consecutive columns.
The headers in the labelmap columns will be used to name the generated atlas labelmaps.

Usage:
atlas.py csv [SWITCHES] csvFile

Switches:
--fusion VALUE:{'avg', 'wavg', 'antsJointFusion'} Also create predicted labelmap(s) by combining the
atlas labelmaps: avg is naive mathematical average,
wavg is weighted average where weights are computed
from MI between the warped atlases and target image,
antsJointFusion is local weighted averaging; the
default is wavg
-n, --nproc VALUE:str number of processes/threads to use (-1 for all
available); the default is 4
-o, --outDir VALUE:str output directory; required
-d Debug mode, saves intermediate labelmaps to atlas-debug-<pid> in output
directory
--fusion VALUE:{'avg', 'wavg', 'antsJointFusion'} Also create predicted labelmap(s) by combining the atlas labelmaps: avg
is naive mathematical average, wavg is weighted average where weights are
computed from MI between the warped atlases and target image,
antsJointFusion is local weighted averaging; the default is wavg
-n, --nproc VALUE:str number of processes/threads to use (-1 for all available); the default is
4
-o, --outPrefix VALUE:str output prefix, output labelmaps are saved as outPrefix-mask.nii.gz,
outPrefix-cingr.nii.gz, ...; required
-t, --target VALUE:ExistingFile target image; required



Example usage:

./atlas.py csv -t t1Nifti -o /tmp/T1labels/ -n 8 ~/pnlpipe/soft_dir/trainingDataT1AHCC-d6e5990/trainingDataT1AHCC-hdr.csv
./atlas.py csv -t t1Nifti -o /tmp/T1-labels -n 8 ~/pnlpipe/soft_dir/trainingDataT1AHCC-d6e5990/trainingDataT1AHCC-hdr.csv

The `csvFile` used here is `trainingDataT1AHCC-d6e5990/trainingDataT1AHCC-hdr.csv` which can be generated by running
[mktrainingfiles.sh](https://github.com/pnlbwh/trainingDataT1AHCC/blob/master/mktrainingfiles.sh) .
Expand Down
182 changes: 105 additions & 77 deletions scripts/atlas.py
View file
Original file line number Diff line number Diff line change
@@ -1,17 +1,29 @@
#!/usr/bin/env python
from __future__ import print_function
from plumbum import local, cli, FG
from plumbum.cmd import ComposeMultiTransform, antsApplyTransforms, MeasureImageSimilarity
from plumbum.cmd import ComposeMultiTransform, antsApplyTransforms, MeasureImageSimilarity, \
head, cut, antsRegistration
from itertools import zip_longest
import pandas as pd
import sys, os
import multiprocessing
from math import exp
from conversion.antsUtil import antsReg
from util import logfmt, save_nifti, TemporaryDirectory, load_nifti, N_CPU, N_PROC
from util import logfmt, save_nifti, TemporaryDirectory, load_nifti, N_CPU, N_PROC, dirname, pjoin

SCRIPTDIR = os.path.dirname(os.path.realpath(__file__))


# determine ANTS_VERSION
# $ antsRegistration --version
# ANTs Version: 2.2.0.dev233-g19285
# Compiled: Sep 2 2018 23:23:33

(antsRegistration['--version'] > '/tmp/ANTS_VERSION') & FG
with open('/tmp/ANTS_VERSION') as f:
content=f.read().split('\n')
ANTS_VERSION= content[0].split()[-1]

import logging
logger = logging.getLogger()
logging.basicConfig(level=logging.INFO, format=logfmt(__file__))
Expand Down Expand Up @@ -70,8 +82,12 @@ def applyWarp(moving, warp, reference, out, interpolation='Linear'):


def computeMI(target, img, miFile):
(MeasureImageSimilarity['-d', '3',
'-m', 'MI[{},{},1,256]'.format(target, img)] > miFile) & FG

if ANTS_VERSION <= '2.1.0':
(MeasureImageSimilarity['3', '2', target, img] | head['-n', '-2'] | cut['-d ', '-f6'] > miFile)()

else:
(MeasureImageSimilarity['-d', '3', '-m', 'MI[{},{},1,256]'.format(target, img)] > miFile) & FG


def weightsFromMIExp(mis, alpha):
Expand Down Expand Up @@ -169,89 +185,94 @@ def train2target(itr):
interpolation='NearestNeighbor')


def makeAtlases(target, trainingTable, outdir, fusion, threads):
def makeAtlases(target, trainingTable, outPrefix, fusion, threads, debug):

outdir = local.path(outdir)
outdir.mkdir()
with TemporaryDirectory() as tmpdir:

L= len(trainingTable)
tmpdir = local.path(tmpdir)

multiDataFrame= pd.concat([trainingTable, pd.DataFrame({'outdir': [outdir]*L, 'target': [str(target)]*L})], axis= 1)
L= len(trainingTable)

logging.info('Create {} atlases: compute transforms from images to target and apply over images'.format(L))
multiDataFrame= pd.concat([trainingTable, pd.DataFrame({'tmpdir': [tmpdir]*L, 'target': [str(target)]*L})], axis= 1)

pool = multiprocessing.Pool(threads) # Use all available cores, otherwise specify the number you want as an argument
logging.info('Create {} atlases: compute transforms from images to target and apply over images'.format(L))

pool.map_async(train2target, multiDataFrame.iterrows())
pool = multiprocessing.Pool(threads) # Use all available cores, otherwise specify the number you want as an argument

pool.close()
pool.join()
pool.map_async(train2target, multiDataFrame.iterrows())

logging.info('Fuse warped labelmaps to compute output labelmaps')
atlasimages = outdir // 'atlas*.nii.gz'
# sorting is required for applying weight to corresponding labelmap
atlasimages.sort()
pool.close()
pool.join()

if fusion.lower() == 'wavg':
logging.info('Fuse warped labelmaps to compute output labelmaps')
atlasimages = tmpdir // 'atlas*.nii.gz'
# sorting is required for applying weight to corresponding labelmap
atlasimages.sort()

ALPHA_DEFAULT= 0.45
if fusion.lower() == 'wavg':

logging.info('Compute MI between warped images and target')
pool = multiprocessing.Pool(threads)
for img in atlasimages:
print('MI between {} and target'.format(img))
miFile= img+'.txt'
pool.apply_async(func= computeMI, args= (target, img, miFile, ))
ALPHA_DEFAULT= 0.45

pool.close()
pool.join()
logging.info('Compute MI between warped images and target')
pool = multiprocessing.Pool(threads)
for img in atlasimages:
print('MI between {} and target'.format(img))
miFile= img+'.txt'
pool.apply_async(func= computeMI, args= (target, img, miFile, ))

mis= []
with open(outdir+'/MI.txt','w') as fw:
pool.close()
pool.join()

for img in atlasimages:
with open(img+'.txt') as f:
mi= f.read().strip()
fw.write(img+','+mi+'\n')
mis.append(float(mi))

weights = weightsFromMIExp(mis, ALPHA_DEFAULT)


pool = multiprocessing.Pool(threads) # Use all available cores, otherwise specify the number you want as an argument
for labelname in list(trainingTable)[1:]: #list(d) gets column names

out = outdir / labelname + '.nii.gz'
if os.path.exists(out):
os.remove(out)
labelmaps = outdir // (labelname + '*')
labelmaps.sort()

if fusion.lower() == 'avg':
print(' ')
# parellelize
# fuseAvg(labelmaps, out)
pool.apply_async(func= fuseAvg, args= (labelmaps, out, ))

elif fusion.lower() == 'antsjointfusion':
print(' ')
# atlasimages are the warped images
# labelmaps are the warped labels
# parellelize
# fuseAntsJointFusion(target, atlasimages, labelmaps, out)
pool.apply_async(func= fuseAntsJointFusion, args= (target, atlasimages, labelmaps, out, ))

elif fusion.lower() == 'wavg':
print(' ')
# parellelize
# fuseWeightedAvg(labelmaps, weights, out)
pool.apply_async(func= fuseWeightedAvg, args= (labelmaps, weights, out, ))
mis= []
with open(tmpdir+'/MI.txt','w') as fw:

else:
print('Unrecognized fusion option: {}. Skipping.'.format(fusion))
for img in atlasimages:
with open(img+'.txt') as f:
mi= f.read().strip()
fw.write(img+','+mi+'\n')
mis.append(float(mi))

weights = weightsFromMIExp(mis, ALPHA_DEFAULT)


pool = multiprocessing.Pool(threads) # Use all available cores, otherwise specify the number you want as an argument
for labelname in list(trainingTable)[1:]: #list(d) gets column names

out = outPrefix+ f'-{labelname}.nii.gz'
if os.path.exists(out):
os.remove(out)
labelmaps = tmpdir // (labelname + '*')
labelmaps.sort()

if fusion.lower() == 'avg':
print(' ')
# parellelize
# fuseAvg(labelmaps, out)
pool.apply_async(func= fuseAvg, args= (labelmaps, out, ))

elif fusion.lower() == 'antsjointfusion':
print(' ')
# atlasimages are the warped images
# labelmaps are the warped labels
# parellelize
# fuseAntsJointFusion(target, atlasimages, labelmaps, out)
pool.apply_async(func= fuseAntsJointFusion, args= (target, atlasimages, labelmaps, out, ))

elif fusion.lower() == 'wavg':
print(' ')
# parellelize
# fuseWeightedAvg(labelmaps, weights, out)
pool.apply_async(func= fuseWeightedAvg, args= (labelmaps, weights, out, ))

else:
print('Unrecognized fusion option: {}. Skipping.'.format(fusion))

pool.close()
pool.join()

if debug:
tmpdir.copy(pjoin(dirname(outPrefix), 'atlas-debug-' + str(os.getpid())))

pool.close()
pool.join()

class Atlas(cli.Application):
"""Makes atlas image/labelmap pairs for a target image. Option to merge labelmaps via averaging
Expand Down Expand Up @@ -282,7 +303,9 @@ class AtlasArgs(cli.Application):
'avg is naive mathematical average, wavg is weighted average where weights are computed from MI '
'between the warped atlases and target image, antsJointFusion is local weighted averaging', default='wavg')
out = cli.SwitchAttr(
['-o', '--out'], help='output directory', mandatory=True)
['-o', '--outPrefix'],
help='output prefix, output labelmaps are saved as outPrefix-mask.nii.gz, outPrefix-cingr.nii.gz, ...',
mandatory=True)

images = cli.SwitchAttr(
['-i', '--images'],
Expand All @@ -299,6 +322,7 @@ class AtlasArgs(cli.Application):
threads= cli.SwitchAttr(['-n', '--nproc'],
help='number of processes/threads to use (-1 for all available)',
default= N_PROC)
debug = cli.Flag('-d', help='Debug mode, saves intermediate labelmaps to atlas-debug-<pid> in output directory')

def main(self):
images = self.images.split()
Expand Down Expand Up @@ -328,8 +352,8 @@ def main(self):
if self.threads==-1 or self.threads>N_CPU:
self.threads= N_CPU

makeAtlases(self.target, trainingTable, self.out, self.fusions, self.threads)
logging.info('Made ' + self.out)
makeAtlases(self.target, trainingTable, self.out, self.fusions, self.threads, self.debug)
logging.info('Made ' + self.out + '-*.nrrd')


@Atlas.subcommand("csv")
Expand All @@ -352,16 +376,20 @@ class AtlasCsv(cli.Application):
'avg is naive mathematical average, wavg is weighted average where weights are computed from MI '
'between the warped atlases and target image, antsJointFusion is local weighted averaging', default='wavg')
out = cli.SwitchAttr(
['-o', '--outDir'], help='output directory', mandatory=True)
['-o', '--outPrefix'],
help='output prefix, output labelmaps are saved as outPrefix-mask.nii.gz, outPrefix-cingr.nii.gz, ...',
mandatory=True)
threads= cli.SwitchAttr(['-n', '--nproc'],
help='number of processes/threads to use (-1 for all available)',
default= N_PROC)
debug = cli.Flag('-d', help='Debug mode, saves intermediate labelmaps to atlas-debug-<pid> in output directory')


@cli.positional(cli.ExistingFile)
def main(self, csvFile):
trainingTable = pd.read_csv(csvFile)
makeAtlases(self.target, trainingTable, self.out, self.fusions, int(self.threads))
logging.info('Made ' + self.out)
makeAtlases(self.target, trainingTable, self.out, self.fusions, int(self.threads), self.debug)
logging.info('Made ' + self.out + '-*.nrrd')


if __name__ == '__main__':
Expand Down
Loading

0 comments on commit c2e7395

Please sign in to comment.