function [IGrad, tSol, k] = ...
    poissonSolverJacobi (IGrad, Mask, B, maxIter, tau)
%
% poissonSolverJacobi -- Jacobi iterative solver for poisson equation
%
% SYNTAX
%
%   [IGrad, tSol, k] = poissonSolverJacobi (IGrad, Mask, B, nMaxIter, tau)
%
% INPUT
%
%   IGrad       Original gradient field                 [M-by-N]
%   Mask        Mask specifies the solving domain       [M-by-N, logical]
%   B           Guidance gradient field                 [M-by-N]
%   maxIter     The maximum iteration of the solver     [integer]
%   tau         The error tolerance level               [scalar]
%
% OUTPUT
%
%   IGrad       Output gradient field                   [M-by-N]
%   tSol        Time for solving the poisson equation   [scalar]
%   k           The maximum iteration number reached    [scalar]
%
% DESCRIPTION
%
%   [IGrad, tSol, k] = poissonSolverJacobi (IGrad, Mask, B, nMaxIter, tau)
%   solves the poisson equation using Jacobi iteration.
%
% NOTES
%
%   Simple and naive Jacobi solver.
%
% REFERENCE(S)
%
%   P?rez, Patrick, Michel Gangnet, and Andrew Blake.
%   "Poisson image editing." ACM Transactions on Graphics (TOG).
%   Vol. 22. No. 3. ACM, 2003.
%
%   Solving the Discrete Poisson Equation using Jacobi, SOR,
%   Conjugate Gradients, and the FFT
%   https://people.eecs.berkeley.edu/~demmel/cs267/lecture24/lecture24.html
%
%
    tic;
    %% Masked domain


    Mask1 = double((Mask));

    mIdx = find(Mask1 == 1);

    Mask2 = mod(mIdx, 2);
    odds = mIdx(Mask2== 1);
    evens = mIdx(Mask2== 0);

    %odd
    %even

    %Find the black nodes
    %Find the red nodes

    fprintf('starting');
    w = 1.97;
    %% Jacobi Iteration
    for k = 1 : maxIter


        % previous-step solution
        IGradPrev = IGrad;

        % Update step
        IGrad(odds) = IGradPrev(odds) + w * ( IGrad(odds+1) + IGrad(odds-1) + ...
            IGrad(odds+size(IGrad, 1)) + ...
            IGrad(odds-size(IGrad, 1)) - B(odds) - 4*IGradPrev(odds) ) / 4;

         IGrad(evens) = IGradPrev(evens)+ w*( IGrad(evens+1) + IGrad(evens-1) + ...
            IGrad(evens+size(IGrad, 1)) + ...
            IGrad(evens-size(IGrad, 1)) - B(evens) - 4*IGradPrev(evens) ) / 4;

        % termination condition (convergence; solution stagnation)
        if (sqrt(sum((IGrad(mIdx) - IGradPrev(mIdx)).^2))) < tau
            break
        end

    end
    tSol = toc;

end


%%------------------------------------------------------------
%
% AUTHORS
%
%   Tiancheng Liu       tcliu@cs.duke.edu
%
% REVISIONS
%
%   0.1 (Spring 2017) - Tiancheng
%
% ------------------------------------------------------------