add ndd algotithm

.gitignore

@@ -1,2 +1,3 @@
 .ipynb_checkpoints
-__pycache__
+__pycache__
+*.pyc

DoseComparison.py

@@ -3,85 +3,103 @@
 DoseComparison
     numerical evaluation of dose distributions comparison
 
-:Date: 2016-11-06
-:Version: 0.0.2
+:Date: 2017-01-02
+:Version: 1.0.0
 :Author: ongchi
 :Copyright: Copyright (c) 2016, ongchi
 :License: BSD 3-Clause License
 '''
 import numpy as np
 
 
-def gamma_index(refimg, tstimg, dta=1, dd=0.05):
-    ''' gamma evaluation of dose image distributions
+def DoseComparison(refimg, tstimg, delta_r=1, delta_d=0.05):
+    ''' gamma evaluation and normalized dose difference of dose image distributions
 
     :param refimg: reference dose image
     :param tstimg: test dose image
-    :param dta: distance-to-agreement criterion (voxels)
-    :param dd: dose difference criterion
+    :param delta_r: spatial criterion (voxels)
+    :param delta_d: dose criterion (percentage)
 
     :type refimg: numpy.ndarray
     :type tstimg: numpy.ndarray
-    :type dta: int
-    :type dd: float
-    :rtype: numpy.ndarray
+    :type delta_r: int
+    :type delta_d: float
+    :rtype: numpy.ndarray, numpy.ndarray
     '''
     # check for valid arguments
     if refimg.shape != tstimg.shape:
         raise Exception("ValueError: shape mismatch: refimg and tstimg must have the same shape")
-    if dta <= 0 or int(dta) != dta:
-        raise Exception("ValueError: dta is an integer greater than zero")
-    if dd <= 0 or dd >= 1:
-        raise Exception("ValueError: dd is a float number between 0 (exclusive) and 1 (exclusive)")
+    if delta_r <= 0 or int(delta_r) != delta_r:
+        raise Exception("ValueError: delta_r is an integer greater than zero")
+    if delta_d <= 0 or delta_d >= 1:
+        raise Exception("ValueError: delta_d is a float number between 0 (exclusive) and 1 (exclusive)")
 
     diff = tstimg - refimg
 
-    gamma = np.ma.empty_like(diff)
-    gamma[:] = np.inf
+    _ = np.empty(np.array(tstimg.shape) + delta_r * 2)
+    exTest = np.ma.array(_, mask=np.ones_like(_, dtype=bool))
+    _ = slice(delta_r, -delta_r)
+    exTest.data[_, _, _], exTest.mask[_, _, _] = tstimg, False
 
-    tmp = np.empty(np.array(tstimg.shape) + dta * 2)
-    exTest = np.ma.array(tmp, mask=np.ones_like(tmp, dtype=bool))
-    exTest.data[dta:-dta, dta:-dta, dta:-dta] = tstimg
-    exTest.mask[dta:-dta, dta:-dta, dta:-dta] = False
+    # distance map
+    distRange = np.arange(-delta_r, delta_r + 1, 1, dtype=float)
+    _ = np.array(np.meshgrid(distRange, distRange, distRange))
+    distMap = np.sqrt(np.sum(_ ** 2, axis=0))
 
-    # distance grid
-    distRange = np.arange(-dta, dta + 1, 1, dtype=float)
-    tmp = np.array(np.meshgrid(distRange, distRange, distRange))
-    distGrid = np.sqrt(np.sum(tmp ** 2, axis=0))
-
-    # mask of distance grid
+    # mask out distance map if center of vexel is in delta_r
     distRange[distRange < 0] += 0.5
     distRange[distRange > 0] -= 0.5
-    tmp = np.array(np.meshgrid(distRange, distRange, distRange))
-    distMask = np.sqrt(np.sum(tmp ** 2, axis=0)) >= dta
+    _ = np.array(np.meshgrid(distRange, distRange, distRange))
+    distMask = np.sqrt(np.sum(_ ** 2, axis=0)) >= delta_r
+
+    # mask distance within delta_r
+    dist = np.ma.array(distMap, mask=distMask)
+    dist[dist > delta_r] = delta_r
 
-    # masked distance within dta
-    dist = np.ma.array(distGrid, mask=distMask)
-    dist[dist > dta] = dta
+    # gamma
+    gamma = np.ma.empty_like(diff)
+    gamma[:] = np.inf
+    _sqDist = (dist / delta_r) ** 2
 
-    # gamma evaluation
-    nz, ny, nx = diff.shape
-    _sqDTA = (dist / dta) ** 2
+    # normalized dose difference
+    madd = np.ma.empty_like(diff)
+    madd[:] = -np.inf
+    l_min_dose = np.ma.empty_like(diff)
+    l_min_dose[:] = np.inf
+    l_min_dist = np.ma.empty_like(diff)
+
+    nx, ny, nz = diff.shape
     it = np.nditer(dist, ("multi_index", ))
     while not it.finished:
         i, j, k = idx = it.multi_index
-        _volSlice = [slice(i, i + nz), slice(j, j + ny), slice(k, k + nx)]
+        _volSlice = [slice(i, i + nx), slice(j, j + ny), slice(k, k + nz)]
 
-        # skip masked voxel
+        # skip masked voxels
         if distMask[idx] or np.alltrue(exTest[_volSlice].mask):
             it.iternext()
             continue
 
-        _sqDD = ((exTest[_volSlice] - refimg) / refimg / dd) ** 2
-        _gamma = np.sqrt(_sqDTA[idx] + _sqDD)
+        # gamma index
+        _sqDose = ((exTest[_volSlice] - refimg) / refimg / delta_d) ** 2
+        _gamma = np.sqrt(_sqDist[idx] + _sqDose)
+        _ = np.bitwise_and(gamma > _gamma, np.bitwise_not(_gamma.mask))
+        gamma[_] = _gamma[_]
 
-        # assign minimum values
-        _mask = np.bitwise_and(gamma > _gamma, np.bitwise_not(_gamma.mask))
-        gamma[_mask] = _gamma[_mask]
+        # madd
+        _ = l_min_dose > exTest[_volSlice]
+        l_min_dose[_] = exTest[_]
+        l_min_dist[_] = dist[idx]
 
         it.iternext()
 
-    return gamma
+    # ndd calculation
+    sr = np.sqrt(delta_d ** 2 - (delta_d * l_min_dist / delta_r) ** 2)
+    madd = ( (np.abs(l_min_dose - refimg) / refimg) < sr
+        ) * ( sr - np.abs(l_min_dose - refimg) + delta_d )
+    madd[madd < delta_d] = delta_d
+    ndd = diff / madd * delta_d
+
+    return gamma, ndd
 
 
 if __name__ == '__main__':
@@ -93,5 +111,4 @@ def wave(x, y, z, v):
     img1 = wave(X, Y, Z, 10)
     img2 = wave(X, Y, Z, 12)
 
-    gamma = gamma_index(img1, img2)
-    print(gamma)
+    gamma, ndd = DoseComparison(img1, img2, delta_r=2, delta_d=0.1)

README.md

@@ -1,5 +1,5 @@
 # DoseComparison
-Radiotherapy dose distributions comparison by gamma index.
+Radiotherapy dose distributions comparison by gamma index and normalized dose difference.
 
 See example.ipynb for example usage.
 
@@ -12,16 +12,23 @@ See example.ipynb for example usage.
 # Change Log
 All notable changes to this project will be documented in this file.
 
-The format is based on [Keep a Changelog](http://keepachangelog.com/) 
+The format is based on [Keep a Changelog](http://keepachangelog.com/)
 and this project adheres to [Semantic Versioning](http://semver.org/).
 
-## 0.0.1 - 2014-05-31
+## [1.0.0] - 2017-01-03
+### Added
+- ndd algorithm
+
+### Changed
+- change function name 'gamma_index' to DoseComparison
+- update example
+
+## [0.0.1] - 2016-11-07
 ### Added
 - gamma_index function
 - example.ipynb for generated wavelet distribution comparison
 
 # TODO
-- Add ndd algorithm
 - Parallel execution support
 - Add tests
 - Add some "real world" cases