Downsampling van zeer grote plaatjes

[DPLuS]

Ik ben in het bezit van satellietfoto's van een gezamenlijke grootte van 500.000 x 500.000 pixels.
Om een soort Google Maps kloon na te kunnen bouwen wil ik dit kaartmateriaal in verschillende zoom-lagen beschikbaar hebben.
Die bovengenoemde afmetingen zijn dan het hoogste zoomniveau.
Ik wil nu 10 lagen creëren die de lagere zoom-niveaus gaan vertegenwoordigen.

Op het laagste niveau moet een plaatje overblijven van ongeveer 5000 x 5000 pixels.

Alleen ben ik nu op zoek naar een algoritme dat me gaat helpen om dit voor deze gigantische afmetingen voorelkaar te krijgen.
Ik kan namelijk niet zomaar wat pixels weggooien bij het berekenen en samenstellen van de hogere zoomniveaus omdat het uiteindelijke plaatje er dan niet uit zal komen te zien.
Ik moet zoiets als dat ze in Photoshop doen bij het downsamplen toepassen, namelijk bicubic resampling.

Ik ben erg bang dat het allemaal te rekenintensief wordt en dat ik meer dan 2 weken moet gaan wachten voordat eindelijk 1 enkele laag berekend is.

Daarom zoek ik naar een algoritme dat als een soort "sliding window" over het plaatje verplaatst kan worden en daar zijn berekeningen doet (verkleinen met bicubic interpolation) en dat deel van het plaatje weer opslaat en verder gaat....

Wie heeft er tips of ervaringen met deze materie?

Alvast dank!

PS: Het bronmateriaal is in JPEG2000, opgeknipt in 3350 tegels van 8000x8000 pixels; voor wie het wil weten...

[leuk_he]

Wat heb je zelf al aan library/tools gevonden.

b.v. zoiets?

http://wiki.panotools.org/

gebruikt voor deze geloof ik:

[google=Largest Digital Photograph in the World ]

[Janoz]

Aangezien je zelf de controle hebt over de zoomlevels zou ik die zo gunstig mogelijk uitzoeken. bicubic resampling is pas ingewikkeld wanneer het te resizen raster nogal afwijkt van het origineel. Dat is de reden waarom je voor de zoomlevels gewoon je resolutie telkens door 2 moet delen. Het bicubic resampling komt dan gewoon neer op het middelen van 4 pixels naar 1 nieuwe pixel.

[Zoijar]

Je kan het in blokken met een read-only rand van een halve kernel width ophakken. Dan kan je al die blokken parallel in CUDA op de GPU berekenen, en tussentijds wegschrijven. Dat zou heel snel moeten gaan; paralleliseerd geweldig, en veel disk I/O kan ook weer asynchroon parallel met de GPU. Ik denk dat een simpele GTX280 daar niet al te lang over doet.

[MSalters]

Als de blokken JPG zijn, dan kun je alles behalve de grondfrequentie weggooien. Dan heb je vrijwel instantaan (in elk geval zonder wiskunde) een 1:8 schaling. Dat scheelt dus een factor 64 in de hoeveelheid pixels die je daadwerkelijk wil samplen.

[DPLuS]

Die zoomlevels kan ik niet zelf uitzoeken, die moeten namelijk namelijk corresponderen met de rasterlagen die al bepaald zijn.
Delen door 2 is dus helaas geen optie.

@MSalters: De bronbestanden zijn ECW-bestanden, JPEG2000 gecomprimeerd.

[ Voor 18% gewijzigd door DPLuS op 16-01-2009 17:20 ]

[Domokoen]

Een eenvoudige manier om dit voor elkaar te krijgen is door een Gaussisch filter toe te passen met bepaalde grootte (sigma). En vervolgens met bilineare interpolatie te subsamplen. Door de interpolatie zit je niet vast aan een bepaald raster.

Je gaat wel tegen wat technische probleempjes oplopen bij de randen van je opgeknipte blokken, omdat het Gaussische filter "bij de buren" wil kijken. Maar als je het implementeert met een sliding window, dan kom je vrij ver.

[matthijsln]

DPLuS schreef op vrijdag 16 januari 2009 @ 15:47:
Ik ben in het bezit van satellietfoto's van een gezamenlijke grootte van 500.000 x 500.000 pixels.
Om een soort Google Maps kloon na te kunnen bouwen wil ik dit kaartmateriaal in verschillende zoom-lagen beschikbaar hebben.
Die bovengenoemde afmetingen zijn dan het hoogste zoomniveau.
Ik wil nu 10 lagen creëren die de lagere zoom-niveaus gaan vertegenwoordigen.

[...]
Wie heeft er tips of ervaringen met deze materie?

PS: Het bronmateriaal is in JPEG2000, opgeknipt in 3350 tegels van 8000x8000 pixels; voor wie het wil weten...

Dit heb ik - ik vermoed met exact dezelfde luchtfoto's - wel eens gedaan. Heel Nederland met 50cm per pixel zeker?

Naast verkleinen zou ik ook de kleinere tiles samenvoegen, want als je heel Nederland wilt bekijken en daarvoor duizenden bestanden moet openen werkt ook niet zo flitsend snel.

Ik heb voor het verkleinen en combineren zelf in Java een tool gemaakt (maar iets van 500 regels). Het verkleinen zelf deed ik gewoon met BufferedImage.getScaledInstance(w, h, Image.SCALE_AREA_AVERAGING). Voor verkleinen werkte averaging beter dan bilinear of bicubic. Ook samenvoegen gewoon met een BufferedImage, voor bestanden van 8k bij 8k werkt dat nog prima. Naast het combineren en resizen berekende mijn tooltje ook nieuwe coordinaten van de hoekpunten van de samengevoegde tile.

De rekentijd valt best mee, iets van een tiental uur ofzo.... Je kan het natuurlijk makkelijk verdelen over meerdere CPU's/computers. Voor een volgende laag die verder uitgezoomd is kan je de vorige laag gebruiken. De brontiles daarvoor gebruiken is niet zo heel zinvol. Daarmee spaar je al heel wat rekentijd uit.

BTW, als je voor je verkleinde plaatjes ook JPEG2000 gebruikt hoef je echt geen 10 lagen te maken (ik heb er zelf maar 3 gemaakt), want omdat JPEG2000 een wavelet compressieformaat is kan al heel efficient een plaatje van lagere resolutie gemaakt worden waarbij niet het gehele bestand hoeft te worden gedecomprimeerd.

Als samenvoegen niet nodig is kan je ook gdaladdo gebruiken.

[matthijsln]

DPLuS schreef op vrijdag 16 januari 2009 @ 17:18:
Die zoomlevels kan ik niet zelf uitzoeken, die moeten namelijk namelijk corresponderen met de rasterlagen die al bepaald zijn.
Delen door 2 is dus helaas geen optie.

@MSalters: De bronbestanden zijn ECW-bestanden, JPEG2000 gecomprimeerd.

Oh, ik lees nu pas dat het ECW bestanden zijn... Die zijn al heel erg goptimaliseerd voor het tonen op een lagere resolutie (wavelet). Daarvoor is het meestal niet nodig om "pyramids/mosaics" te maken. Wat natuurlijk wel kan op mijn eerder genoemde methode.

Daarnaast kan je kan ook Tilecache gebruiken bovenop bijvoorbeeld MapServer. Of Mapnik met mod_tile.

Of met de commerciele tool van de makers van het ECW formaat: ER Mapper Image Compressor.