CImg<T>& _fill(const char *const expression, const bool repeat_values, const unsigned int formula_mode,
                   CImgList<T> *const list_images, const char *const calling_function, const CImg<T> *provides_copy) {
      if (is_empty() || !expression || !*expression) return *this;
      const unsigned int omode = cimg::exception_mode();
      cimg::exception_mode(0);
      CImg<charT> is_error;
      bool is_value_sequence = false;
      cimg_abort_init;

      if (formula_mode) {

        // Try to pre-detect regular value sequence to avoid exception thrown by _cimg_math_parser.
        double value;
        char sep;
        const int err = cimg_sscanf(expression,"%lf %c",&value,&sep);
        if (err==1 || (err==2 && sep==',')) {
          if (err==1) { if (formula_mode==2) return *this; return fill((T)value); }
          else is_value_sequence = true;
        }

        // Try to fill values according to a formula.
        _cimg_abort_init_openmp;
        if (!is_value_sequence) try {
            CImg<T> base = provides_copy?provides_copy->get_shared():get_shared();
            _cimg_math_parser mp(expression + (*expression=='>' || *expression=='<' ||
                                               *expression=='*' || *expression==':'),
                                 calling_function,base,this,list_images,true);
            if (!provides_copy && expression && *expression!='>' && *expression!='<' && *expression!=':' &&
                mp.need_input_copy)
              base.assign().assign(*this,false); // Needs input copy

            // Determine 2nd largest image dimension (used as axis for inner loop in parallelized evaluation).
            unsigned int M;
            if (mp.result_dim) {
              M = cimg::max(_width,_height,_depth);
              M = M==_width?std::max(_height,_depth):M==_height?std::max(_width,_depth):std::max(_width,_height);
            } else {
              M = cimg::max(_width,_height,_depth,_spectrum);
              M = M==_width?cimg::max(_height,_depth,_spectrum):
                M==_height?cimg::max(_width,_depth,_spectrum):
                M==_depth?cimg::max(_width,_height,_spectrum):cimg::max(_width,_height,_depth);
            }

            bool do_in_parallel = false;
#if cimg_use_openmp!=0
            cimg_openmp_if(*expression=='*' || *expression==':' ||
                           (mp.is_parallelizable && M>=(cimg_openmp_sizefactor)*320 && size()/M>=2))
              do_in_parallel = true;
#endif
            if (mp.result_dim) { // Vector-valued expression
              const unsigned int N = std::min(mp.result_dim,_spectrum);
              const ulongT whd = (ulongT)_width*_height*_depth;
              T *ptrd = *expression=='<'?_data + _width*_height*_depth - 1:_data;

              if (*expression=='<') {
                CImg<doubleT> res(1,mp.result_dim);
                mp.begin_t();
                cimg_rofYZ(*this,y,z) {
                  cimg_abort_test;
                  if (formula_mode==2) cimg_rofX(*this,x) mp(x,y,z,0);
                  else cimg_rofX(*this,x) {
                      mp(x,y,z,0,res._data);
                      const double *ptrs = res._data;
                      T *_ptrd = ptrd--; for (unsigned int n = N; n>0; --n) { *_ptrd = (T)(*ptrs++); _ptrd+=whd; }
                    }
                }
                mp.end_t();

              } else if (*expression=='>' || !do_in_parallel) {
                CImg<doubleT> res(1,mp.result_dim);
                mp.begin_t();
                cimg_forYZ(*this,y,z) {
                  cimg_abort_test;
                  if (formula_mode==2) cimg_forX(*this,x) mp(x,y,z,0);
                  else cimg_forX(*this,x) {
                      mp(x,y,z,0,res._data);
                      const double *ptrs = res._data;
                      T *_ptrd = ptrd++; for (unsigned int n = N; n>0; --n) { *_ptrd = (T)(*ptrs++); _ptrd+=whd; }
                    }
                }
                mp.end_t();

              } else {

#if cimg_use_openmp!=0
                cimg_pragma_openmp(parallel)
                {
                  _cimg_math_parser
                    *const _mp = omp_get_thread_num()?new _cimg_math_parser(mp):&mp,
                    &lmp = *_mp;
                  lmp.is_fill = true;
                  cimg_pragma_openmp(barrier)
                  lmp.begin_t();

#define _cimg_fill_openmp_vector(_YZ,_y,_z,_X,_x,_sx,_sy,_sz,_off) \
  cimg_pragma_openmp(for cimg_openmp_collapse(2)) \
  cimg_for##_YZ(*this,_y,_z) _cimg_abort_try_openmp { \
    cimg_abort_test; \
    if (formula_mode==2) cimg_for##_X(*this,_x) lmp(x,y,z,0); \
    else { \
      CImg<doubleT> res(1,lmp.result_dim); \
      T *__ptrd = data(_sx,_sy,_sz,0); \
      const ulongT off = (ulongT)_off; \
      cimg_for##_X(*this,_x) { \
        lmp(x,y,z,0,res._data); \
        const double *ptrs = res._data; \
        T *_ptrd = __ptrd; \
        for (unsigned int n = N; n>0; --n) { *_ptrd = (T)(*ptrs++); _ptrd+=whd; } \
        __ptrd+=off; \
      } \
    } \
  } _cimg_abort_catch_openmp _cimg_abort_catch_fill_openmp

                  if (M==_width) { _cimg_fill_openmp_vector(YZ,y,z,X,x,0,y,z,1) }
                  else if (M==_height) { _cimg_fill_openmp_vector(XZ,x,z,Y,y,x,0,z,_width) }
                  else { _cimg_fill_openmp_vector(XY,x,y,Z,z,x,y,0,_width*_height) }

                  lmp.end_t();
                  cimg_pragma_openmp(barrier) cimg_pragma_openmp(critical) { lmp.merge(mp); }
                  if (&lmp!=&mp) delete &lmp;
                }
#endif
              }

            } else { // Scalar-valued expression
              T *ptrd = *expression=='<'?end() - 1:_data;
              if (*expression=='<') {
                mp.begin_t();
                if (formula_mode==2) cimg_rofYZC(*this,y,z,c) { cimg_abort_test; cimg_rofX(*this,x) mp(x,y,z,c); }
                else cimg_rofYZC(*this,y,z,c) { cimg_abort_test; cimg_rofX(*this,x) *(ptrd--) = (T)mp(x,y,z,c); }
                mp.end_t();

              } else if (*expression=='>' || !do_in_parallel) {
                mp.begin_t();
                if (formula_mode==2) cimg_forYZC(*this,y,z,c) { cimg_abort_test; cimg_forX(*this,x) mp(x,y,z,c); }
                else cimg_forYZC(*this,y,z,c) { cimg_abort_test; cimg_forX(*this,x) *(ptrd++) = (T)mp(x,y,z,c); }
                mp.end_t();

              } else {

#if cimg_use_openmp!=0
                cimg_pragma_openmp(parallel)
                {
                  _cimg_math_parser
                    *const _mp = omp_get_thread_num()?new _cimg_math_parser(mp):&mp,
                    &lmp = *_mp;
                  lmp.is_fill = true;
                  cimg_pragma_openmp(barrier)
                  lmp.begin_t();

#define _cimg_fill_openmp_scalar(_YZC,_y,_z,_c,_X,_x,_sx,_sy,_sz,_sc,_off) \
  cimg_pragma_openmp(for cimg_openmp_collapse(3)) \
  cimg_for##_YZC(*this,_y,_z,_c) _cimg_abort_try_openmp { \
    cimg_abort_test; \
    if (formula_mode==2) cimg_for##_X(*this,_x) lmp(x,y,z,c); \
    else { \
      T *_ptrd = data(_sx,_sy,_sz,_sc); \
      const ulongT off = (ulongT)_off; \
      cimg_for##_X(*this,_x) { *_ptrd = (T)lmp(x,y,z,c); _ptrd+=off; } \
    } \
  } _cimg_abort_catch_openmp _cimg_abort_catch_fill_openmp

                  if (M==_width) { _cimg_fill_openmp_scalar(YZC,y,z,c,X,x,0,y,z,c,1) }
                  else if (M==_height) { _cimg_fill_openmp_scalar(XZC,x,z,c,Y,y,x,0,z,c,_width) }
                  else if (M==_depth) { _cimg_fill_openmp_scalar(XYC,x,y,c,Z,z,x,y,0,c,_width*_height) }
                  else { _cimg_fill_openmp_scalar(XYZ,x,y,z,C,c,x,y,z,0,_width*_height*_depth) }

                  lmp.end_t();
                  cimg_pragma_openmp(barrier) cimg_pragma_openmp(critical) { lmp.merge(mp); }
                  if (&lmp!=&mp) delete &lmp;
                }
#endif
              }
            }
            mp.end();
          } catch (CImgException& e) { CImg<charT>::string(e._message).move_to(is_error); }
      }

      // Try to fill values according to a value sequence.
      if (!formula_mode || is_value_sequence || is_error) {
        CImg<charT> item(256);
        char sep = 0;
        const char *nexpression = expression;
        ulongT nb = 0;
        const ulongT siz = size();
        T *ptrd = _data;
        for (double val = 0; *nexpression && nb<siz; ++nb) {
          sep = 0;
          const int err = cimg_sscanf(nexpression,"%255[ \n\t0-9.eEinfa+-]%c",item._data,&sep);
          if (err>0 && cimg_sscanf(item,"%lf",&val)==1 && (sep==',' || sep==';' || err==1)) {
            nexpression+=std::strlen(item) + (err>1);
            *(ptrd++) = (T)val;
          } else break;
        }
        cimg::exception_mode(omode);
        if (nb<siz && (sep || *nexpression)) {
          if (is_error) throw CImgArgumentException("%s",is_error._data);
          else throw CImgArgumentException(_cimg_instance
                                           "%s(): Invalid sequence of filling values '%s'.",
                                           cimg_instance,calling_function,expression);
        }
        if (repeat_values && nb && nb<siz)
          for (T *ptrs = _data, *const ptre = _data + siz; ptrd<ptre; ++ptrs) *(ptrd++) = *ptrs;
      }

      cimg::exception_mode(omode);
      cimg_abort_test;
      return *this;
    }