36 #ifndef VIGRA_MULTI_HISTOGRAMM_HXX 37 #define VIGRA_MULTI_HISTOGRAMM_HXX 39 #include "multi_array.hxx" 40 #include "tinyvector.hxx" 41 #include "multi_gridgraph.hxx" 42 #include "multi_convolution.hxx" 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,
56 MultiArrayView<DIM+2 , T_HIST> histogram
59 typedef typename Graph::NodeIt graph_scanner;
60 typedef typename Graph::Node Node;
61 typedef T_HIST ValueType;
62 typedef TinyVector<ValueType,CHANNELS> ChannelsVals;
64 const Graph g(image.shape());
65 const ChannelsVals nBins(bins);
68 for (graph_scanner n(g); n != lemon::INVALID; ++n){
70 ChannelsVals binIndex = image[node];
75 for(
size_t d=0;d<DIM;++d){
78 for(
size_t c=0;c<CHANNELS;++c){
79 const float fi = binIndex[c];
81 histCoord[DIM]=std::min(bi,static_cast<size_t>(bins-1));
83 histogram[histCoord]+=1.0;
88 Kernel1D<float> gauss,gaussBin;
89 gauss.initGaussian(sigma);
90 gaussBin.initGaussian(sigmaBin);
91 for(
size_t c=0;c<CHANNELS;++c){
94 MultiArrayView<DIM+1,T_HIST> histc = histogram.bindOuter(c);
97 ConvolutionOptions<DIM+1> opts;
98 TinyVector<double, DIM+1> sigmaVec(sigma);
99 sigmaVec[DIM] = sigmaBin;
100 opts.stdDev(sigmaVec);
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
120 typedef typename Graph::NodeIt graph_scanner;
121 typedef typename Graph::Node Node;
123 const Graph g(imageA.shape());
126 for (graph_scanner n(g); n != lemon::INVALID; ++n){
129 T_DATA binIndexA = imageA[node];
130 T_DATA binIndexB = imageA[node];
132 binIndexA -=minVals[0];
133 binIndexA /=maxVals[0];
134 binIndexA *=nBins[0];
136 binIndexB -=minVals[1];
137 binIndexB /=maxVals[1];
138 binIndexB *=nBins[1];
141 for(
size_t d=0;d<DIM;++d)
142 histCoord[d]=node[d];
144 histCoord[DIM]=binIndexA;
145 histCoord[DIM+1]=binIndexB;
147 const float fiA = binIndexA;
149 const float fiB = binIndexB;
151 histCoord[DIM]=std::min(biA,static_cast<unsigned int>(nBins[0]-1));
152 histCoord[DIM+1]=std::min(biB,static_cast<unsigned int>(nBins[1]-1));
153 histogram[histCoord]+=1.0;
157 MultiArray<DIM+2 , T_HIST> histogramBuffer(histogram);
158 Kernel1D<float> gaussS,gaussA,gaussB;
159 gaussS.initGaussian(sigma[0]);
160 gaussA.initGaussian(sigma[1]);
161 gaussB.initGaussian(sigma[2]);
175 histogram=histogramBuffer;
186 throw std::runtime_error(
"not yet implemented for arbitrary dimension");
194 template<
unsigned int DIM ,
class T,
class V,
class U>
195 void multiGaussianRankOrder(
196 const MultiArrayView<DIM, T > & image,
200 const TinyVector<double, DIM+1> sigmas,
201 const MultiArrayView<1, V> & ranks,
202 MultiArrayView<DIM+1, U> out
204 typedef MultiArray<DIM, T> ImgType;
205 typedef typename ImgType::difference_type ImgCoord;
207 typedef MultiArray<DIM+1, float> HistType;
208 typedef typename HistType::difference_type HistCoord;
210 typedef MultiArray<DIM+1, U> OutType;
211 typedef typename OutType::difference_type OutCoord;
216 std::copy(image.shape().begin(), image.shape().end(), histShape.begin());
217 histShape[DIM] = bins;
218 HistType histA(histShape);
223 HistCoord histCoord,nextHistCoord;
225 MultiCoordinateIterator<DIM> iter(image.shape());
226 for(std::ptrdiff_t i=0 ;i<image.size(); ++i, ++iter){
227 const ImgCoord imgCoord(*iter);
228 std::copy(imgCoord.begin(),imgCoord.end(),histCoord.begin() );
230 const T value = image[imgCoord];
231 const T fbinIndex = ((value-minVal)/(maxVal-minVal))*bins;
233 const int floorBin =
static_cast<int>(fFloorBin);
234 const int ceilBin =
static_cast<int>(
std::ceil(fbinIndex));
236 if(floorBin==ceilBin){
237 histCoord[DIM] = floorBin;
238 histA[histCoord] += 1.0;
243 const double ceilW = (fbinIndex - fFloorBin);
244 const double floorW = 1.0 - ceilW;
245 histCoord[DIM] = floorBin;
246 histA[histCoord] += floorW;
247 histCoord[DIM] = ceilBin;
248 histA[histCoord] += ceilW;
254 ConvolutionOptions<DIM+1> opts;
263 std::vector<float> histBuffer(bins);
266 MultiCoordinateIterator<DIM> iter(image.shape());
267 for(std::ptrdiff_t i=0 ;i<image.size(); ++i, ++iter){
270 const ImgCoord imgCoord(*iter);
273 std::copy(imgCoord.begin(),imgCoord.end(),histCoord.begin() );
274 nextHistCoord = histCoord;
275 std::copy(imgCoord.begin(),imgCoord.end(),outCoord.begin() );
277 for(
size_t bi=0; bi<bins; ++bi){
279 sum += histA[histCoord];
281 for(
size_t bi=0; bi<bins; ++bi){
283 histA[histCoord] /=
sum;
286 histBuffer[0] = histA[histCoord];
287 for(
size_t bi=1; bi<bins; ++bi){
289 double prevVal = histA[histCoord];
291 histA[histCoord] +=prevVal;
292 histBuffer[bi] = histA[histCoord];
298 for(std::ptrdiff_t r=0; r<ranks.size(); ++r){
300 const V rank = ranks[r];
302 nextHistCoord[DIM] = bi +1;
305 if(rank < histA[histCoord] ||
306 std::abs(rank-histA[histCoord])< 0.0000001 ||
309 out[outCoord] =
static_cast<U
>((maxVal-minVal)*bi*bins + minVal);
314 const size_t upperBinIndex =
315 std::distance(histBuffer.begin(),std::lower_bound(histBuffer.begin()+bi, histBuffer.end(),float(rank)));
316 bi = upperBinIndex - 1;
318 nextHistCoord[DIM] = upperBinIndex;
319 const double rankVal0 =
static_cast<U
>((maxVal-minVal)*bi*bins + minVal);
320 const double rankVal1 =
static_cast<U
>((maxVal-minVal)*(bi+1)*bins + minVal);
321 const double dd = histA[nextHistCoord] - histA[histCoord];
322 const double relD0 = (rank - histA[histCoord])/dd;
323 out[outCoord] = rankVal1 * relD0 + (1.0 - relD0)*rankVal0;
333 #endif //VIGRA_MULTI_HISTOGRAMM_HXX Define a grid graph in arbitrary dimensions.
Definition: multi_fwd.hxx:217
MultiArrayShape< actual_dimension >::type difference_type
Definition: multi_array.hxx:687
void gaussianSmoothMultiArray(...)
Isotropic Gaussian smoothing of a multi-dimensional arrays.
void convolveMultiArrayOneDimension(...)
Convolution along a single dimension of a multi-dimensional arrays.
Definition: accessor.hxx:43
NumericTraits< V >::Promote sum(TinyVectorBase< V, SIZE, D1, D2 > const &l)
sum of the vector's elements
Definition: tinyvector.hxx:2073
FFTWComplex< R >::NormType abs(const FFTWComplex< R > &a)
absolute value (= magnitude)
Definition: fftw3.hxx:1002
int ceil(FixedPoint< IntBits, FracBits > v)
rounding up.
Definition: fixedpoint.hxx:675
int floor(FixedPoint< IntBits, FracBits > v)
rounding down.
Definition: fixedpoint.hxx:667