Speckle as an interference phenomenon

Let’s look an illustrative real-world manifestation of this conceptual framework. SAR images are subject to speckle, a type of image distortion or noise which affects pixels in a stochastic, multiplicative way. Why does this occur?

Consider that the returned wave associated with a given pixel likely isn’t the result of a single interation with a single pointwise object. If multiple reflectors (or similarly, a continuous reflector) lie within the pixel, then the returned wave is the sum of multiple individual waves.

Let’s imagine we have a pixel with two identical point reflectors lying within it. How does the returned wave change as a function of where those reflectors are located? Move the blue and red targets in the pixel below:

Code

chart = {
    const svg = d3.create("svg").attr("viewBox", [0, 0, width, height])

    // define a bounding rectangle
    svg.append("rect")
        .attr("x", width/7)
        .attr("y", stroke_width)
        .attr("width", 5*width/7)
        .attr("height", height - 2*stroke_width)
        .attr("stroke", "black")
        .attr("stroke-width", stroke_width)
        .attr("fill-opacity", 0.087124976) //This is a hard-coded value based on the prescribed physical parameters. It's hardcoded because it's easier than dealing with the resulting reactive dependencies since this cell should never be run again

    // define data used for circles
    const circles = d3.range(2).map(i => ({
        x: (i+2) * width / 5,
        y: (i+1) * height / 3,
        index: i, 
    }));

    // define circles as graphic objects
    svg.selectAll("circle")
        .data(circles)
        .join("circle")
            .attr("cx", d => d.x)
            .attr("cy", d => d.y)
            .attr("r", radius)
            .attr("fill", d => d3.schemeCategory10[d.index * 3])
            .attr("id", function(d,i) {return i})
            .call(drag)

    return svg.node();
}

Code

chart2 = {
    const svg = d3.create("svg").attr("viewBox", [0, 0, width, height])

    //define a function to take the canvas position of the leftside chart to the
    //  corresponding region in the second chart
    //  it's messy because I solved the system by hand and didn't simplify, RIP
    //  to future Jake if I ever have to un-hardcode the object position values
    function lineartransform(xposition){
        return (2*height/5 - radius)/(5/7*width -2*radius-2)*xposition + width/7 + 2/5*height + radius/2 - (2*height/5 - radius)*(width/7+radius+1)/(5/7*width-2*radius-2)
    }

    // define a earth curvature line
    svg.append("rect")
        .attr("x", width/7)
        .attr("y", height - 2*stroke_width)
        .attr("width", 5*width/7)
        .attr("height", stroke_width/2)
        .attr("stroke", "black")
        .attr("stroke-width", stroke_width)
        .attr("fill-opacity", 0)

    // define lines going from the target to the satellites (first, so they appear under them)
    svg.append("line")
        .attr("x1", width/7 + 3*height/20)
        .attr("x2", lineartransform(x1))
        .attr("y1", 3*height/20)
        .attr("y2", height - 2*stroke_width)
        .style("stroke-width", stroke_width)
        .style("stroke", "blue")

    svg.append("line")
        .attr("x1", width/7 + 3*height/20)
        .attr("x2", lineartransform(x2))
        .attr("y1", 3*height/20)
        .attr("y2", height - 2*stroke_width)
        .style("stroke-width", stroke_width)
        .style("stroke", "red")


    // define a rectangle representing the satellite
    svg.append("rect")
        .attr("x", width/7 + height/10)
        .attr("y", height/10)
        .attr("width", height/10)
        .attr("height", height/10)
        .attr("stroke", "black")
        .attr("stroke-width", stroke_width)
        .attr("fill-opacity", .50)
    
    // define notches indicating the bounds of the pixel
    svg.append("line")
        .attr("x1", width/7 + 2*height/5)
        .attr("x2", width/7 + 2*height/5)
        .attr("y1", height - 2*stroke_width)
        .attr("y2", height - 2*stroke_width - height/20)
        .style("stroke-width", stroke_width)
        .style("stroke", "black")

    svg.append("line")
        .attr("x1", width/7 + 4*height/5)
        .attr("x2", width/7 + 4*height/5)
        .attr("y1", height - 2*stroke_width)
        .attr("y2", height - 2*stroke_width - height/20)
        .style("stroke-width", stroke_width)
        .style("stroke", "black")

    // define data used for circles
    const circles = d3.range(2).map(i => ({
        x: (i+2) * width / 5,
        y: (i+1) * height / 3,
        index: i, 
    }));

    // define circles as graphic objects
    svg.selectAll("circle")
        .data(circles)
        .join("circle")
            .attr("cx", function(d,i) {
                if (i == 0) {return lineartransform(x1)}
                else {return lineartransform(x2)}
            })
            .attr("cy", height - 2*stroke_width)
            .attr("r", radius/2)
            .attr("fill", d => d3.schemeCategory10[d.index * 3])
            .attr("id", function(d,i) {return i})

    // add text
    svg.append("text")
        .attr("x", width/7)
        .attr("y", .94*height)
        .style("font-size", "28px")
        .text("side view")

    svg.append("text")
        .attr("x", width/2)
        .attr("y", .15*height)
        .style("font-size", "24px")
        .style("font-style", "italic")
        .text("*wildly not to scale")

    return svg.node();
}

Code

formatted_d1 = (0.5*distance_travelled_w1).toFixed(2).replace(/\B(?<!\.\d*)(?=(\d{3})+(?!\d))/g, ",");
formatted_d2 = (0.5*distance_travelled_w2).toFixed(2).replace(/\B(?<!\.\d*)(?=(\d{3})+(?!\d))/g, ",");

Distance from target one to the satellite: m
Distance from target two to the satellite: m

Code

// defines how circles change appearance and position in response to cursor drag
drag = {

    // move the circles to cursor position when dragging is active
    function dragged(event, d){
        var newx = event.x;
        if (newx < xbounds[0]) newx = xbounds[0];
        if (newx > xbounds[1]) newx = xbounds[1];
        var newy = event.y;
        if (newy < ybounds[0]) newy = ybounds[0];
        if (newy > ybounds[1]) newy = ybounds[1];
        
        if (d3.select(this).attr('id') == "0") {
            mutable x1 = newx;
            mutable y1 = newy;
        }
        if (d3.select(this).attr('id') == "1") {
            mutable x2 = newx;
            mutable y2 = newy;
        }
        
        d3.select(this).raise().attr("cx", d.x = newx).attr("cy", d.y = newy);

        var dt1 = 2*Math.sqrt(sat_z**2 + (mutable x1/width*pixel_width - sat_x)**2 + (mutable y1/width*pixel_width - sat_y)**2);
        var dt2 = 2*Math.sqrt(sat_z**2 + (mutable x2/width*pixel_width - sat_x)**2 + (mutable y2/width*pixel_width - sat_y)**2);
        var p1  =  dt1 / wavelength *2 *Math.PI;
        var p2  =  dt2 / wavelength *2 *Math.PI;

        d3.select("rect").attr("fill-opacity", 1-Math.abs(Math.cos((p1-p2)/2)));
        //d3.select("text").text(1-Math.abs(Math.cos((p1-p2)/2))).style("fill","darkOrange");
    }

    // map the d3 drag event functionality to these custom functions
    return d3.drag()
        .on("drag", dragged)
}

Code

mutable x1= 2*width/5;
mutable y1= height/3;
mutable x2= 3*width/5;
mutable y2= 2*height/3;

Code

stroke_width = 7;
width = 700;
height = 500 + 2*stroke_width;
radius = 25;

Code

xbounds = [width/7 + radius + 1, 6*width/7 - radius - 1]
ybounds = [ radius + 1  + stroke_width, height - radius - 1 - stroke_width]

Code

distance = 4;
NumPoints = 1000;
wavelength = .5;
amplitude = 15;
speed = 10;

Code

// display variables
xscale_factor = 100;

Code

sat_x = -100000;
sat_y = 100000;
sat_z = 600000;
pixel_width = 3;

Code

distance_travelled_w1 = 2*Math.sqrt(sat_z**2 + (x1/width*pixel_width - sat_x)**2 + (y1/width*pixel_width - sat_y)**2);
distance_travelled_w2 = 2*Math.sqrt(sat_z**2 + (x2/width*pixel_width - sat_x)**2 + (y2/width*pixel_width - sat_y)**2);

Code

phase_w1 = distance_travelled_w1 / wavelength *2 *Math.PI;
phase_w2 = distance_travelled_w2 / wavelength *2 *Math.PI;

Code

wave1 = {
    var data = [];
    for (var i = 1; i <= NumPoints; i++) {
        var j = i * distance / NumPoints;
        data.push([j*xscale_factor, amplitude * Math.sin(j * 2*Math.PI / wavelength + time + phase_w2)]);
    }
    return data;
}

Code

wave2 = {
    var data = [];
    for (var i = 1; i <= NumPoints; i++) {
        var j = i * distance / NumPoints;
        data.push([j*xscale_factor , amplitude * Math.sin(j * 2*Math.PI / wavelength + time + phase_w1)]);
    }
    return data;
}

Code

waveSum = {
    var data = [];
    for (var i = 1; i <= NumPoints; i++) {
        var j = i * distance / NumPoints;
        data.push([j*xscale_factor , amplitude * Math.sin(j * 2*Math.PI / wavelength + time + phase_w1) + amplitude * Math.sin(j * 2*Math.PI / wavelength + time + phase_w2)]);
    }
    return data;
}

Code

time = {
  let i = 0;
  while (true) {
    i += speed * 0.01;
    yield i
  }
}

Code

r = d3.line()(wave1);
b = d3.line()(wave2);
svg`<svg viewBox="0 -32 400 64">
  <path d="${r}" stroke="red" fill="none" />
  <path d="${b}" stroke="blue" fill="none" />
</svg>`

Code

p = d3.line()(waveSum);
svg`<svg viewBox="0 -64 400 130">
  <path d="${p}" stroke="black" fill="none" />
</svg>`

Code

d3 = require("d3@7")

Here we’ve considered a SAR satellite 600km above Earth, emitting 50cm wavelength radiation with a pixel resolution of 3m.

The distance from each target to the satellite affects where in their cycle each wave is re-intercepted, and thus the difference determines whether those returned waves interfere constructively or deconstructively. Thus a random distribution of targets creates a random difference in distance to the satellite, creating a random level of interference in each pixel. This means the intensity of each pixel is randomly (though not necessarily uniformly) multiplied by a value from 0 (perfect deconstructive interference) to 1 (perfect constructive interference). This is the origin of speckle, and reason why it manifests as multiplicative noise (as opposed to additive).

Let’s now consider two pixels: one with more, stronger scatterers and one with fewer, weaker scatterers – perhaps pixel one belongs to some dense shrubland, while the other is a freshly plowed field. The distribution of scatterers is likely quite random – what happens to the brightness of these pixels under different configurations?

Code

targetnum_1 = 10;
targetnum_2 = 10;
max_radius_p1 = 35;
min_radius_p1 = 25;
max_radius_p2 = 20;
min_radius_p2 = 15;

radii_p1 = {
    var data = [];
        for (var i = 1; i <= targetnum_1; i++) {
            data.push(Math.random()*(max_radius_p1-min_radius_p1) + min_radius_p1);
        }
    return data;
}

radii_p2 = {
    var data = [];
        for (var i = 1; i <= targetnum_2; i++) {
            data.push(Math.random()*(max_radius_p2-min_radius_p2) + min_radius_p2);
        }
    return data;
}

Code

viewof reroll_toggle = Inputs.form(
    [
        Inputs.button("Redistribute Targets"),
        Inputs.button("Auto Redistribute Targets"),
    ],
    {
        template: (inputs) => htl.html`<div style="display: flex; gap: 1em">
        ${inputs}
        </div>`
    })

vtoggle1    = reroll_toggle[0];
auto_toggle = reroll_toggle[1];

Code

time_auto = {
  let i = 0;
  while (auto_toggle%2==1) {
    i += 1;
    yield i
  }
}

Code

pixel1_positions = {
    var dummy  = vtoggle1;
    var dummy2 = time_auto;
    var data = [];
        for (var i = 1; i <= targetnum_1; i++) {
            var x = Math.random()*(5*width/7 - 2*max_radius_p1) + max_radius_p1 + width/7;
            var y = Math.random()*(height - 2*max_radius_p1) + max_radius_p1;
            data.push([x,y]);
        }
    return data
}

pixel2_positions = {
    var dummy  = vtoggle1;
    var dummy2 = time_auto;
    var data = [];
        for (var i = 1; i <= targetnum_2; i++) {
            var x = Math.random()*(5*width/7 - 2*max_radius_p1) + max_radius_p1 + width/7;
            var y = Math.random()*(height - 2*max_radius_p1) + max_radius_p1;
            data.push([x,y]);
        }
    return data
}

Code

//two pixels with targets in them

chart3 = {
    const svg = d3.create("svg").attr("viewBox", [0, 0, width, height])

    // define a bounding rectangle
    svg.append("rect")
        .attr("x", width/7)
        .attr("y", stroke_width)
        .attr("width", 5*width/7)
        .attr("height", height - 2*stroke_width)
        .attr("stroke", "black")
        .attr("stroke-width", stroke_width)
        .attr("fill-opacity", 1-amplitude_p1/200)

    // define data used for circles
    const circles = d3.range(targetnum_1).map(i => ({
        x: pixel1_positions[i][0],
        y: pixel1_positions[i][1],
        index: i, 
    }));

    // define circles as graphic objects
    svg.selectAll("circle")
        .data(circles)
        .join("circle")
            .attr("cx", d => d.x)
            .attr("cy", d => d.y)
            .attr("r", function(d,i) {return radii_p1[i]})
            .attr("stroke", "black")
            .attr("stroke-width", 4)
            .attr("fill", d => d3.schemeCategory10[0])
            .attr("id", function(d,i) {return i})

    return svg.node();
}

Code

chart4 = {
    const svg = d3.create("svg").attr("viewBox", [0, 0, width, height])

    // define a bounding rectangle
    svg.append("rect")
        .attr("x", width/7)
        .attr("y", stroke_width)
        .attr("width", 5*width/7)
        .attr("height", height - 2*stroke_width)
        .attr("stroke", "black")
        .attr("stroke-width", stroke_width)
        .attr("fill-opacity", 1-amplitude_p2/200)

    // define data used for circles
    const circles = d3.range(targetnum_2).map(i => ({
        x: pixel2_positions[i][0],
        y: pixel2_positions[i][1],
        index: i, 
    }));

    // define circles as graphic objects
    svg.selectAll("circle")
        .data(circles)
        .join("circle")
            .attr("cx", d => d.x)
            .attr("cy", d => d.y)
            .attr("r", function(d,i) {return radii_p2[i]})
            .attr("stroke", "black")
            .attr("stroke-width", 4)
            .attr("fill", d => d3.schemeCategory10[0])
            .attr("id", function(d,i) {return i})

    return svg.node();
}

Code

Plot.plot({
    round: true,
    color: {legend: true},
    x: {label: "Brightness"},
    y: {label: "Frequency"},
    marks: [
        Plot.rectY(brightnesses, Plot.binX({y2: "count"}, {x: "x", fill: "p", mixBlendMode: "multiply"})),
        Plot.ruleY([0])
        ]
})

Code

amplitude_p1 = {

    var dummy = vtoggle1;
    
    var amplitude_naught = radii_p1[0];

    var dt_naught =
        2*Math.sqrt(
            sat_z**2
            + (pixel1_positions[0][0]/width*pixel_width - sat_x)**2
            + (pixel1_positions[0][1]/width*pixel_width - sat_y)**2
        );
    
    var phase_naught =
        dt_naught / wavelength *2 *Math.PI;
    
    //successively add the waves, determining the new phases and amplitudes
    for (var i = 1; i <= targetnum_1-1; i++) {
        //compute the phase of the next subwave
        var dt_next =
            2*Math.sqrt(
                sat_z**2
                + (pixel1_positions[i][0]/width*pixel_width - sat_x)**2
                + (pixel1_positions[i][1]/width*pixel_width - sat_y)**2
        );

        var phase_next = 
            dt_next / wavelength *2 *Math.PI;
        
        //compute the amplitude of the new sum
        var new_amplitude =
            Math.sqrt(
                amplitude_naught**2
                + radii_p1[i]**2
                + 2*amplitude_naught*radii_p1[i]*Math.cos(phase_next-phase_naught)
            );

        //compute the phase of the new sum
        var new_phase =
            Math.atan(
                (amplitude_naught*Math.sin(phase_naught) + radii_p1[i]*Math.sin(phase_next))
                /(amplitude_naught*Math.cos(phase_naught) + radii_p1[i]*Math.cos(phase_next))
            )

        amplitude_naught = new_amplitude;
        phase_naught = new_phase;

    }
    return amplitude_naught
}

amplitude_p2 = {

    var dummy = vtoggle1;
    
    var amplitude_naught = radii_p2[0];

    var dt_naught =
        2*Math.sqrt(
            sat_z**2
            + (pixel2_positions[0][0]/width*pixel_width - sat_x)**2
            + (pixel2_positions[0][1]/width*pixel_width - sat_y)**2
        );
    
    var phase_naught =
        dt_naught / wavelength *2 *Math.PI;
    
    //successively add the waves, determining the new phases and amplitudes
    for (var i = 1; i <= targetnum_1-1; i++) {
        //compute the phase of the next subwave
        var dt_next =
            2*Math.sqrt(
                sat_z**2
                + (pixel2_positions[i][0]/width*pixel_width - sat_x)**2
                + (pixel2_positions[i][1]/width*pixel_width - sat_y)**2
        );

        var phase_next = 
            dt_next / wavelength *2 *Math.PI;
        
        //compute the amplitude of the new sum
        var new_amplitude =
            Math.sqrt(
                amplitude_naught**2
                + radii_p2[i]**2
                + 2*amplitude_naught*radii_p2[i]*Math.cos(phase_next-phase_naught)
            );

        //compute the phase of the new sum
        var new_phase =
            Math.atan(
                (amplitude_naught*Math.sin(phase_naught) + radii_p2[i]*Math.sin(phase_next))
                /(amplitude_naught*Math.cos(phase_naught) + radii_p2[i]*Math.cos(phase_next))
            )

        amplitude_naught = new_amplitude;
        phase_naught = new_phase;

    }
    return amplitude_naught
}

{
mutable brightnesses.push({x: amplitude_p1, p: "Pixel 1"});
mutable brightnesses.push({x: amplitude_p2, p: "Pixel 2"});
mutable brightnesses = mutable brightnesses;
}

Code

//histogram to be populated with the brightness values of the two pixels
mutable brightnesses = [];

Selecting the Auto Redistribute option to generate a few hundred random configurations of each pixel type, we can see in the developing histogram of their brightnesses that while the pixel with more strongly reflecting targets is on average brighter, there are many times when the pixel with weaker targets actually gives a higher return. This is due to the random level of interference we’ve called speckle. We can exploit the fact that the average value remains higher to identify likely high-return areas in speckle filtering. As an aside, notice that the distributions of brightness appear relatively Gaussian – this is a lovely invocation of the Central Limit Theorem.