その他にも当然問題があるのですが、誰もが犯す類の過ちなので、個条書きだけして終りにしま す。リスト４−１に示していない部分に関することも書いておきます。
| リスト４−３ 修正版 |
1 /*****************************************************************************/
2 /* ハッチング処理 */
3 /* by BASIC Programmer */
4 /* */
5 /* 変更 1991/07/31 fuji 大幅に変更いたしました。 */
6 /*****************************************************************************/
7 #include <stdio.h>
8 #include <math.h>
9
10 #define MAXSIZE 100
11
12 typedef struct {
13 long x, y, z;
14 } point3 ;
15
16 typedef struct {
17 double x, y, z;
18 } point3d ;
19
20 typedef double matrix3[3][3];
21
22 typedef struct {
23 point3 origin; /* 基準点 */
24 point3 direction; /* 方向 */
25 int n; /* データ数 */
26 point3 positions[MAXSIZE]; /* 座標 */
27 } hatching;
28
29 typedef struct {
30 double cost, sint, costh[4], sinth[4];
31 } elements;
32
33 /*---------- static関数プロトタイプ宣言 ----------*/
34
35 static void get_cos_sin(double x1, double y1, double x2, double y2,
36 double *costh, double *sinth);
37 static void convert1(hatching *hin, elements *hw);
38 static double takasa(long gpx[], long gpy[], long gpz[], double x0, double y0);
39
40 #define f_swap(a,b) { double tmp=a; a=b; b=tmp }
41 #define sgn(a) ((a)==0.0 ? 0.0 : (a)<0.0 ? -1.0 : 1.0)
42
43 /*---------------------------------------------------------------------------*/
44 /* 回転角度の cos, sin を求める */
45 /*---------------------------------------------------------------------------*/
46 static void get_cos_sin(double x1, double y1, double x2, double y2,
47 double *costh, double *sinth)
48 {
49 double d, w1, w2;
50
51 d = hypot( w1 = x2 - x1, w2 = y2 - y1 );
52 if (1.0 > d) {
53 *costh = 1.0;
54 *sinth = 0.0;
55 } else {
56 *costh = w1 / d;
57 *sinth = w2 / d;
58 }
59 }
60
61 /*---------------------------------------------------------------------------*/
62 /* 行列＝行列 */
63 /*---------------------------------------------------------------------------*/
64 static void matcpy( matrix3 out, matrix3 in )
65 {
66 int i, j;
67
68 for( i=0 ; i<3 ; ++i )
69 for( j=0 ; j<3 ; ++j )
70 out[i][j] = in[i][j];
71 }
72
73 /*---------------------------------------------------------------------------*/
74 /* 行列＝行列・行列 */
75 /*---------------------------------------------------------------------------*/
76 static void matmul( matrix3 out, matrix3 m1, matrix3 m2 )
77 {
78 int i, j, k;
79 matrix3 m;
80 double s;
81
82 for( i=0 ; i<3 ; ++i )
83 for( j=0 ; j<3 ; ++j ) {
84 s = 0;
85 for( k=0 ; k<3 ; ++k )
86 s += m1[i][k] * m2[k][j];
87 m[i][j] = s;
88 }
89
90 matcpy( out, m );
91 }
92
93 /*---------------------------------------------------------------------------*/
94 /* 実数ベクトル＝行列・実数ベクトル */
95 /*---------------------------------------------------------------------------*/
96 vecmat( point3d *outp, matrix3 mat, point3d *inp)
97 {
98 point3d wk;
99
100 wk.x = mat[0][0]*inp->x + mat[0][1]*inp->y + mat[0][2]*inp->z;
101 wk.y = mat[1][0]*inp->x + mat[1][1]*inp->y + mat[1][2]*inp->z;
102 wk.z = mat[2][0]*inp->x + mat[2][1]*inp->y + mat[2][2]*inp->z;
103
104 *outp = wk;
105 }
106
107 /*---------------------------------------------------------------------------*/
108 /* 整数ベクトル＝行列・整数ベクトル */
109 /*---------------------------------------------------------------------------*/
110 vecmat_i( point3 *outp, matrix3 mat, point3 *inp)
111 {
112 point3 wk;
113
114 wk.x = mat[0][0]*inp->x + mat[0][1]*inp->y + mat[0][2]*inp->z;
115 wk.y = mat[1][0]*inp->x + mat[1][1]*inp->y + mat[1][2]*inp->z;
116 wk.z = mat[2][0]*inp->x + mat[2][1]*inp->y + mat[2][2]*inp->z;
117
118 *outp = wk;
119 }
120
121 /*---------------------------------------------------------------------------*/
122 /* convert1 */
123 /*---------------------------------------------------------------------------*/
124 /* fuji : 作成者の本当の「意図」が分からないので、無理矢理同じ動作する */
125 /* ようにしました。本当は、もっと、はるかにすっきりなるでしょう。 */
126 /*---------------------------------------------------------------------------*/
127 static void convert1(hatching *hin, elements *hw)
128 {
129 matrix3 rot, r;
130 static matrix3 unitmat = { {1,0,0}, {0,1,0}, {0,0,1} };
131 point3d pos[3];
132 int i, td;
133 double sinth, costh;
134 point3 *p;
135
136 /* 先頭３点の座標値をdoubleで取り出す */
137 for( i=0 ; i<3 ; ++i ) {
138 pos[i].x = hin->positions[i].x;
139 pos[i].y = hin->positions[i].y;
140 pos[i].z = hin->positions[i].z;
141 }
142
143 /* ｚ軸まわりの回転を行なう */
144 td = 0;
145 get_cos_sin( pos[0].x, pos[0].y, pos[1].x, pos[1].y, &costh, &sinth );
146 matcpy( rot, unitmat );
147 rot[0][0] = rot[1][1] = hw->costh[td] = costh ;
148 rot[1][0] = -(rot[0][1] = hw->sinth[td] = sinth);
149 for( i=0 ; i<3 ; ++i )
150 vecmat( &pos[i], rot, &pos[i] );
151
152 /* ｙ軸まわりの回転を行なう */
153 td = 1;
154 get_cos_sin( pos[0].x, pos[0].z, pos[1].x, pos[1].z, &costh, &sinth );
155 matcpy( r, unitmat );
156 r[0][0] = r[2][2] = hw->costh[td] = costh ;
157 r[2][0] = -(r[0][2] = hw->sinth[td] = sinth);
158 for( i=0 ; i<3 ; ++i )
159 vecmat( &pos[i], r, &pos[i] );
160 matmul( rot, r, rot );
161
162 /* ｘ軸まわりの回転を行なう */
163 td = 2;
164 get_cos_sin( pos[1].y, pos[1].z, pos[2].y, pos[2].z, &costh, &sinth );
165 matcpy( r, unitmat );
166 r[1][1] = r[2][2] = hw->costh[td] = costh ;
167 r[2][1] = -(r[1][2] = hw->sinth[td] = sinth);
168 matmul( rot, r, rot );
169
170 /* 求まった回転行列を、回転すべき全対象に適用する */
171 p = &hin->positions[0];
172 for (i = 0; i < hin->n; i++, ++p )
173 vecmat_i( p, rot, p );
174 vecmat_i( &hin->origin, rot, &hin->origin );
175 vecmat_i( &hin->direction, rot, &hin->direction );
176 }
177
178 /*---------------------------------------------------------------------------*/
179 /* 高さ */
180 /*---------------------------------------------------------------------------*/
181 /* fuji : 内容が良く分からないので、最低限だけ直しています。ごめんね。 */
182 /*---------------------------------------------------------------------------*/
183 static double takasa(long gpx[], long gpy[], long gpz[], double x0, double y0)
184 {
185 int tr, i1, i2, k, j, i;
186 double iy1, iy2, dd, d1, sinr, z, gcost, gsint,
187 xd, yd, kx1, kx2, kz1, kz2, gminx;
188 double xdd[3], ydd[3], zdd[3], ixd[3], iyd[3], izd[3];
189 double d;
190
191 iy1 = 10000.0;
192 iy2 = -10000.0;
193 for( tr=0; ; ++tr ) {
194 i1 = tr;
195 i2 = (tr+1) % 3;
196 d = hypot((double)(gpx[i1]-gpy[i1]), (double)(gpx[i2]-gpy[i2]));
197 if( 1.0 <= d ) {
198 gcost = (gpx[i2] - gpx[i1]) / d;
199 gsint = (gpy[i2] - gpy[i1]) / d;
200 gminx = 99999.0;
201 for (j = 0; j < 3; j++ ) {
202 k = (tr + j) % 3;
203 xdd[j] = gcost * gpx[k] + gsint * gpy[k];
204 ydd[j] = -gsint * gpx[k] + gcost * gpy[k];
205 zdd[j] = (double)gpz[k];
206 if(xdd[j] < gminx)
207 gminx = xdd[j];
208 }
209 }
210 xd = gcost * x0 + gsint * y0;
211 yd = -gsint * x0 + gcost * y0;
212 if( (int)gminx == (int)xd )
213 continue;
214 d1 = (iy1 - iy2);
215 sinr = (iy2 - iy1) / d1;
216 for (i = 0; i < 3; i++) {
217 ixd[i] = ((zdd[i] - iy1) / d1) * sinr;
218 iyd[i] = -((xdd[i] - xd ) / d1) * sinr;
219 izd[i] = ydd[i] - yd;
220 }
221 if( iyd[1]==iyd[0] || iyd[2]==iyd[0] )
222 continue;
223
224 kx1 = ixd[0] - (ixd[1]-ixd[0]) / (iyd[1]-iyd[0]) * iyd[0];
225 kz1 = izd[0] - (izd[1]-izd[0]) / (iyd[1]-iyd[0]) * iyd[0];
226 kx2 = ixd[0] - (ixd[2]-ixd[0]) / (iyd[2]-iyd[0]) * iyd[0];
227 kz2 = izd[0] - (izd[2]-izd[0]) / (iyd[2]-iyd[0]) * iyd[0];
228 if (kz2 == kz1 )
229 continue;
230 dd = kx1 - (kx2 - kx1) / (kz2 - kz1) * kz1;
231 z = dd * d1 * sinr + iy1;
232 return(z);
233 }
234 }
235
236 /*---------------------------------------------------------------------------*/
237 /* 途中までです。ごめんなさーい。 by fuji */
238 /*---------------------------------------------------------------------------*/
|