function heatflux_vectorfield(image_data)

    % first step is to flatten it to 1 dimension
    [X,map] = rgb2ind(image_data, 65536);
    T = ind2gray(X,map);
    %T = mean(image_data, 3);
    flat_size = size(T);
    width = flat_size(1);
    height = flat_size(2);
    
    % second step is to calculate a vector field
    % the vector field is comprised of U, V which 
    % U = X axis run, V = Y axis rise.  
    u = zeros(width,height);
    v = zeros(width,height);
    
    for i=2:(width-1)
        for j=2:(height-1)
            kp = 0.49;
        
            qx = -kp*( cast(T(i+1,j),'double') - cast(T(i-1,j), 'double'))/2 ;
            qy = -kp*( cast(T(i,j+1),'double') - cast(T(i,j-1), 'double'))/2;
            qn = sqrt(power(qx,2) + power(qy,2));
            qtheta = atan(qy/qx);
        
            if qx < 0
                qtheta = qtheta + pi;
            end
            
            u(width - i,j) = qn*cos(qtheta);
            v(width - i,j) = qn*sin(qtheta);
        end
    end
    
    % Display the vector field.
    [x, y] = meshgrid(1:1:height, 1:1:width);
    quiver(x,y,u,v);
end