Updated f2py section.

Author: jonaslindemann

f2py/.vscode/settings.json

@@ -0,0 +1,4 @@
+{
+    "fortran.fortls.disabled": true,
+    "python.terminal.activateEnvInCurrentTerminal": true
+}

f2py/arr1.py

@@ -2,15 +2,10 @@
 
 import os
 import numpy as np
-import fortran_utils as utils
+import mf_mixed
 
-# Add the path to the directory containing the Fortran runtime DLL to the system path
-
-utils.f2py_setup()
-
-if utils.f2py_compile("fortmod", ["arr1.f90"])!=0:
-    print("Compilation failed")
-    exit()
+mf_mixed.f2py_setup()
+mf_mixed.f2py_compile('fortmod', ['arr1.f90'])
 
 from fortmod import *
 

f2py/arr2.py

@@ -1,5 +1,7 @@
-import sys
-sys.path.append("./build") 
+import mf_mixed
+
+mf_mixed.f2py_setup()
+mf_mixed.f2py_compile('fortmod', ['arr2.f90'])
 
 from numpy import *
 from fortmod import *

f2py/filters.f90

@@ -0,0 +1,51 @@
+module filters 
+    use iso_c_binding
+
+    implicit none
+  
+  contains
+  
+  subroutine mean_filter(image, kernel_size, filtered_image)
+    use, intrinsic :: iso_c_binding, only: c_double
+    implicit none
+    real(c_double), intent(in) :: image(:,:)
+    integer, intent(in) :: kernel_size
+    real(c_double), intent(inout) :: filtered_image(:,:)
+    integer :: y, x, i, j, height, width, y_start, y_end, x_start, x_end
+    real(c_double) :: sum_value, divisor
+    integer :: window_size
+  
+    height = size(image, 1)
+    width = size(image, 2)
+  
+    ! Precompute divisor
+    divisor = 1.0_c_double / (kernel_size * kernel_size)
+  
+     !$OMP PARALLEL DO PRIVATE(x, i, j, sum_value, y_start, y_end, x_start, x_end, window_size) SCHEDULE(static)
+    do y = 1, height
+      do x = 1, width
+        ! Calculate window boundaries
+        y_start = max(1, y - kernel_size/2)
+        y_end = min(height, y + kernel_size/2)
+        x_start = max(1, x - kernel_size/2)
+        x_end = min(width, x + kernel_size/2)
+  
+        ! Calculate the sum value of the window
+        sum_value = 0.0_c_double
+        window_size = 0
+        
+        do j = y_start, y_end
+          do i = x_start, x_end
+            sum_value = sum_value + image(j, i)
+            window_size = window_size + 1
+          end do
+        end do
+  
+        ! Assign the mean value to the corresponding pixel in the filtered image
+        filtered_image(y, x) = sum_value / window_size
+      end do
+    end do
+    !$OMP END PARALLEL DO
+  
+  end subroutine mean_filter  
+  end module filters

f2py/fortran_python.png

@@ -0,0 +1,104 @@
+import os
+import time
+import cv2
+import sys
+import numpy as np
+import matplotlib.pyplot as plt
+
+import mf_mixed
+
+
+
+import numpy as np
+from numpy.lib.stride_tricks import as_strided
+
+def blur_filter_vectorized(image, kernel_size=3):
+    # Create the kernel
+    kernel = np.ones((kernel_size, kernel_size)) / (kernel_size ** 2)
+    
+    # Get dimensions
+    i_height, i_width = image.shape
+    k_height, k_width = kernel.shape
+    
+    # Pad the image
+    pad_height = k_height // 2
+    pad_width = k_width // 2
+    padded_image = np.pad(image, ((pad_height, pad_height), (pad_width, pad_width)), mode='edge')
+    
+    # Create a view of the padded image with sliding windows
+    windows = as_strided(padded_image,
+                         shape=(i_height, i_width, k_height, k_width),
+                         strides=padded_image.strides + padded_image.strides)
+    
+    # Perform the convolution
+    return np.einsum('ijkl,kl->ij', windows, kernel)
+
+
+if __name__ == "__main__":
+
+    mf_mixed.f2py_setup()
+    mf_mixed.f2py_compile('filters', ['filters.f90'])
+
+    from filters import *   # Import the Fortran module
+
+    # Load the image in grayscale mode
+
+    print("Loading image...")
+    image = cv2.imread('fortran_python.png', cv2.IMREAD_GRAYSCALE)
+    print(image.dtype)
+
+    print("Creating fortran array...")
+    image_f = np.zeros_like(image, order='F', dtype=float)
+    image_f[:] = image
+
+    # Define the kernel size for the mean filter
+
+    kernel_size = 7
+
+    # Apply the mean filter to the image
+
+    print("Applying mean filter (python)...")
+
+    t0 = time.time()
+    filtered_image1 = blur_filter_vectorized(image_f, kernel_size)
+    t1 = time.time()
+
+    print("Time (python): ", t1 - t0)
+
+    print("Applying mean filter (fortran)...")
+    filtered_image2 = np.zeros_like(image_f, order='F', dtype=float)
+
+    n_samples = 100
+
+    t0 = time.time()
+    for i in range(n_samples):
+        filters.mean_filter(image_f, kernel_size, filtered_image2)
+    t1 = time.time()
+
+    print("Time (fortran): ", (t1 - t0) / n_samples)
+
+    print("Applying mean filter (cv2)...")
+
+    t0 = time.time()
+    for i in range(1):
+        filtered_image3 = cv2.blur(image_f, (kernel_size, kernel_size))
+    t1 = time.time()
+
+    print("Time (python cv2): ", t1 - t0)
+
+    plt.figure()
+    plt.title('Python')
+    plt.imshow(filtered_image1, cmap='gray')
+    plt.figure()
+    plt.title('Fortran')
+    plt.imshow(filtered_image2, cmap='gray')
+    plt.figure()
+    plt.title('OpenCV')
+    plt.imshow(filtered_image3, cmap='gray')
+    plt.show()
+
+    # Display the original and filtered images
+    #cv2.imshow('Original Image', image)
+    #cv2.imshow('Mean Filtered Image', filtered_image)
+    #cv2.waitKey(0)
+    #cv2.destroyAllWindows()

f2py/fortran_utils.py

@@ -0,0 +1,21 @@
+# -*- coding: utf-8 -*-
+
+import os
+
+# ---- 1. Replace the following path with corresponding path in your system ----
+f2py_mingw_path = "C:\\Qt\\Tools\\mingw1120_64\\bin"
+
+# ---- 2. Replace the following flags with corresponding flags in your system ----
+#f2py_opt_flags = '-Wall -Wextra -Wimplicit-interface -fPIC -O3 -march=native -ffast-math -funroll-loops -fopenmp'
+f2py_opt_flags = '-Wall -Wextra -Wimplicit-interface -fPIC -O3 -march=native -ffast-math -funroll-loops -fopenmp -ffree-line-length-none'
+#f2py_opt_flags = '-Wall -Wextra -Wimplicit-interface -O3 -march=native'
+f2py_link_flags = '-lgomp'
+
+def f2py_setup():
+    if os.name == 'nt':
+        os.add_dll_directory(f2py_mingw_path)
+
+def f2py_compile(f2py_module_name, f2py_source_files):
+    source_files_string = ' '.join(f2py_source_files)
+    print(f'Running: f2py -m {f2py_module_name} -c {source_files_string} --f90flags="{f2py_opt_flags}" {f2py_link_flags}')
+    return os.system(f'f2py -m {f2py_module_name} -c {source_files_string} --f90flags="{f2py_opt_flags}" {f2py_link_flags}')

f2py/image_processing.py

@@ -1,5 +1,8 @@
-import sys
-sys.path.append("./build") 
+
+import mf_mixed
+
+mf_mixed.f2py_setup()
+mf_mixed.f2py_compile('fortmod', ['simple.f90'])
 
 from numpy import *
 from fortmod import *

f2py/mf_mixed.py

@@ -1,19 +1,21 @@
-import sys
-sys.path.append("./build") 
+import mf_mixed
+
+mf_mixed.f2py_setup()
+mf_mixed.f2py_compile('fortmod', ['matrix.f90'])
 
 from numpy import *
 from fortmod import *
 
-print matrix_multiply2.__doc__
+print(matrix.matrix_multiply2.__doc__)
 
 A = ones((10,10), 'float32', order='F') * 10.0
 B = ones((10,10), 'float32', order='F') * 20.0
 C = zeros((10,10), 'float32', order='F')
 
 print("id of C before multiply =",id(C))
 
-matrix_multiply2(A, B, C)
+matrix.matrix_multiply2(A, B, C)
 
 print("id of C after multiply =",id(C))
 
-print C
+print(C)

f2py/simple.py

@@ -0,0 +1,7 @@
+program myfirst
+
+    implicit none
+
+    print*, 'Hello, Fortran'
+
+end program myfirst

f2py/test.py

@@ -29,9 +29,11 @@ def __init__(self, x=0.0, y=0.0):
         super().__init__()
         self.x = x
         self.y = y
-
         self.vx = 0.0
-        self.vy = 0.0
+        self.vy = 0.
+        self.stroke_color = [255, 0, 0]
+
+        self.stroke_width = 2
 
     def update(self, dt):
         self.x += self.vx*dt
@@ -41,38 +43,43 @@ def on_draw(self):
         sketch.point(self.x, self.y)
 
 class RoundParticle(Particle):
-    def __init__(self, x=0.0, y=0.0, r=1.0):
+    def __init__(self, x=0.0, y=0.0, r=50.0):
         super().__init__(x, y)
         self.r = r
 
     def on_draw(self):
         sketch.ellipse(self.x, self.y, self.r*2, self.r*2)
 
+class SquareParticle(Particle):
+    def __init__(self, x=0.0, y=0.0, w=50.0, h=50.0):
+        super().__init__(x, y)
+        self.w = w
+        self.h = h
+        self.r = w/2
+
+    def on_draw(self):
+        sketch.rect(self.x, self.y, self.w, self.h)
+
 class BoxBoundary:
     def __init__(self):
         self.xmin = 0.0
         self.xmax = 600.0
         self.ymin = 0.0
         self.ymax = 600.0
 
-    def is_inside(self, p):
-        return (p.x - p.r > self.xmin) and (p.x + p.r < self.xmax) and (p.y - p.r > self.ymin) and (p.y + p.r < self.ymin)
-
     def check(self, p):
-        if not self.is_inside(p):
-            if (p.x - p.r) < self.xmin:
-                p.vx = -p.vx
-                p.x = self.xmin + p.r
-            if (p.x + p.r) > self.xmax:
-                p.vx = -p.vx
-                p.x = self.xmax - p.r
-            if (p.y - p.r) < self.ymin:
-                p.vy = -p.vy
-                p.y = self.ymin + p.r
-            if (p.y + p.r) > self.ymax:
-                p.vy = -p.vy
-                p.y = self.ymax - p.r
-
+        if p.x - p.r < self.xmin:
+            p.x = self.xmin + p.r
+            p.vx = -p.vx
+        elif p.x + p.r > self.xmax:
+            p.x = self.xmax - p.r
+            p.vx = -p.vx
+        if p.y - p.r < self.ymin:
+            p.y = self.ymin + p.r
+            p.vy = -p.vy
+        elif p.y + p.r > self.ymax:
+            p.y = self.ymax - p.r
+            p.vy = -p.vy
 
 class ParticleSketch(Sketch):
 
@@ -84,21 +91,28 @@ def setup(self):
         self.ellipse_mode(self.CENTER)
 
         self.particles = []
+
         self.boundary = BoxBoundary()
 
         for i in range(100):
-            p = RoundParticle(self.random(0,600), self.random(0,600), self.random(30, 70))
+            if self.random(1) < 0.5:
+                p = SquareParticle(self.random(0,600), self.random(0, 600), self.random(30, 70), self.random(30, 70))
+            else:
+                p = RoundParticle(self.random(0,600), self.random(0, 600), self.random(30, 70))
+
             p.vx = self.random(-60.0, 60.0)
             p.vy = self.random(-60.0, 60.0)
             p.fill_color = [self.random(255), self.random(255), self.random(255)]
             p.fill_alpha = self.random(50, 255)
+
             self.particles.append(p)
 
+
     def draw(self):
         self.background(40)
 
         for p in self.particles:
-            p.update(1.0/60.0)
+            p.update(1/60)
             self.boundary.check(p)
             p.draw()