Home - General / All posts - Grid maths in M9 - Transfering lidar texture to a coarser DEM
 danb1,697 post(s) #30-Jul-18 03:23 I have a task which I really want to undertake in Manifold 9 due to the volume of data involved. Unfortunately however it will involve using grid maths and conditional logic to manipulate and combine grids (as would be achieved via Surface > Transform in Mfd8).I am attempting to apply the surface ‘texture’ of more detailed, but geographically restricted lidar data to a coarser but geographically extensive DEM. This might sound like an odd thing to wish to do, but the additional information encoded in the lidar ‘texture’ is very useful for fine guidance of covariates used in digital soil mapping of low lying areas. To this end, I have created a simple proof of concept in Mfd8 which I have outlined below, but ideally would like to complete the process in Mfd9 due to the volume of lidar data rather than processing tile by tile in Mfd8. Is this sort of processing currently available in Mfd9 given that while the two grids are expressed in the same projection, the pixel dimensions and offsets are different?BTW Merge in Mfd9 won’t work for me in this case due to height differences between the DEM and Lidar data creating a large height step along the join line.Mfd8 Processing steps (please feel free to make any suggestions for improvements as this is really just a proof of concept):Harvest the lidar texture (from [LIDAR_DEM])This I have done as simply as possible, by generating a highly smoothed 'average' lidar surface using the median pixel value over an arbitrary wide window here 5 pixels deep (11x11 window). I subtract this 'average' surface from the unmodified lidar to retain the texture without the heights.[LIDAR_DEM] - MedValue([LIDAR_DEM], 5)The resulting grid is called [LIDAR_TEXTURE].Combine the lidar texture and the coarse DEM data ([BASE_DEM])Combination is achieved using conditional statements. If the lidar texture is null (here -9999), use the coarser base dataelse use the coarse DEM plus the lidar texture. Because the range of values in the texture is small, this doesn't really change the overall height of the base DEM, rather only imposes the local texture which is more important than the height.Combination by this method will also not create a significant step in heights across the lidar/coarse DEM boundary.IIf([LIDAR_TEXTURE] = -9999, [BASE_DEM], [BASE_DEM] + [LIDAR_TEXTURE])The resulting grid is called [BASE RESULT].The result of this conditional statement does leave some pixel artefacts along the join line between the datasets (not sure why). I have smoothed these using something like:[BASE RESULT FINAL]:IIf([BASE RESULT]>5000, MedValue([BASE RESULT], 4), [BASE RESULT])Attachments: Clipboard-1.jpg Landsystems Ltd ... Know your land | www.landsystems.co.nz
 adamw8,568 post(s) #01-Aug-18 13:46 It all looks fairly simple except for the first step with the median filter which is doable, but more complex than it needs to be because we don't have the median filter as a built-in function (coming).Step 1: Do median filter on the more detailed but covering a smaller area LIDAR DEM. This currently requires a script. Tell me if you want it. Write the result into the second tile field in the same table as the original image.Step 2: Subtract the result of the median filter from the LIDAR DEM, leaving only the texture. This is just a Subtract transform between the tile fields. If we were writing the result of the median filter into a separate image, this would have required a join between the two images on X and Y.Step 3: Fill invisible pixels in the texture with 0. This is me getting rid of IIf - let's just convert it into an addition. The filling of invisible pixels is done using TileCombine(tile, TileMakeNew(cx, cy, 0)), invisible pixels in the first tile get filled from the second tile.Step 4: Project the texture to the exact coordinate system of the less detailed but covering a bigger area base DEM, using the exact same tile size as the base DEM. The coordinate system, including offsets and scales, and the tile size have to be exactly the same, this will make sure the X and Y values in the base DEM and the projected texture correspond to the same tiles.Step 5: Add two images together. This requires a join on X and Y. Adding tiles is done using just tile1+tile2.Step 6: Smooth the result using the median filter again.We are planning to make it easier to combine images dynamically automatically projecting them during transforms, this should make a fair bit of the above automatic.
 danb1,697 post(s) #04-Aug-18 00:34 Apologies for the delay in responding. Many thanks for this Adam, I suspected that the infrastructure would be largely there but not the GUI implementation. Very glad also to hear that it is coming as I believe M9 will ultimately transform much of my work with lidar and point clouds.With regards your kind offer for a median filter script, if a builtin filter is coming, I don't wish to distract you, though I am interested if it will demonstrate a way to generate custom filters further down the track. Landsystems Ltd ... Know your land | www.landsystems.co.nz
 adamw8,568 post(s) #04-Aug-18 08:57 Let's wait until the built-in filter then (yes, it is coming). We will provide examples of custom filters regardless.
 tjhb8,743 post(s) #13-Aug-18 02:54 I'm going to have a decent go at writing a median filter in the next week or so. I need to learn to write Manifold 9 filters in a standardised way for my own purposes.I would be really grateful if Dan could test the particular filter, perhaps against Manifold 8; and if Adam could cast an eye over the eventual code and suggest areas for improvement (including as to speed and efficient data access).
 danb1,697 post(s) #13-Aug-18 05:05 That would be great Tim and a really useful learning exercise for me. I look forward to testing your filter in due course. Landsystems Ltd ... Know your land | www.landsystems.co.nz
 adamw8,568 post(s) #13-Aug-18 07:29 Of course I'll cast an eye, I love things like this.
 tjhb8,743 post(s) #13-Aug-18 08:45 Thank you both. This will be fun.
 tjhb8,743 post(s) #27-Aug-18 07:15 P.s. I have been wondering whether the Manifold 8 transforms (MedianSquare, MedianSquare5 and MedValue; and possibly MedianCross which I haven't tested here) in fact use a weighted median filter, with weights according to radial distance or value difference from the central pixel--or both together, treating each pixel pair as a 3D vector. (The central pixel itself would need special treatment in all these cases.)I haven't tried that.I would like to be able to reproduce whatever 8 does, because from experience, it is very good.
 adamw8,568 post(s) #27-Aug-18 15:21 The code looks fine except it does not handle invisible pixels (as you say in the post and as you note in one of the TODOs in the script).You can rework the innermost loop to keep the number of valid values in 'window' (eg, put new values into it using something like 'window[valueCount++] = ...'), then sort just the valid values in the array (System.Array.Sort allows specifying the range to sort) and pick the median from there.I re-checked and 8 does not do any weighting, it just skips invisible pixels. I suspect this is the reason for the differences, at least on edges.Added: are you saying the results of the above function differ from output of M8 mid-image with no invisible pixels in sight? If so, I'll check closer.
 adamw8,568 post(s) #27-Aug-18 15:33 I think I know where the difference comes from: the index of the final value should be pixels/2, not (pixels+1)/2, because it is zero-based. Ie, for a radius of 1 = 3x3 window = 9 pixels, with all pixels visible, the zero-based index of the median value should be 4, not 5. Because pixels are: 0-1-2-3---4---5-6-7-8.
 tjhb8,743 post(s) #27-Aug-18 16:03 Thanks! I bet I would have compounded my error when adding invisible pixels to the mix. I didn't spot it.I'll do a fresh comparison after fixing that and confirm I get the same results as 8 now, then add in the pixel masks.By the way using outer and inner functions in the script was supposedly to help in reusing the same code/structure for other data types and other filters. That might change as I learn more.Very kind of you to check Adam, thanks again.Oh, and I also have to thank you for ending my irrational fear of spiders C#.
 tjhb8,743 post(s) #27-Aug-18 23:29 Confirming that Adam's fix has fixed it. Same result as Manifold 8--almost identical.The only differences now are in (some) marginal pixels. I am not at all inclined to worry about this. For anyone interested see attached image. The maximum (absolute) difference is just short of 10 units (metres). Roughly 2/3 of differences are > 1 unit (selected). I can reproduce if necessary.This is using the above query and script (with median derived per Adam's correction), window/radius 1, against transform MedianSquare in Manifold 8.Different pixels filtered in Manifold 8 using transform "Abs([ML 2] - [Median])" then drawn as points in SQL8.Otherwise identical.By the way, query+script in Manifold 9 is noticeably faster than prep+CUDA in Manifold 8.Attachments: diff image.png
 tjhb8,743 post(s) #28-Aug-18 00:28 In case anyone is still interested in the detail (I barely am myself):The only difference between MedianSquare(s) and MedValue(s, 1) in Manifold 8 is in their treatment of marginal/edge pixels. Otherwise their results are identical.Since marginal pixel results are nominal for window functions (and normally should be ignored), it might have been better if only one of those functions had existed in 8. Failing that, that they were literally identical. I expect that the same applies for MedianSquare5(s) and MedValue(s, 2), but I haven't tested.But this doesn't matter now.To repeat, there is absolutely no reason for concern here. The different functions handle (redundant) edge pixels slightly differently, but that has no analytic importance.Finally, the Manifold 9 treatment above is exactly like MedValue(s, w) in Manifold 8. (And also exactly like MedianSquare or MedianSquare5, if edge pixels are ignored.)
 adamw8,568 post(s) #28-Aug-18 06:43 By the way, query+script in Manifold 9 is noticeably faster than prep+CUDA in Manifold 8.This is because the window is small (so computations per pixel aren't big and GPGPU is of limited help) and because 9 uses multiple threads and 8 does not. If the window size parameter could grow higher than 5x5, GPGPU with a single thread would have eventually outperformed CPU with four threads. We are going to add median filter functions, including GPGPU variants, with extensions, as a built-in option to 9. But it is good that the performance of a query+script is pretty competitive. There will always be some function that we don't have built-in, being able to write a custom function and have it perform reasonably well is great.
 tjhb8,743 post(s) #27-Aug-18 23:35 I will post a version allowing for invisible pixels shortly (with perhaps other changes, taking more of Adam's advice).