diff --git a/ndreg/ndreg.py b/ndreg/ndreg.py index 3b180f6..fed854a 100755 --- a/ndreg/ndreg.py +++ b/ndreg/ndreg.py @@ -28,7 +28,7 @@ def register_brain(atlas, img, modality, outdir=None): """ if outdir is None: outdir = './' - final_transform = register_affine(sitk.Normalize(atlas), + final_transform = register_affine(sitk.Normalize(atlas), img, learning_rate=1e-1, grad_tol=4e-6, @@ -51,9 +51,9 @@ def register_brain(atlas, img, modality, outdir=None): # then run lddmm e = 5e-3 s = 0.1 - atlas_lddmm, field, inv_field = register_lddmm(sitk.Normalize(atlas_affine_w), + atlas_lddmm, field, inv_field = register_lddmm(sitk.Normalize(atlas_affine_w), sitk.Normalize(img_w), - alpha_list=[0.05], + alpha_list=[0.05], scale_list = [0.0625, 0.125, 0.25, 0.5, 1.0], epsilon_list=e, sigma=s, min_epsilon_list=e*1e-6, @@ -84,7 +84,7 @@ def register_affine(atlas, img, learning_rate=1e-2, iters=200, min_step=1e-10, s numberOfIterations=iters) registration_method.SetOptimizerScalesFromPhysicalShift() - # Setup for the multi-resolution framework. + # Setup for the multi-resolution framework. registration_method.SetShrinkFactorsPerLevel(shrinkFactors=shrink_factors) registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=sigmas) registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn() @@ -101,14 +101,14 @@ def register_affine(atlas, img, learning_rate=1e-2, iters=200, min_step=1e-10, s if verbose: registration_method.AddCommand(sitk.sitkStartEvent, util.start_plot) registration_method.AddCommand(sitk.sitkEndEvent, util.end_plot) - registration_method.AddCommand(sitk.sitkMultiResolutionIterationEvent, util.update_multires_iterations) + registration_method.AddCommand(sitk.sitkMultiResolutionIterationEvent, util.update_multires_iterations) registration_method.AddCommand(sitk.sitkIterationEvent, lambda: util.plot_values(registration_method)) final_transform = registration_method.Execute(sitk.Cast(img, sitk.sitkFloat32), sitk.Cast(atlas, sitk.sitkFloat32)) return final_transform -def register_lddmm(affine_img, target_img, alpha_list=0.05, scale_list=[0.0625, 0.125, 0.25, 0.5, 1.0], +def register_lddmm(affine_img, target_img, alpha_list=0.05, scale_list=[0.0625, 0.125, 0.25, 0.5, 1.0], epsilon_list=1e-4, min_epsilon_list=1e-10, sigma=0.1, use_mi=False, iterations=200, inMask=None, refMask=None, verbose=True, out_dir=''): if sigma == None: @@ -123,12 +123,12 @@ def register_lddmm(affine_img, target_img, alpha_list=0.05, scale_list=[0.0625, useMI=use_mi, inMask=inMask, refMask=refMask, - iterations=iterations, + iterations=iterations, verbose=verbose, outDirPath=out_dir) - source_lddmm = registerer.imgApplyField(affine_img, field, - size=target_img.GetSize(), + source_lddmm = registerer.imgApplyField(affine_img, field, + size=target_img.GetSize(), spacing=target_img.GetSpacing()) return source_lddmm, field, invField @@ -155,7 +155,7 @@ def register_rigid(atlas, img, learning_rate=1e-2, iters=200, min_step=1e-10, sh numberOfIterations=iters) registration_method.SetOptimizerScalesFromPhysicalShift() - # Setup for the multi-resolution framework. + # Setup for the multi-resolution framework. registration_method.SetShrinkFactorsPerLevel(shrinkFactors=shrink_factors) registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=sigmas) registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn() @@ -172,9 +172,488 @@ def register_rigid(atlas, img, learning_rate=1e-2, iters=200, min_step=1e-10, sh if verbose: registration_method.AddCommand(sitk.sitkStartEvent, util.start_plot) registration_method.AddCommand(sitk.sitkEndEvent, util.end_plot) - registration_method.AddCommand(sitk.sitkMultiResolutionIterationEvent, util.update_multires_iterations) + registration_method.AddCommand(sitk.sitkMultiResolutionIterationEvent, util.update_multires_iterations) registration_method.AddCommand(sitk.sitkIterationEvent, lambda: util.plot_values(registration_method)) final_transform = registration_method.Execute(sitk.Cast(img, sitk.sitkFloat32), sitk.Cast(atlas, sitk.sitkFloat32)) return final_transform + + +def register_rigid_n_way(image_list, + epsilon = 1e-2, + medians = 10, + learning_rate=1e-2, + iters=200, + min_step=1e-10, + shrink_factors=[1], + sigmas=[.150], + use_mi=False, + grad_tol=1e-6, + verbose=False): + + """ + Register N 3-D images with rigid transformation. + https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3389460/ + + 1. Compute median of input images + 2. Register and resample images with median as atlas + 3. Repeat until transforms converge + + + Parameters: + n-way: + image_list, - input images (array of SimpleITK.SimpleITK.Image s) + epsilon = 1e-2 - convergence threshold: MSE between consecutive median atlases + medians = 10, - maximum number of medians to calculate (upper limit for running time) + + register_rigid: + learning_rate=1e-2, + iters=200, + min_step=1e-10, + shrink_factors=[1], + sigmas=[.150], + use_mi=False, + grad_tol=1e-6, + verbose=False + + """ + + depth = min([img.GetDepth() for img in image_list]) + source_images = [img[:,:,:depth] for img in image_list] + + atlas = sitk.GetImageFromArray(np.median([sitk.GetArrayFromImage(img) for img in image_list], axis=0)) + + errors = [] + #init tranforms to identity and compose new transforms on top of it + final_transforms = [sitk.Transform(sitk.TranslationTransform(3, (0,0,0))) for i in range(len(images))] + + #repeat until convergence + for _ in range(medians): + new_images = [] + for i in range(len(source_images)): + print("Image {} of {}".format(i+1, len(source_images))) + img = source_images[i] + transform = registerer.register_rigid(atlas, img, + learning_rate=learning_rate, + iters=iters, + min_step=min_step, + shrink_factors=shrink_factors, + sigmas=sigmas, + use_mi=use_mi, + grad_tol=grad_tol, + verbose=verbose) + new_images.append(registerer.resample(atlas, transform, img)) + final_transforms[i].AddTransform(transform) + + #images for next pass or output + source_images = new_images + new_atlas = sitk.GetImageFromArray(np.median(np.array( + [sitk.GetArrayFromImage(img) for img in source_images]), axis=0)) + + total_error = np.sum(calculate_error(atlas, new_atlas)) + errors.append(total_error) + + #Convergence check (early stopping) + if total_error < epsilon: + break + return final_transforms, errors + + + +def imgShow(img, vmin=None, vmax=None, cmap=None, alpha=None, + newFig=True, flip=[0, 0, 0], numSlices=3, useNearest=False): + """ + Displays an image. Only 2D images are supported for now + """ + if newFig: + fig = plt.figure() + + if (vmin is None) or (vmax is None): + stats = sitk.StatisticsImageFilter() + stats.Execute(img) + if vmin is None: + vmin = stats.GetMinimum() + if vmax is None: + vmax = stats.GetMaximum() + + if cmap is None: + cmap = plt.cm.gray + if alpha is None: + alpha = 1.0 + + interpolation = ['bilinear', 'none'][useNearest] + + if img.GetDimension() == 2: + plt.axis('off') + ax = plt.imshow(sitk.GetArrayFromImage(img), cmap=cmap, vmin=vmin, + vmax=vmax, alpha=alpha, interpolation=interpolation) + + elif img.GetDimension() == 3: + size = img.GetSize() + for i in range(img.GetDimension()): + start = size[2 - i] / (numSlices + 1) + sliceList = np.linspace(start, size[2 - i] - start, numSlices) + sliceSize = list(size) + sliceSize[2 - i] = 0 + + for (j, slice) in enumerate(sliceList): + sliceIndex = [0] * img.GetDimension() + sliceIndex[2 - i] = int(slice) + sliceImg = sitk.Extract(img, sliceSize, sliceIndex) + sliceArray = sitk.GetArrayFromImage(sliceImg) + if flip[i]: + sliceArray = np.transpose(sliceArray) + + plt.subplot(numSlices, img.GetDimension(), + i + img.GetDimension() * j + 1) + ax = plt.imshow(sliceArray, cmap=cmap, vmin=vmin, + vmax=vmax, alpha=alpha, interpolation=interpolation) + plt.axis('off') + else: + raise Exception("Image dimension must be 2 or 3.") + + if newFig: + plt.show() + + +def imgShowResults(inImg, refImg, field, logPath=""): + numRows = 5 + numCols = 3 + defInImg = registerer.imgApplyField(inImg, field, size=refImg.GetSize()) + checker = registerer.imgChecker(defInImg, refImg) + + sliceList = [] + for i in range(inImg.GetDimension()): + step = [5] * inImg.GetDimension() + step[2 - i] = None + grid = imgGrid(inImg.GetSize(), inImg.GetSpacing(), + step=step, field=field) + + sliceList.append(imgSlices(grid, flip=[0, 1, 1])[i]) + fig = plt.figure() + imgShowResultsRow(inImg, numRows, numCols, 0, title="$I_0$") + imgShowResultsRow(defInImg, numRows, numCols, 1, + title="$I_0 \circ \phi_{10}$") + imgShowResultsRow(checker, numRows, numCols, 2, + title="$I_0$ and $I_1$\n Checker") + imgShowResultsRow(refImg, numRows, numCols, 3, title="$I_1$") + imgShowResultsRow(sliceList, numRows, numCols, 4, title="$\phi_{10}$") + fig.subplots_adjust(hspace=0.05, wspace=0) + plt.show() + + +def imgShowResultsRow(img, numRows=1, numCols=3, rowIndex=0, title=""): + if isinstance(img, list): + sliceImgList = img + else: + sliceImgList = imgSlices(img, flip=[0, 1, 1]) + + for (i, sliceImg) in enumerate(sliceImgList): + ax = plt.subplot(numRows, numCols, rowIndex * numCols + i + 1) + plt.imshow(sitk.GetArrayFromImage(sliceImg), + cmap=plt.cm.gray, aspect='auto') + ax.set_yticklabels([]) + ax.set_xticklabels([]) + if i == 0: + plt.ylabel(title, rotation=0, labelpad=30) + # plt.axis('off') + + +def imgGrid(size, spacing, step=[10, 10, 10], field=None): + """ + Creates a grid image using with specified size and spacing with distance between lines defined by step. + If step is None along a dimention no grid lines will be plotted. + For example step=[5,5,None] will result in a grid image with grid lines every 5 voxels in the x and y directions but no grid lines in the z direction. + An optinal displacement field can be applied to the grid as well. + """ + + if not(util.is_iterable(size)): + raise Exception("size must be a list.") + if not(util.is_iterable(spacing)): + raise Exception("spacing must be a list.") + if not(util.is_iterable(step)): + raise Exception("step must be a list.") + if len(size) != len(spacing): + raise Exception("len(size) != len(spacing)") + if len(size) != len(step): + raise Exception("len(size) != len(step)") + + dimension = len(size) + offset = [0] * dimension + + for i in range(dimension): + if step[i] is None: + step[i] = size[i] + 2 + offset[i] = -1 + + gridSource = sitk.GridImageSource() + gridSource.SetSpacing(spacing) + gridSource.SetGridOffset(np.array(offset) * np.array(spacing)) + gridSource.SetOrigin([0] * dimension) + gridSource.SetSize(np.array(size)) + gridSource.SetGridSpacing(np.array(step) * np.array(spacing)) + gridSource.SetScale(255) + gridSource.SetSigma(1 * np.array(spacing)) + grid = gridSource.Execute() + + if not(field is None): + grid = sitk.WrapPad(grid, [20] * dimension, [20] * dimension) + grid = registerer.imgApplyField(grid, field, size=size) + + return grid + + +def imgSlices(img, flip=[0, 0, 0], numSlices=1): + size = img.GetSize() + sliceImgList = [] + for i in range(img.GetDimension()): + start = size[2 - i] / (numSlices + 1) + sliceList = np.linspace(start, size[2 - i] - start, numSlices) + sliceSize = list(size) + sliceSize[2 - i] = 0 + + for (j, slice) in enumerate(sliceList): + sliceIndex = [0] * img.GetDimension() + sliceIndex[2 - i] = int(slice) + sliceImg = sitk.Extract(img, sliceSize, sliceIndex) + + if flip[i]: + sliceImgDirection = sliceImg.GetDirection() + sliceImg = sitk.PermuteAxesImageFilter().Execute( + sliceImg, range(sliceImg.GetDimension() - 1, -1, -1)) + sliceImg.SetDirection(sliceImgDirection) + sliceImgList.append(sliceImg) + + return sliceImgList + + +def imgPercentile(img, percentile): + if percentile < 0.0 or percentile > 1.0: + raise Exception("Percentile should be between 0.0 and 1.0") + + (values, bins) = np.histogram(sitk.GetArrayFromImage(img), bins=255) + cumValues = np.cumsum(values).astype(float) + cumValues = (cumValues - cumValues.min()) / cumValues.ptp() + + index = np.argmax(cumValues > percentile) - 1 + value = bins[index] + return value + + +def imgMetamorphosisSlicePlotterRow( + img, numRows=1, numCols=3, rowIndex=0, title="", vmin=None, vmax=None): + if isinstance(img, list): + sliceImgList = img + else: + if vmax is None or (vmin is None): + stats = sitk.StatisticsImageFilter() + stats.Execute(img) + if vmin is None: + vmin = stats.GetMinimum() + if vmax is None: + vmax = stats.GetMaximum() + sliceImgList = imgSlices(img, flip=[0, 1, 1]) + + for (i, sliceImg) in enumerate(sliceImgList): + ax = plt.subplot(numRows, numCols, rowIndex * numCols + i + 1) + plt.imshow( + sitk.GetArrayFromImage(sliceImg), + cmap=plt.cm.gray, + aspect='auto', + vmax=vmax, + vmin=vmin) + ax.set_yticks([]) + ax.set_xticks([]) + if i == 0: + plt.ylabel(title, rotation=0, labelpad=40) + + +def imgMetamorphosisSlicePlotter(inImg, refImg, field): + numRows = 5 + numCols = 3 + defInImg = registerer.imgApplyField(inImg, field, size=refImg.GetSize()) + inImg = registerer.imgApplyAffine( + inImg, [ + 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0], size=refImg.GetSize()) + checker = registerer.imgChecker(defInImg, refImg) + + sliceList = [] + for i in range(inImg.GetDimension()): + step = [20] * inImg.GetDimension() + step[2 - i] = None + grid = imgGrid( + inImg.GetSize(), + inImg.GetSpacing(), + step=step, + field=field) + + sliceList.append(imgSlices(grid, flip=[0, 1, 1])[i]) + + imgMetamorphosisSlicePlotterRow( + inImg, + numRows, + numCols, + 0, + title="$I_0$", + vmax=imgPercentile( + inImg, + 0.99)) + imgMetamorphosisSlicePlotterRow( + defInImg, + numRows, + numCols, + 1, + title="$I(1)$", + vmax=imgPercentile( + defInImg, + 0.99)) + imgMetamorphosisSlicePlotterRow( + checker, + numRows, + numCols, + 2, + title="$I(1)$ and $I_1$\n Checker", + vmax=imgPercentile( + checker, + 0.99)) + imgMetamorphosisSlicePlotterRow( + refImg, + numRows, + numCols, + 3, + title="$I_1$", + vmax=imgPercentile( + refImg, + 0.99)) + imgMetamorphosisSlicePlotterRow( + sliceList, numRows, numCols, 4, title="$\phi_{10}$") + plt.gcf().subplots_adjust(hspace=0.1, wspace=0.025) + + +def imgMetamorphosisLogPlotter( + logPathList, labelList=None, useLog=False, useTime=False): + if not(util.is_iterable(logPathList)): + raise Exception("logPathList should be a list.") + + if labelList is None: + labelList = ["Step {0}".format(i) + for i in range(1, len(logPathList) + 1)] + else: + if not(util.is_iterable(labelList)): + raise Exception("labelList should be a list.") + if len(labelList) != len(logPathList): + raise Exception( + "Number of labels should equal number of log files.") + + initialPercent = 1.0 + initialX = 0 + levelXList = [] + levelPercentList = [] + for (i, logPath) in enumerate(logPathList): + percentList = imgMetamorphosisLogParser(logPath)[:, 1] * initialPercent + numIterations = len(percentList) + if useTime: + # Parse run time from log and convert to minutes + time = float(util.txt_read(logPath).split( + "Time = ")[1].split("s ")[0]) / 60.0 + xList = np.linspace(0, time, numIterations + 1)[1:] + initialX + else: + xList = np.arange(0, numIterations) + initialX + + if not useLog: + if i == 0: + xList = np.array([initialX] + list(xList)) + percentList = np.array([initialPercent] + list(percentList)) + + levelXList += [xList] + levelPercentList += [percentList] + + initialPercent = percentList[-1] + initialX = xList[-1] + + for i in range(len(levelXList)): + if i > 0: + xList = np.concatenate((levelXList[i - 1][-1:], levelXList[i])) + percentList = np.concatenate( + (levelPercentList[i - 1][-1:], levelPercentList[i])) + + else: + xList = levelXList[i] + percentList = levelPercentList[i] + + plt.plot(xList, percentList, label=labelList[i], linewidth=1.5) + + # Add plot annotations + if useTime: + plt.xlabel("Time (Minutes)") + else: + plt.xlabel("Iteration") + + plt.ylabel("Normalized $M(I(1), I_1)$") + plt.legend(loc="center left", bbox_to_anchor=(1, 0.5)) + if useLog: + plt.xscale("log") + ax = plt.gca() + ax.xaxis.set_major_formatter(ScalarFormatter()) + plt.autoscale(enable=True, axis='x', tight=True) + + # Fix maximum y to 1.0 + ylim = list(ax.get_ylim()) + ylim[1] = 1.0 + ax.set_ylim(ylim) + + +def imgMetamorphosisLogParser(logPath): + logText = util.txt_read(logPath) + lineList = logText.split("\n") + + for (lineIndex, line) in enumerate(lineList): + if "E, E_velocity, E_rate, E_image" in line: + break + + dataArray = np.empty((0, 5), float) + for line in lineList[lineIndex:]: + if "E =" in line: + break + + try: + (iterationString, dataString) = line.split(".\t") + except BaseException: + continue + + (energyString, velocityEnergyString, rateEnergyString, + imageEnergyString, learningRateString) = (dataString.split(",")) + (energy, velocityEnergy, rateEnergy, learningRate) = map(float, [ + energyString, velocityEnergyString, rateEnergyString, learningRateString]) + (imageEnergy, imageEnergyPercent) = map( + float, imageEnergyString.replace("(", "").replace("%)", "").split()) + + imageEnergy = float(imageEnergyString.split(" (")[0]) + dataRow = np.array( + [[energy, imageEnergyPercent / 100, velocityEnergy, learningRate, rateEnergy]]) + dataArray = np.concatenate((dataArray, dataRow), axis=0) + + return dataArray + +def imgChecker(inImg, refImg, useHM=True, pattern=None): + """ + Checkerboards input image with reference image + """ + inImg = sitk.Cast(inImg, refImg.GetPixelID()) + inSize = list(inImg.GetSize()) + refSize = list(refImg.GetSize()) + if pattern is None: pattern = [4]*inImg.GetDimension() + + if(inSize != refSize): + sourceSize = np.array([inSize, refSize]).min(0) + # Empty image with same size as reference image + tmpImg = sitk.Image(refSize, refImg.GetPixelID()) + tmpImg.CopyInformation(refImg) + inImg = sitk.Paste(tmpImg, inImg, sourceSize) + + if useHM: + inImg = preprocessor.imgHM(inImg, refImg) + + return sitk.CheckerBoardImageFilter().Execute(inImg, refImg, pattern) + diff --git a/ndreg/preprocessor.py b/ndreg/preprocessor.py index 8c0ef8e..aa0b9ba 100644 --- a/ndreg/preprocessor.py +++ b/ndreg/preprocessor.py @@ -11,7 +11,7 @@ def preprocess_brain(img, spacing, modality, image_orientation, atlas_orientation='pir'): """Perform all preprocessing steps associated with a given imaging modality. - + Parameters: ---------- img : {SimpleITK.SimpleITK.Image} @@ -24,7 +24,7 @@ def preprocess_brain(img, spacing, modality, image_orientation, atlas_orientatio A 3-letter string describing the orientation of the brain along the x, y, and z axes. See (http://www.grahamwideman.com/gw/brain/orientation/orientterms.htm) for more information atlas_orientation : {str}, optional Orientation of the atlas you are using. (the default is 'pir', which is the orientation for the Allen Reference Atlas.) - + Returns ------- SimpleITK.SimpleITK.Image @@ -42,19 +42,19 @@ def preprocess_brain(img, spacing, modality, image_orientation, atlas_orientatio def create_mask(img, use_triangle=False): """Creates a mask of the image to separate brain from background using triangle or otsu thresholding. Otsu thresholding is the default. - + Parameters: ---------- img : {SimpleITK.SimpleITK.Image} Image to compute the mask on. use_triangle : {bool}, optional Set to True if you want to use triangle thresholding. (the default is False, which results in Otsu thresholding) - + Returns ------- SimpleITK.SimpleITK.Image Binary mask with 1s as the foreground and 0s as the background. - """ + """ test_mask = None if use_triangle: @@ -78,7 +78,7 @@ def create_mask(img, use_triangle=False): def correct_bias_field(img, mask=None, scale=0.25, niters=[50, 50, 50, 50]): """Correct bias field in image using the N4ITK algorithm (http://bit.ly/2oFwAun) - + Parameters: ---------- img : {SimpleITK.SimpleITK.Image} @@ -89,7 +89,7 @@ def correct_bias_field(img, mask=None, scale=0.25, niters=[50, 50, 50, 50]): Scale at which to compute the bias correction. (the default is 0.25, which results in bias correction computed on an image downsampled to 1/4 of it's original size) niters : {list}, optional Number of iterations per resolution. Each additional entry in the list adds an additional resolution at which the bias is estimated. (the default is [50, 50, 50, 50] which results in 50 iterations per resolution at 4 resolutions) - + Returns ------- SimpleITK.SimpleITK.Image @@ -103,11 +103,11 @@ def correct_bias_field(img, mask=None, scale=0.25, niters=[50, 50, 50, 50]): stats.Execute(img) std = math.sqrt(stats.GetVariance()) img_rescaled = sitk.Cast(img, sitk.sitkFloat32) + 0.1*std - + spacing = np.array(img_rescaled.GetSpacing())/scale img_ds = imgResample(img_rescaled, spacing=spacing) img_ds = sitk.Cast(img_ds, sitk.sitkFloat32) - + # Calculate bias if mask is None: @@ -119,10 +119,10 @@ def correct_bias_field(img, mask=None, scale=0.25, niters=[50, 50, 50, 50]): mask_sitk.CopyInformation(img) mask = mask_sitk mask = imgResample(mask, spacing=spacing) - + img_ds_bc = sitk.N4BiasFieldCorrection(img_ds, mask, 0.001, niters) bias_ds = img_ds_bc / sitk.Cast(img_ds, img_ds_bc.GetPixelID()) - + # Upsample bias bias = imgResample(bias_ds, spacing=img.GetSpacing(), size=img.GetSize()) @@ -132,7 +132,7 @@ def correct_bias_field(img, mask=None, scale=0.25, niters=[50, 50, 50, 50]): def remove_grid_artifact(img, z_axis=1, sigma=10, mask=None): """Remove the gridding artifact from COLM images. - + Parameters: ---------- img : {SimpleITK.SimpleITK.Image} @@ -143,7 +143,7 @@ def remove_grid_artifact(img, z_axis=1, sigma=10, mask=None): The variance of the gaussian used to blur the image. Larger sigma means more grid correction but stronger edge artifacts. (the default is 10, which empirically works well our data at 50 um) mask : {SimpleITK.SimpleITK.Image}, optional An image with 1s representing the foreground (brain) and 0s representing the background. (the default is None, which will use otsu thresholding to create the brain mask.) - + Returns ------- SimpleITK.SimpleITK.Image @@ -169,7 +169,7 @@ def remove_grid_artifact(img, z_axis=1, sigma=10, mask=None): def imgReorient(img, in_orient, out_orient): """Reorients input image to match out_orient. - + Parameters: ---------- img : {SimpleITK.SimpleITK.Image} @@ -178,7 +178,7 @@ def imgReorient(img, in_orient, out_orient): 3-letter string indicating orientation of brain. out_orient : {str} 3-letter string indicating desired orientation of input. - + Returns ------- SimpleITK.SimpleITK.Image @@ -230,7 +230,7 @@ def imgReorient(img, in_orient, out_orient): def imgResample(img, spacing, size=[], useNearest=False, origin=[0,0,0], outsideValue=0): """Resample image to certain spacing and size. - + Parameters: ---------- img : {SimpleITK.SimpleITK.Image} @@ -245,14 +245,14 @@ def imgResample(img, spacing, size=[], useNearest=False, origin=[0,0,0], outside The location in physical space representing the [0,0,0] voxel in the input image. (the default is [0,0,0]) outsideValue : {int}, optional value used to pad are outside image (the default is 0) - + Returns ------- SimpleITK.SimpleITK.Image Resampled input image. """ - + if len(spacing) != img.GetDimension(): raise Exception( "len(spacing) != " + str(img.GetDimension())) @@ -285,6 +285,26 @@ def imgResample(img, spacing, size=[], useNearest=False, origin=[0,0,0], outside identityDirection, outsideValue) +def normalize(img, percentile=0.99): + if percentile < 0.0 or percentile > 1.0: + raise Exception("Percentile should be between 0.0 and 1.0") + + #Accept ndarray images or sitk images + if type(img) is np.ndarray: + sitk_img = sitk.GetImageFromArray(img) + else: + sitk_img = img + + #just taking code from ndreg.py.... + (values, bins) = np.histogram(sitk.GetArrayFromImage(sitk_img), bins=255) + cumValues = np.cumsum(values).astype(float) + cumValues = (cumValues - cumValues.min()) / cumValues.ptp() + + index = np.argmax(cumValues > percentile) - 1 + max_val = bins[index] + + return sitk.Clamp(sitk_img, upperBound=max_val) / max_val + def downsample_and_reorient(atlas, target, atlas_orient, target_orient, spacing, size=[], set_origin=True, dv_atlas=0.0, dv_target=0.0): """ make sure img1 is the source and img2 is the target. @@ -306,7 +326,7 @@ def downsample_and_reorient(atlas, target, atlas_orient, target_orient, spacing, out_target = resampler.Execute(target_r) resampler.SetDefaultPixelValue(dv_atlas) out_atlas = resampler.Execute(atlas) - + assert(out_target.GetOrigin() == out_atlas.GetOrigin()) assert(out_target.GetSize() == out_atlas.GetSize()) assert(out_target.GetSpacing() == out_atlas.GetSpacing()) @@ -314,7 +334,7 @@ def downsample_and_reorient(atlas, target, atlas_orient, target_orient, spacing, def imgHM(img, ref_img, numMatchPoints=64, numBins=256): """Performs histogram matching on two images. - + Parameters: ---------- img : {SimpleITK.SimpleITK.Image} @@ -325,7 +345,7 @@ def imgHM(img, ref_img, numMatchPoints=64, numBins=256): number of quantile values to be matched. (the default is 64) numBins : {int}, optional Number of bins used in computation of the histogram(the default is 256) - + Returns ------- SimpleITK.SimpleITK.Image @@ -334,4 +354,4 @@ def imgHM(img, ref_img, numMatchPoints=64, numBins=256): img = sitk.Cast(img, ref_img.GetPixelID()) return sitk.HistogramMatchingImageFilter().Execute( - img, ref_img, numBins, numMatchPoints, False) \ No newline at end of file + img, ref_img, numBins, numMatchPoints, False) diff --git a/preprocess.py b/preprocess.py new file mode 100644 index 0000000..6c2b97b --- /dev/null +++ b/preprocess.py @@ -0,0 +1,127 @@ +import matplotlib +matplotlib.use('Agg') + +import matplotlib.gridspec as gridspec +import matplotlib.pyplot as plt + +import ndreg +from ndreg import preprocessor, registerer, util +import SimpleITK as sitk +import numpy as np + +import time + + + +def normalize(img, percentile=0.99): + #Accept ndarray images or sitk images + if type(img) is np.ndarray: + sitk_img = sitk.GetImageFromArray(img) + else: + sitk_img = img + max_val = ndreg.imgPercentile(sitk_img, percentile) + return sitk.Clamp(sitk_img, upperBound=max_val) / max_val + +def process(img): + """ + Input: SITK Image + Output: Clipped and normalized SITK Image + """ + temp = np.clip(sitk.GetArrayFromImage(img), -3000, 12000) + return normalize(temp) + +def overlay(img1, img2, title, save=None): + fig = plt.figure(figsize=(5, 7)) + plt.imshow(sitk.GetArrayViewFromImage(img1), cmap='Purples', alpha=0.5) + plt.imshow(sitk.GetArrayViewFromImage(img2), cmap='Greens', alpha=0.5) + plt.axis('off') + plt.title(title) + if save is not None: + plt.savefig(save) +# plt.show() + plt.close() + +print("Reading data") +tp1 = util.imgRead('./data/R04_tp1.tif') +tp2 = util.imgRead('./data/R04_tp2.tif') +tp3 = util.imgRead('./data/R04_tp3.tif') +tp4 = util.imgRead('./data/R04_tp4.tif') + +tp1_slices = [tp1[:,:,i] for i in range(50)] +tp2_slices = [tp2[:,:,i] for i in range(50)] +tp3_slices = [tp3[:,:,i] for i in range(50)] +tp4_slices = [tp4[:,:,i] for i in range(50)] + +print("Processing images") +tp1_processed = sitk.GetImageFromArray([sitk.GetArrayFromImage(process(tp1_slices[i])) for i in range(50)]) +tp2_processed = sitk.GetImageFromArray([sitk.GetArrayFromImage(process(tp2_slices[i])) for i in range(50)]) +tp3_processed = sitk.GetImageFromArray([sitk.GetArrayFromImage(process(tp3_slices[i])) for i in range(50)]) +tp4_processed = sitk.GetImageFromArray([sitk.GetArrayFromImage(process(tp4_slices[i])) for i in range(50)]) + +tp1_processed1 = [sitk.GetArrayFromImage(process(tp1_slices[i])) for i in range(50)] +tp2_processed1 = [sitk.GetArrayFromImage(process(tp2_slices[i])) for i in range(50)] +tp3_processed1 = [sitk.GetArrayFromImage(process(tp3_slices[i])) for i in range(50)] +tp4_processed1 = [sitk.GetArrayFromImage(process(tp4_slices[i])) for i in range(50)] + + +final_transform_12 = registerer.register_rigid(tp1_processed, tp2_processed, learning_rate=1e-1, iters=25) +print("Computing corrected image (of timepoint 2)") +corrected_img_12 = registerer.resample(tp1_processed, final_transform_12, tp2_processed) + +errors_12 = [] +for i in range(50): + error = registerer.imgMSE(normalize(tp1_processed[:,:,i]), normalize(corrected_img_12[:,:,i])) + errors_12.append(error) + + overlay(tp1_processed[:,:,i], corrected_img_12[:,:,i], "Timepoint 1 and Timepoint 2 Registered Overlay\n(With Processing)\nSlice {}".format(i), save='output/process/tp1tp2/tp1tp2_slice{:02d}'.format(i)) + print("Slice {}: Registration error is: {} voxels^2".format(i, error)) + +plt.title("Square Voxel Error vs Z-Slice\nTimepoint 1 vs Timepoint 2") +plt.xlabel('Log Square Voxel Error') +plt.ylabel('Z-Slice') +plt.gca().invert_yaxis() +plt.plot(np.log(errors_12), range(50)) +plt.savefig('output/voxelerror12.png') +plt.close() + + +final_transform_13 = registerer.register_rigid(tp1_processed, tp3_processed, learning_rate=1e-1, iters=25) +print("Computing corrected image (of timepoint 3)") +corrected_img_13 = registerer.resample(tp1_processed, final_transform_13, tp3_processed) + +errors_13 = [] +for i in range(50): + error = registerer.imgMSE(normalize(tp1_processed[:,:,i]), normalize(corrected_img_13[:,:,i])) + errors_13.append(error) + + overlay(tp1_processed[:,:,i], corrected_img_13[:,:,i], "Timepoint 1 and Timepoint 3 Registered Overlay\n(With Processing)\nSlice {}".format(i), save='output/process/tp1tp3/tp1tp3_slice{:02d}'.format(i)) + print("Slice {}: Registration error is: {} voxels^2".format(i, error)) + +plt.title("Square Voxel Error vs Z-Slice\nTimepoint 1 vs Timepoint 3") +plt.xlabel('Log Square Voxel Error') +plt.ylabel('Z-Slice') +plt.gca().invert_yaxis() +plt.plot(np.log(errors_13), range(50)) +plt.savefig('output/voxelerror13.png') +plt.close() + + +final_transform_14 = registerer.register_rigid(tp1_processed, tp4_processed, learning_rate=1e-1, iters=25) +print("Computing corrected image (of timepoint 4)") +corrected_img_14 = registerer.resample(tp1_processed, final_transform_14, tp4_processed) + +errors_14 = [] +for i in range(50): + error = registerer.imgMSE(normalize(tp1_processed[:,:,i]), normalize(corrected_img_14[:,:,i])) + errors_14.append(error) + + overlay(tp1_processed[:,:,i], corrected_img_14[:,:,i], "Timepoint 1 and Timepoint 2 Registered Overlay\n(With Processing)\nSlice {}".format(i), save='output/process/tp1tp4/tp1tp4_slice{:02d}'.format(i)) + print("Slice {}: Registration error is: {} voxels^2".format(i, error)) + +plt.title("Square Voxel Error vs Z-Slice\nTimepoint 1 vs Timepoint 4") +plt.xlabel('Log Square Voxel Error') +plt.ylabel('Z-Slice') +plt.gca().invert_yaxis() +plt.plot(np.log(errors_14), range(50)) +plt.savefig('output/voxelerror14.png') +plt.close()