add freak kernel init

This commit is contained in:
VahdPourAmin
2021-08-26 09:35:15 +02:00
parent 50192e3b04
commit c8cebce678
133 changed files with 378 additions and 14079 deletions
+8 -2
View File
@@ -19,6 +19,9 @@
#include <vio/global.h>
#include <vio/matcher.h>
#include <CL/cl.h>
#include <vio/cl_class.h>
#include <vio/initialization.h>
namespace vk {
class AbstractCamera;
@@ -44,7 +47,7 @@ public:
size_t max_n_kfs; //!< max number of keyframes to reproject from
bool find_match_direct;
Options()
: max_n_kfs(10),
: max_n_kfs(30),
find_match_direct(true)
{}
} options_;
@@ -52,7 +55,7 @@ public:
size_t n_matches_;
size_t n_trials_;
Reprojector(vk::AbstractCamera* cam, Map& map);
Reprojector(vk::AbstractCamera* cam, Map& map, opencl* kernel);
~Reprojector();
@@ -94,6 +97,9 @@ private:
void resetGrid();
bool reprojectCell(Cell& cell, FramePtr frame);
bool reprojectPoint(FramePtr frame, Point* point);
protected:
opencl* gpu_freak;
};
} // namespace vio
+129
View File
@@ -0,0 +1,129 @@
void meanIntensity( image2d_t _image, image2d_t _integral, float kp_x, float kp_y,
int scale, int rot, int* point, float3* patternLookupPtr, uchar* output)
{
// Prepare a suitable OpenCL image sampler.
sampler_t const sampler = CLK_ADDRESS_CLAMP | CLK_FILTER_NEAREST;
// get point position in image
float3* FreakPoint = patternLookupPtr + (scale * FREAK_NB_ORIENTATION * FREAK_NB_POINTS + rot * FREAK_NB_POINTS + (*point));
float xf = FreakPoint->x + kp_x;
float yf = FreakPoint->y + kp_y;
int x = (int)xf;
int y = (int)yf;
// get the sigma:
float radius = FreakPoint->z;//sigma
// calculate output:
if( radius < 0.5 )
{
// interpolation multipliers:
const int r_x = (int)((xf - x) * 1024);
const int r_y = (int)((yf - y) * 1024);
const int r_x_1 = (1024 - r_x);
const int r_y_1 = (1024 - r_y);
uint ret_val;
// linear interpolation:
ret_val = (r_x_1 * r_y_1 * (int)(read_imageui(_image, sampler, (int2)(y , x ))).x)
+ (r_x * r_y_1 * (int)(read_imageui(_image, sampler, (int2)(y , x+1))).x)
+ (r_x_1 * r_y * (int)(read_imageui(_image, sampler, (int2)(y+1, x ))).x)
+ (r_x * r_y * (int)(read_imageui(_image, sampler, (int2)(y+1, x+1))).x);
// return the rounded mean
ret_val += 2 * 1024 * 1024;
output[*point] = (uchar)(ret_val / (4 * 1024 * 1024));
}
// calculate borders
const int x_left = round(xf - radius);
const int y_top = round(yf - radius);
const int x_right = round(xf + radius + 1); // integral image is 1px wider
const int y_bottom = round(yf + radius + 1); // integral image is 1px higher
int ret_val;
ret_val = read_imageui(_integral, sampler, (int2)(y_bottom, x_right)).x; // bottom right corner
ret_val -= read_imageui(_integral, sampler, (int2)(y_bottom, x_left)).x; // bottom left corner
ret_val += read_imageui(_integral, sampler, (int2)(y_top, x_left)).x; // top left corner
ret_val -= read_imageui(_integral, sampler, (int2)(y_top, x_right)).x; // top right corner
const int area = (x_right - x_left) * (y_bottom - y_top);
ret_val = (ret_val + area / 2) / area;
output[*point] = (uchar)(ret_val);
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
__kernel void freak(
__global image2d_t* image, // grayescale image
__global image2d_t* imgIntegral,
__global float4* keypoints, // ptx, pty, angle
__global uchar* descriptors, // kp_size * (64 = NB_PAIRS/8) - CV_8U with zero values
__global int* patternSizes, // array size : NB_SCALES = 64
__global int4* orientationPairs, // array size : FREAK_NB_ORIENTATION = 256
__global float3* patternLookup, // array size : NB_SCALES(64) * FREAK_NB_ORIENTATION(256) * FREAK_NB_POINTS(43)
__global uchar2* descriptionPairs, // array size : NB_PAIRS = 512
__global int* nOctaves,
__global int* kpScaleIdx //array size = kp_size - used to save pattern scale index corresponding to each keypoints
)
{
const float sizeCst = (float)(NB_SCALES / (FREAK_LOG2 * *nOctaves));
uchar pointsValue[FREAK_NB_POINTS];
int thetaIdx = 0;
int direction0;
int direction1;
int k = get_global_id(0);
kpScaleIdx[k] = max( (int)(log(7.0 / FREAK_SMALLEST_KP_SIZE) * sizeCst + 0.5) ,0);
if( kpScaleIdx[k] >= NB_SCALES )
kpScaleIdx[k] = NB_SCALES - 1;
if( keypoints[k].x <= patternSizes[kpScaleIdx[k]] || // check if the description at this specific position and scale fits inside the image
keypoints[k].y <= patternSizes[kpScaleIdx[k]] ||
keypoints[k].x >= get_image_dim(image[(int)keypoints[k].w]).y - patternSizes[kpScaleIdx[k]] ||
keypoints[k].y >= get_image_dim(image[(int)keypoints[k].w]).x - patternSizes[kpScaleIdx[k]])
return;
// get the points intensity value in the un-rotated pattern
for( int i = FREAK_NB_POINTS; i--; )
{
meanIntensity(image[(int)keypoints[k].w], imgIntegral[(int)keypoints[k].w], keypoints[k].x, keypoints[k].y,
*(kpScaleIdx + k), 0, &i, patternLookup, pointsValue);
}
direction0 = 0;
direction1 = 0;
for( int m = 45; m--; )
{
// iterate through the orientation pairs
const int delta = (pointsValue[ orientationPairs[m].w ] - pointsValue[ orientationPairs[m].x ]);
direction0 += delta * (orientationPairs[m].y) / 2048; // weight_dx
direction1 += delta * (orientationPairs[m].z) / 2048; // weight_dy
}
thetaIdx = round(FREAK_NB_ORIENTATION *
(float)(atan2((float)direction1,(float)direction0) * (180.0 / 3.141592653589793238462643383279502884197169399375)) *
(1 / 360.0));
if( thetaIdx < 0 )
thetaIdx += FREAK_NB_ORIENTATION;
if( thetaIdx >= FREAK_NB_ORIENTATION )
thetaIdx -= FREAK_NB_ORIENTATION;
// extract descriptor at the computed orientation
for( int i = FREAK_NB_POINTS; i--; )
{
meanIntensity(image[(int)keypoints[k].w], imgIntegral[(int)keypoints[k].w], keypoints[k].x, keypoints[k].y,
*(kpScaleIdx + k), thetaIdx, &i, patternLookup, pointsValue);
}
// Extract descriptor
int cnt = 0;
for( int n = 0; n<4; n++ )
{
for (int m = 8; m >= 1; m--)
{
for(int kk = (n*16)+15; kk >= (n*16); --kk, ++cnt)
{
descriptors[(int)keypoints[k].w*300*64 + (k * 64) + kk] *= 2;
if (pointsValue[(int)descriptionPairs[cnt].x] >= pointsValue[(int)descriptionPairs[cnt].y])
descriptors[(int)keypoints[k].w*300*64 + (k * 64) + kk] += 1;
}
}
}
}
+14
View File
@@ -0,0 +1,14 @@
__kernel void integral_image(
__read_only image2d_t* int_img, // grayescale image - height X width
__write_only image2d_t* out_img // (height+1) X (width+1) with zero values
)
{
sampler_t const sampler = CLK_ADDRESS_CLAMP | CLK_FILTER_NEAREST;
int k = get_global_id(0)
int2 coord = (int2)(get_global_id(1)+1, get_global_id(2)+1);
int res = 0;
for (int i = 0; i < coord.x; i++)
for (int j = 0; j < coord.y; j++)
res += read_imageui(int_img[k], sampler, (int2)(i, j)); // bottom right corner
write_imageui(out_img[k], coord, res);
}
+3 -1
View File
@@ -41,7 +41,9 @@ opencl::opencl(vk::AbstractCamera* cam):cam(cam) {
" -DF_Y="+std::to_string(camera[1])+
" -DC_X="+std::to_string(camera[2])+
" -DC_Y="+std::to_string(camera[3])+
" -DS="+std::to_string(camera[4]);
" -DS="+std::to_string(camera[4])+
" -DFREAK_LOG2=0.693147180559945 -DFREAK_NB_ORIENTATION=256 -DFREAK_NB_POINTS=43"+
" -DFREAK_SMALLEST_KP_SIZE=7 -DNB_PAIRS=512 -DNB_SCALES=64";
if(program->build({ *device },options.c_str()) !=0)
std::cout << " Error building: " << program->getBuildInfo<CL_PROGRAM_BUILD_STATUS>(*device)<<'\n'
<< " Binary type: " << program->getBuildInfo<CL_PROGRAM_BINARY_TYPE>(*device)<<'\n'
@@ -171,6 +171,7 @@ FrameHandlerBase::UpdateResult FrameHandlerMono::processFrame()
SparseImgAlignGpu img_align(Config::kltMaxLevel(), Config::kltMinLevel(),30, SparseImgAlignGpu::GaussNewton, false,gpu_fast_);
if(img_align.run(last_frame_, new_frame_)==0)return RESULT_FAILURE;
reprojector_.reprojectMap(new_frame_, overlap_kfs_);
// assert // TODO
/*
std::cout<<"Reprojection Map nPoint: "<<overlap_kfs_.back().second
<<"\tnCell: "<<reprojector_.n_trials_<<"\t nMatches: "<<reprojector_.n_matches_
+223 -11
View File
@@ -30,8 +30,8 @@
namespace vio {
Reprojector::Reprojector(vk::AbstractCamera* cam, Map& map) :
map_(map)
Reprojector::Reprojector(vk::AbstractCamera* cam, Map& map, opencl* kernel) :
map_(map), gpu_freak(kernel)
{
initializeGrid(cam);
}
@@ -130,15 +130,227 @@ void Reprojector::reprojectMap(
// Now we go through each grid cell and select one point to match.
// At the end, we should have at maximum one reprojected point per cell.
//SVO_START_TIMER("feature_align");
for(size_t i=0; i<grid_.cells.size(); ++i)
{
// we prefer good quality points over unkown quality (more likely to match)
// and unknown quality over candidates (position not optimized)
if(reprojectCell(*grid_.cells.at(grid_.cell_order[i]), frame))
++n_matches_;
if(n_matches_ > (size_t) Config::maxFts())
break;
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
std::vector<cl_float4> kps; // keypoints
cl_float4 data;
int frame_number = 0;
for(auto it_frame = overlap_kfs.begin(), ite_frame = overlap_kfs.end(); it_frame != ite_frame; it_frame++){
for(auto it_ftr = it_frame->first->fts_.begin(), ite_ftr = it_frame->first->fts_.end(); it_ftr != ite_ftr; it_ftr++){
data.w = frame_number;
data.x = (*it_ftr)->px.x();
data.y = (*it_ftr)->px.y();
data.z = 0.0;
kps.push_back(data);
}
frame_number++;
}
cl_float4 keypoints[kps.size()];
size_t counter = 0;
for(auto it = kps.begin(), ite = kps.end(); it != ite; it++) {
cl_float4 ftr;
ftr.w = it->w;
ftr.x = it->x;
ftr.y = it->y;
ftr.z = it->z;
keypoints[counter++] = ftr;
}
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
/*cl_float3 keypoints[overlap_kfs.back().second];
size_t counter = 0;
for(size_t i = 0; i < grid_.cells.size(); ++i)
{
Cell* cell = grid_.cells.at(grid_.cell_order[i]);
for(auto it = cell->begin(), ite = cell->end(); it != ite; it++, counter++)
{
keypoints[counter].x = it->px.x(); //x
keypoints[counter].y = it->px.y(); //y
keypoints[counter].z = 0.0; //angle
}
}*/
// opencl_initializer();
// cv::Mat img = overlap_kfs.back().first->img();
// cv::Mat integral_img=cv::Mat(cv::Size(701, 501),CV_8UC1);
// gpu_freak->reload_buf(0,0,img);
// gpu_freak->run(0);
// gpu_freak->read(0,1,1,&integral_img);
// FREAK preprocessing
static const float FREAK_LOG2 = 0.693147180559945;
static const int FREAK_NB_ORIENTATION = 256;
static const int NB_ORIENPAIRS = 45;
static const int FREAK_NB_POINTS = 43;
static const int FREAK_SMALLEST_KP_SIZE = 7; // smallest size of keypoints
static const int NB_PAIRS = 512;
static const int NB_SCALES = 64;
bool orientationNormalized = true;
bool scaleNormalized = true;
float patternScale = 22.0f;
int nOctaves = 4;
bool extAll = false; // true if all pairs need to be extracted for pairs selection
cl_float3 patternLookup[NB_SCALES * FREAK_NB_ORIENTATION * FREAK_NB_POINTS]; // look-up table for the pattern points (position+sigma of all points at all scales and orientation)
int patternSizes[NB_SCALES]; // size of the pattern at a specific scale (used to check if a point is within image boundaries)
cl_uchar2 descriptionPairs[NB_PAIRS];
cl_int4 orientationPairs[NB_ORIENPAIRS];
// default pairs
static const int FREAK_DEF_PAIRS[NB_PAIRS] = {
404,431,818,511,181,52,311,874,774,543,719,230,417,205,11,
560,149,265,39,306,165,857,250,8,61,15,55,717,44,412,
592,134,761,695,660,782,625,487,549,516,271,665,762,392,178,
796,773,31,672,845,548,794,677,654,241,831,225,238,849,83,
691,484,826,707,122,517,583,731,328,339,571,475,394,472,580,
381,137,93,380,327,619,729,808,218,213,459,141,806,341,95,
382,568,124,750,193,749,706,843,79,199,317,329,768,198,100,
466,613,78,562,783,689,136,838,94,142,164,679,219,419,366,
418,423,77,89,523,259,683,312,555,20,470,684,123,458,453,833,
72,113,253,108,313,25,153,648,411,607,618,128,305,232,301,84,
56,264,371,46,407,360,38,99,176,710,114,578,66,372,653,
129,359,424,159,821,10,323,393,5,340,891,9,790,47,0,175,346,
236,26,172,147,574,561,32,294,429,724,755,398,787,288,299,
769,565,767,722,757,224,465,723,498,467,235,127,802,446,233,
544,482,800,318,16,532,801,441,554,173,60,530,713,469,30,
212,630,899,170,266,799,88,49,512,399,23,500,107,524,90,
194,143,135,192,206,345,148,71,119,101,563,870,158,254,214,
276,464,332,725,188,385,24,476,40,231,620,171,258,67,109,
844,244,187,388,701,690,50,7,850,479,48,522,22,154,12,659,
736,655,577,737,830,811,174,21,237,335,353,234,53,270,62,
182,45,177,245,812,673,355,556,612,166,204,54,248,365,226,
242,452,700,685,573,14,842,481,468,781,564,416,179,405,35,
819,608,624,367,98,643,448,2,460,676,440,240,130,146,184,
185,430,65,807,377,82,121,708,239,310,138,596,730,575,477,
851,797,247,27,85,586,307,779,326,494,856,324,827,96,748,
13,397,125,688,702,92,293,716,277,140,112,4,80,855,839,1,
413,347,584,493,289,696,19,751,379,76,73,115,6,590,183,734,
197,483,217,344,330,400,186,243,587,220,780,200,793,246,824,
41,735,579,81,703,322,760,720,139,480,490,91,814,813,163,
152,488,763,263,425,410,576,120,319,668,150,160,302,491,515,
260,145,428,97,251,395,272,252,18,106,358,854,485,144,550,
131,133,378,68,102,104,58,361,275,209,697,582,338,742,589,
325,408,229,28,304,191,189,110,126,486,211,547,533,70,215,
670,249,36,581,389,605,331,518,442,822
};
float scaleStep = pow(2.0, (float)(nOctaves)/NB_SCALES ); // 2 ^ ( (nOctaves-1) /nbScales)
const int nn[8] = {6, 6, 6, 6, 6, 6, 6, 1}; // number of points on each concentric circle (from outer to inner)
const float bigR = (float)(2.0 / 3.0); // bigger radius
const float smallR = (float)(2.0 / 24.0); // smaller radius
const float unitSpace = (float)((bigR - smallR) / 21.0 ); // define spaces between concentric circles (from center to outer: 1,2,3,4,5,6)
// radii of the concentric cirles (from outer to inner)
const float radius[8] = {bigR , bigR - 6 * unitSpace, bigR - 11 * unitSpace, bigR - 15 * unitSpace,
bigR - 18 * unitSpace, bigR - 20 * unitSpace, smallR , 0.0};
// sigma of pattern points (each group of 6 points on a concentric cirle has the same sigma)
const float sigma[8] = {float(radius[0] / 2.0), float(radius[1] / 2.0), float(radius[2] / 2.0), float(radius[3] / 2.0),
float(radius[4] / 2.0), float(radius[5] / 2.0), float(radius[6] / 2.0), float(radius[6] / 2.0)};
for( int scaleIdx=0; scaleIdx < NB_SCALES; ++scaleIdx )
{
patternSizes[scaleIdx] = 0; // proper initialization
float scalingFactor = pow(scaleStep,scaleIdx); //scale of the pattern, scaleStep ^ scaleIdx
for( int orientationIdx = 0; orientationIdx < FREAK_NB_ORIENTATION; ++orientationIdx )
{
int pointIdx = 0;
for( size_t i = 0; i < 8; ++i )
{
for( int k = 0 ; k < nn[i]; ++k )
{
float alpha = (float)(k) * 2 * (float)3.141592653589793238462643383279502884197169399375 / (float)(nn[i])
+ ((float)3.141592653589793238462643383279502884197169399375 / nn[i] * ( i % 2 )/*beta*/)
+ ((float)(orientationIdx * 2 * (float)3.141592653589793238462643383279502884197169399375 / (float)(FREAK_NB_ORIENTATION)) /* theta*/);
// add the point to the look-up table
cl_float3* point = &patternLookup[scaleIdx * FREAK_NB_ORIENTATION * FREAK_NB_POINTS + orientationIdx * FREAK_NB_POINTS + pointIdx];
point->x = (float)(radius[i] * cos(alpha) * scalingFactor * patternScale); //x
point->y = (float)(radius[i] * sin(alpha) * scalingFactor * patternScale); //y
point->z = (float)(sigma[i] * scalingFactor * patternScale); //sigma
// adapt the sizeList if necessary
const int sizeMax = (int)(ceil((radius[i] + sigma[i]) * scalingFactor * patternScale)) + 1;
if( patternSizes[scaleIdx] < sizeMax )
patternSizes[scaleIdx] = sizeMax;
++pointIdx;
}
}
}
}
// build the list of orientation pairs
orientationPairs[0].w=0; orientationPairs[0].x=3; orientationPairs[1].w=1; orientationPairs[1].x=4; orientationPairs[2].w=2; orientationPairs[2].x=5;
orientationPairs[3].w=0; orientationPairs[3].x=2; orientationPairs[4].w=1; orientationPairs[4].x=3; orientationPairs[5].w=2; orientationPairs[5].x=4;
orientationPairs[6].w=3; orientationPairs[6].x=5; orientationPairs[7].w=4; orientationPairs[7].x=0; orientationPairs[8].w=5; orientationPairs[8].x=1;
orientationPairs[9].w=6; orientationPairs[9].x=9; orientationPairs[10].w=7; orientationPairs[10].x=10; orientationPairs[11].w=8; orientationPairs[11].x=11;
orientationPairs[12].w=6; orientationPairs[12].x=8; orientationPairs[13].w=7; orientationPairs[13].x=9; orientationPairs[14].w=8; orientationPairs[14].x=10;
orientationPairs[15].w=9; orientationPairs[15].x=11; orientationPairs[16].w=10; orientationPairs[16].x=6; orientationPairs[17].w=11; orientationPairs[17].x=7;
orientationPairs[18].w=12; orientationPairs[18].x=15; orientationPairs[19].w=13; orientationPairs[19].x=16; orientationPairs[20].w=14; orientationPairs[20].x=17;
orientationPairs[21].w=12; orientationPairs[21].x=14; orientationPairs[22].w=13; orientationPairs[22].x=15; orientationPairs[23].w=14; orientationPairs[23].x=16;
orientationPairs[24].w=15; orientationPairs[24].x=17; orientationPairs[25].w=16; orientationPairs[25].x=12; orientationPairs[26].w=17; orientationPairs[26].x=13;
orientationPairs[27].w=18; orientationPairs[27].x=21; orientationPairs[28].w=19; orientationPairs[28].x=22; orientationPairs[29].w=20; orientationPairs[29].x=23;
orientationPairs[30].w=18; orientationPairs[30].x=20; orientationPairs[31].w=19; orientationPairs[31].x=21; orientationPairs[32].w=20; orientationPairs[32].x=22;
orientationPairs[33].w=21; orientationPairs[33].x=23; orientationPairs[34].w=22; orientationPairs[34].x=18; orientationPairs[35].w=23; orientationPairs[35].x=19;
orientationPairs[36].w=24; orientationPairs[36].x=27; orientationPairs[37].w=25; orientationPairs[37].x=28; orientationPairs[38].w=26; orientationPairs[38].x=29;
orientationPairs[39].w=30; orientationPairs[39].x=33; orientationPairs[40].w=31; orientationPairs[40].x=34; orientationPairs[41].w=32; orientationPairs[41].x=35;
orientationPairs[42].w=36; orientationPairs[42].x=39; orientationPairs[43].w=37; orientationPairs[43].x=40; orientationPairs[44].w=38; orientationPairs[44].x=41;
for( unsigned m = NB_ORIENPAIRS; m--; )
{
float dx = patternLookup[orientationPairs[m].w].x - patternLookup[orientationPairs[m].x].x;
float dy = patternLookup[orientationPairs[m].w].y - patternLookup[orientationPairs[m].x].y;
float norm_sq = (dx * dx + dy * dy);
orientationPairs[m].y = round((dx / (norm_sq)) * 4096.0); //weight_dx
orientationPairs[m].z = round((dy / (norm_sq)) * 4096.0); //weight_dy
}
// build the list of description pairs
cl_uchar2 allPairs[(int)((FREAK_NB_POINTS * (FREAK_NB_POINTS + 1)) / 2)];
for( int i = 1, idx = 0; i < FREAK_NB_POINTS; ++i )
{
// (generate all the pairs)
for( int j = 0; j < i; ++j )
{
allPairs[idx++] = (cl_uchar2){(uchar)i, (uchar)j};
}
}
for( int i = 0; i < NB_PAIRS; ++i )
descriptionPairs[i] = allPairs[FREAK_DEF_PAIRS[i]];
int kpScaleIdx[overlap_kfs.back().second];
// gpu_freak->reload_buf(1,0,img);
// gpu_freak->reload_buf(1,1,integral_img);
// gpu_freak->reload_buf(1,2, keypoints);
// gpu_freak->reload_buf(1,4,patternSizes);
// gpu_freak->reload_buf(1,5,orientationPairs);
// gpu_freak->reload_buf(1,6,patternLookup);
// gpu_freak->reload_buf(1,7,descriptionPairs);
// gpu_freak->write_buf(1,8,nOctaves);
// gpu_freak->reload_buf(1,9,kpScaleIdx);
// gpu_freak->run(1);
// cl_uchar descriptor[frame_number*300*64]={0};
// gpu_freak->read(1,3,300*64,descriptor);
//////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// for(size_t i=0; i<grid_.cells.size(); ++i)
// {
// // we prefer good quality points over unkown quality (more likely to match)
// // and unknown quality over candidates (position not optimized)
// if(reprojectCell(*grid_.cells.at(grid_.cell_order[i]), frame))
// ++n_matches_;
// if(n_matches_ > (size_t) Config::maxFts())
// break;
// }
// SVO_STOP_TIMER("feature_align");
}
@@ -151,7 +363,7 @@ bool Reprojector::pointQualityComparator(Candidate& lhs, Candidate& rhs)
bool Reprojector::reprojectCell(Cell& cell, FramePtr frame)
{
cell.sort(boost::bind(&Reprojector::pointQualityComparator, _1, _2));
cell.sort(boost::bind(&Reprojector::pointQualityComparator, _1, _2));
Cell::iterator it=cell.begin();
while(it!=cell.end())
{