[ VIGRA Homepage | Function Index | Class Index | Namespaces | File List | Main Page ]

multi_histogram.hxx VIGRA

1 /************************************************************************/
2 /* */
3 /* Copyright 2002-2003 by Ullrich Koethe, Hans Meine */
4 /* */
5 /* This file is part of the VIGRA computer vision library. */
6 /* The VIGRA Website is */
7 /* http://hci.iwr.uni-heidelberg.de/vigra/ */
8 /* Please direct questions, bug reports, and contributions to */
9 /* ullrich.koethe@iwr.uni-heidelberg.de or */
10 /* vigra@informatik.uni-hamburg.de */
11 /* */
12 /* Permission is hereby granted, free of charge, to any person */
13 /* obtaining a copy of this software and associated documentation */
14 /* files (the "Software"), to deal in the Software without */
15 /* restriction, including without limitation the rights to use, */
16 /* copy, modify, merge, publish, distribute, sublicense, and/or */
17 /* sell copies of the Software, and to permit persons to whom the */
18 /* Software is furnished to do so, subject to the following */
19 /* conditions: */
20 /* */
21 /* The above copyright notice and this permission notice shall be */
22 /* included in all copies or substantial portions of the */
23 /* Software. */
24 /* */
25 /* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND */
26 /* EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES */
27 /* OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND */
28 /* NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT */
29 /* HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, */
30 /* WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING */
31 /* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR */
32 /* OTHER DEALINGS IN THE SOFTWARE. */
33 /* */
34 /************************************************************************/
35 
36 #ifndef VIGRA_MULTI_HISTOGRAMM_HXX
37 #define VIGRA_MULTI_HISTOGRAMM_HXX
38 
39 #include "multi_array.hxx"
40 #include "tinyvector.hxx"
41 #include "multi_gridgraph.hxx"
42 #include "multi_convolution.hxx"
43 
44 
45 namespace vigra{
46 
47 
48  template< unsigned int DIM , class T_DATA ,unsigned int CHANNELS, class T_HIST >
49  void multiGaussianHistogram(
50  const MultiArrayView<DIM, TinyVector<T_DATA,CHANNELS> > & image,
51  const TinyVector<T_DATA,CHANNELS> minVals,
52  const TinyVector<T_DATA,CHANNELS> maxVals,
53  const size_t bins,
54  const float sigma,
55  const float sigmaBin,
56  MultiArrayView<DIM+2 , T_HIST> histogram
57  ){
58  typedef vigra::GridGraph< DIM , undirected_tag> Graph;
59  typedef typename Graph::NodeIt graph_scanner;
60  typedef typename Graph::Node Node;
61  typedef T_HIST ValueType;
62  typedef TinyVector<ValueType,CHANNELS> ChannelsVals;
63  typedef typename MultiArrayView<DIM+2 , T_HIST>::difference_type HistCoord;
64  const Graph g(image.shape());
65  const ChannelsVals nBins(bins);
66  histogram.init(1.0);
67  // iterate over all nodes (i.e. pixels)
68  for (graph_scanner n(g); n != lemon::INVALID; ++n){
69  const Node node(*n);
70  ChannelsVals binIndex = image[node];
71  binIndex -=minVals;
72  binIndex /=maxVals;
73  binIndex *=nBins;
74  HistCoord histCoord;
75  for(size_t d=0;d<DIM;++d){
76  histCoord[d]=node[d];
77  }
78  for(size_t c=0;c<CHANNELS;++c){
79  const float fi = binIndex[c];
80  const size_t bi = std::floor(fi+0.5);
81  histCoord[DIM]=std::min(bi,static_cast<size_t>(bins-1));
82  histCoord[DIM+1]=c;
83  histogram[histCoord]+=1.0;
84  }
85  }
86 
87  //MultiArray<DIM+2 , T_HIST> histogramBuffer(histogram);
88  Kernel1D<float> gauss,gaussBin;
89  gauss.initGaussian(sigma);
90  gaussBin.initGaussian(sigmaBin);
91  for(size_t c=0;c<CHANNELS;++c){
92 
93  // histogram for one channel
94  MultiArrayView<DIM+1,T_HIST> histc = histogram.bindOuter(c);
95  //MultiArrayView<DIM+1,T_HIST> histcBuffer = histogram.bindOuter(c);
96 
97  ConvolutionOptions<DIM+1> opts;
98  TinyVector<double, DIM+1> sigmaVec(sigma);
99  sigmaVec[DIM] = sigmaBin;
100  opts.stdDev(sigmaVec);
101 
102  // convolve spatial dimensions
103  gaussianSmoothMultiArray(histc, histc, opts);
104 
105  }
106 
107  }
108 
109  template< unsigned int DIM , class T_DATA, class T_HIST >
110  void multiGaussianCoHistogram(
111  const MultiArrayView<DIM, T_DATA > & imageA,
112  const MultiArrayView<DIM, T_DATA > & /*imageB*/,
113  const TinyVector<T_DATA,2> & minVals,
114  const TinyVector<T_DATA,2> & maxVals,
115  const TinyVector<int,2> & nBins,
116  const TinyVector<float,3> & sigma,
117  MultiArrayView<DIM+2, T_HIST> histogram
118  ){
119  typedef vigra::GridGraph< DIM , undirected_tag> Graph;
120  typedef typename Graph::NodeIt graph_scanner;
121  typedef typename Graph::Node Node;
122  typedef typename MultiArrayView<DIM+2 , T_HIST>::difference_type HistCoord;
123  const Graph g(imageA.shape());
124  histogram.init(0.0);
125  // iterate over all nodes (i.e. pixels)
126  for (graph_scanner n(g); n != lemon::INVALID; ++n){
127 
128  const Node node(*n);
129  T_DATA binIndexA = imageA[node];
130  T_DATA binIndexB = imageA[node];
131 
132  binIndexA -=minVals[0];
133  binIndexA /=maxVals[0];
134  binIndexA *=nBins[0];
135 
136  binIndexB -=minVals[1];
137  binIndexB /=maxVals[1];
138  binIndexB *=nBins[1];
139 
140  HistCoord histCoord;
141  for(size_t d=0;d<DIM;++d)
142  histCoord[d]=node[d];
143 
144  histCoord[DIM]=binIndexA;
145  histCoord[DIM+1]=binIndexB;
146 
147  const float fiA = binIndexA;
148  const unsigned int biA = std::floor(fiA+0.5);
149  const unsigned int biB = std::floor(fiA+0.5);
150  histCoord[DIM]=std::min(biA,static_cast<unsigned int>(nBins[0]-1));
151  histCoord[DIM+1]=std::min(biB,static_cast<unsigned int>(nBins[1]-1));
152  histogram[histCoord]+=1.0;
153 
154  }
155 
156  MultiArray<DIM+2 , T_HIST> histogramBuffer(histogram);
157  Kernel1D<float> gaussS,gaussA,gaussB;
158  gaussS.initGaussian(sigma[0]);
159  gaussA.initGaussian(sigma[1]);
160  gaussB.initGaussian(sigma[2]);
161 
162  if(DIM==2){
163  convolveMultiArrayOneDimension(srcMultiArrayRange(histogram), destMultiArray(histogramBuffer), 0, gaussS);
164  convolveMultiArrayOneDimension(srcMultiArrayRange(histogramBuffer), destMultiArray(histogram), 1, gaussS);
165  convolveMultiArrayOneDimension(srcMultiArrayRange(histogram), destMultiArray(histogramBuffer), 2, gaussA);
166  convolveMultiArrayOneDimension(srcMultiArrayRange(histogramBuffer), destMultiArray(histogram), 3, gaussB);
167  }
168  else if(DIM==3){
169  convolveMultiArrayOneDimension(srcMultiArrayRange(histogram), destMultiArray(histogramBuffer), 0, gaussS);
170  convolveMultiArrayOneDimension(srcMultiArrayRange(histogramBuffer), destMultiArray(histogram), 1, gaussS);
171  convolveMultiArrayOneDimension(srcMultiArrayRange(histogram), destMultiArray(histogramBuffer), 2, gaussS);
172  convolveMultiArrayOneDimension(srcMultiArrayRange(histogramBuffer), destMultiArray(histogram), 3, gaussA);
173  convolveMultiArrayOneDimension(srcMultiArrayRange(histogram), destMultiArray(histogramBuffer), 4, gaussB);
174  histogram=histogramBuffer;
175  }
176  else if(DIM==4){
177  convolveMultiArrayOneDimension(srcMultiArrayRange(histogram), destMultiArray(histogramBuffer), 0, gaussS);
178  convolveMultiArrayOneDimension(srcMultiArrayRange(histogramBuffer), destMultiArray(histogram), 1, gaussS);
179  convolveMultiArrayOneDimension(srcMultiArrayRange(histogram), destMultiArray(histogramBuffer), 2, gaussS);
180  convolveMultiArrayOneDimension(srcMultiArrayRange(histogramBuffer), destMultiArray(histogram), 3, gaussS);
181  convolveMultiArrayOneDimension(srcMultiArrayRange(histogram), destMultiArray(histogramBuffer), 4, gaussA);
182  convolveMultiArrayOneDimension(srcMultiArrayRange(histogramBuffer), destMultiArray(histogram), 5, gaussA);
183  }
184  else{
185  throw std::runtime_error("not yet implemented for arbitrary dimension");
186  }
187 
188  }
189 
190 
191 
192 
193  template< unsigned int DIM , class T, class V, class U>
194  void multiGaussianRankOrder(
195  const MultiArrayView<DIM, T > & image,
196  const T minVal,
197  const T maxVal,
198  const size_t bins,
199  const TinyVector<double, DIM+1> sigmas,
200  const MultiArrayView<1, V> & ranks,
201  MultiArrayView<DIM+1, U> out
202  ){
203  typedef MultiArray<DIM, T> ImgType;
204  typedef typename ImgType::difference_type ImgCoord;
205 
206  typedef MultiArray<DIM+1, float> HistType;
207  typedef typename HistType::difference_type HistCoord;
208 
209  typedef MultiArray<DIM+1, U> OutType;
210  typedef typename OutType::difference_type OutCoord;
211 
212  // FIXME: crashes on Python3
213 
214  HistCoord histShape;
215  std::copy(image.shape().begin(), image.shape().end(), histShape.begin());
216  histShape[DIM] = bins;
217  HistType histA(histShape);
218 
219  histA = 0.0;
220 
221  // collect values
222  HistCoord histCoord,nextHistCoord;
223  {
224  MultiCoordinateIterator<DIM> iter(image.shape());
225  for(std::ptrdiff_t i=0 ;i<image.size(); ++i, ++iter){
226  const ImgCoord imgCoord(*iter);
227  std::copy(imgCoord.begin(),imgCoord.end(),histCoord.begin() );
228 
229  const T value = image[imgCoord];
230  const T fbinIndex = ((value-minVal)/(maxVal-minVal))*bins;
231  const T fFloorBin = std::floor(fbinIndex);
232  const int floorBin = static_cast<int>(fFloorBin);
233  const int ceilBin = static_cast<int>(std::ceil(fbinIndex));
234 
235  if(floorBin==ceilBin){
236  histCoord[DIM] = floorBin;
237  histA[histCoord] += 1.0;
238  }
239  else{
240  const T floorBin = std::floor(fbinIndex);
241  const T ceilBin = std::ceil(fbinIndex);
242  const double ceilW = (fbinIndex - fFloorBin);
243  const double floorW = 1.0 - ceilW;
244  histCoord[DIM] = floorBin;
245  histA[histCoord] += floorW;
246  histCoord[DIM] = ceilBin;
247  histA[histCoord] += ceilW;
248  }
249 
250  }
251  }
252  //
253  ConvolutionOptions<DIM+1> opts;
254  opts.stdDev(sigmas);
255 
256  // convolve spatial dimensions
257  gaussianSmoothMultiArray(histA, histA, opts);
258 
259  OutCoord outCoord;
260 
261 
262  std::vector<float> histBuffer(bins);
263  //std::cout<<"normalize and compute ranks\n";
264  {
265  MultiCoordinateIterator<DIM> iter(image.shape());
266  for(std::ptrdiff_t i=0 ;i<image.size(); ++i, ++iter){
267 
268  // normalize
269  const ImgCoord imgCoord(*iter);
270  //std::cout<<"at pixel "<<imgCoord<<"\n";
271 
272  std::copy(imgCoord.begin(),imgCoord.end(),histCoord.begin() );
273  nextHistCoord = histCoord;
274  std::copy(imgCoord.begin(),imgCoord.end(),outCoord.begin() );
275  double sum = 0;
276  for(size_t bi=0; bi<bins; ++bi){
277  histCoord[DIM] = bi;
278  sum += histA[histCoord];
279  }
280  for(size_t bi=0; bi<bins; ++bi){
281  histCoord[DIM] = bi;
282  histA[histCoord] /= sum;
283  }
284  histCoord[DIM] = 0;
285  histBuffer[0] = histA[histCoord];
286  for(size_t bi=1; bi<bins; ++bi){
287 
288  double prevVal = histA[histCoord];
289  histCoord[DIM] = bi;
290  histA[histCoord] +=prevVal;
291  histBuffer[bi] = histA[histCoord];
292  }
293 
294 
295 
296  size_t bi=0;
297  for(std::ptrdiff_t r=0; r<ranks.size(); ++r){
298  outCoord[DIM] = r;
299  const V rank = ranks[r];
300  histCoord[DIM] = bi;
301  nextHistCoord[DIM] = bi +1;
302  //std::cout<<" bi "<<bi<<" rank "<<rank<<" "<<histA[histCoord]<<"\n";
303  // corner cases
304  if(rank < histA[histCoord] ||
305  std::abs(rank-histA[histCoord])< 0.0000001 ||
306  bi==bins-1
307  ){
308  out[outCoord] = static_cast<U>((maxVal-minVal)*bi*bins + minVal);
309  break;
310  }
311  else{
312  // with binary search
313  const size_t upperBinIndex =
314  std::distance(histBuffer.begin(),std::lower_bound(histBuffer.begin()+bi, histBuffer.end(),float(rank)));
315  bi = upperBinIndex - 1;
316  histCoord[DIM] = bi;
317  nextHistCoord[DIM] = upperBinIndex;
318  const double rankVal0 = static_cast<U>((maxVal-minVal)*bi*bins + minVal);
319  const double rankVal1 = static_cast<U>((maxVal-minVal)*(bi+1)*bins + minVal);
320  const double dd = histA[nextHistCoord] - histA[histCoord];
321  const double relD0 = (rank - histA[histCoord])/dd;
322  out[outCoord] = rankVal1 * relD0 + (1.0 - relD0)*rankVal0;
323  break;
324  }
325  }
326  }
327  }
328  }
329 }
330 //end namespace vigra
331 
332 #endif //VIGRA_MULTI_HISTOGRAMM_HXX
vigra::GridGraph
Define a grid graph in arbitrary dimensions.
Definition: multi_fwd.hxx:217
vigra::sum
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:2073
vigra::abs
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
vigra::convolveMultiArrayOneDimension
void convolveMultiArrayOneDimension(...)
Convolution along a single dimension of a multi-dimensional arrays.
vigra::gaussianSmoothMultiArray
void gaussianSmoothMultiArray(...)
Isotropic Gaussian smoothing of a multi-dimensional arrays.
vigra::ceil
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition: fixedpoint.hxx:675
vigra::floor
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition: fixedpoint.hxx:667

© Ullrich Köthe (ullrich.koethe@iwr.uni-heidelberg.de)
Heidelberg Collaboratory for Image Processing, University of Heidelberg, Germany

html generated using doxygen and Python
vigra 1.11.1 (Fri May 19 2017)