Using Fourrier-Mellin transform in PostGIS Part III
Intégration de la transformée de Fourier-Mellin dans PostgreSQL
Nous avons vu comment intégrer du code C++ dans PostgreSQL lorsque celui-ci est simplissime. Dans cette nouvelle série d’articles, nous allons intégrer le travail de @Thoduka2017.
Les codes dont nous partirons sont ici sur GitHub.
Quatre librairies sont nécessaires pour faire tourner le code référencé:
- fftw3 pour la partie mathématiques
- OpenCV pour la gestion des images
- Eigen pour les calculs matriciels
- Boost pour la gestion des fichiers
Je vous recommande d’installer les codes du Repo dans un premier temps avant de penser à les faire tourner par Postgres. Une vérification rapide de l’existence des libariries dans votre système peut-être effectuée via:
ldconfig -p|grep fftw
Objectif
Nous allons utliser cette librairie pour créer une fonction
St_AlignRaster(Orast,Drast)
qui aligne Orast sur Drast en utilisant l’algorithme de Reddy et Chatterji [@506761]. Après une brève présentation de l’extension Raster de PostGIS, nous verrons comment obtenir de PostgreSQL les pointeurs vers les rasters afin de les fournir à OpenCV dans un premier temps, puis à la libraire de Thoduka dans un second temps.
Présentation du type RASTER de PostGIS
Le modèle RASTER, dans sa version sérialisée est exposé ici. La lecture de ce document permet de comprendre certains points du code de l’extension Raster de PostGIS.
Deux structures permettent de gérer les rasters dans Postgres. La version désérialisée qui contient un pointeur vers les bandes et qui est appellée “rt_raster” dans le code:
struct rt_raster_t {
uint32_t size;
uint16_t version;
/* Number of bands, all share the same dimension
* and georeference */
uint16_t numBands;
/* Georeference (in projection units) */
double scaleX; /* pixel width */
double scaleY; /* pixel height */
double ipX; /* geo x ordinate of the corner of upper-left pixel */
double ipY; /* geo y ordinate of the corner of bottom-right pixel */
double skewX; /* skew about the X axis*/
double skewY; /* skew about the Y axis */
int32_t srid; /* spatial reference id */
uint16_t width; /* pixel columns - max 65535 */
uint16_t height; /* pixel rows - max 65535 */
rt_band *bands; /* actual bands */
};
typedef struct rt_raster_t* rt_raster
La version sérialisée est modélisée par une structure d’entête:
struct rt_raster_serialized_t {
/*---[ 8 byte boundary ]---{ */
uint32_t size; /* required by postgresql: 4 bytes */
uint16_t version; /* format version (this is version 0): 2 bytes */
uint16_t numBands; /* Number of bands: 2 bytes */
/* }---[ 8 byte boundary ]---{ */
double scaleX; /* pixel width: 8 bytes */
/* }---[ 8 byte boundary ]---{ */
double scaleY; /* pixel height: 8 bytes */
/* }---[ 8 byte boundary ]---{ */
double ipX; /* insertion point X: 8 bytes */
/* }---[ 8 byte boundary ]---{ */
double ipY; /* insertion point Y: 8 bytes */
/* }---[ 8 byte boundary ]---{ */
double skewX; /* skew about the X axis: 8 bytes */
/* }---[ 8 byte boundary ]---{ */
double skewY; /* skew about the Y axis: 8 bytes */
/* }---[ 8 byte boundary ]--- */
int32_t srid; /* Spatial reference id: 4 bytes */
uint16_t width; /* pixel columns: 2 bytes */
uint16_t height; /* pixel rows: 2 bytes */
};
Lors de la sérialisation cette structure précéde la donnée des bandes sur le modèle:
[HEADER] [BAND0] [BAND1] [BAND2]..
^Alignée ^Alignée ^Alignée
Le header est sur 8 octets, Chaque bande est alignée sur 8 octets. Le respect de l’alignement engendre quelques tracas mais finalement ceux-ci sont transparents lorsque nous cherchons à accéder à la donnée.
A l’usage, c’est la version désérialisée qui est utilisée. Nous aurons l’occasion, un peu plus loin, d’illustrer la façon de désérialiser un raster lorsqu’il est appellé par une fonction SQL.
Pour terminer avec les objets centraux de l’extension Raster, voici la structure rt_band_t (Noter le typedef: typedef struct rt_band_t rt_band*):
struct rt_band_t {
rt_pixtype pixtype;
int32_t offline;
uint16_t width;
uint16_t height;
int32_t hasnodata; /* a flag indicating if this band contains nodata values */
int32_t isnodata; /* a flag indicating if this band is filled only with
nodata values. flag CANNOT be TRUE if hasnodata is FALSE */
double nodataval; /* int will be converted ... */
int8_t ownsdata; /* 0, externally owned. 1, internally owned. only applies to data.mem */
rt_raster raster; /* reference to parent raster */
union {
void* mem; /* actual data, externally owned */
struct rt_extband_t offline;
} data;
};
Ce qui nous importe donc, c’est de récupérer un pointeur sur la donnée d’une bande afin, par exemple de le passer à une librairie comme OpenCV car in fine, pour Postgres, la donnée se trouve dans les bandes.
Voici une ébauche de code illustrant l’usage de la librairie Raster de PostGIS (PgCv2.c). Notre fonction reçoit le nom de deux rasters monobandes et de pixels codés sur un octet chacun, récupère un pointeur vers la donnée de chaque bande avant de les confier à notre Wrapper OpenCV:
PG_FUNCTION_INFO_V1(RASTER_asCV);
Datum RASTER_asCV(PG_FUNCTION_ARGS)
{
rt_pgraster *pgraster = NULL;
rt_raster raster;
rt_band band = NULL; /* Pointer to rt_band_t */
uint32_t numBands;
void* pband=NULL;
rt_pixtype px_type;
uint8_t type;
/* pgraster is null, return null */
if (PG_ARGISNULL(0)) PG_RETURN_NULL();
/* Get serialised raster by name*/
pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(0));
/* Deserialise pgraster into raster */
raster = rt_raster_deserialize(pgraster, FALSE);
if (!raster) {
PG_FREE_IF_COPY(pgraster, 0);
elog(ERROR, "RASTER_asCV: Could not deserialize raster");
PG_RETURN_NULL();
}
/* Getting the num. of bands:
uint16_t
rt_raster_get_num_bands(rt_raster raster) {
assert(NULL != raster);
return raster->numBands;
}
*/
numBands = rt_raster_get_num_bands(raster);
if(numbands! = 1)
{
PG_FREE_IF_COPY(pgraster, 0);
rt_raster_destroy(raster);
elog(ERROR, "RASTER_asCV: One band raster required");
PG_RETURN_NULL();
}
band = rt_raster_get_band(raster,0);
/* Here's the code of that last function:
rt_band
rt_raster_get_band(rt_raster raster, int n) {
assert(NULL != raster);
if (n >= raster->numBands || n < 0)
return NULL;
return raster->bands[n];
} */
px_type = band->pixtype;
int pixbytes = rt_pixtype_size(px_type);
pband = band->data.mem; /* Pointer to the band, to be transfered to OpenCV */
/* Second band */
pgraster = (rt_pgraster *) PG_DETOAST_DATUM(PG_GETARG_DATUM(1));
/* Deserialise pgraster into raster */
rt_raster To_rast = rt_raster_deserialize(pgraster, FALSE);
if (!To_rast) {
PG_FREE_IF_COPY(pgraster, 0);
elog(ERROR, "RASTER_asCV: Could not deserialize raster");
PG_RETURN_NULL();
}
numBands = rt_raster_get_num_bands(To_rast);
if(numbands! = 1)
{
PG_FREE_IF_COPY(pgraster, 0);
rt_raster_destroy(To_rast);
elog(ERROR, "RASTER_asCV: One band raster required");
PG_RETURN_NULL();
}
rt_band To_band = rt_raster_get_band(To_rast,0);
rt_pixtype To_px_type = To_band->pixtype;
int To_pixbytes = rt_pixtype_size(To_px_type);
void* To_pband = To_band->data.mem; /* Pointer to the band, to be transfered to OpenCV */
/*For this simple example, images have CV_8U pixel's type, which is defined as 0 */
type=0 /*CV_8UC1*/;
/* Now we can pass pband and To_pband (and size) to openCV's wrapped mat constructor */
uint_8* img=ST_register(raster->height,raster->width,type,pband,To_pband,raster->width * pixbytes);
/* Rest is for destroying pointers */
}
Lien PostgreSQL-OpenCV
Dans cette première étape, nous allons créer une fonction ST_Raster2CV(rast) dont l’objectif est juste de voir comment passer un pointeur d’une bande du raster à OpenCV pour l’afficher à l’écran.
Code initial OpenCV
Nous partons, là encore d’un code, le plus simple possible, trouvé sur le site d’OpenCV
//! [includes]
#include <opencv2/core.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/highgui.hpp>
#include <iostream>
using namespace cv;
//! [includes]
int main()
{
//! [imread]
std::string image_path = samples::findFile("sample1.png");
Mat img = imread(image_path, IMREAD_GRAYSCALE);
//! [imread]
//! [empty]
if(img.empty())
{
std::cout << "Could not read the image: " << image_path << std::endl;
return 1;
}
//! [empty]
//! [imshow]
imshow("Display window", img);
//! [imshow]
//! [imshow]
int k = waitKey(0); // Wait for a keystroke in the window
//! [imshow]
//! [release memory on keystroke]
if(k)
{
img.release() ;
}
//! [release]
return 0;
}
Makefile OpenCV
Le premier petit soucis concerne la compilation de ce code qui est proposée nativement via CMake. Voici un Makefile qui fait cependant le travail. Cela nous permettra de compiler ultérieurement OpenCV avec Postgres grâce à un seul Makefile.
CC = g++
CFLAGS = -g -Wall
SRCS = DisplayImage.cpp
PROG = DisplayImage
OPENCV = `pkg-config opencv4 --cflags --libs`
LIBS = $(OPENCV)
$(PROG):$(SRCS)
$(CC) $(CFLAGS) -o $(PROG) $(SRCS) $(LIBS)
Compatibilité du type Raster
Ce qui nous intéresse dans le code ci-desssus, c’est qu’un objet Mat, img, est créé qui reçoit les données d’un fichier JPEG. Or il est possible d’instancier un tel objet en lui passant un pointeur vers une image. La signature de ce constructeur est:
Mat (int rows, int cols, int type, void *data, size_t step=AUTO_STEP)
A ce stade, il est donc nécessaire de comprendre comment est organisée la donnée pointée par data.
De manière générale, OpenCV stocke les images dans un tableau multidimensionnel M tel que l’on passe d’une bande à une autre grâce aux valeurs d’un autre tableau, M.step[]. Les valeurs de ce tableau sont décroissantes.
Ainsi l’élément de coordonnées
(i_0,...,i_{M.dims-1})
où
0 \leq i_k<M.size[k]
a pour adresse:
addr(M_{i_0,...,i_{M.dims−1}})=M.data+M.step[0]∗i_0+M.step[1]∗i_1+...+M.step[M.dims−1]∗i_{M.dims−1}
En conséquence, une bande simple verra un pixel de coordonnée (i,j) stocké à l’adresse suivante:
addr(M_{i,j})=M.data+M.step[0]∗i+M.step[1]∗j
La décroissance de M.step implique que les pixels sont stockés ligne par ligne et donc M.step[1] est égal à la taille (en octets) d’un pixel.
Pour un raster bi-bande, les valeurs de chaque bande d’un pixel se suivent contrairement à ce qui se passe dans Postgres où, comme on l’a vu, les bandes sont stockées les unes après les autres. Ceci implique que le passage d’un raster multibande de Postgres à OpenCV aura un coût important puisqu’il faudra réorganiser le raster en mémoire.
Nous nous contenterons donc par la suite de travailler avec des rasters mono-canaux pour lesquels les accès seront directs.
Vers un usage d’OpenCV dans PostgreSQL
Une fois ces contraintes comprises, on peut commencer à drapper les codes dans du C comme indiqué dans notre premier article:
Nous créons CvWrapper.cc:
//! [includes]
#include <opencv2/core.hpp>
#include <opencv2/imgcodecs.hpp>
#include <opencv2/highgui.hpp>
#include <unistd.h>
#include <iostream>
#include "CvWrapper.h"
using namespace cv;
//! [includes]
extern "C" {
int Cv_W(){
//! [imread]
std::string image_path = samples::findFile("sample3.png");
Mat img = imread(image_path, IMREAD_GRAYSCALE);
//! [imread]
//! [empty]
if(img.empty())
{
std::cout << "Could not read the image: " << image_path << std::endl;
return 1;
}
//! [empty]
std::cout << "Could read the image: " << image_path << std::endl;
//! [imshow]
imshow("Display window", img);
//! [imshow]
int k = waitKey(0); // Wait for a keystroke in the window
//! [imshow]
//! [release memory on keystroke]
if(k)
{
img.release() ;
}
//! [release]
return 0;
}
}
Puis CvWrapper.h:
#ifndef __CVWRAPPER_H
#define __CVWRAPPER_H
#ifdef __cplusplus
extern "C" {
#endif
int Cv_W();
#ifdef __cplusplus
}
#endif
#endif
Enfin, MyMain_C.c pour tester le bon fonctionnement de ce wrappage:
#include "CvWrapper.h"
#include <stdio.h>
int main(int argc, char* argv[]){
int reponse;
reponse=Cv_W();
return reponse;
}
Reste à compiler et linker (je renvoie le lecteur au premier article de cette série pour les explications):
g++ -c -Wall -fPIC -o CvWrapper.o CvWrapper.cc `pkg-config opencv4 --cflags --libs`
g++ -shared -o libCvWrapper.so CvWrapper.o `pkg-config opencv4 --cflags --libs`
gcc -c MyMain_C.c -o MyMain_C.o
gcc MyMain_C.o -lCvWrapper -L. -Wl,-rpath=. -o Main_lib
Par curiosité, je me suis essayé à installer en l’état ce code comme extension postgreSQL. Voici le Makefile qui permet ceci.
L’usage conjoint des deux librairies fonctionne bien, en revanche postgres ne peut se connecter à aucun serveur X, et j’avoue avoir à ce stade été pressé de poursuivre vers la mise en route de notre algo. Il est juste rassurant de constater que pour compiler et linker correctement, l’ajout dans le Makefile de pkg-config opencv4 --cflags --libs suffit:
MODULES = PgCv
EXTENSION = PgCv
DATA = PgCv--0.0.1.sql
PG_CONFIG = pg_config
PGXS := $(shell $(PG_CONFIG) --pgxs)
PG_LDFLAGS := -lCvWrapper -L. -Wl,-rpath=/home/dac/Bureau/OpenGeo/blog/codes/PgOpenCV `pkg-config opencv4 --cflags --libs`
include $(PGXS)