ndimage.hpp 16.7 KB
Newer Older
1
2
#ifndef MLN_CORE_IMAGE_NDIMAGE_HH
# define MLN_CORE_IMAGE_NDIMAGE_HH
edwin's avatar
edwin committed
3

4
5

# include <mln/core/image_base.hpp>
edwin's avatar
edwin committed
6
7
8
9
10
11
12
13
14

# include <mln/core/domain/box.hpp>
# include <mln/core/memory.hpp>
# include <mln/core/assert.hpp>
# include <mln/core/image_traits.hpp>
# include <mln/core/image_category.hpp>

# include <boost/range/iterator_range.hpp>

15
16
17
# include <mln/core/image/ndimage_iter.hpp>
# include <mln/core/image/ndimage_pixel_iterator.hpp>

edwin's avatar
edwin committed
18
19
namespace mln
{
20
21
22
23
24
25
26
  // FWD
  //template <typename I, typename T> struct ndimage_iter;
  //template <typename I, typename T> struct ndimage_rev_iter;
  //template <typename T, unsigned dim, typename I> struct ndimage_pixel_iterator;
  //template <typename T, unsigned dim, typename I> struct ndimage_rev_pixel_iterator;
  template <typename T, unsigned dim, typename I> struct ndimage_pixel;
  template <typename T, unsigned dim, typename E> struct ndimage_base;
edwin's avatar
edwin committed
27

28
29
30
31
32
33
34
/******************************************/
/****              Traits              ****/
/******************************************/


  template <typename T, unsigned dim, typename E>
  struct image_traits< ndimage_base<T, dim, E> >
edwin's avatar
edwin committed
35
  {
36
37
38
    typedef raw_image_tag               category;
    typedef std::true_type              accessible;
  };
edwin's avatar
edwin committed
39

40
41
42
43
44
45
/******************************************/
/****            Definition            ****/
/******************************************/

  namespace internal
  {
edwin's avatar
edwin committed
46
47
48
    template <typename T, unsigned dim>
    struct ndimage_data
    {
49
      ndimage_data(size_t* shape_, unsigned border, T v = T());
edwin's avatar
edwin committed
50
51
52
53
54
55
56
      ~ndimage_data();

      size_t shape[dim];
      size_t strides[dim];
      size_t nbytes;
      char*  buffer;

57
58
59
    private:
      ndimage_data(const ndimage_data&);
    };
edwin's avatar
edwin committed
60
61
62
63
  }


  template <typename T, unsigned dim, typename E>
64
65
66
67
  struct ndimage_base : image_base<E, point<short, dim>, T,
                                   ndimage_pixel<T, dim, E>,
                                   ndimage_pixel<const T, dim, const E>
                                   >
edwin's avatar
edwin committed
68
  {
69
70
71
    BOOST_CONCEPT_ASSERT((IterableImage<E>));
    BOOST_CONCEPT_ASSERT((AccessibleImage<E>));

edwin's avatar
edwin committed
72
73
74
75
  private:
    typedef ndimage_base<T, dim, E>                             this_type;
  public:
    // As an Image
76

edwin's avatar
edwin committed
77
    typedef box<short, dim>					domain_type;
78
79
    typedef point<short,dim>                                    site_type;
    typedef point<short,dim>                                    point_type;
edwin's avatar
edwin committed
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
    typedef ndimage_pixel<const T, dim, const E>                const_pixel_type;
    typedef ndimage_pixel<T, dim, E>                            pixel_type;
    typedef T							value_type;
    enum { ndim = dim};

    // As a ContainerImage
    typedef T&			reference;
    typedef const T&		const_reference;
    typedef T*			pointer;
    typedef const T*		const_pointer;
    typedef ptrdiff_t		difference_type;
    typedef size_t		size_type;


    // As an Image
95
96
    // \group Constructors
    // \{
edwin's avatar
edwin committed
97
98
    explicit ndimage_base(unsigned border = 3);
    explicit ndimage_base(const domain_type& domain, unsigned border = 3);
99
    // \}
edwin's avatar
edwin committed
100
101
102
103

    const domain_type&  domain() const;

    // As an ContainerImage
104
105
    // \group Point-wise access
    // \{
edwin's avatar
edwin committed
106
107
    reference operator() (site_type p);
    const_reference operator() (site_type p) const;
108
    // \}
edwin's avatar
edwin committed
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

    // FIXME move to base
    pixel_type pix_at(site_type p)
    {
      mln_precondition(domain_.has(p));
      pixel_type pix;
      pix.point_ = p;
      pix.ima_  = (E*)this;
      pix.ptr_ = (char*) & (operator () (p));
      return pix;
    }

    const_pixel_type pix_at(site_type p) const
    {
      mln_precondition(domain_.has(p));
124
      const_pixel_type pix((const E*) this);
edwin's avatar
edwin committed
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
      pix.point_ = p;
      pix.ptr_ = (char*) & (operator () (p));
      return pix;
    }


    // As a Resizable Image
    void resize(const domain_type& domain, unsigned border = 3);

    // As a RandomAccessImage
    reference at (difference_type n);
    const_reference at (difference_type n) const;

    // As an IterableImage
    typedef ndimage_iter<this_type, T>                                  value_iterator;
140
    typedef ndimage_iter<const this_type, const T>                      const_value_iterator;
edwin's avatar
edwin committed
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
    typedef boost::reverse_iterator<value_iterator>                     reverse_value_iterator;
    typedef boost::reverse_iterator<const_value_iterator>               const_reverse_value_iterator;
    typedef ndimage_pixel_iterator<T, dim, E>                           pixel_iterator;
    typedef ndimage_pixel_iterator<const T, dim, const E>               const_pixel_iterator;
    typedef boost::reverse_iterator<pixel_iterator>                     reverse_pixel_iterator;
    typedef boost::reverse_iterator<const_pixel_iterator>               const_reverse_pixel_iterator;

    typedef boost::iterator_range<value_iterator>                                   value_range;
    typedef boost::iterator_range<const_value_iterator>                       const_value_range;
    typedef boost::iterator_range<pixel_iterator>                                   pixel_range;
    typedef boost::iterator_range<const_pixel_iterator>                       const_pixel_range;


    value_range         values();
    const_value_range   values() const;
156

edwin's avatar
edwin committed
157
158
159

    pixel_range                 pixels();
    const_pixel_range           pixels() const;
160

edwin's avatar
edwin committed
161
162
163
164
165
166
167
168
169
170
171
172




    // As a Raw Image
    const size_t*       strides() const;
    int border() const { return border_; }

  protected:
    friend struct ndimage_pixel<T, dim, E>;
    friend struct ndimage_pixel<const T, dim, const E>;

173
174
    //friend struct ndimage_iter<this_type, T>;
    //friend struct ndimage_iter<this_type, const T>;
edwin's avatar
edwin committed
175

176
177
    //friend struct ndimage_pixel_iterator<T, dim, E>;
    //friend struct ndimage_pixel_iterator<const T, dim, const E>;
edwin's avatar
edwin committed
178
179
180


    domain_type	domain_;	///< Domain of image
181
    std::array<size_t, dim>	strides_;	///< Strides in bytes
edwin's avatar
edwin committed
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
    std::shared_ptr< internal::ndimage_data<T, dim> > data_;
    int		border_;
    char*	ptr_;           ///< Pointer to the first element
    char*	last_;          ///< Pointer to the last element
  };


  template <typename T>
  struct image3d : ndimage_base<T, 3, image3d<T> >
  {
  protected:
    typedef ndimage_base<T, 3, image3d<T> > base;
    typedef typename base::domain_type domain_type;

  public:
    explicit image3d (unsigned border = 3) : base (border) {}

    explicit image3d(const domain_type& domain, unsigned border = 3)
      : base(domain, border)
    {
    }

    image3d(short nslices, short nrows, short ncols, unsigned border = 3)
      : base( (box<short,3>) {{{0,0,0}},{{nslices, nrows, ncols}}}, border)
    {
    }
  };


  namespace internal
  {
    template <typename T, unsigned dim>
214
    ndimage_data<T, dim>::ndimage_data(size_t* shape_, unsigned border, T v)
edwin's avatar
edwin committed
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
    {
      for (unsigned i = 0; i < dim; ++i)
	shape[i] = shape_[i] + 2 * border;

      strides[dim-1] = sizeof(T);

      // Each row / page ... are 16 bytes aligned
      unsigned ndim = dim;
      if (ndim >= 2)
	{
	  strides[dim-2] = ((shape[dim-1] * sizeof(T)) & ~(size_t)15) + (size_t) 16;
	  for (int i = dim-3; i >= 0; --i)
	    strides[i] = shape[i+1] * strides[i+1];
	}

      // Allocate data
      nbytes = shape[0] * strides[0];
      buffer = (char*) mln::aligned_malloc(nbytes, border * sizeof(T));

      // Construct
      {
	char* ptr = buffer;
	unsigned nlines = 1;
	for (unsigned i = 0; i < dim-1; ++i)
	  nlines *= shape[i];
	unsigned nelements = shape[dim-1];

        if (dim > 1) {
          for (unsigned i = 0; i < nlines; ++i, ptr += strides[dim-2]) {
            T* p_ = (T*) ptr;
            for (unsigned j = 0; j < nelements; ++j, ++p_)
246
              new (p_) T(v);
edwin's avatar
edwin committed
247
248
249
250
          }
        } else {
          T* p_ = (T*) ptr;
          for (unsigned j = 0; j < nelements; ++j, ++p_)
251
            new (p_) T(v);
edwin's avatar
edwin committed
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
        }
      }
    }


    template <typename T, unsigned dim>
    ndimage_data<T, dim>::~ndimage_data()
    {
      char* ptr = buffer;
      unsigned nlines = 1;
      for (unsigned i = 0; i < dim-1; ++i)
	nlines *= shape[i];
      unsigned nelements = shape[dim-1];
      if (dim == 1) {
        for (unsigned k = 0; k < nelements; ++k)
          ((T*)ptr + k)->~T();
        //    std::destroy( (T*)ptr, (T*) ptr + nelements);
      } else {
        for (unsigned i = 0; i < nlines; ++i, ptr += strides[dim-2])
          for (unsigned k = 0; k < nelements; ++k)
            ((T*)ptr + k)->~T();
        //std::destroy((T*)ptr, (T*)ptr + nelements);
      }
      mln::aligned_free(buffer);
    }
  }


  template <typename T, unsigned dim, typename E>
  ndimage_base<T,dim,E>::ndimage_base(unsigned border)
    : domain_ (), border_ (border), ptr_ (NULL)
  {
    for (unsigned i = 0; i < dim; ++i){
      mln_postcondition(domain_.pmin[i] == 0);
      mln_postcondition(domain_.pmax[i] == 0);
    }
  }

  template <typename T, unsigned dim, typename E>
  ndimage_base<T,dim,E>::ndimage_base(const domain_type& domain, unsigned border)
    : domain_ (domain),
      border_ (border)
  {
    //mln_precondition(domain.size() > 0);
    site_type shp = domain.shape();
297
    point<size_t, dim> sz;
edwin's avatar
edwin committed
298
299
300
301
302
    sz = shp;

    // Compute strides size (in bytes)
    // The row stride is 16 bytes aligned
    data_.reset(new internal::ndimage_data<T, dim>(&(sz[0]), border));
303
    std::copy(data_->strides, data_->strides + dim, strides_.begin());
edwin's avatar
edwin committed
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322

    // Compute pointer at (0,0)
    ptr_ = data_->buffer;
    last_ = data_->buffer;
    for (unsigned i = 0; i < dim; ++i) {
      ptr_ += border * strides_[i];
      last_ += (border + sz[i] - 1) * strides_[i];
    }

  }

  template <typename T, unsigned dim, typename E>
  void
  ndimage_base<T,dim,E>::resize(const domain_type& domain, unsigned border)
  {
    domain_ = domain;
    border_ = border;

    site_type shp = domain.shape();
323
    point<size_t, dim> sz;
edwin's avatar
edwin committed
324
325
326
327
328
    sz = shp;

    // Compute strides size (in bytes)
    // The row stride is 16 bytes aligned
    data_.reset(new internal::ndimage_data<T, dim>(&(sz[0]), border));
329
    std::copy(data_->strides, data_->strides + dim, strides_.begin());
edwin's avatar
edwin committed
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389

    // Compute pointer at (0,0)
    ptr_ = data_->buffer;
    last_ = data_->buffer;
    for (unsigned i = 0; i < dim; ++i) {
      ptr_ += border * strides_[i];
      last_ += (border + sz[i] - 1) * strides_[i];
    }
  }

  template <typename T, unsigned dim, typename E>
  inline
  T&
  ndimage_base<T,dim,E>::operator() (site_type p)
  {
    mln_precondition(domain_.has(p));
    site_type q = p - domain_.pmin;

    char* ptr = ptr_;
    for (unsigned i = 0; i < dim; ++i)
      ptr += q[i] * strides_[i];
    return * (T*)ptr;
  }

  template <typename T, unsigned dim, typename E>
  inline
  const T&
  ndimage_base<T,dim,E>::operator() (site_type p) const
  {
    mln_precondition(domain_.has(p));
    site_type q = p - domain_.pmin;

    const char* ptr = ptr_;
    for (unsigned i = 0; i < dim; ++i)
      ptr += q[i] * strides_[i];
    return * (const T*)ptr;
  }

  template <typename T, unsigned dim, typename E>
  inline
  T&
  ndimage_base<T,dim,E>::at (difference_type n)
  {
    mln_precondition(0 <= n && n < data_->nbytes);
    return *reinterpret_cast<T*>(ptr_+n);
  }

  template <typename T, unsigned dim, typename E>
  inline
  const T&
  ndimage_base<T,dim,E>::at (difference_type n) const
  {
    mln_precondition(0 <= n && n < data_->nbytes);
    return *reinterpret_cast<const T*>(ptr_+n);
  }

  template <typename T, unsigned dim, typename E>
  const size_t*
  ndimage_base<T,dim,E>::strides () const
  {
390
    return &strides_[0];
edwin's avatar
edwin committed
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
  }

  // template <typename T, unsigned dim, typename E>
  // ptrdiff_t
  // ndimage_base<T,dim,E>::offset (point_type dp) const
  // {
  //   ptrdiff_t x = 0;
  //   for (unsigned i = 0; i < dim; ++i)
  //     x += p[i] * strides_[i];
  //   return x;
  // }

  template <typename T, unsigned dim, typename E>
  inline
  const typename ndimage_base<T,dim,E>::domain_type&
  ndimage_base<T,dim,E>::domain () const
  {
    return domain_;
  }

  /* -- Value range -- */

  template <typename T, unsigned dim, typename E>
  inline
  typename ndimage_base<T,dim,E>::const_value_range
  ndimage_base<T,dim,E>::values () const
  {
418
419
420
421
422
423
424
425
426
427
428
429
    point_type shp = domain_.shape();
    size_t sz = shp[0] * strides_[0];
    const_pixel_type pix_begin_(exact(this)), pix_end_(exact(this));
    pix_begin_.ptr_ = ptr_;
    pix_begin_.point_ = (point_type) {{0,}} ;
    pix_end_.ptr_ = ptr_ + sz;
    pix_end_.point_ = (point_type) {{0,}};
    pix_end_.point_[0] = shp[0];

    const_value_iterator begin_( pix_begin_,
                                 internal::make_point_visitor(shp),
                                 internal::make_strided_pointer_value_visitor(ptr_, shp, strides_.begin()));
edwin's avatar
edwin committed
430

431
432
433
434
435
436

    const_value_iterator end_(   pix_end_,
                                 internal::make_point_visitor(shp),
                                 internal::make_strided_pointer_value_visitor(ptr_ + sz, shp, strides_.begin()));

    return boost::make_iterator_range(begin_, end_);
edwin's avatar
edwin committed
437
438
439
440
441
442
443
  }

  template <typename T, unsigned dim, typename E>
  inline
  typename ndimage_base<T,dim,E>::value_range
  ndimage_base<T,dim,E>::values ()
  {
444
445
446
447
448
449
450
451
452
453
454
455
    point_type shp = domain_.shape();
    size_t sz = shp[0] * strides_[0];
    pixel_type pix_begin_(exact(this)), pix_end_(exact(this));
    pix_begin_.ptr_ = ptr_;
    pix_begin_.point_ = point_type ();
    pix_end_.ptr_ = ptr_ + sz;
    pix_end_.point_ = point_type ();
    pix_end_.point_[0] = shp[0];

    value_iterator begin_( pix_begin_,
                           internal::make_point_visitor(shp),
                           internal::make_strided_pointer_value_visitor(ptr_, shp, strides_.begin()));
edwin's avatar
edwin committed
456
457


458
459
460
461
462
463
464
465
    value_iterator end_(   pix_end_,
                           internal::make_point_visitor(shp),
                           internal::make_strided_pointer_value_visitor(ptr_ + sz, shp, strides_.begin()));

    return boost::make_iterator_range(begin_, end_);
  }

  /* -- Value range -- */
edwin's avatar
edwin committed
466
467

  template <typename T, unsigned dim, typename E>
468
  inline
edwin's avatar
edwin committed
469
  typename ndimage_base<T,dim,E>::const_pixel_range
470
  ndimage_base<T,dim,E>::pixels () const
edwin's avatar
edwin committed
471
  {
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
    point_type shp = domain_.shape();
    size_t sz = shp[0] * strides_[0];
    const_pixel_type pix_begin_(exact(this)), pix_end_(exact(this));
    pix_begin_.ptr_ = ptr_;
    pix_begin_.point_ = domain_.pmin;
    pix_end_.ptr_ = ptr_ + sz;
    pix_end_.point_ = domain_.pmin;
    pix_end_.point_[0] = domain_.pmax[0];

    const_pixel_iterator begin_( pix_begin_,
                                 internal::make_point_visitor(domain_.pmin, domain_.pmax),
                                 internal::make_strided_pointer_value_visitor(ptr_, shp, strides_.begin()));


    const_pixel_iterator end_(   pix_end_,
                                 internal::make_point_visitor(domain_.pmin, domain_.pmax),
                                 internal::make_strided_pointer_value_visitor(ptr_ + sz, shp, strides_.begin()));

    return boost::make_iterator_range(begin_, end_);
edwin's avatar
edwin committed
491
492
493
  }

  template <typename T, unsigned dim, typename E>
494
  inline
edwin's avatar
edwin committed
495
  typename ndimage_base<T,dim,E>::pixel_range
496
  ndimage_base<T,dim,E>::pixels ()
edwin's avatar
edwin committed
497
  {
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
    point_type shp = domain_.shape();
    size_t sz = shp[0] * strides_[0];
    pixel_type pix_begin_(exact(this)), pix_end_(exact(this));
    pix_begin_.ptr_ = ptr_;
    pix_begin_.point_ = domain_.pmin;
    pix_end_.ptr_ = ptr_ + sz;
    pix_end_.point_ = domain_.pmin;
    pix_end_.point_[0] = domain_.pmax[0];

    pixel_iterator begin_( pix_begin_,
                           internal::make_point_visitor(domain_.pmin, domain_.pmax),
                           internal::make_strided_pointer_value_visitor(ptr_, shp, strides_.begin()));

    pixel_iterator end_(   pix_end_,
                           internal::make_point_visitor(domain_.pmin, domain_.pmax),
                           internal::make_strided_pointer_value_visitor(ptr_ + sz, shp, strides_.begin()));

    return boost::make_iterator_range(begin_, end_);
edwin's avatar
edwin committed
516
517
518
  }


519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
  // /* -- Pixel range  */
  // template <typename T, unsigned dim, typename E>
  // inline
  // typename ndimage_base<T,dim,E>::const_pixel_range
  // ndimage_base<T,dim,E>::pixels() const
  // {
  //   size_t sz = (domain_.pmax[0] - domain_.pmin[0]) * strides_[0];
  //   point_type pend = domain_.pmin; pend[0] = domain_.pmax[0];
  //   return boost::make_iterator_range(const_pixel_iterator(reinterpret_cast<const E*>(this), domain_.pmin, ptr_),
  //                                     const_pixel_iterator(reinterpret_cast<const E*>(this), pend, ptr_ + sz));
  // }

  // template <typename T, unsigned dim, typename E>
  // inline
  // typename ndimage_base<T,dim,E>::pixel_range
  // ndimage_base<T,dim,E>::pixels()
  // {
  //   size_t sz = (domain_.pmax[0] - domain_.pmin[0]) * strides_[0];
  //   point_type pend = domain_.pmin; pend[0] = domain_.pmax[0];
  //   return boost::make_iterator_range(pixel_iterator(reinterpret_cast<E*>(this), domain_.pmin, ptr_),
  //                                     pixel_iterator(reinterpret_cast<E*>(this), pend, ptr_ + sz));
  // }


edwin's avatar
edwin committed
543
544
545
} // end of namespace mln


546
547

#endif // !MLN_CORE_IMAGE_NDIMAGE_HH