Allink  v0.1
DrImage.cpp
1 #include "DrImage.h"
2 
3 #ifndef USE_PNG
4 int main(int agrc,char **argv){
5  printf("No pngwriter libs used.\nExit.\n");
6 }
7 #else
8 int main(int argc, char **argv){
9  DrImage *DrI = new DrImage(argc,argv);
10  for(int i=1;i<argc;i++){
11  if(!strncmp(argv[i],"-ga",3)){
12  DrI->Gauss();
13  }
14  else if(!strncmp(argv[i],"-c",2)){
15  DrI->Cut();
16  }
17  else if(!strncmp(argv[i],"-g",2)){
18  DrI->Gravity();
19  }
20  else if(!strncmp(argv[i],"-f",2)){
21  DrI->Fourier();
22  }
23  else if(!strncmp(argv[i],"-i",2)){
24  DrI->Ising();
25  }
26  else if(!strncmp(argv[i],"-l",2)){
27  DrI->LennardJones();
28  }
29  else if(!strncmp(argv[i],"-m",2)){
30  DrI->ConvMatrix();
31  }
32  else if(!strncmp(argv[i],"-s",2)){
33  DrI->Shear();
34  }
35  else if(!strncmp(argv[i],"-r",2)){
36  DrI->Rotor();
37  }
38  else if(!strncmp(argv[i],"-d",2)){
39  DrI->Difference();
40  }
41  else if(!strncmp(argv[i],"-sl",2)){
42  DrI->SlitScan();
43  }
44  else if(!strncmp(argv[i],"-y",2)){
45  }
46  }
47  return 0;
48 }
49 DrImage::DrImage(int argc, char *argv){
50  NFile = 0;
51  NWidth = 0;
52  NHeight = 0;
53  for(int c=0;c<argc;c++){
54  if(!strncmp(argv[i],"-",1)){
55  continue;
56  }
57 
58  NFile++;
59  }
60  FileList = (char **)calloc(NFile,sizeof(char));
61  for(int c=0,f=0;c<argc;c++){
62  if(!strncmp(argv[i],"-",1)){
63  continue;
64  }
65  FileList[f] = (char *)calloc(120,sizeof(char));
66  sprintf(FileList[f],argv[c]);
67  f++;
68  }
69  Load(FileList[0]);
70  Load2(FileList[NFile-1]);
71 }
72 DrImage::~DrImage(){
73  for(int l=0;l<3;l++){
74  free(data[l]);
75  free(data1[l]);
76  free(data);
77  free(data1);
78  for(int f=0;f<NFile;f++){
79  free(FileList[f]);
80  }
81  free(FileList);
82 }
83 void DrImage::Load(char *FName){
84  printf("Load %s\n",FName);
85  ImageIn = new pngwriter(0,0,1.0,FName);
86  ImageIn->readfromfile(FName);
87  NWidth = ImageIn->getwidth();
88  NHeight = ImageIn->getheight();
89  if(NWidth*NHeight == 0){
90  printf("The image is empty\n");
91  exit(0);
92  }
93  data = (double **)calloc(3,sizeof(double));
94  for(int l=0;l<3;l++){
95  data[l] = (double *)calloc(NWidth*NHeight,sizeof(double));
96  }
97  for(int l=0;l<3;l++){
98  for(int w=0;w<NWidth;w++){
99  for(int h=0;h<NHeight;h++){
100  data[l][h*NWidth + w] = ImageIn->dread(w,h,l+1);
101  }
102  }
103  }
104  //ImageIn.close();
105 }
106 void DrImage::ReLoad(char *FName){
107  ImageIn->readfromfile(FName);
108  for(int l=0;l<3;l++){
109  for(int w=0;w<NWidth;w++){
110  for(int h=0;h<NHeight;h++){
111  data[l][h*NWidth + w] = ImageIn->dread(w,h,l+1);
112  }
113  }
114  }
115 }
116 void DrImage::Load2(char *FName){
117  printf("Load %s\n",FName);
118  ImageIn->readfromfile(FName);
119  data1 = (double **)calloc(3,sizeof(double));
120  for(int l=0;l<3;l++){
121  data1[l] = (double *)calloc(NWidth*NHeight,sizeof(double));
122  }
123  for(int l=0;l<3;l++){
124  for(int w=0;w<NWidth;w++){
125  for(int h=0;h<NHeight;h++){
126  data1[l][h*NWidth + w] = ImageIn->dread(w,h,l+1);
127  }
128  }
129  }
130  //ImageIn.close();
131 }
132 void DrImage::Gravity(){
133  int NStep = (int)(25.*60.*.2);
134  double InvNStep = 1./(double)NStep;
135  double InvNHei = 1./(double)NHeight;
136  double Acc = -1.;
137  double Vel = fabs(Acc)*NStep*.75;
138  int Pos = 0.;
139  for(int s=0;s<NStep;s++){
140  fprintf(stderr,"Elaborating %lf %%\r",s*InvNStep*100.);
141  char cImage[160];
142  sprintf(cImage,"Image%05u.png",s);
143  pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
144  for(int h=0;h<NHeight;h++){
145  int h1 = h + Pos;
146  if(h1 >= NHeight) h1 -= NHeight;
147  if(h1 < 0) h1 += NHeight;
148  for(int w=0;w<NWidth;w++){
149  ImageOut.plot(w,h1,data[0][h*NWidth + w],data[1][h*NWidth + w],data[2][h*NWidth + w]);
150  }
151  }
152  Vel += Acc;
153  if(Vel > NHeight) Vel = NHeight;
154  Pos += (int)Vel;
155  Pos -= (int)(floor(Pos*InvNHei)*NHeight);
156  ImageOut.close();
157  }
158  printf("\n");
159 }
160 void DrImage::Ising(){
161  int NStill = 50;
162  int NStep = (int)(25.*60.*.2) - NStill;
163  int NChange = 60;
164  int NSquare = 4;
165  double InvNStep = 1./(double)NStep;
166  Matematica *Mate = new Matematica();
167  for(int s=NStill+NStep;s>=NStep;s--){
168  char cImage[160];
169  sprintf(cImage,"Image%05u.png",s);
170  pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
171  for(int h=0;h<NHeight;h++){
172  for(int w=0;w<NWidth;w++){
173  ImageOut.plot(w,h,data[0][h*NWidth+w],data[1][h*NWidth+w],data[2][h*NWidth+w]);
174  }
175  }
176  ImageOut.close();
177  }
178  for(int s=NStep;s>=0;s--){
179  fprintf(stderr,"Elaborating %lf %%\r",s*InvNStep*100.);
180  char cImage[160];
181  sprintf(cImage,"Image%05u.png",s);
182  pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
183  for(int c=0;c<NChange;c++){
184  int w = (int)(NWidth*Mate->Casuale());
185  int h = (int)(NHeight*Mate->Casuale());
186  w = w - w%NSquare;
187  h = h - h%NSquare;
188  if(Mate->Casuale() < .5){
189  for(int ws=0;ws<NSquare;ws++){
190  if(ws+w > NWidth) continue;
191  for(int hs=0;hs<NSquare;hs++){
192  if(hs+h > NHeight) continue;
193  for(int l=0;l<3;l++){
194  data[l][(h+hs)*NWidth+(w+ws)] = 0.;
195  }
196  }
197  }
198  }
199  else {
200  for(int ws=0;ws<NSquare;ws++){
201  if(ws+w > NWidth) continue;
202  for(int hs=0;hs<NSquare;hs++){
203  if(hs+h > NHeight) continue;
204  for(int l=0;l<3;l++){
205  data[l][(h+hs)*NWidth+(w+ws)] = 1.;
206  }
207  }
208  }
209  }
210  }
211  for(int h=0;h<NHeight;h++){
212  for(int w=0;w<NWidth;w++){
213  ImageOut.plot(w,h,data[0][h*NWidth+w],data[1][h*NWidth+w],data[2][h*NWidth+w]);
214  }
215  }
216  ImageOut.close();
217  }
218  printf("\n");
219 }
220 void DrImage::ConvMatrix(){
221  double Average = 0.;
222  double Count = 0.;
223  Matrice ImIn(NWidth,NHeight);
224  for(int h=0;h<NHeight;h++){
225  for(int w=0;w<NWidth;w++){
226  double Sum = 1.- .3333*(data[0][h*NWidth+w]+data[1][h*NWidth+w]+data[2][h*NWidth+w]);
227  //double Sum = 1. - (data[0][h*NWidth+w]+data[1][h*NWidth+w]+data[2][h*NWidth+w]);
228  ImIn.Set(w,h,Sum);
229  Average += Sum;
230  Count += 1.;
231  }
232  }
233  Average /= Count;
234  //---------Canny---------------------
235  Matematica *Mat = new Matematica;
236  Matrice Canny(5,5);
237  Canny.FillCanny();
238  Canny.Print();
239  // Mat->ApplyFilter(&ImIn,&Canny);
240  // Mat->ApplyFilter(&ImIn,&Canny);
241  //-----------Edge-------------------
242  SPLINE Weight;
243  Weight.a0 = 0.;
244  Weight.a1 = 1.; Weight.a2 = 0.;
245  Weight.a3 = 0.; Weight.a4 = 0.;
246  Matrice *Mask = new Matrice(Weight,3);
247  Mask->Print();
248  //Mat->ApplyFilter(&ImIn,Mask);
249  // Mask->Transpose();
250  // Mat->ApplyFilter(ImIn,Mask);
251  //-------------Smooth----------------
252  const int NMatSize = 5;
253  Matrice GaussEdge(NMatSize,NMatSize);
254  GaussEdge.FillGaussian(.5,1.);
255  //double LatPos[5] = {.125,-.25,0.,.25,-.125};
256  // double LatPos[3] = {-1.,0.,1.};
257  // for(int w=0;w<NMatSize;w++){
258  // for(int h=0;h<NMatSize;h++){
259  // GaussEdge.Set(w,h,GaussEdge.Val(w,h)*LatPos[w]);
260  // }
261  // }
262  GaussEdge.Print();
263  //Mat->ApplyFilter(ImIn,&GaussEdge);
264  Matrice GaussSmooth(5,5);
265  GaussSmooth.FillGaussian(.5,3.);
266  // Mat->ApplyFilter(ImIn,&GaussSmooth);
267  //------------PixDev------------------
268  for(int h=0;h<NHeight;h++){
269  for(int w=0;w<NWidth;w++){
270  ImIn.Set(w,h,Average-ImIn.Val(w,h));
271  }
272  }
273  int PixelDev = 5;
274  double ValStep[3] = {0.,0.,0.};
275  int ValNStep[3] = {0,0,0};
276  for(int h=0;h<NHeight;h++){
277  for(int w=0;w<NWidth;w++){
278  ValNStep[0] = w - PixelDev;
279  if(ValNStep[0] < 0 ) continue;
280  if(ValNStep[0] >= NWidth) continue;
281  ValNStep[1] = w;
282  ValNStep[2] = w + PixelDev;
283  if(ValNStep[2] < 0 ) continue;
284  if(ValNStep[2] >= NWidth) continue;
285  for(int d=0;d<3;d++){
286  ValStep[d] = ImIn.Val(ValNStep[d],h);
287  if(d == 1){
288  ValStep[d] = ValStep[d] > 0. ? ValStep[d] : 0.;
289  continue;
290  }
291  ValStep[d] = ValStep[d] < 0. ? -ValStep[d] : 0.;
292  }
293  double Resp = ValStep[0]*ValStep[1]*ValStep[2];
294  //double Resp = ValStep[1];
295  //ImIn.Set(w,h,Resp);
296  }
297  }
298  char cImage[160];
299  sprintf(cImage,"Canny.png");
300  pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
301  FILE *Ciccia = fopen("Pos3d.dat","w");
302  double NormH = 1./(double)NHeight;
303  double NormW = 1./(double)NWidth;
304  for(int h=0;h<NHeight;h++){
305  for(int w=0;w<NWidth;w++){
306  fprintf(Ciccia,"%lf %lf %lf\n",w*NormH,h*NormH,ImIn.Val(w,h));
307  ImageOut.plot(w,h,ImIn.Val(w,h),ImIn.Val(w,h),ImIn.Val(w,h));
308  }
309  }
310  fclose(Ciccia);
311  ImageOut.close();
312 }
313 void DrImage::Cut(){
314  char cImage[160];
315  sprintf(cImage,"Small.png");
316  double pixel[4];
317  double xBound[2] = {.5,.75};
318  double yBound[2] = {.2,.45};
319  int xNBound[2];
320  int yNBound[2];
321  for(int d=0;d<2;d++){
322  xNBound[d] = (int)(xBound[d]*NWidth);
323  yNBound[d] = (int)(yBound[d]*NHeight);
324  printf("%d %d %d %d\n",xNBound[d],NWidth,yNBound[d],NHeight);
325  }
326  int Dx = xNBound[1]-xNBound[0];
327  int Dy = yNBound[1]-yNBound[0];
328  double *Plot = new double[Dx*Dy];
329  Matematica *Mat = new Matematica;
330  //VarData *Var = new VarData;
331  Matrice *ImIn = new Matrice(Dx,Dy);
332  pngwriter ImageOut(xNBound[1]-xNBound[0],yNBound[1]-yNBound[0],1.0,cImage);
333  double Average = 0.;
334  double Count = 0.;
335  for(int h=yNBound[0];h<yNBound[1];h++){
336  for(int w=xNBound[0];w<xNBound[1];w++){
337  int hh = h -yNBound[0];
338  int ww = w -xNBound[0];
339  double Sum = data[0][h*NWidth+w]+data[1][h*NWidth+w]+data[2][h*NWidth+w];
340  //Sum *= .5;
341  // ImageOut.plot(ww,hh,data[0][h*NWidth+w],data[1][h*NWidth+w],data[2][h*NWidth+w]);
342  ImIn->Set(ww,hh,Sum);
343  Plot[hh*Dx+ww] = Sum;
344  }
345  }
346 
347 
348  for(int h=0;h<Dy;h++){
349  for(int w=0;w<Dx;w++){
350  //if(ImIn->Val(w,h) >= 0.)
351  ImageOut.plot(w,h,10.*ImIn->Val(w,h),ImIn->Val(w,h),ImIn->Val(w,h));
352  }
353  }
354  ImageOut.close();
355  delete [] Plot;
356 }
357 void DrImage::Difference(){
358  char cImage[160];
359  sprintf(cImage,"Difference.png");
360  pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
361  double pixel[4];
362  for(int h=0;h<NHeight;h++){
363  for(int w=0;w<NWidth;w++){
364  for(int l=0;l<3;l++){
365  pixel[l] = data[l][h*NWidth+w] - data1[l][h*NWidth+w];
366  }
367  ImageOut.plot(w,h,pixel[0],pixel[1],pixel[2]);
368  }
369  }
370  ImageOut.close();
371  printf("\n");
372  return;
373  printf("Difference\n");
374  //char cImage[60];
375  double **data2 = (double **)calloc(3,sizeof(double));
376  for(int l=0;l<3;l++){
377  data2[l] = (double *)calloc(NHeight*NWidth,sizeof(double));
378  }
379  for(int h=0;h<NHeight;h++){
380  for(int w=0;w<NWidth;w++){
381  for(int l=0;l<4;l++){
382  data2[l][h*NWidth+w] = data[l][h*NWidth+w];// - data1[l][h*NWidth+w];
383  }
384  }
385  }
386  sprintf(cImage,"Difference.png");
387  //pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
388  for(int h=0;h<NHeight;h++){
389  for(int w=0;w<NWidth;w++){
390  ImageOut.plot(w,h,data2[0][h*NWidth+w],data2[1][h*NWidth+w],data2[2][h*NWidth+w]);
391  }
392  }
393  ImageOut.close();
394 }
395 void DrImage::Shear(){
396  int NStill = 50;
397  int NStep = (int)(25.*60.*.2) - NStill;
398  double InvNStep = 1./(double)NStep;
399  double InvNHei = 1./(double)NHeight;
400  int NStrip = 50;
401  double InvNStrip = 1./(double)NStrip;
402  double *Acc = (double *)calloc(NStrip,sizeof(double));
403  double *Vel = (double *)calloc(NStrip,sizeof(double));
404  int *Pos = (int *)calloc(NStrip,sizeof(int));
405  double NotRational = .1*exp(1)*sqrt(M_PI)/(double)(NStrip);
406  for(int s=0;s<NStrip;s++){
407  Acc[s] = (s*NotRational);
408  Vel[s] = 0.;
409  Pos[s] = 0;
410  }
411  for(int s=NStill+NStep;s>=NStep;s--){
412  char cImage[160];
413  sprintf(cImage,"Image%05u.png",s);
414  pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
415  for(int h=0;h<NHeight;h++){
416  for(int w=0;w<NWidth;w++){
417  ImageOut.plot(w,h,data[0][h*NWidth+w],data[1][h*NWidth+w],data[2][h*NWidth+w]);
418  }
419  }
420  ImageOut.close();
421  }
422  for(int s=NStep;s>=0;s--){
423  //for(int s=0;s<NStep;s++){
424  fprintf(stderr,"Elaborating %lf %%\r",s*InvNStep*100.);
425  char cImage[160];
426  sprintf(cImage,"Image%05u.png",s);
427  pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
428  for(int h=0;h<NHeight;h++){
429  int st = (int)(h*InvNHei*NStrip);
430  for(int w=0;w<NWidth;w++){
431  int w1 = w + Pos[st];
432  if(w1 >= NWidth) w1 -= NWidth;
433  if(w1 < 0) w1 += NWidth;
434  ImageOut.plot(w1,h,data[0][h*NWidth+w],data[1][h*NWidth+w],data[2][h*NWidth+w]);
435  }
436  }
437  for(int st=0;st<NStrip;st++){
438  Vel[st] += Acc[st];
439  //if(Vel[st] > NHeight) Vel[st] = NHeight;
440  Pos[st] += (int)Vel[st];
441  Pos[st] -= (int)(floor(Pos[st]*InvNHei)*NHeight);
442  }
443  ImageOut.close();
444  }
445  printf("\n");
446 }
447 void DrImage::EffectOnDataR(double *data2,int l,int w,int h,int ws,int hs,int NSquare){
448  data[l][(h+hs)*NWidth+(w+ws)] = data2[(NSquare-1-ws)*NSquare+hs];
449 }
450 void DrImage::EffectOnDataT(double *data2,int l,int w,int h,int ws,int hs,int NSquare){
451  data[l][(h+hs)*NWidth+(w+ws)] = data2[ws*NSquare+hs];
452 }
453 void DrImage::EffectOnDataM(double *data2,int l,int w,int h,int ws,int hs,int NSquare){
454  data[l][(h+hs)*NWidth+(w+ws)] = data2[(NSquare-1-hs)*NSquare+ws];
455 }
456 void DrImage::Rotor(){
457  int NStill = 50;//50;
458  int NStep = (int)(25.*60.*.2) - NStill;
459  const int NSquare = 6;
460  int NDecr = (int)((NStep)/(double)NSquare);
461  int NChange = 4;
462  int NSq[7] = {64,49,36,25,16,9,4};
463  int NIncrement = 2;
464  int NSkip = 15;
465  double *data2 = (double *)calloc(NSq[0]*NSq[0],sizeof(double));
466  double InvNStep = 1./(double)NStep;
467  Matematica *Mate = new Matematica();
468  for(int s=NStill+NStep;s>=NStep;s--){
469  char cImage[160];
470  sprintf(cImage,"Image%05u.png",s);
471  pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
472  for(int h=0;h<NHeight;h++){
473  for(int w=0;w<NWidth;w++){
474  ImageOut.plot(w,h,data[0][h*NWidth+w],data[1][h*NWidth+w],data[2][h*NWidth+w]);
475  }
476  }
477  ImageOut.close();
478  }
479  for(int s=NStep,sq=0;s>=0;s--){
480  // if( (s%NDecr) == 0 ) sq++;
481  // if(sq < 0) sq = 0;
482  // if(sq > NSquare) sq = NSquare-1;
483  fprintf(stderr,"Elaborating %lf %%\r",s*InvNStep*100.);
484  char cImage[160];
485  sprintf(cImage,"Image%05u.png",s);
486  pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
487  if((s%NSkip) == 0){
488  for(int c=0;c<NChange;c++){
489  int w = (int)(NWidth*Mate->Casuale());
490  int h = (int)(NHeight*Mate->Casuale());
491  for(int l=0;l<3;l++){
492  for(int ws=0;ws<NSq[sq];ws++){
493  if(ws+w > NWidth) continue;
494  for(int hs=0;hs<NSq[sq];hs++){
495  if(hs+h > NHeight) continue;
496  data2[hs*NSq[sq]+ws] = data[l][(h+hs)*NWidth+(w+ws)];
497  }
498  }
499  for(int ws=0;ws<NSq[sq];ws++){
500  if(ws+w > NWidth) continue;
501  for(int hs=0;hs<NSq[sq];hs++){
502  if(hs+h > NHeight) continue;
503  EffectOnDataR(data2,l,w,h,ws,hs,NSq[sq]);
504  //data[l][(h+hs)*NWidth+(w+ws)] = data2[(NMin-hs)*NSq[sq]+ws];
505  }
506  }
507  }
508  }
509  }
510  for(int h=0;h<NHeight;h++){
511  for(int w=0;w<NWidth;w++){
512  ImageOut.plot(w,h,data[0][h*NWidth+w],data[1][h*NWidth+w],data[2][h*NWidth+w]);
513  }
514  }
515  ImageOut.close();
516  }
517  free(data2);
518  printf("\n");
519 }
520 void DrImage::LennardJones(){
521  //---------------------alloc----------------------
522  int NStill = 50;
523  int NStep = (int)(25.*60.*.2) - NStill;
524  int NRow = 20;
525  int NCol = 20;
526  double NPRow = 1.;//(double)NRow;
527  double NPCol = 1.;//(double)NCol;
528  double NStepPRow = NRow/(double)NStep;
529  double NStepPCol = NCol/(double)NStep;
530  double Eps = 0.1;
531  double Sigma = 3.0;
532  double Dt = 2.;//.1;
533  double InvNStep = 1./(double)NStep;
534  double InvWid = 1./(double)NWidth;
535  double InvHei = 1./(double)NHeight;
536  double InvRow = 1./(double)NRow;
537  double InvCol = 1./(double)NCol;
538  //int HeiSize = (int)(NHeight*InvRow);
539  //int WidSize = (int)(NWidth*InvCol);
540  double **data2 = (double **)calloc(3,sizeof(double));
541  for(int l=0;l<3;l++){
542  data2[l] = (double *)calloc(NHeight*NWidth,sizeof(double));
543  }
544  double *PPos = (double *)calloc(2*NRow*NCol,sizeof(double));
545  int *PBound = (int *)calloc(4*NRow*NCol,sizeof(int));
546  double *PVel = (double *)calloc(2*NRow*NCol,sizeof(double));
547  double *PFor = (double *)calloc(2*NRow*NCol,sizeof(double));
548  double Edge[2] = {NWidth,NHeight};
549  double InvEdge[2] = {1./(double)NWidth,1./(double)NHeight};
550  Matematica *Mate = new Matematica();
551  PERMUTE *Perm = (PERMUTE *)calloc(NRow*NCol,sizeof(PERMUTE));
552  //-------------------------initializing-----------------------
553  for(int l=0;l<3;l++){
554  for(int h=0;h<NHeight;h++){
555  for(int w=0;w<NWidth;w++){
556  data2[l][h*NWidth+w] = data[l][h*NWidth+w];
557  }
558  }
559  }
560  for(int r=0;r<NRow;r++){
561  for(int c=0;c<NCol;c++){
562  double x = (c)*InvCol*NWidth;
563  double y = (r)*InvRow*NHeight;
564  int p1 = (r*NCol+c)*2;
565  PPos[p1 ] = x;
566  PPos[p1+1] = y;
567  PVel[p1 ] = (2.*Mate->Casuale()-1.);
568  PVel[p1+1] = (2.*Mate->Casuale()-1.);
569  PBound[(r*NCol+c)*4 ] = (int)(c*InvCol*NWidth);
570  PBound[(r*NCol+c)*4+1] = (int)((c+1)*InvCol*NWidth);
571  PBound[(r*NCol+c)*4+2] = (int)(r*InvRow*NHeight);
572  PBound[(r*NCol+c)*4+3] = (int)((r+1)*InvRow*NHeight);
573  Perm[r*NCol+c].n = r*NCol+c;
574  Perm[r*NCol+c].m = r*NCol+c;
575  }
576  }
577  //squares of different sizes
578  if(1==1){
579  for(int c=0;c<NCol-1;c++){
580  int NBorder = (int)((2.*Mate->Casuale()-1.)*10.);
581  for(int r=0;r<NRow;r++){
582  int p1 = (r*NCol+c)*4;
583  int p2 = (r*NCol+(c+1))*4;
584  PBound[p1+1] += NBorder;
585  PBound[p2 ] += NBorder;
586  }
587  }
588  for(int r=0;r<NRow-1;r++){
589  int NBorder = (int)((2.*Mate->Casuale()-1.)*10.);
590  for(int c=0;c<NCol;c++){
591  int p1 = (r*NCol+c)*4;
592  int p2 = ((r+1)*NCol+c)*4;
593  PBound[p1+3] += NBorder;
594  PBound[p2+2] += NBorder;
595  }
596  }
597  for(int r=0;r<NRow;r++){
598  for(int c=0;c<NCol;c++){
599  int p1 = (r*NCol+c)*4;
600  double x = PBound[p1 ];
601  double y = PBound[p1+2];
602  int p2 = (r*NCol+c)*2;
603  PPos[p2 ] = x;
604  PPos[p2+1] = y;
605  PPos[p2 ] -= floor(PPos[p2 ]*InvEdge[0])*Edge[0];
606  PPos[p2+1] -= floor(PPos[p2+1]*InvEdge[1])*Edge[1];
607  }
608  }
609  }
610  // for(int c=0;c<NCol;c++){
611  // int p1 = (0*NCol+c)*4;
612  // printf("x %d) %d %d %d\n",c,PBound[p1],PBound[p1+1],PBound[p1+1]-PBound[p1],NWidth);
613  // }
614  // for(int r=0;r<NRow;r++){
615  // int p1 = (r*NCol+0)*4;
616  // printf("y %d) %d %d %d\n",r,PBound[p1+2],PBound[p1+3],PBound[p1+3]-PBound[p1+2],NHeight);
617  // }
618  Mate->PermuteRandomAll(Perm,NRow*NCol);
619  //------------------------loop----------------------------
620  for(int s=NStill+NStep;s>=NStep;s--){
621  char cImage[160];
622  sprintf(cImage,"Image%05u.png",s);
623  pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
624  for(int h=0;h<NHeight;h++){
625  for(int w=0;w<NWidth;w++){
626  ImageOut.plot(w,h,data[0][h*NWidth+w],data[1][h*NWidth+w],data[2][h*NWidth+w]);
627  }
628  }
629  ImageOut.close();
630  }
631  for(int s=NStep;s>=0;s--){
632  NPCol += NStepPCol;
633  if(NPCol > NCol) NPCol = NCol;
634  NPRow += NStepPRow;
635  if(NPRow > NRow) NPRow = NRow;
636  fprintf(stderr,"Elaborating %lf %%\r",s*InvNStep*100.);
637  for(int r=0;r<NRow;r++){
638  for(int c=0;c<NCol;c++){
639  int p1 = (r*NCol+c)*2;
640  PFor[p1 ] = 0.;
641  PFor[p1+1] = 0.;
642  }
643  }
644  for(int r=0;r<NRow;r++){
645  for(int c=0;c<NCol;c++){
646  if(c > (int)NPCol) continue;
647  if(r > (int)NPRow) continue;
648  int p1 = Perm[r*NCol+c].m*2;
649  if(p1 < 0 || p1 >= NRow*NCol*2){
650  continue;
651  }
652  for(int r1=r+1;r1<NRow;r1++){
653  for(int c1=c+1;c1<NCol;c1++){
654  //int p2 = (r1*NCol+c1)*2;
655  int p2 = Perm[r1*NCol+c1].m*2;
656  if(p2 < 0 || p2 >= NRow*NCol) continue;
657  double Dist[3];
658  for(int d=0;d<2;d++){
659  Dist[d] = PPos[p1+d] - PPos[p2+d];
660  Dist[d] -= floor(Dist[d]/Edge[d] + .5)*Edge[d];
661  }
662  double Dist2 = SQR(Dist[0]) + SQR(Dist[1]);
663  if(Dist2 > SQR(3.*Sigma))continue;
664  Dist[2] = sqrt(Dist2);
665  double a = (Sigma/sqrt(Dist2));
666  //double Force = Eps*12.*pow(a,13.) - Eps*6.*pow(a,7.);
667  // if(Force > 20.) Force = 20.;
668  double Force = -Eps/Dist2;
669  for(int d=0;d<2;d++){
670  PFor[p1+d] -= Force*Dist[d]/Dist[2];
671  PFor[p2+d] += Force*Dist[d]/Dist[2];
672  }
673  }
674  }
675  }
676  }
677  //-------------integration--------------
678  for(int r=0;r<NRow;r++){
679  for(int c=0;c<NCol;c++){
680  if(c > (int)NPCol) continue;
681  if(r > (int)NPRow) continue;
682  int p1 = Perm[r*NCol+c].m*2;
683  if(p1 < 0 || p1 >= NRow*NCol*2){
684  continue;
685  }
686  for(int d=0;d<2;d++){
687  PVel[p1+d] += PFor[p1+d]*Dt;
688  PPos[p1+d] += PVel[p1+d]*Dt;
689  PPos[p1+d] -= floor(PPos[p1+d]*InvEdge[d])*Edge[d];
690  }
691  // printf("%d %lf %lf\n",r*NCol+c,PPos[p1],PPos[p1+1]);
692  }
693  }
694  for(int l=0;l<3;l++){
695  for(int h=0;h<NHeight;h++){
696  for(int w=0;w<NWidth;w++){
697  data[l][h*NWidth+w] = 0.;
698  }
699  }
700  }
701  //---------------updating position-------------
702  for(int r=0;r<NRow;r++){
703  for(int c=0;c<NCol;c++){
704  int p1 = (r*NCol+c)*2;
705  int xn = (int)(PPos[p1 ]);
706  int yn = (int)(PPos[p1+1]);
707  int xo = PBound[(r*NCol+c)*4 ];
708  int yo = PBound[(r*NCol+c)*4+2];
709  int WidSize = PBound[(r*NCol+c)*4+1]-PBound[(r*NCol+c)*4 ];
710  int HeiSize = PBound[(r*NCol+c)*4+3]-PBound[(r*NCol+c)*4+2];
711  //printf("%d %d) %d %d -> %d %d |%d %d| %f %f\n",s,r*NCol+c,xo,yo,xn,yn,WidSize,HeiSize,PPos[p1],PPos[p1+1]);
712  for(int l=0;l<3;l++){
713  for(int ws=0;ws<WidSize;ws++){
714  int wo = xo+ws;
715  int wn = xn+ws;
716  if(wo >= NWidth){
717  printf("x out of bound %d > %d | %d\n",wo,NWidth,WidSize);
718  continue;
719  }
720  if(wn >= NWidth) wn -= NWidth;
721  for(int hs=0;hs<HeiSize;hs++){
722  int ho = yo+hs;
723  int hn = yn+hs;
724  if(ho >= NHeight){
725  printf("y out of bound %d+%d = %d > %d |%d\n",yo,hs,ho,NHeight,HeiSize);
726  continue;
727  }
728  if(hn > NHeight) hn -= NHeight;
729  //printf("%d= %d %d) %d->%d (%d %d) %d->%d (%d %d)\n",p1,r,c,xo,xn,WidSize,NWidth,yo,yn,HeiSize,NHeight);
730  if(hn*NWidth+wn < 0 || hn*NWidth+wn >= NWidth*NHeight){
731  printf("input out of border hn %d wn %d > %d %d\n",hn,wn,yn,xn,hn*NWidth+wn,NWidth,NHeight);
732  continue;
733  }
734  if(ho*NWidth+wo < 0 || ho*NWidth+wo >= NWidth*NHeight){
735  printf("output out of border %d %d %d > %d %d\n",hn,wn,hn*NWidth+wn,NWidth,NHeight);
736  continue;
737  }
738  data[l][hn*NWidth+wn] = data2[l][ho*NWidth+wo];
739  }
740  }
741  }
742  }
743  }
744  //-------------saving images
745  char cImage[160];
746  sprintf(cImage,"Image%05u.png",s);
747  pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
748  for(int h=0;h<NHeight;h++){
749  for(int w=0;w<NWidth;w++){
750  ImageOut.plot(w,h,data[0][h*NWidth+w],data[1][h*NWidth+w],data[2][h*NWidth+w]);
751  }
752  }
753  ImageOut.close();
754  }
755  printf("\n");
756  for(int l=0;l<3;l++){
757  free(data2[l]);
758  }
759  free(data2);
760  free(PPos);
761  free(PFor);
762  free(PVel);
763  free(PBound);
764 }
765 #include <fftw3.h>
766 void DrImage::Fourier(){
767  int NStill = 0;//50;
768  int NStep = 10;//(int)(25.*60.*.2) - NStill;
769  int NChange = 60;
770  int NSquare = 2;
771  double Norm = 1./(double)(NWidth*NHeight);
772  double InvNStep = 1./(double)NStep;
773  Matematica *Mate = new Matematica();
774  fftw_complex *out = (fftw_complex *)fftw_malloc(NWidth*NHeight*sizeof(fftw_complex));
775  fftw_complex *in = (fftw_complex *)fftw_malloc(NWidth*NHeight*sizeof(fftw_complex));
776  fftw_plan direct = fftw_plan_dft_2d(NWidth,NHeight,
777  in,out,FFTW_FORWARD,FFTW_PATIENT);
778  fftw_plan reverse = fftw_plan_dft_2d(NWidth,NHeight,
779  out,in,FFTW_BACKWARD,FFTW_PATIENT);
780  int NhInit = NHeight;
781  int NwInit = NWidth;
782  double **data2 = (double **)calloc(3,sizeof(double));
783  for(int l=0;l<3;l++){
784  data2[l] = (double *)calloc(NHeight*NWidth,sizeof(double));
785  }
786  for(int l=0;l<3;l++){
787  for(int h=0;h<NHeight;h++){
788  for(int w=0;w<NWidth;w++){
789  data2[l][h*NWidth+w] = data[l][h*NWidth+w];
790  }
791  }
792  }
793  // for(int s=NStill+NStep;s>=NStep;s--){
794  // char cImage[160];
795  // sprintf(cImage,"Image%05u.png",s);
796  // pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
797  // for(int h=0;h<NHeight;h++){
798  // for(int w=0;w<NWidth;w++){
799  // ImageOut.plot(w,h,data[0][h*NWidth+w],data[1][h*NWidth+w],data[2][h*NWidth+w]);
800  // }
801  // }
802  // ImageOut.close();
803  // }
804  for(int s=NStep;s>0;s--){
805  // NhInit--;
806  // if(NhInit < 0) NhInit = 0;
807  // NwInit--;
808  // if(NwInit < 0) NwInit = 0;
809  fprintf(stderr,"Elaborating %lf %%\r",s*InvNStep*100.);
810  for(int l=0;l<3;l++){
811  for(int h=0;h<NHeight;h++){
812  for(int w=0;w<NWidth;w++){
813  in[h*NWidth+w][0] = data2[l][h*NWidth+w];
814  out[h*NWidth+w][0] = 0.;
815  out[h*NWidth+w][1] = 0.;
816  }
817  }
818  fftw_execute(direct);
819  // for(int h=NhInit;h<NHeight;h++){
820  // for(int w=NwInit;w<NWidth;w++){
821  // printf("%d %d\n",h,w);
822  // out[h*NWidth+w][0] = 0.;
823  // out[h*NWidth+w][1] = 0.;
824  // }
825  // }
826  fftw_execute(reverse);
827  for(int h=0;h<NHeight;h++){
828  for(int w=0;w<NWidth;w++){
829  data[l][h*NWidth+w] = in[h*NWidth+w][0]*Norm;
830  }
831  }
832  }
833  char cImage[160];
834  sprintf(cImage,"Image%05u.png",s);
835  pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
836  for(int h=0;h<NHeight;h++){
837  for(int w=0;w<NWidth;w++){
838  ImageOut.plot(w,h,data[0][h*NWidth+w],data[1][h*NWidth+w],data[2][h*NWidth+w]);
839  }
840  }
841  ImageOut.close();
842  // FILE *FOut = fopen("Result1.dat","w");
843  // for(int h=0;h<NHeight;h++){
844  // for(int w=0;w<NWidth;w++){
845  // double x = w/(double)NWidth;
846  // double y = h/(double)NHeight;
847  // fprintf(FOut,"%lf %lf %lf\n",x,y,data[2][h*NWidth+w]);
848  // }
849  // }
850  // fclose(FOut);
851  // FOut = fopen("Result2.dat","w");
852  // for(int h=0;h<NHeight;h++){
853  // for(int w=0;w<NWidth;w++){
854  // double x = w/(double)NWidth;
855  // double y = h/(double)NHeight;
856  // fprintf(FOut,"%lf %lf %lf\n",x,y,in[h*NWidth+w][0]);
857  // }
858  // }
859  // fclose(FOut);
860  }
861  fftw_destroy_plan(direct);
862  fftw_destroy_plan(reverse);
863  fftw_free(out);
864  fftw_free(in);
865  printf("\n");
866 }
867 void DrImage::Gauss(){
868  int NStill = 0;//50;
869  int NStep = (int)(25.*60.*.2) - NStill;
870  int NSkip = 15;
871  double *Plot2 = (double *)calloc(NWidth*NHeight,sizeof(double));
872  double InvNStep = 1./(double)NStep;
873  Matrice Mask(3,3);
874  Mask.FillGaussian(.5,3.);
875  Mask.Print();
876  int NMat2 = (int)Mask.pNCol()/2;
877  for(int s=NStill+NStep;s>=NStep;s--){
878  char cImage[160];
879  sprintf(cImage,"Image%05u.png",s);
880  pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
881  for(int h=0;h<NHeight;h++){
882  for(int w=0;w<NWidth;w++){
883  ImageOut.plot(w,h,data[0][h*NWidth+w],data[1][h*NWidth+w],data[2][h*NWidth+w]);
884  }
885  }
886  ImageOut.close();
887  }
888  for(int s=NStep,sq=0;s>=0;s--){
889  // if( (s%NDecr) == 0 ) sq++;
890  // if(sq < 0) sq = 0;
891  // if(sq > NSquare) sq = NSquare-1;
892  fprintf(stderr,"Elaborating %lf %%\r",s*InvNStep*100.);
893  char cImage[160];
894  sprintf(cImage,"Image%05u.png",s);
895  pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
896  //if((s%NSkip) == 0)
897  {
898  //convolute filter
899  for(int l=0;l<3;l++){
900  for(int gx=0;gx<NWidth;gx++){
901  for(int gy=0;gy<NHeight;gy++){
902  Plot2[gy*NWidth+gx] = 0.;
903  for(int mx=0;mx<Mask.pNRow();mx++){
904  int g1x = gx + mx - NMat2;
905  if(g1x >= NWidth) continue;//g1x -= NGrid;
906  if(g1x < 0) continue;//g1x + NGrid;
907  for(int my=0;my<Mask.pNCol();my++){
908  int g1y = gy + my - NMat2;
909  if(g1y >= NHeight) continue;//g1y -= NGrid;
910  if(g1y < 0) continue;//g1y + NGrid;
911  Plot2[gy*NWidth+gx] += data[l][g1y*NWidth+g1x]*Mask.Val(mx,my);
912  }
913  }
914  }
915  }
916  for(int gx=0;gx<NWidth;gx++){
917  for(int gy=0;gy<NHeight;gy++){
918  data[l][gy*NWidth+gx] = Plot2[gy*NWidth+gx];
919  }
920  }
921  }
922  }
923  for(int h=0;h<NHeight;h++){
924  for(int w=0;w<NWidth;w++){
925  ImageOut.plot(w,h,data[0][h*NWidth+w],data[1][h*NWidth+w],data[2][h*NWidth+w]);
926  }
927  }
928  ImageOut.close();
929  }
930  free(Plot2);
931  printf("\n");
932 }
933 void DrImage::SlitScan(){
934  int NStep = 1;
935  pngwriter ImageOut[NStep];
936  for(int f=0;f<NStep;f++){
937  char cImage[160];
938  sprintf(cImage,"Image%05u.png",f);
939  ImageOut[f] = new pngwriter(NWidth,NHeight,1.0,cImage);
940  }
941  for(int f=0;f<NFile;f++){
942  fprintf(stderr,"Elaborating %lf %%\r",f/(double)NFile*100.);
943  ReLoad(FileList[f]);
944  int wInit = f;
945  if(wInit >= NWidth) break;//wInit -= NWidth;
946  for(int w=wInit;w<NWidth;w++){
947  int f1 = 0;
948  for(int h=0;h<NHeight;h++){
949  //int set = f + h*NFile;
950  ImageOut[f1].plot(w,h,data[0][h*NWidth + w],data[1][h*NWidth + w],data[2][h*NWidth + w]);
951  }
952  }
953  }
954  for(int s=0;s<NStep;s++){
955  ImageOut[s].close();
956  }
957  printf("\n");
958 }
959 #endif //USE_PNG
double a4
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
int pNCol()
Number of columns.
int m
Second number of the couple.
Indices for the permutation.
int PermuteRandomAll(int *Sequence, int NMass)
Permutes a sequence without repeating.
int n
First number of the couple.
double a3
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
double Casuale()
Random uniform number.
double Val(int row)
Returns a value in 1d.
void FillGaussian(double Sigma, double CutOff)
Fill the entries for the Gauss blur.
double a2
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
double a1
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
Matrice computes the algebric operations on matrices.
bool Set(int row, int column, double Val)
Set a coefficient.
Implementation of useful algorythms.
Definition: Matematica.h:76
int pNRow()
Number of rows.
double a0
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
Coefficient of a spline.
void Print()
Print the entries.