Cone of influence

[davidjorg]

Check whether long kymographs work and avoid the 'cone of influence' problem.

[dscho]

I just had an idea: We essentially infer the most likely wavelet weight in the phase() method by evaluating an integral of the form ∫ ∫ f(x)gabor(x - t) dx dt where the assumption is that ∫ gabor(x - t) dx is constant. However, that assumption is incorrect because we only have partial data and therefore integrate only within a specific interval which changes depending on the value of t... Therefore, we should compensate by inferring ∫ (∫ f(x)gabor(x - t) dx / ∫ gabor(x - t) dx) dt instead, where both inner integrals are calculated in the interval for which f(x) is defined.

This is the diff I came up with (to apply on top of 0a507c3):

diff --git a/src/main/java/sc/fiji/timelapse/Phase_Map.java b/src/main/java/sc/fiji/timelapse/Phase_Map.java
index 0855fbe..7360ee4 100644
--- a/src/main/java/sc/fiji/timelapse/Phase_Map.java
+++ b/src/main/java/sc/fiji/timelapse/Phase_Map.java
@@ -26,12 +26,13 @@ public class Phase_Map implements PlugInFilter {

    private double octaveNumber = 4, voicesPerOctave = 50;
    private double gaussSigma = 2, x0 = 100, x1 = 400, sigma0 = 1, sigma1 = 1, subtractionPoint = 50;
-   private boolean plotWaveCounts, showProfileStack, anchorProfileStack, cutTailsFromProfileStack, showPhaseProfileMap;
+   private boolean plotWaveCounts, showProfileStack, anchorProfileStack, cutTailsFromProfileStack, showPhaseProfileMap,
+       compensateBoundaryEffect;

    private ImagePlus imp;

    private double phase(double[] data, int dataSize, double s, int tau) {
-       double wR = 0, wI = 0;
+       double wR = 0, wI = 0, totalR = 0, totalI = 0;

        for (int i = 0; i < dataSize; i++) {
            double u = (i - tau) / s;
@@ -44,8 +45,20 @@ public class Phase_Map implements PlugInFilter {

            wR += data[i] * gaborR;
            wI += data[i] * -gaborI;
+
+           totalR += gaborR;
+           totalI += gaborI;
        }

+       if (compensateBoundaryEffect) {
+           // divide by integral of gabor over the specified interval
+           final double denominator = totalR * totalR + totalI * totalI;
+           if (denominator != 0) {
+               final double r = wR, i = wI;
+               wR = (r * totalR + i * totalI) / denominator;
+               wI = (-r * totalI + i * totalR) / denominator;
+           }
+       }
        // normalization unnecessary:
        // wR /= Math.sqrt(s);
        // wI /= Math.sqrt(s);
@@ -421,6 +434,7 @@ public class Phase_Map implements PlugInFilter {
        gd.addCheckbox("Cut_tails_from_profile_stack i.e. skip spurious signal at tail", cutTailsFromProfileStack);
        gd.addCheckbox("Show_phase_profile_map", showPhaseProfileMap);
        gd.addNumericField("Subtraction_point", subtractionPoint, 0);
+       gd.addCheckbox("Compensate_for_boundary_effect", compensateBoundaryEffect);
        gd.showDialog();
        if (gd.wasCanceled())
            return;
@@ -438,6 +452,7 @@ public class Phase_Map implements PlugInFilter {
        cutTailsFromProfileStack = gd.getNextBoolean();
        showPhaseProfileMap = gd.getNextBoolean();
        subtractionPoint = gd.getNextNumber();
+       compensateBoundaryEffect = gd.getNextBoolean();

        final int width = ip.getWidth(), height = ip.getHeight();
        final Calibration calibration = imp.getCalibration();
diff --git a/src/test/java/sc/fiji/timelapse/PhaseMapTestDrive.java b/src/test/java/sc/fiji/timelapse/PhaseMapTestDrive.java
index 70738fb..5c2cbdb 100644
--- a/src/test/java/sc/fiji/timelapse/PhaseMapTestDrive.java
+++ b/src/test/java/sc/fiji/timelapse/PhaseMapTestDrive.java
@@ -7,6 +7,7 @@ public class PhaseMapTestDrive {
        final String testImage = PhaseMapTestDrive.class.getResource("/kymograph_140816_yfp_s0017.tif").getPath();
        Debug.runFilter(testImage, "Phase Map", " octave_number=3 "
                + "plot_wave_counts "
+               + "compensate_for_boundary_effect "
                + "show_profile_stack anchor_profile_stack cut_tails_from_profile_stack "
                + "show_phase_profile_map ", false);
    }

However, the results are less than convincing: where we had this boundary effect without my change:

it looks like this with my change:

It appears that my attempt to compensate for the boundary artifact introduced another, even worse artifact... 😞