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.
--Chris On Tue, Mar 11, 2014 at 2:44 PM, Schregle Roland HSLU T&A <[email protected]> wrote: > Hi all, > > I'm new to PDL and semi-proficient with Perl at best. I'm posting this a > last resort after *lots* of frustration while trying to combine a set of RGB > images of dimensionality [3, N, M] using a selective averaging technique > often used in astrophotography called sigma clip (see > http://www.ccdware.com/support/presentations/advanced_image_combine_gralak20081116.pdf). > > The images are generated iteratively and the set size is not known in > advance. > > I append the images into a 4D array on every iteration using glue(). Of > course I'm open to alternatives... > > I'm seriously stuck on actually combining the [3, N, M, n] image stack after > n iterations using sigma clipping. This requires getting the [N, M] mean and > standard deviation (sigma) from the pixel luminance (average over RGB), then > averaging the stacked images using only those pixels whose deviation from > the mean is within some factor of sigma. This effectively rejects outliers. > > Here's what I have so far... > > ---------------------------------------------------------------------------------------------------- > use constant SIGMACLIP => 0.5; > my $imgStack; > > do { > ... > # Append new image $img to stack, adding 4th dimension > $imgStack = defined $imgStack ? $imgStack -> glue(3, $img) : $img; > ... 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. > # Get mean and std deviation for RGB average (= luminance) of all images > my $imgStackLum = average($imgStack); > # [N, M, n] --> [n, N, M] > my ($mean, undef, undef, undef, undef, undef, $sigma) = > statsover($imgStackLum -> reorder(2,0,1)); I find it more readable to use mv() rather than reorder() especially when you are bringing an axis from the end to the beginning of the dimensions since you can use -1 instead of the actual value. An array slice assignment clarifies the LHS to capture the desired 2 outputs from statsover: ($mean,$sigma) = ( $imgStackLum->mv(-1,0)->statsover )[0,6]; > # Add 3rd dim to mean for subtraction > my $imgStackDiff = $imgStackLum - $mean -> dummy(2); > > # Get mask array containing non-clipped pixels > $sigma *= SIGMACLIP; > my $mask = (abs($imgStackDiff) < $sigma); > > # Selectively combine $imgStack into selectively averaged [3, N, M] image > based on $mask > $imgComb = where($imgStack, $mask)... -> average; ??? 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); NOTE: I used PDL::NiceSlice syntax for the $imgTotal calculation to get the dummy dimensions on the $mask piddle. The above is equivalent to $mask->dummy(0,3); > ... > } until (converged); > ---------------------------------------------------------------------------------------------------- > > Note that I have to redo this for every iteration, which is anything but > efficient. Again, I'm open to suggestions for optimisation... once it's > working. 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). > This is all pretty standard stuff, so nothing exotic, and I'll bet > somebody's done it in PDL. I couldn't find any hints in the archives tho. > (Apparently the search function in the recent archives is broken, btw). > > Any help appreciated. Many thanks & best regards, 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. > --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 _______________________________________________ Perldl mailing list [email protected] http://mailman.jach.hawaii.edu/mailman/listinfo/perldl
