tried choosing only relevant voxels but unsuccessful in implementation

main.jl

@@ -530,8 +530,21 @@ function construct_cost(cdata::contourdata, vd::voxdata, mat::material)
 
     # considered voxels
     not_empty_voxels = length.([m.seglen for m in vd.maps]) .>0
-    considered_voxels = (1:length(vd.below))[(vd.below.!=0) .& not_empty_voxels]
-    considered_voxels = considered_voxels[not_empty_voxels[vd.below[considered_voxels]]]
+    valid_voxels = (1:length(vd.below))[(vd.below.!=0) .& not_empty_voxels]
+    valid_voxels = valid_voxels[not_empty_voxels[vd.below[valid_voxels]]]
+
+    #stress_size = .√( vd.voxels.Sx.^2 + vd.voxels.Sy.^2 + vd.voxels.Sz.^2)
+    #top_voxes = Set(sortperm(stress_size)[1:Int(length(stress_size)/4)]) # top stressors only
+    #top_valid_voxels = sort(collect(Set(valid_voxels) ∩ top_voxes))
+    #rel_voxels = vd.voxels[top_valid_voxels,:]
+
+    #relmaps = sort(collect(Set(vd.below[top_valid_voxels]) ∪ Set(top_valid_voxels)))
+    #considered_voxels = indexin(top_valid_voxels, relmaps)
+    # relbelows = indexin(vd.below[top_valid_voxels], relmaps)
+    considered_voxels = valid_voxels
+    relbelows = vd.below[valid_voxels]
+    rel_voxels = vd.voxels[considered_voxels,:]
+    relmaps = vd.maps
 
     # calculate stresses at each voxel
     F = (2/mat.σ̄[1,1]^2 - 1/mat.σ̄[3,3]^2)/2
@@ -539,18 +552,13 @@ function construct_cost(cdata::contourdata, vd::voxdata, mat::material)
     L = 1/(2*mat.σ̄[1,2]^2)
     M = 1/(2*mat.σ̄[1,3]^2)
 
-    # voxtimes = [:(sum($(v.seglen).*($(v.segoffset) + timestart[$(v.segcontours)]))/$(sum(v.seglen))) for v in vd.maps]
-    # voxtimes = voxtimes[considered_voxels]
-    relmaps = vd.maps[considered_voxels]
-    relbelows = vd.below[considered_voxels]
-    rel_voxels = vd.voxels[considered_voxels,:]
     return quote
         function cost_f(rl::Vector{Int})
             # println("Cumulative Sum") # 100k
             timestart = cumsum([$contour_times[c] for c in rl])
 
             #println("Vox times") # 600
-            voxtimes = [v.seglen ⋅ timestart[v.segcontours] + v.c for v in vd.maps]
+            voxtimes = [v.seglen ⋅ timestart[v.segcontours] + v.c for v in relmaps]
 
             # println("time diff") # 28k
             Δt = voxtimes[considered_voxels] - voxtimes[relbelows]
@@ -572,7 +580,6 @@ function construct_cost(cdata::contourdata, vd::voxdata, mat::material)
                 $G * ((σ33 - σ11).^2 + (σ33 - σ22).^2) +
                 $(2 * L) * (σ12).^2 +
                 $(2 * M) * (σ23 + σ31).^2)
-
             return cost
         end
     end
@@ -645,4 +652,11 @@ o=0
         swap!(rl, i, j)
     end
 end
-println(o)
\ No newline at end of file
+println(o)
+
+
+### Visualise distribution of stresses
+data = .√( vd.voxels.Sx.^2 + vd.voxels.Sy.^2 + vd.voxels.Sz.^2)
+data = .√(vd.voxels.Txy.^2 + vd.voxels.Tyz.^2 + vd.voxels.Txz.^2)
+histogram(data) # vast majority of voxels near 0 stress - can ignore
+histogram!(sort(data, rev=true)[1:4000]) # vast majority of voxels near 0 stress - can ignore
\ No newline at end of file