polygon.cpp 3.7 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165
  1. #include "polygon.hpp"
  2. //---------------
  3. // Polygon::init
  4. //---------------
  5. void
  6. Polygon::init(){
  7. for(size_t i=0;i<grid_size;++i){
  8. grid[i]=empty;
  9. grid_histo[i]=0;
  10. }
  11. for(size_t i=0;i<grid_width;++i){
  12. grid[pos(i,0)]=border;
  13. grid[pos(i,grid_height-1)]=border;
  14. }
  15. for(size_t i=0;i<grid_height;++i){
  16. grid[pos(0,i)]=border;
  17. grid[pos(grid_width-1,i)]=border;
  18. }
  19. for(size_t i=0;i<xinit;++i){
  20. grid[pos(i,1)]=border;
  21. }
  22. head=pos(xinit,yinit);
  23. grid[head]=1;
  24. }
  25. //--------------------
  26. // Polygon::operator=
  27. //--------------------
  28. Polygon&
  29. Polygon::operator=(const Polygon& P){
  30. head=P.head;
  31. length=P.length;
  32. memcpy(grid,P.grid,grid_size*sizeof(size_t));
  33. memcpy(grid_histo,P.grid_histo,grid_size*sizeof(size_t));
  34. memcpy(step_histo,P.step_histo,(max_len+1)*sizeof(size_t));
  35. return *this;
  36. }
  37. //------------------
  38. // Polygon::display
  39. //------------------
  40. void
  41. Polygon::display() const{
  42. for(size_t z=0;z<grid_height;++z){
  43. size_t y=grid_height-z-1;
  44. for(size_t x=0;x<grid_width;++x){
  45. size_t v=grid[pos(x,y)];
  46. if(v==border) cout<<"\033[47m \033[0m";
  47. else if(v==empty or step_histo[v]!=grid_histo[pos(x,y)] or v>length) cout<<" ";
  48. else cout<<"\033[42m\033[30m"<<" "<<"\033[0m";
  49. }
  50. cout<<endl;
  51. }
  52. }
  53. //-------------
  54. // Polygon::fp
  55. //-------------
  56. Reel
  57. Polygon::fp() const{
  58. map<size_t,InfoVertex> infos;
  59. size_t id=0;
  60. int64 x=xinit;
  61. int64 y=yinit;
  62. size_t c=pos(x,y);
  63. for(size_t t=1;t<=length;++t){
  64. size_t l=left(c);
  65. size_t r=right(c);
  66. size_t u=up(c);
  67. size_t d=down(c);
  68. auto it=infos.find(c);
  69. if(it==infos.end()){
  70. infos[c]=InfoVertex(id++,x,y,l,r,u,d);
  71. }
  72. else{
  73. InfoVertex& info=it->second;
  74. info.v[0]=l;
  75. info.v[1]=r;
  76. info.v[2]=u;
  77. info.v[3]=d;
  78. info.deg=4;
  79. }
  80. it=infos.find(l);
  81. if(it==infos.end()) infos[l]=InfoVertex(id++,x-1,y);
  82. it=infos.find(r);
  83. if(it==infos.end()) infos[r]=InfoVertex(id++,x+1,y);
  84. it=infos.find(u);
  85. if(it==infos.end()) infos[u]=InfoVertex(id++,x,y+1);
  86. it=infos.find(d);
  87. if(it==infos.end()) infos[d]=InfoVertex(id++,x,y-1);
  88. size_t s=step_histo[t+1];
  89. if(grid[l]==t+1 and s==grid_histo[l]){
  90. c=l;
  91. --x;
  92. }
  93. else if(grid[r]==t+1 and s==grid_histo[r]){
  94. c=r;
  95. ++x;
  96. }
  97. else if(grid[u]==t+1 and s==grid_histo[u]){
  98. c=u;
  99. ++y;
  100. }
  101. else if(grid[d]==t+1 and s==grid_histo[d]){
  102. c=d;
  103. --y;
  104. }
  105. }
  106. size_t ne=infos.size();
  107. AvxMatrix AvxB;
  108. AvxMatrix AvxC;
  109. AvxB.clear();
  110. AvxC.clear();
  111. for(auto it=infos.begin();it!=infos.end();++it){
  112. InfoVertex& info=it->second;
  113. size_t i=info.id;
  114. for(size_t k=0;k<info.deg;++k){
  115. size_t j=infos[info.v[k]].id;
  116. AvxB.get(i,j)=1;
  117. AvxB.get(j,i)=1;
  118. }
  119. for(auto it2=infos.begin();it2!=infos.end();++it2){
  120. InfoVertex& info2=it2->second;
  121. size_t j=info2.id;
  122. int64 dx=abs(info.x-info2.x);
  123. int64 dy=abs(info.y-info2.y);
  124. AvxC.get(i,j)=get_coeff(dx,dy);
  125. }
  126. }
  127. AvxMatrix AvxM;
  128. AvxM.clear();
  129. AvxM.from_C_B(AvxC,AvxB,ne);
  130. Reel det=AvxM.Gauss(ne,ne+1);
  131. Reel res=0;
  132. for(size_t i=0;i<ne;++i){
  133. res+=AvxB.get_diag_square_sym(i,ne)*AvxM.get(i,ne);
  134. }
  135. res*=(det*0.25);
  136. return res;
  137. }
  138. //----------------------
  139. // Polygon::next_vertex
  140. //----------------------
  141. size_t
  142. Polygon::next_vertex(size_t t,size_t c) const{
  143. size_t l=left(c);
  144. size_t r=right(c);
  145. size_t u=up(c);
  146. size_t d=down(c);
  147. size_t s=step_histo[t+1];
  148. if(grid[l]==t+1 and s==grid_histo[l]) return l;
  149. if(grid[r]==t+1 and s==grid_histo[r]) return r;
  150. if(grid[u]==t+1 and s==grid_histo[u]) return u;
  151. if(grid[d]==t+1 and s==grid_histo[d]) return d;
  152. assert(false);
  153. return 0;
  154. }