Sign up ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Q1: I have a page of a printed text file, but it was cut by a paper shredder (from top to bottom). You can download all the paper fragments here. Here are some of them. How can I reconstruct it?

shredded text

Q2: This time, the situation is even worse - the page was cut by a shredder twice (from top to bottom and from left to right)! (Get the fragments here).

really shredded text

My code for Q1:

filenames = FileNames["*.bmp", "C:\\Users\\Liao\\Desktop\\B\\file1"];
pics = Import /@ filenames;
binarypic = Binarize[#, 0.9] & /@ pics;
binarydata = ImageData /@ binarypic;
leftandrightdata = 
MapIndexed[{#2[[1]], #1[[All, 1]], #1[[All, -1]]} &, binarydata];
result = {};
firstPic = Select[leftandrightdata, ! MemberQ[#[[2]], 0] &][[1, 1]];
start = leftandrightdata[[firstPic]];
DeleteCases[leftandrightdata, start] //. 
      x_ /; Length[x] > 0 :> (sort = 
   SortBy[x, HammingDistance[#[[2]], start[[3]]] &];
     AppendTo[result, start = sort[[1]]]; Rest[sort]);
Prepend[List /@ result[[All, 1]], {firstPic}]] // ImageAssemble

This code maybe dizzy for you, but it doesn't matter. The core idea is simple and straightforward, which use HammingDistance to determine whether two images can be combined together.

But for Q2, I can't really find a good way to solve it.

share|improve this question
how did you manage to put that stripes to the shredder? :p – Kuba Sep 19 '13 at 4:58
@Kuba It's a tough work! Well, actually it's a problem from a Mathematical Contest in Modeling. Now the contest have finished, so I post it here, looking forward to funny and exellent anwsers. – mm.Jang Sep 19 '13 at 5:27
You can try ImmageAssemble at each permutation, then TextRecognize -> DictionaryLookup -> Length. I think it would be fun-compact-neverending :) – Kuba Sep 19 '13 at 10:13
@Kuba It's too horrible to imagine what will happen :p – mm.Jang Sep 19 '13 at 14:54
Can you upload the files in a non .rar format? Perhaps a .zip or .tar.gz? – R. M. Sep 20 '13 at 3:17

3 Answers 3


pics = Binarize /@ (Import /@ FileNames["*.bmp"]);

(*extracting edges*)
left = ImageData[#][[All, 1]] & /@ pics; 
right = ImageData[#][[All, -1]] & /@ pics;

Let's find the left edge, here I'm assuming that it is all white. So if

pics // First // ImageDimensions
{72, 1980}

and an image is binarized then the left edge is:

Total /@ left // Position[#, 1980, 1] &

Instead of HammingDistance I'm calculating Norm:

cons = {4};
   With[{dists = Norm[right[[cons[[-1]]]] - #] & /@ left}, 
        AppendTo[cons, Position[dists, Min@dists, 1][[1, 1]] ]
  , {Length@pics - 1}]


enter image description here

share|improve this answer
Scroll down, my answer only solves simple case. – Kuba Apr 30 at 5:25

Greedy algorithm (1D & 2D)

Data loading:

data = N@ImageData@Import[#] & /@ FileNames["*.bmp"] // Developer`ToPackedArray;
n = Length[data];

The matching of two edges is defined by the property

$$ \frac{\mathop{\rm Var} [e_1+e_2]}{\mathop{\rm Var} [e_1-e_2]} $$

where $\mathop{\rm Var}[X]$ is the variance of $X$

r = 7;
var[list_] := 
  Total[#, {3, 4}] &[(Transpose[list, {1, 2, 4, 3}] - 
     Mean@Transpose[list, {2, 3, 1, 4}])^2];
corr[list1_, list2_] := 
  var@GaussianFilter[list1 + list2, {{0, 0, r, 0}}]/(0.0001 + 
     var@GaussianFilter[list1 - list2, {{0, 0, r, 0}}]);

Here r=7 is the radius of smoothing. It is very helpfull for real images (see below). 0.0001 is a small number to prevent division by zero.

The comparison each tile with each other:

top = corr[ConstantArray[data[[All, 1]], n], 
   Transpose@ConstantArray[data[[All, -1]], n]];
bottom = Transpose[top];
left = corr[ConstantArray[data[[All, All, 1]], n], 
   Transpose@ConstantArray[data[[All, All, -1]], n]];
right = Transpose[left];

Tiles in the resulting mosaic and possible positions of the next tiles:

tiles = SparseArray[{}, {2 n, 2 n}];
next = {};

The addition of the tile num to the position {i,j} of the resulting mosaic and the update of next tiles:

add[num_, {i_, j_}] := (next = DeleteCases[next, {i, j}];
  tiles[[i, j]] = num;
  left[[num, All]] = right[[num, All]] = top[[num, All]] = bottom[[num, All]] = 0.0; 
  If[tiles[[##]] == 0 && Not@MatchQ[next, {##}], AppendTo[next, {##}];] & 
      @@@ {{i, j + 1}, {i, j - 1}, {i + 1, j}, {i - 1, j}})

Quality of the matching of all possible tiles at the position {i,j} is just the sum of the matching with neighbors:

quality[i_, j_] := 
  If[#2 != 0, #[[All, #2]], 0 #[[All, 1]]] & @@@ {{left, 
      tiles[[i, j + 1]]}, {right, tiles[[i, j - 1]]}, {top, 
      tiles[[i + 1, j]]}, {bottom, tiles[[i - 1, j]]}} // Total;

Addition of two tiles with the best matching:

(add[#1, {n, n}]; add[#2, {n, n + 1}]) & @@@ Position[left, Max[left]];

Addition of n-2 remaining tiles to the best positions:

Do[add[#2, next[[#]]] & @@ First@Position[#, Max[#]] &[quality @@@ next], {n - 2}];

The result for Q1:

view[margins__] := 
 Image@Flatten[#, {{1, 3}, {2, 4}}] &@
  Map[If[# <= 0, 1 + 0 data[[1]], data[[#]]] &, 
   tiles[[#1 ;; #2, #3 ;; #4]] &[margins], {2}]

view[Min[#1], Max[#1], Min[#2], Max[#2]] & @@ Transpose@tiles["NonzeroPositions"]

enter image description here

Unfortunately Q2 have too many blank edges to reconstruction so I use Lena photo

data = N@ImageData /@ #[[PermutationList@RandomPermutation@Length[#]]] &@
  Flatten@ImagePartition[ExampleData[{"TestImage", "Lena"}], 64];

After repeating the steps above we got the result

enter image description here

Without smoothing of edges (see r above) the matching fails because of noise and fine structure (especially in the hat).

share|improve this answer
Cool! But my bounty has expired......:-/ – mm.Jang Oct 10 '13 at 6:11
@mm.Jang If the bounty is expired the reputation is losted? I have never started a bounty. – ybeltukov Oct 10 '13 at 19:21
Very nice! My suggestion for the blank edges may be to give them huge (unnatural) weights, that way it avoids connecting any white edges, if there is no way to move across a image with blank edges (in the case of the letter) it may be a natural gap. – lalmei Nov 1 '13 at 16:13
I'm trying to understand your weights between edges better (specifically how those transpose's are acting on the image data), since I'm not very familiar with image analysis. Are you taking a variance over each color channel ? then summing ? – lalmei Nov 4 '13 at 17:15
@lalmei Yes, I sum squares of every edge pixel and each color channel (Total with parameter {3,4}). Transpose is a bit complicated, but it helps me to subtract the mean value (also with each channel separately). – ybeltukov Nov 4 '13 at 17:23

Q2: This answer has been sitting in my trunk of files for a year and a half now. I was hoping to improve it and use it to solve the original problem but never got around to it, at the moment it solves images/photos with noisy data. Might as well post it before I forget how it even works.

So I wanted to solve this problem with a Minimum Spanning Tree Algorithm. I decided to use Kruskal's algorithm, the idea being that people solve a jigsaw puzzle as they find pieces that match (maybe there is a way to parallelize this). This way you build a Forest of Trees made up of the best matches. So you can stop at any moment, and in the worst case you have a bunch of larger matched pieces.

At the time my idea was to use this algorithm for the main questions and get it to a point where matching using text recognition wasn't too time consuming since now you would have less pieces with larger pieces of words on them (make a greedy algorithm that matched the pieces to maximize the number of english words it found). But even then using text recognition at the time took a very long time.

The main problem I came into was trying to make sure that when you merge two Trees the images weren't overlapping. I made a modification to make it work, I call this the Minimum Spaning Geometrical Tree. There are two functions, one that gives the relative positions of pieces in a tree Chargeds. And one that checks if there is a path to a specific relative position for a piece, ChargedPath. Here I call the positions charges, since I envisioned this for a more generalized use (e.g. particle tracking. )

The main part of the code is to build an adjacency list like this $$ \{ node1, \{\{"up", node2\}, \{"down", node24\}, \{"left", node64\}, \cdots \}\}, \\ \hspace{-11.5cm} \{node2, \cdots\} $$ Except instead of "up","down","left","right", I use 1,2,3,4. (maybe not in that order)


For weights I used the same as @ybeltukov, since it works best for noisy images.

r = 12;
var2[list_] := Total@(Variance[#] & /@ (Thread[list]));
corr[list1_, list2_] := 
   var2@GaussianFilter[list1 - list2, {{0, 0, r, 0}}]/
   (var2@GaussianFilter[list1 + list2, {{0, 0, r, 0}}]);

And I organize the edge pixels as the OP.

testdata = #[[PermutationList@RandomPermutation@Length[#]]] &@
 Flatten@ImagePartition[ExampleData[{"TestImage", "Lena"}], 64];
pics = testdata;

test2 = ImageData[#] & /@ pics;
test3 = MapIndexed[{#2[[1]], #1[[All, 1]], #1[[All, -1]], 
     #1[[1, All]], #1[[-1, All]]} &, test2];

Then I create a distance matrix, it checks $n^2$ distance between the edges, for left/right and up/down. It would be better organized as a bipartite graph but so it is for now. Amazingly, this is the slowest part of the code.

n = Length[pics]
 leftrightF = 
   Table[If[i == j, 64^2, N[corr[test3[[i, 3]], test3[[j, 2]]]]], {i, 
    1, n}, {j, 1, n}];
 updownF = 
   Table[If[i == j, 64^2, N[corr[test3[[i, 5]], test3[[j, 4]]]]], {i, 
   1, n}, {j, 1, n}];

Minimum Spaning Geometrical tree

The main ingredient of the Kruskal Algorithm is that when you connect an edge you have to make you create no loops. If there is already a path between two points adding a connection between them will lead to a loop, so this function checks if there is a path between two points. I use Throw and Catch with a While loop for a depth-first search, and pack it all inside the Return. ( There might be a more efficient way to do this with Scan)

ITAPath[INadjgraph_List, point1_Integer, point2_Integer, label_: 1] :=
        Module[{adjgraph = fillin[Function[{y}, 
        {y[[1]], Function[{x}, x[[2]] ] /@ (y[[2]])}] /@ INadjgraph], 
        stack = {}, currentpos, marked},

  marked = Table[False, {i, 1, Length[adjgraph]}];
  marked[[point1]] = True;
  stack = Join[stack, adjgraph[[point1, 2]]];

           If[point1 == point2, Throw[True];];
           While[Length[stack] != 0,
            currentpos = Last@stack;
            stack = Drop[stack, {-1}];
              If[currentpos == point2, Throw[True];];
              marked[[currentpos]] = True;
              stack = Join[stack, adjgraph[[currentpos, 2]]];

This is essentially a Depth-First Search Algorithm to see if there is a path (ITAP) between two points with the adjancecy list. I use a fill function to make the edges undirected in the adjacency list. For some reason I keep the adjacency list directed. ( Actually I think this is because I think the minimum weight between two edges is not always back to the same node)

The second ingredient is to make sure the new point added to the tree does not overlap ( in 2D space) with other points already in the tree. I do this by looking up the "charge" or position of each point in the tree relative to the new point. This is done with the ChargePath function

ChargedPath[INadjgraph_List, point1_Integer, charge_] :=
   Module[{adjgraph = fillin2@INadjgraph, stack = {}, currentpos, 
           marked, chargem, chargeadd, side, parent},

   marked = Table[False, {i, 1, Length[adjgraph]}];
   parent = Table[Null, {i, 1, Length[adjgraph]}];
   marked[[point1]] = True;
   chargem = {0, 0};
   stack = Join[stack, adjgraph[[point1, 2]]];
        parent[[#[[2]]]] = {chargem, point1}; &, 
        adjgraph[[point1, 2]]];

     If[chargem == charge, Throw[True];];
     While[Length[stack] != 0,
           {side, currentpos} = Last@stack;
           stack = Drop[stack, {-1}];

        chargeadd = 
        Which[side == 1, {0, -1}, 
              side == 2, {0, 1}, 
              side == 3, {-1, 0},
              side == 4, {1, 0}];
        chargem = parent[[currentpos, 1]];
        chargem = chargeadd + chargem;
        If[chargem == charge, Throw[True];];
        marked[[currentpos]] = True;
        stack = Join[stack, adjgraph[[currentpos, 2]]];
             parent[[#[[2]]]] = {chargem, currentpos}; &, 
             adjgraph[[currentpos, 2]]];

Again it is a Depth-First Search algorithm to quickly find if there is a point with the same "charge" or position.

These are some custodial functions to filling the adjacency list with undirected edges for traversing the graph. A function to list the charges (positions) of each node in the tree.

fillin[adj_] := Module[{adjtest = adj},
    If[! MemberQ[adjtest[[x, 2]], y[[1]]],
           adjtest[[x, 2]] = Append[adjtest[[x, 2]], y[[1]]];
    ], y[[2]]];
  ], adjtest];
reve[y_] := Which[y == 1, 2, y == 2, 1, y == 3, 4, y == 4, 3];

fillin2[adj_] := Module[{adjtest = adj},
    If[! MemberQ[adjtest[[x[[2]], 2]], {reve[x[[1]]], y[[1]]}],

      adjtest[[x[[2]], 2]] = 
        Append[adjtest[[x[[2]], 2]], {reve[x[[1]]], y[[1]]}];
    ], y[[2]]];
 ], adjtest];

 Chargeds[INadjgraph_List, point1_Integer] :=
 Module[{adjgraph = fillin2@INadjgraph, stack = {}, currentpos, 
         marked, chargem, chargeadd, side, parent},

  marked = Table[False, {i, 1, Length[adjgraph]}];
  parent = Table[Null, {i, 1, Length[adjgraph]}];
  marked[[point1]] = True;

         chargem = {0, 0};
         Sow[{point1, chargem}, 1];
         stack = Join[stack, adjgraph[[point1, 2]]];

              parent[[#[[2]]]] = {chargem, point1}; &, 
              adjgraph[[point1, 2]]];

     While[Length[stack] != 0,
           {side, currentpos} = Last@stack;
           stack = Drop[stack, {-1}];
           If[! marked[[currentpos]],
           chargeadd = 
               Which[side == 1, {-1, 0}, 
                     side == 2, {1, 0}, 
                     side == 3, {0, 1},
                     side == 4, {0, -1}];

      chargem = parent[[currentpos, 1]];
      chargem = chargeadd + chargem;

      Sow[{currentpos, chargem}, 1];

      marked[[currentpos]] = True;
      stack = Join[stack, adjgraph[[currentpos, 2]]];
          parent[[#[[2]]]] = {chargem, currentpos}; &, 
          adjgraph[[currentpos, 2]]];

Building the Adjacency List

Update: simplified the code significantly

Building the spanning tree is the following code.

The main if statements are to separate which labels to give the new edges, since we are looking at global minima of edges(up down left right). Once we go to the correct label that's where we check for "no loops" (kruskal algorithm for mst) and the fact that the tree's don't overlap with due to the geometry (the new part of the algorithm), on top of a bunch of checks to make sure we only have 4 edges and that edges has not been connected already to some other part of the tree.

You can also set a max weight, so that if the weight gets larger than a certain value, then the algorithm stops and just gives the remaining Forest, which you could at that point apply some additional weights to how to match up those remaining Trees.

MSGT[leftright_, updown_, n_, minset_: 50., LRsidesize_: 180,UDsidesize_: 72] := 
   Module[{leftrightweights = leftright, updownweights = updown,
        Adjs = Table[{i, {}}, {i, 1, n}], minlr, minud, pos, 
        location2, location1, i, k, directions, labels, 
        Adjs2, oldimage, newimage, newcount, oldcount, t},

     (*Since the result is a tree, I can only add n-1 edges*)       

     k = n - 1;

     While[k > 0,

     (*find global minimum*) 

        minlr = Min[leftrightweights];
        minud = Min[updownweights];

     (* check if the minimum is an up down edge or left right *)
      (* Once I pass by a minimum I replace it with a large weight, It might be easier or faster to just keep track of *)
        If[minlr < minud,
        pos = Position[leftrightweights, minlr];
        leftrightweights [[First@(First@pos), Last@(First@pos)]] =LRsidesize^2;
        directions = {-1, 0};
        labels = {1, 2};
        pos = Position[updownweights, minud];
        updownweights [[First@(First@pos), Last@(First@pos)]] =UDsidesize^2;
        directions = {0, 1};
        labels = {3, 4};

         location1 = pos[[1, 1]];
         location2 =  pos[[1, 2]];

      If[location1 < location2, 
         {location1, location2} = {location1, location2} /. {location1 ->location2,location2 -> location1};
         labels = RotateLeft[labels];
         directions = -directions;

       (* if the weights are less the pre-set minimum quit *)
       If[(minlr > minset && minud > minset), Break[];];

 (* Main meat of the algorithm *)    
(* I group the checks in two sets *)   
(* first set is to check if the adjacency list is being build correctly, i.e. is it full ? does another node points to it from the opposite label etc. *)

If[(MemberQ[Adjs[[location1, 2]], {Last@labels, _}]), Continue[];];
If[(MemberQ[Adjs, {Last@labels, location2}, {3}]), Continue[];];  
If[(MemberQ[Adjs, {First@labels, location1}, {3}]), Continue[];];
If[(MemberQ[Adjs[[location2, 2]], {First@labels, _}]), 

(* The second set checks is the Path check and the charge path check *)

If[(ITAPath[Adjs, location1, location2, 1]), Continue[];];
        Or[##] & @@ (ChargedPath[Adjs, 
          location1, -directions + #] & /@ (Last@
          Thread@Last@Chargeds[Adjs, location2])) && (! 
         Or[##] & @@ (ChargedPath[Adjs, location2, 
           directions + #] & /@ (Last@
           Thread@Last@Chargeds[Adjs, location1])))))),    

 Adjs[[location1, 2]] = 
  Append[Adjs[[location1, 2]], {Last@labels, location2}];
  (* I don't add the incoming edge, but you could, I don't think it affects the algorithm *)
 (*   Adjs[[location2,2]]=Append[Adjs[[location2,2]],{First@


Now we can finally build the adjacency tree:

{ttt, res3} = AbsoluteTiming@MSGT[leftrightF, updownF, n, 3];


adjtest = 
 Function[{y}, {y[[1]], Function[{x}, x[[2]] ] /@ (y[[2]])}] /@ 
d = DeleteDuplicates[
    Sort[#] & /@ ((#[[1]] -> #[[2]]) & /@ (Flatten[
   Thread[#] & /@ adjtest, 1]))];
GraphPlot[d, VertexCoordinateRules -> Rule @@@ Last@Chargeds[res3, 2]]

enter image description here

Building the image from the Adjacency List

Using VertexRenderingFunction this is rather simple,

GraphPlot[d, VertexCoordinateRules -> Rule @@@ Last@Chargeds[res3, 2],
  VertexRenderingFunction -> (Inset[Image[pics[[#2]], 
     ImageSize -> {41, 41}], #1] &)]

enter image description here

Here is an animation of the building of the puzzle for every step in MSGT,

enter image description here

share|improve this answer
No judgements, but it will take me a while to go through it. – rcollyer Apr 29 at 15:47
@rcollyer: I wonder if anyone can find a method to reconstruct the entire Lenna image from just the piece Mathematica knows... ;) – J. M. is back. May 4 at 3:58
Hey lalmei, I a had just started editing your answer when I saw you were doing the same. Perhaps you can fix the indentation in the first three code blocks? By the way I love the .gif :). – Jacob Akkerboom May 9 at 11:34
@JacobAkkerboom I fixed edited most of the code, but it might need more work, I want to do another update tomorrow. btw the frames for the gif took 6 hours to render ( I used heike's packing algorithm to produce each frame), but I then I ran into the size limit of gifs and had to throw away a third of the frames. Don't know why the exported gif has these edge effects in some of the frames though. – lalmei May 9 at 13:12
@lalmei I made a few more tiny edits. I added AbsoluteTiming somewhere, so that the code now completely works when copy pasting. It looks great! I don't know much about gifs unfortunately. But I often hear that in many cases gifs can be made a lot smaller. The gif looks as if it shouldn't be too big. Great that you had the perseverance to let it finish if it took 6 hours :). – Jacob Akkerboom May 9 at 13:27

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.