using PyPlot using NtToolBox # using Autoreload # arequire("NtToolBox") gamma0 = [.78, .14, .42, .18, .32, .16, .75, .83, .57, .68, .46, .40, .72, .79, .91, .90] + 1im*[.87, .82, .75, .63, .34, .17, .08, .46, .50, .25, .27, .57, .73, .57, .75, .79]; periodize = gamma -> [gamma; [gamma[1]]] function cplot(gamma, s = "b", lw = 1) plot(real(periodize(gamma)), imag(periodize(gamma)), s, linewidth = lw) axis("equal") axis("off") end cplot(gamma0, "b.-") p = 256; using Dierckx curvabs = gamma -> [0; cumsum(1e-5 .+ abs(gamma[1:end - 1] - gamma[2:end]) )] resample11 = (gamma, d) -> evaluate(Spline1D(d./d[end], real(gamma), k=1), (0:p-1)./p) resample12 = (gamma, d) -> evaluate(Spline1D(d./d[end], imag(gamma), k=1), (0:p-1)./p) resample = gamma -> resample11( [gamma; gamma[1]], curvabs( [gamma; gamma[1]] ) ) + 1im*resample12( [gamma; gamma[1]], curvabs( [gamma; gamma[1]] ) ); gamma1 = resample(gamma0); cplot(gamma1, "k"); BwdDiff = c -> c - [c[end]; c[1:end - 1]] FwdDiff = c -> [c[2:end]; c[1]] - c dotp = (c1, c2) -> real(c1.*conj(c2)); normalize = v -> v./max(abs(v), 1e-10) tangent = gamma -> normalize( FwdDiff(gamma) ) normal = gamma -> -1im*tangent(gamma); delta = .03 gamma2 = gamma1 .+ delta .* normal(gamma1) gamma3 = gamma1 .- delta .* normal(gamma1); cplot(gamma1, "k") cplot(gamma2, "r--") cplot(gamma3, "b--") axis("tight") axis("off") normalC = gamma -> BwdDiff(tangent(gamma)) ./ abs( FwdDiff(gamma) ); dt = 0.001 / 100; Tmax = 3.0 / 100 niter = round(Tmax/dt); gamma = gamma1; gamma = gamma .+ dt .* normalC(gamma); gamma = resample(gamma); include("NtSolutions/segmentation_2_snakes_param/exo1.jl") n = 200; nbumps = 40 theta = rand(nbumps,1).*(2*pi) r = .6*n/2 a = [.62*n,.6*n] x = round( a[1] + r*cos(theta) ) y = round( a[2] + r*sin(theta) ) W = zeros(n, n) for i in 1:nbumps W[Int(x[i]), Int(y[i])] = 1 end W = gaussian_blur(W, 6.0) W = rescale( -min(W, .05), .3,1); imageplot(W) G = grad(W) G = G[:,:,1] .+ 1im*G[:, :, 2]; imageplot(abs(G)) EvalG = gamma -> NtToolBox.bilinear_interpolate(G, imag(gamma), real(gamma)) EvalW = gamma -> NtToolBox.bilinear_interpolate(W, imag(gamma), real(gamma)); r = .98*n/2 # radius p = 128 # number of points on the curve theta = transpose( linspace(0, 2*pi, p + 1) ) theta = theta[1 : end - 1] gamma0 = n/2 * (1 + 1im) + r*(cos(theta) + 1im*sin(theta)); gamma = gamma0; dt = 1; Tmax = 5000 niter = round(Tmax/ dt); lw = 2 clf imageplot(transpose(W)) cplot(gamma, "r", lw) dotp = (c1,c2) -> real(c1.*conj(c2)); N = normal(gamma) g = - EvalW(gamma) .* normalC(gamma) + dotp(EvalG(gamma), N) .* N gamma = (gamma - dt*g); gamma = resample(gamma); include("NtSolutions/segmentation_2_snakes_param/exo2.jl") # gamma = copy(gamma0); # displist = round(linspace(0, niter, 10)) # k = 1; # clf; # imageplot(transpose(W)); # track = gamma # for i in 0 : niter # N = normal(gamma); # g = EvalW(gamma) .* normalC(gamma) - dotp(EvalG(gamma), N) .* N; # gamma = resample( (gamma + dt*g)); # if i == displist[k] # lw = 1 # if (i == 0) | (i == niter) # lw = 4 # end # cplot(gamma, "r", lw); # k = k + 1; # axis("equal"); axis("off"); # end # end n = 256 name = "NtToolBox/src/data/cortex.png" f = load_image(name, n); imageplot(f) G = grad(f) d0 = sqrt(sum(G.^2, 3)) imageplot(d0[:, :, 1]) a = 2 d = gaussian_blur(d0, a) imageplot(d[:, :, 1]) f = copy(d0) n = maximum(size(f)) # t = [collect(0:n/2); collect(-n/2:-2)] # (Y, X) = meshgrid(t, t) # h = exp( -(X.^2 + Y.^2)./(2.0*float(sigma)^2) ) # h = h./sum(h) # real(plan_ifft((plan_fft(f)*f).*(plan_fft(h)*h))*((plan_fft(f)*f).*(plan_fft(h)*h))) d = min(d, .4) W = rescale(-d, .8, 1); imageplot(W[:, :, 1]) p = 128; include("NtSolutions/segmentation_2_snakes_param/exo3.jl") dt = 2; Tmax = 9000 niter = round(Tmax/ dt); include("NtSolutions/segmentation_2_snakes_param/exo4.jl") # G = grad(W); # G = G[:, :, 1] + 1im*G[:, :, 2] # EvalG = gamma -> NtToolBox.bilinear_interpolate(G, imag(gamma), real(gamma)) # EvalW = gamma -> NtToolBox.bilinear_interpolate(W, imag(gamma), real(gamma)) # # # gamma = gamma0 # displist = round(linspace(0, niter, 10)) # k = 1; # clf # imageplot(transpose(f)) # for i in 0 : niter - 1 # n = normal(gamma) # g = EvalW(gamma) .* normalC(gamma) - dotp(EvalG(gamma), n) .* n # gamma = resample( gamma + dt*g ) # if i == displist[k] # lw = 1; # if (i == 0) | (i == niter) # lw = 4; # end # cplot(gamma, "r", lw); # k = k + 1; # axis("equal"); axis("off"); # end # end n = 256 f = load_image(name, n) f = f[45:105, 60:120] n = size(f)[1]; imageplot(f) include("NtSolutions/segmentation_2_snakes_param/exo5.jl") x0 = 4 + 55*1im x1 = 53 + 4*1im; p = 128 t = transpose(linspace(0, 1, p)) gamma0 = t*x1 + (1 - t)*x0; gamma = gamma0; clf imageplot(transpose(W)) cplot(gamma[:],"r", 2) plot(real(gamma[1]), imag(gamma[1]), "b.", markersize = 20) plot(real(gamma[end]), imag(gamma[end]), "b.", markersize = 20); using Dierckx curvabs = gamma -> [0; cumsum(1e-5 .+ abs(gamma[1 : end - 1] - gamma[2 : end]) )] resample11 = (gamma, d) -> evaluate(Spline1D(d./d[end], real(gamma), k=1), (0:p-1)./p) resample12 = (gamma, d) -> evaluate(Spline1D(d./d[end], imag(gamma), k=1), (0:p-1)./p) resample = gamma -> resample11( [gamma; gamma[1]], curvabs( [gamma; gamma[1]] ) ) + 1im*resample12( [gamma; gamma[1]], curvabs( [gamma; gamma[1]] ) ) dt = 1/10; Tmax = 2000*4/ 7 niter = round(Tmax/ dt); include("NtSolutions/segmentation_2_snakes_param/exo6.jl")