JuliaEarth / ImageQuilting.jl

Fast image quilting simulation solver for the GeoStats.jl framework
https://github.com/JuliaEarth/GeoStats.jl
MIT License
41 stars 12 forks source link

hard data issue #14

Closed mariosgeo closed 6 years ago

mariosgeo commented 6 years ago

Hi,

I am using the ImageQuilting toolbox to generate realizations of subsurface. On the attached image, the left figure if the training image, the right figure is the borehole location and the center is one of the realizations.

The windows used in iqsim is 40x40. Apparently, this is much bigger than the actual distance between two boreholes. How can I use more efficiently the simulation ? Reducing the window or am I doing something wrong?

figure_1-1

I am attaching the two files I use as taring image and borehole information. The code looks like this

TI1=npzread("flummy.npy") TI1=TI1[:,:,4] TI1= reshape(TI1,201,201,1) geo=npzread("real_geo2.npy") shape = HardData() for i=1:size(geo,1) push!(shape, (geo[i,1]+1,geo[i,2]+1,1)=>geo[i,3]) end reals= iqsim(TI1, 40, 40, 1, size(TI1)..., hard=shape, nreal=1)

Data_for_Marios.zip

juliohm commented 6 years ago

Hi @mariosgeo, thank you for reaching out.

I wouldn't recommend image quilting for dense hard data configurations. We investigated this issue in our paper, which I am sharing here. It can work if you reduce the template size, but it becomes more computationally intensive. For dense hard data configurations, pixel-based algorithms like direct sampling are recommended. I didn't have the time to implement direct sampling yet, but meanwhile I think there is a Python implementation on the internet that one of my colleagues used.

If you are actually simulating a 3D domain or a very large 2D domain, then direct sampling and other existing pixel-based algorithms become too computationally intensive. I would try to reduce the template size and activate the GPU with image quilting to see if it helps: http://juliohm.github.io/ImageQuilting.jl/stable/gpu-support.html

Another parameter of image quilting that is set by default for the user, and that also impacts performance, is the overlap configuration. By default the overlaps are set to 1/6 of the template size, in your case ceil(1/6 * 40) --> 7. You may want to reduce this overlap region by decreasing the fraction further.

Please let me know if that clarifies the issue.

Best,

mariosgeo commented 6 years ago

Thanks for the reply and advice. I tried various square template sizes (5,8,10,20,40,50, attached the one with template=10), with similar results.

figure_1-2

I tried also several overlap settings with the following loop (the title indicates the factor)

`for i=2:20 reals= iqsim(TI1, 10, 10, 1, size(TI1)..., hard=shape,nreal=1,gpu=true,overlapx=1/i,overlapy=1/i,overlapz=1/i)

end`

figure_1-3

Is there something more I can try?

juliohm commented 6 years ago

I did some quick exploration with your data to see if it makes sense at all to try match it with the given training image. Because the data is dense, you can get a rough sense of the spatial patterns with Kriging. First notice that the variogram is very flat, there is almost no correltation at zero lag:

using GeoStats
using NPZ, Plots

csv = npzread("Data_for_Marios/real_geo2.npy")

hd = PointSetData(Dict(:var => csv[:,3]), csv[:,1:2]')

γ = EmpiricalVariogram(hd, :var, maxlag=100.)
γhat = fit(SphericalVariogram, γ)

plot(γ)
plot!(γhat, maxlag=100.)

screenshot from 2018-08-22 12-14-10

I understand this is a categorical variable, but in any case we can ignore the mathematical formalities for a moment and produce a Kriging map, pretending it is continuous. Global Kriging will take some time to run because there are 1349 data points (need to solve a 1349x1349 linear system for every pixel).

grid = RegularGrid{Float64}(200,200)

problem = EstimationProblem(hd, grid, :var)
solution = solve(problem, Kriging(:var => @NT(variogram=γhat)))

plot(solution, size=(1000,400), color=:viridis)

screenshot from 2018-08-22 12-17-07

Maybe you could ask yourself if by breaking the training image into tiles (like image quilting does) you could ever expect to produce patterns like the above. I still think that pixel-based algorithms are more appropriate, but if you decide to continue with image quilting, I suggest trying different training images.

juliohm commented 6 years ago

I am closing the issue. If you have questions again in the future, please ask in our gitter channel. The issue tracker on GitHub is mainly for bug reports and feature requests.

Thanks,

mariosgeo commented 6 years ago

Thanks for the help and reply