CImg<ulongT> _label(const unsigned int nb, const int *const dx,
                        const int *const dy, const int *const dz,
                        const Tfloat tolerance, const bool is_L2_norm) const {
      CImg<ulongT> res(_width,_height,_depth);
      const Tfloat _tolerance = _spectrum>1 && is_L2_norm?cimg::sqr(tolerance):tolerance;

      // Init label numbers.
      ulongT *ptr = res.data();
      cimg_foroff(res,p) *(ptr++) = p;

      // For each neighbour-direction, label.
      for (unsigned int n = 0; n<nb; ++n) {
        const int _dx = dx[n], _dy = dy[n], _dz = dz[n];
        if (_dx || _dy || _dz) {
          const int
            x0 = _dx<0?-_dx:0,
            x1 = _dx<0?width():width() - _dx,
            y0 = _dy<0?-_dy:0,
            y1 = _dy<0?height():height() - _dy,
            z0 = _dz<0?-_dz:0,
            z1 = _dz<0?depth():depth() - _dz;
          const longT
            wh = (longT)width()*height(),
            whd = (longT)width()*height()*depth(),
            offset = _dz*wh + _dy*width() + _dx;
          for (longT z = z0, nz = z0 + _dz, pz = z0*wh; z<z1; ++z, ++nz, pz+=wh) {
            for (longT y = y0, ny = y0 + _dy, py = y0*width() + pz; y<y1; ++y, ++ny, py+=width()) {
              for (longT x = x0, nx = x0 + _dx, p = x0 + py; x<x1; ++x, ++nx, ++p) {
                Tfloat diff;
                switch (_spectrum) {
                case 1 :
                  diff = cimg::abs((Tfloat)(*this)(x,y,z,0,wh,whd) - (Tfloat)(*this)(nx,ny,nz,0,wh,whd));
                  break;
                case 2 :
                  if (is_L2_norm)
                    diff = cimg::sqr((Tfloat)(*this)(x,y,z,0,wh,whd) - (Tfloat)(*this)(nx,ny,nz,0,wh,whd)) +
                      cimg::sqr((Tfloat)(*this)(x,y,z,1,wh,whd) - (Tfloat)(*this)(nx,ny,nz,1,wh,whd));
                  else
                    diff = cimg::abs((Tfloat)(*this)(x,y,z,0,wh,whd) - (Tfloat)(*this)(nx,ny,nz,0,wh,whd)) +
                      cimg::abs((Tfloat)(*this)(x,y,z,1,wh,whd) - (Tfloat)(*this)(nx,ny,nz,1,wh,whd));
                  break;
                case 3 :
                  if (is_L2_norm)
                    diff = cimg::sqr((Tfloat)(*this)(x,y,z,0,wh,whd) - (Tfloat)(*this)(nx,ny,nz,0,wh,whd)) +
                      cimg::sqr((Tfloat)(*this)(x,y,z,1,wh,whd) - (Tfloat)(*this)(nx,ny,nz,1,wh,whd)) +
                      cimg::sqr((Tfloat)(*this)(x,y,z,2,wh,whd) - (Tfloat)(*this)(nx,ny,nz,2,wh,whd));
                  else
                    diff = cimg::abs((Tfloat)(*this)(x,y,z,0,wh,whd) - (Tfloat)(*this)(nx,ny,nz,0,wh,whd)) +
                      cimg::abs((Tfloat)(*this)(x,y,z,1,wh,whd) - (Tfloat)(*this)(nx,ny,nz,1,wh,whd)) +
                      cimg::abs((Tfloat)(*this)(x,y,z,2,wh,whd) - (Tfloat)(*this)(nx,ny,nz,2,wh,whd));
                  break;
                case 4 :
                  if (is_L2_norm)
                    diff = cimg::sqr((Tfloat)(*this)(x,y,z,0,wh,whd) - (Tfloat)(*this)(nx,ny,nz,0,wh,whd)) +
                      cimg::sqr((Tfloat)(*this)(x,y,z,1,wh,whd) - (Tfloat)(*this)(nx,ny,nz,1,wh,whd)) +
                      cimg::sqr((Tfloat)(*this)(x,y,z,2,wh,whd) - (Tfloat)(*this)(nx,ny,nz,2,wh,whd)) +
                      cimg::sqr((Tfloat)(*this)(x,y,z,3,wh,whd) - (Tfloat)(*this)(nx,ny,nz,3,wh,whd));
                  else
                    diff = cimg::abs((Tfloat)(*this)(x,y,z,0,wh,whd) - (Tfloat)(*this)(nx,ny,nz,0,wh,whd)) +
                      cimg::abs((Tfloat)(*this)(x,y,z,1,wh,whd) - (Tfloat)(*this)(nx,ny,nz,1,wh,whd)) +
                      cimg::abs((Tfloat)(*this)(x,y,z,2,wh,whd) - (Tfloat)(*this)(nx,ny,nz,2,wh,whd)) +
                      cimg::abs((Tfloat)(*this)(x,y,z,3,wh,whd) - (Tfloat)(*this)(nx,ny,nz,3,wh,whd));
                  break;
                default :
                  diff = 0;
                  if (is_L2_norm)
                    cimg_forC(*this,c)
                      diff+=cimg::sqr((Tfloat)(*this)(x,y,z,c,wh,whd) - (Tfloat)(*this)(nx,ny,nz,c,wh,whd));
                  else
                    cimg_forC(*this,c)
                      diff+=cimg::abs((Tfloat)(*this)(x,y,z,c,wh,whd) - (Tfloat)(*this)(nx,ny,nz,c,wh,whd));
                }

                if (diff<=_tolerance) {
                  const longT q = p + offset;
                  ulongT xk, yk;
                  for (xk = (ulongT)(p<q?q:p), yk = (ulongT)(p<q?p:q); xk!=yk && res[xk]!=xk; ) {
                    xk = res[xk]; if (xk<yk) cimg::swap(xk,yk);
                  }
                  if (xk!=yk) res[xk] = (ulongT)yk;
                  for (ulongT _p = (ulongT)p; _p!=yk; ) {
                    const ulongT h = res[_p];
                    res[_p] = (ulongT)yk;
                    _p = h;
                  }
                  for (ulongT _q = (ulongT)q; _q!=yk; ) {
                    const ulongT h = res[_q];
                    res[_q] = (ulongT)yk;
                    _q = h;
                  }
                }
              }
            }
          }
        }
      }

      // Resolve equivalences.
      ulongT counter = 0;
      ptr = res.data();
      cimg_foroff(res,p) { *ptr = *ptr==p?counter++:res[*ptr]; ++ptr; }
      return res;
    }