4 int main(
int agrc,
char **argv){
5 printf(
"No pngwriter libs used.\nExit.\n");
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)){
14 else if(!strncmp(argv[i],
"-c",2)){
17 else if(!strncmp(argv[i],
"-g",2)){
20 else if(!strncmp(argv[i],
"-f",2)){
23 else if(!strncmp(argv[i],
"-i",2)){
26 else if(!strncmp(argv[i],
"-l",2)){
29 else if(!strncmp(argv[i],
"-m",2)){
32 else if(!strncmp(argv[i],
"-s",2)){
35 else if(!strncmp(argv[i],
"-r",2)){
38 else if(!strncmp(argv[i],
"-d",2)){
41 else if(!strncmp(argv[i],
"-sl",2)){
44 else if(!strncmp(argv[i],
"-y",2)){
49 DrImage::DrImage(
int argc,
char *argv){
53 for(
int c=0;c<argc;c++){
54 if(!strncmp(argv[i],
"-",1)){
60 FileList = (
char **)calloc(NFile,
sizeof(
char));
61 for(
int c=0,f=0;c<argc;c++){
62 if(!strncmp(argv[i],
"-",1)){
65 FileList[f] = (
char *)calloc(120,
sizeof(
char));
66 sprintf(FileList[f],argv[c]);
70 Load2(FileList[NFile-1]);
78 for(
int f=0;f<NFile;f++){
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");
93 data = (
double **)calloc(3,
sizeof(
double));
95 data[l] = (
double *)calloc(NWidth*NHeight,
sizeof(
double));
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);
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);
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));
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);
132 void DrImage::Gravity(){
133 int NStep = (int)(25.*60.*.2);
134 double InvNStep = 1./(double)NStep;
135 double InvNHei = 1./(double)NHeight;
137 double Vel = fabs(Acc)*NStep*.75;
139 for(
int s=0;s<NStep;s++){
140 fprintf(stderr,
"Elaborating %lf %%\r",s*InvNStep*100.);
142 sprintf(cImage,
"Image%05u.png",s);
143 pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
144 for(
int h=0;h<NHeight;h++){
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]);
153 if(Vel > NHeight) Vel = NHeight;
155 Pos -= (int)(floor(Pos*InvNHei)*NHeight);
160 void DrImage::Ising(){
162 int NStep = (int)(25.*60.*.2) - NStill;
165 double InvNStep = 1./(double)NStep;
167 for(
int s=NStill+NStep;s>=NStep;s--){
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]);
178 for(
int s=NStep;s>=0;s--){
179 fprintf(stderr,
"Elaborating %lf %%\r",s*InvNStep*100.);
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());
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.;
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.;
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]);
220 void DrImage::ConvMatrix(){
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]);
244 Weight.
a1 = 1.; Weight.
a2 = 0.;
245 Weight.
a3 = 0.; Weight.
a4 = 0.;
252 const int NMatSize = 5;
253 Matrice GaussEdge(NMatSize,NMatSize);
254 GaussEdge.FillGaussian(.5,1.);
265 GaussSmooth.FillGaussian(.5,3.);
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));
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;
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);
288 ValStep[d] = ValStep[d] > 0. ? ValStep[d] : 0.;
291 ValStep[d] = ValStep[d] < 0. ? -ValStep[d] : 0.;
293 double Resp = ValStep[0]*ValStep[1]*ValStep[2];
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));
315 sprintf(cImage,
"Small.png");
317 double xBound[2] = {.5,.75};
318 double yBound[2] = {.2,.45};
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);
326 int Dx = xNBound[1]-xNBound[0];
327 int Dy = yNBound[1]-yNBound[0];
328 double *Plot =
new double[Dx*Dy];
332 pngwriter ImageOut(xNBound[1]-xNBound[0],yNBound[1]-yNBound[0],1.0,cImage);
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];
342 ImIn->
Set(ww,hh,Sum);
343 Plot[hh*Dx+ww] = Sum;
348 for(
int h=0;h<Dy;h++){
349 for(
int w=0;w<Dx;w++){
351 ImageOut.plot(w,h,10.*ImIn->
Val(w,h),ImIn->
Val(w,h),ImIn->
Val(w,h));
357 void DrImage::Difference(){
359 sprintf(cImage,
"Difference.png");
360 pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
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];
367 ImageOut.plot(w,h,pixel[0],pixel[1],pixel[2]);
373 printf(
"Difference\n");
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));
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];
386 sprintf(cImage,
"Difference.png");
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]);
395 void DrImage::Shear(){
397 int NStep = (int)(25.*60.*.2) - NStill;
398 double InvNStep = 1./(double)NStep;
399 double InvNHei = 1./(double)NHeight;
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);
411 for(
int s=NStill+NStep;s>=NStep;s--){
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]);
422 for(
int s=NStep;s>=0;s--){
424 fprintf(stderr,
"Elaborating %lf %%\r",s*InvNStep*100.);
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]);
437 for(
int st=0;st<NStrip;st++){
440 Pos[st] += (int)Vel[st];
441 Pos[st] -= (int)(floor(Pos[st]*InvNHei)*NHeight);
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];
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];
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];
456 void DrImage::Rotor(){
458 int NStep = (int)(25.*60.*.2) - NStill;
459 const int NSquare = 6;
460 int NDecr = (int)((NStep)/(double)NSquare);
462 int NSq[7] = {64,49,36,25,16,9,4};
465 double *data2 = (
double *)calloc(NSq[0]*NSq[0],
sizeof(
double));
466 double InvNStep = 1./(double)NStep;
468 for(
int s=NStill+NStep;s>=NStep;s--){
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]);
479 for(
int s=NStep,sq=0;s>=0;s--){
483 fprintf(stderr,
"Elaborating %lf %%\r",s*InvNStep*100.);
485 sprintf(cImage,
"Image%05u.png",s);
486 pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
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)];
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]);
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]);
520 void DrImage::LennardJones(){
523 int NStep = (int)(25.*60.*.2) - NStill;
528 double NStepPRow = NRow/(double)NStep;
529 double NStepPCol = NCol/(double)NStep;
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;
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));
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};
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];
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;
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;
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;
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;
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;
605 PPos[p2 ] -= floor(PPos[p2 ]*InvEdge[0])*Edge[0];
606 PPos[p2+1] -= floor(PPos[p2+1]*InvEdge[1])*Edge[1];
620 for(
int s=NStill+NStep;s>=NStep;s--){
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]);
631 for(
int s=NStep;s>=0;s--){
633 if(NPCol > NCol) NPCol = NCol;
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;
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){
652 for(
int r1=r+1;r1<NRow;r1++){
653 for(
int c1=c+1;c1<NCol;c1++){
655 int p2 = Perm[r1*NCol+c1].
m*2;
656 if(p2 < 0 || p2 >= NRow*NCol)
continue;
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];
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));
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];
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){
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];
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.;
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];
712 for(
int l=0;l<3;l++){
713 for(
int ws=0;ws<WidSize;ws++){
717 printf(
"x out of bound %d > %d | %d\n",wo,NWidth,WidSize);
720 if(wn >= NWidth) wn -= NWidth;
721 for(
int hs=0;hs<HeiSize;hs++){
725 printf(
"y out of bound %d+%d = %d > %d |%d\n",yo,hs,ho,NHeight,HeiSize);
728 if(hn > NHeight) hn -= 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);
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);
738 data[l][hn*NWidth+wn] = data2[l][ho*NWidth+wo];
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]);
756 for(
int l=0;l<3;l++){
766 void DrImage::Fourier(){
771 double Norm = 1./(double)(NWidth*NHeight);
772 double InvNStep = 1./(double)NStep;
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;
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));
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];
804 for(
int s=NStep;s>0;s--){
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.;
818 fftw_execute(direct);
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;
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]);
861 fftw_destroy_plan(direct);
862 fftw_destroy_plan(reverse);
867 void DrImage::Gauss(){
869 int NStep = (int)(25.*60.*.2) - NStill;
871 double *Plot2 = (
double *)calloc(NWidth*NHeight,
sizeof(
double));
872 double InvNStep = 1./(double)NStep;
876 int NMat2 = (int)Mask.
pNCol()/2;
877 for(
int s=NStill+NStep;s>=NStep;s--){
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]);
888 for(
int s=NStep,sq=0;s>=0;s--){
892 fprintf(stderr,
"Elaborating %lf %%\r",s*InvNStep*100.);
894 sprintf(cImage,
"Image%05u.png",s);
895 pngwriter ImageOut(NWidth,NHeight,1.0,cImage);
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;
906 if(g1x < 0)
continue;
907 for(
int my=0;my<Mask.
pNCol();my++){
908 int g1y = gy + my - NMat2;
909 if(g1y >= NHeight)
continue;
910 if(g1y < 0)
continue;
911 Plot2[gy*NWidth+gx] += data[l][g1y*NWidth+g1x]*Mask.
Val(mx,my);
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];
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]);
933 void DrImage::SlitScan(){
935 pngwriter ImageOut[NStep];
936 for(
int f=0;f<NStep;f++){
938 sprintf(cImage,
"Image%05u.png",f);
939 ImageOut[f] =
new pngwriter(NWidth,NHeight,1.0,cImage);
941 for(
int f=0;f<NFile;f++){
942 fprintf(stderr,
"Elaborating %lf %%\r",f/(
double)NFile*100.);
945 if(wInit >= NWidth)
break;
946 for(
int w=wInit;w<NWidth;w++){
948 for(
int h=0;h<NHeight;h++){
950 ImageOut[f1].plot(w,h,data[0][h*NWidth + w],data[1][h*NWidth + w],data[2][h*NWidth + w]);
954 for(
int s=0;s<NStep;s++){
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.
int pNRow()
Number of rows.
double a0
a0 + a1*x + a2*x^2 + a3*x^3 + a4^4
void Print()
Print the entries.