% Data

% Blocks - point at origin is implicitly included
block([p(1,0,0),p(0,1,0),p(0,0,1)]).
block([p(0,1,0),p(1,1,0),p(1,2,0)]).
block([p(0,1,0),p(1,1,0),p(1,1,1)]).
block([p(0,1,0),p(0,1,1),p(1,1,1)]).
block([p(1,0,0),p(1,1,0),p(2,0,0)]).
block([p(0,1,0),p(1,0,0),p(2,0,0)]).
block([p(1,0,0),p(0,1,0)]).

% 24 posssible permutation of the axes
direction(v(1,2,3)).
direction(v(1,-2,-3)).
direction(v(1,3,-2)).
direction(v(1,-3,2)).
direction(v(-1,2,-3)).
direction(v(-1,-2,3)).
direction(v(-1,3,2)).
direction(v(-1,-3,-2)).
direction(v(2,1,-3)).
direction(v(2,-1,3)).
direction(v(2,3,1)).
direction(v(2,-3,-1)).
direction(v(-2,1,3)).
direction(v(-2,-1,-3)).
direction(v(-2,3,-1)).
direction(v(-2,-3,1)).
direction(v(3,1,2)).
direction(v(3,-1,-2)).
direction(v(3,2,-1)).
direction(v(3,-2,1)).
direction(v(-3,1,-2)).
direction(v(-3,-1,2)).
direction(v(-3,2,1)).
direction(v(-3,-2,-1)).


% Generating rotations

select_axis(1, p(X,_,_), X).
select_axis(-1, p(X,_,_), X1) :- X1 is -X.
select_axis(2, p(_,Y,_), Y).
select_axis(-2, p(_,Y,_), Y1) :- Y1 is -Y.
select_axis(3, p(_,_,Z), Z).
select_axis(-3, p(_,_,Z), Z1) :- Z1 is -Z.

permute_axes(v(I,J,K), Point, p(X,Y,Z)) :-
    select_axis(I, Point, X),
    select_axis(J, Point, Y),
    select_axis(K, Point, Z).

rotation(Block, Rotated) :-
    direction(D),
    maplist(permute_axes(D), [p(0,0,0)|Block], Rotated).

minimum([Point], Point).
minimum([p(X,Y,Z)|Points], p(X1,Y1,Z1)) :-
    minimum(Points, p(X0,Y0,Z0)),
    X1 is min(X, X0),
    Y1 is min(Y, Y0),
    Z1 is min(Z, Z0).

subtract(p(X1,Y1,Z1), p(X2,Y2,Z2), p(X,Y,Z)) :-
    X is X2 - X1,
    Y is Y2 - Y1,
    Z is Z2 - Z1.

normalize(Block, Norm) :-
    minimum(Block, Min),
    maplist(subtract(Min), Block, Translated),
    sort(Translated, Norm).

all_rotations(Block, Rotations) :-
    findall(Norm, (rotation(Block, Rot), normalize(Rot, Norm)), All),
    sort(All, Rotations).


% Model loader

read_points(0, _, []).
read_points(N, Fd, [p(X,Y,Z)|Model]) :-
    N > 0, N1 is N - 1,
    read(Fd, (X, Y, Z)),
    read_points(N1, Fd, Model).

load_model(Filename, Model) :-
    open(Filename, read, Fd, []),
    read_points(27, Fd, Model),
    close(Fd).


% Solver

solve(Model, Solution) :-
    findall(B, block(B), Blocks),
    solve(Blocks, Model, Solution).

solve([], [], []).
solve([Block|Rest], Model, [Transformed|Solution]) :-
    all_rotations(Block, Rotations),
    member(Rot, Rotations),
    translate(Rot, Model, Transformed),
    remove(Transformed, Model, Model1),
    solve(Rest, Model1, Solution).

translate([Pos|Rest], Model, [Pos1|Trans]) :-
    member(Pos1, Model),
    subtract(Pos, Pos1, Delta),
    maplist(add(Delta), Rest, Trans).

add(p(X1,Y1,Z1), p(X2,Y2,Z2), p(X,Y,Z)) :-
    X is X2 + X1,
    Y is Y2 + Y1,
    Z is Z2 + Z1.

remove([], Y, Y).
remove([X|Xs], Y, Z) :- select(X, Y, Y1), remove(Xs, Y1, Z).


% Output

writeVTK(Solution, Filename) :-
    open(Filename, write, Fd),
    write(Fd, '# vtk DataFile Version 2.0'), nl(Fd),
    write(Fd, 'Block by Block solution'), nl(Fd),
    write(Fd, 'ASCII'), nl(Fd),
    write(Fd, 'DATASET POLYDATA'), nl(Fd),
    write(Fd, 'POINTS 27 int'), nl(Fd),
    write_points(Fd, Solution),
    write(Fd, 'POINT_DATA 27'), nl(Fd),
    write(Fd, 'SCALARS group int 1'), nl(Fd),
    write(Fd, 'LOOKUP_TABLE default'), nl(Fd),
    write(Fd, '1 1 1 1 2 2 2 2 3 3 3 3 4 4 4 4 5 5 5 5 6 6 6 6 7 7 7'), nl(Fd),
    close(Fd).

write_points(_, []).
write_points(Fd, [[]|Rest]) :- write_points(Fd, Rest).
write_points(Fd, [[p(X,Y,Z)|Points]|Rest]) :-
    write(Fd, X), write(Fd, ' '),
    write(Fd, Y), write(Fd, ' '),
    write(Fd, Z), nl(Fd),
    write_points(Fd, [Points|Rest]).