for (int pi = 0; pi < _px.size(); pi++)
    {
        if (_weights[pi] == 0.f) continue;

        //pos from array
        float px = _px.at(pi);
        float py = _py.at(pi);
        float pz = _pz.at(pi);

        //index of cell particle belongs to
        int i = (_px.at(pi)-_offset[0])/_h[0];
        int j = (_py.at(pi)-_offset[1])/_h[1];
        int k = (_pz.at(pi)-_offset[2])/_h[2];

        //mass to distribute
        float p_mass = _h_part_density[pi];
        float weight = _weights[pi];

        //bounds
        int low_i = std::max(0,i-1);
        int low_j = std::max(0,j-1);
        int low_k = std::max(0,k-1);

        int high_i = std::min(i+1,_nx-1);
        int high_j = std::min(j+1,_ny-1);
        int high_k = std::min(k+1,_nz-1);

        //+1/-1 support
        for (int ii = low_i; ii <= high_i; ii++)
            for (int jj = low_j; jj <= high_j; jj++)
                for (int kk = low_k; kk <= high_k; kk++)
                {
                    //target cell center
                    float posii = ii * _h[0] + _offset[0] + 0.5f*_h[0];
                    float posjj = jj * _h[1] + _offset[1] + 0.5f*_h[1];
                    float poskk = kk * _h[2] + _offset[2] + 0.5f*_h[2];

                    float mass = p_mass * tent(px-posii,_h[0])*tent(py-posjj,_h[1])*tent(pz-poskk,_h[2]) / weight;

                    //sum mass to target particle
                    _h_density->at(ii,jj,kk) += mass;
                }

    }