1 /*
2 * Medical Image Registration ToolKit (MIRTK)
3 *
4 * Copyright 2008-2017 Imperial College London
5 * Copyright 2008-2015 Daniel Rueckert, Julia Schnabel
6 * Copyright 2015-2017 Andreas Schuh
7 *
8 * Licensed under the Apache License, Version 2.0 (the "License");
9 * you may not use this file except in compliance with the License.
10 * You may obtain a copy of the License at
11 *
12 * http://www.apache.org/licenses/LICENSE-2.0
13 *
14 * Unless required by applicable law or agreed to in writing, software
15 * distributed under the License is distributed on an "AS IS" BASIS,
16 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
17 * See the License for the specific language governing permissions and
18 * limitations under the License.
19 */
20
21 #ifndef MIRTK_Histogram2D_H
22 #define MIRTK_Histogram2D_H
23
24 #include "mirtk/Math.h"
25 #include "mirtk/Object.h"
26 #include "mirtk/Histogram1D.h"
27
28
29 namespace mirtk {
30
31
32 /**
33 * Class for 2D histograms.
34 */
35 template <class HistogramType>
36 class Histogram2D : public Object
37 {
38 mirtkObjectMacro(Histogram2D);
39
40 public:
41
42 typedef HistogramType BinType;
43
44 private:
45
46 /// Number of bins in x-direction
47 int _nbins_x;
48
49 /// Number of bins in x-direction
50 int _nbins_y;
51
52 /// Number of samples
53 HistogramType _nsamp;
54
55 /// Min. value for x-samples, everything below is ignored
56 double _min_x;
57
58 /// Min. value for y-samples, everything below is ignored
59 double _min_y;
60
61 /// Max. value for x-samples, everything above is ignored
62 double _max_x;
63
64 /// Max. value for y-samples, everything above is ignored
65 double _max_y;
66
67 /// Width of bin in x-direction
68 double _width_x;
69
70 /// Width of bin in y-direction
71 double _width_y;
72
73 /// Dynamic memory for bins
74 HistogramType **_bins;
75
76 public:
77
78 /// Construct a histogram from another histogram
79 Histogram2D(const Histogram2D &);
80
81 /// Construct a histogram with 256 bins and samples ranging from 0 to 255
82 Histogram2D(int nbins_x = 256, int nbins_y = 256);
83
84 /// Construct a histogram for samples ranging from min to max and width
85 Histogram2D(double min_x, double max_x, double width_x,
86 double min_y, double max_y, double width_y);
87
88 /// Destructor
89 ~Histogram2D();
90
91 /// Assignment operator
92 Histogram2D &operator =(const Histogram2D &);
93
94 /// Construct a histogram for samples ranging from min to max and width
95 void Initialize(double min_x, double max_x, double width_x,
96 double min_y, double max_y, double width_y);
97
98 /// Clear and reset histogram
99 void Reset();
100
101 /// Get raw pointer to histogram bins
102 HistogramType *RawPointer();
103
104 /// Get raw pointer to histogram bins
105 const HistogramType *RawPointer() const;
106
107 /// Clear and copy histogram
108 void Reset(const Histogram2D &);
109
110 /// Transpose histogram in place with x and y direction exchanged
111 Histogram2D<HistogramType> &Transpose();
112
113 /// Get transposed histogram with x and y direction exchanged
114 Histogram2D<HistogramType> Transposed() const;
115
116 /// Get total number of bins
117 int NumberOfBins() const;
118
119 /// Get number of bins in x-direction
120 int NumberOfBinsX() const;
121
122 /// Put number of bins in x-direction
123 void PutNumberOfBinsX(int);
124
125 /// Get number of bins in x-direction
126 int NumberOfBinsY() const;
127
128 /// Put number of bins in x-direction
129 void PutNumberOfBinsY(int);
130
131 /// Get number of bins in x- and y-direction
132 void GetNumberOfBins(int *, int *) const;
133
134 /// Get number of bins in x- and y-direction
135 void GetNumberOfBins(int &, int &) const;
136
137 /// Put number of bins in x- and y-direction
138 void PutNumberOfBins(int, int);
139
140 /// Get minimum value in histogram in x-direction
141 double MinX() const;
142
143 /// Get maximum value in histogram in x-direction
144 double MaxX() const;
145
146 /// Get minimum value in histogram in y-direction
147 double MinY() const;
148
149 /// Get maximum value in histogram in y-direction
150 double MaxY() const;
151
152 /// Get width of bins in histogram in x-direction
153 double WidthX() const;
154
155 /// Get width of bins in histogram in y-direction
156 double WidthY() const;
157
158 /// Get minimum value in histogram
159 void GetMin(double *, double *) const;
160
161 /// Get minimum value in histogram
162 void GetMin(double &, double &) const;
163
164 /// Put minimum value in histogram
165 void PutMin(double, double);
166
167 /// Get maximum value in histogram
168 void GetMax(double *, double *) const;
169
170 /// Get maximum value in histogram
171 void GetMax(double &, double &) const;
172
173 /// Put maximum value in histogram
174 void PutMax(double, double);
175
176 /// Get width of bins in histogram
177 void GetWidth(double *, double *) const;
178
179 /// Get width of bins in histogram
180 void GetWidth(double &, double &) const;
181
182 /// Put width of bins in histogram
183 void PutWidth(double, double);
184
185 /// Get number of samples in histogram
186 HistogramType NumberOfSamples() const;
187
188 /// Set number of samples in histogram
189 void NumberOfSamples(HistogramType);
190
191 /// Get number of samples in bin(i, j)
192 HistogramType &operator()(int, int);
193
194 /// Get number of samples in bin(i, j)
195 const HistogramType &operator()(int, int) const;
196
197 /// Add counts to bins
198 void Add(int, int, HistogramType = 1);
199
200 /// Delete counts from bins
201 void Delete(int, int, HistogramType = 1);
202
203 /// Add samples
204 void AddSample(double, double, HistogramType = 1);
205
206 /// Delete samples
207 void DelSample(double, double, HistogramType = 1);
208
209 /// Convert sample value to bin index
210 int ValToBinX(double val) const;
211
212 /// Convert bin index to sample value
213 double BinToValX(int bin) const;
214
215 /// Convert sample value to bin index
216 int ValToBinY(double val) const;
217
218 /// Convert bin index sample value
219 double BinToValY(int bin) const;
220
221 /// Compute marginal histogram of X
222 void HistogramX(Histogram1D<HistogramType> &) const;
223
224 /// Compute marginal histogram of X
225 Histogram1D<HistogramType> HistogramX() const;
226
227 /// Compute marginal histogram of Y
228 void HistogramY(Histogram1D<HistogramType> &) const;
229
230 /// Compute marginal histogram of Y
231 Histogram1D<HistogramType> HistogramY() const;
232
233 /// Log transform histogram
234 void Log();
235
236 /// Smooth histogram
237 void Smooth();
238
239 /// Get smoothed histogram, optionally with padded boundaries
240 Histogram2D<HistogramType> Smoothed(bool = false);
241
242 /// Calculate joint probability p(x, y)
243 double JointProbability(int, int) const;
244
245 /// Calculate marginal probability p(x)
246 double MarginalProbabilityX(int) const;
247
248 /// Calculate marginal probability p(y)
249 double MarginalProbabilityY(int) const;
250
251 /// Calculate conditional probability p(x|y)
252 double ConditionalProbabilityXY(int, int) const;
253
254 /// Calculate conditional probability p(y|x)
255 double ConditionalProbabilityYX(int, int) const;
256
257 /// Calculate mean
258 double MeanX() const;
259
260 /// Calculate mean
261 double MeanY() const;
262
263 /// Calculate conditional mean
264 double ConditionalMeanXY(int) const;
265
266 /// Calculate conditional mean
267 double ConditionalMeanYX(int) const;
268
269 /// Calculate variance
270 double VarianceX() const;
271
272 /// Calculate variance
273 double VarianceY() const;
274
275 /// Calculate conditional variance
276 double ConditionalVarianceXY(int) const;
277
278 /// Calculate conditional variance
279 double ConditionalVarianceYX(int) const;
280
281 /// Calculate covariance
282 double Covariance() const;
283
284 /// Calculate standard deviation
285 double StandardDeviationX() const;
286
287 /// Calculate standard deviation
288 double StandardDeviationY() const;
289
290 /// Calculate marginal entropy
291 double EntropyX() const;
292
293 /// Calculate marginal entropy
294 double EntropyY() const;
295
296 /// Calculate joint entropy
297 double JointEntropy() const;
298
299 /// Calculate mutual information
300 double MutualInformation() const;
301
302 /// Calculate normalized mutual information
303 double NormalizedMutualInformation() const;
304
305 /// Calculate cross correlation
306 double CrossCorrelation() const;
307
308 /// Calculate correlation ratio
309 double CorrelationRatioXY() const;
310
311 /// Calculate correlation ratio
312 double CorrelationRatioYX() const;
313
314 /// Calculate sums of squared differences
315 double SumsOfSquaredDifferences() const;
316
317 /// Calcualate label consistency
318 double LabelConsistency() const;
319
320 /// Calcualate kappa statistic
321 double Kappa() const;
322
323 /// Read histogram
324 void Read(const char *);
325
326 /// Write histogram
327 void Write(const char *) const;
328
329 /// Write histogram as 2D image
330 void WriteAsImage(const char *) const;
331
332 /// Print histogram
333 void Print() const;
334 };
335
336 ////////////////////////////////////////////////////////////////////////////////
337 // Inline definitions
338 ////////////////////////////////////////////////////////////////////////////////
339
340 // -----------------------------------------------------------------------------
341 template <class HistogramType>
342 inline Histogram2D<HistogramType> &
343 Histogram2D<HistogramType>::operator =(const Histogram2D &h)
344 {
345 Reset(h);
346 return *this;
347 }
348
349 // -----------------------------------------------------------------------------
350 template <class HistogramType>
NumberOfBins()351 inline int Histogram2D<HistogramType>::NumberOfBins() const
352 {
353 return _nbins_x * _nbins_y;
354 }
355
356 // -----------------------------------------------------------------------------
357 template <class HistogramType>
NumberOfBinsX()358 inline int Histogram2D<HistogramType>::NumberOfBinsX() const
359 {
360 return _nbins_x;
361 }
362
363 // -----------------------------------------------------------------------------
364 template <class HistogramType>
NumberOfBinsY()365 inline int Histogram2D<HistogramType>::NumberOfBinsY() const
366 {
367 return _nbins_y;
368 }
369
370 // -----------------------------------------------------------------------------
371 template <class HistogramType>
RawPointer()372 inline HistogramType *Histogram2D<HistogramType>::RawPointer()
373 {
374 return _bins[0];
375 }
376
377 // -----------------------------------------------------------------------------
378 template <class HistogramType>
RawPointer()379 inline const HistogramType *Histogram2D<HistogramType>::RawPointer() const
380 {
381 return _bins[0];
382 }
383
384 // -----------------------------------------------------------------------------
385 template <class HistogramType>
NumberOfSamples()386 inline HistogramType Histogram2D<HistogramType>::NumberOfSamples() const
387 {
388 return _nsamp;
389 }
390
391 // -----------------------------------------------------------------------------
392 template <class HistogramType>
NumberOfSamples(HistogramType n)393 inline void Histogram2D<HistogramType>::NumberOfSamples(HistogramType n)
394 {
395 _nsamp = n;
396 }
397
398 // -----------------------------------------------------------------------------
399 template <class HistogramType>
MinX()400 inline double Histogram2D<HistogramType>::MinX() const
401 {
402 return _min_x;
403 }
404
405 // -----------------------------------------------------------------------------
406 template <class HistogramType>
MaxX()407 inline double Histogram2D<HistogramType>::MaxX() const
408 {
409 return _max_x;
410 }
411
412 // -----------------------------------------------------------------------------
413 template <class HistogramType>
MinY()414 inline double Histogram2D<HistogramType>::MinY() const
415 {
416 return _min_y;
417 }
418
419 // -----------------------------------------------------------------------------
420 template <class HistogramType>
MaxY()421 inline double Histogram2D<HistogramType>::MaxY() const
422 {
423 return _max_y;
424 }
425
426 // -----------------------------------------------------------------------------
427 template <class HistogramType>
WidthX()428 inline double Histogram2D<HistogramType>::WidthX() const
429 {
430 return _width_x;
431 }
432
433 // -----------------------------------------------------------------------------
434 template <class HistogramType>
WidthY()435 inline double Histogram2D<HistogramType>::WidthY() const
436 {
437 return _width_y;
438 }
439
440 // -----------------------------------------------------------------------------
441 template <class HistogramType>
operator()442 inline HistogramType &Histogram2D<HistogramType>::operator()(int i, int j)
443 {
444 return _bins[j][i];
445 }
446
447 // -----------------------------------------------------------------------------
448 template <class HistogramType>
operator()449 inline const HistogramType &Histogram2D<HistogramType>::operator()(int i, int j) const
450 {
451 return _bins[j][i];
452 }
453
454 // -----------------------------------------------------------------------------
455 template <class HistogramType>
Add(int i,int j,HistogramType n)456 inline void Histogram2D<HistogramType>::Add(int i, int j, HistogramType n)
457 {
458 _bins[j][i] += n;
459 _nsamp += n;
460 }
461
462 // -----------------------------------------------------------------------------
463 template <class HistogramType>
Delete(int i,int j,HistogramType n)464 inline void Histogram2D<HistogramType>::Delete(int i, int j, HistogramType n)
465 {
466 _bins[j][i] -= n;
467 _nsamp -= n;
468 }
469
470 // -----------------------------------------------------------------------------
471 template <class HistogramType>
ValToBinX(double val)472 inline int Histogram2D<HistogramType>::ValToBinX(double val) const
473 {
474 int index = static_cast<int>(round(_nbins_x * (val - _min_x - 0.5*_width_x) / (_max_x - _min_x)));
475 if (index < 0) index = 0;
476 if (index > _nbins_x-1) index = _nbins_x - 1;
477 return index;
478 }
479
480 // -----------------------------------------------------------------------------
481 template <class HistogramType>
ValToBinY(double val)482 inline int Histogram2D<HistogramType>::ValToBinY(double val) const
483 {
484 int index = static_cast<int>(round(_nbins_y * (val - _min_y - 0.5*_width_y) / (_max_y - _min_y)));
485 if (index < 0) index = 0;
486 if (index > _nbins_y-1) index = _nbins_y - 1;
487 return index;
488 }
489
490 // -----------------------------------------------------------------------------
491 template <class HistogramType>
BinToValX(int i)492 inline double Histogram2D<HistogramType>::BinToValX(int i) const
493 {
494 return (i*(_max_x - _min_x)/_nbins_x + _min_x) + _width_x/2.0;
495 }
496
497 // -----------------------------------------------------------------------------
498 template <class HistogramType>
BinToValY(int i)499 inline double Histogram2D<HistogramType>::BinToValY(int i) const
500 {
501 return (i*(_max_y - _min_y)/_nbins_y + _min_y) + _width_y/2.0;
502 }
503
504 // -----------------------------------------------------------------------------
505 template <class HistogramType>
HistogramX()506 inline Histogram1D<HistogramType> Histogram2D<HistogramType>::HistogramX() const
507 {
508 Histogram1D<HistogramType> hx(_nbins_x);
509 HistogramX(hx);
510 return hx;
511 }
512
513 // -----------------------------------------------------------------------------
514 template <class HistogramType>
HistogramY()515 inline Histogram1D<HistogramType> Histogram2D<HistogramType>::HistogramY() const
516 {
517 Histogram1D<HistogramType> hy(_nbins_y);
518 HistogramY(hy);
519 return hy;
520 }
521
522 // -----------------------------------------------------------------------------
523 template <class HistogramType>
JointProbability(int i,int j)524 inline double Histogram2D<HistogramType>::JointProbability(int i, int j) const
525 {
526 return _bins[j][i] / (double) _nsamp;
527 }
528
529 // -----------------------------------------------------------------------------
530 template <class HistogramType>
MarginalProbabilityX(int i)531 inline double Histogram2D<HistogramType>::MarginalProbabilityX(int i) const
532 {
533 HistogramType n = 0;
534 for (int j = 0; j < _nbins_y; j++) {
535 n += _bins[j][i];
536 }
537 return n / (double) _nsamp;
538 }
539
540 // -----------------------------------------------------------------------------
541 template <class HistogramType>
MarginalProbabilityY(int i)542 inline double Histogram2D<HistogramType>::MarginalProbabilityY(int i) const
543 {
544 HistogramType n = 0;
545 for (int j = 0; j < _nbins_x; j++) {
546 n += _bins[i][j];
547 }
548 return static_cast<double>(n) / _nsamp;
549 }
550
551 // -----------------------------------------------------------------------------
552 template <class HistogramType>
ConditionalProbabilityXY(int i,int j)553 inline double Histogram2D<HistogramType>::ConditionalProbabilityXY(int i, int j) const
554 {
555 double p = this->MarginalProbabilityY(j);
556 if (p > 0) {
557 return this->JointProbability(i, j) / p;
558 } else {
559 return 0;
560 }
561 }
562
563 // -----------------------------------------------------------------------------
564 template <class HistogramType>
ConditionalProbabilityYX(int i,int j)565 inline double Histogram2D<HistogramType>::ConditionalProbabilityYX(int i, int j) const
566 {
567 double p = this->MarginalProbabilityX(j);
568 if (p > 0) {
569 return this->JointProbability(j, i) / p;
570 } else {
571 return 0;
572 }
573 }
574
575
576 } // namespace mirtk
577
578 #endif // MIRTK_Histogram2D_H
579