Fine shift density

src/c_modules/fine_movement/libimagedeform/include/image_deform.h

@@ -123,12 +123,18 @@ void image_deform_apply_point(const struct ImageDeform *deform,
  * \param deform deformation structure
  * \param input_image original image
  * \param output_image resulting image
- * \param subpixels upscaling coefficient of image during deform
  */
 void image_deform_apply_image(const struct ImageDeform *deform,
                               const struct ImageGrid *input_image,
-                              struct ImageGrid *output_image,
-                              int subpixels);
+                              struct ImageGrid *output_image);
+
+/**
+ * \brief Generate image deform density
+ * \param deform deformation sructure
+ * \param divergence output image for divergence
+ */
+void image_deform_calculate_divergence(const struct ImageDeform *deform,
+                                       struct ImageGrid *divergence);
 
 #ifdef __cplusplus
 }

src/c_modules/fine_movement/libimagedeform/src/image_deform.c

@@ -157,17 +157,78 @@ void image_deform_apply_point(const struct ImageDeform *deform,
 
 void image_deform_apply_image(const struct ImageDeform *deform,
                               const struct ImageGrid *input_image,
-                              struct ImageGrid *output_image,
-                              int subpixels)
+                              struct ImageGrid *output_image)
 {
     int y, x;
+    double kx = (double)input_image->w / output_image->w;
+    double ky = (double)input_image->h / output_image->h;
     for (y = 0; y < output_image->h; y++)
         for (x = 0; x < output_image->w; x++)
         {
             double orig_y, orig_x;
-            image_deform_apply_point(deform, (double)x/subpixels, (double)y/subpixels, &orig_x, &orig_y);
+            image_deform_apply_point(deform, x*kx, y*ky, &orig_x, &orig_y);
 
             double val = image_grid_get_pixel(input_image, orig_x, orig_y);
             image_grid_set_pixel(output_image, x, y, val);
         }
 }
+
+void image_deform_calculate_divergence(const struct ImageDeform *deform,
+                                       struct ImageGrid *divergence)
+{
+    int y, x;
+    double kx = (double)deform->image_w / divergence->w;
+    double ky = (double)deform->image_h / divergence->h;
+    // Calculate density from original coordinates
+    for (y = 0; y < divergence->h; y++)
+    {
+        double dy1, dy2;
+        if (y == 0)
+        {
+            dy1 = 0;
+            dy2 = 1;
+        }
+        else if (y == divergence->h-1)
+        {
+            dy1 = -1;
+            dy2 = 0;
+        }
+        else
+        {
+            dy1 = -1;
+            dy2 = 1;
+        }
+
+        for (x = 0; x < divergence->w; x++)
+        {
+            double dx1, dx2;
+            if (x == 0)
+            {
+                dx1 = 0;
+                dx2 = 1;
+            }
+            else if (x == divergence->w-1)
+            {
+                dx1 = -1;
+                dx2 = 0;
+            }
+            else
+            {
+                dx1 = -1;
+                dx2 = 1;
+            }
+
+            double vx1 = image_deform_get_shift(deform, (x+dx1)*kx, y*ky, 0);
+            double vx2 = image_deform_get_shift(deform, (x+dx2)*kx, y*ky, 0);
+            double ddx = (vx2 - vx1) / (dx2-dx1);
+
+            double vy1 = image_deform_get_shift(deform, x*kx, (y+dy1)*ky, 0);
+            double vy2 = image_deform_get_shift(deform, x*kx, (y+dy2)*ky, 0);
+            double ddy = (vy2 - vy1) / (dy2-dy1);
+
+            double div = ddx + ddy;
+
+            image_grid_set_pixel(divergence, x, y, div);
+        }
+    }
+}

src/c_modules/fine_movement/module/imagedeform.c

@@ -138,10 +138,34 @@ static PyObject *ImageDeform_apply_image(PyObject *_self,
         (struct ImageGridObject *)PyObject_CallObject((PyObject *)&ImageGrid, argList);
     Py_DECREF(argList);
 
-    image_deform_apply_image(&self->deform, &in_img->grid, &out_img->grid, subpixels);
+    image_deform_apply_image(&self->deform, &in_img->grid, &out_img->grid);
     return (PyObject *)out_img;
 }
 
+static PyObject *ImageDeform_divergence(PyObject *_self,
+                                         PyObject *args,
+                                         PyObject *kwds)
+{
+    struct ImageDeformObject *self = (struct ImageDeformObject *)_self;
+    static char *kwlist[] = {"subpixels", NULL};
+    int subpixels;
+    if (!PyArg_ParseTupleAndKeywords(args, kwds, "i", kwlist, &subpixels))
+    {
+        PyErr_SetString(PyExc_ValueError, "invalid function arguments");
+        Py_INCREF(Py_None);
+        return Py_None;
+    }
+
+    PyObject *argList = Py_BuildValue("ii", self->deform.image_w*subpixels, self->deform.image_w*subpixels);
+    struct ImageGridObject *out_img =
+        (struct ImageGridObject *)PyObject_CallObject((PyObject *)&ImageGrid, argList);
+    Py_DECREF(argList);
+
+    image_deform_calculate_divergence(&self->deform, &out_img->grid);
+    return (PyObject *)out_img;
+}
+
+
 static PyObject *ImageDeform_apply_point(PyObject *_self,
                                          PyObject *args,
                                          PyObject *kwds)
@@ -168,7 +192,9 @@ static PyMethodDef ImageDeform_methods[] = {
      "Return image deform content as numpy array"},
     {"apply_image", (PyCFunction)ImageDeform_apply_image, METH_VARARGS | METH_KEYWORDS,
      "Apply ImageDeform to ImageGrid"},
-     {"apply_point", (PyCFunction)ImageDeform_apply_point, METH_VARARGS | METH_KEYWORDS,
+    {"divergence", (PyCFunction)ImageDeform_divergence, METH_VARARGS | METH_KEYWORDS,
+     "Divergence of ImageDeform"},
+    {"apply_point", (PyCFunction)ImageDeform_apply_point, METH_VARARGS | METH_KEYWORDS,
      "Apply ImageDeform to point"},
     {NULL} /* Sentinel */
 };

tests/test_fine_shift.py

@@ -267,3 +267,58 @@ def test_memory_leak():
         deltas.append(memory2 - memory1 - prevd)
         prevd = memory2 - memory1
     assert deltas[-1] == 0
+
+def test_divergence_1():
+    df1 = next(readjpeg(os.path.join(dir_path, "fine_shift/image3.png")))
+    df2 = next(readjpeg(os.path.join(dir_path, "fine_shift/image3.png")))
+
+    image = df1.get_channel("L")[0].astype('double')
+    image_ref = df2.get_channel("L")[0].astype('double')
+
+    w = image.shape[1]
+    h = image.shape[0]
+
+    grid = ImageGrid(image_w=w, image_h=h)
+    grid.fill(image)
+    grid_ref = ImageGrid(image_w=w, image_h=h)
+    grid_ref.fill(image_ref)
+
+    lc = ImageDeformLC(image_w=w, image_h=h, pixels=1)
+    deform = lc.find(grid, None, grid_ref, None, 5, 5, 1)
+    assert deform is not None
+
+    divergence = deform.divergence(subpixels=1)
+    assert divergence is not None
+    divergence = divergence.content()
+    assert len(divergence.shape) == 2
+    assert divergence.shape[0] == h
+    assert divergence.shape[1] == w
+    assert np.amin(divergence) == 0
+    assert np.amax(divergence) == 0
+
+def test_divergence_2():
+    df1 = next(readjpeg(os.path.join(dir_path, "fine_shift/image3.png")))
+    df2 = next(readjpeg(os.path.join(dir_path, "fine_shift/image4.png")))
+
+    image = df1.get_channel("L")[0].astype('double')
+    image_ref = df2.get_channel("L")[0].astype('double')
+
+    w = image.shape[1]
+    h = image.shape[0]
+
+    grid = ImageGrid(image_w=w, image_h=h)
+    grid.fill(image)
+    grid_ref = ImageGrid(image_w=w, image_h=h)
+    grid_ref.fill(image_ref)
+
+    lc = ImageDeformLC(image_w=w, image_h=h, pixels=1)
+    deform = lc.find(grid, None, grid_ref, None, 5, 5, 1)
+    assert deform is not None
+
+    divergence = deform.divergence(subpixels=1)
+    assert divergence is not None
+    divergence = divergence.content()
+
+    assert len(divergence.shape) == 2
+    assert divergence.shape[0] == h
+    assert divergence.shape[1] == w