Beefy Boxes and Bandwidth Generously Provided by pair Networks
Perl: the Markov chain saw
 
PerlMonks  

Map grid to boundary conversion algorithm

by Willard B. Trophy (Hermit)
on Apr 16, 2004 at 14:52 UTC ( [id://345780]=perlquestion: print w/replies, xml ) Need Help??

Willard B. Trophy has asked for the wisdom of the Perl Monks concerning the following question:

Despite what I said in 341327, I find myself prototyping in Perl. Fortran's still too weird (again).

I've got a grid of values, a tiny trivial example being:

P P B R B P R B B B B R R R B B R R R B B G G B B

(For the interested, these are terrain "clutter" (roughness) maps from satellite survey, as managed by Golden Software's Surfer mapping package.)

Unfortunately, the industry standard package WaSP doesn't accept grid input, but needs the locations of the boundaries where terrain types change, something like (with values):

P P|B|R|B +-+ +-+ P|R|B B B -+ +---+ B|R R R|B | | B|R R R|B +---+-+ B|G G|B B

or maybe clearer without the cell values, just showing boundaries:

| | | +-+ +-+ | | -+ +---+ | | | | | | +---+-+ | |

This looks to me like a multicolour version of Moore Neighbourhood Tracing. I've bounced some ideas around the office, but implementing them has proved problematic.

I probably haven't defined the problem very well, but ... how would you go about designing a routine to solve this?

UPDATE: Many thanks to all who responded. I should maybe have been a little less figurative with the ASCII art; it was the actual roughness contour lines I was after, though I'm sure I could write an ASCII to vector converter.

--
bowling trophy thieves, die!

Replies are listed 'Best First'.
Re: Map grid to boundary conversion algorithm
by kvale (Monsignor) on Apr 16, 2004 at 16:04 UTC
    Here is a solution to the problem:
    # read in matrix while (<DATA>) { next if /^\s*$/; push @r, [split]; } # find the boundaries for my $r (0..@r-1) { for my $c (0..@{$r[0]}-2) { $vert[$r][$c] = $r[$r][$c] ne $r[$r][$c+1] ? 1 : 0; } } for my $r (0..@r-2) { for my $c (0..@{$r[0]}-1) { $horz[$r][$c] = $r[$r][$c] ne $r[$r+1][$c] ? 1 : 0; } } for my $r (0..@r-2) { for my $c (0..@{$r[0]}-2) { $cross[$r][$c] = $vert[$r][$c] + $vert[$r+1][$c] + $horz[$r][$c] + $horz[$r][$c+1] > 1 ? 1 : 0; } } # Print it out for my $r (0..@r-1) { for my $c (0..@{$r[0]}-1) { print $r[$r][$c]; print $vert[$r][$c] ? '|' : ' ' if $c < @{$r[0]}-1; } print "\n"; next unless $r < @r-1; for my $c (0..@{$r[0]}-1) { print $horz[$r][$c] ? '-' : ' ' ; print $cross[$r][$c] ? '+' : ' ' if $c < @{$r[0]}-1; } print "\n"; } __DATA__ P P B R B P R B B B B R R R B B R R R B B G G B B
    I have simplified the printing to just put a cross at all intersections, you can easily generalize to vertical and horizontal connections, etc. I'll leave the conversion to Fortran 1-based matrices to you :-)

    Update: If efficiency is a concern, the three passes to find the boundary can be combined into one with a little extra logic. For large matrices, this change may blow the cache fewer times.

    -Mark

      Here is my stab. As kvale notes you can decrease the number of passes.

      #!/usr/bin/perl; use strict; use warnings; my @grid = ( [qw(P P B R B)], [qw(P R B B B)], [qw(B R R R B)], [qw(B R R R B)], [qw(B G G B B)], ); $" = ''; my @new; # pass 1 for my $row(0..$#grid) { my $last = $grid[$row][0]; push @{$new[$row*2]}, $grid[$row][0]; for my $pt(0..$#{$grid[$row]}) { next if $pt == 0; if ( $grid[$row][$pt] eq $last ) { push @{$new[$row*2]}, (' ', $grid[$row][$pt]); } else { push @{$new[$row*2]}, ('|', $grid[$row][$pt]); $last = $grid[$row][$pt]; } last if $pt == $#{$grid[$row]}; } last if $row == $#grid; my $width = $#{$grid[$row]}; $new[$row*2+1] = [(' ') x ($width * 2 + 1)]; } # pass 2 for ( my $row=2; $row < @new; $row+=2 ) { for ( my $pt=0; $pt < @{$new[$row]}; $pt +=2 ) { if ( $new[$row][$pt] eq $new[$row-2][$pt] ) { $new[$row-1][$pt] = ' '; } else { $new[$row-1][$pt] = '-'; } } } # find the corners for ( my $row=1; $row < @new; $row+=2 ) { for ( my $pt=1; $pt < @{$new[$row]}; $pt +=2 ) { # do we have corners ? if ( $new[$row][$pt-1].$new[$row-1][$pt] eq '-|' or $new[$row][$pt-1].$new[$row+1][$pt] eq '-|' or $new[$row][$pt+1].$new[$row-1][$pt] eq '-|' or $new[$row][$pt+1].$new[$row+1][$pt] eq '-|' ) { $new[$row][$pt] = '+'; } } } # final pass to complete now we have corners for ( my $row=1; $row < @new; $row+=2 ) { for ( my $pt=1; $pt < @{$new[$row]}; $pt +=2 ) { if ( $new[$row][$pt-1].$new[$row][$pt].$new[$row][$pt+1] eq '- + -' ) { $new[$row][$pt] = '-'; } elsif ( $new[$row-1][$pt].$new[$row][$pt].$new[$row+1][$pt] e +q '| |' ) { $new[$row][$pt] = '|'; } } } print 'Wanted: P P|B|R|B +-+ +-+ P|R|B B B -+ +---+ B|R R R|B | | B|R R R|B +---+-+ B|G G|B B Got: '; print "@$_\n" for @new;

      cheers

      tachyon

Re: Map grid to boundary conversion algorithm
by simonm (Vicar) on Apr 16, 2004 at 16:16 UTC

    Assuming there's no performance pressure, I think a straightforward search would be sufficient.

    Here's an untested code outline:

    use strict; # Build grid array my @grid; foreach ( <DATA> ) { s/^\s*\b// or next; push @grid, [ split ' ', $_ ], } # Assign arbitrary regions my @region; foreach my $x ( 0 .. $#grid ) { foreach my $y ( 0 .. $#{ $grid[0] } ) { $region[$x][$y] = "$x.$y"; } } # Merge neighboring regions foreach my $x ( 0 .. $#grid ) { foreach my $y ( 0 .. $#{ $grid[0] } ) { my @neighbors = ( [ $x-1, $y ], [$x, $y-1] ); foreach my $coord ( @neighbors ) { my ( $nx, $ny ) = @$coord; next if ( $nx < 0 or $ny < 0 ); if ( $grid[$x][$y] eq $grid[ $nx ][ $ny ] and $region[$x][$y] ne $region[ $nx ][ $ny ] ) { # Cell is the same color as its neighbor, join it merge_regions( $x, $y, $nx, $ny ) } } } } sub merge_regions { my ( $tx, $ty, $nx, $ny ) = @_; my $tr = $region[$tx][$ty]; my $nr = $region[$nx][$ny]; foreach my $x ( 0 .. $#grid ) { foreach my $y ( 0 .. $#{ $grid[0] } ) { if ( $region[$x][$y] eq $tr ) { $region[$x][$y] = $nr } } } } my $count = 0; my %regions; foreach my $x ( 0 .. $#grid ) { foreach my $y ( 0 .. $#{ $grid[0] } ) { my $num = ( $regions{ $region[$x][$y] } ||= ++ $count ); print $num . " "; } print "\n" } __DATA__ P P B R B P R B B B B R R R B B R R R B B G G B B

    The output needs to be adjusted to your application; for now it just prints a grid of contiguous regions:

    1 1 2 3 2 
    1 4 2 2 2 
    5 4 4 4 2 
    5 4 4 4 2 
    5 6 6 2 2 

    P.S. I've only seen "Moore Neighbourhood" used to refer to the immediate lateral and diagonal neighbors of a cell, but perhaps there's more than one meaning? Do you want to allow diagonal connections between neighbors, or only lateral ones?

    Update: Looking at kvale's solution, I noticed that I implicitly assumed that the OP was looking for "contiguous regions" rather than "boundaries" -- it works out to the same thing, but influenced my choice of approach.

      I might be wrong but to me this looks very similar as the input ... ?!

      pelagic
        I might be wrong but to me this looks very similar as the input

        Yes, but the regions have been assigned distinct numbers. (Looking closely, you'll see that the clump of "B"s in the top right corner is given a different number than the clump of "B"s in the bottom left.)

Re: Map grid to boundary conversion algorithm
by BrowserUk (Patriarch) on Apr 16, 2004 at 18:58 UTC

    A slightly different approach.

    #! perl -slw use strict; my $g = qq[PPBRB PRBBB BRRRB BRRRB BGGBB]; my $w = 5 -1; my( $o1, $o2 )= ( $w * 2, $w *2 +1 ); $g =~ s[\n]["\n" . '-+' x $w . "-\n"]eg; $g =~ s[([PRGB])(?=[PRGB])][$1|]g; $g =~ s[(?<=([PRGB]))\|(?=\1)] [ ]g; $g =~ s[(?<=(.).{$o2})[-](?=.{$o2}\1)] [ ]gs; $g =~ s[(?<= .{$o1}[-])\+(?=[-].{$o1} )] [-]gs; $g =~ s[(?<= .{$o1}[ -])\+(?=[ -].{$o1} )][ ]gs; $g =~ s[(?<=\|.{$o1} )\+(?= .{$o1}\|)] [|]gs; print $g; __END__ P:\test>345780 P P|B|R|B +-+ +-+ P|R|B B B -+ +---+ B|R R R|B | | B|R R R|B +---+-+ B|G G|B B

    Examine what is said, not who speaks.
    "Efficiency is intelligent laziness." -David Dunham
    "Think for yourself!" - Abigail
Re: Map grid to boundary conversion algorithm
by pelagic (Priest) on Apr 16, 2004 at 15:24 UTC
    Hi,
    I'd first define a 2-dimensional system to identify the terrain-data and to identify all distinct border pieces within the same terrain. You should be able of identifying all 4 borders of any terrain-piece. You should also be able of identifying the crossings in the same way.
    So if you start with a piece of terrain e.g. $t[x][y] you should have an algorithm to find its borders and border-crossings.
    One possibility might be to define that every piece of terrain has it's own borders and crossings and make the identification dependant like that: $t[x][y] would lead to $b[x][y][E], $b[x][y][S], $c[x][y][NE] and $c[x][y][SE].

    pelagic
Re: Map grid to boundary conversion algorithm
by jdporter (Paladin) on Apr 16, 2004 at 15:33 UTC
    The WAsP web site mentions a map conversion utility called Surfer. Does it not do what you need?

    jdporter
    The 6th Rule of Perl Club is -- There is no Rule #6.

      Thanks, JD, but this is output from Surfer I'm dealing with. WaSP doesn't read roughness data in Surfer format.

      --
      bowling trophy thieves, die!

Re: Map grid to boundary conversion algorithm
by BrowserUk (Patriarch) on Apr 16, 2004 at 19:31 UTC
    UPDATE: Many thanks to all who responded. I should maybe have been a little less figurative with the ASCII art; it was the actual roughness contour lines I was after, though I'm sure I could write an ASCII to vector converter.

    Perhaps you could give an example of the format you need, for your sample?


    Examine what is said, not who speaks.
    "Efficiency is intelligent laziness." -David Dunham
    "Think for yourself!" - Abigail
      Sorry for the delay. With your comments in mind, I ended up coding it this way:
      1. define a zero-value one cell border around the whole array.
      2. walk through each row horizontally, looking for vertical changes in value. Store each one in a structure containing the left and right values, and the coordinates of the end points
      3. walk through each column vertically, looking for horizontal changes in value, as above.
      4. search for adjoining vertical segments, and merge them.
      5. search for adjoining horizontal segments, and merge them.
      6. Write out the lists of line segments.

      It appears that the mapping tool we use is clever enough to join up adjacent lines into areas, so lines are all we need. The reason for the zero-value border isn't initially obvious, but it's necessary for the contouring routine to get its senses right.

      I put an example of what a tiny (300m x 300m) test file in my blog, if anyone cares to see what it looks like.

      I know my algorithm isn't that efficient; it took nearly 8 hours to process a real 73 x 58 km dataset on a P4-2800. My Fortran version (compiled with the excellent free-for-personal use Intel Fortran Compiler for Linux) looks like it will run at least an order of magnitude quicker.

      I'm beginning to warm to Fortran (90, or 95; F77 I cannot love). It has a nifty WHERE command which is somewhere between map and grep.

      Yes, at this point, you can say 'He's got away from us Jack.' ...

      --
      bowling trophy thieves, die!

        Glad you have a solution that works for you. It is an interesting problem that I think I would have enjoyed having a crack at, but it still isn't clear to me what format you want the output in? What would the line segments for your original example be?


        Examine what is said, not who speaks.
        "Efficiency is intelligent laziness." -David Dunham
        "Think for yourself!" - Abigail

Log In?
Username:
Password:

What's my password?
Create A New User
Domain Nodelet?
Node Status?
node history
Node Type: perlquestion [id://345780]
Approved by broquaint
Front-paged by Enlil
help
Chatterbox?
and the web crawler heard nothing...

How do I use this?Last hourOther CB clients
Other Users?
Others imbibing at the Monastery: (4)
As of 2024-03-29 15:51 GMT
Sections?
Information?
Find Nodes?
Leftovers?
    Voting Booth?

    No recent polls found