amazer.m

function [solution,solution_image] = amazer(mz,p,q)
% Solve mazes using image processing!
% mz is a binary image of a maze (1s represent walls)
% p and q are the start and end points

% walls = bwmorph(mz,'dilate',2);
% walls = bwmorph(mz,'thicken',5);
walls = mz;

DA = bwdistgeodesic(~walls,p(1),p(2),'quasi-euclidean');
DB = bwdistgeodesic(~walls,q(1),q(2),'quasi-euclidean');

D = DA + DB;
D = round(D * 8) / 8;

D(isnan(D)) = inf;

paths = imregionalmin(D);
solution = bwmorph(paths,'thin',Inf);
thick_solution = bwmorph(solution,'dilate');

solution_image = imoverlay(mz,thick_solution,[1 0 0]);