diff --git a/2014/2014_08/ada_server/.gitignore b/2014/2014_08/ada_server/.gitignore new file mode 100644 index 0000000..7c35ed5 --- /dev/null +++ b/2014/2014_08/ada_server/.gitignore @@ -0,0 +1,3 @@ +*.ali +*.o +server diff --git a/2014/2014_08/ada_server/constants.ads b/2014/2014_08/ada_server/constants.ads new file mode 100644 index 0000000..1be98f1 --- /dev/null +++ b/2014/2014_08/ada_server/constants.ads @@ -0,0 +1,7 @@ +package constants is + + cr : Character := Character'Val(13); + lf : Character := Character'Val(10); + newline: String := (cr, lf); + +end constants; diff --git a/2014/2014_08/ada_server/dispatchers.adb b/2014/2014_08/ada_server/dispatchers.adb new file mode 100644 index 0000000..162b68b --- /dev/null +++ b/2014/2014_08/ada_server/dispatchers.adb @@ -0,0 +1,75 @@ +package body dispatchers is + + task body dispatcher is + HANDLERS_COUNT : constant Integer := 2; + MAX_SOCKETS_COUNT : constant Integer := 1000; + UNDISPATCHED_SOCKET_RETRY_TIME : constant Duration := 0.5; + + handlers : array (1 .. HANDLERS_COUNT) of Handler; + + type Socket_Index is range 1..MAX_SOCKETS_COUNT; + package Socket_Vector is new Ada.Containers.Vectors ( + Element_Type => Socket_Type, + Index_Type => Socket_Index); + + function Find_Free_Handler(s: Socket_Type) return Boolean is + begin + for id in handlers'Range loop + select + handlers(id).handle(s); + return True; + else + null; + end select; + end loop; + return False; + end Find_Free_Handler; + + result : Boolean; + undispatched : Socket_Vector.Vector; + + begin + accept start; + + for id in handlers'Range loop + handlers(id).start(id); + end loop; + + loop + select + accept dispatch (s: Socket_Type) do + Put_Line ("dispatch command"); + result := Find_Free_Handler(s); + if not result then + Put_Line("All handlers are busy"); + undispatched.Append(s); + end if; + end dispatch; + or + accept stop do + Put_Line("stop command, stopping dispatcher"); + end stop; + exit; + or + delay UNDISPATCHED_SOCKET_RETRY_TIME; + if Integer(undispatched.Length) > 0 then + declare + s: Socket_Type := undispatched.Element(1); + begin + result := Find_Free_Handler(s); + if result then + undispatched.Delete(1); + Put_Line("Handled undispatched client"); + end if; + end; + end if; + end select; + end loop; + + for id in handlers'Range loop + handlers(id).stop; + end loop; + + + end dispatcher; +end dispatchers; diff --git a/2014/2014_08/ada_server/dispatchers.ads b/2014/2014_08/ada_server/dispatchers.ads new file mode 100644 index 0000000..551d3cf --- /dev/null +++ b/2014/2014_08/ada_server/dispatchers.ads @@ -0,0 +1,17 @@ +with Ada.Strings.Unbounded; use Ada.Strings.Unbounded; +with Ada.Text_IO; use Ada.Text_IO; +with Gnat.Sockets; use Gnat.Sockets; +with Ada.Calendar; use Ada.Calendar; +with Ada.Containers.Vectors; +with Handlers; use Handlers; + +package dispatchers is + task dispatcher is + entry start; + -- signal to start + entry dispatch(s: Socket_Type); + -- we got accepted socket from the listener + entry stop; + -- signal to stop + end dispatcher; +end dispatchers; diff --git a/2014/2014_08/ada_server/handlers.adb b/2014/2014_08/ada_server/handlers.adb new file mode 100644 index 0000000..ddc8f77 --- /dev/null +++ b/2014/2014_08/ada_server/handlers.adb @@ -0,0 +1,247 @@ +package body handlers is + task body handler is + socket : Socket_Type; + sel : access Gnat.Sockets.Selector_Type := new + Gnat.Sockets.Selector_Type; + working : Boolean := True; + subtype Socket_Index is Integer range 1 .. 1000; + root_path : constant String := "/var/www"; + + package Socket_Vector is new Ada.Containers.Vectors ( + Element_Type => Socket_Type, + Index_Type => Socket_Index + ); + + v : Socket_Vector.Vector; + + task watcher is + entry start; + end watcher; + + task body watcher is + R : Socket_Set_Type; + W : Socket_Set_Type; + E : Socket_Set_Type; + status : Selector_Status; + active : Socket_Type; + unparsed: Unbounded_String; + begin + accept start; + Put_Line("start watcher"); + Create_Selector (sel.all); + loop + Put_Line("Watch loop start"); + Empty(R); + Empty(W); + Empty(E); + + Put("There is "); + Put(v.Last_Index, Width => 0); + Put_Line(" sockets to watch"); + + Put_Line("setting to read sockets"); + + for id in 1..v.Last_Index loop + set(R, v.Element(id)); + end loop; + + Put_Line("start of check_selector"); + + Check_Selector( + sel.all, + R, + W, + E, + status, + Gnat.Sockets.Forever); + + Put_Line("end of check_selector"); + + -- wychodzimy jedynie wtedy, gdy zmienna working + -- jest False. Inaczej to jedynie przeładowanie + -- socketów np. dodano nam socket do obserwacji + exit when status = Gnat.Sockets.Aborted + and not working; + + if status /= Gnat.Sockets.Aborted then + Put_Line("Sending command to client handler"); + Get(R, active); + Handler.command(active); + end if; + + Put_Line("Watch loop end"); + + end loop; + Close_Selector (sel.all); + Put_Line("exit watcher"); + end watcher; + + procedure Handle_File(s: Socket_Type; + path: String; + finished: out Boolean) is + + sz : Natural := Natural(Size(path)); + subtype Content_Type is String (1 .. sz); + package dio is new Ada.Direct_IO(Content_Type); + File : dio.File_Type; + contents : Content_Type; + channel : Stream_Access := Stream(s); + begin + Put_Line("Handle_File start"); + Put_Line("Path: " & path); + Put_Line("File size: " & Natural'Image(sz)); + + dio.Open( + File => File, + Mode => dio.In_File, + Name => path); + + Put_Line("File opened"); + + String'Write(channel, "HTTP/1.0 200 OK" & newline); + String'Write(channel, "Server: Ada Server 0.1" & + newline); + + Put_Line("Root path: " & root_path & path); + Put_Line("File size: " & Natural'Image(sz)); + + dio.Read (File, Item => contents); + Put_Line("Contents: " & contents & newline); + String'Write(channel, newline); + String'Write(channel, contents); + + dio.Close(File); + Put_Line("Closed file"); + + Close_Socket(s); + Put_Line("Closed socket"); + finished := True; + end Handle_File; + + procedure Handle_Socket(s: Socket_Type; finished: out Boolean) is + subtype Line_Type is String (1 .. 4096); + channel : Stream_Access := Stream (s); + line : Line_Type; + elem_line : Stream_Element_Array (1 .. 4096); + last : Stream_Element_Offset; + + package af renames Ada.Strings.Fixed; + + function Convert is new Ada.Unchecked_Conversion ( + Source => Stream_Element_Array, + Target => Line_Type + ); + begin + Put_Line("Got command to execute"); + channel := Stream (s); + Receive_Socket(s, elem_line, last); + line := Convert(elem_line); + Put("Odczytano: "); + Put(Integer(last), Width => 0); + Put_Line(" znakow"); + if last >= 4 and line (1 .. 4) = "GET " then + Put_Line("GET Command: " & line); + + declare + EOL : Integer := af.Index(line, newline); + HTTP_Version_Offset : Integer := af.Index(line, "HTTP/1."); + path : String := line (5..HTTP_Version_Offset - 2); + begin + if Exists (root_path & path) then + Put_Line("Kind: " & File_Kind'Image(Kind(root_path & path))); + if Kind(root_path & path) = Ordinary_File then + Put_Line("Full path: '" & root_path & path & "'"); + handle_file(s, + root_path & path, + finished); + return; + else + -- pewnie katalog + if Exists(root_path & path & "index.html") then + handle_file(s, + root_path & path & + "index.html", + finished); + + return; + end if; + + if Exists(root_path & path & "index.htm") then + handle_file(s, + root_path & path & + "index.htm", + finished); + + return; + end if; + + String'Write(channel, "HTTP/1.0 200 OK" & newline); + String'Write(channel, "Server: Ada Server 0.1" & + newline & newline); + String'Write(channel, "" & + "directory" & newline); + end if; + else + String'Write (channel, "HTTP/1.0 404 NOT FOUND" & + newline); + String'Write (channel, "Server: Ada Server 0.1" & + newline & newline); + String'Write (channel, "" & + "404 Error"); + end if; + end; + + else + String'Write(channel, "Unknown" & newline); + end if; + Put_Line("Closing socket"); + Close_Socket(s); + finished := True; + end Handle_Socket; + + idx : Integer; + + begin + accept start(id: Integer) do + idx := id; + end start; + + Watcher.Start; + + Put("start client handler with id: "); + Put(idx, Width => 0); + New_Line; + loop + select + accept handle(s: Socket_Type) do + socket := s; + v.Append(s); + end handle; + Abort_Selector(sel.all); + Put_Line("Socket consumed by client handler"); + or + accept command(s: Socket_Type) do + declare + finished: Boolean; + begin + Handle_Socket(s, finished); + if finished then + v.Delete(v.Find_Index(s)); + end if; + end; + end command; + or + accept stop do + Put_Line("Got command to stop"); + working := False; + Abort_Selector(sel.all); + end stop; + exit; + or + terminate; + end select; + end loop; + Put_Line("end client handler"); + end handler; + +end handlers; diff --git a/2014/2014_08/ada_server/handlers.ads b/2014/2014_08/ada_server/handlers.ads new file mode 100644 index 0000000..d3d7200 --- /dev/null +++ b/2014/2014_08/ada_server/handlers.ads @@ -0,0 +1,37 @@ +with Ada.Containers.Vectors; +with Ada.Integer_Text_IO; use Ada.Integer_Text_IO; +with Ada.Strings.Unbounded; use Ada.Strings.Unbounded; +with Ada.Text_IO; use Ada.Text_IO; +with GNAT.Sockets; use GNAT.Sockets; +with Constants; use Constants; +with Ada.Calendar; use Ada.Calendar; +with Ada.Streams; use Ada.Streams; +with Ada.Unchecked_Conversion; +with Ada.Strings.Fixed; +with Ada.Directories; use Ada.Directories; +with Ada.Direct_IO; +with Ada.Exceptions; use Ada.Exceptions; +with Ada.IO_Exceptions; use Ada.IO_Exceptions; + +package handlers is + + type Client_Type is record + socket : Socket_Type; + read : Unbounded_String; + unread : Unbounded_String; + last_active : Time; + end record; + + type Client_Index_Type is range 1..100; + package Client_Vector_Type is new Ada.Containers.Vectors ( + Element_Type => Client_Type, + Index_Type => Client_Index_Type); + + task type handler is + entry start(id: Integer); + entry handle(s: Socket_Type); + entry command(s: Socket_Type); + entry stop; + end handler; + +end handlers; diff --git a/2014/2014_08/ada_server/server.adb b/2014/2014_08/ada_server/server.adb new file mode 100644 index 0000000..7c95dab --- /dev/null +++ b/2014/2014_08/ada_server/server.adb @@ -0,0 +1,67 @@ +with Ada.Text_IO; use Ada.Text_IO; +with Ada.Integer_Text_IO; use Ada.Integer_Text_IO; +with GNAT.Sockets; use GNAT.Sockets; +with signals; use signals; +with Ada.Containers.Vectors; +with Handlers; use Handlers; +with Constants; use Constants; +with Dispatchers; use Dispatchers; + +procedure server is + + HOST : constant String := "localhost"; + PORT : constant Port_Type := 6666; + work : Boolean := True; + address : Sock_Addr_Type; + server : Socket_Type; + socket : Socket_Type; + h2 : Sigint_Handler; + status : Gnat.Sockets.Selector_Status; + selector: signals.selector_Access; + req : Request_Type(Non_Blocking_IO); +begin + Put_Line("start of server"); + Dispatcher.Start; + selector := new Gnat.Sockets.Selector_Type; + h2.selector(selector); + Gnat.Sockets.Create_Selector(selector.all); + address.Addr := Addresses (Get_Host_By_Name(host), 1); + address.Port := PORT; + + Create_Socket(server); + Set_Socket_Option ( + server, + Socket_Level, + (Reuse_Address, True)); + + Bind_Socket (server, address); + Put_Line("Binded"); + Listen_Socket(server); + Put_Line("Listening"); + + loop + Put_Line("Waiting to accept"); + Accept_Socket ( + Server => server, + Socket => socket, + Address => address, + Timeout => GNAT.Sockets.Forever, + Selector => selector, + Status => status); + exit when status = Gnat.Sockets.Aborted; + Put_Line("Accepted"); + + Control_Socket(socket, req); + + Dispatcher.dispatch(socket); + Put_Line("Handled by listener"); + + end loop; + + Gnat.Sockets.Close_Selector(selector.all); + Close_Socket (server); + + Dispatcher.Stop; + + Put_Line("end of server"); +end server; diff --git a/2014/2014_08/ada_server/signals.adb b/2014/2014_08/ada_server/signals.adb new file mode 100644 index 0000000..d172195 --- /dev/null +++ b/2014/2014_08/ada_server/signals.adb @@ -0,0 +1,16 @@ +package body signals is + protected body Sigint_Handler is + + procedure selector(s: Selector_Access) is + begin + sel := s; + end; + + procedure Handle is + begin + Call_Count := Call_Count + 1; + Put_Line("SIGINT handled"); + Gnat.Sockets.Abort_Selector(sel.all); + end Handle; + end Sigint_Handler; +end signals; diff --git a/2014/2014_08/ada_server/signals.ads b/2014/2014_08/ada_server/signals.ads new file mode 100644 index 0000000..32d68d3 --- /dev/null +++ b/2014/2014_08/ada_server/signals.ads @@ -0,0 +1,23 @@ +with Ada.Interrupts; use Ada.Interrupts; +with Ada.Interrupts.Names; use Ada.Interrupts.Names; +with Ada.Text_IO; use Ada.Text_IO; +with Gnat.Sockets; use Gnat.Sockets; + +package signals is + + pragma Unreserve_All_Interrupts; + + type Selector_Access is access all Gnat.Sockets.Selector_Type; + + protected type Sigint_Handler is + procedure selector(s: Selector_Access); + procedure Handle; + + pragma Interrupt_Handler(Handle); + pragma Attach_Handler(Handle, Sigint); + private + Call_Count : Natural := 0; + sel : Selector_Access; + end Sigint_Handler; + +end signals; diff --git a/2015/2015_03/shoe-soles/README.md b/2015/2015_03/shoe-soles/README.md new file mode 100644 index 0000000..85e492b --- /dev/null +++ b/2015/2015_03/shoe-soles/README.md @@ -0,0 +1 @@ +Projekt rozpoznawania obuwia. Bardzo wczesna faza. diff --git a/2015/2015_03/shoe-soles/dect.py b/2015/2015_03/shoe-soles/dect.py new file mode 100644 index 0000000..72904c2 --- /dev/null +++ b/2015/2015_03/shoe-soles/dect.py @@ -0,0 +1,75 @@ +from scipy import ndimage +from skimage import filters,feature +import sys +import matplotlib.pyplot as plt +import math +import numpy as np +import pdb + +def calculate_edges(img): + # TODO fix the performance - added ad hoc + (n,m) = img.shape + edge = np.zeros((n,m), dtype=np.bool) +# pdb.set_trace() + for i in range(0, n): + for j in range(0, m): + edge_point = False + for k in range(-1, 2): + if edge_point: + break + for l in range(-1,2): + if edge_point: + break + + y = (i+k)%n + x = (j+l)%m + + if img[y,x] != img[i,j]: + edge_point = True + break + edge[i,j] = edge_point + return edge + +def calc_shape_corners(img_mask2, range_ext=8): + (my,mx) = ndimage.center_of_mass(img_mask2) +# edge = feature.canny(img_mask2) + edge = calculate_edges(img_mask2) + + (n,m) = img_mask2.shape + + print("Center of mass {}x{}".format(mx,my)) + distances = np.zeros((n,m), dtype=np.float) + distances[:,:] = 0 + for i in range(0, n): + for j in range(0, m): + if edge[i,j]: + dx = j - mx + dy = i - my + d = math.sqrt(dx*dx+dy*dy) + distances[i,j] = d + + max_distances = ndimage.maximum_filter(distances, range_ext) + diff = np.abs(max_distances - distances) + corners = edge * (diff < 1e-3) + return corners + +if __name__ == "__main__": + img = ndimage.imread(sys.argv[1]) + img_gray = 0.21 * img[:,:,0] + 0.72 * img[:,:,1] + 0.07 * img[:,:,2] + otsu_lvl = filters.threshold_otsu(img_gray) + img_mask = img_gray <= otsu_lvl + img_mask2 = ndimage.median_filter(img_mask, 3) + (n,m) = img_mask2.shape + + + corners = calc_shape_corners(img_mask2) + + mat_x = np.tile(np.arange(0,m),n).reshape((n,m)) + mat_y = np.tile(np.arange(0,n),m).reshape((m,n)).T + + xs = mat_x[corners] + ys = mat_y[corners] + + plt.imshow(img_gray, cmap="Greys_r") + plt.plot(xs, ys, "r+") + plt.show() diff --git a/2015/2015_03/shoe-soles/filtering-print.py b/2015/2015_03/shoe-soles/filtering-print.py new file mode 100644 index 0000000..5db0234 --- /dev/null +++ b/2015/2015_03/shoe-soles/filtering-print.py @@ -0,0 +1,8 @@ +from scipy import ndimage +import matplotlib.pyplot as plt + +img = ndimage.imread("/home/tpolgrabia/Pobrane/CSFID/tracks_cropped/00003.jpg") +img_filtered = ndimage.median_filter(img, 5) + +plt.imshow(img_filtered, cmap="Greys_r") +plt.show() diff --git a/2015/2015_03/shoe-soles/image_boundaries.py b/2015/2015_03/shoe-soles/image_boundaries.py new file mode 100644 index 0000000..826fb3f --- /dev/null +++ b/2015/2015_03/shoe-soles/image_boundaries.py @@ -0,0 +1,50 @@ +# -*- coding: UTF-8 -*- +import numpy as np + +def find_minx(mat,val = 1): + mat2 = mat.reshape((-1,1), order='F') + (n,m) = mat.shape + n2 = mat2.shape[0] + for (el,idx) in zip(np.nditer(mat2),range(0,n2)): + if el == val: + return idx / n + + return None + +def find_miny(mat,val = 1): + mat2 = mat.reshape((1,-1)) + (n,m) = mat.shape + n2 = mat2.shape[1] + for (el,idx) in zip(np.nditer(mat2),range(0,n2)): + if el == val: + return idx / m + + return None + +def find_maxx(mat,val = 1): + mat2 = np.rot90(np.rot90(mat)).reshape((-1,1), order='F') + (n,m) = mat.shape + n2 = mat2.shape[0] + for (el,idx) in zip(np.nditer(mat2),range(0,n2)): + if el == val: + return m - 1 - idx / n + return None + +def find_maxy(mat,val = 1): + mat_rotated = np.rot90(np.rot90(mat)) + mat2 = mat_rotated.reshape((1,-1)) + (n,m) = mat.shape + n2 = mat2.shape[1] + for (el,idx) in zip(np.nditer(mat2),range(0,n2)): + if el == val: + return n - 1 - idx / m + + return None + +def calc_image_crop(img, val = 255): + minx = find_minx(img,val) + maxx = find_maxx(img,val) + miny = find_miny(img,val) + maxy = find_maxy(img,val) + return (minx,maxx,miny,maxy) + diff --git a/2015/2015_03/shoe-soles/image_boundaries2.py b/2015/2015_03/shoe-soles/image_boundaries2.py new file mode 100644 index 0000000..840d24d --- /dev/null +++ b/2015/2015_03/shoe-soles/image_boundaries2.py @@ -0,0 +1,66 @@ +# -*- coding: UTF-8 -*- +import numpy as np +from scipy import ndimage + +def find_miny(mat,val = 1): + (n,m) = mat.shape + for i in range(0,n): + for j in range(0,m): + if mat[i,j] == val: + return i + + return None + +def find_minx(mat,val = 1): + (n,m) = mat.shape + for j in range(0,m): + for i in range(0,n): + if mat[i,j] == val: + return j + + return None + + +def find_maxy(mat,val = 1): + (n,m) = mat.shape + for i in range(n-1,-1,-1): + for j in range(m-1,-1,-1): + if mat[i,j] == val: + return i + + return None + +def find_maxx(mat,val = 1): + (n,m) = mat.shape + for j in range(m-1,-1,-1): + for i in range(n-1,-1,-1): + if mat[i,j] == val: + return j + + return None + + +def calc_image_crop(img, val = 255): + minx = find_minx(img,val) + maxx = find_maxx(img,val) + miny = find_miny(img,val) + maxy = find_maxy(img,val) + return (minx,maxx,miny,maxy) + +def crop_image(img,val = 255): + (minx,maxx,miny,maxy) = calc_image_crop(img,val) + return img[miny:maxy,minx:maxx] + +def crop_image_prefilled(img,val = 255): + (minx,maxx,miny,maxy) = calc_image_crop(img,val) + img_filled = ndimage.binary_fill_holes(img) + return img_filled[miny:maxy,minx:maxx] + +def dimension_rate(img): + (n,m) = img.shape + return float(m)/float(n) + +def image_to_mask(img,threshold = 128): + m = img >= threshold + return 1*m + diff --git a/2015/2015_03/shoe-soles/measures.py b/2015/2015_03/shoe-soles/measures.py new file mode 100644 index 0000000..fcdb18a --- /dev/null +++ b/2015/2015_03/shoe-soles/measures.py @@ -0,0 +1,85 @@ +#!/usr/bin/python + +from scipy import ndimage +import sys +import math +import matplotlib.pyplot as plt +import numpy as np + +print("Name: {}".format(__name__)) +if __name__ != "__main__": + print("not runned") + sys.exit(0) + +if len(sys.argv) <= 1: + print("Too small arguments") + sys.exit(1) + +path = sys.argv[1] + +img = ndimage.imread(path) +img_gray = img[:,:,0] + +(n,m) = img_gray.shape +mask = img_gray >= 128 + +sx = 0.0 +sy = 0.0 +c = 0 + +for i in range(0, n): + for j in range(0, m): + if mask[i,j]: + sx += j + sy += i + c += 1 + +ax = sx / c +ay = sy / c + +print("Position: {} x {}".format(ax,ay)) + +grav_measure = 0.0 + +for i in range(0, n): + for j in range(0, m): + if mask[i,j]: + grav_measure += (j - ax) * (j - ax) + (i - ay) * (i - ay) + +grav_measure2 = math.sqrt(grav_measure / c) + +print("Grav measure: {}".format(grav_measure)) + +# vector of offset + +ox = 0.0 +oy = 0.0 + +for i in range(0, n): + for j in range(0, m): + if mask[i,j]: + ox += math.pow(j - ax, 1.0) + oy += math.pow(i - ay, 1.0) + +print("Vector of offset (measure of beaing simetrical): {}x{}".format(ox,oy)) + +img_filled = ndimage.binary_fill_holes(img_gray) +s = np.sum(img_filled) +unit = max(n,m) +c2 = c / float(unit) +s2 = s / (float(unit)*float(unit)) +print("l: {}, s: {}".format(c2,s2)) +# plt.imshow(img_filled, cmap="Greys_r") +# plt.show() + +circularity1 = 2.0 * math.sqrt(s2 / np.pi) +print("Circularity1: {}".format(circularity1)) + +circularity2 = c2 / math.pi +print("Circularity2: {}".format(circularity2)) + +w_measure = c2 / (2.0 * math.sqrt(np.pi * s2)) - 1.0 +print("W measure: {}".format(w_measure)) + +w9_measure = 2.0 * math.sqrt(np.pi * s2) / c2 +print("W9 measure (malinkowskiej) {}".format(w9_measure)) diff --git a/2015/2015_03/shoe-soles/old/ipython_log.py b/2015/2015_03/shoe-soles/old/ipython_log.py new file mode 100644 index 0000000..fb4205c --- /dev/null +++ b/2015/2015_03/shoe-soles/old/ipython_log.py @@ -0,0 +1,437 @@ +# IPython log file + +from dect import calc_shape_corners +calc_shape_corners(img) +from scipy import ndimage +img = ndimage.imread("/home/tpolgrabia/Pobrane/CSFID/tracks_cropped/00003.jpg") +img.shape +#[Out]# (570, 193) +import matplotlib.pyplot as plt +plt.imshow(img, cmap="Greys_r") +#[Out]# +plt.show() +from skimage import filters +otsu_lvl = filters.threshold_otsu(img) +otsu_lvl +#[Out]# 139 +plt.imshow(img, cmap="Greys_r") +#[Out]# +plt.show() +plt.imshow(img >= otsu_lvl, cmap="Greys_r") +#[Out]# +plt.show() +m = img <= otsu_lvl +plt.imshow(m, cmap="Greys_r") +#[Out]# +plt.show() +img_label, nb_labels = ndimage.label(m) +nb_labels +#[Out]# 571 +plt.imshow(nb_labels) +plt.imshow(img_label)) +plt.imshow(img_label) +#[Out]# +plt.show() +sizes = ndimage.sum(m, img_label, range(0, nb_labels+1)) +sizes +#[Out]# array([ 0.00000000e+00, 1.00000000e+00, 4.50000000e+01, +#[Out]# 1.70000000e+01, 6.20000000e+01, 1.50000000e+01, +#[Out]# 5.00000000e+00, 2.00000000e+00, 2.00000000e+00, +#[Out]# 3.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 3.00000000e+00, 9.00000000e+00, +#[Out]# 2.90000000e+01, 1.00000000e+00, 5.00000000e+00, +#[Out]# 2.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 4.00000000e+00, 5.00000000e+00, 1.00000000e+00, +#[Out]# 3.00000000e+00, 2.40000000e+01, 6.00000000e+00, +#[Out]# 5.00000000e+00, 3.10000000e+01, 1.00000000e+00, +#[Out]# 5.30000000e+01, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 2.00000000e+00, +#[Out]# 2.00000000e+00, 2.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 2.00000000e+00, 6.00000000e+00, +#[Out]# 3.80000000e+01, 1.00000000e+00, 3.00000000e+00, +#[Out]# 3.00000000e+00, 1.20000000e+01, 1.10000000e+01, +#[Out]# 1.00000000e+01, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 5.60000000e+01, 3.00000000e+00, 1.00000000e+00, +#[Out]# 2.00000000e+00, 1.00000000e+00, 2.50000000e+01, +#[Out]# 1.00000000e+00, 3.00000000e+00, 4.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 3.00000000e+00, 3.00000000e+00, +#[Out]# 3.00000000e+00, 2.00000000e+00, 1.00000000e+00, +#[Out]# 3.00000000e+00, 3.00000000e+00, 2.20000000e+01, +#[Out]# 1.00000000e+00, 3.00000000e+00, 3.00000000e+00, +#[Out]# 1.00000000e+00, 8.00000000e+00, 1.00000000e+01, +#[Out]# 1.00000000e+00, 2.00000000e+00, 1.46000000e+02, +#[Out]# 2.00000000e+00, 3.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.40000000e+01, 1.00000000e+00, +#[Out]# 1.00000000e+00, 4.00000000e+00, 1.15550000e+04, +#[Out]# 1.00000000e+00, 1.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 5.80000000e+01, 2.00000000e+00, +#[Out]# 1.00000000e+00, 2.00000000e+00, 8.00000000e+00, +#[Out]# 1.00000000e+01, 2.00000000e+00, 3.00000000e+00, +#[Out]# 1.00000000e+00, 3.00000000e+00, 5.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 5.00000000e+00, +#[Out]# 1.00000000e+00, 2.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 2.00000000e+00, 3.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 3.00000000e+00, +#[Out]# 2.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 5.40000000e+01, 5.00000000e+00, +#[Out]# 8.00000000e+00, 1.90000000e+01, 1.00000000e+00, +#[Out]# 1.00000000e+00, 2.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+01, 3.00000000e+00, 2.00000000e+00, +#[Out]# 2.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 3.30000000e+01, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 4.00000000e+00, 1.00000000e+00, +#[Out]# 4.00000000e+00, 1.30000000e+01, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 2.00000000e+00, +#[Out]# 2.00000000e+00, 2.00000000e+00, 1.00000000e+00, +#[Out]# 2.00000000e+00, 1.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 2.00000000e+00, 3.00000000e+00, +#[Out]# 1.40000000e+01, 2.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 2.00000000e+00, 3.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 3.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 5.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 5.00000000e+00, 4.00000000e+00, +#[Out]# 4.00000000e+00, 2.20000000e+01, 2.00000000e+00, +#[Out]# 3.00000000e+00, 1.00000000e+00, 3.00000000e+00, +#[Out]# 3.00000000e+00, 2.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 2.00000000e+00, 1.00000000e+00, 4.00000000e+00, +#[Out]# 4.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 5.00000000e+00, 1.10000000e+01, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 3.00000000e+00, 1.00000000e+00, 3.00000000e+00, +#[Out]# 3.00000000e+00, 4.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 2.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 4.00000000e+00, 5.00000000e+00, +#[Out]# 3.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 2.00000000e+00, 1.20000000e+01, 1.01400000e+03, +#[Out]# 4.00000000e+00, 1.00000000e+00, 3.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 5.00000000e+00, +#[Out]# 4.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 2.00000000e+00, 1.00000000e+00, +#[Out]# 1.90000000e+01, 2.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+01, 3.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 2.00000000e+00, 4.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 2.00000000e+00, +#[Out]# 1.10000000e+01, 5.00000000e+00, 1.00000000e+00, +#[Out]# 7.00000000e+00, 2.00000000e+00, 1.00000000e+00, +#[Out]# 3.00000000e+00, 1.00000000e+00, 2.00000000e+00, +#[Out]# 6.40000000e+01, 1.00000000e+00, 1.00000000e+00, +#[Out]# 2.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 2.00000000e+00, 1.00000000e+00, 3.00000000e+00, +#[Out]# 4.00000000e+00, 5.80000000e+01, 9.00000000e+00, +#[Out]# 9.30000000e+01, 8.10000000e+01, 6.00000000e+00, +#[Out]# 2.00000000e+00, 1.00000000e+00, 2.91000000e+02, +#[Out]# 3.00000000e+00, 1.00000000e+00, 4.00000000e+00, +#[Out]# 1.30000000e+01, 1.35000000e+02, 3.10000000e+01, +#[Out]# 4.00000000e+00, 2.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 3.00000000e+01, 1.00000000e+01, +#[Out]# 5.20000000e+01, 2.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 4.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 2.62000000e+02, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 4.00000000e+00, 2.28000000e+03, 1.00000000e+00, +#[Out]# 1.39000000e+02, 2.00000000e+00, 2.00000000e+00, +#[Out]# 2.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 3.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 3.00000000e+00, 1.00000000e+00, 2.75000000e+02, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 8.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.13000000e+02, 2.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 2.00000000e+00, 2.27000000e+02, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.40000000e+01, 2.00000000e+00, +#[Out]# 2.26000000e+02, 2.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 3.54000000e+02, 1.00000000e+00, 1.10000000e+01, +#[Out]# 1.00000000e+00, 1.00000000e+01, 1.00000000e+00, +#[Out]# 1.00000000e+00, 9.00000000e+00, 1.00000000e+00, +#[Out]# 1.40000000e+01, 1.00000000e+00, 2.00000000e+00, +#[Out]# 1.11000000e+02, 2.00000000e+00, 1.00000000e+00, +#[Out]# 5.00000000e+00, 6.00000000e+00, 1.89000000e+02, +#[Out]# 1.00000000e+00, 8.00000000e+00, 1.00000000e+00, +#[Out]# 7.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 7.90000000e+01, 3.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 5.00000000e+00, 1.00000000e+00, +#[Out]# 8.00000000e+00, 2.00000000e+00, 1.94000000e+02, +#[Out]# 1.20000000e+01, 1.00000000e+00, 4.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 2.00000000e+00, +#[Out]# 8.70000000e+01, 1.55000000e+02, 1.68400000e+03, +#[Out]# 3.00000000e+00, 1.30000000e+01, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 2.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 3.00000000e+00, +#[Out]# 6.00000000e+00, 1.87900000e+03, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 6.00000000e+00, +#[Out]# 3.15000000e+02, 1.59000000e+02, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 5.00000000e+00, +#[Out]# 6.80000000e+01, 2.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 2.00000000e+00, +#[Out]# 1.59000000e+02, 3.95000000e+02, 6.00000000e+00, +#[Out]# 5.00000000e+00, 6.00000000e+00, 6.80000000e+01, +#[Out]# 4.00000000e+01, 5.10000000e+01, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 3.69000000e+02, 1.17000000e+02, 1.00000000e+00, +#[Out]# 2.00000000e+00, 3.80000000e+01, 5.61000000e+02, +#[Out]# 3.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 9.00000000e+00, +#[Out]# 1.00000000e+00, 2.00000000e+00, 1.00000000e+00, +#[Out]# 2.89000000e+02, 1.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 2.00000000e+00, +#[Out]# 1.70000000e+01, 2.00000000e+00, 1.00000000e+00, +#[Out]# 4.00000000e+00, 3.00000000e+00, 4.00000000e+00, +#[Out]# 3.00000000e+00, 1.00000000e+00, 4.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00, 3.00000000e+00, +#[Out]# 3.00000000e+00, 2.00000000e+00, 6.00000000e+00, +#[Out]# 1.00000000e+00, 3.00000000e+00, 1.00000000e+00, +#[Out]# 2.00000000e+00, 2.00000000e+00, 3.10000000e+01, +#[Out]# 7.00000000e+00, 2.00000000e+00, 9.00000000e+00, +#[Out]# 6.10000000e+01, 1.00000000e+00, 6.00000000e+00, +#[Out]# 3.00000000e+00, 1.00000000e+00, 7.00000000e+00, +#[Out]# 2.00000000e+00, 1.00000000e+00, 6.00000000e+00, +#[Out]# 2.40000000e+01, 5.00000000e+00, 2.00000000e+00, +#[Out]# 1.00000000e+00, 2.00000000e+00, 9.00000000e+00, +#[Out]# 9.40000000e+01, 1.00000000e+00, 1.00000000e+00, +#[Out]# 2.00000000e+00, 1.00000000e+00, 8.00000000e+00, +#[Out]# 2.60000000e+01, 1.00000000e+00, 2.50000000e+01, +#[Out]# 1.00000000e+00, 1.00000000e+00, 4.00000000e+00, +#[Out]# 1.80000000e+01, 2.00000000e+00, 7.00000000e+00, +#[Out]# 3.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.40000000e+01, 1.40000000e+01, 4.90000000e+01, +#[Out]# 1.00000000e+00, 2.00000000e+00, 4.00000000e+00, +#[Out]# 1.00000000e+00, 2.50000000e+01, 4.00000000e+00, +#[Out]# 4.00000000e+00, 7.00000000e+00, 2.30000000e+01, +#[Out]# 3.00000000e+00, 5.00000000e+00, 8.00000000e+00, +#[Out]# 6.60000000e+01, 8.00000000e+00, 1.00000000e+01, +#[Out]# 3.00000000e+00, 6.00000000e+00, 2.00000000e+00, +#[Out]# 4.00000000e+00, 1.60000000e+01, 1.00000000e+00, +#[Out]# 6.70000000e+01, 7.00000000e+00, 1.50000000e+02, +#[Out]# 4.10000000e+01, 1.00000000e+00, 2.20000000e+01, +#[Out]# 3.00000000e+00, 1.00000000e+00, 1.00000000e+00, +#[Out]# 1.00000000e+00, 1.00000000e+00]) +size_min_limit = 50 +remove_segments = sizes < size_min_limit +remove_segments[img_label] +#[Out]# array([[ True, True, True, ..., True, True, True], +#[Out]# [ True, True, True, ..., True, True, True], +#[Out]# [ True, True, True, ..., True, True, True], +#[Out]# ..., +#[Out]# [False, False, False, ..., True, True, True], +#[Out]# [False, False, False, ..., True, True, True], +#[Out]# [False, False, False, ..., True, True, True]], dtype=bool) +remove_mask = remove_segments[img_label] +img_label[remove_mask] = 0 +plt.imshow(img_label) +#[Out]# +plt.show() +import numpy as np +get_ipython().magic(u'pinfo np.unique') +labels = np.unique(img_label) +labels +#[Out]# array([ 0, 4, 30, 54, 83, 92, 97, 127, 239, 276, 289, 291, 292, +#[Out]# 296, 301, 309, 319, 325, 327, 338, 351, 356, 363, 369, 381, 386, +#[Out]# 393, 401, 408, 409, 410, 430, 435, 436, 441, 447, 448, 452, 454, +#[Out]# 459, 460, 464, 474, 504, 519, 552, 561, 563], dtype=int32) +get_ipython().magic(u'pinfo np.searchsorted') +img_label = np.searchsorted(labels, img_label) +nb_labels = len(labels) +plt.imshow(img_label) +#[Out]# +plt.show() +ndimage.center_of_mass(img_label, labels, nb_labels+1) +ndimage.center_of_mass(img_label, labels, nb_labels+1) +get_ipython().magic(u'pinfo ndimage.center_of_mass') +ndimage.center_of_mass(m, img_label, labels) +#[Out]# [(274.99289099526067, 99.926540284360186), (60.863013698630134, 45.006849315068493), (413.25593824228031, 114.4667458432304), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan)] +centers = ndimage.center_of_mass(m, img_label, labels) +centers +#[Out]# [(274.99289099526067, 99.926540284360186), (60.863013698630134, 45.006849315068493), (413.25593824228031, 114.4667458432304), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan), (nan, nan)] +centers[0] +#[Out]# (274.99289099526067, 99.926540284360186) +centers[1] +#[Out]# (60.863013698630134, 45.006849315068493) +centers[2] +#[Out]# (413.25593824228031, 114.4667458432304) +centers[3] +#[Out]# (nan, nan) +centers[4] +#[Out]# (nan, nan) +centers[5] +#[Out]# (nan, nan) +labels] +labels +#[Out]# array([ 0, 4, 30, 54, 83, 92, 97, 127, 239, 276, 289, 291, 292, +#[Out]# 296, 301, 309, 319, 325, 327, 338, 351, 356, 363, 369, 381, 386, +#[Out]# 393, 401, 408, 409, 410, 430, 435, 436, 441, 447, 448, 452, 454, +#[Out]# 459, 460, 464, 474, 504, 519, 552, 561, 563], dtype=int32) +labels = np.unique(img_label) +labels +#[Out]# array([ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, +#[Out]# 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, +#[Out]# 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47]) +centers = ndimage.center_of_mass(m, img_label, labels) +centers +#[Out]# [(274.99289099526067, 99.926540284360186), (7.112903225806452, 30.677419354838708), (28.69811320754717, 103.75471698113208), (42.089285714285715, 112.85714285714286), (60.863013698630134, 45.006849315068493), (154.04898312418865, 112.49026395499784), (71.5, 32.413793103448278), (89.5, 60.388888888888886), (184.49112426035504, 34.280078895463511), (222.828125, 47.25), (238.68965517241378, 186.10344827586206), (241.40860215053763, 52.473118279569896), (241.8641975308642, 121.19753086419753), (250.85910652920961, 163.26116838487974), (259.62222222222221, 56.814814814814817), (259.84615384615387, 121.51923076923077), (269.1641221374046, 156.87022900763358), (445.33201754385965, 185.51271929824563), (278.58992805755395, 67.2158273381295), (286.66181818181821, 153.59272727272727), (297.10619469026551, 74.13274336283186), (305.29955947136563, 151.12775330396477), (318.63274336283183, 77.853982300884951), (327.37853107344631, 144.00564971751413), (342.48648648648651, 78.090090090090087), (350.68783068783068, 140.42328042328043), (360.32911392405066, 78.455696202531641), (368.29381443298968, 140.70103092783506), (377.85057471264366, 74.402298850574709), (386.11612903225807, 141.27741935483871), (413.25593824228031, 114.4667458432304), (495.27833954230977, 6.3432676955827567), (443.28888888888889, 145.54920634920634), (442.47169811320754, 99.345911949685529), (448.38235294117646, 61.676470588235297), (460.94339622641508, 100.12578616352201), (462.9594936708861, 144.1113924050633), (466.64705882352939, 60.441176470588232), (476.21568627450978, 175.94117647058823), (481.65853658536588, 143.49864498644988), (479.21367521367523, 99.282051282051285), (507.28163992869878, 96.748663101604279), (500.46020761245677, 144.54325259515571), (523.13114754098365, 145.36065573770492), (529.47872340425533, 126.72340425531915), (548.36363636363637, 83.015151515151516), (554.2388059701492, 81.880597014925371), (563.0333333333333, 103.40000000000001)] +centers[:,0] +np.array(centers) +#[Out]# array([[ 274.992891 , 99.92654028], +#[Out]# [ 7.11290323, 30.67741935], +#[Out]# [ 28.69811321, 103.75471698], +#[Out]# [ 42.08928571, 112.85714286], +#[Out]# [ 60.8630137 , 45.00684932], +#[Out]# [ 154.04898312, 112.49026395], +#[Out]# [ 71.5 , 32.4137931 ], +#[Out]# [ 89.5 , 60.38888889], +#[Out]# [ 184.49112426, 34.2800789 ], +#[Out]# [ 222.828125 , 47.25 ], +#[Out]# [ 238.68965517, 186.10344828], +#[Out]# [ 241.40860215, 52.47311828], +#[Out]# [ 241.86419753, 121.19753086], +#[Out]# [ 250.85910653, 163.26116838], +#[Out]# [ 259.62222222, 56.81481481], +#[Out]# [ 259.84615385, 121.51923077], +#[Out]# [ 269.16412214, 156.87022901], +#[Out]# [ 445.33201754, 185.5127193 ], +#[Out]# [ 278.58992806, 67.21582734], +#[Out]# [ 286.66181818, 153.59272727], +#[Out]# [ 297.10619469, 74.13274336], +#[Out]# [ 305.29955947, 151.1277533 ], +#[Out]# [ 318.63274336, 77.8539823 ], +#[Out]# [ 327.37853107, 144.00564972], +#[Out]# [ 342.48648649, 78.09009009], +#[Out]# [ 350.68783069, 140.42328042], +#[Out]# [ 360.32911392, 78.4556962 ], +#[Out]# [ 368.29381443, 140.70103093], +#[Out]# [ 377.85057471, 74.40229885], +#[Out]# [ 386.11612903, 141.27741935], +#[Out]# [ 413.25593824, 114.46674584], +#[Out]# [ 495.27833954, 6.3432677 ], +#[Out]# [ 443.28888889, 145.54920635], +#[Out]# [ 442.47169811, 99.34591195], +#[Out]# [ 448.38235294, 61.67647059], +#[Out]# [ 460.94339623, 100.12578616], +#[Out]# [ 462.95949367, 144.11139241], +#[Out]# [ 466.64705882, 60.44117647], +#[Out]# [ 476.21568627, 175.94117647], +#[Out]# [ 481.65853659, 143.49864499], +#[Out]# [ 479.21367521, 99.28205128], +#[Out]# [ 507.28163993, 96.7486631 ], +#[Out]# [ 500.46020761, 144.5432526 ], +#[Out]# [ 523.13114754, 145.36065574], +#[Out]# [ 529.4787234 , 126.72340426], +#[Out]# [ 548.36363636, 83.01515152], +#[Out]# [ 554.23880597, 81.88059701], +#[Out]# [ 563.03333333, 103.4 ]]) +centers = np.array(centers) +centers[:,0] +#[Out]# array([ 274.992891 , 7.11290323, 28.69811321, 42.08928571, +#[Out]# 60.8630137 , 154.04898312, 71.5 , 89.5 , +#[Out]# 184.49112426, 222.828125 , 238.68965517, 241.40860215, +#[Out]# 241.86419753, 250.85910653, 259.62222222, 259.84615385, +#[Out]# 269.16412214, 445.33201754, 278.58992806, 286.66181818, +#[Out]# 297.10619469, 305.29955947, 318.63274336, 327.37853107, +#[Out]# 342.48648649, 350.68783069, 360.32911392, 368.29381443, +#[Out]# 377.85057471, 386.11612903, 413.25593824, 495.27833954, +#[Out]# 443.28888889, 442.47169811, 448.38235294, 460.94339623, +#[Out]# 462.95949367, 466.64705882, 476.21568627, 481.65853659, +#[Out]# 479.21367521, 507.28163993, 500.46020761, 523.13114754, +#[Out]# 529.4787234 , 548.36363636, 554.23880597, 563.03333333]) +xs = centers[:,0] +ys = centers[:,0] +xs = centers[:,1] +plt.imshow(m, cmap="Greys_r") +#[Out]# +plt.plot(xs, ys, "r+") +#[Out]# [] +plt.show() +get_ipython().magic(u'logstart') +get_ipython().magic(u'pwd ') +#[Out]# u'/home/tpolgrabia/Dokumenty/prace/signal-processing/shoe-soles' +help %logstart +get_ipython().magic(u'pinfo %logstart') +get_ipython().magic(u'logstart -o') +get_ipython().magic(u'logstop') +get_ipython().magic(u'logstart -o') +get_ipython().magic(u'logstop') +# nie, 10 lip 2016 18:10:54 +get_ipython().magic(u'pwd ') +#[Out]# u'/home/tpolgrabia/Dokumenty/prace/signal-processing/shoe-soles' +# nie, 10 lip 2016 18:10:59 +get_ipython().system(u'cat ipython_log.py') +# nie, 10 lip 2016 18:11:52 +pw +# nie, 10 lip 2016 18:12:19 +get_ipython().magic(u'edit') +#[Out]# 'corners = np.zeros(img.shape, dtype=np.bool)\n\nfor i in range(0, nb_labels):\n mseg = img_label == i\n seg_corner = calc_shape_corners(mseg)\n corners += seg_corner\n\n(n,m) = img.shape\nmat_x = np.tile(np.arange(0, m),n).reshape((n,m))\nmat_y = np.tile(np.arange(0, n),m).reshape((m,n)).T\n\nxs = mat_x[corners]\nys = mat_y[corners]\n\nplt.imshow(img, cmap="Greys_r")\nplt.plot(xs, ys, "b+")\nplt.show()\n' +# nie, 10 lip 2016 18:16:31 +get_ipython().magic(u'edit _') +#[Out]# 'corners = np.zeros(img.shape, dtype=np.bool)\n\nfor i in range(0, nb_labels):\n mseg = img_label == i\n seg_corner = calc_shape_corners(mseg)\n corners += seg_corner\n\n(n,m) = img.shape\nmat_x = np.tile(np.arange(0, m),n).reshape((n,m))\nmat_y = np.tile(np.arange(0, n),m).reshape((m,n)).T\n\nxs = mat_x[corners]\nys = mat_y[corners]\n\nplt.imshow(img, cmap="Greys_r")\nplt.plot(xs, ys, "b+")\nplt.show()\n' +# nie, 10 lip 2016 18:17:03 +get_ipython().magic(u'edit dect.py') +# nie, 10 lip 2016 18:17:28 +get_ipython().magic(u'edit __') +#[Out]# 'corners = np.zeros(img.shape, dtype=np.bool)\n\nfor i in range(0, nb_labels):\n mseg = img_label == i\n seg_corner = calc_shape_corners(mseg)\n corners += seg_corner\n\n(n,m) = img.shape\nmat_x = np.tile(np.arange(0, m),n).reshape((n,m))\nmat_y = np.tile(np.arange(0, n),m).reshape((m,n)).T\n\nxs = mat_x[corners]\nys = mat_y[corners]\n\nplt.imshow(img, cmap="Greys_r")\nplt.plot(xs, ys, "b+")\nplt.show()\n' +# nie, 10 lip 2016 18:18:12 +get_ipython().magic(u'edit _') +#[Out]# 'corners = np.zeros(img.shape, dtype=np.bool)\n\nfor i in range(0, nb_labels):\n mseg = img_label == i\n seg_corner = calc_shape_corners(mseg)\n corners = corners + seg_corner\n\n(n,m) = img.shape\nmat_x = np.tile(np.arange(0, m),n).reshape((n,m))\nmat_y = np.tile(np.arange(0, n),m).reshape((m,n)).T\n\nxs = mat_x[corners]\nys = mat_y[corners]\n\nplt.imshow(img, cmap="Greys_r")\nplt.plot(xs, ys, "b+")\nplt.show()\n' +# nie, 10 lip 2016 18:18:50 +get_ipython().magic(u'edit _') +#[Out]# 'corners = np.zeros(img.shape, dtype=np.bool)\n\nfor i in range(0, nb_labels):\n mseg = img_label == i\n seg_corner = calc_shape_corners(mseg)\n plt.imshow(seg_corner, cmap="Greys_r")\n plt.show()\n corners = corners + seg_corner\n\n(n,m) = img.shape\nmat_x = np.tile(np.arange(0, m),n).reshape((n,m))\nmat_y = np.tile(np.arange(0, n),m).reshape((m,n)).T\n\nxs = mat_x[corners]\nys = mat_y[corners]\n\nplt.imshow(img, cmap="Greys_r")\nplt.plot(xs, ys, "b+")\nplt.show()\n' +# nie, 10 lip 2016 18:20:23 +get_ipython().magic(u'edit __') +#[Out]# 'corners = np.zeros(img.shape, dtype=np.bool)\n\nfor i in range(0, nb_labels):\n mseg = img_label == i\n plt.imshow(mseg, cmap="Greys_r")\n seg_corner = calc_shape_corners(mseg)\n corners = corners + seg_corner\n\n(n,m) = img.shape\nmat_x = np.tile(np.arange(0, m),n).reshape((n,m))\nmat_y = np.tile(np.arange(0, n),m).reshape((m,n)).T\n\nxs = mat_x[corners]\nys = mat_y[corners]\n\nplt.imshow(img, cmap="Greys_r")\nplt.plot(xs, ys, "b+")\nplt.show()\n' +# nie, 10 lip 2016 18:20:44 +get_ipython().magic(u'edit __') +#[Out]# 'corners = np.zeros(img.shape, dtype=np.bool)\n\nfor i in range(0, nb_labels):\n mseg = img_label == i\n seg_corner = calc_shape_corners(mseg)\n plt.imshow(seg_corner, cmap="Greys_r")\n plt.show()\n corners = corners + seg_corner\n\n(n,m) = img.shape\nmat_x = np.tile(np.arange(0, m),n).reshape((n,m))\nmat_y = np.tile(np.arange(0, n),m).reshape((m,n)).T\n\nxs = mat_x[corners]\nys = mat_y[corners]\n\nplt.imshow(img, cmap="Greys_r")\nplt.plot(xs, ys, "b+")\nplt.show()\n' +# nie, 10 lip 2016 18:21:07 +get_ipython().magic(u'edit _') +#[Out]# 'corners = np.zeros(img.shape, dtype=np.bool)\n\nfor i in range(0, nb_labels):\n mseg = img_label == i\n plt.imshow(mseg, cmap="Greys_r")\n plt.show()\n seg_corner = calc_shape_corners(mseg)\n corners = corners + seg_corner\n\n(n,m) = img.shape\nmat_x = np.tile(np.arange(0, m),n).reshape((n,m))\nmat_y = np.tile(np.arange(0, n),m).reshape((m,n)).T\n\nxs = mat_x[corners]\nys = mat_y[corners]\n\nplt.imshow(img, cmap="Greys_r")\nplt.plot(xs, ys, "b+")\nplt.show()\n' +# nie, 10 lip 2016 18:22:01 +get_ipython().magic(u'edit test.py') +# nie, 10 lip 2016 18:22:10 +get_ipython().magic(u'save 1-92 try.py') +# nie, 10 lip 2016 18:22:16 +get_ipython().magic(u'save try.py 1-92') +# nie, 10 lip 2016 18:22:22 +get_ipython().magic(u'edit test.py') +# nie, 10 lip 2016 18:22:33 +get_ipython().magic(u'edit ___') +#[Out]# 'corners = np.zeros(img.shape, dtype=np.bool)\n\nfor i in range(0, nb_labels):\n mseg = img_label == i\n# plt.imshow(mseg, cmap="Greys_r")\n seg_corner = calc_shape_corners(mseg)\n corners = corners + seg_corner\n\n(n,m) = img.shape\nmat_x = np.tile(np.arange(0, m),n).reshape((n,m))\nmat_y = np.tile(np.arange(0, n),m).reshape((m,n)).T\n\nxs = mat_x[corners]\nys = mat_y[corners]\n\nplt.imshow(img, cmap="Greys_r")\nplt.plot(xs, ys, "b+")\nplt.show()\n' +# nie, 10 lip 2016 18:22:58 +get_ipython().magic(u'edit ___') +#[Out]# 'corners = np.zeros(img.shape, dtype=np.bool)\n\nfor i in range(0, nb_labels):\n mseg = img_label == i\n seg_corner = calc_shape_corners(mseg)\n plt.imshow(seg_corner, cmap="Greys_r")\n plt.show()\n corners = corners + seg_corner\n\n(n,m) = img.shape\nmat_x = np.tile(np.arange(0, m),n).reshape((n,m))\nmat_y = np.tile(np.arange(0, n),m).reshape((m,n)).T\n\nxs = mat_x[corners]\nys = mat_y[corners]\n\nplt.imshow(img, cmap="Greys_r")\nplt.plot(xs, ys, "b+")\nplt.show()\n' +# nie, 10 lip 2016 18:23:05 +get_ipython().magic(u'edit __') +#[Out]# 'corners = np.zeros(img.shape, dtype=np.bool)\n\nfor i in range(0, nb_labels):\n mseg = img_label == i\n# plt.imshow(mseg, cmap="Greys_r")\n seg_corner = calc_shape_corners(mseg)\n corners = corners + seg_corner\n\n(n,m) = img.shape\nmat_x = np.tile(np.arange(0, m),n).reshape((n,m))\nmat_y = np.tile(np.arange(0, n),m).reshape((m,n)).T\n\nxs = mat_x[corners]\nys = mat_y[corners]\n\nplt.imshow(img, cmap="Greys_r")\nplt.plot(xs, ys, "b+")\nplt.show()\n' +# nie, 10 lip 2016 18:23:13 +get_ipython().magic(u'edit _') +#[Out]# 'corners = np.zeros(img.shape, dtype=np.bool)\n\nfor i in range(0, nb_labels):\n mseg = img_label == i\n# plt.imshow(mseg, cmap="Greys_r")\n seg_corner = calc_shape_corners(mseg)\n corners = corners + seg_corner\n\n(n,m) = img.shape\nmat_x = np.tile(np.arange(0, m),n).reshape((n,m))\nmat_y = np.tile(np.arange(0, n),m).reshape((m,n)).T\n\nxs = mat_x[corners]\nys = mat_y[corners]\n\nplt.imshow(img, cmap="Greys_r")\nplt.plot(xs, ys, "b+")\nplt.show()\n' +# nie, 10 lip 2016 18:24:06 +get_ipython().magic(u'reload_ext') +# nie, 10 lip 2016 18:24:50 +get_ipython().magic(u'load_ext %autoreload') +# nie, 10 lip 2016 18:25:30 +get_ipython().magic(u'edit _') +#[Out]# 'corners = np.zeros(img.shape, dtype=np.bool)\n\nfor i in range(0, nb_labels):\n mseg = img_label == i\n# plt.imshow(mseg, cmap="Greys_r")\n seg_corner = calc_shape_corners(mseg)\n corners = corners + seg_corner\n\n(n,m) = img.shape\nmat_x = np.tile(np.arange(0, m),n).reshape((n,m))\nmat_y = np.tile(np.arange(0, n),m).reshape((m,n)).T\n\nxs = mat_x[corners]\nys = mat_y[corners]\n\nplt.imshow(img, cmap="Greys_r")\nplt.plot(xs, ys, "b+")\nplt.show()\n' +# nie, 10 lip 2016 18:26:10 +with open("processing.py", "w+") as f: + f.write(Out[102]) + +# nie, 10 lip 2016 18:26:14 +get_ipython().system(u'cat processing.py') diff --git a/2015/2015_03/shoe-soles/old/try.py b/2015/2015_03/shoe-soles/old/try.py new file mode 100644 index 0000000..3ea07b8 --- /dev/null +++ b/2015/2015_03/shoe-soles/old/try.py @@ -0,0 +1,93 @@ +# coding: utf-8 +from dect import calc_shape_corners +calc_shape_corners(img) +from scipy import ndimage +img = ndimage.imread("/home/tpolgrabia/Pobrane/CSFID/tracks_cropped/00003.jpg") +img.shape +import matplotlib.pyplot as plt +plt.imshow(img, cmap="Greys_r") +plt.show() +from skimage import filters +otsu_lvl = filters.threshold_otsu(img) +otsu_lvl +plt.imshow(img, cmap="Greys_r") +plt.show() +plt.imshow(img >= otsu_lvl, cmap="Greys_r") +plt.show() +m = img <= otsu_lvl +plt.imshow(m, cmap="Greys_r") +plt.show() +img_label, nb_labels = ndimage.label(m) +nb_labels +plt.imshow(nb_labels) +plt.imshow(img_label)) +plt.imshow(img_label) +plt.show() +sizes = ndimage.sum(m, img_label, range(0, nb_labels+1)) +sizes +size_min_limit = 50 +remove_segments = sizes < size_min_limit +remove_segments[img_label] +remove_mask = remove_segments[img_label] +img_label[remove_mask] = 0 +plt.imshow(img_label) +plt.show() +import numpy as np +get_ipython().magic(u'pinfo np.unique') +labels = np.unique(img_label) +labels +get_ipython().magic(u'pinfo np.searchsorted') +img_label = np.searchsorted(labels, img_label) +nb_labels = len(labels) +plt.imshow(img_label) +plt.show() +ndimage.center_of_mass(img_label, labels, nb_labels+1) +ndimage.center_of_mass(img_label, labels, nb_labels+1) +get_ipython().magic(u'pinfo ndimage.center_of_mass') +ndimage.center_of_mass(m, img_label, labels) +centers = ndimage.center_of_mass(m, img_label, labels) +centers +centers[0] +centers[1] +centers[2] +centers[3] +centers[4] +centers[5] +labels] +labels +labels = np.unique(img_label) +labels +centers = ndimage.center_of_mass(m, img_label, labels) +centers +centers[:,0] +np.array(centers) +centers = np.array(centers) +centers[:,0] +xs = centers[:,0] +ys = centers[:,0] +xs = centers[:,1] +plt.imshow(m, cmap="Greys_r") +plt.plot(xs, ys, "r+") +plt.show() +get_ipython().magic(u'logstart') +get_ipython().magic(u'pwd ') +help %logstart +get_ipython().magic(u'pinfo %logstart') +get_ipython().magic(u'logstart -o') +get_ipython().magic(u'logstop') +get_ipython().magic(u'logstart -o') +get_ipython().magic(u'logstop') +get_ipython().magic(u'logstart -o -t') +get_ipython().magic(u'pwd ') +get_ipython().system(u'cat ipython_log.py') +pw +get_ipython().magic(u'edit') +get_ipython().magic(u'edit _') +get_ipython().magic(u'edit dect.py') +get_ipython().magic(u'edit __') +get_ipython().magic(u'edit _') +get_ipython().magic(u'edit _') +get_ipython().magic(u'edit __') +get_ipython().magic(u'edit __') +get_ipython().magic(u'edit _') +get_ipython().magic(u'edit test.py') diff --git a/2015/2015_03/shoe-soles/poi_shoe_detection.py b/2015/2015_03/shoe-soles/poi_shoe_detection.py new file mode 100644 index 0000000..564e56f --- /dev/null +++ b/2015/2015_03/shoe-soles/poi_shoe_detection.py @@ -0,0 +1,43 @@ +# coding: utf-8 +import sys +import matplotlib.pyplot as plt +from scipy import ndimage +path = sys.argv[1] +print("I am reading {}".format(path)) +img = ndimage.imread(path) +print("I have read") +img_gray = 0.21 * img[:,:,0] + 0.72 * img[:,:,1] + 0.07 * img[:,:,2] +# plt.imshow(img_gray, cmap="Greys_r") +# plt.show() +from sklearn.cluster import KMeans +k = KMeans(n_clusters=5) +img_reshaped = img_gray.reshape((-1,1)) +k.fit(img_reshaped) +values = k.cluster_centers_.squeeze() +labels = k.labels_ +# labels +img_labels = labels.reshape(img_gray.shape) +# img_labels +# plt.imshow(img_labels) +# plt.show() +# get_ipython().magic(u'pinfo ndimage.median_filter') +# get_ipython().magic(u'pinfo ndimage.median_filter') +# get_ipython().magic(u'pinfo ndimage.median_filter') +img_labels2 = ndimage.median_filter(img_labels, 3) +# plt.imshow(img_labels2) +# plt.show() +from skimage.feature import corner_harris, corner_subpix, corner_peaks +coords = corner_peaks(corner_harris(img_labels2), min_distance=5) +coords_subpix = corner_subpix(img_labels2, coords, window_size=13) +fig, ax = plt.subplots() +ax.imshow(img_labels2, interpolation='nearest', cmap=plt.cm.gray) +ax.plot(coords[:, 1], coords[:, 0], '.b', markersize=3) +ax.plot(coords_subpix[:, 1], coords_subpix[:, 0], '+r', markersize=15) +# ax.axis((0, 350, 350, 0)) +plt.show() +# fig, ax = plt.subplots() +# ax.imshow(img_labels2, interpolation='nearest', cmap=plt.cm.gray) +# ax.plot(coords[:, 1], coords[:, 0], '.b', markersize=3) +# ax.plot(coords_subpix[:, 1], coords_subpix[:, 0], '+r', markersize=15) +# plt.show() +# get_ipython().magic(u'save poi_shoe_detection.py 1-39') diff --git a/2015/2015_03/shoe-soles/poi_shoe_detection2.py b/2015/2015_03/shoe-soles/poi_shoe_detection2.py new file mode 100644 index 0000000..47b0d33 --- /dev/null +++ b/2015/2015_03/shoe-soles/poi_shoe_detection2.py @@ -0,0 +1,242 @@ +import numpy as np +import math +from scipy import ndimage,spatial +import matplotlib.pyplot as plt +from skimage import feature +from skimage.feature import corner_harris, corner_subpix, corner_peaks +from sklearn.cluster import KMeans +from skimage.filters import threshold_otsu +from skimage.draw import line_aa +import math +import sys + +def grayscale_luminosity(img): + return 0.21*img[:,:,0] + 0.72 * img[:,:,1] + 0.07 * img[:,:,2] + +def kmeans_cluster(img_gray): + img_gray_reshaped = img_gray.reshape((-1,1)) + print("About to start k-means clustering") + k = KMeans(n_clusters = 5) + print("Starting k-means clustring") + k.fit(img_gray_reshaped) + print("Finished k-means clustering") + values = k.cluster_centers_.squeeze() + labels = k.labels_ + return (values, labels) + +def corners_seg(img_labels2): + coords = corner_peaks(corner_harris(img_labels2), min_distance=5) + coords_subpix = corner_subpix(img_labels2, coords, window_size=13) + return (coords, coords_subpix) + +def reassign_labels(labels1, values1, labels2, values2): + # TODO not efficient computionally + if values1.length != values2.length: + return None + + n = values1.length + c_min = 1e8 # TODO max value of float + reassignment = np.array((n, 2)) + for i in range(0, n): + e1 = values1[i] + for j in range(0, n): + e2 = values2[j] + d = np.sum(np.abs(e2-e1)) + if d < c_min: + c_min = d + reassignment[i,0] = labels1[i] + reassignment[i,1] = labels2[j] + + return reassignment + +def compare_by_poi(img1, img2, xwindow = 5, ywindow = 5): + # both images are matrix of 0s or 1s + k = np.ones((xwindow,ywindow)) + middlex = (xwindow-1) / 2 + middley = (ywindow-1) / 2 + k[middlex,middley] = 0 + img2_neighbour_count_int = ndimage.convolve(img2) + img2_neighbour_count_bool = img2_neighbour_count_int > 0 + img2_neighbour_count_int2 = np.array(img2_neighbour_count_bool, dtype=np.uint32) + img_comparision = img2_neighbour_count_int2 * img1 + n1 = np.sum(img1) + n2 = np.sum(img_comparision) + return float(n2)/float(n1) + +def find_first_marked_point(crossing_mask): + (n,m) = crossing_mask.shape + for i in range(0,n): + for j in range(0,m): + if crossing_mask[i,j]: + return (j,i) + + return None + +def MSE(a1, a2): + c = 1 + for n in a1.shape: + c *= n + + d = a1 - a2 + return np.sum(d*d) / n + +nsampling = 100 + +def calculate_shape_vector( + path, + nsampling, + median_filter_size = 10, + sweep_line_length = 1024.0): + + img = ndimage.imread(path) + + print("Image read") + img_gray_not_filtered = grayscale_luminosity(img) + img_gray = ndimage.median_filter(img_gray_not_filtered, + size=median_filter_size) + # fig = plt.figure(figsize=(15,15)) + otsu_lvl = threshold_otsu(img_gray) + print(otsu_lvl) + mask = img_gray <= otsu_lvl + + # plt.imshow(mask, cmap="Greys_r") + # plt.show() + + otsu_sub_mask_lvl = threshold_otsu(img_gray[mask]) + mask2 = img_gray <= otsu_sub_mask_lvl + + sub_mask = mask * mask2 + sub_mask2 = mask * (True ^ mask2) + + # plt.imshow(sub_mask, cmap="Greys_r") + # plt.show() + + # plt.imshow(sub_mask2, cmap="Greys_r") + # plt.show() + + label_im, nb_labels = ndimage.label(sub_mask2) + print("Nr of regions: {}".format(nb_labels)) + + sizes = ndimage.sum(sub_mask2, label_im, range(nb_labels + 1)) + + mask_size = sizes < 500 + print("Liczba usunietych: {}".format(np.sum(mask_size))) + remove_pixel = mask_size[label_im] + label_im[remove_pixel] = 0 + # plt.imshow(label_im) + # plt.show() + + centers = np.zeros((nb_labels+1-np.sum(mask_size),2)) + idx = 0 + for i in range(0, nb_labels+1): + if mask_size[i]: + continue + + # processing that means calculating segment parameters + + mask_segment = label_im == i + centers[idx,:] = ndimage.center_of_mass(mask_segment) + idx += 1 + + print("Centers: {}".format(centers)) + centers2 = np.copy(centers) + (n,m) = img_gray.shape + centers2[:,0] = centers[:,0] / float(n) + centers2[:,1] = centers[:,0] / float(m) + + print("Centers 2: {}".format(centers2)) + + mask_int = np.array(mask, dtype=np.uint32) + label_im, nb_labels = ndimage.label(mask_int) + sizes = ndimage.sum(mask, label_im, range(nb_labels + 1)) + mask_size = sizes < 0 + remove_pixel = mask_size[label_im] + label_im[remove_pixel] = 0 + + edges2 = feature.canny(label_im > 0) + + # plt.imshow(edges2, cmap="Greys_r") + + sx = 0.0 + sy = 0.0 + c = 0 + + (n,m) = edges2.shape + print("Width: {}, height: {}".format(n, m)) + for i in range(0, n): + for j in range(0, m): + if edges2[i,j]: + c += 1 + sx += j + sy += i + + sx /= c + sy /= c + + print("Middle ({}x{})".format(sx,sy)) + # plt.plot([sx], [sy], "r+") + + line_length = sweep_line_length + angle_step = 2.0 * np.pi / nsampling + dist_vector = np.zeros((nsampling), dtype=np.float32) + dist_vector[:] = 1e8 # TODO add here infinitium + for i in range(0,nsampling): + line_arr = np.zeros((n,m), dtype=np.bool) + rr,cc,val = line_aa(int(math.floor(sy)), + int(math.floor(sx)), + int(math.floor(sy+line_length*math.sin(angle_step*i))), + int(math.floor(sx+line_length*math.cos(angle_step*i)))) + m1 = rr >= 0 + m2 = rr < n + m3 = cc >= 0 + m4 = cc < m + mall = m1 * m2 * m3 * m4 + cc2 = cc[mall] + rr2 = rr[mall] + val2 = val[mall] + + line_arr[rr2,cc2] = val2 > 0 + crossing_mask = (edges2 * line_arr) + # TODO find one (first - best) crossing point and mark him + p = find_first_marked_point(crossing_mask) + if p != None: + plt.plot([p[0]], [p[1]], "r+") + dist_vector[i] = math.sqrt((p[0] - sx)*(p[0] - sx) + (p[1] - sy)*(p[1] - sy)) + + # plt.show() + return (mask, dist_vector / np.min(dist_vector), centers2) + +if len(sys.argv) < 2: + print("Too small arguments") + sys.exit(1) + +(img_gray1, v1, centers1) = calculate_shape_vector(sys.argv[1],nsampling) +if len(sys.argv) < 3: + print("To small arguments to compare") + sys.exit(0) + +(img_gray2, v2, centers2) = calculate_shape_vector(sys.argv[2],nsampling) +# (img_gray2, v2) = calculate_shape_vector("italian-sole-new-800_1_2_1.jpg",nsampling) + +def compare_sets(c1, c2): + k = spatial.KDTree(c1) + s = 0.0 + (n,m) = c2.shape + for i in range(0, n): + (d,idx) = k.query(c2[i,:]) + s += d*d + + s /= float(n) + return math.sqrt(s) + +print("v1") +print(v1) +print("v2") +print(v2) +diff = v1-v2 +print("Diff: {}".format(diff)) +mse = np.sum(diff*diff)/nsampling +print("MSE: {}".format(MSE(v1,v2))) + +print("Diff between sets of points: {}" + .format(compare_sets(centers1, centers2))) diff --git a/2015/2015_03/shoe-soles/print-processing.py b/2015/2015_03/shoe-soles/print-processing.py new file mode 100644 index 0000000..42a7a80 --- /dev/null +++ b/2015/2015_03/shoe-soles/print-processing.py @@ -0,0 +1,73 @@ +from scipy import ndimage +import math +import matplotlib.pyplot as plt +from skimage import filters +import numpy as np +from skimage import feature + +img = ndimage.imread("/home/tpolgrabia/Pobrane/CSFID/tracks_cropped/00003.jpg") +otsu_lvl = filters.threshold_otsu(img) +img_mask = img >= otsu_lvl +img_mask_fil = True ^ ndimage.minimum_filter(img_mask, 3) +img_label, nb_labels = ndimage.label(img_mask_fil) +sizes = ndimage.sum(img_mask_fil, img_label, range(0, nb_labels+1)) +remove_labels = sizes < 100 +lbls = True & remove_labels +lbls[0] = False + +print("Remove labels: {}".format(ndimage.sum(remove_labels))) + +remove_mask = remove_labels[img_label] +img_label[remove_mask] = 0 + +idx = 0 +n2 = nb_labels + 1 - np.sum(remove_labels) +centers = np.zeros((nb_labels+1,2)) +(n,m) = img.shape +corners = np.zeros((n,m), dtype=np.bool) + +for i in range(0, nb_labels+1): + if remove_labels[i]: + continue + + seg = img_label == i + slice_y, slice_x = ndimage.find_objects(seg)[0] + seg_window = seg[slice_y, slice_x] + edge = feature.canny(seg_window) + plt.imshow(edge, cmap="Greys_r") + plt.show() + (en,em) = edge.shape + distances = np.zeros((en,em)) + + centers[idx,:] = ndimage.center_of_mass(seg) + for j in range(0, en): + for k in range(0, em): + if edge[j,k]: + diff = centers[i,:] - np.array([slice_y.start, slice_x.start]) + dx = k - diff[1] + dy = j - diff[0] + d = math.sqrt(dx*dx+dy*dy) + distances[j,k] = d + # max_distances = ndimage.maximum_filter(distances, 20) + # t_corners = np.abs(max_distances - distances) < 1e-6 + # print(t_corners) + # plt.imshow(t_corners, cmap="Greys_r") + # plt.show() +# corners[slice_y,slice_x] = corners[slice_y,slice_x] + edge + + idx += 1 + +(n,m) = img.shape + +mx = np.tile(np.arange(0,m),n).reshape((n,m)) +my = np.tile(np.arange(0,n),m).reshape((m,n)).T +xs2 = mx[corners] +ys2 = my[corners] + +xs = centers[:,1] +ys = centers[:,0] + +plt.imshow(img_label, cmap="Greys_r") +# plt.plot(xs, ys, "r+") +plt.plot(xs2,ys2,"b+") +plt.show() diff --git a/2015/2015_03/shoe-soles/processing.py b/2015/2015_03/shoe-soles/processing.py new file mode 100644 index 0000000..c8a24ce --- /dev/null +++ b/2015/2015_03/shoe-soles/processing.py @@ -0,0 +1,63 @@ +from dect import calc_shape_corners +import math +from scipy import ndimage +from skimage import filters +import numpy as np +import matplotlib.pyplot as plt + +img = ndimage.imread("/home/tpolgrabia/Pobrane/CSFID/tracks_cropped/00003.jpg") +(n,m) = img.shape +otsu_lvl = filters.threshold_otsu(img) +m = img <= otsu_lvl +img_label, nb_labels = ndimage.label(m) +sizes = ndimage.sum(m, img_label, range(0, nb_labels+1)) +min_seg_size = 100 +to_be_removed_segs = sizes < min_seg_size +remove_mask = to_be_removed_segs[img_label] +img_label[remove_mask] = 0 +labels = np.unique(img_label) +nb_labels = len(labels) +img_label = np.searchsorted(labels, img_label) +labels = np.unique(img_label) +centers = ndimage.center_of_mass(m, img_label, labels) +centers = np.array(centers, dtype=np.float) + +(nc,who_cares) = centers.shape + +def distance_point(m1, m2): + df = m1-m2 + df2 = df*df + sd2 = np.sum(df2) + d = math.sqrt(sd2) + return d + +# TODO check if generates all rates of distance pairs +for i in range(0, nc): + for j in range(0, i): + for k in range(0, j): + m1 = centers[i,:] + m2 = centers[j,:] + m3 = centers[k,:] + d1 = distance_point(m1, m2) + d2 = distance_point(m2, m3) + d = d1 / d2 + print("{}-{}/{}-{}: {}".format(i, j, j, k, d)) + +corners = np.zeros(img.shape, dtype=np.bool) + +for i in range(0, nb_labels): + mseg = img_label == i +# plt.imshow(mseg, cmap="Greys_r") + seg_corner = calc_shape_corners(mseg) + corners = corners + seg_corner + +(n,m) = img.shape +mat_x = np.tile(np.arange(0, m),n).reshape((n,m)) +mat_y = np.tile(np.arange(0, n),m).reshape((m,n)).T + +xs = mat_x[corners] +ys = mat_y[corners] + +plt.imshow(img, cmap="Greys_r") +plt.plot(xs, ys, "b+") +plt.show() diff --git a/2015/2015_03/shoe-soles/sole-detection.py b/2015/2015_03/shoe-soles/sole-detection.py new file mode 100644 index 0000000..832a77e --- /dev/null +++ b/2015/2015_03/shoe-soles/sole-detection.py @@ -0,0 +1,114 @@ +#!/usr/bin/python3 + +import numpy as np +import matplotlib.pyplot as plt +import sys +from scipy import ndimage +from skimage import feature +from skimage import filters +from skimage.feature import corner_harris, corner_subpix, corner_peaks + +def detect_corner(img, windowx = 2, windowy = 2): + (n,m) = img.shape + corners = np.zeros((n,m), dtype=np.bool) + nr_cuts = 0 + for i in range(0, n): + print("Row: {}".format(i)) + for j in range(0, m): + corner = False + for k in range(-windowy, windowy): + nr_cuts = 0 + for l in range(-windowx, windowx): + y = (i+k)%n + x1 = (j+l)%m + x2 = (j+l+1)%m + if img[y,x1] ^ img[y,x2]: + nr_cuts += 1 + if nr_cuts > 1: + corner = True +# print("{}x{} is corner".format(j,i)) + break + corners[i,j] = corner + return corners + +def produce_poi_list(img): + r = [] + (n,m) = img.shape + for i in range(0, n): + row = [] + for j in range(0, m): + if img[i,j]: + row.append((float(j)/m,float(i)/n)) + r.append(row) + + return r + +# not time-effective comparision +# TODO compare only rows in the r window +def compare_poi_list(pl1, pl2, r): + matches = [] + for pr1 in pl1: + for (x1,y1) in pr1: + for pr2 in pl2: + for (x2,y2) in pr2: + if r * r > ((x1 - x2)*(x1 - x2) + (y1-y2)*(y1-y2)): + matches.append(((x1,y1),(x2,y2))) + return matches + +file_path = sys.argv[1] +print("File {} to be analysed".format(file_path)) + +img = ndimage.imread(file_path) +img_gray = 0.21 * img[:,:,0] + 0.72 * img[:,:,1] + 0.07 * img[:,:,2] + +n,m = img_gray.shape + +print("Width {}, Height {}".format(m,n)) + +otsu_lvl = filters.threshold_otsu(img_gray) +seg_shoe = img_gray <= otsu_lvl +img_shoe_selected = seg_shoe * img_gray +data_selected = img_gray[seg_shoe] +shoe_otsu_lvl = filters.threshold_otsu(data_selected) +img_sole_segs1 = img_shoe_selected <= shoe_otsu_lvl +img_sole_segs = ndimage.median_filter(img_sole_segs1, 5) +img_sole_segs_nr = np.array(img_sole_segs, np.int) +k = np.ones((11, 11)) +img_sole_seg_count = ndimage.convolve(img_sole_segs_nr, k) +m1 = img_sole_seg_count >= 80 +m2 = True ^ img_sole_segs +# m2 = img_sole_seg_count <= 1 +m3 = m1 * m2 +# m3i = np.array(m3, np.int) + +corners = detect_corner(img_sole_segs) +poi_list = produce_poi_list(corners) +print(poi_list) + +m3i = np.array(corners, np.int) + +# color_red = np.ones((n,m,3)) * np.array([255, 0, 0]) +color_red = 255 +img_annotated = img_gray * (1 - m3i) + m3i * color_red +# img_annotated = img_gray + + +# plt.imshow(m1 * m2, cmap="Greys_r") +# plt.imshow(img_annotated, cmap="Greys_r") + +# coords = corner_peaks(corner_harris(img_gray), min_distance=5) +# coords_subpix = corner_subpix(img_gray, coords, window_size=13) + +fig = plt.figure() +ax1 = fig.add_subplot(2,1,1) +ax1.imshow(img_sole_segs, cmap="Greys_r") + +ax2 = fig.add_subplot(2,1,2) +ax2.imshow(img_annotated, cmap="Greys_r") +# ax2.plot(coords[:, 1], coords[:, 0], '.b', markersize=3) +# ax2.plot(coords_subpix[:, 1], coords_subpix[:, 0], '+r', markersize=15) + +plt.show() + +plt.imsave("sole_seg.png", img_sole_segs) +