Regardless of the type of data taken in an experiment (voltages versus time, or photographic images, or whatever), there is often a “data processing” stage that precedes a “data analysis” stage. When our data is taken in the form of an image, we have “image processing” (which might involve, say, adjusting the brightness and contrast and, perhaps, a low-pass filter aimed at reducing extraneous noise), and only after the image is deemed ready do we  proceed onward to “image  analysis,” where we try to extract key bits of information that help us to construct a physical story or model.

Image Processing  (expand to see details):

To import an image into Mathematica, one way is just to drop and drag it into the notebook, and then you can surround it with the syntax needed to help Mathematica to properly handle and identify this object:

Im135 = Image[]

This photo, taken from a student’s (Eugene' s) phone, was taken when “grayscale” values applied to the top and bottom halves of the SLM differed by 135 “grayscale units.” [Your computer thinks that the SLM is a computer monitor, a plain old HDMI display, and so the way that we apply voltages to the distinct pixels on the SLM is by sending them each a grayscale value in the range from zero to 255. Inside the SLM itself is an 8-bit CMOS circuit, meaning that it can apply (2^8) distinct voltages to the liquid crystal solution floating above each pixel.]

I find it preferable to separate the 3 - color dataset into the constituent simple arrays that describe the Red, Green, and Blue parts of the image :

RGBIm135 = ColorSeparate[Im135]

Unsurprisingly, since our data was taken using a green laser, the information is essentially captured by the single array describing the green portion of the original three-color image.

GIm135 = RGBIm135[[2]]

Our  original expectation was that our Michelson interferometer would produce uniform "fringes" -- but there is significant “pincushion distortion” in this student’s image. -- Still, if we trim this to only show the central region, these results aren' t too terribly far away from our expectation :

CrIm135 = ImageTrim[GIm135, { {107, 56}, {170, 160} } ]

There are many ways to apply a tilt correction to an image. Given everything that’s going on here (the fringes on the top half are shifted with respect to those on the bottom half), I’m just going to use a Fourier transform, so as to gain a statistical sense of the region shown above

Unfortunately, when we take a Fourier transform,  zero frequency ends up, by default, in the upper left-hand corner of the “map of spatial frequencies” that make up our trimmed image. Although you probably weren’t expecting the way Mathematica wraps the frequencies by default, this is common (LabVIEW and Python and many other packages do the same thing). I prefer to move the zero frequency back to the center of the ensuing map, either by using RotateLeft on the spatial frequency map, ...or by operating on the image above (i.e., the spatial domain), simply multiplying by + /-1^(i, j) as the following:

data135 = ImageData[CrIm135];
{nRow,nCol}=Dimensions[data135];
data135=data135 *(-1)^Table[i+j,{i,nRow},{j,nCol}];
fw135=Fourier[data135,FourierParameters->{1,1}];

Just as an aside, including FourierParameters->{1,1} just means that no normalization pre-factors are added

fw135 contains both  the phases and amplitudes of the various spatial frequencies making up our image. For now, I just want to look at the relative magnitudes of the different spatial frequencies present in the Fourier transform.  If we take a logarithm, this helps us to see small “peaks” in the map of spatial frequencies -- without that, we will be find it difficult to see anything more than just the strongest (zeroth-order) peak. -- Note by adding one to each magnitude, before taking the logarithm, we ensure that all of the resulting values will be positive:

abs135=Log[(1+Abs@fw135)];

Next, we re - scale the output, so that we are, in fact, looking at RELATIVE magnitudes :

RelMagFFT135 = Image[abs135/Max[abs135],ImageSize->300]

Before asking Mathematica to find highlighted features within this figure, I decided to reduce the noise a bit by using a Low - Pass Filter :

SmFFT135=LowpassFilter[RelMagFFT135,1.1]

Now, I' ll ask Mathematica to give me the coordinates of any points it considers to be significant:

ImageKeypoints[SmFFT135]

Clearly, Mathematica finds more points of interest than I do. So, in order to assess which of these highlighted points is of interest to me, I simply overlay them atop the figure :

HighlightImage[SmFFT135,ImageKeypoints[SmFFT135]]

Because of the uniformity of the periodic fringes in our trimmed image, in the Fourier Transform we only expect to see just a couple of peaks (in addition to the “zeroth order” peak in the center. (We expect to see two peaks since a "standing wave" can be envisioned as a superposition of two "counter-propagating" traveling waves.)

Right-clicking allows me to see cursor coordinates, which let’s me sort  out that it is the first two of Mathematica’s highlighted points that correspond to the first-order peaks in this Fourier spectrum. When I then copy the coordinates of those points from the list given by ImageKeypoints[SmFFT135] , Mathematica reports them, to many more significant digits, as:
{20.656536102294922`, 53.146148681640625`}, and 
{45.30823516845703`, 51.8387565612793`}

At this point, you could go so far as to mask all other spectral frequencies and then re - construct the image, which would then be free of noise. See : 
https : // mathematica.stackexchange.com/questions/110914/how - to - use - 2 d - fourier - analysis - to - clean - the - noise - in - an - image 
(Or, as a quick and dirty way of masking other spatial frequencies, you could simple threshold the image, using syntax like: Fred = Threshold[RelMagFFT0, 0.5] )

For now, I’ll just use those points to apply a tilt correction to the original image:

TIm135 = ImageRotate[,ArcTan[(((53.146148681640625`)-(51.8387565612793`))/((45.30823516845703`)-(20.656536102294922`)))]]

Now I' ll zoom back in to a region showing the shift between fringes coming from the top half of the SLM and those coming from the bottom half:

CrIm135 = ImageTrim[TIm135, { {107, 88}, {180, 146} } ]

Again, I' ll only need the green part of the original three - color image :

RGB2Im135 = ColorSeparate[CrIm135]

GrIm135 = RGB2Im135[[2]]

Because the top and bottom halves of the SLM were addressed differently, the fringes from our  Michelson interferometer are displaced, for the top half relative to the bottom half.

Noting that the illumination is not uniform (but is brighter on the far right side), I’ll trim a bit further, but must make sure that when I make separate image sections highlighting the fringes on the top and bottom that they are spanning the exact same set of horizontal coordinates.

ImageDimensions[GrIm135]

Top135 = ImageTrim[GrIm135, { {0, 34}, {53, 59} } ]

Bot135 = ImageTrim[GrIm135, { {0, 0}, {53, 26} } ]

Image Analysis :

Our aim is to produce a plot of the phase shift produced (as a fraction of 2 π) for each grayscale difference applied between the top and bottom halves of the SLM. We could also produce a Look - Up - Table (LUT) to apply to each individual SLM, for automatically converting any desired phase distribution into the corresponding grayscale distribution to use in addressing that particular device.

Now that we' ve adjusted the tilt and confined our attention to an appropriate region of interest, it might seem that we could simply sum the data in each column of these images, and then fit that, in each case, to a sine squared (it is squared because the detected intensity goes as the square of the net local field).  We could then use the extracted fitting parameters to determine the phase offset in each case, thereby producing the desired calibration from grayscale offset to phase difference.

Let' s look at the top section first (expand to see details):

First, we have to convert the image into something that Mathematica recognizes as an array :

TopArr135 = ImageData[Top135];

This is simple array :

ListPlot3D[TopArr135,PlotStyle->RGBColor[0.11`,1.`,0.`]]

Adjusting our eye point, we can see that over a span of 40 pixels there are approximately 9 fringe periods, meaning that the fringe period is something of order (40/9) pixels, ...but since the measured intensity goes as the square of a sine function, our initial estimates of the spatial frequency are 0.5*(9/40) and for the wave number 2π*0.5*(9/40) = 0.707

In Mathematica, we have a row-first convention, and so the columns are associated with the second index. Armed with this information, we see that in order to sum along columns we need to be summing as follows:

TopFringe135 = Total[TopArr135, {1}]

Before I can fit this, I want to be able to estimate more of the fitting parameters. For that, I’ll find the max and min:

Min[TopFringe135]

Max[TopFringe135]

Now I have enough information to provide some ballpark initial estimates on each of the fitting parameters in our simple model of the fringe intensity (I’ve tweaked these a bit until the result looked okay in the graph below):

nlmTop135 = NonlinearModelFit[(TopFringe135), {d+(a2*(Sin[k*x+ϕ])^2)},{ {d,11.3},{a2, 7.06},{k,0.636},{ϕ,2.718}} , x]

Dimensions[TopFringe135]

Show[ListPlot[(TopFringe135)], Plot[nlmTop135 [x], {x, 0, 54}], Frame->True]

It must be said that this approach does not seem terribly robust, in part because there are so very many fitting parameters, and in part because the fringes are not perfectly straight, but MOSTLY because there just aren' t enough pixels per period to provide a clear enough sampling of the waveform

Nevertheless, I' ll demonstrate the method and then, at the end of this notebook, will make recommendations for improvement. What follows are the “best fit” parameters extracted from our nonlinear modeling:

nlmTop135["ParameterTable"]

The phase offset for the top section is best fit (given our initial  estimates on the fitting parameters) to be ϕ = 2.73171  radians, which I would re-express in multiples of π:

pTop135 = ((2.7317057662829254`)/π)

Thus, using this method resulted in an estimated top-half phase offset of (0.869529) *  π

Next, we just have to do the same things to the bottom section (expand to see details) :

BotArr135 = ImageData[Bot135];

ListPlot3D[BotArr135,PlotStyle->RGBColor[0.11`,1.`,0.`]]

BotFringe135 = Total[BotArr135, {1}]

Min[BotFringe135]

Max[BotFringe135]

Using the best-fit k from the top half as our starting estimate for the bottom half, and expecting to add around π to the phase shift, we have :

nlmBot135 = NonlinearModelFit[(BotFringe135), {d+(a2*(Sin[k*x+ϕ])^2)},{ {d,9.9},{a2,5.3},{k,0.635346},{ϕ, 5.86}} , x]

Show[ListPlot[(BotFringe135)], Plot[nlmBot135 [x], {x, 0, 75}], Frame->True]

nlmBot135["ParameterTable"]

For the bottom half, the phase offset (given our initial  estimates on the fitting parameters) is estimated to be ϕ = 6.00698  radians, which I would re-express in multiples of π:

pBot135 = ((6.006981305421901`)/π)

That is, while the top half saw an estimated top-half phase offset of (0.869529) *  π, we find the bottom half has a phase offset of (1.91208)*  π. The phase difference between top and bottom halves is then determined to be 1.04 π

(pBot135 - pTop135)

Concluding remarks/suggestions :

Such methods would yield far more reliable estimates if we (a) remove as much pincushion distortion as possible, while still in the lab (...or, if you must, by applying a mathematical transform to the raw image) and, more importantly, (b) move the imaging screen further away, so as to let the fringe pattern expand out by another factor of two (which would give us more pixels per fringe, and therefore more reliable sampling of the fringe pattern, which would better ensure consistent convergence of our non-linear fitting algorithm,  with less sensitivity to initial estimates of the fitting parameters).