<div dir="ltr"><div class="gmail_default" style="font-family:monospace,monospace">Hey Tom,<br></div><div class="gmail_default" style="font-family:monospace,monospace">glad your query is faster,<br></div><div class="gmail_default" style="font-family:monospace,monospace">and thanks for your answers (helps my curiosity)!<br>Your <a href="https://github.com/catmaid/CATMAID">CATMAID</a> project is very interesting (I come from a computer vision/graphics background).<br></div><div class="gmail_default" style="font-family:monospace,monospace"><br>Index on 2D + index on range is the way to go then!<br><br></div><div class="gmail_default" style="font-family:monospace,monospace">It is faster because you manually indicate that your data can be explained in 2+1 dimension<br></div><div class="gmail_default" style="font-family:monospace,monospace">which in your case fits the data well (from what I understood, medical imaging is often physically sensed in layers along a main axis, like in .).<br></div><div class="gmail_default" style="font-family:monospace,monospace">3D index are based on R-Tree, which are difficult to balance & build in 2D, and even more in 3D.<br><br></div><div class="gmail_default" style="font-family:monospace,monospace"> * Your server conf/hard looks super fine, I had to ask because sometime people forget to tune ^^<br></div><div class="gmail_default" style="font-family:monospace,monospace"> * If you get tired of the cast to numeric, you could create a custom "floatrange" as explained in <a href="https://www.postgresql.org/docs/current/static/rangetypes.html" target="_blank">the doc</a>.<br></div><div class="gmail_default" style="font-family:monospace,monospace"> * the dual would be interesting because when looking for the original edges, you could look for points in the dual, simplifying partitioning/indexing a lot. This was just a remark I don't think you need it.<br></div><div class="gmail_default" style="font-family:monospace,monospace"> * introducing new nodes is a bad idea indeed, I put it for completeness<br></div><div class="gmail_default" style="font-family:monospace,monospace"> * Because of Rtree, your grid approach can have some cells that overlaps (you relax the strict partitioning of space), which would make this easy to write without locking considerations (which are a super pain to deal with.<br></div><div class="gmail_default" style="font-family:monospace,monospace"> * You could partition by numrange(Z). The key idea would be inheritance and partial partitioning on Z <br></div><div class="gmail_default" style="font-family:monospace,monospace">  - you use table inheritance : treenode_edge would have child tables treenode_edge_N. Each children table contains a defined numrange on Z:  table1(z=0 to z=10), table2(z=8 to 20). You can overlap (or not)<br></div><div class="gmail_default" style="font-family:monospace,monospace">  - you don't perform a full partition : <br>   * you put your edges into the best partition if they entirely fit in it (edge(2,6) would go in table1 )<br></div><div class="gmail_default" style="font-family:monospace,monospace">   * the edges that don't fit entirely in one table (for instance edge(5,22)), you keep them in treenode_edge<br></div><div class="gmail_default" style="font-family:monospace,monospace">  - this could actually also help your slice by X and Y <br></div><div class="gmail_default" style="font-family:monospace,monospace">  - inserting would be easy, updating quite easy.<br></div><div class="gmail_default" style="font-family:monospace,monospace">  - You wouldn't have complex lock issues.<br></div><div class="gmail_default" style="font-family:monospace,monospace"><br></div><div class="gmail_default" style="font-family:monospace,monospace">I don't think you need partitioning for the moment, but if you want to increase the number of nodes from 50 millions to several hundreds or thousands of millions, you'll need it (or switch to MonetDB)!<br></div><div class="gmail_default" style="font-family:monospace,monospace"><br>All the best for this project,<br></div><div class="gmail_default" style="font-family:monospace,monospace">Cheers,<br></div><div class="gmail_default" style="font-family:monospace,monospace">Rémi-C<br></div><div class="gmail_default" style="font-family:monospace,monospace"><br> <br></div><div class="gmail_default" style="font-family:monospace,monospace"><br></div></div><div class="gmail_extra"><br><div class="gmail_quote">2017-01-21 1:18 GMT+01:00 Tom Kazimiers <span dir="ltr"><<a href="mailto:tom@voodoo-arts.net" target="_blank">tom@voodoo-arts.net</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi Rémi,<br>
<br>
Thanks a lot for your thorough reply. Especially scenario A.1 has helped me to improve performance by over an order of magnitude for our most common use case. My reply got a bit longer, sorry, I wanted to provide a clearer picture of my setup. For further tests I changed the system, which is an exact mirror of my previous test setup, but this one has SSDs. The generated query plan is the same, but the timing is of course a little bit bitter. But for completeness, this is my original query now:<span class=""><br>
<br>
 SELECT <a href="http://te.id" rel="noreferrer" target="_blank">te.id</a><br>
 FROM treenode_edge te<br></span>
 WHERE te.edge &&& 'LINESTRINGZ( -537284.0 912464.0 131705.0,<br>
                                 1456444.0  -340912.0 131670.0)'<br>
  AND ST_3DDWithin(te.edge, ST_MakePolygon(ST_GeomFromText<wbr>(<br>
   'LINESTRING( -537284.0 -340912.0 131687.5, 1456444.0 -340912.0 131687.5,<br>
                1456444.0  912464.0 131687.5, -537284.0  912464.0 131687.5,<br>
                -537284.0 -340912.0 131687.5)')), 17.5)<br>
  AND te.project_id = 1;<br>
<br>
Like before, I use &&& to get all edges that have a bounding box that<br>
intersects with the query bounding box. ST_3DDWithin is used to lower<br>
number of false positives, i.e. edges with an intersecting BB, but that<br>
don't intersect the query BB. The query returns 6799 rows and runs in<br>
about 600ms:<br>
<br>
 Bitmap Heap Scan on treenode_edge te  (cost=77.04..6555.06 rows=112 width=8) (actual time=595.209..614.640 rows=6799 loops=1)<br>
   Recheck Cond: (edge &&& '01020000800200000000000000886<wbr>520C100000000A0D82B4100000000C<wbr>8130041000000003C3936410000000<wbr>0C0CE14C100000000B0120041'::ge<wbr>ometry)<br>
   Filter: ((edge && '01030000800100000005000000000<wbr>00000AB6520C10000000006CF14C10<wbr>0000000B012004100000000AB6520C<wbr>100000000C3D82B4100000000B0120<wbr>041000000804D39364100000000C3D<wbr>82B4100000000C8130041000000804<wbr>D3936410000000006CF14C10000000<wbr>0C813004100000000AB6520C100000<wbr>00006CF14C100000000B0120041'::<wbr>geometry) AND (project_id = 1) AND ('0103000080010000000500000000<wbr>000000886520C100000000C0CE14C1<wbr>000000003C130041000000003C3936<wbr>4100000000C0CE14C1000000003C13<wbr>0041000000003C39364100000000A0<wbr>D82B41000000003C13004100000000<wbr>886520C100000000A0D82B41000000<wbr>003C13004100000000886520C10000<wbr>0000C0CE14C1000000003C130041':<wbr>:geometry && st_expand(edge, '17.5'::double precision)) AND _st_3ddwithin(edge, '01030000800100000005000000000<wbr>00000886520C100000000C0CE14C10<wbr>00000003C130041000000003C39364<wbr>100000000C0CE14C1000000003C130<wbr>041000000003C39364100000000A0D<wbr>82B41000000003C130041000000008<wbr>86520C100000000A0D82B410000000<wbr>03C13004100000000886520C100000<wbr>000C0CE14C1000000003C130041'::<wbr>geometry, '17.5'::double precision))<br>
   Rows Removed by Filter: 26<br>
   Heap Blocks: exact=4120<br>
   Buffers: shared hit=24703<br>
   ->  Bitmap Index Scan on treenode_edge_gix  (cost=0.00..77.01 rows=1679 width=0) (actual time=594.444..594.444 rows=6825 loops=1)<br>
         Index Cond: (edge &&& '01020000800200000000000000886<wbr>520C100000000A0D82B4100000000C<wbr>8130041000000003C3936410000000<wbr>0C0CE14C100000000B0120041'::ge<wbr>ometry)<br>
         Buffers: shared hit=20583<br>
 Planning time: 0.661 ms<br>
 Execution time: 615.102 ms<br>
<br>
The bitmap index scan on the ND GIST index is still the major problem.<span class=""><br>
<br>
On Thu, Jan 19, 2017 at 03:40:47PM +0100, Rémi Cura wrote:<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
- hardware : are your tables on SSD, does your index fit in RAM<br>
</blockquote>
<br></span>
Our production and test system have the following setup, it is mainly used for the database and a web-server (no other heavy processes, barely any load).<br>
<br>
- 24 x Intel Xeon CPU E5-2630 v2 @ 2.60GHz<br>
- 128GB RAM (all tables and indices fit into it)<br>
- LSI MegaRAID SAS 2208, RAID 1, 1TB total, built from 2 x 1 TB Samsung  SSD 850 PRO 1TB (this is where the Postgres data directory is)<br>
<br>
The Postgres data directory has a size of 40GB and the size of one complete dump of the database of interest is about 5GB at the moment. Of course we made sure we mounted the partition with the noatime option and we use the ext4 filesystem.<span class=""><br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
- postgres tuning : you have customised postgresql.conf to fit your<br>
hardware<br>
</blockquote>
<br></span>
Yes, we modified it a bit. This is mainly what we changed:<br>
<br>
max_connections = 250<br>
shared_buffers = 16GB<br>
work_mem = 256MB<br>
maintenance_work_mem = 2GB<br>
wal_buffers = 16MB<br>
effective_cache_size = 16GB<br>
default_statistics_target = 1000<br>
<br>
These settings worked well so far for us on most queries. But I wonder if less connections and more work_mem might be useful.<br>
<br>
For the sake of completeness, we also tried the parallel worker feature of<br>
Postgres 9.6, but couldn't see any query using it (since no sequential were used in the first place, I believe).<span class=""><br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
- query writing : your query is optimally written (looks fine here)<br>
</blockquote>
<br></span>
Thanks, this improves my confidence that I didn't miss anything obvious on the query level. :-)<span class=""><br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
- data model (detailed after)<br>
</blockquote>
<br></span>
I add some more details below, but first I want to add some more<br>
context: typically we deal with individual nodes, stored in a separate<br>
treenode table (the result of the PostGIS query is then joined into it).<br>
Each node has a parent (except the root node) and the edge between them<br>
is an actual part of the neuron (matches underlying image data). It will always be a straight line, though. This is also why we need to respect edges during queries: we need to know if the edge between two nodes intersect with a query bounding box, even if the nodes themselves aren't. Edges often cross multiple query bounding boxes/sections, i.e.  nodes are multiple query bounding boxes/sections away. Users create new nodes and modify or delete old ones all the time.<br>
<br>
We query typically in sections of Z, but support also queries in sections of X or Y. Ideally, all are equally fast. Z-sections are, however, the clear focus.<span class=""><br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
IF you want fast results, you need a good indexing, and this is possible if you exploit the sparsness of your data. So where is the sparsness?<br>
</blockquote>
<br></span>
The neurons/trees can cover quite a bit of space. There are many really big ones that cross hundreds of slices. In my current test data set we have about 10.000 neurons, with an average of 6000 nodes each. About 50.000 nodes get added per day. Most neurons span between 1000 and 2500 sections. An typical section we want to improve performance for might involve 5000-10000 edge intersections which might belong to 1000-2000 different neurons/trees, but this also varies a lot by region.<span class=""><br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
* Scenario A.1: your usual slices affects a majority of trees<br>
 - you don't really need an index on X,Y, but mainly on Z.<br>
   The idea is to separate your index into index(X,Y) + index(Z)<br>
    Try adding an index on GIST(numrange(ST_Zmin,ST_ZMax)<wbr>) (you could add<br>
it anyway)<br>
    and an index on GIST(geom) (aka only 2D).<br>
</blockquote>
<br></span>
As long as I can also query from orthogonal perspectives with reasonable performance, I don't need separate X and Y indices, it seems.<br>
<br>
First, I adding the following indices:<br>
<br>
 CREATE INDEX treenode_edge_z_range_gist ON treenode_edge<br>
 USING GIST(numrange(ST_ZMin(edge)::n<wbr>umeric, ST_ZMax(edge)::numeric), '[]');<br>
<br>
 CREATE INDEX treenode_edge_2d_gist ON treenode_edge USING GIST(edge);<br>
<br>
Without the type casts, Postgres wouldn't allow it and inclusive bounds<br>
are the current behavior. My regular query of course doesn't make use of<br>
this index in its current form, but replacing the &&& operator test with<br>
a 2D bound box check plus a Z index condition turns out to be really<br>
fast:<br>
<br>
 EXPLAIN (ANALYZE, BUFFERS) SELECT <a href="http://te.id" rel="noreferrer" target="_blank">te.id</a><br>
 FROM treenode_edge te<br>
 WHERE te.edge && ST_MakeEnvelope(-537284.0, -340912.0, 1456444.0, 912464.0)<br>
   AND numrange(ST_ZMin(te.edge)::num<wbr>eric, ST_ZMax(te.edge)::numeric) && '[131670.0,131705.0]'::numrang<wbr>e<br>
   AND ST_3DDWithin(te.edge, ST_MakePolygon(ST_GeomFromText<wbr>(<br>
     'LINESTRING( -537284.0 -340912.0 131687.5, 1456444.0 -340912.0 131687.5,<br>
                  1456444.0  912464.0 131687.5, -537284.0  912464.0 131687.5,<br>
                  -537284.0 -340912.0 131687.5)')), 17.5)<br>
   AND te.project_id = 1;<br>
 ---<br>
 Bitmap Heap Scan on treenode_edge te  (cost=3398.97..123260.14 rows=5007 width=8) (actual time=12.919..33.256 rows=6799 loops=1)<br>
   Recheck Cond: (numrange((st_zmin((edge)::box<wbr>3d))::numeric, (st_zmax((edge)::box3d))::nume<wbr>ric, '[]'::text) && '[131670.0,131705.0]'::numrang<wbr>e)<br>
   Filter: ((edge && '01030000000100000005000000000<wbr>00000886520C100000000C0CE14C10<wbr>0000000886520C100000000A0D82B4<wbr>1000000003C39364100000000A0D82<wbr>B41000000003C39364100000000C0C<wbr>E14C100000000886520C100000000C<wbr>0CE14C1'::geometry) AND (edge && '01030000800100000005000000000<wbr>00000AB6520C10000000006CF14C10<wbr>0000000B012004100000000AB6520C<wbr>100000000C3D82B4100000000B0120<wbr>041000000804D39364100000000C3D<wbr>82B4100000000C8130041000000804<wbr>D3936410000000006CF14C10000000<wbr>0C813004100000000AB6520C100000<wbr>00006CF14C100000000B0120041'::<wbr>geometry) AND (project_id = 1) AND ('0103000080010000000500000000<wbr>000000886520C100000000C0CE14C1<wbr>000000003C130041000000003C3936<wbr>4100000000C0CE14C1000000003C13<wbr>0041000000003C39364100000000A0<wbr>D82B41000000003C13004100000000<wbr>886520C100000000A0D82B41000000<wbr>003C13004100000000886520C10000<wbr>0000C0CE14C1000000003C130041':<wbr>:geometry && st_expand(edge, '17.5'::double precision)) AND _st_3ddwithin(edge, '01030000800100000005000000000<wbr>00000886520C100000000C0CE14C10<wbr>00000003C130041000000003C39364<wbr>100000000C0CE14C1000000003C130<wbr>041000000003C39364100000000A0D<wbr>82B41000000003C130041000000008<wbr>86520C100000000A0D82B410000000<wbr>03C13004100000000886520C100000<wbr>000C0CE14C1000000003C130041'::<wbr>geometry, '17.5'::double precision))<br>
   Rows Removed by Filter: 26<br>
   Heap Blocks: exact=4120<br>
   Buffers: shared hit=4428<br>
   ->  Bitmap Index Scan on treenode_edge_range_gist  (cost=0.00..3397.72 rows=75374 width=0) (actual time=12.351..12.351 rows=6825 loops=1)<br>
         Index Cond: (numrange((st_zmin((edge)::box<wbr>3d))::numeric, (st_zmax((edge)::box3d))::nume<wbr>ric, '[]'::text) && '[131670.0,131705.0]'::numrang<wbr>e)<br>
         Buffers: shared hit=308<br>
 Planning time: 0.856 ms<br>
 Execution time: 33.731 ms<br>
<br>
This is a very impressive speed-up, thanks for this suggestion! For all<br>
examples I tried the results were exactly the same as with the &&&<br>
operator. I only stumped across the fact the current semantic is to<br>
have inclusive bounds, but numrange()'s behavior is '[)' by default.  Being able to change this is also a nice feature.<br>
<br>
I also tested creating range indices for X and Y so that I could replace<br>
the &&& operator with three separate range checks. This, however, was<br>
much slower with >2500ms. This setup would have been nicer, because it could work equally well for "slice" queries in each dimension.<br>
<br>
Anyhow, using your suggested indices and query works very well and fast even for differently oriented slices. It is really suprising that the two 2D index is faster than a 3D index for us.<span class=""><br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
* Scenario A.2 : your usual slices affect few trees<br>
 - you could easily perform a 2 part query, where the first part use an<br>
index on tree bounding box, and the second perform the query on edge. This<br>
require to create and maintain a tree_bounding_box table with appropriate<br>
indexes, which is quite easy with triggers Ex (not optimised):<br>
  WITH filtering_by_tree AS (<br>
    SELECT distinct tree_id<br>
    FROM my_sync_table_of_tree_bounding<wbr>_boxes as tree<br>
    WHERE tree.bbox &&& your_slice<br>
  ) SELECT edge_id<br>
    FROM treenode_edge  AS te<br>
    WHERE te.geom &&& your_slice<br>
     AND EXISTS (SELECT 1 FROM filtering_by_tree As ft WHERE ft.tree_id =<br>
te.tree_id)<br>
 - you could use pgpointcloud to similar effect by construction a patch<br>
per tree.<br>
</blockquote>
<br></span>
This is certainly also an appealing idea. I quickly tested this by<br>
creating a neuron summary table with its current bounding boxes as box3d<br>
and then queried in a similar way like you suggested above. I couldn't<br>
get this very fast though, regardless of whether &&& or a 2D-approach is<br>
in use. I typically ended up with about 500ms.<span class=""><br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Now something is very ambiguous in your data, and it could have major<br>
consequences :<br>
* Are the edges between nodes that represent neuron __only logical__ (they<br>
are straight lines, always the same Zmax-Zmin, and represent a conceptual<br>
graph), or are they polylines that represent the __actual__ neurone path?<br>
</blockquote>
<br></span>
Edges reflect the morphology of the actual neuron and their length<br>
(even only in Z) can vary a lot.<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class="">
* Scenario B.1.1 : edges are only logical and have a constant Zmax-Zmin,<br>
i.e. nodes have quantified Z :<br>
 - you don't need edges at all for indexing, you should work only on nodes<br>
(i.e. index nodes IN Z + (X,Y), find nodes, then find relevant edges). The<br>
great advantage in this case is that you can use the full power of<br>
pgpointcloud to create meaningfull patches (of variables size for instance,<br>
see this<br></span>
<<a href="http://www.sciencedirect.com/science/article/pii/S092427161630123X" rel="noreferrer" target="_blank">http://www.sciencedirect.com/<wbr>science/article/pii/S092427161<wbr>630123X</a>>, 3.1,<span class=""><br>
Adapt patch grouping rules ), which means scalling in the range of billions<br>
of points possible.<br>
</span></blockquote>
<br>
Thanks, that's good to know, but our Zmax-Zmin can vary.<span class=""><br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
 - you could try to think of using the dual of your current graph (edges<br>
become nodes, nodes become edges).<br>
</blockquote>
<br></span>
Why could this be an advantage here? I guess if there are helpful data structures for a dual representation, search time can be improved.<span class=""><br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
* Scenario B.1.2 : edges are only logical, but nodes don't have a constant<br>
Zmax-Zmin. : you could insert (virtual) nodes so you are again in the same<br>
case as scenario B.1.1<br>
</blockquote>
<br></span>
I would prefer not doing this. We did this in the past, but wanted to<br>
reduce the number of nodes needed and improve usability with large distances.<span class=""><br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
* Scenario B.2 : edges are actual neurone path.<br>
 - you could construct intermediary objects (groups of neurone) of<br>
variable size (those within cubes), so as to have again a 2 steps query :<br>
first look for groups of neuron, then for neurons inside<br>
</blockquote>
<br></span>
This sounds similar to the grid I have implemented at the moment. So far<br>
I was not able to test the update performance of this. In the worst case this requires asynchronous updates that need to be clever about locking.<br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class="">
* Scenario C : your data don't have obvious sparsity.<br>
 - use the same strategy you currently use, but scale by partitionning<br>
your data into numerous tables.<br>
  You can see the partitionning entry of the postgres manual.<br></span>
<<a href="https://www.postgresql.org/docs/current/static/ddl-partitioning.html" rel="noreferrer" target="_blank">https://www.postgresql.org/do<wbr>cs/current/static/ddl-partitio<wbr>ning.html</a>><span class=""><br>
  This would be much easier to create and maintain if you could partition<br>
on only Z dimension.<br>
  Note that you will need to duplicate (and then de-duplicate) some edges.<br>
</span></blockquote>
<br>
Thanks for the pointer. Partitioning sounds interesting, but similarly to the current grid based approach, getting the CHECK constraints right seems not so easy for different use cases. With my grid I could see performance changes when I changed cell size. But I might try and set this up anyway.<span class=""><br>
<br>
<blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
It's hard to recommend anything without further information on your<br>
data/your requirement/your ressources.<br>
</blockquote>
<br></span>
I know and I hope I could provide enough information to make some things clearer. This E-mail got a bit lengthy, but I wanted to address all the points you listed. You have certainly helped me a lot already and I am<br>
currently testing the promising looking 2D based approach (2D GIST and range GIST) further.<br>
<br>
Thanks again, also for the nice explanatory graphics!<span class="HOEnZb"><font color="#888888"><br>
Tom<br>
</font></span><br>______________________________<wbr>_________________<br>
postgis-users mailing list<br>
<a href="mailto:postgis-users@lists.osgeo.org">postgis-users@lists.osgeo.org</a><br>
<a href="https://lists.osgeo.org/mailman/listinfo/postgis-users" rel="noreferrer" target="_blank">https://lists.osgeo.org/<wbr>mailman/listinfo/postgis-users</a><br></blockquote></div><br></div>