Seg3D  2.4
Seg3D is a free volume segmentation and processing tool developed by the NIH Center for Integrative Biomedical Computing at the University of Utah Scientific Computing and Imaging (SCI) Institute.
threshold.hxx
1 /*
2  For more information, please see: http://software.sci.utah.edu
3 
4  The MIT License
5 
6  Copyright (c) 2016 Scientific Computing and Imaging Institute,
7  University of Utah.
8 
9 
10  Permission is hereby granted, free of charge, to any person obtaining a
11  copy of this software and associated documentation files (the "Software"),
12  to deal in the Software without restriction, including without limitation
13  the rights to use, copy, modify, merge, publish, distribute, sublicense,
14  and/or sell copies of the Software, and to permit persons to whom the
15  Software is furnished to do so, subject to the following conditions:
16 
17  The above copyright notice and this permission notice shall be included
18  in all copies or substantial portions of the Software.
19 
20  THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
21  OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
22  FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
23  THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
24  LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
25  FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
26  DEALINGS IN THE SOFTWARE.
27 */
28 
29 // File : threshold.hxx
30 // Author : Pavel A. Koshevoy
31 // Created : 2006/04/05 12:36
32 // Copyright : (C) 2004-2008 University of Utah
33 // Description : Helper functions for thresholding 2D image to isolate
34 // peaks/clusters in the data.
35 
36 #ifndef THRESHOLD_HXX_
37 #define THRESHOLD_HXX_
38 
39 // local includes:
40 #include <Core/ITKCommon/common.hxx>
41 
42 
43 //----------------------------------------------------------------
44 // threshold_by_intensity_inplace
45 //
46 // The image is thresholded, intensities below/above threshold are set
47 // to a background value, specified by bg_offset.
48 //
49 // Depending on whether the image is being thresholded above or below
50 // the threshold, the background value is calculated as:
51 //
52 // a. keep pixels above threshold,
53 // bg = clip - bg_offset * (max - clip)
54 //
55 // b. keep pixels below threshold,
56 // bg = clip + bg_offset * (clip_max - min)
57 //
58 // The function returns the background value of the thresholded image.
59 //
60 template <class image_t>
61 double
62 threshold_by_intensity_inplace(typename image_t::Pointer & image,
63  const double clip_param,
64  const bool keep_pixels_above_threshold = true,
65  const unsigned int bins = 256,
66  const double bg_offset = 1e-3)
67 {
68  typedef itk::ImageRegionConstIterator<image_t> iter_t;
69  typedef typename image_t::IndexType index_t;
70 
71  // first find minima/maxima of the image:
72  double v_min = std::numeric_limits<double>::max();
73  double v_max = -v_min;
74 
75  iter_t iter(image, image->GetLargestPossibleRegion());
76  for (iter.GoToBegin(); !iter.IsAtEnd(); ++iter)
77  {
78  const double & v = iter.Get();
79  v_min = std::min(v_min, v);
80  v_max = std::max(v_max, v);
81  }
82 
83  // calculate the min/max range:
84  const double v_rng = v_max - v_min;
85 
86  if (v_rng == 0.0)
87  {
88  // there are no peaks in this image:
89  return v_min;
90  }
91 
92  double clip = v_min + clip_param * v_rng;
93 
94  // threshold the peaks:
95  double background =
96  keep_pixels_above_threshold
97  ? clip - bg_offset * v_rng
98  : clip + bg_offset * v_rng;
99 
100  if (keep_pixels_above_threshold)
101  {
102  image = threshold<image_t>(image, clip, v_max, background, v_max);
103  }
104  else
105  {
106  image = threshold<image_t>(image, v_min, clip, v_min, background);
107  }
108 
109  return background;
110 }
111 
112 //----------------------------------------------------------------
113 // threshold_by_intensity
114 //
115 template <class image_t>
116 typename image_t::Pointer
117 threshold_by_intensity(const image_t * image,
118  const double clip_param,
119  const bool keep_pixels_above_threshold = true,
120  const unsigned int bins = 256,
121  const double bg_offset = 1e-3)
122 {
123  typename image_t::Pointer peaks = cast<image_t, image_t>(image);
124  threshold_by_intensity_inplace(peaks,
125  clip_param,
126  keep_pixels_above_threshold,
127  bins,
128  bg_offset);
129  return peaks;
130 }
131 
132 //----------------------------------------------------------------
133 // threshold_by_area_inplace
134 //
135 // The clip_param below refers to the number of pixels that fall
136 // below/above the threshold. Thus, the number of pixels above
137 // the maxima is 1 - clip_param.
138 //
139 // Thus, it is possible to specify a threshold value without
140 // knowing anything about the image.
141 //
142 // The image is thresholded, intensities below/above threshold are set
143 // to a background value, specified by bg_offset.
144 //
145 // Depending on whether the image is being thresholded above or below
146 // the threshold, the background value is calculated as:
147 //
148 // a. keep pixels above threshold,
149 // bg = clip - bg_offset * (max - clip)
150 //
151 // b. keep pixels below threshold,
152 // bg = clip + bg_offset * (clip_max - min)
153 //
154 // The function returns the background value of the thresholded image
155 //
156 template <class image_t>
157 double
158 threshold_by_area_inplace(typename image_t::Pointer & image,
159  const double clip_param,
160  const bool keep_pixels_above_threshold = true,
161  const unsigned int bins = 4096,
162  const double bg_offset = 1e-3)
163 {
164  typedef itk::ImageRegionConstIterator<image_t> iter_t;
165  typedef typename image_t::IndexType index_t;
166 
167  // first find minima/maxima of the image:
168  double v_min = std::numeric_limits<double>::max();
169  double v_max = -v_min;
170 
171  iter_t iter(image, image->GetLargestPossibleRegion());
172  for (iter.GoToBegin(); !iter.IsAtEnd(); ++iter)
173  {
174  const double v = iter.Get();
175  v_min = std::min(v_min, v);
176  v_max = std::max(v_max, v);
177  }
178 
179  // calculate the min/max range:
180  const double v_rng = v_max - v_min;
181 
182  if (v_rng == 0.0)
183  {
184  // there are no peaks in this image:
185  return v_min;
186  }
187 
188  // build a histogram:
189  std::vector<unsigned int> pdf;
190  pdf.assign(bins, 0);
191 
192  for (iter.GoToBegin(); !iter.IsAtEnd(); ++iter)
193  {
194  const double v = iter.Get();
195  const unsigned int bin =
196  (unsigned int)(double((v - v_min) / v_rng) * double(bins - 1));
197  pdf[bin]++;
198  }
199 
200  // build the cumulative histogram:
201  std::vector<unsigned int> cdf(bins);
202  cdf[0] = pdf[0];
203  for (unsigned int i = 1; i < bins; i++)
204  {
205  cdf[i] = cdf[i - 1] + pdf[i];
206  }
207 
208  // calculate the total number of pixels in the image:
209  typename image_t::SizeType size =
210  image->GetLargestPossibleRegion().GetSize();
211 
212  double total_number_of_pixels = 1.0;
213  for (unsigned int i = 0; i < size.GetSizeDimension(); i++)
214  {
215  total_number_of_pixels *= size[i];
216  }
217 
218  // find the CDF bin that contains a given percentage of the total image:
219  double clip = 0.0;
220  for (unsigned int i = 1; i < bins; i++)
221  {
222  clip = v_min + (double(i) / double(bins - 1)) * v_rng;
223  if (double(cdf[i]) >= clip_param * total_number_of_pixels)
224  {
225  break;
226  }
227  }
228 
229  // threshold the peaks:
230  double background =
231  keep_pixels_above_threshold
232  ? clip - bg_offset * v_rng
233  : clip + bg_offset * v_rng;
234 
235  if (keep_pixels_above_threshold)
236  {
237  image = threshold<image_t>(image, clip, v_max, background, v_max);
238  }
239  else
240  {
241  image = threshold<image_t>(image, v_min, clip, v_min, background);
242  }
243 
244  return background;
245 }
246 
247 //----------------------------------------------------------------
248 // threshold_by_area
249 //
250 template <class image_t>
251 typename image_t::Pointer
252 threshold_by_area(const image_t * image,
253  const double clip_param,
254  const bool keep_pixels_above_threshold = true,
255  const unsigned int bins = 256,
256  const double bg_offset = 1e-3)
257 {
258  typename image_t::Pointer peaks = cast<image_t, image_t>(image);
259  threshold_by_area_inplace(peaks,
260  clip_param,
261  keep_pixels_above_threshold,
262  bins,
263  bg_offset);
264  return peaks;
265 }
266 
267 
268 #endif // THRESHOLD_HXX_