Note: MaTeX is a LATEX typesetting package for Mathematica which can be installed from the sitehttps://github.com/szhorvat/MaTeX.

<<MaTeX` (* only used in labeling some figures *)

Naive difference-of-projections flow

Let’s consider the simple example where the two sets A and B represent two intersecting straight lines: y = m x and y = 0, respectively (m ≠ 0). The two projections (Subscript[P, A]) and (Subscript[P, B]) map a point (x,y) to the nearest point on lines A and B. Thus, (Subscript[P, A])(x,y)=(((x+m y)/((m^2)+1)),m((x+m y)/((m^2)+1)))  and  (Subscript[P, B])(x,y)=(x,0).

PA[{x_,y_},m_]={1,m}*((x+m*y)/((m^2)+1))

PB[{x_,y_}]={x,0}

Check that the two projections yield sensible results:

With[{m=0.8},
Plot[{m*x,0},{x,-2,2},AspectRatio->1,Frame->True,Axes->False,PlotStyle->Black,PlotRange->{-2,2},FrameTicksStyle->Directive[Black,16],FrameLabel->{MaTeX["x",Magnification->2],MaTeX["y",Magnification->2]},ImageSize->400,Epilog->Join[{Red,PointSize[0.018]},{Point[{#,PA[#,m],PB[#]}],Line[{PA[#,m],#,PB[#]}]}&@RandomReal[{-2,2},2]]]
]

Let’s see what happens if we evolve the dynamics with the difference of projections:

With[{m=0.8},
StreamPlot[PA[{x,y},m]-PB[{x,y}],{x,-2,2},{y,-2,2},PlotRangePadding->None,FrameTicksStyle->Directive[Black,16],StreamColorFunction->"Rainbow",StreamStyle->AbsoluteThickness[1.],ImageSize->450,FrameLabel->{MaTeX["x",Magnification->2],MaTeX["y",Magnification->2]},StreamStyle->AbsoluteThickness[1.1],PlotLegends->Placed[BarLegend[Automatic,LabelStyle->{Black,FontSize->14},LegendMarkerSize->220,LegendLabel->MaTeX["\\left\\lVert \\dot{\\vec{r}} \\right\\rVert",Magnification->1.6]],{1.04,0.56}],Epilog->{Line[{{-2,0},{2,0}}],Line[{{-2,-2*m},{2,2*m}}],PointSize[0.024],Point[{{0,0}}]},ImagePadding->{{50,2},{50,2}}]]

We see that although the desired solution (0,0) is a fixed point, it is a saddle. Hence, it is unstable and the system will fail to find the solution.

On the other hand, with the RRR algorithm, the system spirals in toward (0,0):

With[{m=0.8},
StreamPlot[PA[2*PB[{x,y}]-{x,y},m]-PB[{x,y}],{x,-2,2},{y,-2,2},PlotRangePadding->None,FrameTicksStyle->Directive[Black,16],StreamColorFunction->"Rainbow",StreamStyle->AbsoluteThickness[1.],ImageSize->450,FrameLabel->{MaTeX["x",Magnification->2],MaTeX["y",Magnification->2]},StreamStyle->AbsoluteThickness[1.1],PlotLegends->Placed[BarLegend[Automatic,LabelStyle->{Black,FontSize->14},LegendMarkerSize->220,LegendLabel->MaTeX["\\left\\lVert \\dot{\\vec{r}} \\right\\rVert",Magnification->1.6]],{1.04,0.56}],Epilog->{Line[{{-2,0},{2,0}}],Line[{{-2,-2*m},{2,2*m}}],PointSize[0.024],Point[{{0,0}}]},ImagePadding->{{50,2},{50,2}}]]

Square and Line

Projection (Subscript[P, A]) maps (x,y) to the nearest of the 4 points (1,1), (-1,1), (-1,-1), and (1,-1). Thus, we simply have (Subscript[P, A])(x,y)=(Sign(x), Sign(y)). Conversely, projection (Subscript[P, B]) maps (x,y) to the nearest point on the line y = x. Thus, (Subscript[P, B])(x,y)=(((x+y)/2),((x+y)/2)) .

Pa[{x_,y_}]=Sign[{x,y}]

Pb[{x_,y_}]=Table[((x+y)/2),2]

The first term in the RRR flow can be simplified as follows:

Pa[2*Pb[{x,y}]-{x,y}]

Pa[{y,x}]

Let’s see what the RRR flow looks like:

StreamPlot[10*(Pa[{y,x}]-Pb[{x,y}]),{x,-2.6,2.6},{y,-2.6,2.6},StreamStyle->AbsoluteThickness[1.05],StreamPoints->{Join[Flatten[(0.1*#*{-1,1}&/@Complement[Range[-10,10],{0}])+#&/@{1.001,0.999,-0.999,-1.001},1],{{0,-0.001},{0,0.001},{1.001,1.001},{-1.001,-1.001}}],0.1,1000},PlotRangePadding->None,StreamColorFunction->"Rainbow",FrameTicksStyle->Directive[Black,16],ImageSize->450,FrameLabel->{MaTeX["x",Magnification->2],MaTeX["y",Magnification->2]},PlotLegends->Placed[BarLegend[Automatic,LabelStyle->{Black,FontSize->14},LegendMarkerSize->220,LegendLabel->MaTeX["10\\left\\lVert \\dot{\\vec{r}} \\right\\rVert",Magnification->1.6]],{1.04,0.55}],ImagePadding->{{50,2},{50,2}},Epilog->{GrayLevel[0.8],Line[{{-2,0},{2,0}}],Line[{{0,-2},{0,2}}],Black,AbsoluteThickness[1.02],Line[{{-3,-3},{3,3}}],PointSize[0.023],Point[{{-1,-1},{1,1},{1,-1},{-1,1}}],AbsoluteThickness[1.4],Dashing[{0.025,0.012}],Line[{{0,2},{2,0}}],Line[{{0,-2},{-2,0}}]}]

The fixed points are indicated by dashed black lines. They are given by the set S=(Subscript[S, 1])⋃(Subscript[S, 2]) where (Subscript[S, 1])={(x,1-x):0<x<1} and (Subscript[S, 2])={(x,-1-x):-1<x<0}. The points in (Subscript[S, 1]) represent the solution at (1,1), whereas the points in (Subscript[S, 2]) represent the solution at (-1,-1). We find that all trajectories lead to one of the fixed points (except possibly those on the separatrix x+y = 0, which constitute a set of measure zero).

(*StreamPlot[100*(Pa[{y,x}]-Pb[{x,y}]),{x,-2,2},{y,-2,2},StreamPoints->{Join[Tuples[{Range[-1.8,1.8,0.3],{-0.001,0.001}}],Tuples[{{-0.001,0.001},Range[-1.8,1.8,0.3]}],{{0,0.001}}],0.1,1000},PlotRangePadding->None,StreamColorFunction->"Rainbow"]*)

Histogram projection

1st part:

The constraint set S is composed of all permutations of the final sorted pixel values (OverVector[p])≡{(Subscript[p, 1]),(Subscript[p, 2]),..., (Subscript[p, N])}. We need to show that the Euclidean distance ‖(OverVector[x])-(OverVector[q])‖ : (OverVector[q])∈S is minimized for (OverVector[q])=r((OverVector[p])) where r denotes the ordering of (OverVector[x]). This is equivalent to showing that the Euclidean distance between two given arrays of equal length N is minimized by sorting them in the same order.

We will prove the above proposition by induction. First, consider N=2. Take two arrays (OverVector[a])={(Subscript[a, 1]),(Subscript[a, 2])} and (OverVector[b])={(Subscript[b, 1]),(Subscript[b, 2])} with (Subscript[a, 1])<(Subscript[a, 2]) and (Subscript[b, 1])<(Subscript[b, 2]). We can compare the distances for the two different orderings as (Subscript[d, {(Subscript[a, 1]),(Subscript[a, 2])},{(Subscript[b, 1]),(Subscript[b, 2])}]^2)-(Subscript[d, {(Subscript[a, 1]),(Subscript[a, 2])},{(Subscript[b, 2]),(Subscript[b, 1])}]^2)=((((Subscript[a, 1])-(Subscript[b, 1])))^2)+((((Subscript[a, 2])-(Subscript[b, 2])))^2)-((((Subscript[a, 1])-(Subscript[b, 2])))^2)+((((Subscript[a, 2])-(Subscript[b, 1])))^2)=2((Subscript[a, 1])-(Subscript[a, 2]))((Subscript[b, 2])-(Subscript[b, 1]))<0. Hence, (Subscript[d, {(Subscript[a, 1]),(Subscript[a, 2])},{(Subscript[b, 1]),(Subscript[b, 2])}])<(Subscript[d, {(Subscript[a, 1]),(Subscript[a, 2])},{(Subscript[b, 2]),(Subscript[b, 1])}]). Thus, the proposition is true for N=2.

Now consider N=3. Take (OverVector[a])={(Subscript[a, 1]),(Subscript[a, 2]),(Subscript[a, 3])} and (OverVector[b])={(Subscript[b, 1]),(Subscript[b, 2]),(Subscript[b, 3])}, where (Subscript[a, 1])<(Subscript[a, 2])<(Subscript[a, 3]) and (Subscript[b, 1])<(Subscript[b, 2])<(Subscript[b, 3]). There are 6 permutations of (OverVector[b]). Any such permutation P((OverVector[b])) is either a cyclic rotation of (OverVector[b]) [such as {(Subscript[b, 3]),(Subscript[b, 1]),(Subscript[b, 2])}] or differs from (OverVector[b]) in at most 2 positions [such as {(Subscript[b, 2]),(Subscript[b, 1]),(Subscript[b, 3])}]. For the latter case, the problem reduces to comparing distances between N=2 arrays. Let’s consider the former situation: consider a cyclic rotation of (OverVector[b]) by n sites [1≤n≤N-1, here N=3] and call it (Subscript[C, n])((OverVector[b])). Then

MaTeX["\\frac{1}{2}\\Big(d^2_{\\vec{a},\\vec{b}}-d^2_{\\vec{a},C_n(\\vec{b})}\\Big) = \\sum_{j=1}^{N-n} a_j \\underbrace{(b_{j+n}-b_j)}_{>0}+\\sum_{j=N-n+1}^{N} a_j \\underbrace{(b_{j+n-N}-b_j)}_{<0} \\;<\\; a_{N-n+1}\\left[\\sum_{j=1}^{N-n} b_{j+n}-b_j+\\sum_{j=N-n+1}^{N} b_{j+n-N}-b_j\\right]=0\\;.",Magnification->1.7]

Hence, (Subscript[d, (OverVector[a]),(OverVector[b])])<(Subscript[d, (OverVector[a]),(Subscript[C, n])((OverVector[b]))]). Thus, we have proved the proposition for N=3. These arguments readily generalize to larger values of N. Any permutation of (OverVector[b])={(Subscript[b, 1]),(Subscript[b, 2]),...,(Subscript[b, N])} is either a cyclic rotation or differs from (OverVector[b]) in at most N-1 positions. Either case yields a larger Euclidean distance from (OverVector[a])={(Subscript[a, 1]),(Subscript[a, 2]),...,(Subscript[a, N])}  [where (Subscript[a, 1])<(Subscript[a, 2])<...<(Subscript[a, N]) and (Subscript[b, 1])<(Subscript[b, 2])<...<(Subscript[b, N])]. Therefore, the proposition holds for all integers N≥2 (QED).

2nd part

Import the Escher image:

Escher=Import[NotebookDirectory[]<>"DrawingHands.jpg"]

Let’s see the pixel distribution:

Histogram[Flatten[ImageData[Escher]],50,Frame->True,FrameTicksStyle->Directive[Black,14],ImageSize->420,FrameLabel->{Style["pixel values",Black,18],Style["number of pixels",Black,18]},PlotRange->{-250,7000}]

We use the projection algorithm to change the pixel distribution to a uniform one:

Block[{pixels,ordering,newsortedpixelvalues,rowlength},
pixels=ImageData[Escher]; (* pixel values in the original image *)
ordering=Ordering[Ordering[Flatten[pixels]]]; (* ordering of pixels: positions when sorted *)
newsortedpixelvalues=Subdivide[0.,1.,Length[Flatten[pixels]]-1]; (* sorted pixel values in the projected image  *)
rowlength=Last[Dimensions[pixels]]; (* number of pixels in a horozontal slice *)
projectedimage=Image[Partition[newsortedpixelvalues[[ordering]],rowlength]] (* construction of the projected image *)
]

Let’s check the new pixel distribution:

Histogram[Flatten[ImageData[projectedimage]],50,Frame->True,FrameTicksStyle->Directive[Black,14],ImageSize->420,FrameLabel->{Style["pixel values",Black,18],Style["number of pixels",Black,18]},PlotRange->{{-0.02,1.02},{-65,2050}},PlotRangePadding->None]

Escher and his uniform grayscale projection:

Row[{Column[{Style["Original",18],Show[Escher,ImageSize->320]},Alignment->Center],Column[{Style["Projected",18],Show[projectedimage,ImageSize->320]},Alignment->Center]},"         "]

---------------The End---------------