3 ms·

Need more details on what you're doing, but many simple operations used in satellite image analysis are embarrassingly parallel. You could, for example, read in

by ovis 11y ago

Need more details on what you're doing, but many simple operations used in satellite image analysis are embarrassingly parallel. You could, for example, read in chunks or strips using GDAL, distribute them across multiple threads or cores, and get some speed-up. If much of the time is spent in numpy rather than Python, I think it should be non-blocking so the GIL won't get in your way.

porker 11y ago
It's a friend's research project I'm tinkering with, so if I don't sound like I know what I'm doing I don't... :) The first step is masking out clouds (morphological erosion + more), and then applying a haze optimized transform to the image data. I'm confused about making code parallelisable when the transform of a pixel relies on neighbouring pixels - how do you remove boundary errors (take extra pixels, then crop off each strip before recombining I guess)? There must be something in-place to deal with boundary errors round the edge of the image, but so far I haven't grasped it! Numpy, Scipy, and particularly scikit-image are where the time is spent - my dabbling was originally to try and speed things up, but I've also become interested in understanding what is going on.
- ovis 11y ago
  Yeah, when I need neighbouring pixels I extract overlapping chunks. As far as the edge of the scene - there aren't any perfect solutions that I know of. You could try to fudge it by modifying the algorithm to handle the edges with fewer pixels, or perhaps use data from different scenes as an approximation. Or you may simply lose the data there.
- gh02t 11y ago
  This is where you get into the nitty-gritty (and stupidly difficult) type of parallel algorithms. When you have e.g. cells that are dependent on their neighbors you can no longer use the trivial divide-and-conquer strategy. Usually the strategy is to break the domain into chunks, with each chunk being multiple cells. Inside the chunk, the processor already has the data for adjacent cells so it can just read it. On the boundaries where some of the adjacent cells are in another chunk and on a different processor, it has to ask for the data... the vast majority of scientific codes here use OpenMP (not so bad, shared memory) or MPI (absolutely miserable to use IMO, but powerful distributed memory) to communicate. There's other tools for the job that are nicer than MPI, but they don't get much attention from HPC people (who tend to be very conservative with new technology, especially "high-level" abstractions). Inter-processor communication is orders of magnitude slower than querying local memory (especially if that processor is say on another machine), so on top of that there's lots of clever methods of what is essentially figuring out how to slice up the problem domain; examples include Red-Black schemes, domain decomposition and Wavefront propagation. To answer your question specifically on boundary conditions (and here I'm talking about the edges of the whole image where there's literally nothing in some cells nearby, not boundaries between partitioned regions like I was above)... you just have to choose something. There are three common approaches, which I'll illustrate by imagining you're doing a blur that decides how to blur a pixel by looking at the 8 pixels that surround it and that you're on the far right edge of the image (so the right three neighboring pixels would lie outside the image): 1. Zero - treat the three pixels as black/zero when calculating. Equivalently, treat them as whatever arbitrary value you find works best (e.g. white/1). 2. Reflective - treat the three pixels as whatever value the pixels on the left side at the appropriate positions would be. That is to say, wrap around to the other side. 3. Symmetric - treat the right three pixels as if they were the same as the left three (or the center three if you like). This is a common problem also when solving something like a differential equation on a discretized spatial grid. In this case, the edge values are determined based on your boundary conditions- method 1 corresponds to vacuum boundaries, 2. corresponds to periodic boundaries and 3. is a specific type of extrapolated boundary.

by ovis 11y ago

porker 11y ago

It's a friend's research project I'm tinkering with, so if I don't sound like I know what I'm doing I don't... :) The first step is masking out clouds (morphological erosion + more), and then applying a haze optimized transform to the image data. I'm confused about making code parallelisable when the transform of a pixel relies on neighbouring pixels - how do you remove boundary errors (take extra pixels, then crop off each strip before recombining I guess)? There must be something in-place to deal with boundary errors round the edge of the image, but so far I haven't grasped it! Numpy, Scipy, and particularly scikit-image are where the time is spent - my dabbling was originally to try and speed things up, but I've also become interested in understanding what is going on.

ovis 11y ago
Yeah, when I need neighbouring pixels I extract overlapping chunks. As far as the edge of the scene - there aren't any perfect solutions that I know of. You could try to fudge it by modifying the algorithm to handle the edges with fewer pixels, or perhaps use data from different scenes as an approximation. Or you may simply lose the data there.
gh02t 11y ago
This is where you get into the nitty-gritty (and stupidly difficult) type of parallel algorithms. When you have e.g. cells that are dependent on their neighbors you can no longer use the trivial divide-and-conquer strategy. Usually the strategy is to break the domain into chunks, with each chunk being multiple cells. Inside the chunk, the processor already has the data for adjacent cells so it can just read it. On the boundaries where some of the adjacent cells are in another chunk and on a different processor, it has to ask for the data... the vast majority of scientific codes here use OpenMP (not so bad, shared memory) or MPI (absolutely miserable to use IMO, but powerful distributed memory) to communicate. There's other tools for the job that are nicer than MPI, but they don't get much attention from HPC people (who tend to be very conservative with new technology, especially "high-level" abstractions). Inter-processor communication is orders of magnitude slower than querying local memory (especially if that processor is say on another machine), so on top of that there's lots of clever methods of what is essentially figuring out how to slice up the problem domain; examples include Red-Black schemes, domain decomposition and Wavefront propagation. To answer your question specifically on boundary conditions (and here I'm talking about the edges of the whole image where there's literally nothing in some cells nearby, not boundaries between partitioned regions like I was above)... you just have to choose something. There are three common approaches, which I'll illustrate by imagining you're doing a blur that decides how to blur a pixel by looking at the 8 pixels that surround it and that you're on the far right edge of the image (so the right three neighboring pixels would lie outside the image): 1. Zero - treat the three pixels as black/zero when calculating. Equivalently, treat them as whatever arbitrary value you find works best (e.g. white/1). 2. Reflective - treat the three pixels as whatever value the pixels on the left side at the appropriate positions would be. That is to say, wrap around to the other side. 3. Symmetric - treat the right three pixels as if they were the same as the left three (or the center three if you like). This is a common problem also when solving something like a differential equation on a discretized spatial grid. In this case, the edge values are determined based on your boundary conditions- method 1 corresponds to vacuum boundaries, 2. corresponds to periodic boundaries and 3. is a specific type of extrapolated boundary.