diff --git a/code/external/deepmatching_1.2.2_c++_linux/Makefile b/code/external/deepmatching_1.2.2_c++_linux/Makefile
index fa644e7..3f31795 100755
--- a/code/external/deepmatching_1.2.2_c++_linux/Makefile
+++ b/code/external/deepmatching_1.2.2_c++_linux/Makefile
@@ -1,46 +1,46 @@
CC=g++
OS_NAME=$(shell uname -s)
ifeq ($(OS_NAME),Linux)
#old# LAPACKLDFLAGS=/usr/lib64/atlas/libsatlas.so # single-threaded blas
- LAPACKLDFLAGS=/usr/lib/atlas-base/libatlas.so /usr/lib/libblas/libblas.so
- #LAPACKLDFLAGS=/usr/lib64/atlas/libtatlas.so # multi-threaded blas
- #BLAS_THREADING=-D MULTITHREADED_BLAS # remove this if wrong
+ #LAPACKLDFLAGS=/usr/lib/atlas-base/libatlas.so /usr/lib/libblas/libblas.so
+ LAPACKLDFLAGS=/usr/lib64/atlas/libtatlas.so.3.10 # multi-threaded blas
+ BLAS_THREADING=-D MULTITHREADED_BLAS # remove this if wrong
endif
ifeq ($(OS_NAME),Darwin) # Mac OS X
LAPACKLDFLAGS=-framework Accelerate # for OS X
endif
LAPACKCFLAGS=-Dinteger=int $(BLAS_THREADING)
#old#
#STATICLAPACKLDFLAGS=-fPIC -Wall -g -fopenmp -static -static-libstdc++ /home/lear/douze/tmp/jpeg-6b/libjpeg.a /usr/lib64/libpng.a /usr/lib64/libz.a /usr/lib64/libblas.a /usr/lib/gcc/x86_64-redhat-linux/4.9.2/libgfortran.a /usr/lib/gcc/x86_64-redhat-linux/4.9.2/libquadmath.a # statically linked version
STATICLAPACKLDFLAGS=-fPIC -Wall -g -fopenmp -static -static-libstdc++ /home/lear/douze/tmp/jpeg-6b/libjpeg.a /usr/lib/x86_64-linux-gnu/libpng.a /usr/lib/x86_64-linux-gnu/libz.a /usr/lib/libblas/libblas.a /usr/lib/gcc/x86_64-linux-gnu/4.8/libgfortran.a /usr/lib/gcc/x86_64-linux-gnu/4.8/libquadmath.a # statically linked version
CFLAGS= -fPIC -Wall -g -std=c++11 $(LAPACKCFLAGS) -fopenmp -DUSE_OPENMP -O3
LDFLAGS=-fPIC -Wall -g -ljpeg -lpng -fopenmp
CPYTHONFLAGS=-I/usr/include/python2.7
SOURCES := $(shell find . -name '*.cpp' ! -name 'deepmatching_matlab.cpp')
OBJ := $(SOURCES:%.cpp=%.o)
HEADERS := $(shell find . -name '*.h')
all: deepmatching
.cpp.o: %.cpp %.h
$(CC) -o $@ $(CFLAGS) -c $+
deepmatching: $(HEADERS) $(OBJ)
$(CC) -o $@ $^ $(LDFLAGS) $(LAPACKLDFLAGS)
deepmatching-static: $(HEADERS) $(OBJ)
$(CC) -o $@ $^ $(STATICLAPACKLDFLAGS)
python: $(HEADERS) $(OBJ)
# swig -python $(CPYTHONFLAGS) deepmatching.i # not necessary, only do if you have swig compiler
g++ $(CFLAGS) -c deepmatching_wrap.c $(CPYTHONFLAGS)
g++ -shared $(LDFLAGS) $(LAPACKLDFLAGS) deepmatching_wrap.o $(OBJ) -o _deepmatching.so $(LIBFLAGS)
clean:
rm -f $(OBJ) deepmatching *~ *.pyc .gdb_history deepmatching_wrap.o _deepmatching.so deepmatching.mex???
diff --git a/code/external/deepmatching_1.2.2_c++_linux/deep_matching.cpp b/code/external/deepmatching_1.2.2_c++_linux/deep_matching.cpp
index 0f9a07c..e9a8297 100755
--- a/code/external/deepmatching_1.2.2_c++_linux/deep_matching.cpp
+++ b/code/external/deepmatching_1.2.2_c++_linux/deep_matching.cpp
@@ -1,943 +1,943 @@
/*
Copyright (C) 2014 Jerome Revaud
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see
*/
#include "deep_matching.h"
#include "std.h"
#include "conv.h"
#include "maxfilter.h"
-
+#include "iostream"
// return size of atomic patches
int get_atomic_patch_size( const dm_params_t* params )
{
int upsize = (1 << params->prior_img_downscale);
return 4*upsize;
}
// crop dimensions to a multiple of patch_size
void get_source_shape( const int width, const int height, const int patch_size, int* res ) {
// crop the reference image to a multiple of patch size
res[0] = patch_size * int(width / patch_size);
res[1] = patch_size * int(height / patch_size);
}
// extract pixel descriptor for both images
void extract_image_desc( image_t* img0, image_t* img1, const dm_params_t* params,
float_layers** desc0, float_layers** desc1 )
{
// slightly reduce img0 size to fit the patch tiling
int patch_size = get_atomic_patch_size( params );
int size[2]; // = {width, height}
get_source_shape( img0->width, img0->height, patch_size, size );
image_crop(img0, size[0], size[1]);
// extract gradient-based information
*desc0 = extract_desc( img0, ¶ms->desc_params, params->n_thread );
*desc1 = extract_desc( img1, ¶ms->desc_params, params->n_thread );
}
void avgpool2( float_layers* hog, const dm_params_t* params )
{
int niter = params->prior_img_downscale;
while(niter--) {
float_layers res = empty_layers(float,hog->tx/2,hog->ty/2,hog->tz);
_avgpool2(hog,&res,params->n_thread);
// replace hog by res
free(hog->pixels);
*hog = res;
}
}
/* compute the grid of parent cell position, and their connection to children cells
cells can be half-overlapping if =1
forces the grid spacing if >0
*/
void prepare_big_cells( const int imshape[2], int cell_size, int overlap, int child_overlap,
int_cube* child_grid, float_image* child_norms, int dense_step,
int_cube* grid, int_cube* children, float_image* norms )
{
int offset, step, gtx, gty;
if( dense_step ) {
step = dense_step;
offset = 0;
// we do not care if the patches are overlapping outside the image
#define grid_size(imsize) (1+imsize/step)
gtx = grid_size(imshape[0]);
gty = grid_size(imshape[1]);
#undef grid_size
} else {
// we want patches fully included in the image
offset = cell_size/2;
step = cell_size/(overlap+1);
#define grid_size(imsize) (1+MAX(0,imsize-2*offset)/step)
gtx = grid_size(imshape[0]);
gty = grid_size(imshape[1]);
#undef grid_size
}
assert(!grid->pixels);
*grid = empty_cube(int,gtx,gty,2);
assert(0<=overlap && overlap<=1);
int nc = pow2(2+child_overlap); // number of children per cell
if(child_grid) {
assert(!norms->pixels);
*norms = image_like(float,grid);
assert(!children->pixels);
*children = empty_cube(int,gtx,gty,nc);
}
_prepare_big_cells( cell_size, offset, step, child_grid, child_norms, grid, children, norms );
}
void sample_patches( float_layers* hog, int_cube* pos, int patch_size, int f, float norm, int n_thread,
float_image* patches, float_array* norms )
{
assert(norm>0);
const int npos = pos->tx*pos->ty;
int_image new_pos = empty_image(int,2,npos);
for(int i=0; i<2*npos; i++)
new_pos.pixels[i] = (pos->pixels[i]-patch_size/2)/f;
patch_size /= f;
const int nh = get_patch_desc_dim(hog,patch_size);
assert(!patches->pixels);
*patches = empty_image(float,nh,npos);
assert(norms->tx==npos);
_sample_patches( hog, NULL, &new_pos, patch_size, norm, patches, norms, n_thread );
free(new_pos.pixels);
}
const float trans_inv = 0.9f;
void convolve_atomic_patches( float_layers* source, float_layers* target,
const dm_params_t* params, res_scale* first_level )
{
const int extend = 1; // slightly spatially extend response maps
const float norm = 1; // renorm patches
const int f = first_level->f; // scale factor w.r.t. original image
const int psize = first_level->patch_size; // current patch size
// first, sample patches
float_image patches = {0};
assert(!first_level->norms.pixels);
first_level->norms = image_like(float, &first_level->grid);
float_array norms_arr = {first_level->norms.pixels, (int)IMG_SIZE(&first_level->norms)};
sample_patches( source, &first_level->grid, psize, f, norm, params->n_thread, &patches, &norms_arr );
//hash_image(&patches)
// rectify the norm to a boolean (0 or 1) (useless ?)
first_level->assign = empty_array(int,norms_arr.tx);
int n=0, tx = patches.tx;
for(int i=0; i0;
// eliminate zero-norm patches
if( norms_arr.pixels[i] ) {
if( n < i ) // copy
memcpy( patches.pixels + n*tx, patches.pixels + i*tx, tx*sizeof(float));
first_level->assign.pixels[i] = n++;
} else
first_level->assign.pixels[i] = -1;
// convolution is not fully invariant to the image border:
// blank cells outside the image are a bit disadvantageous
if( norms_arr.pixels[i] == 0 )
norms_arr.pixels[i] = 1-trans_inv;
}
patches.ty = n; // update new number of valid patches
//hash_image(&first_level->norms)
//hash_image(&patches)
// compute the first level convolutions
fastconv( &patches, target, psize/f, params->ngh_rad/f, extend, norm, params->n_thread, first_level );
free(patches.pixels);
}
int_image* maxpool3_and_subsample2( float_layers* hog, int true_shape[2], int_image* offsets, float_layers* res, int nt )
{
assert(!res->pixels);
if ( offsets->pixels == NULL )
assert( hog->tx == true_shape[0] && hog->ty == true_shape[1] );
// set downsampled size
true_shape[0] = (true_shape[0]+1)/2;
true_shape[1] = (true_shape[1]+1)/2;
assert( true_shape[0]>0 && true_shape[1]>0 );
if ( offsets->pixels == NULL ) {
// joint max-pooling and subsampling
*res = empty_layers(float, true_shape[0], true_shape[1], hog->tz);
_max_filter_3_and_subsample_layers( hog, res, nt );
return NULL;
} else {
// with offsets
float_layers maxpooled_hog = layers_like(float, hog);
_max_filter_3_layers( hog, &maxpooled_hog, nt );
//CHECK_MAPS(&maxpooled_hog);
// slightly bigger, so that mininum size always >= 2
int width = (hog->tx+2)/2;
int height = (hog->ty+2)/2;
*res = empty_layers(float, width, height, hog->tz);
_subsample2_offset( &maxpooled_hog, offsets, res, nt );
free(maxpooled_hog.pixels);
// compute new offsets
int_image* res_offsets = NEW(int_image);
*res_offsets = image_like(int, offsets);
for(long i=0; ipixels[i] = (int)floor( offsets->pixels[i]/2.f );
return res_offsets;
}
}
#define CHECK_MAPS(rmaps) assert(min_array_f((rmaps)->pixels,LAYERS_SIZE(rmaps))>=0 && \
max_array_f((rmaps)->pixels,LAYERS_SIZE(rmaps))<=1.001)
/* aggregate response maps of children patches to form response maps of parent patches */
int sparse_conv( int_cube* children, int_array* children_assign, float_image* child_norms,
int true_patch_size, float_layers* map, int_image* offsets, int nt,
res_scale* res )
{
float_layers ext_map;
if( MIN(map->tx,map->ty) < 5 ) {
ext_map = zeros_layers(float,MAX(5,map->tx),MAX(5,map->ty),map->tz);
for(int l=0; ltz; l++)
for(int j=0; jty; j++)
for(int i=0; itx; i++)
ext_map.pixels[(l*ext_map.ty + j)*ext_map.tx + i] = map->pixels[(l*map->ty + j)*map->tx + i];
map = &ext_map;
res->true_shape[0] = ext_map.tx;
res->true_shape[1] = ext_map.ty;
}
int_image _children = reshape_z_xy(int, &res->children);
if( offsets )
res->offsets = empty_image(int, 2, _children.ty);
assert(!res->res_map.pixels);
res->res_map = empty_layers(float, map->tx, map->ty, _children.ty);
int gap = true_patch_size / 4;
assert(gap>0);
float_array _norms = reshape_xy(float, &res->norms);
float_array _child_norms = reshape_xy(float, child_norms);
// allocate useless assign
res->assign = empty_array(int, res->res_map.tz);
for(int i=0; iassign.tx; i++) res->assign.pixels[i] = i;
int_array* _assign = NULL;
int_array* _ch_assign = children_assign->pixels ? children_assign : NULL;
int n = _sparse_conv( &_children, _ch_assign, gap, trans_inv, map, offsets,
&_child_norms, &_norms, _assign, &res->res_map, &res->offsets, nt );
//CHECK_MAPS(res);
if(map==&ext_map) free(ext_map.pixels);
return n;
}
res_scale new_pyramid_level(int f, int psize)
{
res_scale res = {0}; // initialize everything to 0/NULL
res.f = f; // subsampling factor with respect to original image size
res.patch_size = psize; // patch size in original image coordinates
return res;
}
// Compute the multi-scale pyramid response
void compute_matching_pyr( float_layers* source, float_layers* target, const dm_params_t* params,
matching_pyramid_t& res_maps )
{
const int src_shape[2] = {source->tx, source->ty};
int L = 0; // current pyramid level
const int atomic_psize = get_atomic_patch_size( params );
int psize = atomic_psize; // will grow by a factor 2 at each level
int f = psize/4; // initial scaling factor
// subsample if needed
avgpool2( source, params );
avgpool2( target, params );
//hash_layers(source)
//hash_layers(target)
res_maps.clear();
res_maps.push_back(new_pyramid_level(f,psize));
res_scale *child, *last = &res_maps[res_maps.size()-1];
// compute the initial patches in source image
if( params->verbose ) std_printf("layer %d, patch_size = %dx%d\n", L, psize, psize);
prepare_big_cells( src_shape, psize, params->overlapgrid, NULL, NULL );
//hash_cube(&last->grid)
//hash_layers(source)
convolve_atomic_patches( source, target, params, last );
//hash_layers(&last->res_map)
if( params->verbose )
std_printf("remaining %ld big cells (actually, %d are unique)\n", IMG_SIZE(&last->grid), last->res_map.tz);
// non-linear correction
if( params->nlpow>0 )
fastipow( &last->res_map, params->nlpow, params->n_thread );
//hash_layers(&last->res_map)
const int dense_step = params->subsample_ref ? 0 : psize/(1+(params->overlap<1));
// aggregate patches for all subsequent levels
while( 2*psize <= MIN(params->max_psize, MAX(src_shape[0], src_shape[1])) ) {
L++;
f *= 2;
psize *= 2;
res_maps.push_back(new_pyramid_level(f,psize));
child = &res_maps[res_maps.size()-2]; // previous level
last = &res_maps[res_maps.size()-1]; // current level
if( params->verbose ) std_printf("layer %d, patch_size = %dx%d\n", L, psize, psize);
// max pooling + subsampling
//CHECK_MAPS(&child->res_map);
last->true_shape[0] = child->true_shape[0]; // will be modified in subsampled2()
last->true_shape[1] = child->true_shape[1];
float_layers subs_res_map = {0};
int_image* offsets = maxpool3_and_subsample2( &child->res_map, last->true_shape, &child->offsets,
&subs_res_map, params->n_thread );
//CHECK_MAPS(&subs_res_map);
// build the set of patches at this scale
prepare_big_cells( src_shape, psize, params->overlapoverlapgrid, &child->norms, dense_step, &last->grid, &last->children, &last->norms );
//DA(last->true_shape,2)
//hash_cube(&last->grid)
//hash_image(&last->norms)
//hash_cube(&last->children)
// aggregate children response maps to form parent response maps
sparse_conv( &last->children, &child->assign, &child->norms, psize/f, &subs_res_map, offsets,
params->n_thread, last );
free(subs_res_map.pixels);
free_image(offsets);
//CHECK_MAPS(&last->res_map);
if( params->verbose )
std_printf("remaining %ld big cells (actually, %d are unique)\n", IMG_SIZE(&last->grid), last->res_map.tz);
// non-linear correction
if( params->nlpow>0 )
fastipow(&last->res_map, params->nlpow, params->n_thread );
//hash_layers(&last->res_map)
}
}
void free_matching_pyramid( matching_pyramid_t& res_maps ) {
unsigned int i;
for(i=0; i0); // descending order
}
#else
static int arg_sort_maxima(const void* a, const void* b, void* arr) {
float diff = ((float*)arr)[5*(*(int*)a)+4] - ((float*)arr)[5*(*(int*)b)+4];
return (diff<0) - (diff>0); // descending order
}
#endif
void reorder_rows( int_image* img, int_array* order )
{
assert(order->tx==img->ty);
const int tx = img->tx;
int_image res = image_like(int, img);
for(int i=0; itx; i++)
memcpy(res.pixels + i*tx, img->pixels+order->pixels[i]*tx, tx*sizeof(int));
free(img->pixels);
*img = res;
}
// return points corresponding to patch matches
int_image* find_optimal_matchings( matching_pyramid_t& mp, const dm_params_t* params )
{
const int nobordure = 0;
int_image* maxima = NEW(int_image);
int_array order = {0};
if( params->maxima_mode ) { // normal process: maxima detection
float th=0;
int check_parents=false, check_children=false;
float_array sc_maxima = empty_array(float,int(mp.size()));
for(unsigned int i=0; imin_level, params->nlpow,
check_parents, check_children, nobordure, maxima, params->n_thread );
free(sc_maxima.pixels);
order = empty_array(int,maxima->ty);
for(int i=0; ity; i++) order.pixels[i] = maxima->ty-1-i; // last first
} else { // we just analyse all cells at the top level
const float_layers* rmap = &mp[mp.size()-1].res_map;
const int tx = rmap->tx, txy=tx*rmap->ty;
*maxima = empty_image(int, 5, (int)LAYERS_SIZE(rmap));
for(int i=0; ity; i++) {
int* row = maxima->pixels + 5*i;
row[0] = mp.size()-1; // pyramid level
row[1] = i/txy; // layer number
row[2] = i%tx; // x position
row[3] = (i%txy)/tx; // y position
((float*)row)[4] = rmap->pixels[i];
}
//hash_image(maxima)
order = empty_array(int,maxima->ty);
for(int i=0; ity; i++) order.pixels[i] = i;
#ifdef __APPLE__
qsort_r(order.pixels, maxima->ty, sizeof(int), maxima->pixels, arg_sort_maxima);
#else
qsort_r(order.pixels, maxima->ty, sizeof(int), arg_sort_maxima, maxima->pixels);
#endif
}
if( params->verbose>0 )
std_printf("found %d local matches\n",maxima->ty);
// reorder maxima
reorder_rows( maxima, &order );
free(order.pixels);
return maxima;
}
static inline float ptdot( const float* m, float x, float y ) {
return x*m[0] + y*m[1] + m[2];
}
void apply_rot( float_cube* corres, float rot[6] ) {
assert( corres->tz == 6 );
const int nb = IMG_SIZE(corres);
float* p = corres->pixels;
for(int i=0; itx;
const int n_maxima = maxima->ty;
float_cube corres0 = zeros_cube(float, (src_shape[0]+step-1)/step, (src_shape[1]+step-1)/step,6);
float_cube corres1 = zeros_cube(float, (target_shape[0]+step-1)/step, (target_shape[1]+step-1)/step,6);
int i;
// allocate temporary optimization maps
for(i=0; ilow_mem && size > 1000003 ) size = 1000003; // big prime
assert( size <= 2147483647 || !"try using -mem parameter");
scales[i].passed = zeros_array(float, (int)size);
}
#if defined(USE_OPENMP)
#pragma omp parallel for schedule(static,1) num_threads(params->n_thread)
#endif
for(i=0; iverbose && i%100==0) std_printf("\rgathering correspondences %d%%...",100*i/n_maxima);
int* m = maxima->pixels + tx*i;
int level = m[0], num_map = m[1];
int x = m[2], y = m[3];
assert(levelscoring_mode ) // new mode
_argmax_correspondences( scales.data(), level, num_map, x, y, ((float*)m)[4],
&corres0, step, &corres1, step, i );
else // old iccv mode
_argmax_correspondences_v1( scales.data(), level, num_map, x, y, m[0]*((float*)m)[4],
&corres0, step, &corres1, step, i );
}
// free optimization maps
for(i=0; iverbose) std_printf("\n");
if( params->rot45 ) { // rectify correspondences
assert( corres_out );
apply_rot( &corres0, corres_out->rot );
apply_rot( &corres1, corres_out->rot );
}
// keep only reciprocal matches
int nres;
float* corres = _intersect_corres( &corres0, &corres1, &nres );
float_image* res = NEW(float_image);
*res = (float_image){corres, 6, nres};
if( corres_out == NULL ) {
free(corres0.pixels);
free(corres1.pixels);
}
else { // save unfiltered correspondences
corres_out->corres0 = corres0;
corres_out->corres1 = corres1;
}
return res;
}
void eye_rot3x3( float rot[6] ) {
memset( rot, 0, 6*sizeof(float));
rot[0] = rot[4] = 1;
}
inline float bilinear_interp(const float* img, const int tx, const int ty, float x, float y ) {
if( x < 0 || x+1.001 >= tx ) return 0; // outside
if( y < 0 || y+1.001 >= ty ) return 0; // outside
int ix = int(x);
int iy = int(y);
img += ix + iy*tx; // move pointer
float rx = x - ix;
float ry = y - iy;
return (1-ry)*((1-rx)*img[0] + rx*img[1]) +
ry *((1-rx)*img[tx]+ rx*img[tx+1]);
}
void scale_rot3x3( float rot[6], float sc ) {
for(int i=0; i<6; i++)
rot[i] *= sc;
}
void inv_rot3x3( float rot[6], float res[6] ) {
assert( fabs((rot[0]*rot[4] - rot[1]*rot[3]) - 1) < 1e-6 );
// because rot is unitary, invert == transpose
res[0] = rot[0];
res[1] = rot[3];
res[3] = rot[1];
res[4] = rot[4];
res[2] = -rot[2]*rot[0] - rot[5]*rot[3];
res[5] = -rot[2]*rot[1] - rot[5]*rot[4];
}
// rotate a descriptor HOG image by a given angle
float_layers* rotate45( float_layers* hog, const dm_params_t* params, full_corres_t* corres_out ) {
assert( corres_out ); // we need it to write rot !
const int patch_size = get_atomic_patch_size( params );
const int n_rot45 = params->rot45;
if( (n_rot45 % 8) == 0 ) { // nothing to do
eye_rot3x3( corres_out->rot );
return hog;
}
const int tx = hog->tx;
const int ty = hog->ty;
// rotation matrix
float angle = n_rot45 * M_PI / 4;
float c = cos(angle), s = sin(angle);
float rot[6] = {c, -s, 0, s, c, 0};
// pt_in_original_image = rot * pt_in_rotated_image
// determine center of rotation before
float cx_before = tx/2.0;
float cy_before = ty/2.0;
// determine center of rotation after
float corners[2][4] = {{0, (float)tx, (float)tx, 0}, {0, 0, (float)ty, (float)ty}};
for(int i=0; i<4; i++) { // rotate corners
float x = corners[0][i], y = corners[1][i];
corners[0][i] = ptdot(rot+0, x, y);
corners[1][i] = ptdot(rot+3, x, y);
}
int rot_size[2] = {int(0.5 + max_array_f(corners[0], 4) - min_array_f(corners[0], 4)),
int(0.5 + max_array_f(corners[1], 4) - min_array_f(corners[1], 4)) };
get_source_shape( rot_size[0], rot_size[1], patch_size, rot_size );
float cx_after = rot_size[0]/2.0;
float cy_after = rot_size[1]/2.0;
// compute the translation
rot[2] = cx_before - ptdot(rot+0, cx_after, cy_after);
rot[5] = cy_before - ptdot(rot+3, cx_after, cy_after);
// create result
assert( hog->tz == 9 );
float_layers* rot_hog = NEW(float_layers);
*rot_hog = empty_layers(float, rot_size[0], rot_size[1], 9);
for(int c=0; ctz; c++) {
const int src_c = (c<8) ? int((c+n_rot45+256)%8) : c; // roll channels except for last one (see hog.h)
const float* f = hog->pixels + src_c * IMG_SIZE(hog);
float* p = rot_hog->pixels + c * IMG_SIZE(rot_hog);
for(int y=0; yrot, rot, 6*sizeof(float) );
return rot_hog;
}
// set default parameters
void set_default_dm_params( dm_params_t* params )
{
// pixel descriptor params
set_default_desc_params( ¶ms->desc_params );
// general parameters
params->prior_img_downscale = 1; // resolution R = 1/2^downscale, default = 1/2
params->rot45 = 0; // don't rotate the first image
params->overlap = 999; // don't use overlapping patches
params->subsample_ref = false; // don't subsample patches in reference image (=first image)
params->nlpow = 1.4;
params->ngh_rad = 0; // no limit by default
params->maxima_mode = 0; // don't use maxima, just start from all top patches
params->min_level = 2; // useless
params->max_psize = 999; // maximum patch size
params->low_mem = true; // optimize mem but then results are slightly unstable/non-reproducible
params->verbose = 0;
params->scoring_mode = 1; // improved scoring scheme
params->n_thread = 1; // no multithreading by default
}
// main function
float_image* deep_matching( image_t* img0, image_t* img1, const dm_params_t* params, full_corres_t* corres_out )
{
cout << "###########################################################" << endl;
// verify parameters
assert(between(0,params->prior_img_downscale,3));
assert(between(0,params->overlap,999));
assert(between(0,params->subsample_ref,1));
assert(between(0.1,params->nlpow,10));
assert(between(0,params->ngh_rad,1<<16));
assert(between(0,params->maxima_mode,1));
assert(between(0,params->min_level,4));
assert(between(0,params->low_mem,1));
assert(between(0,params->scoring_mode,1));
assert(between(0,params->verbose,10));
assert(between(1,params->n_thread,128));
cout << "###########################################################" << endl;
// extract pixel descriptors
float_layers *source, *target;
extract_image_desc( img0, img1, params, &source, &target );
if( corres_out ) // the first image is rotated
source = rotate45( source, params, corres_out );
int src_shape[2] = {source->tx, source->ty};
assert( LAYERS_SIZE(source) > 0 );
int target_shape[2] = {target->tx, target->ty};
assert( LAYERS_SIZE(target) > 0 );
cout << "###########################################################" << endl;
//hash_layers(source)
//hash_layers(target)
// compute local matchings
matching_pyramid_t matching_pyr;
compute_matching_pyr( source, target, params, matching_pyr );
free_layers(source);
free_layers(target);
cout << "###########################################################" << endl;
//hash_layers(&matching_pyr[matching_pyr.size()-1].res_map);
// find optmial matchings (maxima)
int_image* maxima = find_optimal_matchings(matching_pyr, params);
cout << "###########################################################" << endl;
//hash_image(maxima);
// select the best displacements (maxpool merge)
float_image* corres = gather_correspondences( src_shape, target_shape, matching_pyr, maxima, params, corres_out );
cout << "###########################################################" << endl;
//hash_image(corres);
// free everything
free_matching_pyramid(matching_pyr);
free_layers(maxima);
return corres;
}
void swap_first_second_img( float_cube* corres ) {
assert( corres->tz == 6 );
const int nb = IMG_SIZE(corres);
float* p = corres->pixels;
for(int i = 0; i < nb; i++) {
float a = p[0];
float b = p[1];
float c = p[2];
float d = p[3];
*p++ = c;
*p++ = d;
*p++ = a;
*p++ = b;
p += 2;
}
}
void rescale_corres( float_cube* corres, float f0, float f1, int code ) {
assert( corres->tz == 6 );
const int nb = IMG_SIZE(corres);
float* p = corres->pixels;
for(int i = 0; i < nb; i++) {
p[0] *= f0;
p[1] *= f0;
p[2] *= f1;
p[3] *= f1;
p[5] = code;
p += 6;
}
}
// set default parameters
void set_default_scalerot_params( scalerot_params_t* params ) {
params->fast = true;
params->min_sc0 = 0; // scale = 2^(-0/2) = 1
params->max_sc0 = 5; // scale = 2^(-5/2) = 0.176
params->min_sc1 = 0;
params->max_sc1 = 5;
params->min_rot = 0; // rot = 0*45 = 0
params->max_rot = 8; // rot = 8*45 = 360
}
// main function for scale/rotation invariant version
float_image* deep_matching_scale_rot( image_t* img0, image_t* img1, dm_params_t* params,
const scalerot_params_t* sr_params ) {
// verify parameters
assert(sr_params->min_sc0 < sr_params->max_sc0);
assert(sr_params->min_sc1 < sr_params->max_sc1);
assert(between(0, sr_params->min_sc0, 5));
assert(between(0, sr_params->max_sc0, 5));
assert(between(0, sr_params->min_sc1, 5));
assert(between(0, sr_params->max_sc1, 5));
assert(sr_params->min_rot >= 0);
assert(between(1,sr_params->max_rot - sr_params->min_rot, 8));
// init shape
const int psize = get_atomic_patch_size(params);
int imshape0[2];
get_source_shape( img0->width, img0->height, psize, imshape0 );
int imshape1[2] = {img1->width, img1->height};
// check dm params to ensure everything goes fine from now on
#define mean_dim(shape) ((shape[0] + shape[1])/2)
params->max_psize = MIN(mean_dim(imshape0), mean_dim(imshape1));
const int verbose = params->verbose;
params->verbose = MAX(0, verbose - 1); // decrease for inner deepmatchings
// prepare output
const int step0 = psize/2;
const int step1 = psize/2;
float_cube all_corres0 = zeros_cube(float, (imshape0[0]+step0/2-1)/step0, (imshape0[1]+step0/2-1)/step0, 6);
float_cube all_corres1 = zeros_cube(float, (imshape1[0]+step1/2-1)/step1, (imshape1[1]+step1/2-1)/step1, 6);
full_corres_t out;
const int NS = 5;
image_t *scaled_images1[NS] = {NULL};
// loop over all scale*rot combinations
for(int sc0 = sr_params->min_sc0;
sc0 < sr_params->max_sc0;
sc0++) {
const float scale0 = pow(2, -0.5*sc0 ); // scale factor for img0
assert( scale0<=1 && sc0<5 );
image_t* scaled_img0 = ( scale0 >= 1 ) ? img0 :
image_resize_bilinear_scale( img0, scale0 );
for(int sc1 = sr_params->min_sc1;
sc1 < sr_params->max_sc1;
sc1++) {
const float scale1 = pow(2, -0.5*sc1 ); // scale factor for img1
assert( scale1<=1 && sc1<5 );
// optimization, deactivate only if eg. both images are blurry
if( sr_params->fast && !(scale0==1 || scale1==1)) continue;
image_t* scaled_img1 = scaled_images1[sc1 - sr_params->min_sc1];
if( scaled_img1 == NULL ) {
scaled_img1 = ( scale1 >= 1 ) ? img1 :
image_resize_bilinear_scale( img1, scale1 );
// remember result
scaled_images1[sc1 - sr_params->min_sc1] = scaled_img1;
}
for(int rotation = sr_params->min_rot;
rotation < sr_params->max_rot;
rotation++) {
assert( rotation >= 0 );
const int rot_scale_code = 8*(sc1*5+sc0) + (rotation%8); // cannot be negative, because of bin count
if( verbose )
std_printf( "processing scale = (x%g, x%g) + rotation = %d deg (code %d)...\n",
scale0, scale1, 45*rotation, rot_scale_code);
float rot0[6], rot1[6];
// compute correspondences with rotated+scaled image
#define max_dim(img) MAX(img->width, img->height)
if( max_dim(scaled_img0) >= max_dim(scaled_img1) ) { // first image is always the largest
params->rot45 = rotation;
float_image* corres = deep_matching(scaled_img0, scaled_img1, params, &out );
free_image( corres ); // we don't care
inv_rot3x3(out.rot, rot0);
eye_rot3x3(rot1);
} else { // scaled_img1 is larger
params->rot45 = -rotation;
float_image* corres = deep_matching(scaled_img1, scaled_img0, params, &out );
free_image( corres ); // we don't care
// swap first and second image coordinates
memswap( &out.corres0, &out.corres1, sizeof(float_cube) );
swap_first_second_img( &out.corres0 );
swap_first_second_img( &out.corres1 );
inv_rot3x3(out.rot, rot1);
eye_rot3x3(rot0);
}
// change scale of correspondences
rescale_corres( &out.corres0, 1/scale0, 1/scale1, rot_scale_code );
rescale_corres( &out.corres1, 1/scale0, 1/scale1, rot_scale_code );
scale_rot3x3(rot0, scale0);
scale_rot3x3(rot1, scale1);
// merge correspondences in the reference frame
merge_corres( rot0, rot1,
psize, psize, &out.corres0, &out.corres1, 2,
step0, step1, &all_corres0, &all_corres1 ); // finer grid for merge
free(out.corres0.pixels);
free(out.corres1.pixels);
}
}
// free memory
if( img0 != scaled_img0 )
image_delete( scaled_img0 );
}
// final intersection
int nres;
float* corres = _intersect_corres( &all_corres0, &all_corres1, &nres );
float_image* res = NEW(float_image);
*res = (float_image){corres, 6, nres};
// free memory
for(int i=0; i
-*/
-#include "deep_matching.h"
-#include "std.h"
-#include "conv.h"
-#include "maxfilter.h"
-
-
-
-// return size of atomic patches
-int get_atomic_patch_size( const dm_params_t* params )
-{
- int upsize = (1 << params->prior_img_downscale);
- return 4*upsize;
-}
-
-// crop dimensions to a multiple of patch_size
-void get_source_shape( const int width, const int height, const int patch_size, int* res ) {
- // crop the reference image to a multiple of patch size
- res[0] = patch_size * int(width / patch_size);
- res[1] = patch_size * int(height / patch_size);
-}
-
-// extract pixel descriptor for both images
-void extract_image_desc( image_t* img0, image_t* img1, const dm_params_t* params,
- float_layers** desc0, float_layers** desc1 )
-{
- // slightly reduce img0 size to fit the patch tiling
- int patch_size = get_atomic_patch_size( params );
-
- int size[2]; // = {width, height}
- get_source_shape( img0->width, img0->height, patch_size, size );
- image_crop(img0, size[0], size[1]);
-
- // extract gradient-based information
- *desc0 = extract_desc( img0, ¶ms->desc_params, params->n_thread );
- *desc1 = extract_desc( img1, ¶ms->desc_params, params->n_thread );
-}
-
-
-void avgpool2( float_layers* hog, const dm_params_t* params )
-{
- int niter = params->prior_img_downscale;
- while(niter--) {
- float_layers res = empty_layers(float,hog->tx/2,hog->ty/2,hog->tz);
- _avgpool2(hog,&res,params->n_thread);
-
- // replace hog by res
- free(hog->pixels);
- *hog = res;
- }
-}
-
-
-/* compute the grid of parent cell position, and their connection to children cells
- cells can be half-overlapping if =1
- forces the grid spacing if >0
-*/
-void prepare_big_cells( const int imshape[2], int cell_size, int overlap, int child_overlap,
- int_cube* child_grid, float_image* child_norms, int dense_step,
- int_cube* grid, int_cube* children, float_image* norms )
-{
- int offset, step, gtx, gty;
- if( dense_step ) {
- step = dense_step;
- offset = 0;
- // we do not care if the patches are overlapping outside the image
- #define grid_size(imsize) (1+imsize/step)
- gtx = grid_size(imshape[0]);
- gty = grid_size(imshape[1]);
- #undef grid_size
- } else {
- // we want patches fully included in the image
- offset = cell_size/2;
- step = cell_size/(overlap+1);
- #define grid_size(imsize) (1+MAX(0,imsize-2*offset)/step)
- gtx = grid_size(imshape[0]);
- gty = grid_size(imshape[1]);
- #undef grid_size
- }
-
- assert(!grid->pixels);
- *grid = empty_cube(int,gtx,gty,2);
-
- assert(0<=overlap && overlap<=1);
- int nc = pow2(2+child_overlap); // number of children per cell
- if(child_grid) {
- assert(!norms->pixels);
- *norms = image_like(float,grid);
- assert(!children->pixels);
- *children = empty_cube(int,gtx,gty,nc);
- }
-
- _prepare_big_cells( cell_size, offset, step, child_grid, child_norms, grid, children, norms );
-}
-
-
-void sample_patches( float_layers* hog, int_cube* pos, int patch_size, int f, float norm, int n_thread,
- float_image* patches, float_array* norms )
-{
- assert(norm>0);
- const int npos = pos->tx*pos->ty;
- int_image new_pos = empty_image(int,2,npos);
- for(int i=0; i<2*npos; i++)
- new_pos.pixels[i] = (pos->pixels[i]-patch_size/2)/f;
-
- patch_size /= f;
- const int nh = get_patch_desc_dim(hog,patch_size);
-
- assert(!patches->pixels);
- *patches = empty_image(float,nh,npos);
- assert(norms->tx==npos);
-
- _sample_patches( hog, NULL, &new_pos, patch_size, norm, patches, norms, n_thread );
-
- free(new_pos.pixels);
-}
-
-
-const float trans_inv = 0.9f;
-
-void convolve_atomic_patches( float_layers* source, float_layers* target,
- const dm_params_t* params, res_scale* first_level )
-{
- const int extend = 1; // slightly spatially extend response maps
- const float norm = 1; // renorm patches
-
- const int f = first_level->f; // scale factor w.r.t. original image
- const int psize = first_level->patch_size; // current patch size
-
- // first, sample patches
- float_image patches = {0};
- assert(!first_level->norms.pixels);
- first_level->norms = image_like(float, &first_level->grid);
- float_array norms_arr = {first_level->norms.pixels, (int)IMG_SIZE(&first_level->norms)};
- sample_patches( source, &first_level->grid, psize, f, norm, params->n_thread, &patches, &norms_arr );
- //hash_image(&patches)
-
- // rectify the norm to a boolean (0 or 1) (useless ?)
- first_level->assign = empty_array(int,norms_arr.tx);
- int n=0, tx = patches.tx;
- for(int i=0; i0;
-
- // eliminate zero-norm patches
- if( norms_arr.pixels[i] ) {
- if( n < i ) // copy
- memcpy( patches.pixels + n*tx, patches.pixels + i*tx, tx*sizeof(float));
- first_level->assign.pixels[i] = n++;
- } else
- first_level->assign.pixels[i] = -1;
-
- // convolution is not fully invariant to the image border:
- // blank cells outside the image are a bit disadvantageous
- if( norms_arr.pixels[i] == 0 )
- norms_arr.pixels[i] = 1-trans_inv;
- }
- patches.ty = n; // update new number of valid patches
-
- //hash_image(&first_level->norms)
- //hash_image(&patches)
-
- // compute the first level convolutions
- fastconv( &patches, target, psize/f, params->ngh_rad/f, extend, norm, params->n_thread, first_level );
-
- free(patches.pixels);
-}
-
-int_image* maxpool3_and_subsample2( float_layers* hog, int true_shape[2], int_image* offsets, float_layers* res, int nt )
-{
- assert(!res->pixels);
- if ( offsets->pixels == NULL )
- assert( hog->tx == true_shape[0] && hog->ty == true_shape[1] );
-
- // set downsampled size
- true_shape[0] = (true_shape[0]+1)/2;
- true_shape[1] = (true_shape[1]+1)/2;
- assert( true_shape[0]>0 && true_shape[1]>0 );
-
- if ( offsets->pixels == NULL ) {
- // joint max-pooling and subsampling
- *res = empty_layers(float, true_shape[0], true_shape[1], hog->tz);
- _max_filter_3_and_subsample_layers( hog, res, nt );
- return NULL;
-
- } else {
- // with offsets
- float_layers maxpooled_hog = layers_like(float, hog);
- _max_filter_3_layers( hog, &maxpooled_hog, nt );
- //CHECK_MAPS(&maxpooled_hog);
-
- // slightly bigger, so that mininum size always >= 2
- int width = (hog->tx+2)/2;
- int height = (hog->ty+2)/2;
- *res = empty_layers(float, width, height, hog->tz);
- _subsample2_offset( &maxpooled_hog, offsets, res, nt );
- free(maxpooled_hog.pixels);
-
- // compute new offsets
- int_image* res_offsets = NEW(int_image);
- *res_offsets = image_like(int, offsets);
- for(long i=0; ipixels[i] = (int)floor( offsets->pixels[i]/2.f );
- return res_offsets;
- }
-}
-
-#define CHECK_MAPS(rmaps) assert(min_array_f((rmaps)->pixels,LAYERS_SIZE(rmaps))>=0 && \
- max_array_f((rmaps)->pixels,LAYERS_SIZE(rmaps))<=1.001)
-
-/* aggregate response maps of children patches to form response maps of parent patches */
-int sparse_conv( int_cube* children, int_array* children_assign, float_image* child_norms,
- int true_patch_size, float_layers* map, int_image* offsets, int nt,
- res_scale* res )
-{
- float_layers ext_map;
- if( MIN(map->tx,map->ty) < 5 ) {
- ext_map = zeros_layers(float,MAX(5,map->tx),MAX(5,map->ty),map->tz);
- for(int l=0; ltz; l++)
- for(int j=0; jty; j++)
- for(int i=0; itx; i++)
- ext_map.pixels[(l*ext_map.ty + j)*ext_map.tx + i] = map->pixels[(l*map->ty + j)*map->tx + i];
- map = &ext_map;
- res->true_shape[0] = ext_map.tx;
- res->true_shape[1] = ext_map.ty;
- }
-
- int_image _children = reshape_z_xy(int, &res->children);
-
- if( offsets )
- res->offsets = empty_image(int, 2, _children.ty);
-
- assert(!res->res_map.pixels);
- res->res_map = empty_layers(float, map->tx, map->ty, _children.ty);
- int gap = true_patch_size / 4;
- assert(gap>0);
- float_array _norms = reshape_xy(float, &res->norms);
- float_array _child_norms = reshape_xy(float, child_norms);
-
- // allocate useless assign
- res->assign = empty_array(int, res->res_map.tz);
- for(int i=0; iassign.tx; i++) res->assign.pixels[i] = i;
-
- int_array* _assign = NULL;
- int_array* _ch_assign = children_assign->pixels ? children_assign : NULL;
- int n = _sparse_conv( &_children, _ch_assign, gap, trans_inv, map, offsets,
- &_child_norms, &_norms, _assign, &res->res_map, &res->offsets, nt );
- //CHECK_MAPS(res);
-
- if(map==&ext_map) free(ext_map.pixels);
- return n;
-}
-
-res_scale new_pyramid_level(int f, int psize)
-{
- res_scale res = {0}; // initialize everything to 0/NULL
- res.f = f; // subsampling factor with respect to original image size
- res.patch_size = psize; // patch size in original image coordinates
- return res;
-}
-
-// Compute the multi-scale pyramid response
-void compute_matching_pyr( float_layers* source, float_layers* target, const dm_params_t* params,
- matching_pyramid_t& res_maps )
-{
- const int src_shape[2] = {source->tx, source->ty};
- int L = 0; // current pyramid level
- const int atomic_psize = get_atomic_patch_size( params );
- int psize = atomic_psize; // will grow by a factor 2 at each level
- int f = psize/4; // initial scaling factor
-
- // subsample if needed
- avgpool2( source, params );
- avgpool2( target, params );
-
- //hash_layers(source)
- //hash_layers(target)
-
- res_maps.clear();
- res_maps.push_back(new_pyramid_level(f,psize));
- res_scale *child, *last = &res_maps[res_maps.size()-1];
-
- // compute the initial patches in source image
- if( params->verbose ) std_printf("layer %d, patch_size = %dx%d\n", L, psize, psize);
- prepare_big_cells( src_shape, psize, params->overlapgrid, NULL, NULL );
- //hash_cube(&last->grid)
-
- //hash_layers(source)
- convolve_atomic_patches( source, target, params, last );
- //hash_layers(&last->res_map)
- if( params->verbose )
- std_printf("remaining %ld big cells (actually, %d are unique)\n", IMG_SIZE(&last->grid), last->res_map.tz);
-
- // non-linear correction
- if( params->nlpow>0 )
- fastipow( &last->res_map, params->nlpow, params->n_thread );
-
- //hash_layers(&last->res_map)
-
- const int dense_step = params->subsample_ref ? 0 : psize/(1+(params->overlap<1));
-
- // aggregate patches for all subsequent levels
- while( 2*psize <= MIN(params->max_psize, MAX(src_shape[0], src_shape[1])) ) {
- L++;
- f *= 2;
- psize *= 2;
- res_maps.push_back(new_pyramid_level(f,psize));
- child = &res_maps[res_maps.size()-2]; // previous level
- last = &res_maps[res_maps.size()-1]; // current level
- if( params->verbose ) std_printf("layer %d, patch_size = %dx%d\n", L, psize, psize);
-
- // max pooling + subsampling
- //CHECK_MAPS(&child->res_map);
- last->true_shape[0] = child->true_shape[0]; // will be modified in subsampled2()
- last->true_shape[1] = child->true_shape[1];
- float_layers subs_res_map = {0};
- int_image* offsets = maxpool3_and_subsample2( &child->res_map, last->true_shape, &child->offsets,
- &subs_res_map, params->n_thread );
- //CHECK_MAPS(&subs_res_map);
-
- // build the set of patches at this scale
- prepare_big_cells( src_shape, psize, params->overlapoverlapgrid, &child->norms, dense_step, &last->grid, &last->children, &last->norms );
- //DA(last->true_shape,2)
- //hash_cube(&last->grid)
- //hash_image(&last->norms)
- //hash_cube(&last->children)
-
- // aggregate children response maps to form parent response maps
- sparse_conv( &last->children, &child->assign, &child->norms, psize/f, &subs_res_map, offsets,
- params->n_thread, last );
- free(subs_res_map.pixels);
- free_image(offsets);
- //CHECK_MAPS(&last->res_map);
- if( params->verbose )
- std_printf("remaining %ld big cells (actually, %d are unique)\n", IMG_SIZE(&last->grid), last->res_map.tz);
-
- // non-linear correction
- if( params->nlpow>0 )
- fastipow(&last->res_map, params->nlpow, params->n_thread );
- //hash_layers(&last->res_map)
- }
-}
-
-
-void free_matching_pyramid( matching_pyramid_t& res_maps ) {
- unsigned int i;
- for(i=0; i0); // descending order
-}
-#else
-static int arg_sort_maxima(const void* a, const void* b, void* arr) {
- float diff = ((float*)arr)[5*(*(int*)a)+4] - ((float*)arr)[5*(*(int*)b)+4];
- return (diff<0) - (diff>0); // descending order
-}
-#endif
-
-void reorder_rows( int_image* img, int_array* order )
-{
- assert(order->tx==img->ty);
- const int tx = img->tx;
- int_image res = image_like(int, img);
-
- for(int i=0; itx; i++)
- memcpy(res.pixels + i*tx, img->pixels+order->pixels[i]*tx, tx*sizeof(int));
-
- free(img->pixels);
- *img = res;
-}
-
-// return points corresponding to patch matches
-int_image* find_optimal_matchings( matching_pyramid_t& mp, const dm_params_t* params )
-{
- const int nobordure = 0;
- int_image* maxima = NEW(int_image);
- int_array order = {0};
-
- if( params->maxima_mode ) { // normal process: maxima detection
-
- float th=0;
- int check_parents=false, check_children=false;
-
- float_array sc_maxima = empty_array(float,int(mp.size()));
- for(unsigned int i=0; imin_level, params->nlpow,
- check_parents, check_children, nobordure, maxima, params->n_thread );
- free(sc_maxima.pixels);
-
- order = empty_array(int,maxima->ty);
- for(int i=0; ity; i++) order.pixels[i] = maxima->ty-1-i; // last first
-
- } else { // we just analyse all cells at the top level
- const float_layers* rmap = &mp[mp.size()-1].res_map;
- const int tx = rmap->tx, txy=tx*rmap->ty;
- *maxima = empty_image(int, 5, (int)LAYERS_SIZE(rmap));
-
- for(int i=0; ity; i++) {
- int* row = maxima->pixels + 5*i;
- row[0] = mp.size()-1; // pyramid level
- row[1] = i/txy; // layer number
- row[2] = i%tx; // x position
- row[3] = (i%txy)/tx; // y position
- ((float*)row)[4] = rmap->pixels[i];
- }
- //hash_image(maxima)
-
- order = empty_array(int,maxima->ty);
- for(int i=0; ity; i++) order.pixels[i] = i;
- #ifdef __APPLE__
- qsort_r(order.pixels, maxima->ty, sizeof(int), maxima->pixels, arg_sort_maxima);
- #else
- qsort_r(order.pixels, maxima->ty, sizeof(int), arg_sort_maxima, maxima->pixels);
- #endif
- }
-
- if( params->verbose>0 )
- std_printf("found %d local matches\n",maxima->ty);
-
- // reorder maxima
- reorder_rows( maxima, &order );
- free(order.pixels);
- return maxima;
-}
-
-
-static inline float ptdot( const float* m, float x, float y ) {
- return x*m[0] + y*m[1] + m[2];
-}
-
-void apply_rot( float_cube* corres, float rot[6] ) {
- assert( corres->tz == 6 );
- const int nb = IMG_SIZE(corres);
- float* p = corres->pixels;
-
- for(int i=0; itx;
- const int n_maxima = maxima->ty;
-
- float_cube corres0 = zeros_cube(float, (src_shape[0]+step-1)/step, (src_shape[1]+step-1)/step,6);
- float_cube corres1 = zeros_cube(float, (target_shape[0]+step-1)/step, (target_shape[1]+step-1)/step,6);
-
- int i;
- // allocate temporary optimization maps
- for(i=0; ilow_mem && size > 1000003 ) size = 1000003; // big prime
- assert( size <= 2147483647 || !"try using -mem parameter");
- scales[i].passed = zeros_array(float, (int)size);
- }
-
- #if defined(USE_OPENMP)
- #pragma omp parallel for schedule(static,1) num_threads(params->n_thread)
- #endif
- for(i=0; iverbose && i%100==0) std_printf("\rgathering correspondences %d%%...",100*i/n_maxima);
- int* m = maxima->pixels + tx*i;
- int level = m[0], num_map = m[1];
- int x = m[2], y = m[3];
- assert(levelscoring_mode ) // new mode
- _argmax_correspondences( scales.data(), level, num_map, x, y, ((float*)m)[4],
- &corres0, step, &corres1, step, i );
- else // old iccv mode
- _argmax_correspondences_v1( scales.data(), level, num_map, x, y, m[0]*((float*)m)[4],
- &corres0, step, &corres1, step, i );
- }
-
- // free optimization maps
- for(i=0; iverbose) std_printf("\n");
-
- if( params->rot45 ) { // rectify correspondences
- assert( corres_out );
- apply_rot( &corres0, corres_out->rot );
- apply_rot( &corres1, corres_out->rot );
- }
-
- // keep only reciprocal matches
- int nres;
- float* corres = _intersect_corres( &corres0, &corres1, &nres );
- float_image* res = NEW(float_image);
- *res = (float_image){corres, 6, nres};
-
- if( corres_out == NULL ) {
- free(corres0.pixels);
- free(corres1.pixels);
- }
- else { // save unfiltered correspondences
- corres_out->corres0 = corres0;
- corres_out->corres1 = corres1;
- }
- return res;
-}
-
-
-void eye_rot3x3( float rot[6] ) {
- memset( rot, 0, 6*sizeof(float));
- rot[0] = rot[4] = 1;
-}
-
-inline float bilinear_interp(const float* img, const int tx, const int ty, float x, float y ) {
- if( x < 0 || x+1.001 >= tx ) return 0; // outside
- if( y < 0 || y+1.001 >= ty ) return 0; // outside
- int ix = int(x);
- int iy = int(y);
- img += ix + iy*tx; // move pointer
- float rx = x - ix;
- float ry = y - iy;
- return (1-ry)*((1-rx)*img[0] + rx*img[1]) +
- ry *((1-rx)*img[tx]+ rx*img[tx+1]);
-}
-
-void scale_rot3x3( float rot[6], float sc ) {
- for(int i=0; i<6; i++)
- rot[i] *= sc;
-}
-
-void inv_rot3x3( float rot[6], float res[6] ) {
- assert( fabs((rot[0]*rot[4] - rot[1]*rot[3]) - 1) < 1e-6 );
- // because rot is unitary, invert == transpose
- res[0] = rot[0];
- res[1] = rot[3];
- res[3] = rot[1];
- res[4] = rot[4];
- res[2] = -rot[2]*rot[0] - rot[5]*rot[3];
- res[5] = -rot[2]*rot[1] - rot[5]*rot[4];
-}
-
-
-
-// rotate a descriptor HOG image by a given angle
-float_layers* rotate45( float_layers* hog, const dm_params_t* params, full_corres_t* corres_out ) {
- assert( corres_out ); // we need it to write rot !
- const int patch_size = get_atomic_patch_size( params );
- const int n_rot45 = params->rot45;
-
- if( (n_rot45 % 8) == 0 ) { // nothing to do
- eye_rot3x3( corres_out->rot );
- return hog;
- }
- const int tx = hog->tx;
- const int ty = hog->ty;
-
- // rotation matrix
- float angle = n_rot45 * M_PI / 4;
- float c = cos(angle), s = sin(angle);
- float rot[6] = {c, -s, 0, s, c, 0};
- // pt_in_original_image = rot * pt_in_rotated_image
-
- // determine center of rotation before
- float cx_before = tx/2.0;
- float cy_before = ty/2.0;
- // determine center of rotation after
- float corners[2][4] = {{0, (float)tx, (float)tx, 0}, {0, 0, (float)ty, (float)ty}};
- for(int i=0; i<4; i++) { // rotate corners
- float x = corners[0][i], y = corners[1][i];
- corners[0][i] = ptdot(rot+0, x, y);
- corners[1][i] = ptdot(rot+3, x, y);
- }
- int rot_size[2] = {int(0.5 + max_array_f(corners[0], 4) - min_array_f(corners[0], 4)),
- int(0.5 + max_array_f(corners[1], 4) - min_array_f(corners[1], 4)) };
- get_source_shape( rot_size[0], rot_size[1], patch_size, rot_size );
- float cx_after = rot_size[0]/2.0;
- float cy_after = rot_size[1]/2.0;
- // compute the translation
- rot[2] = cx_before - ptdot(rot+0, cx_after, cy_after);
- rot[5] = cy_before - ptdot(rot+3, cx_after, cy_after);
-
- // create result
- assert( hog->tz == 9 );
- float_layers* rot_hog = NEW(float_layers);
- *rot_hog = empty_layers(float, rot_size[0], rot_size[1], 9);
-
- for(int c=0; ctz; c++) {
- const int src_c = (c<8) ? int((c+n_rot45+256)%8) : c; // roll channels except for last one (see hog.h)
- const float* f = hog->pixels + src_c * IMG_SIZE(hog);
- float* p = rot_hog->pixels + c * IMG_SIZE(rot_hog);
-
- for(int y=0; yrot, rot, 6*sizeof(float) );
-
- return rot_hog;
-}
-
-
-// set default parameters
-void set_default_dm_params( dm_params_t* params )
-{
- // pixel descriptor params
- set_default_desc_params( ¶ms->desc_params );
-
- // general parameters
- params->prior_img_downscale = 1; // resolution R = 1/2^downscale, default = 1/2
- params->rot45 = 0; // don't rotate the first image
- params->overlap = 999; // don't use overlapping patches
- params->subsample_ref = false; // don't subsample patches in reference image (=first image)
- params->nlpow = 1.4;
- params->ngh_rad = 0; // no limit by default
- params->maxima_mode = 0; // don't use maxima, just start from all top patches
- params->min_level = 2; // useless
- params->max_psize = 999; // maximum patch size
- params->low_mem = true; // optimize mem but then results are slightly unstable/non-reproducible
- params->verbose = 0;
- params->scoring_mode = 1; // improved scoring scheme
- params->n_thread = 1; // no multithreading by default
-}
-
-
-// main function
-float_image* deep_matching( image_t* img0, image_t* img1, const dm_params_t* params, full_corres_t* corres_out )
-{
- // verify parameters
- assert(between(0,params->prior_img_downscale,3));
- assert(between(0,params->overlap,999));
- assert(between(0,params->subsample_ref,1));
- assert(between(0.1,params->nlpow,10));
- assert(between(0,params->ngh_rad,1<<16));
- assert(between(0,params->maxima_mode,1));
- assert(between(0,params->min_level,4));
- assert(between(0,params->low_mem,1));
- assert(between(0,params->scoring_mode,1));
- assert(between(0,params->verbose,10));
- assert(between(1,params->n_thread,128));
-
- cout << "###########################################################" << endl;
-
- // extract pixel descriptors
- float_layers *source, *target;
- extract_image_desc( img0, img1, params, &source, &target );
- if( corres_out ) // the first image is rotated
- source = rotate45( source, params, corres_out );
- int src_shape[2] = {source->tx, source->ty};
- assert( LAYERS_SIZE(source) > 0 );
- int target_shape[2] = {target->tx, target->ty};
- assert( LAYERS_SIZE(target) > 0 );
- cout << "###########################################################" << endl;
-
- //hash_layers(source)
- //hash_layers(target)
-
- // compute local matchings
- matching_pyramid_t matching_pyr;
- compute_matching_pyr( source, target, params, matching_pyr );
- free_layers(source);
- free_layers(target);
- cout << "###########################################################" << endl;
-
- //hash_layers(&matching_pyr[matching_pyr.size()-1].res_map);
-
- // find optmial matchings (maxima)
- int_image* maxima = find_optimal_matchings(matching_pyr, params);
- cout << "###########################################################" << endl;
-
- //hash_image(maxima);
-
- // select the best displacements (maxpool merge)
- float_image* corres = gather_correspondences( src_shape, target_shape, matching_pyr, maxima, params, corres_out );
- cout << "###########################################################" << endl;
-
- //hash_image(corres);
-
- // free everything
- free_matching_pyramid(matching_pyr);
- free_layers(maxima);
-
- return corres;
-}
-
-
-void swap_first_second_img( float_cube* corres ) {
- assert( corres->tz == 6 );
- const int nb = IMG_SIZE(corres);
- float* p = corres->pixels;
-
- for(int i = 0; i < nb; i++) {
- float a = p[0];
- float b = p[1];
- float c = p[2];
- float d = p[3];
- *p++ = c;
- *p++ = d;
- *p++ = a;
- *p++ = b;
- p += 2;
- }
-}
-
-void rescale_corres( float_cube* corres, float f0, float f1, int code ) {
- assert( corres->tz == 6 );
- const int nb = IMG_SIZE(corres);
- float* p = corres->pixels;
-
- for(int i = 0; i < nb; i++) {
- p[0] *= f0;
- p[1] *= f0;
- p[2] *= f1;
- p[3] *= f1;
- p[5] = code;
- p += 6;
- }
-}
-
-// set default parameters
-void set_default_scalerot_params( scalerot_params_t* params ) {
- params->fast = true;
- params->min_sc0 = 0; // scale = 2^(-0/2) = 1
- params->max_sc0 = 5; // scale = 2^(-5/2) = 0.176
- params->min_sc1 = 0;
- params->max_sc1 = 5;
- params->min_rot = 0; // rot = 0*45 = 0
- params->max_rot = 8; // rot = 8*45 = 360
-}
-
-
-// main function for scale/rotation invariant version
-float_image* deep_matching_scale_rot( image_t* img0, image_t* img1, dm_params_t* params,
- const scalerot_params_t* sr_params ) {
- // verify parameters
- assert(sr_params->min_sc0 < sr_params->max_sc0);
- assert(sr_params->min_sc1 < sr_params->max_sc1);
- assert(between(0, sr_params->min_sc0, 5));
- assert(between(0, sr_params->max_sc0, 5));
- assert(between(0, sr_params->min_sc1, 5));
- assert(between(0, sr_params->max_sc1, 5));
- assert(sr_params->min_rot >= 0);
- assert(between(1,sr_params->max_rot - sr_params->min_rot, 8));
-
- // init shape
- const int psize = get_atomic_patch_size(params);
- int imshape0[2];
- get_source_shape( img0->width, img0->height, psize, imshape0 );
- int imshape1[2] = {img1->width, img1->height};
-
- // check dm params to ensure everything goes fine from now on
- #define mean_dim(shape) ((shape[0] + shape[1])/2)
- params->max_psize = MIN(mean_dim(imshape0), mean_dim(imshape1));
- const int verbose = params->verbose;
- params->verbose = MAX(0, verbose - 1); // decrease for inner deepmatchings
-
- // prepare output
- const int step0 = psize/2;
- const int step1 = psize/2;
- float_cube all_corres0 = zeros_cube(float, (imshape0[0]+step0/2-1)/step0, (imshape0[1]+step0/2-1)/step0, 6);
- float_cube all_corres1 = zeros_cube(float, (imshape1[0]+step1/2-1)/step1, (imshape1[1]+step1/2-1)/step1, 6);
- full_corres_t out;
-
- const int NS = 5;
- image_t *scaled_images1[NS] = {NULL};
-
- // loop over all scale*rot combinations
- for(int sc0 = sr_params->min_sc0;
- sc0 < sr_params->max_sc0;
- sc0++) {
- const float scale0 = pow(2, -0.5*sc0 ); // scale factor for img0
- assert( scale0<=1 && sc0<5 );
- image_t* scaled_img0 = ( scale0 >= 1 ) ? img0 :
- image_resize_bilinear_scale( img0, scale0 );
-
- for(int sc1 = sr_params->min_sc1;
- sc1 < sr_params->max_sc1;
- sc1++) {
- const float scale1 = pow(2, -0.5*sc1 ); // scale factor for img1
- assert( scale1<=1 && sc1<5 );
- // optimization, deactivate only if eg. both images are blurry
- if( sr_params->fast && !(scale0==1 || scale1==1)) continue;
-
- image_t* scaled_img1 = scaled_images1[sc1 - sr_params->min_sc1];
- if( scaled_img1 == NULL ) {
- scaled_img1 = ( scale1 >= 1 ) ? img1 :
- image_resize_bilinear_scale( img1, scale1 );
- // remember result
- scaled_images1[sc1 - sr_params->min_sc1] = scaled_img1;
- }
-
- for(int rotation = sr_params->min_rot;
- rotation < sr_params->max_rot;
- rotation++) {
- assert( rotation >= 0 );
- const int rot_scale_code = 8*(sc1*5+sc0) + (rotation%8); // cannot be negative, because of bin count
-
- if( verbose )
- std_printf( "processing scale = (x%g, x%g) + rotation = %d deg (code %d)...\n",
- scale0, scale1, 45*rotation, rot_scale_code);
-
- float rot0[6], rot1[6];
-
- // compute correspondences with rotated+scaled image
- #define max_dim(img) MAX(img->width, img->height)
- if( max_dim(scaled_img0) >= max_dim(scaled_img1) ) { // first image is always the largest
- params->rot45 = rotation;
-
- float_image* corres = deep_matching(scaled_img0, scaled_img1, params, &out );
- free_image( corres ); // we don't care
-
- inv_rot3x3(out.rot, rot0);
- eye_rot3x3(rot1);
-
- } else { // scaled_img1 is larger
- params->rot45 = -rotation;
-
- float_image* corres = deep_matching(scaled_img1, scaled_img0, params, &out );
- free_image( corres ); // we don't care
-
- // swap first and second image coordinates
- memswap( &out.corres0, &out.corres1, sizeof(float_cube) );
- swap_first_second_img( &out.corres0 );
- swap_first_second_img( &out.corres1 );
-
- inv_rot3x3(out.rot, rot1);
- eye_rot3x3(rot0);
- }
-
- // change scale of correspondences
- rescale_corres( &out.corres0, 1/scale0, 1/scale1, rot_scale_code );
- rescale_corres( &out.corres1, 1/scale0, 1/scale1, rot_scale_code );
- scale_rot3x3(rot0, scale0);
- scale_rot3x3(rot1, scale1);
-
- // merge correspondences in the reference frame
- merge_corres( rot0, rot1,
- psize, psize, &out.corres0, &out.corres1, 2,
- step0, step1, &all_corres0, &all_corres1 ); // finer grid for merge
-
- free(out.corres0.pixels);
- free(out.corres1.pixels);
- }
- }
-
- // free memory
- if( img0 != scaled_img0 )
- image_delete( scaled_img0 );
- }
-
- // final intersection
- int nres;
- float* corres = _intersect_corres( &all_corres0, &all_corres1, &nres );
- float_image* res = NEW(float_image);
- *res = (float_image){corres, 6, nres};
-
- // free memory
- for(int i=0; i
-*/
-#include
-#include
-#include
-#include
-#include
-
-#include
-#include
-#include
-#include
-
-void std_printf(const char* format, ... ) {
- va_list arglist;
- va_start( arglist, format );
- char buffer[1024];
- vsprintf( buffer, format, arglist );
- va_end(arglist);
-
- mexPrintf(buffer);
-}
-
-void err_printf(const char* format, ... ) {
- va_list arglist;
- va_start( arglist, format );
- char buffer[1024];
- vsprintf( buffer, format, arglist );
- va_end(arglist);
-
- mexErrMsgTxt(buffer);
-}
-
-
-#include "image.h"
-#include "deep_matching.h"
-#include "io.h"
-#include "main.h"
-
-
-static inline bool ispowerof2( long n ) {
- return (n & (n-1))==0;
-}
-
-color_image_t *input3darray_to_color_image(const mxArray *p){
- const int *dims = mxGetDimensions(p);
- const int h = dims[0], w = dims[1];
- assert( dims[2]==3 );
- float *in = (float*) mxGetData(p);
- color_image_t *out = color_image_new(w, h);
- for(int c=0 ; c<3 ; c++){
- float *inptr = in + c*w*h;
- float *outptr = out->c1 + c*w*h;
- for( int j=0 ; jty, w = corres->tx;
- float *data = (float*) mxGetData(p);
- for( int j=0 ; jpixels[j*w+i];
- }
- }
-}
-
-void mexFunction( int nl, mxArray *pl[], int nr, const mxArray *pr[] ) {
-
- if( nr==0 ) {
- usage(MATLAB_OPTIONS);
- return;
- }
-
- if ( nl != 1){
- usage(MATLAB_OPTIONS);
- mexErrMsgTxt("error: returns one output");
- }
- if( nr < 2 || nr > 3){
- usage(MATLAB_OPTIONS);
- mexErrMsgTxt("error: takes two to four inputs");
- }
- // The code is originally written for C-order arrays.
- // We thus transpose all arrays in this mex-function which is not efficient...
-
- const int *pDims;
- if(mxGetNumberOfDimensions(pr[0]) != 3) mexErrMsgTxt("input images must have 3 dimensions");
- if(!mxIsClass(pr[0], "single")) mexErrMsgTxt("input images must be single");
- pDims = mxGetDimensions(pr[0]);
- if( pDims[2]!=3 ) mexErrMsgTxt("input images must have 3 channels");
- const int h = pDims[0], w = pDims[1];
- color_image_t *cim1 = input3darray_to_color_image( pr[0] );
-
- if(mxGetNumberOfDimensions(pr[1]) != 3) mexErrMsgTxt("input images must have 3 dimensions");
- if(!mxIsClass(pr[1], "single")) mexErrMsgTxt("input images must be single");
- pDims = mxGetDimensions(pr[1]);
- if( pDims[2]!=3) mexErrMsgTxt("input images must have 3 channels");
- color_image_t *cim2 = input3darray_to_color_image( pr[1] );
-
- // convert images to gray
- image_t *im1=image_gray_from_color(cim1), *im2=image_gray_from_color(cim2);;
- color_image_delete(cim1);
- color_image_delete(cim2);
-
- // set params to default
- dm_params_t params;
- set_default_dm_params(¶ms);
- scalerot_params_t sr_params;
- set_default_scalerot_params(&sr_params);
- bool use_scalerot = false;
- float fx=1, fy=1;
-
- // read options
- if( nr == 3 ){
- char *options = mxArrayToString(pr[2]);
- if( !options ) mexErrMsgTxt("Third parameter must be a string");
- int argc=0;
- const char* argv[256];
- argv[argc] = strtok(options," ");
- while(argv[argc]!=NULL)
- argv[++argc] = strtok(NULL," ");
-
- parse_options(¶ms, &sr_params, &use_scalerot, &fx, &fy, argc, argv, MATLAB_OPTIONS, &im1, &im2);
- }
-
- if( use_scalerot )
- assert( params.ngh_rad == 0 || !"max trans cannot be used in full scale and rotation mode");
- else
- if( params.subsample_ref && (!ispowerof2(im1->width) || !ispowerof2(im1->height)) ) {
- std_printf("WARNING: first image has dimension which are not power-of-2\n");
- std_printf("For improved results, you should consider resizing the images with '-resize '\n");
- }
-
- // compute deep matching
-// float_image* corres = use_scalerot ?
-// deep_matching_scale_rot( im1, im2, ¶ms, &sr_params ) :
-// deep_matching ( im1, im2, ¶ms, NULL ); // standard call
- cout << "###########################################################" << endl;
- float_image* corres = deep_matching ( im1, im2, ¶ms, NULL ); // standard call
-
- // output
- pl[0] = mxCreateNumericMatrix(corres->ty, corres->tx, mxSINGLE_CLASS, mxREAL);
- corres_to_output(corres, pl[0]);
-
- image_delete(im1);
- image_delete(im2);
- free_image(corres);
- return;
-}
diff --git a/fiji_script_vertical_combination.py b/fiji_script_vertical_combination.py
new file mode 100644
index 0000000..e6d21b4
--- /dev/null
+++ b/fiji_script_vertical_combination.py
@@ -0,0 +1,61 @@
+#@String name
+
+from ij import IJ, ImagePlus
+from ij.io import Opener, FileSaver
+from ij.plugin import RGBStackMerge, StackCombiner, Resizer
+from ij.process import ImageConverter
+from ij.plugin import Resizer
+
+import json
+
+#import file paths
+json_file = open('files.json','r')
+text = json_file.read()
+files = json.loads(text)
+json_file.close()
+
+#Instanciate necessary classes
+opener = Opener()
+RGB_stack_merge = RGBStackMerge()
+stack_combiner = StackCombiner()
+resizer = Resizer()
+
+
+#Open images
+unregistered_GC6s_info = opener.getTiffFileInfo(files['GC6s_input'])
+unregistered_GC6s_imp = opener.openTiffStack(unregistered_GC6s_info)
+
+unregistered_tdTom_info = opener.getTiffFileInfo(files['tdTom_input'])
+unregistered_tdTom_imp = opener.openTiffStack(unregistered_tdTom_info)
+
+registered_GC6s_info = opener.getTiffFileInfo(files['GC6s_registered'])
+registered_GC6s_imp = opener.openTiffStack(registered_GC6s_info)
+
+registered_tdTom_info = opener.getTiffFileInfo(files['tdTom_registered'])
+registered_tdTom_imp = opener.openTiffStack(registered_tdTom_info)
+
+vector_info = opener.getTiffFileInfo(files['vectors_out'])
+vector_imp = opener.openTiffStack(vector_info)
+
+
+#Upsample to resolution of vector image
+IJ.run(unregistered_GC6s_imp, "Size...", "width=%d height=%d depth=%d interpolation=None" % (vector_imp.getWidth(),vector_imp.getHeight(),vector_imp.getStackSize()))
+IJ.run(registered_GC6s_imp, "Size...", "width=%d height=%d depth=%d interpolation=None" % (vector_imp.getWidth(),vector_imp.getHeight(),vector_imp.getStackSize()))
+IJ.run(unregistered_tdTom_imp, "Size...", "width=%d height=%d depth=%d interpolation=None" % (vector_imp.getWidth(),vector_imp.getHeight(),vector_imp.getStackSize()))
+IJ.run(registered_tdTom_imp, "Size...", "width=%d height=%d depth=%d interpolation=None" % (vector_imp.getWidth(),vector_imp.getHeight(),vector_imp.getStackSize()))
+
+
+#Convert all images to RGB
+unregistered_merged_channels_imp = RGBStackMerge.mergeChannels([unregistered_tdTom_imp, unregistered_GC6s_imp, None],False)
+registered_merged_channels_imp = RGBStackMerge.mergeChannels([registered_tdTom_imp, registered_GC6s_imp, None],False)
+IJ.run(unregistered_merged_channels_imp, "RGB Color", "slices")
+IJ.run(registered_merged_channels_imp, "RGB Color", "slices")
+IJ.run(vector_imp, "RGB Color", "slices")
+
+#Combine image into one (stacked on top of each other)
+two_combined_stack = stack_combiner.combineVertically(unregistered_merged_channels_imp.getStack(), vector_imp.getStack())
+three_combined_stack = stack_combiner.combineVertically(two_combined_stack, registered_merged_channels_imp.getStack())
+
+#Save result
+file_saver_combined_image = FileSaver(ImagePlus('combined',three_combined_stack))
+file_saver_combined_image.saveAsTiffStack('/home/aymanns/registration/results/test/combined.tif')
diff --git a/files.json b/files.json
new file mode 100644
index 0000000..337f6b2
--- /dev/null
+++ b/files.json
@@ -0,0 +1,9 @@
+{
+ "tdTom_input": "/home/aymanns/data/1_Horizontal_VNC/tdTom.tif",
+ "GC6s_input": "/home/aymanns/data/1_Horizontal_VNC/GC6s.tif",
+ "tdTom_registered": "/home/aymanns/motion-correction-paper/results/test/warped1.tif",
+ "GC6s_registered": "/home/aymanns/motion-correction-paper/results/test/warped2.tif",
+ "color_flow_out" : "/home/aymanns/motion-correction-paper/results/test/color_flow.tif",
+ "vectors_out": "/home/aymanns/motion-correction-paper/results/test/vectors.tif",
+ "column_out": "/home/aymanns/motion-correction-paper/results/test/column.tif"
+}
diff --git a/motion_compensate.m b/motion_compensate.m
index 696a3dd..fbe102d 100755
--- a/motion_compensate.m
+++ b/motion_compensate.m
@@ -1,74 +1,77 @@
-function motion_compensate(fnInput1,fnInput2,fnMatch,fnDeepMatching,fnOut1,fnOut2,fnColor,N,param)
+function motion_compensate(fnInput1,fnInput2,fnMatch,fnDeepMatching,fnOut1,fnOut2,fnColor,fnVector,N,param)
%% Description
% Motion is estimated with a brigthness constancy data term defined on fnInput1
% and a feature matching similarity term defined on fnMatch. The sequences fnIput1 and
% fnInput2 are warped according to the estimated motion field.
% For more details see the paper: "Imaging neural activity in the ventral
% nerve cord of behaving adult Drosophila", bioRxiv
%
%% Input
% fnInput1: filename of the sequence used for the brightness constancy term, in TIF format
% fnInput2: filename of the sequence warped with the motion field estimated from fnInput1 and fnMatch, in TIF format
% fnMatch: filename of the sequence used for feature matching, in TIF format
% fnDeepMatching: filename of the deepmatching code
% fnOut1: filename used to save the warped version of fnInput1, in TIF format
% fnOut2: filename used to save the warped version of fnInput2, in TIF format
% fnOut2: filename used to save the color visualization of the estimated motion, in TIF format
% N: number of frames to process
% param: parameters of the algorithm (see 'default_parameters.m')
%
% Copyright (C) 2017 D. Fortun, denis.fortun@epfl.ch
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%% Input data
seq=mijread_stack(fnInput1);
seqW=mijread_stack(fnInput2);
seqMatch=mijread_stack(fnMatch);
if N==-1
N=size(seq,3);
end
seq=double(seq(:,:,1:N));
seqRescale=(seq-min(seq(:)))/(max(seq(:))-min(seq(:)))*255;
seqW=double(seqW(:,:,1:N));
seqMatch=double(seqMatch(:,:,1:N));
seqMatch=(seqMatch-min(seqMatch(:)))/(max(seqMatch(:))-min(seqMatch(:)))*255;
%% Motion estimation
w=zeros(size(seqRescale,1),size(seqRescale,2),2,size(seqRescale,3)-1);
colorFlow=zeros(size(seqRescale,1),size(seqRescale,2),3,size(seqRescale,3)-1);
+vectorFlow=zeros(512,512,size(seqRescale,3)-1);
i1=seqRescale(:,:,1);
i1Match=seqMatch(:,:,1);
for t=1:N-1 % Replace parfor by for if you don't want to parallelize
fprintf('Frame %i\n',t);
i2=seqRescale(:,:,t+1);
i2Match=seqMatch(:,:,t+1);
[i10,i2]=midway(i1,i2);
w(:,:,:,t) = compute_motion(i10,i2,i1Match,i2Match,fnDeepMatching,param,t);
colorFlow(:,:,:,t)=flowToColor(w(:,:,:,t));
+ vectorFlow(:,:,t)=flowToVectorImage(w(:,:,:,t),15,[512,512]);
end
%% Registration
seqWarped=seq;
seqwWarped=seqW;
parfor t=1:N-1
seqWarped(:,:,t+1)=warpImg(seq(:,:,t+1),w(:,:,1,t),w(:,:,2,t));
seqwWarped(:,:,t+1)=warpImg(seqW(:,:,t+1),w(:,:,1,t),w(:,:,2,t));
end
%% Save
mijwrite_stack(single(seqWarped),fnOut1);
mijwrite_stack(single(seqwWarped),fnOut2);
mijwrite_stack(single(colorFlow),fnColor,1);
+mijwrite_stack(single(vectorFlow),fnVector);
diff --git a/pipeline.run b/pipeline.run
new file mode 100644
index 0000000..ccb642b
--- /dev/null
+++ b/pipeline.run
@@ -0,0 +1,15 @@
+#!/bin/bash
+#SBATCH -p debug
+#SBATCH --workdir /home/aymanns/motion-correction-paper
+#SBATCH --nodes 1
+#SBATCH --ntasks 1
+#SBATCH --cpus-per-task 1
+#SBATCH --mem 64000
+#SBATCH --time 00:10:00
+module load matlab
+module load jdk
+echo STARTING AT `date`
+matlab -nodesktop -nodisplay -r test
+echo FINISHED Matlab at `date`
+/home/aymanns/Fiji.app/ImageJ-linux64 --ij2 --headless --console --run /home/aymanns/motion-correction-paper/fiji_script_vertical_combination.py
+echo FINISHED Fiji at `date`
diff --git a/test.m b/test.m
index 2378044..9811827 100755
--- a/test.m
+++ b/test.m
@@ -1,37 +1,53 @@
close all;
clear all;
warn = warning ('off','all');
poolobj = gcp('nocreate');
if isempty(poolobj)
parpool;
end
pctRunOnAll warning('off','all')
+%% set path
+fnIJ='/home/aymanns/registration/ij.jar'; % Path to "ij.jar"
+fnMIJ='/home/aymanns/registration/mij.jar'; % Path to "mij.jar"
+javaaddpath(fnIJ);
+javaaddpath(fnMIJ);
+addpath('code');
+addpath('code/external/utils');
+addpath(genpath('code/external/InvPbLib'));
+
+%% load file names
+fname = 'files.json';
+fid = fopen(fname);
+raw = fread(fid,inf);
+str = char(raw');
+fclose(fid);
+files = parse_json(str);
+
%% Parameters
param = default_parameters();
param.lambda = 1000; % Regularization parameter
param.gamma = 100; % Sets the strength of the feature matching constraint
%% Input - Output
-type=0; % 0 for Mac users, 1 for Linux users
-fnMatch='data/1_Horizontal_VNC/tdTom.tif'; % Sequence used for the feature matching similarity term
-fnIn1='data/1_Horizontal_VNC/GC6s.tif'; % Sequence used for the brightness constancy term
-fnIn2='data/1_Horizontal_VNC/tdTom.tif'; % Sequence warped with the motion field estimated from fnIn1 and fnMatch
+type=1; % 0 for Mac users, 1 for Linux users
+fnMatch=files.tdTom_input;%'/home/aymanns/data/1_Horizontal_VNC/tdTom.tif'; % Sequence used for the feature matching similarity term
+fnIn1=files.GC6s_input; % Sequence used for the brightness constancy term
+fnIn2=files.tdTom_input; % Sequence warped with the motion field estimated from fnIn1 and fnMatch
N=3; % Motion is estimated on frames 1 to N. If N=-1, motion is estimated on
% the whole sequence
-fnSave='results/test/'; % Folder where the results are saved
-mkdir(fnSave);
-fnOut1=[fnSave,'warped1.tif']; % Sequence fnIn1 warped
-fnOut2=[fnSave,'warped2.tif']; % Sequence fnIn2 warped
-fnColor=[fnSave,'colorFlow.tif']; % Color visualization of the motion field
+fnOut1=files.GC6s_registered; % Sequence fnIn1 warped
+fnOut2=files.tdTom_registered; % Sequence fnIn2 warped
+fnColor=files.color_flow_out; % Color visualization of the motion field
+fnVector=files.vectors_out;
if type==0
fnDeepMatching='code/external/deepmatching_1.2.2_c++_mac';
elseif type==1
fnDeepMatching='code/external/deepmatching_1.2.2_c++_linux';
end
addpath(fnDeepMatching);
%% Perform motion compensation
-motion_compensate(fnIn1,fnIn2,fnMatch,fnDeepMatching,fnOut1,fnOut2,fnColor,N,param);
+motion_compensate(fnIn1,fnIn2,fnMatch,fnDeepMatching,fnOut1,fnOut2,fnColor,fnVector,N,param);