Intrinsic image decomposition faces the long-standing challenge from the coupling of the components of the image-the surface albedo, direct illumination, and ambient illumination in the observed image. Without knowing the absolute values of the image components, we propose inferring shading by ordering pixels by relative brightness. The pairwise shading orders are estimated in two ways: brightness order and low-order fittings of the local shading field. The brightness order is a nonlocal metric that can be used to compare any two pixels, including those with different reflectances and shadings. Low-order fittings are used for pixel pairs within local regions of smooth shading. They can capture the global order structure and local variations in the shading when used together. To integrate the pairwise orders into a globally consistent order, we propose a Consistency-aware Selective Fusion method. The iterative selection process solves the inconsistencies between pairwise orders obtained using different estimation methods. To avoid polluting the global order, inconsistent or unreliable pairwise orders will be automatically excluded from the fusion. Experimental results show that the proposed model effectively recovers the shading, including deep shadows, on the MIT Intrinsic Image dataset. Moreover, our model works well on natural images from the IIW, UIUC Shadow, and NYU-depth datasets, where the colors of direct lights and ambient lights are quite different.