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);