    Tfloat _linear_atXYZ_p(const float fx, const float fy=0, const float fz=0, const int c=0) const {
      const float
        nfx = cimg::mod(fx,_width - 0.5f),
        nfy = cimg::mod(fy,_height - 0.5f),
        nfz = cimg::mod(fz,_depth - 0.5f);
      const unsigned int
        x = (unsigned int)nfx,
        y = (unsigned int)nfy,
        z = (unsigned int)nfz;
      const float
        dx = nfx - x,
        dy = nfy - y,
        dz = nfz - z;
      const unsigned int
        nx = cimg::mod(x + 1,_width),
        ny = cimg::mod(y + 1,_height),
        nz = cimg::mod(z + 1,_depth);
      const Tfloat
        Iccc = (Tfloat)(*this)(x,y,z,c), Incc = (Tfloat)(*this)(nx,y,z,c),
        Icnc = (Tfloat)(*this)(x,ny,z,c), Innc = (Tfloat)(*this)(nx,ny,z,c),
        Iccn = (Tfloat)(*this)(x,y,nz,c), Incn = (Tfloat)(*this)(nx,y,nz,c),
        Icnn = (Tfloat)(*this)(x,ny,nz,c), Innn = (Tfloat)(*this)(nx,ny,nz,c);
      return Iccc +
        (Incc - Iccc +
         (Iccc + Innc - Icnc - Incc +
          (Iccn + Innn + Icnc + Incc - Icnn - Incn - Iccc - Innc)*dz)*dy +
         (Iccc + Incn - Iccn - Incc)*dz)*dx +
        (Icnc - Iccc +
         (Iccc + Icnn - Iccn - Icnc)*dz)*dy +
        (Iccn - Iccc)*dz;
    }
    Tfloat _linear_atXYZ_p(const float fx, const float fy=0, const float fz=0, const int c=0) const {
      const float
        nfx = cimg::mod(fx,_width - 0.5f),
        nfy = cimg::mod(fy,_height - 0.5f),
        nfz = cimg::mod(fz,_depth - 0.5f);
      const unsigned int
        x = (unsigned int)nfx,
        y = (unsigned int)nfy,
        z = (unsigned int)nfz;
      const float
        dx = nfx - x,
        dy = nfy - y,
        dz = nfz - z;
      const unsigned int
        nx = cimg::mod(x + 1,_width),
        ny = cimg::mod(y + 1,_height),
        nz = cimg::mod(z + 1,_depth);
      const Tfloat
        Iccc = (Tfloat)(*this)(x,y,z,c), Incc = (Tfloat)(*this)(nx,y,z,c),
        Icnc = (Tfloat)(*this)(x,ny,z,c), Innc = (Tfloat)(*this)(nx,ny,z,c),
        Iccn = (Tfloat)(*this)(x,y,nz,c), Incn = (Tfloat)(*this)(nx,y,nz,c),
        Icnn = (Tfloat)(*this)(x,ny,nz,c), Innn = (Tfloat)(*this)(nx,ny,nz,c);
      return Iccc +
        (Incc - Iccc +
         (Iccc + Innc - Icnc - Incc +
          (Iccn + Innn + Icnc + Incc - Icnn - Incn - Iccc - Innc)*dz)*dy +
         (Iccc + Incn - Iccn - Incc)*dz)*dx +
        (Icnc - Iccc +
         (Iccc + Icnn - Iccn - Icnc)*dz)*dy +
        (Iccn - Iccc)*dz;
    }