Forgot to CC this to the list... hope it comes in handy.
--Roland ------- Forwarded message ------- From: "Schregle Roland HSLU T&A" <[email protected]> To: "Chris Marshall" <[email protected]> Cc: Subject: Re: [Perldl] Implementing sigma clipping on RGB image stack Date: Wed, 12 Mar 2014 15:33:03 +0100 On Tue, 11 Mar 2014 23:28:20 +0100, Chris Marshall <[email protected]> wrote:
Hi Roland- It looks like you've already gotten everything working but the final selective average. I've put comments inline with your text below on ways I would do things for efficiency if performance (speed or memory) were important.
Hi Chris, megathanks up front for your feedback -- you just catapulted our project forward not insignificantly. ;^)
Since glue() and append() all reallocate the entire pdl object to add one entry, it is more efficient to allocate by bigger chunks to reduce the overhead. Then just use .= slice assignment to add new frame data to the stack.
The images are generated in p parallel threads, so I could allocate a [3,M,N,p]-dim PDL and drop them in there as the threads complete, then glue/append them to the stack in one go before averaging.
An array slice assignment clarifies the LHS to capture the desired 2 outputs from statsover: ($mean,$sigma) = ( $imgStackLum->mv(-1,0)->statsover )[0,6];
How efficient is it to use statsover() here, considering I just need the mean and stddev and ignore the remaining output? Does it actually skip calculating the unwanted stats or are they simply thrown away? Would it be faster to get mean and stddev manually?
You could use bad values to perform the computation but another approach is to calculated the average by hand. The sum of the values is: $imgTotal = ( $imgStack * $mask(*3) )->mv(-1,0)->sumover; and the number of values is just the same sum but just over the mask with one subtlety: where the N==1 you need to avoid a division: $imgNum = $mask->mv(-1,0)->sumover; And the resulting average over the participating pixels is then: $imgComb = $imgTotal / $imgNum->lclip(1);
This is great -- and it works! I was tempted to use average() after multiplying with $mask, until I realised it uses uniform rather than per-pixel weights; needless to say, the results were even noisier than without sigma clip! ;^)
The usual things for iterative calculation of statistics is to keep the raw sum and sum-of-squares totals so you can update the values for each pixel in O(1) time rather than O(n+1).
Updating a raw sum for the mean is straightforward, but I can't see a way around (mostly) redoing the stddev every time. Similarly, $imgTotal and $imgNum depend on $mask, which also changes with every iteration (albeit less frequently as more images are added).
I'm assuming (hoping) that you are using a current version of PDL but the above should be correct for the last several releases. Maybe someone else will post with the example using bad values.
I'm using PDL v2.007 with Perl 5.10 and updated all the relevant modules via CPAN. Your suggestions work fine! Craig posted an example using bad values; will look into that next. Many thanks again, much obliged! Regards, --Roland -- Dr. Roland Schregle Senior Research Associate T direct: +41 41 349 39 77 [email protected] Lucerne University of Applied Sciences and Arts School of Engineering and Architecture CC Envelopes and Solar Energy (EASE) Technikumstrasse 21, CH-6048 Horw T +41 41 349 33 11, F +41 41 349 39 60 www.hslu.ch/ccease _______________________________________________ Perldl mailing list [email protected] http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
