.....죽는줄 알았다
화상 클러스터링이랑 이미지에 있는 꽃 개수 체크하는 프로그램보단 훨씬 쉽지만..
프로그래밍이 약한 나에게는 엄청 고역이였다 ....ㅠㅠㅠ
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include"/home1/prof/kudo/jikken/data_inout.h"
#define N 100
#define pi 3.141592
#define fmin 0.98
#define fmax 1.06
//8p (32)式
double a2theta(double A,double B,double theta){
double x;
x = theta * pi / 180;
return A * A * cos(x) * cos(x) + B * B * sin(x) * sin(x);
}
//8p (31)式
double p1(double A,double B,double lo,double r,double theta){
double aa2theta,r2,sqrted,result;
if(r < 0)
r = r * -1;
r2 = r * r;
aa2theta = a2theta(A,B,theta);
sqrted = sqrt(aa2theta - r2);
if(r2 <= aa2theta)
result = 2 * lo * A * B / aa2theta * sqrted;
else
result = 0;
return result;
}
//8p (30)式
double p(double x0,double y0,double A,double B,double alpha,double lo,double r,double theta){
double x,r_new,theta_new,result;
x = theta * 3.1416 / 180;
r_new = r - x0 * cos(x) - y0 * sin(x);
theta_new = theta - alpha;
result = p1(A,B,lo,r_new,theta_new);
return result;
}
//8p (29)式 Shepp-Loganの畳み込み関数
double h(double m,double delta){
return 4/( pi*delta*delta * (1-4*m*m));
}
/*
//Ramachandran-Lakshminarayanan
double h(double m,double delta){
int temp_m;
temp_m = (int)m;
if(temp_m == 0){
return (double)pi / (2*delta*delta);
}
else if(temp_m%2 != 0){
return (double)-1*2/(m*m*pi*delta*delta);
}
else
return 0;
}
*/
//7p (27)式
int m0(double k,double l,int n){
int result;
result = (int)floor(k*cos(n*pi/N) + l*sin(n*pi/N));
//printf("%d ",result);
return result;
}
//7p (26)式
double gamm(double k,double l,int n){
return k * cos(n * pi / N) + l * sin(n * pi / N) - m0(k,l,n);
}
//7p (25)式
double q1(double k,double l,int n,double q_data[][N]){
int am0;
double agamma,result;
agamma = gamm(k,l,n);
am0 = m0(k,l,n);
if(am0 > 127)
am0 = 127;
else if(am0 < -127)
am0 = -127;
am0 = am0 + 127;
result = (1-agamma) * q_data[am0][n] + agamma * q_data[am0+1][n];
return result;
}
//7p (24)式
double f(double k,double l,double q_data[][N]){
int n;
double f_data_temp,result;
f_data_temp = 0;
for(n=0;n<N;n++){
f_data_temp += q1(k,l,n,q_data);
}
result = (double)f_data_temp/(2*N);
// printf("%f ",result);
return result;
}
//10p (33)式
int to_seisu(double f){
if(f >= fmax){
return 255;
}
else if(f <= fmin){
return 0;
}
else{
return (int)floor((256 * (f - fmin) / (fmax - fmin)));
}
}
int main(void){
FILE *fo;
char *fname="dataforSL.csv";
unsigned char imgout[256][256];
int size = 256*256;
char filename[] = "outputforSL";
double ar,atheta,p_data_temp,delta;
int i,j,k,l,kk,ll;
double mmdash,temp2,mm,sumforq;
int m,n,mdash;
double p_data[256][256];
double q_data[256][256];
double f_data[256][256];
int i_data[256][256];
//Shepp-Logan head phantom data
double SL[10][6] = {
// {x0,y0,A,B,alpha,lo}
{0,0,0.92,0.69,90,2.0},
{0,-0.0184,0.874,0.6624,90,-0.98},
{0.22,0,0.31,0.11,72,-0.02},
{-0.22,0,0.41,0.16,108,-0.02},
{0,0.35,0.25,0.21,90,0.01},
{0,0.1,0.046,0.046,0,0.01},
{0,-0.1,0.046,0.046,0,0.01},
{-0.08,-0.605,0.046,0.023,0,0.01},
{0,-0.605,0.023,0.023,0,0.01},
{0.06,-0.605,0.046,0.023,90,0.01}
};
// initialization for p(theta,r)
for(i=0;i<256;i++){
for(j=0;j<N;j++){
p_data[i][j] = 0;
q_data[i][j] = 0;
}
}
// for p(m,n)
for(i=0;i<10;i++){
for(j=0;j<256;j++){
ar = (double)(j - 127) / 127;
for(k=0;k<N;k++){
atheta = (double) k / N * 180;
p_data_temp = p(SL[i][0],SL[i][1],SL[i][2],SL[i][3],SL[i][4],SL[i][5],ar,atheta);
p_data[j][k] += p_data_temp;
}
}
}
// for q(m',n)
delta = (double)2.0/256;
for(mdash=0;mdash<256;mdash++){
mmdash = mdash-127;
for(n=0;n<N;n++){
sumforq = 0;
for(m=0;m<256;m++){
mm = m-127;
temp2 = (double)mmdash-mm;
sumforq += (double)h(temp2,delta)*p_data[m][n];
}
q_data[mdash][n] = (double)sumforq * delta;
}
}
//f(k,l) >> imgout
for(k=0;k<256;k++){
for(l=0;l<256;l++){
kk = (k - 127);
ll = (l - 127);
f_data[k][l] = f(kk,ll,q_data);
i_data[k][l] = to_seisu(f_data[k][l]);
imgout[k][l] = i_data[k][l];
}
}
write_image(imgout,size,filename);
if((fo=fopen(fname,"w"))==NULL){
printf("File[%s] does not open!\n",fname);
exit(0);
}