[postgis-users] Improve 3D BB query on large set of lines

Tom Kazimiers tom at voodoo-arts.net
Fri Jan 20 16:18:24 PST 2017

Hi Rémi,

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 

  SELECT te.id
  FROM treenode_edge te
  WHERE te.edge &&& 'LINESTRINGZ( -537284.0 912464.0 131705.0,
                                  1456444.0  -340912.0 131670.0)'
   AND ST_3DDWithin(te.edge, ST_MakePolygon(ST_GeomFromText(
    'LINESTRING( -537284.0 -340912.0 131687.5, 1456444.0 -340912.0 131687.5,
                 1456444.0  912464.0 131687.5, -537284.0  912464.0 131687.5,
                 -537284.0 -340912.0 131687.5)')), 17.5)
   AND te.project_id = 1;

Like before, I use &&& to get all edges that have a bounding box that
intersects with the query bounding box. ST_3DDWithin is used to lower
number of false positives, i.e. edges with an intersecting BB, but that
don't intersect the query BB. The query returns 6799 rows and runs in
about 600ms:

  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)
    Recheck Cond: (edge &&& '01020000800200000000000000886520C100000000A0D82B4100000000C8130041000000003C39364100000000C0CE14C100000000B0120041'::geometry)
    Filter: ((edge && '0103000080010000000500000000000000AB6520C10000000006CF14C100000000B012004100000000AB6520C100000000C3D82B4100000000B0120041000000804D39364100000000C3D82B4100000000C8130041000000804D3936410000000006CF14C100000000C813004100000000AB6520C10000000006CF14C100000000B0120041'::geometry) AND (project_id = 1) AND ('0103000080010000000500000000000000886520C100000000C0CE14C1000000003C130041000000003C39364100000000C0CE14C1000000003C130041000000003C39364100000000A0D82B41000000003C13004100000000886520C100000000A0D82B41000000003C13004100000000886520C100000000C0CE14C1000000003C130041'::geometry && st_expand(edge, '17.5'::double precision)) AND _st_3ddwithin(edge, '0103000080010000000500000000000000886520C100000000C0CE14C1000000003C130041000000003C39364100000000C0CE14C1000000003C130041000000003C39364100000000A0D82B41000000003C13004100000000886520C100000000A0D82B41000000003C13004100000000886520C100000000C0CE14C1000000003C130041'::geometry, '17.5'::double precision))
    Rows Removed by Filter: 26
    Heap Blocks: exact=4120
    Buffers: shared hit=24703
    ->  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)
          Index Cond: (edge &&& '01020000800200000000000000886520C100000000A0D82B4100000000C8130041000000003C39364100000000C0CE14C100000000B0120041'::geometry)
          Buffers: shared hit=20583
  Planning time: 0.661 ms
  Execution time: 615.102 ms

The bitmap index scan on the ND GIST index is still the major problem.

On Thu, Jan 19, 2017 at 03:40:47PM +0100, Rémi Cura wrote:
> - hardware : are your tables on SSD, does your index fit in RAM

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).

- 24 x Intel Xeon CPU E5-2630 v2 @ 2.60GHz
- 128GB RAM (all tables and indices fit into it)
- 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)

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.

> - postgres tuning : you have customised postgresql.conf to fit your

Yes, we modified it a bit. This is mainly what we changed:

max_connections = 250
shared_buffers = 16GB
work_mem = 256MB
maintenance_work_mem = 2GB
wal_buffers = 16MB
effective_cache_size = 16GB
default_statistics_target = 1000

These settings worked well so far for us on most queries. But I wonder 
if less connections and more work_mem might be useful.

For the sake of completeness, we also tried the parallel worker feature of
Postgres 9.6, but couldn't see any query using it (since no sequential 
were used in the first place, I believe).

> - query writing : your query is optimally written (looks fine here)

Thanks, this improves my confidence that I didn't miss anything obvious 
on the query level. :-)

> - data model (detailed after)

I add some more details below, but first I want to add some more
context: typically we deal with individual nodes, stored in a separate
treenode table (the result of the PostGIS query is then joined into it).
Each node has a parent (except the root node) and the edge between them
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.

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.

>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 

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.

> * Scenario A.1: your usual slices affects a majority of trees
>  - you don't really need an index on X,Y, but mainly on Z.
>    The idea is to separate your index into index(X,Y) + index(Z)
>     Try adding an index on GIST(numrange(ST_Zmin,ST_ZMax)) (you could add
>it anyway)
>     and an index on GIST(geom) (aka only 2D).

As long as I can also query from orthogonal perspectives with reasonable 
performance, I don't need separate X and Y indices, it seems.

First, I adding the following indices:

  CREATE INDEX treenode_edge_z_range_gist ON treenode_edge
  USING GIST(numrange(ST_ZMin(edge)::numeric, ST_ZMax(edge)::numeric), '[]');

  CREATE INDEX treenode_edge_2d_gist ON treenode_edge USING GIST(edge);

Without the type casts, Postgres wouldn't allow it and inclusive bounds
are the current behavior. My regular query of course doesn't make use of
this index in its current form, but replacing the &&& operator test with
a 2D bound box check plus a Z index condition turns out to be really

  FROM treenode_edge te
  WHERE te.edge && ST_MakeEnvelope(-537284.0, -340912.0, 1456444.0, 912464.0)
    AND numrange(ST_ZMin(te.edge)::numeric, ST_ZMax(te.edge)::numeric) && '[131670.0,131705.0]'::numrange
    AND ST_3DDWithin(te.edge, ST_MakePolygon(ST_GeomFromText(
      'LINESTRING( -537284.0 -340912.0 131687.5, 1456444.0 -340912.0 131687.5,
                   1456444.0  912464.0 131687.5, -537284.0  912464.0 131687.5,
                   -537284.0 -340912.0 131687.5)')), 17.5)
    AND te.project_id = 1;
  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)
    Recheck Cond: (numrange((st_zmin((edge)::box3d))::numeric, (st_zmax((edge)::box3d))::numeric, '[]'::text) && '[131670.0,131705.0]'::numrange)
    Filter: ((edge && '0103000000010000000500000000000000886520C100000000C0CE14C100000000886520C100000000A0D82B41000000003C39364100000000A0D82B41000000003C39364100000000C0CE14C100000000886520C100000000C0CE14C1'::geometry) AND (edge && '0103000080010000000500000000000000AB6520C10000000006CF14C100000000B012004100000000AB6520C100000000C3D82B4100000000B0120041000000804D39364100000000C3D82B4100000000C8130041000000804D3936410000000006CF14C100000000C813004100000000AB6520C10000000006CF14C100000000B0120041'::geometry) AND (project_id = 1) AND ('0103000080010000000500000000000000886520C100000000C0CE14C1000000003C130041000000003C39364100000000C0CE14C1000000003C130041000000003C39364100000000A0D82B41000000003C13004100000000886520C100000000A0D82B41000000003C13004100000000886520C100000000C0CE14C1000000003C130041'::geometry && st_expand(edge, '17.5'::double precision)) AND _st_3ddwithin(edge, '0103000080010000000500000000000000886520C100000000C0CE14C1000000003C130041000000003C39364100000000C0CE14C1000000003C130041000000003C39364100000000A0D82B41000000003C13004100000000886520C100000000A0D82B41000000003C13004100000000886520C100000000C0CE14C1000000003C130041'::geometry, '17.5'::double precision))
    Rows Removed by Filter: 26
    Heap Blocks: exact=4120
    Buffers: shared hit=4428
    ->  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)
          Index Cond: (numrange((st_zmin((edge)::box3d))::numeric, (st_zmax((edge)::box3d))::numeric, '[]'::text) && '[131670.0,131705.0]'::numrange)
          Buffers: shared hit=308
  Planning time: 0.856 ms
  Execution time: 33.731 ms

This is a very impressive speed-up, thanks for this suggestion! For all
examples I tried the results were exactly the same as with the &&&
operator. I only stumped across the fact the current semantic is to
have inclusive bounds, but numrange()'s behavior is '[)' by default.  
Being able to change this is also a nice feature.

I also tested creating range indices for X and Y so that I could replace
the &&& operator with three separate range checks. This, however, was
much slower with >2500ms. This setup would have been nicer, because it 
could work equally well for "slice" queries in each dimension.

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.

> * Scenario A.2 : your usual slices affect few trees
>  - you could easily perform a 2 part query, where the first part use an
>index on tree bounding box, and the second perform the query on edge. This
>require to create and maintain a tree_bounding_box table with appropriate
>indexes, which is quite easy with triggers Ex (not optimised):
>   WITH filtering_by_tree AS (
>     SELECT distinct tree_id
>     FROM my_sync_table_of_tree_bounding_boxes as tree
>     WHERE tree.bbox &&& your_slice
>   ) SELECT edge_id
>     FROM treenode_edge  AS te
>     WHERE te.geom &&& your_slice
>      AND EXISTS (SELECT 1 FROM filtering_by_tree As ft WHERE ft.tree_id =
>  - you could use pgpointcloud to similar effect by construction a patch
>per tree.

This is certainly also an appealing idea. I quickly tested this by
creating a neuron summary table with its current bounding boxes as box3d
and then queried in a similar way like you suggested above. I couldn't
get this very fast though, regardless of whether &&& or a 2D-approach is
in use. I typically ended up with about 500ms.

>Now something is very ambiguous in your data, and it could have major
>consequences :
> * Are the edges between nodes that represent neuron __only logical__ (they
>are straight lines, always the same Zmax-Zmin, and represent a conceptual
>graph), or are they polylines that represent the __actual__ neurone path?

Edges reflect the morphology of the actual neuron and their length
(even only in Z) can vary a lot.

> * Scenario B.1.1 : edges are only logical and have a constant Zmax-Zmin,
>i.e. nodes have quantified Z :
>  - you don't need edges at all for indexing, you should work only on nodes
>(i.e. index nodes IN Z + (X,Y), find nodes, then find relevant edges). The
>great advantage in this case is that you can use the full power of
>pgpointcloud to create meaningfull patches (of variables size for instance,
>see this
><http://www.sciencedirect.com/science/article/pii/S092427161630123X>, 3.1,
>Adapt patch grouping rules ), which means scalling in the range of billions
>of points possible.

Thanks, that's good to know, but our Zmax-Zmin can vary.

>  - you could try to think of using the dual of your current graph (edges
>become nodes, nodes become edges).

Why could this be an advantage here? I guess if there are helpful data 
structures for a dual representation, search time can be improved.

> * Scenario B.1.2 : edges are only logical, but nodes don't have a constant
>Zmax-Zmin. : you could insert (virtual) nodes so you are again in the same
>case as scenario B.1.1

I would prefer not doing this. We did this in the past, but wanted to
reduce the number of nodes needed and improve usability with large distances.

> * Scenario B.2 : edges are actual neurone path.
>  - you could construct intermediary objects (groups of neurone) of
>variable size (those within cubes), so as to have again a 2 steps query :
>first look for groups of neuron, then for neurons inside

This sounds similar to the grid I have implemented at the moment. So far
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 

> * Scenario C : your data don't have obvious sparsity.
>  - use the same strategy you currently use, but scale by partitionning
>your data into numerous tables.
>   You can see the partitionning entry of the postgres manual.
>   This would be much easier to create and maintain if you could partition
>on only Z dimension.
>   Note that you will need to duplicate (and then de-duplicate) some edges.

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.

>It's hard to recommend anything without further information on your
>data/your requirement/your ressources.

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
currently testing the promising looking 2D based approach (2D GIST and 
range GIST) further.

Thanks again, also for the nice explanatory graphics!
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 833 bytes
Desc: not available
URL: <http://lists.osgeo.org/pipermail/postgis-users/attachments/20170120/5a320c94/attachment.sig>

More information about the postgis-users mailing list