function histog = hod(x, y, z, config)   
    %computes HOD for one joints. x, y, and z are the locations of this
    %joint in nFrames. nBins is the length of the histogram
    nFrames = length(x);
    nBins = config.nBins;
    allDegrees = 360;
    histog = zeros(1, 3 * nBins);
    
    centers = ((0:nBins-1) * (allDegrees / nBins)) + (allDegrees / (2 * nBins));
    
    for j=1:nFrames-1
        % First, we make sure the frame is not corrupted
        if (x(j+1) ~= 0 || y(j+1) ~= 0 || z(j+1) ~= 0) && ...
                (x(j) ~= 0 || y(j) ~= 0 || z(j) ~= 0)
%             && floor(conf(j+1) + 0.1) && floor(conf(j) + 0.1)
            dx = x(j+1) - x(j);
            dy = y(j+1) - y(j);
            dz = z(j+1) - z(j);
            
            angxy = getangle(dx, dy);
            addVote(0, angxy, getIndex(angxy), sqrt(dx^2 + dy^2));
            
            angzy = getangle(dz, dy);
            addVote(1, angzy, getIndex(angzy), sqrt(dz^2 + dy^2));
            
            angzx = getangle(dz, dx);
            addVote(2, angzx, getIndex(angzx), sqrt(dz^2 + dx^2));
        end
    end
    function addVote(projection, angle, index, value)
        if config.softBinning
            p0 = centers(r(index));
            p1 = centers(r(index - 1));
            p2 = centers(r(index + 1));
            
            a = abs(angle - p0);
            b = abs(angle - p1);
            c = abs(angle - p2);

            nom = a * b * c / (a + b + c);

            histog(projection * nBins + index) = histog(projection * nBins + index) + ((nom / a) * value);
            
            rind = r(index - 1);
            histog(projection * nBins + rind) = histog(projection * nBins + rind) + ((nom / b) * value);
            
            rind = r(index + 1);
            histog(projection * nBins + rind) = histog(projection * nBins + rind) + ((nom / c) * value);
        else
            histog(projection * nBins + index) = histog(projection * nBins + index) + value;
        end
    end
    function ind = getIndex(angle)
        ind = floor(angle*nBins/allDegrees) + 1;
    end
    function ri = r(ind)
        if ind > nBins
            ri = ind - nBins;
        elseif ind < 1
            ri = ind + nBins;
        else
            ri = ind;    
        end 
    end
end