1 /*
2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
3 % %
4 % %
5 % %
6 % FFFFF EEEEE AAA TTTTT U U RRRR EEEEE %
7 % F E A A T U U R R E %
8 % FFF EEE AAAAA T U U RRRR EEE %
9 % F E A A T U U R R E %
10 % F EEEEE A A T UUU R R EEEEE %
11 % %
12 % %
13 % MagickCore Image Feature Methods %
14 % %
15 % Software Design %
16 % Cristy %
17 % July 1992 %
18 % %
19 % %
20 % Copyright 1999-2020 ImageMagick Studio LLC, a non-profit organization %
21 % dedicated to making software imaging solutions freely available. %
22 % %
23 % You may not use this file except in compliance with the License. You may %
24 % obtain a copy of the License at %
25 % %
26 % https://imagemagick.org/script/license.php %
27 % %
28 % Unless required by applicable law or agreed to in writing, software %
29 % distributed under the License is distributed on an "AS IS" BASIS, %
30 % WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. %
31 % See the License for the specific language governing permissions and %
32 % limitations under the License. %
33 % %
34 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35 %
36 %
37 %
38 */
39
40 /*
41 Include declarations.
42 */
43 #include "MagickCore/studio.h"
44 #include "MagickCore/animate.h"
45 #include "MagickCore/artifact.h"
46 #include "MagickCore/blob.h"
47 #include "MagickCore/blob-private.h"
48 #include "MagickCore/cache.h"
49 #include "MagickCore/cache-private.h"
50 #include "MagickCore/cache-view.h"
51 #include "MagickCore/channel.h"
52 #include "MagickCore/client.h"
53 #include "MagickCore/color.h"
54 #include "MagickCore/color-private.h"
55 #include "MagickCore/colorspace.h"
56 #include "MagickCore/colorspace-private.h"
57 #include "MagickCore/composite.h"
58 #include "MagickCore/composite-private.h"
59 #include "MagickCore/compress.h"
60 #include "MagickCore/constitute.h"
61 #include "MagickCore/display.h"
62 #include "MagickCore/draw.h"
63 #include "MagickCore/enhance.h"
64 #include "MagickCore/exception.h"
65 #include "MagickCore/exception-private.h"
66 #include "MagickCore/feature.h"
67 #include "MagickCore/gem.h"
68 #include "MagickCore/geometry.h"
69 #include "MagickCore/list.h"
70 #include "MagickCore/image-private.h"
71 #include "MagickCore/magic.h"
72 #include "MagickCore/magick.h"
73 #include "MagickCore/matrix.h"
74 #include "MagickCore/memory_.h"
75 #include "MagickCore/module.h"
76 #include "MagickCore/monitor.h"
77 #include "MagickCore/monitor-private.h"
78 #include "MagickCore/morphology-private.h"
79 #include "MagickCore/option.h"
80 #include "MagickCore/paint.h"
81 #include "MagickCore/pixel-accessor.h"
82 #include "MagickCore/profile.h"
83 #include "MagickCore/property.h"
84 #include "MagickCore/quantize.h"
85 #include "MagickCore/quantum-private.h"
86 #include "MagickCore/random_.h"
87 #include "MagickCore/resource_.h"
88 #include "MagickCore/segment.h"
89 #include "MagickCore/semaphore.h"
90 #include "MagickCore/signature-private.h"
91 #include "MagickCore/string_.h"
92 #include "MagickCore/thread-private.h"
93 #include "MagickCore/timer.h"
94 #include "MagickCore/utility.h"
95 #include "MagickCore/version.h"
96
97 /*
98 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
99 % %
100 % %
101 % %
102 % C a n n y E d g e I m a g e %
103 % %
104 % %
105 % %
106 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
107 %
108 % CannyEdgeImage() uses a multi-stage algorithm to detect a wide range of
109 % edges in images.
110 %
111 % The format of the CannyEdgeImage method is:
112 %
113 % Image *CannyEdgeImage(const Image *image,const double radius,
114 % const double sigma,const double lower_percent,
115 % const double upper_percent,ExceptionInfo *exception)
116 %
117 % A description of each parameter follows:
118 %
119 % o image: the image.
120 %
121 % o radius: the radius of the gaussian smoothing filter.
122 %
123 % o sigma: the sigma of the gaussian smoothing filter.
124 %
125 % o lower_percent: percentage of edge pixels in the lower threshold.
126 %
127 % o upper_percent: percentage of edge pixels in the upper threshold.
128 %
129 % o exception: return any errors or warnings in this structure.
130 %
131 */
132
133 typedef struct _CannyInfo
134 {
135 double
136 magnitude,
137 intensity;
138
139 int
140 orientation;
141
142 ssize_t
143 x,
144 y;
145 } CannyInfo;
146
IsAuthenticPixel(const Image * image,const ssize_t x,const ssize_t y)147 static inline MagickBooleanType IsAuthenticPixel(const Image *image,
148 const ssize_t x,const ssize_t y)
149 {
150 if ((x < 0) || (x >= (ssize_t) image->columns))
151 return(MagickFalse);
152 if ((y < 0) || (y >= (ssize_t) image->rows))
153 return(MagickFalse);
154 return(MagickTrue);
155 }
156
TraceEdges(Image * edge_image,CacheView * edge_view,MatrixInfo * canny_cache,const ssize_t x,const ssize_t y,const double lower_threshold,ExceptionInfo * exception)157 static MagickBooleanType TraceEdges(Image *edge_image,CacheView *edge_view,
158 MatrixInfo *canny_cache,const ssize_t x,const ssize_t y,
159 const double lower_threshold,ExceptionInfo *exception)
160 {
161 CannyInfo
162 edge,
163 pixel;
164
165 MagickBooleanType
166 status;
167
168 register Quantum
169 *q;
170
171 register ssize_t
172 i;
173
174 q=GetCacheViewAuthenticPixels(edge_view,x,y,1,1,exception);
175 if (q == (Quantum *) NULL)
176 return(MagickFalse);
177 *q=QuantumRange;
178 status=SyncCacheViewAuthenticPixels(edge_view,exception);
179 if (status == MagickFalse)
180 return(MagickFalse);
181 if (GetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
182 return(MagickFalse);
183 edge.x=x;
184 edge.y=y;
185 if (SetMatrixElement(canny_cache,0,0,&edge) == MagickFalse)
186 return(MagickFalse);
187 for (i=1; i != 0; )
188 {
189 ssize_t
190 v;
191
192 i--;
193 status=GetMatrixElement(canny_cache,i,0,&edge);
194 if (status == MagickFalse)
195 return(MagickFalse);
196 for (v=(-1); v <= 1; v++)
197 {
198 ssize_t
199 u;
200
201 for (u=(-1); u <= 1; u++)
202 {
203 if ((u == 0) && (v == 0))
204 continue;
205 if (IsAuthenticPixel(edge_image,edge.x+u,edge.y+v) == MagickFalse)
206 continue;
207 /*
208 Not an edge if gradient value is below the lower threshold.
209 */
210 q=GetCacheViewAuthenticPixels(edge_view,edge.x+u,edge.y+v,1,1,
211 exception);
212 if (q == (Quantum *) NULL)
213 return(MagickFalse);
214 status=GetMatrixElement(canny_cache,edge.x+u,edge.y+v,&pixel);
215 if (status == MagickFalse)
216 return(MagickFalse);
217 if ((GetPixelIntensity(edge_image,q) == 0.0) &&
218 (pixel.intensity >= lower_threshold))
219 {
220 *q=QuantumRange;
221 status=SyncCacheViewAuthenticPixels(edge_view,exception);
222 if (status == MagickFalse)
223 return(MagickFalse);
224 edge.x+=u;
225 edge.y+=v;
226 status=SetMatrixElement(canny_cache,i,0,&edge);
227 if (status == MagickFalse)
228 return(MagickFalse);
229 i++;
230 }
231 }
232 }
233 }
234 return(MagickTrue);
235 }
236
CannyEdgeImage(const Image * image,const double radius,const double sigma,const double lower_percent,const double upper_percent,ExceptionInfo * exception)237 MagickExport Image *CannyEdgeImage(const Image *image,const double radius,
238 const double sigma,const double lower_percent,const double upper_percent,
239 ExceptionInfo *exception)
240 {
241 #define CannyEdgeImageTag "CannyEdge/Image"
242
243 CacheView
244 *edge_view;
245
246 CannyInfo
247 element;
248
249 char
250 geometry[MagickPathExtent];
251
252 double
253 lower_threshold,
254 max,
255 min,
256 upper_threshold;
257
258 Image
259 *edge_image;
260
261 KernelInfo
262 *kernel_info;
263
264 MagickBooleanType
265 status;
266
267 MagickOffsetType
268 progress;
269
270 MatrixInfo
271 *canny_cache;
272
273 ssize_t
274 y;
275
276 assert(image != (const Image *) NULL);
277 assert(image->signature == MagickCoreSignature);
278 if (image->debug != MagickFalse)
279 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
280 assert(exception != (ExceptionInfo *) NULL);
281 assert(exception->signature == MagickCoreSignature);
282 /*
283 Filter out noise.
284 */
285 (void) FormatLocaleString(geometry,MagickPathExtent,
286 "blur:%.20gx%.20g;blur:%.20gx%.20g+90",radius,sigma,radius,sigma);
287 kernel_info=AcquireKernelInfo(geometry,exception);
288 if (kernel_info == (KernelInfo *) NULL)
289 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
290 edge_image=MorphologyImage(image,ConvolveMorphology,1,kernel_info,exception);
291 kernel_info=DestroyKernelInfo(kernel_info);
292 if (edge_image == (Image *) NULL)
293 return((Image *) NULL);
294 if (TransformImageColorspace(edge_image,GRAYColorspace,exception) == MagickFalse)
295 {
296 edge_image=DestroyImage(edge_image);
297 return((Image *) NULL);
298 }
299 (void) SetImageAlphaChannel(edge_image,OffAlphaChannel,exception);
300 /*
301 Find the intensity gradient of the image.
302 */
303 canny_cache=AcquireMatrixInfo(edge_image->columns,edge_image->rows,
304 sizeof(CannyInfo),exception);
305 if (canny_cache == (MatrixInfo *) NULL)
306 {
307 edge_image=DestroyImage(edge_image);
308 return((Image *) NULL);
309 }
310 status=MagickTrue;
311 edge_view=AcquireVirtualCacheView(edge_image,exception);
312 #if defined(MAGICKCORE_OPENMP_SUPPORT)
313 #pragma omp parallel for schedule(static) shared(status) \
314 magick_number_threads(edge_image,edge_image,edge_image->rows,1)
315 #endif
316 for (y=0; y < (ssize_t) edge_image->rows; y++)
317 {
318 register const Quantum
319 *magick_restrict p;
320
321 register ssize_t
322 x;
323
324 if (status == MagickFalse)
325 continue;
326 p=GetCacheViewVirtualPixels(edge_view,0,y,edge_image->columns+1,2,
327 exception);
328 if (p == (const Quantum *) NULL)
329 {
330 status=MagickFalse;
331 continue;
332 }
333 for (x=0; x < (ssize_t) edge_image->columns; x++)
334 {
335 CannyInfo
336 pixel;
337
338 double
339 dx,
340 dy;
341
342 register const Quantum
343 *magick_restrict kernel_pixels;
344
345 ssize_t
346 v;
347
348 static double
349 Gx[2][2] =
350 {
351 { -1.0, +1.0 },
352 { -1.0, +1.0 }
353 },
354 Gy[2][2] =
355 {
356 { +1.0, +1.0 },
357 { -1.0, -1.0 }
358 };
359
360 (void) memset(&pixel,0,sizeof(pixel));
361 dx=0.0;
362 dy=0.0;
363 kernel_pixels=p;
364 for (v=0; v < 2; v++)
365 {
366 ssize_t
367 u;
368
369 for (u=0; u < 2; u++)
370 {
371 double
372 intensity;
373
374 intensity=GetPixelIntensity(edge_image,kernel_pixels+u);
375 dx+=0.5*Gx[v][u]*intensity;
376 dy+=0.5*Gy[v][u]*intensity;
377 }
378 kernel_pixels+=edge_image->columns+1;
379 }
380 pixel.magnitude=hypot(dx,dy);
381 pixel.orientation=0;
382 if (fabs(dx) > MagickEpsilon)
383 {
384 double
385 slope;
386
387 slope=dy/dx;
388 if (slope < 0.0)
389 {
390 if (slope < -2.41421356237)
391 pixel.orientation=0;
392 else
393 if (slope < -0.414213562373)
394 pixel.orientation=1;
395 else
396 pixel.orientation=2;
397 }
398 else
399 {
400 if (slope > 2.41421356237)
401 pixel.orientation=0;
402 else
403 if (slope > 0.414213562373)
404 pixel.orientation=3;
405 else
406 pixel.orientation=2;
407 }
408 }
409 if (SetMatrixElement(canny_cache,x,y,&pixel) == MagickFalse)
410 continue;
411 p+=GetPixelChannels(edge_image);
412 }
413 }
414 edge_view=DestroyCacheView(edge_view);
415 /*
416 Non-maxima suppression, remove pixels that are not considered to be part
417 of an edge.
418 */
419 progress=0;
420 (void) GetMatrixElement(canny_cache,0,0,&element);
421 max=element.intensity;
422 min=element.intensity;
423 edge_view=AcquireAuthenticCacheView(edge_image,exception);
424 #if defined(MAGICKCORE_OPENMP_SUPPORT)
425 #pragma omp parallel for schedule(static) shared(status) \
426 magick_number_threads(edge_image,edge_image,edge_image->rows,1)
427 #endif
428 for (y=0; y < (ssize_t) edge_image->rows; y++)
429 {
430 register Quantum
431 *magick_restrict q;
432
433 register ssize_t
434 x;
435
436 if (status == MagickFalse)
437 continue;
438 q=GetCacheViewAuthenticPixels(edge_view,0,y,edge_image->columns,1,
439 exception);
440 if (q == (Quantum *) NULL)
441 {
442 status=MagickFalse;
443 continue;
444 }
445 for (x=0; x < (ssize_t) edge_image->columns; x++)
446 {
447 CannyInfo
448 alpha_pixel,
449 beta_pixel,
450 pixel;
451
452 (void) GetMatrixElement(canny_cache,x,y,&pixel);
453 switch (pixel.orientation)
454 {
455 case 0:
456 default:
457 {
458 /*
459 0 degrees, north and south.
460 */
461 (void) GetMatrixElement(canny_cache,x,y-1,&alpha_pixel);
462 (void) GetMatrixElement(canny_cache,x,y+1,&beta_pixel);
463 break;
464 }
465 case 1:
466 {
467 /*
468 45 degrees, northwest and southeast.
469 */
470 (void) GetMatrixElement(canny_cache,x-1,y-1,&alpha_pixel);
471 (void) GetMatrixElement(canny_cache,x+1,y+1,&beta_pixel);
472 break;
473 }
474 case 2:
475 {
476 /*
477 90 degrees, east and west.
478 */
479 (void) GetMatrixElement(canny_cache,x-1,y,&alpha_pixel);
480 (void) GetMatrixElement(canny_cache,x+1,y,&beta_pixel);
481 break;
482 }
483 case 3:
484 {
485 /*
486 135 degrees, northeast and southwest.
487 */
488 (void) GetMatrixElement(canny_cache,x+1,y-1,&beta_pixel);
489 (void) GetMatrixElement(canny_cache,x-1,y+1,&alpha_pixel);
490 break;
491 }
492 }
493 pixel.intensity=pixel.magnitude;
494 if ((pixel.magnitude < alpha_pixel.magnitude) ||
495 (pixel.magnitude < beta_pixel.magnitude))
496 pixel.intensity=0;
497 (void) SetMatrixElement(canny_cache,x,y,&pixel);
498 #if defined(MAGICKCORE_OPENMP_SUPPORT)
499 #pragma omp critical (MagickCore_CannyEdgeImage)
500 #endif
501 {
502 if (pixel.intensity < min)
503 min=pixel.intensity;
504 if (pixel.intensity > max)
505 max=pixel.intensity;
506 }
507 *q=0;
508 q+=GetPixelChannels(edge_image);
509 }
510 if (SyncCacheViewAuthenticPixels(edge_view,exception) == MagickFalse)
511 status=MagickFalse;
512 }
513 edge_view=DestroyCacheView(edge_view);
514 /*
515 Estimate hysteresis threshold.
516 */
517 lower_threshold=lower_percent*(max-min)+min;
518 upper_threshold=upper_percent*(max-min)+min;
519 /*
520 Hysteresis threshold.
521 */
522 edge_view=AcquireAuthenticCacheView(edge_image,exception);
523 for (y=0; y < (ssize_t) edge_image->rows; y++)
524 {
525 register ssize_t
526 x;
527
528 if (status == MagickFalse)
529 continue;
530 for (x=0; x < (ssize_t) edge_image->columns; x++)
531 {
532 CannyInfo
533 pixel;
534
535 register const Quantum
536 *magick_restrict p;
537
538 /*
539 Edge if pixel gradient higher than upper threshold.
540 */
541 p=GetCacheViewVirtualPixels(edge_view,x,y,1,1,exception);
542 if (p == (const Quantum *) NULL)
543 continue;
544 status=GetMatrixElement(canny_cache,x,y,&pixel);
545 if (status == MagickFalse)
546 continue;
547 if ((GetPixelIntensity(edge_image,p) == 0.0) &&
548 (pixel.intensity >= upper_threshold))
549 status=TraceEdges(edge_image,edge_view,canny_cache,x,y,lower_threshold,
550 exception);
551 }
552 if (image->progress_monitor != (MagickProgressMonitor) NULL)
553 {
554 MagickBooleanType
555 proceed;
556
557 #if defined(MAGICKCORE_OPENMP_SUPPORT)
558 #pragma omp atomic
559 #endif
560 progress++;
561 proceed=SetImageProgress(image,CannyEdgeImageTag,progress,image->rows);
562 if (proceed == MagickFalse)
563 status=MagickFalse;
564 }
565 }
566 edge_view=DestroyCacheView(edge_view);
567 /*
568 Free resources.
569 */
570 canny_cache=DestroyMatrixInfo(canny_cache);
571 return(edge_image);
572 }
573
574 /*
575 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
576 % %
577 % %
578 % %
579 % G e t I m a g e F e a t u r e s %
580 % %
581 % %
582 % %
583 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
584 %
585 % GetImageFeatures() returns features for each channel in the image in
586 % each of four directions (horizontal, vertical, left and right diagonals)
587 % for the specified distance. The features include the angular second
588 % moment, contrast, correlation, sum of squares: variance, inverse difference
589 % moment, sum average, sum varience, sum entropy, entropy, difference variance,% difference entropy, information measures of correlation 1, information
590 % measures of correlation 2, and maximum correlation coefficient. You can
591 % access the red channel contrast, for example, like this:
592 %
593 % channel_features=GetImageFeatures(image,1,exception);
594 % contrast=channel_features[RedPixelChannel].contrast[0];
595 %
596 % Use MagickRelinquishMemory() to free the features buffer.
597 %
598 % The format of the GetImageFeatures method is:
599 %
600 % ChannelFeatures *GetImageFeatures(const Image *image,
601 % const size_t distance,ExceptionInfo *exception)
602 %
603 % A description of each parameter follows:
604 %
605 % o image: the image.
606 %
607 % o distance: the distance.
608 %
609 % o exception: return any errors or warnings in this structure.
610 %
611 */
612
MagickLog10(const double x)613 static inline double MagickLog10(const double x)
614 {
615 #define Log10Epsilon (1.0e-11)
616
617 if (fabs(x) < Log10Epsilon)
618 return(log10(Log10Epsilon));
619 return(log10(fabs(x)));
620 }
621
GetImageFeatures(const Image * image,const size_t distance,ExceptionInfo * exception)622 MagickExport ChannelFeatures *GetImageFeatures(const Image *image,
623 const size_t distance,ExceptionInfo *exception)
624 {
625 typedef struct _ChannelStatistics
626 {
627 PixelInfo
628 direction[4]; /* horizontal, vertical, left and right diagonals */
629 } ChannelStatistics;
630
631 CacheView
632 *image_view;
633
634 ChannelFeatures
635 *channel_features;
636
637 ChannelStatistics
638 **cooccurrence,
639 correlation,
640 *density_x,
641 *density_xy,
642 *density_y,
643 entropy_x,
644 entropy_xy,
645 entropy_xy1,
646 entropy_xy2,
647 entropy_y,
648 mean,
649 **Q,
650 *sum,
651 sum_squares,
652 variance;
653
654 PixelPacket
655 gray,
656 *grays;
657
658 MagickBooleanType
659 status;
660
661 register ssize_t
662 i,
663 r;
664
665 size_t
666 length;
667
668 unsigned int
669 number_grays;
670
671 assert(image != (Image *) NULL);
672 assert(image->signature == MagickCoreSignature);
673 if (image->debug != MagickFalse)
674 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
675 if ((image->columns < (distance+1)) || (image->rows < (distance+1)))
676 return((ChannelFeatures *) NULL);
677 length=MaxPixelChannels+1UL;
678 channel_features=(ChannelFeatures *) AcquireQuantumMemory(length,
679 sizeof(*channel_features));
680 if (channel_features == (ChannelFeatures *) NULL)
681 ThrowFatalException(ResourceLimitFatalError,"MemoryAllocationFailed");
682 (void) memset(channel_features,0,length*
683 sizeof(*channel_features));
684 /*
685 Form grays.
686 */
687 grays=(PixelPacket *) AcquireQuantumMemory(MaxMap+1UL,sizeof(*grays));
688 if (grays == (PixelPacket *) NULL)
689 {
690 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
691 channel_features);
692 (void) ThrowMagickException(exception,GetMagickModule(),
693 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
694 return(channel_features);
695 }
696 for (i=0; i <= (ssize_t) MaxMap; i++)
697 {
698 grays[i].red=(~0U);
699 grays[i].green=(~0U);
700 grays[i].blue=(~0U);
701 grays[i].alpha=(~0U);
702 grays[i].black=(~0U);
703 }
704 status=MagickTrue;
705 image_view=AcquireVirtualCacheView(image,exception);
706 #if defined(MAGICKCORE_OPENMP_SUPPORT)
707 #pragma omp parallel for schedule(static) shared(status) \
708 magick_number_threads(image,image,image->rows,1)
709 #endif
710 for (r=0; r < (ssize_t) image->rows; r++)
711 {
712 register const Quantum
713 *magick_restrict p;
714
715 register ssize_t
716 x;
717
718 if (status == MagickFalse)
719 continue;
720 p=GetCacheViewVirtualPixels(image_view,0,r,image->columns,1,exception);
721 if (p == (const Quantum *) NULL)
722 {
723 status=MagickFalse;
724 continue;
725 }
726 for (x=0; x < (ssize_t) image->columns; x++)
727 {
728 grays[ScaleQuantumToMap(GetPixelRed(image,p))].red=
729 ScaleQuantumToMap(GetPixelRed(image,p));
730 grays[ScaleQuantumToMap(GetPixelGreen(image,p))].green=
731 ScaleQuantumToMap(GetPixelGreen(image,p));
732 grays[ScaleQuantumToMap(GetPixelBlue(image,p))].blue=
733 ScaleQuantumToMap(GetPixelBlue(image,p));
734 if (image->colorspace == CMYKColorspace)
735 grays[ScaleQuantumToMap(GetPixelBlack(image,p))].black=
736 ScaleQuantumToMap(GetPixelBlack(image,p));
737 if (image->alpha_trait != UndefinedPixelTrait)
738 grays[ScaleQuantumToMap(GetPixelAlpha(image,p))].alpha=
739 ScaleQuantumToMap(GetPixelAlpha(image,p));
740 p+=GetPixelChannels(image);
741 }
742 }
743 image_view=DestroyCacheView(image_view);
744 if (status == MagickFalse)
745 {
746 grays=(PixelPacket *) RelinquishMagickMemory(grays);
747 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
748 channel_features);
749 return(channel_features);
750 }
751 (void) memset(&gray,0,sizeof(gray));
752 for (i=0; i <= (ssize_t) MaxMap; i++)
753 {
754 if (grays[i].red != ~0U)
755 grays[gray.red++].red=grays[i].red;
756 if (grays[i].green != ~0U)
757 grays[gray.green++].green=grays[i].green;
758 if (grays[i].blue != ~0U)
759 grays[gray.blue++].blue=grays[i].blue;
760 if (image->colorspace == CMYKColorspace)
761 if (grays[i].black != ~0U)
762 grays[gray.black++].black=grays[i].black;
763 if (image->alpha_trait != UndefinedPixelTrait)
764 if (grays[i].alpha != ~0U)
765 grays[gray.alpha++].alpha=grays[i].alpha;
766 }
767 /*
768 Allocate spatial dependence matrix.
769 */
770 number_grays=gray.red;
771 if (gray.green > number_grays)
772 number_grays=gray.green;
773 if (gray.blue > number_grays)
774 number_grays=gray.blue;
775 if (image->colorspace == CMYKColorspace)
776 if (gray.black > number_grays)
777 number_grays=gray.black;
778 if (image->alpha_trait != UndefinedPixelTrait)
779 if (gray.alpha > number_grays)
780 number_grays=gray.alpha;
781 cooccurrence=(ChannelStatistics **) AcquireQuantumMemory(number_grays,
782 sizeof(*cooccurrence));
783 density_x=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
784 sizeof(*density_x));
785 density_xy=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
786 sizeof(*density_xy));
787 density_y=(ChannelStatistics *) AcquireQuantumMemory(2*(number_grays+1),
788 sizeof(*density_y));
789 Q=(ChannelStatistics **) AcquireQuantumMemory(number_grays,sizeof(*Q));
790 sum=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(*sum));
791 if ((cooccurrence == (ChannelStatistics **) NULL) ||
792 (density_x == (ChannelStatistics *) NULL) ||
793 (density_xy == (ChannelStatistics *) NULL) ||
794 (density_y == (ChannelStatistics *) NULL) ||
795 (Q == (ChannelStatistics **) NULL) ||
796 (sum == (ChannelStatistics *) NULL))
797 {
798 if (Q != (ChannelStatistics **) NULL)
799 {
800 for (i=0; i < (ssize_t) number_grays; i++)
801 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
802 Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
803 }
804 if (sum != (ChannelStatistics *) NULL)
805 sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
806 if (density_y != (ChannelStatistics *) NULL)
807 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
808 if (density_xy != (ChannelStatistics *) NULL)
809 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
810 if (density_x != (ChannelStatistics *) NULL)
811 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
812 if (cooccurrence != (ChannelStatistics **) NULL)
813 {
814 for (i=0; i < (ssize_t) number_grays; i++)
815 cooccurrence[i]=(ChannelStatistics *)
816 RelinquishMagickMemory(cooccurrence[i]);
817 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(
818 cooccurrence);
819 }
820 grays=(PixelPacket *) RelinquishMagickMemory(grays);
821 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
822 channel_features);
823 (void) ThrowMagickException(exception,GetMagickModule(),
824 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
825 return(channel_features);
826 }
827 (void) memset(&correlation,0,sizeof(correlation));
828 (void) memset(density_x,0,2*(number_grays+1)*sizeof(*density_x));
829 (void) memset(density_xy,0,2*(number_grays+1)*sizeof(*density_xy));
830 (void) memset(density_y,0,2*(number_grays+1)*sizeof(*density_y));
831 (void) memset(&mean,0,sizeof(mean));
832 (void) memset(sum,0,number_grays*sizeof(*sum));
833 (void) memset(&sum_squares,0,sizeof(sum_squares));
834 (void) memset(density_xy,0,2*number_grays*sizeof(*density_xy));
835 (void) memset(&entropy_x,0,sizeof(entropy_x));
836 (void) memset(&entropy_xy,0,sizeof(entropy_xy));
837 (void) memset(&entropy_xy1,0,sizeof(entropy_xy1));
838 (void) memset(&entropy_xy2,0,sizeof(entropy_xy2));
839 (void) memset(&entropy_y,0,sizeof(entropy_y));
840 (void) memset(&variance,0,sizeof(variance));
841 for (i=0; i < (ssize_t) number_grays; i++)
842 {
843 cooccurrence[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,
844 sizeof(**cooccurrence));
845 Q[i]=(ChannelStatistics *) AcquireQuantumMemory(number_grays,sizeof(**Q));
846 if ((cooccurrence[i] == (ChannelStatistics *) NULL) ||
847 (Q[i] == (ChannelStatistics *) NULL))
848 break;
849 (void) memset(cooccurrence[i],0,number_grays*
850 sizeof(**cooccurrence));
851 (void) memset(Q[i],0,number_grays*sizeof(**Q));
852 }
853 if (i < (ssize_t) number_grays)
854 {
855 for (i--; i >= 0; i--)
856 {
857 if (Q[i] != (ChannelStatistics *) NULL)
858 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
859 if (cooccurrence[i] != (ChannelStatistics *) NULL)
860 cooccurrence[i]=(ChannelStatistics *)
861 RelinquishMagickMemory(cooccurrence[i]);
862 }
863 Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
864 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
865 sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
866 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
867 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
868 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
869 grays=(PixelPacket *) RelinquishMagickMemory(grays);
870 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
871 channel_features);
872 (void) ThrowMagickException(exception,GetMagickModule(),
873 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
874 return(channel_features);
875 }
876 /*
877 Initialize spatial dependence matrix.
878 */
879 status=MagickTrue;
880 image_view=AcquireVirtualCacheView(image,exception);
881 for (r=0; r < (ssize_t) image->rows; r++)
882 {
883 register const Quantum
884 *magick_restrict p;
885
886 register ssize_t
887 x;
888
889 ssize_t
890 offset,
891 u,
892 v;
893
894 if (status == MagickFalse)
895 continue;
896 p=GetCacheViewVirtualPixels(image_view,-(ssize_t) distance,r,image->columns+
897 2*distance,distance+2,exception);
898 if (p == (const Quantum *) NULL)
899 {
900 status=MagickFalse;
901 continue;
902 }
903 p+=distance*GetPixelChannels(image);;
904 for (x=0; x < (ssize_t) image->columns; x++)
905 {
906 for (i=0; i < 4; i++)
907 {
908 switch (i)
909 {
910 case 0:
911 default:
912 {
913 /*
914 Horizontal adjacency.
915 */
916 offset=(ssize_t) distance;
917 break;
918 }
919 case 1:
920 {
921 /*
922 Vertical adjacency.
923 */
924 offset=(ssize_t) (image->columns+2*distance);
925 break;
926 }
927 case 2:
928 {
929 /*
930 Right diagonal adjacency.
931 */
932 offset=(ssize_t) ((image->columns+2*distance)-distance);
933 break;
934 }
935 case 3:
936 {
937 /*
938 Left diagonal adjacency.
939 */
940 offset=(ssize_t) ((image->columns+2*distance)+distance);
941 break;
942 }
943 }
944 u=0;
945 v=0;
946 while (grays[u].red != ScaleQuantumToMap(GetPixelRed(image,p)))
947 u++;
948 while (grays[v].red != ScaleQuantumToMap(GetPixelRed(image,p+offset*GetPixelChannels(image))))
949 v++;
950 cooccurrence[u][v].direction[i].red++;
951 cooccurrence[v][u].direction[i].red++;
952 u=0;
953 v=0;
954 while (grays[u].green != ScaleQuantumToMap(GetPixelGreen(image,p)))
955 u++;
956 while (grays[v].green != ScaleQuantumToMap(GetPixelGreen(image,p+offset*GetPixelChannels(image))))
957 v++;
958 cooccurrence[u][v].direction[i].green++;
959 cooccurrence[v][u].direction[i].green++;
960 u=0;
961 v=0;
962 while (grays[u].blue != ScaleQuantumToMap(GetPixelBlue(image,p)))
963 u++;
964 while (grays[v].blue != ScaleQuantumToMap(GetPixelBlue(image,p+offset*GetPixelChannels(image))))
965 v++;
966 cooccurrence[u][v].direction[i].blue++;
967 cooccurrence[v][u].direction[i].blue++;
968 if (image->colorspace == CMYKColorspace)
969 {
970 u=0;
971 v=0;
972 while (grays[u].black != ScaleQuantumToMap(GetPixelBlack(image,p)))
973 u++;
974 while (grays[v].black != ScaleQuantumToMap(GetPixelBlack(image,p+offset*GetPixelChannels(image))))
975 v++;
976 cooccurrence[u][v].direction[i].black++;
977 cooccurrence[v][u].direction[i].black++;
978 }
979 if (image->alpha_trait != UndefinedPixelTrait)
980 {
981 u=0;
982 v=0;
983 while (grays[u].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p)))
984 u++;
985 while (grays[v].alpha != ScaleQuantumToMap(GetPixelAlpha(image,p+offset*GetPixelChannels(image))))
986 v++;
987 cooccurrence[u][v].direction[i].alpha++;
988 cooccurrence[v][u].direction[i].alpha++;
989 }
990 }
991 p+=GetPixelChannels(image);
992 }
993 }
994 grays=(PixelPacket *) RelinquishMagickMemory(grays);
995 image_view=DestroyCacheView(image_view);
996 if (status == MagickFalse)
997 {
998 for (i=0; i < (ssize_t) number_grays; i++)
999 cooccurrence[i]=(ChannelStatistics *)
1000 RelinquishMagickMemory(cooccurrence[i]);
1001 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1002 channel_features=(ChannelFeatures *) RelinquishMagickMemory(
1003 channel_features);
1004 (void) ThrowMagickException(exception,GetMagickModule(),
1005 ResourceLimitError,"MemoryAllocationFailed","`%s'",image->filename);
1006 return(channel_features);
1007 }
1008 /*
1009 Normalize spatial dependence matrix.
1010 */
1011 for (i=0; i < 4; i++)
1012 {
1013 double
1014 normalize;
1015
1016 register ssize_t
1017 y;
1018
1019 switch (i)
1020 {
1021 case 0:
1022 default:
1023 {
1024 /*
1025 Horizontal adjacency.
1026 */
1027 normalize=2.0*image->rows*(image->columns-distance);
1028 break;
1029 }
1030 case 1:
1031 {
1032 /*
1033 Vertical adjacency.
1034 */
1035 normalize=2.0*(image->rows-distance)*image->columns;
1036 break;
1037 }
1038 case 2:
1039 {
1040 /*
1041 Right diagonal adjacency.
1042 */
1043 normalize=2.0*(image->rows-distance)*(image->columns-distance);
1044 break;
1045 }
1046 case 3:
1047 {
1048 /*
1049 Left diagonal adjacency.
1050 */
1051 normalize=2.0*(image->rows-distance)*(image->columns-distance);
1052 break;
1053 }
1054 }
1055 normalize=PerceptibleReciprocal(normalize);
1056 for (y=0; y < (ssize_t) number_grays; y++)
1057 {
1058 register ssize_t
1059 x;
1060
1061 for (x=0; x < (ssize_t) number_grays; x++)
1062 {
1063 cooccurrence[x][y].direction[i].red*=normalize;
1064 cooccurrence[x][y].direction[i].green*=normalize;
1065 cooccurrence[x][y].direction[i].blue*=normalize;
1066 if (image->colorspace == CMYKColorspace)
1067 cooccurrence[x][y].direction[i].black*=normalize;
1068 if (image->alpha_trait != UndefinedPixelTrait)
1069 cooccurrence[x][y].direction[i].alpha*=normalize;
1070 }
1071 }
1072 }
1073 /*
1074 Compute texture features.
1075 */
1076 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1077 #pragma omp parallel for schedule(static) shared(status) \
1078 magick_number_threads(image,image,number_grays,1)
1079 #endif
1080 for (i=0; i < 4; i++)
1081 {
1082 register ssize_t
1083 y;
1084
1085 for (y=0; y < (ssize_t) number_grays; y++)
1086 {
1087 register ssize_t
1088 x;
1089
1090 for (x=0; x < (ssize_t) number_grays; x++)
1091 {
1092 /*
1093 Angular second moment: measure of homogeneity of the image.
1094 */
1095 channel_features[RedPixelChannel].angular_second_moment[i]+=
1096 cooccurrence[x][y].direction[i].red*
1097 cooccurrence[x][y].direction[i].red;
1098 channel_features[GreenPixelChannel].angular_second_moment[i]+=
1099 cooccurrence[x][y].direction[i].green*
1100 cooccurrence[x][y].direction[i].green;
1101 channel_features[BluePixelChannel].angular_second_moment[i]+=
1102 cooccurrence[x][y].direction[i].blue*
1103 cooccurrence[x][y].direction[i].blue;
1104 if (image->colorspace == CMYKColorspace)
1105 channel_features[BlackPixelChannel].angular_second_moment[i]+=
1106 cooccurrence[x][y].direction[i].black*
1107 cooccurrence[x][y].direction[i].black;
1108 if (image->alpha_trait != UndefinedPixelTrait)
1109 channel_features[AlphaPixelChannel].angular_second_moment[i]+=
1110 cooccurrence[x][y].direction[i].alpha*
1111 cooccurrence[x][y].direction[i].alpha;
1112 /*
1113 Correlation: measure of linear-dependencies in the image.
1114 */
1115 sum[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1116 sum[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1117 sum[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1118 if (image->colorspace == CMYKColorspace)
1119 sum[y].direction[i].black+=cooccurrence[x][y].direction[i].black;
1120 if (image->alpha_trait != UndefinedPixelTrait)
1121 sum[y].direction[i].alpha+=cooccurrence[x][y].direction[i].alpha;
1122 correlation.direction[i].red+=x*y*cooccurrence[x][y].direction[i].red;
1123 correlation.direction[i].green+=x*y*
1124 cooccurrence[x][y].direction[i].green;
1125 correlation.direction[i].blue+=x*y*
1126 cooccurrence[x][y].direction[i].blue;
1127 if (image->colorspace == CMYKColorspace)
1128 correlation.direction[i].black+=x*y*
1129 cooccurrence[x][y].direction[i].black;
1130 if (image->alpha_trait != UndefinedPixelTrait)
1131 correlation.direction[i].alpha+=x*y*
1132 cooccurrence[x][y].direction[i].alpha;
1133 /*
1134 Inverse Difference Moment.
1135 */
1136 channel_features[RedPixelChannel].inverse_difference_moment[i]+=
1137 cooccurrence[x][y].direction[i].red/((y-x)*(y-x)+1);
1138 channel_features[GreenPixelChannel].inverse_difference_moment[i]+=
1139 cooccurrence[x][y].direction[i].green/((y-x)*(y-x)+1);
1140 channel_features[BluePixelChannel].inverse_difference_moment[i]+=
1141 cooccurrence[x][y].direction[i].blue/((y-x)*(y-x)+1);
1142 if (image->colorspace == CMYKColorspace)
1143 channel_features[BlackPixelChannel].inverse_difference_moment[i]+=
1144 cooccurrence[x][y].direction[i].black/((y-x)*(y-x)+1);
1145 if (image->alpha_trait != UndefinedPixelTrait)
1146 channel_features[AlphaPixelChannel].inverse_difference_moment[i]+=
1147 cooccurrence[x][y].direction[i].alpha/((y-x)*(y-x)+1);
1148 /*
1149 Sum average.
1150 */
1151 density_xy[y+x+2].direction[i].red+=
1152 cooccurrence[x][y].direction[i].red;
1153 density_xy[y+x+2].direction[i].green+=
1154 cooccurrence[x][y].direction[i].green;
1155 density_xy[y+x+2].direction[i].blue+=
1156 cooccurrence[x][y].direction[i].blue;
1157 if (image->colorspace == CMYKColorspace)
1158 density_xy[y+x+2].direction[i].black+=
1159 cooccurrence[x][y].direction[i].black;
1160 if (image->alpha_trait != UndefinedPixelTrait)
1161 density_xy[y+x+2].direction[i].alpha+=
1162 cooccurrence[x][y].direction[i].alpha;
1163 /*
1164 Entropy.
1165 */
1166 channel_features[RedPixelChannel].entropy[i]-=
1167 cooccurrence[x][y].direction[i].red*
1168 MagickLog10(cooccurrence[x][y].direction[i].red);
1169 channel_features[GreenPixelChannel].entropy[i]-=
1170 cooccurrence[x][y].direction[i].green*
1171 MagickLog10(cooccurrence[x][y].direction[i].green);
1172 channel_features[BluePixelChannel].entropy[i]-=
1173 cooccurrence[x][y].direction[i].blue*
1174 MagickLog10(cooccurrence[x][y].direction[i].blue);
1175 if (image->colorspace == CMYKColorspace)
1176 channel_features[BlackPixelChannel].entropy[i]-=
1177 cooccurrence[x][y].direction[i].black*
1178 MagickLog10(cooccurrence[x][y].direction[i].black);
1179 if (image->alpha_trait != UndefinedPixelTrait)
1180 channel_features[AlphaPixelChannel].entropy[i]-=
1181 cooccurrence[x][y].direction[i].alpha*
1182 MagickLog10(cooccurrence[x][y].direction[i].alpha);
1183 /*
1184 Information Measures of Correlation.
1185 */
1186 density_x[x].direction[i].red+=cooccurrence[x][y].direction[i].red;
1187 density_x[x].direction[i].green+=cooccurrence[x][y].direction[i].green;
1188 density_x[x].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1189 if (image->alpha_trait != UndefinedPixelTrait)
1190 density_x[x].direction[i].alpha+=
1191 cooccurrence[x][y].direction[i].alpha;
1192 if (image->colorspace == CMYKColorspace)
1193 density_x[x].direction[i].black+=
1194 cooccurrence[x][y].direction[i].black;
1195 density_y[y].direction[i].red+=cooccurrence[x][y].direction[i].red;
1196 density_y[y].direction[i].green+=cooccurrence[x][y].direction[i].green;
1197 density_y[y].direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1198 if (image->colorspace == CMYKColorspace)
1199 density_y[y].direction[i].black+=
1200 cooccurrence[x][y].direction[i].black;
1201 if (image->alpha_trait != UndefinedPixelTrait)
1202 density_y[y].direction[i].alpha+=
1203 cooccurrence[x][y].direction[i].alpha;
1204 }
1205 mean.direction[i].red+=y*sum[y].direction[i].red;
1206 sum_squares.direction[i].red+=y*y*sum[y].direction[i].red;
1207 mean.direction[i].green+=y*sum[y].direction[i].green;
1208 sum_squares.direction[i].green+=y*y*sum[y].direction[i].green;
1209 mean.direction[i].blue+=y*sum[y].direction[i].blue;
1210 sum_squares.direction[i].blue+=y*y*sum[y].direction[i].blue;
1211 if (image->colorspace == CMYKColorspace)
1212 {
1213 mean.direction[i].black+=y*sum[y].direction[i].black;
1214 sum_squares.direction[i].black+=y*y*sum[y].direction[i].black;
1215 }
1216 if (image->alpha_trait != UndefinedPixelTrait)
1217 {
1218 mean.direction[i].alpha+=y*sum[y].direction[i].alpha;
1219 sum_squares.direction[i].alpha+=y*y*sum[y].direction[i].alpha;
1220 }
1221 }
1222 /*
1223 Correlation: measure of linear-dependencies in the image.
1224 */
1225 channel_features[RedPixelChannel].correlation[i]=
1226 (correlation.direction[i].red-mean.direction[i].red*
1227 mean.direction[i].red)/(sqrt(sum_squares.direction[i].red-
1228 (mean.direction[i].red*mean.direction[i].red))*sqrt(
1229 sum_squares.direction[i].red-(mean.direction[i].red*
1230 mean.direction[i].red)));
1231 channel_features[GreenPixelChannel].correlation[i]=
1232 (correlation.direction[i].green-mean.direction[i].green*
1233 mean.direction[i].green)/(sqrt(sum_squares.direction[i].green-
1234 (mean.direction[i].green*mean.direction[i].green))*sqrt(
1235 sum_squares.direction[i].green-(mean.direction[i].green*
1236 mean.direction[i].green)));
1237 channel_features[BluePixelChannel].correlation[i]=
1238 (correlation.direction[i].blue-mean.direction[i].blue*
1239 mean.direction[i].blue)/(sqrt(sum_squares.direction[i].blue-
1240 (mean.direction[i].blue*mean.direction[i].blue))*sqrt(
1241 sum_squares.direction[i].blue-(mean.direction[i].blue*
1242 mean.direction[i].blue)));
1243 if (image->colorspace == CMYKColorspace)
1244 channel_features[BlackPixelChannel].correlation[i]=
1245 (correlation.direction[i].black-mean.direction[i].black*
1246 mean.direction[i].black)/(sqrt(sum_squares.direction[i].black-
1247 (mean.direction[i].black*mean.direction[i].black))*sqrt(
1248 sum_squares.direction[i].black-(mean.direction[i].black*
1249 mean.direction[i].black)));
1250 if (image->alpha_trait != UndefinedPixelTrait)
1251 channel_features[AlphaPixelChannel].correlation[i]=
1252 (correlation.direction[i].alpha-mean.direction[i].alpha*
1253 mean.direction[i].alpha)/(sqrt(sum_squares.direction[i].alpha-
1254 (mean.direction[i].alpha*mean.direction[i].alpha))*sqrt(
1255 sum_squares.direction[i].alpha-(mean.direction[i].alpha*
1256 mean.direction[i].alpha)));
1257 }
1258 /*
1259 Compute more texture features.
1260 */
1261 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1262 #pragma omp parallel for schedule(static) shared(status) \
1263 magick_number_threads(image,image,number_grays,1)
1264 #endif
1265 for (i=0; i < 4; i++)
1266 {
1267 register ssize_t
1268 x;
1269
1270 for (x=2; x < (ssize_t) (2*number_grays); x++)
1271 {
1272 /*
1273 Sum average.
1274 */
1275 channel_features[RedPixelChannel].sum_average[i]+=
1276 x*density_xy[x].direction[i].red;
1277 channel_features[GreenPixelChannel].sum_average[i]+=
1278 x*density_xy[x].direction[i].green;
1279 channel_features[BluePixelChannel].sum_average[i]+=
1280 x*density_xy[x].direction[i].blue;
1281 if (image->colorspace == CMYKColorspace)
1282 channel_features[BlackPixelChannel].sum_average[i]+=
1283 x*density_xy[x].direction[i].black;
1284 if (image->alpha_trait != UndefinedPixelTrait)
1285 channel_features[AlphaPixelChannel].sum_average[i]+=
1286 x*density_xy[x].direction[i].alpha;
1287 /*
1288 Sum entropy.
1289 */
1290 channel_features[RedPixelChannel].sum_entropy[i]-=
1291 density_xy[x].direction[i].red*
1292 MagickLog10(density_xy[x].direction[i].red);
1293 channel_features[GreenPixelChannel].sum_entropy[i]-=
1294 density_xy[x].direction[i].green*
1295 MagickLog10(density_xy[x].direction[i].green);
1296 channel_features[BluePixelChannel].sum_entropy[i]-=
1297 density_xy[x].direction[i].blue*
1298 MagickLog10(density_xy[x].direction[i].blue);
1299 if (image->colorspace == CMYKColorspace)
1300 channel_features[BlackPixelChannel].sum_entropy[i]-=
1301 density_xy[x].direction[i].black*
1302 MagickLog10(density_xy[x].direction[i].black);
1303 if (image->alpha_trait != UndefinedPixelTrait)
1304 channel_features[AlphaPixelChannel].sum_entropy[i]-=
1305 density_xy[x].direction[i].alpha*
1306 MagickLog10(density_xy[x].direction[i].alpha);
1307 /*
1308 Sum variance.
1309 */
1310 channel_features[RedPixelChannel].sum_variance[i]+=
1311 (x-channel_features[RedPixelChannel].sum_entropy[i])*
1312 (x-channel_features[RedPixelChannel].sum_entropy[i])*
1313 density_xy[x].direction[i].red;
1314 channel_features[GreenPixelChannel].sum_variance[i]+=
1315 (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1316 (x-channel_features[GreenPixelChannel].sum_entropy[i])*
1317 density_xy[x].direction[i].green;
1318 channel_features[BluePixelChannel].sum_variance[i]+=
1319 (x-channel_features[BluePixelChannel].sum_entropy[i])*
1320 (x-channel_features[BluePixelChannel].sum_entropy[i])*
1321 density_xy[x].direction[i].blue;
1322 if (image->colorspace == CMYKColorspace)
1323 channel_features[BlackPixelChannel].sum_variance[i]+=
1324 (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1325 (x-channel_features[BlackPixelChannel].sum_entropy[i])*
1326 density_xy[x].direction[i].black;
1327 if (image->alpha_trait != UndefinedPixelTrait)
1328 channel_features[AlphaPixelChannel].sum_variance[i]+=
1329 (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1330 (x-channel_features[AlphaPixelChannel].sum_entropy[i])*
1331 density_xy[x].direction[i].alpha;
1332 }
1333 }
1334 /*
1335 Compute more texture features.
1336 */
1337 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1338 #pragma omp parallel for schedule(static) shared(status) \
1339 magick_number_threads(image,image,number_grays,1)
1340 #endif
1341 for (i=0; i < 4; i++)
1342 {
1343 register ssize_t
1344 y;
1345
1346 for (y=0; y < (ssize_t) number_grays; y++)
1347 {
1348 register ssize_t
1349 x;
1350
1351 for (x=0; x < (ssize_t) number_grays; x++)
1352 {
1353 /*
1354 Sum of Squares: Variance
1355 */
1356 variance.direction[i].red+=(y-mean.direction[i].red+1)*
1357 (y-mean.direction[i].red+1)*cooccurrence[x][y].direction[i].red;
1358 variance.direction[i].green+=(y-mean.direction[i].green+1)*
1359 (y-mean.direction[i].green+1)*cooccurrence[x][y].direction[i].green;
1360 variance.direction[i].blue+=(y-mean.direction[i].blue+1)*
1361 (y-mean.direction[i].blue+1)*cooccurrence[x][y].direction[i].blue;
1362 if (image->colorspace == CMYKColorspace)
1363 variance.direction[i].black+=(y-mean.direction[i].black+1)*
1364 (y-mean.direction[i].black+1)*cooccurrence[x][y].direction[i].black;
1365 if (image->alpha_trait != UndefinedPixelTrait)
1366 variance.direction[i].alpha+=(y-mean.direction[i].alpha+1)*
1367 (y-mean.direction[i].alpha+1)*
1368 cooccurrence[x][y].direction[i].alpha;
1369 /*
1370 Sum average / Difference Variance.
1371 */
1372 density_xy[MagickAbsoluteValue(y-x)].direction[i].red+=
1373 cooccurrence[x][y].direction[i].red;
1374 density_xy[MagickAbsoluteValue(y-x)].direction[i].green+=
1375 cooccurrence[x][y].direction[i].green;
1376 density_xy[MagickAbsoluteValue(y-x)].direction[i].blue+=
1377 cooccurrence[x][y].direction[i].blue;
1378 if (image->colorspace == CMYKColorspace)
1379 density_xy[MagickAbsoluteValue(y-x)].direction[i].black+=
1380 cooccurrence[x][y].direction[i].black;
1381 if (image->alpha_trait != UndefinedPixelTrait)
1382 density_xy[MagickAbsoluteValue(y-x)].direction[i].alpha+=
1383 cooccurrence[x][y].direction[i].alpha;
1384 /*
1385 Information Measures of Correlation.
1386 */
1387 entropy_xy.direction[i].red-=cooccurrence[x][y].direction[i].red*
1388 MagickLog10(cooccurrence[x][y].direction[i].red);
1389 entropy_xy.direction[i].green-=cooccurrence[x][y].direction[i].green*
1390 MagickLog10(cooccurrence[x][y].direction[i].green);
1391 entropy_xy.direction[i].blue-=cooccurrence[x][y].direction[i].blue*
1392 MagickLog10(cooccurrence[x][y].direction[i].blue);
1393 if (image->colorspace == CMYKColorspace)
1394 entropy_xy.direction[i].black-=cooccurrence[x][y].direction[i].black*
1395 MagickLog10(cooccurrence[x][y].direction[i].black);
1396 if (image->alpha_trait != UndefinedPixelTrait)
1397 entropy_xy.direction[i].alpha-=
1398 cooccurrence[x][y].direction[i].alpha*MagickLog10(
1399 cooccurrence[x][y].direction[i].alpha);
1400 entropy_xy1.direction[i].red-=(cooccurrence[x][y].direction[i].red*
1401 MagickLog10(density_x[x].direction[i].red*density_y[y].direction[i].red));
1402 entropy_xy1.direction[i].green-=(cooccurrence[x][y].direction[i].green*
1403 MagickLog10(density_x[x].direction[i].green*
1404 density_y[y].direction[i].green));
1405 entropy_xy1.direction[i].blue-=(cooccurrence[x][y].direction[i].blue*
1406 MagickLog10(density_x[x].direction[i].blue*density_y[y].direction[i].blue));
1407 if (image->colorspace == CMYKColorspace)
1408 entropy_xy1.direction[i].black-=(
1409 cooccurrence[x][y].direction[i].black*MagickLog10(
1410 density_x[x].direction[i].black*density_y[y].direction[i].black));
1411 if (image->alpha_trait != UndefinedPixelTrait)
1412 entropy_xy1.direction[i].alpha-=(
1413 cooccurrence[x][y].direction[i].alpha*MagickLog10(
1414 density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1415 entropy_xy2.direction[i].red-=(density_x[x].direction[i].red*
1416 density_y[y].direction[i].red*MagickLog10(density_x[x].direction[i].red*
1417 density_y[y].direction[i].red));
1418 entropy_xy2.direction[i].green-=(density_x[x].direction[i].green*
1419 density_y[y].direction[i].green*MagickLog10(density_x[x].direction[i].green*
1420 density_y[y].direction[i].green));
1421 entropy_xy2.direction[i].blue-=(density_x[x].direction[i].blue*
1422 density_y[y].direction[i].blue*MagickLog10(density_x[x].direction[i].blue*
1423 density_y[y].direction[i].blue));
1424 if (image->colorspace == CMYKColorspace)
1425 entropy_xy2.direction[i].black-=(density_x[x].direction[i].black*
1426 density_y[y].direction[i].black*MagickLog10(
1427 density_x[x].direction[i].black*density_y[y].direction[i].black));
1428 if (image->alpha_trait != UndefinedPixelTrait)
1429 entropy_xy2.direction[i].alpha-=(density_x[x].direction[i].alpha*
1430 density_y[y].direction[i].alpha*MagickLog10(
1431 density_x[x].direction[i].alpha*density_y[y].direction[i].alpha));
1432 }
1433 }
1434 channel_features[RedPixelChannel].variance_sum_of_squares[i]=
1435 variance.direction[i].red;
1436 channel_features[GreenPixelChannel].variance_sum_of_squares[i]=
1437 variance.direction[i].green;
1438 channel_features[BluePixelChannel].variance_sum_of_squares[i]=
1439 variance.direction[i].blue;
1440 if (image->colorspace == CMYKColorspace)
1441 channel_features[BlackPixelChannel].variance_sum_of_squares[i]=
1442 variance.direction[i].black;
1443 if (image->alpha_trait != UndefinedPixelTrait)
1444 channel_features[AlphaPixelChannel].variance_sum_of_squares[i]=
1445 variance.direction[i].alpha;
1446 }
1447 /*
1448 Compute more texture features.
1449 */
1450 (void) memset(&variance,0,sizeof(variance));
1451 (void) memset(&sum_squares,0,sizeof(sum_squares));
1452 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1453 #pragma omp parallel for schedule(static) shared(status) \
1454 magick_number_threads(image,image,number_grays,1)
1455 #endif
1456 for (i=0; i < 4; i++)
1457 {
1458 register ssize_t
1459 x;
1460
1461 for (x=0; x < (ssize_t) number_grays; x++)
1462 {
1463 /*
1464 Difference variance.
1465 */
1466 variance.direction[i].red+=density_xy[x].direction[i].red;
1467 variance.direction[i].green+=density_xy[x].direction[i].green;
1468 variance.direction[i].blue+=density_xy[x].direction[i].blue;
1469 if (image->colorspace == CMYKColorspace)
1470 variance.direction[i].black+=density_xy[x].direction[i].black;
1471 if (image->alpha_trait != UndefinedPixelTrait)
1472 variance.direction[i].alpha+=density_xy[x].direction[i].alpha;
1473 sum_squares.direction[i].red+=density_xy[x].direction[i].red*
1474 density_xy[x].direction[i].red;
1475 sum_squares.direction[i].green+=density_xy[x].direction[i].green*
1476 density_xy[x].direction[i].green;
1477 sum_squares.direction[i].blue+=density_xy[x].direction[i].blue*
1478 density_xy[x].direction[i].blue;
1479 if (image->colorspace == CMYKColorspace)
1480 sum_squares.direction[i].black+=density_xy[x].direction[i].black*
1481 density_xy[x].direction[i].black;
1482 if (image->alpha_trait != UndefinedPixelTrait)
1483 sum_squares.direction[i].alpha+=density_xy[x].direction[i].alpha*
1484 density_xy[x].direction[i].alpha;
1485 /*
1486 Difference entropy.
1487 */
1488 channel_features[RedPixelChannel].difference_entropy[i]-=
1489 density_xy[x].direction[i].red*
1490 MagickLog10(density_xy[x].direction[i].red);
1491 channel_features[GreenPixelChannel].difference_entropy[i]-=
1492 density_xy[x].direction[i].green*
1493 MagickLog10(density_xy[x].direction[i].green);
1494 channel_features[BluePixelChannel].difference_entropy[i]-=
1495 density_xy[x].direction[i].blue*
1496 MagickLog10(density_xy[x].direction[i].blue);
1497 if (image->colorspace == CMYKColorspace)
1498 channel_features[BlackPixelChannel].difference_entropy[i]-=
1499 density_xy[x].direction[i].black*
1500 MagickLog10(density_xy[x].direction[i].black);
1501 if (image->alpha_trait != UndefinedPixelTrait)
1502 channel_features[AlphaPixelChannel].difference_entropy[i]-=
1503 density_xy[x].direction[i].alpha*
1504 MagickLog10(density_xy[x].direction[i].alpha);
1505 /*
1506 Information Measures of Correlation.
1507 */
1508 entropy_x.direction[i].red-=(density_x[x].direction[i].red*
1509 MagickLog10(density_x[x].direction[i].red));
1510 entropy_x.direction[i].green-=(density_x[x].direction[i].green*
1511 MagickLog10(density_x[x].direction[i].green));
1512 entropy_x.direction[i].blue-=(density_x[x].direction[i].blue*
1513 MagickLog10(density_x[x].direction[i].blue));
1514 if (image->colorspace == CMYKColorspace)
1515 entropy_x.direction[i].black-=(density_x[x].direction[i].black*
1516 MagickLog10(density_x[x].direction[i].black));
1517 if (image->alpha_trait != UndefinedPixelTrait)
1518 entropy_x.direction[i].alpha-=(density_x[x].direction[i].alpha*
1519 MagickLog10(density_x[x].direction[i].alpha));
1520 entropy_y.direction[i].red-=(density_y[x].direction[i].red*
1521 MagickLog10(density_y[x].direction[i].red));
1522 entropy_y.direction[i].green-=(density_y[x].direction[i].green*
1523 MagickLog10(density_y[x].direction[i].green));
1524 entropy_y.direction[i].blue-=(density_y[x].direction[i].blue*
1525 MagickLog10(density_y[x].direction[i].blue));
1526 if (image->colorspace == CMYKColorspace)
1527 entropy_y.direction[i].black-=(density_y[x].direction[i].black*
1528 MagickLog10(density_y[x].direction[i].black));
1529 if (image->alpha_trait != UndefinedPixelTrait)
1530 entropy_y.direction[i].alpha-=(density_y[x].direction[i].alpha*
1531 MagickLog10(density_y[x].direction[i].alpha));
1532 }
1533 /*
1534 Difference variance.
1535 */
1536 channel_features[RedPixelChannel].difference_variance[i]=
1537 (((double) number_grays*number_grays*sum_squares.direction[i].red)-
1538 (variance.direction[i].red*variance.direction[i].red))/
1539 ((double) number_grays*number_grays*number_grays*number_grays);
1540 channel_features[GreenPixelChannel].difference_variance[i]=
1541 (((double) number_grays*number_grays*sum_squares.direction[i].green)-
1542 (variance.direction[i].green*variance.direction[i].green))/
1543 ((double) number_grays*number_grays*number_grays*number_grays);
1544 channel_features[BluePixelChannel].difference_variance[i]=
1545 (((double) number_grays*number_grays*sum_squares.direction[i].blue)-
1546 (variance.direction[i].blue*variance.direction[i].blue))/
1547 ((double) number_grays*number_grays*number_grays*number_grays);
1548 if (image->colorspace == CMYKColorspace)
1549 channel_features[BlackPixelChannel].difference_variance[i]=
1550 (((double) number_grays*number_grays*sum_squares.direction[i].black)-
1551 (variance.direction[i].black*variance.direction[i].black))/
1552 ((double) number_grays*number_grays*number_grays*number_grays);
1553 if (image->alpha_trait != UndefinedPixelTrait)
1554 channel_features[AlphaPixelChannel].difference_variance[i]=
1555 (((double) number_grays*number_grays*sum_squares.direction[i].alpha)-
1556 (variance.direction[i].alpha*variance.direction[i].alpha))/
1557 ((double) number_grays*number_grays*number_grays*number_grays);
1558 /*
1559 Information Measures of Correlation.
1560 */
1561 channel_features[RedPixelChannel].measure_of_correlation_1[i]=
1562 (entropy_xy.direction[i].red-entropy_xy1.direction[i].red)/
1563 (entropy_x.direction[i].red > entropy_y.direction[i].red ?
1564 entropy_x.direction[i].red : entropy_y.direction[i].red);
1565 channel_features[GreenPixelChannel].measure_of_correlation_1[i]=
1566 (entropy_xy.direction[i].green-entropy_xy1.direction[i].green)/
1567 (entropy_x.direction[i].green > entropy_y.direction[i].green ?
1568 entropy_x.direction[i].green : entropy_y.direction[i].green);
1569 channel_features[BluePixelChannel].measure_of_correlation_1[i]=
1570 (entropy_xy.direction[i].blue-entropy_xy1.direction[i].blue)/
1571 (entropy_x.direction[i].blue > entropy_y.direction[i].blue ?
1572 entropy_x.direction[i].blue : entropy_y.direction[i].blue);
1573 if (image->colorspace == CMYKColorspace)
1574 channel_features[BlackPixelChannel].measure_of_correlation_1[i]=
1575 (entropy_xy.direction[i].black-entropy_xy1.direction[i].black)/
1576 (entropy_x.direction[i].black > entropy_y.direction[i].black ?
1577 entropy_x.direction[i].black : entropy_y.direction[i].black);
1578 if (image->alpha_trait != UndefinedPixelTrait)
1579 channel_features[AlphaPixelChannel].measure_of_correlation_1[i]=
1580 (entropy_xy.direction[i].alpha-entropy_xy1.direction[i].alpha)/
1581 (entropy_x.direction[i].alpha > entropy_y.direction[i].alpha ?
1582 entropy_x.direction[i].alpha : entropy_y.direction[i].alpha);
1583 channel_features[RedPixelChannel].measure_of_correlation_2[i]=
1584 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].red-
1585 entropy_xy.direction[i].red)))));
1586 channel_features[GreenPixelChannel].measure_of_correlation_2[i]=
1587 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].green-
1588 entropy_xy.direction[i].green)))));
1589 channel_features[BluePixelChannel].measure_of_correlation_2[i]=
1590 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].blue-
1591 entropy_xy.direction[i].blue)))));
1592 if (image->colorspace == CMYKColorspace)
1593 channel_features[BlackPixelChannel].measure_of_correlation_2[i]=
1594 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].black-
1595 entropy_xy.direction[i].black)))));
1596 if (image->alpha_trait != UndefinedPixelTrait)
1597 channel_features[AlphaPixelChannel].measure_of_correlation_2[i]=
1598 (sqrt(fabs(1.0-exp(-2.0*(double) (entropy_xy2.direction[i].alpha-
1599 entropy_xy.direction[i].alpha)))));
1600 }
1601 /*
1602 Compute more texture features.
1603 */
1604 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1605 #pragma omp parallel for schedule(static) shared(status) \
1606 magick_number_threads(image,image,number_grays,1)
1607 #endif
1608 for (i=0; i < 4; i++)
1609 {
1610 ssize_t
1611 z;
1612
1613 for (z=0; z < (ssize_t) number_grays; z++)
1614 {
1615 register ssize_t
1616 y;
1617
1618 ChannelStatistics
1619 pixel;
1620
1621 (void) memset(&pixel,0,sizeof(pixel));
1622 for (y=0; y < (ssize_t) number_grays; y++)
1623 {
1624 register ssize_t
1625 x;
1626
1627 for (x=0; x < (ssize_t) number_grays; x++)
1628 {
1629 /*
1630 Contrast: amount of local variations present in an image.
1631 */
1632 if (((y-x) == z) || ((x-y) == z))
1633 {
1634 pixel.direction[i].red+=cooccurrence[x][y].direction[i].red;
1635 pixel.direction[i].green+=cooccurrence[x][y].direction[i].green;
1636 pixel.direction[i].blue+=cooccurrence[x][y].direction[i].blue;
1637 if (image->colorspace == CMYKColorspace)
1638 pixel.direction[i].black+=cooccurrence[x][y].direction[i].black;
1639 if (image->alpha_trait != UndefinedPixelTrait)
1640 pixel.direction[i].alpha+=
1641 cooccurrence[x][y].direction[i].alpha;
1642 }
1643 /*
1644 Maximum Correlation Coefficient.
1645 */
1646 if ((fabs(density_x[z].direction[i].red) > MagickEpsilon) &&
1647 (fabs(density_y[x].direction[i].red) > MagickEpsilon))
1648 Q[z][y].direction[i].red+=cooccurrence[z][x].direction[i].red*
1649 cooccurrence[y][x].direction[i].red/density_x[z].direction[i].red/
1650 density_y[x].direction[i].red;
1651 if ((fabs(density_x[z].direction[i].green) > MagickEpsilon) &&
1652 (fabs(density_y[x].direction[i].red) > MagickEpsilon))
1653 Q[z][y].direction[i].green+=cooccurrence[z][x].direction[i].green*
1654 cooccurrence[y][x].direction[i].green/
1655 density_x[z].direction[i].green/density_y[x].direction[i].red;
1656 if ((fabs(density_x[z].direction[i].blue) > MagickEpsilon) &&
1657 (fabs(density_y[x].direction[i].blue) > MagickEpsilon))
1658 Q[z][y].direction[i].blue+=cooccurrence[z][x].direction[i].blue*
1659 cooccurrence[y][x].direction[i].blue/
1660 density_x[z].direction[i].blue/density_y[x].direction[i].blue;
1661 if (image->colorspace == CMYKColorspace)
1662 if ((fabs(density_x[z].direction[i].black) > MagickEpsilon) &&
1663 (fabs(density_y[x].direction[i].black) > MagickEpsilon))
1664 Q[z][y].direction[i].black+=cooccurrence[z][x].direction[i].black*
1665 cooccurrence[y][x].direction[i].black/
1666 density_x[z].direction[i].black/density_y[x].direction[i].black;
1667 if (image->alpha_trait != UndefinedPixelTrait)
1668 if ((fabs(density_x[z].direction[i].alpha) > MagickEpsilon) &&
1669 (fabs(density_y[x].direction[i].alpha) > MagickEpsilon))
1670 Q[z][y].direction[i].alpha+=
1671 cooccurrence[z][x].direction[i].alpha*
1672 cooccurrence[y][x].direction[i].alpha/
1673 density_x[z].direction[i].alpha/
1674 density_y[x].direction[i].alpha;
1675 }
1676 }
1677 channel_features[RedPixelChannel].contrast[i]+=z*z*
1678 pixel.direction[i].red;
1679 channel_features[GreenPixelChannel].contrast[i]+=z*z*
1680 pixel.direction[i].green;
1681 channel_features[BluePixelChannel].contrast[i]+=z*z*
1682 pixel.direction[i].blue;
1683 if (image->colorspace == CMYKColorspace)
1684 channel_features[BlackPixelChannel].contrast[i]+=z*z*
1685 pixel.direction[i].black;
1686 if (image->alpha_trait != UndefinedPixelTrait)
1687 channel_features[AlphaPixelChannel].contrast[i]+=z*z*
1688 pixel.direction[i].alpha;
1689 }
1690 /*
1691 Maximum Correlation Coefficient.
1692 Future: return second largest eigenvalue of Q.
1693 */
1694 channel_features[RedPixelChannel].maximum_correlation_coefficient[i]=
1695 sqrt((double) -1.0);
1696 channel_features[GreenPixelChannel].maximum_correlation_coefficient[i]=
1697 sqrt((double) -1.0);
1698 channel_features[BluePixelChannel].maximum_correlation_coefficient[i]=
1699 sqrt((double) -1.0);
1700 if (image->colorspace == CMYKColorspace)
1701 channel_features[BlackPixelChannel].maximum_correlation_coefficient[i]=
1702 sqrt((double) -1.0);
1703 if (image->alpha_trait != UndefinedPixelTrait)
1704 channel_features[AlphaPixelChannel].maximum_correlation_coefficient[i]=
1705 sqrt((double) -1.0);
1706 }
1707 /*
1708 Relinquish resources.
1709 */
1710 sum=(ChannelStatistics *) RelinquishMagickMemory(sum);
1711 for (i=0; i < (ssize_t) number_grays; i++)
1712 Q[i]=(ChannelStatistics *) RelinquishMagickMemory(Q[i]);
1713 Q=(ChannelStatistics **) RelinquishMagickMemory(Q);
1714 density_y=(ChannelStatistics *) RelinquishMagickMemory(density_y);
1715 density_xy=(ChannelStatistics *) RelinquishMagickMemory(density_xy);
1716 density_x=(ChannelStatistics *) RelinquishMagickMemory(density_x);
1717 for (i=0; i < (ssize_t) number_grays; i++)
1718 cooccurrence[i]=(ChannelStatistics *)
1719 RelinquishMagickMemory(cooccurrence[i]);
1720 cooccurrence=(ChannelStatistics **) RelinquishMagickMemory(cooccurrence);
1721 return(channel_features);
1722 }
1723
1724 /*
1725 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1726 % %
1727 % %
1728 % %
1729 % H o u g h L i n e I m a g e %
1730 % %
1731 % %
1732 % %
1733 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
1734 %
1735 % Use HoughLineImage() in conjunction with any binary edge extracted image (we
1736 % recommand Canny) to identify lines in the image. The algorithm accumulates
1737 % counts for every white pixel for every possible orientation (for angles from
1738 % 0 to 179 in 1 degree increments) and distance from the center of the image to
1739 % the corner (in 1 px increments) and stores the counts in an accumulator
1740 % matrix of angle vs distance. The size of the accumulator is 180x(diagonal/2).
1741 % Next it searches this space for peaks in counts and converts the locations
1742 % of the peaks to slope and intercept in the normal x,y input image space. Use
1743 % the slope/intercepts to find the endpoints clipped to the bounds of the
1744 % image. The lines are then drawn. The counts are a measure of the length of
1745 % the lines.
1746 %
1747 % The format of the HoughLineImage method is:
1748 %
1749 % Image *HoughLineImage(const Image *image,const size_t width,
1750 % const size_t height,const size_t threshold,ExceptionInfo *exception)
1751 %
1752 % A description of each parameter follows:
1753 %
1754 % o image: the image.
1755 %
1756 % o width, height: find line pairs as local maxima in this neighborhood.
1757 %
1758 % o threshold: the line count threshold.
1759 %
1760 % o exception: return any errors or warnings in this structure.
1761 %
1762 */
1763
MagickRound(double x)1764 static inline double MagickRound(double x)
1765 {
1766 /*
1767 Round the fraction to nearest integer.
1768 */
1769 if ((x-floor(x)) < (ceil(x)-x))
1770 return(floor(x));
1771 return(ceil(x));
1772 }
1773
RenderHoughLines(const ImageInfo * image_info,const size_t columns,const size_t rows,ExceptionInfo * exception)1774 static Image *RenderHoughLines(const ImageInfo *image_info,const size_t columns,
1775 const size_t rows,ExceptionInfo *exception)
1776 {
1777 #define BoundingBox "viewbox"
1778
1779 DrawInfo
1780 *draw_info;
1781
1782 Image
1783 *image;
1784
1785 MagickBooleanType
1786 status;
1787
1788 /*
1789 Open image.
1790 */
1791 image=AcquireImage(image_info,exception);
1792 status=OpenBlob(image_info,image,ReadBinaryBlobMode,exception);
1793 if (status == MagickFalse)
1794 {
1795 image=DestroyImageList(image);
1796 return((Image *) NULL);
1797 }
1798 image->columns=columns;
1799 image->rows=rows;
1800 draw_info=CloneDrawInfo(image_info,(DrawInfo *) NULL);
1801 draw_info->affine.sx=image->resolution.x == 0.0 ? 1.0 : image->resolution.x/
1802 DefaultResolution;
1803 draw_info->affine.sy=image->resolution.y == 0.0 ? 1.0 : image->resolution.y/
1804 DefaultResolution;
1805 image->columns=(size_t) (draw_info->affine.sx*image->columns);
1806 image->rows=(size_t) (draw_info->affine.sy*image->rows);
1807 status=SetImageExtent(image,image->columns,image->rows,exception);
1808 if (status == MagickFalse)
1809 return(DestroyImageList(image));
1810 if (SetImageBackgroundColor(image,exception) == MagickFalse)
1811 {
1812 image=DestroyImageList(image);
1813 return((Image *) NULL);
1814 }
1815 /*
1816 Render drawing.
1817 */
1818 if (GetBlobStreamData(image) == (unsigned char *) NULL)
1819 draw_info->primitive=FileToString(image->filename,~0UL,exception);
1820 else
1821 {
1822 draw_info->primitive=(char *) AcquireMagickMemory((size_t)
1823 GetBlobSize(image)+1);
1824 if (draw_info->primitive != (char *) NULL)
1825 {
1826 (void) memcpy(draw_info->primitive,GetBlobStreamData(image),
1827 (size_t) GetBlobSize(image));
1828 draw_info->primitive[GetBlobSize(image)]='\0';
1829 }
1830 }
1831 (void) DrawImage(image,draw_info,exception);
1832 draw_info=DestroyDrawInfo(draw_info);
1833 (void) CloseBlob(image);
1834 return(GetFirstImageInList(image));
1835 }
1836
HoughLineImage(const Image * image,const size_t width,const size_t height,const size_t threshold,ExceptionInfo * exception)1837 MagickExport Image *HoughLineImage(const Image *image,const size_t width,
1838 const size_t height,const size_t threshold,ExceptionInfo *exception)
1839 {
1840 #define HoughLineImageTag "HoughLine/Image"
1841
1842 CacheView
1843 *image_view;
1844
1845 char
1846 message[MagickPathExtent],
1847 path[MagickPathExtent];
1848
1849 const char
1850 *artifact;
1851
1852 double
1853 hough_height;
1854
1855 Image
1856 *lines_image = NULL;
1857
1858 ImageInfo
1859 *image_info;
1860
1861 int
1862 file;
1863
1864 MagickBooleanType
1865 status;
1866
1867 MagickOffsetType
1868 progress;
1869
1870 MatrixInfo
1871 *accumulator;
1872
1873 PointInfo
1874 center;
1875
1876 register ssize_t
1877 y;
1878
1879 size_t
1880 accumulator_height,
1881 accumulator_width,
1882 line_count;
1883
1884 /*
1885 Create the accumulator.
1886 */
1887 assert(image != (const Image *) NULL);
1888 assert(image->signature == MagickCoreSignature);
1889 if (image->debug != MagickFalse)
1890 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
1891 assert(exception != (ExceptionInfo *) NULL);
1892 assert(exception->signature == MagickCoreSignature);
1893 accumulator_width=180;
1894 hough_height=((sqrt(2.0)*(double) (image->rows > image->columns ?
1895 image->rows : image->columns))/2.0);
1896 accumulator_height=(size_t) (2.0*hough_height);
1897 accumulator=AcquireMatrixInfo(accumulator_width,accumulator_height,
1898 sizeof(double),exception);
1899 if (accumulator == (MatrixInfo *) NULL)
1900 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1901 if (NullMatrix(accumulator) == MagickFalse)
1902 {
1903 accumulator=DestroyMatrixInfo(accumulator);
1904 ThrowImageException(ResourceLimitError,"MemoryAllocationFailed");
1905 }
1906 /*
1907 Populate the accumulator.
1908 */
1909 status=MagickTrue;
1910 progress=0;
1911 center.x=(double) image->columns/2.0;
1912 center.y=(double) image->rows/2.0;
1913 image_view=AcquireVirtualCacheView(image,exception);
1914 for (y=0; y < (ssize_t) image->rows; y++)
1915 {
1916 register const Quantum
1917 *magick_restrict p;
1918
1919 register ssize_t
1920 x;
1921
1922 if (status == MagickFalse)
1923 continue;
1924 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
1925 if (p == (Quantum *) NULL)
1926 {
1927 status=MagickFalse;
1928 continue;
1929 }
1930 for (x=0; x < (ssize_t) image->columns; x++)
1931 {
1932 if (GetPixelIntensity(image,p) > (QuantumRange/2.0))
1933 {
1934 register ssize_t
1935 i;
1936
1937 for (i=0; i < 180; i++)
1938 {
1939 double
1940 count,
1941 radius;
1942
1943 radius=(((double) x-center.x)*cos(DegreesToRadians((double) i)))+
1944 (((double) y-center.y)*sin(DegreesToRadians((double) i)));
1945 (void) GetMatrixElement(accumulator,i,(ssize_t)
1946 MagickRound(radius+hough_height),&count);
1947 count++;
1948 (void) SetMatrixElement(accumulator,i,(ssize_t)
1949 MagickRound(radius+hough_height),&count);
1950 }
1951 }
1952 p+=GetPixelChannels(image);
1953 }
1954 if (image->progress_monitor != (MagickProgressMonitor) NULL)
1955 {
1956 MagickBooleanType
1957 proceed;
1958
1959 #if defined(MAGICKCORE_OPENMP_SUPPORT)
1960 #pragma omp atomic
1961 #endif
1962 progress++;
1963 proceed=SetImageProgress(image,CannyEdgeImageTag,progress,image->rows);
1964 if (proceed == MagickFalse)
1965 status=MagickFalse;
1966 }
1967 }
1968 image_view=DestroyCacheView(image_view);
1969 if (status == MagickFalse)
1970 {
1971 accumulator=DestroyMatrixInfo(accumulator);
1972 return((Image *) NULL);
1973 }
1974 /*
1975 Generate line segments from accumulator.
1976 */
1977 file=AcquireUniqueFileResource(path);
1978 if (file == -1)
1979 {
1980 accumulator=DestroyMatrixInfo(accumulator);
1981 return((Image *) NULL);
1982 }
1983 (void) FormatLocaleString(message,MagickPathExtent,
1984 "# Hough line transform: %.20gx%.20g%+.20g\n",(double) width,
1985 (double) height,(double) threshold);
1986 if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1987 status=MagickFalse;
1988 (void) FormatLocaleString(message,MagickPathExtent,
1989 "viewbox 0 0 %.20g %.20g\n",(double) image->columns,(double) image->rows);
1990 if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1991 status=MagickFalse;
1992 (void) FormatLocaleString(message,MagickPathExtent,
1993 "# x1,y1 x2,y2 # count angle distance\n");
1994 if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
1995 status=MagickFalse;
1996 line_count=image->columns > image->rows ? image->columns/4 : image->rows/4;
1997 if (threshold != 0)
1998 line_count=threshold;
1999 for (y=0; y < (ssize_t) accumulator_height; y++)
2000 {
2001 register ssize_t
2002 x;
2003
2004 for (x=0; x < (ssize_t) accumulator_width; x++)
2005 {
2006 double
2007 count;
2008
2009 (void) GetMatrixElement(accumulator,x,y,&count);
2010 if (count >= (double) line_count)
2011 {
2012 double
2013 maxima;
2014
2015 SegmentInfo
2016 line;
2017
2018 ssize_t
2019 v;
2020
2021 /*
2022 Is point a local maxima?
2023 */
2024 maxima=count;
2025 for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++)
2026 {
2027 ssize_t
2028 u;
2029
2030 for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++)
2031 {
2032 if ((u != 0) || (v !=0))
2033 {
2034 (void) GetMatrixElement(accumulator,x+u,y+v,&count);
2035 if (count > maxima)
2036 {
2037 maxima=count;
2038 break;
2039 }
2040 }
2041 }
2042 if (u < (ssize_t) (width/2))
2043 break;
2044 }
2045 (void) GetMatrixElement(accumulator,x,y,&count);
2046 if (maxima > count)
2047 continue;
2048 if ((x >= 45) && (x <= 135))
2049 {
2050 /*
2051 y = (r-x cos(t))/sin(t)
2052 */
2053 line.x1=0.0;
2054 line.y1=((double) (y-(accumulator_height/2.0))-((line.x1-
2055 (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
2056 sin(DegreesToRadians((double) x))+(image->rows/2.0);
2057 line.x2=(double) image->columns;
2058 line.y2=((double) (y-(accumulator_height/2.0))-((line.x2-
2059 (image->columns/2.0))*cos(DegreesToRadians((double) x))))/
2060 sin(DegreesToRadians((double) x))+(image->rows/2.0);
2061 }
2062 else
2063 {
2064 /*
2065 x = (r-y cos(t))/sin(t)
2066 */
2067 line.y1=0.0;
2068 line.x1=((double) (y-(accumulator_height/2.0))-((line.y1-
2069 (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
2070 cos(DegreesToRadians((double) x))+(image->columns/2.0);
2071 line.y2=(double) image->rows;
2072 line.x2=((double) (y-(accumulator_height/2.0))-((line.y2-
2073 (image->rows/2.0))*sin(DegreesToRadians((double) x))))/
2074 cos(DegreesToRadians((double) x))+(image->columns/2.0);
2075 }
2076 (void) FormatLocaleString(message,MagickPathExtent,
2077 "line %g,%g %g,%g # %g %g %g\n",line.x1,line.y1,line.x2,line.y2,
2078 maxima,(double) x,(double) y);
2079 if (write(file,message,strlen(message)) != (ssize_t) strlen(message))
2080 status=MagickFalse;
2081 }
2082 }
2083 }
2084 (void) close(file);
2085 /*
2086 Render lines to image canvas.
2087 */
2088 image_info=AcquireImageInfo();
2089 image_info->background_color=image->background_color;
2090 (void) FormatLocaleString(image_info->filename,MagickPathExtent,"%s",path);
2091 artifact=GetImageArtifact(image,"background");
2092 if (artifact != (const char *) NULL)
2093 (void) SetImageOption(image_info,"background",artifact);
2094 artifact=GetImageArtifact(image,"fill");
2095 if (artifact != (const char *) NULL)
2096 (void) SetImageOption(image_info,"fill",artifact);
2097 artifact=GetImageArtifact(image,"stroke");
2098 if (artifact != (const char *) NULL)
2099 (void) SetImageOption(image_info,"stroke",artifact);
2100 artifact=GetImageArtifact(image,"strokewidth");
2101 if (artifact != (const char *) NULL)
2102 (void) SetImageOption(image_info,"strokewidth",artifact);
2103 lines_image=RenderHoughLines(image_info,image->columns,image->rows,exception);
2104 artifact=GetImageArtifact(image,"hough-lines:accumulator");
2105 if ((lines_image != (Image *) NULL) &&
2106 (IsStringTrue(artifact) != MagickFalse))
2107 {
2108 Image
2109 *accumulator_image;
2110
2111 accumulator_image=MatrixToImage(accumulator,exception);
2112 if (accumulator_image != (Image *) NULL)
2113 AppendImageToList(&lines_image,accumulator_image);
2114 }
2115 /*
2116 Free resources.
2117 */
2118 accumulator=DestroyMatrixInfo(accumulator);
2119 image_info=DestroyImageInfo(image_info);
2120 (void) RelinquishUniqueFileResource(path);
2121 return(GetFirstImageInList(lines_image));
2122 }
2123
2124 /*
2125 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2126 % %
2127 % %
2128 % %
2129 % M e a n S h i f t I m a g e %
2130 % %
2131 % %
2132 % %
2133 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2134 %
2135 % MeanShiftImage() delineate arbitrarily shaped clusters in the image. For
2136 % each pixel, it visits all the pixels in the neighborhood specified by
2137 % the window centered at the pixel and excludes those that are outside the
2138 % radius=(window-1)/2 surrounding the pixel. From those pixels, it finds those
2139 % that are within the specified color distance from the current mean, and
2140 % computes a new x,y centroid from those coordinates and a new mean. This new
2141 % x,y centroid is used as the center for a new window. This process iterates
2142 % until it converges and the final mean is replaces the (original window
2143 % center) pixel value. It repeats this process for the next pixel, etc.,
2144 % until it processes all pixels in the image. Results are typically better with
2145 % colorspaces other than sRGB. We recommend YIQ, YUV or YCbCr.
2146 %
2147 % The format of the MeanShiftImage method is:
2148 %
2149 % Image *MeanShiftImage(const Image *image,const size_t width,
2150 % const size_t height,const double color_distance,
2151 % ExceptionInfo *exception)
2152 %
2153 % A description of each parameter follows:
2154 %
2155 % o image: the image.
2156 %
2157 % o width, height: find pixels in this neighborhood.
2158 %
2159 % o color_distance: the color distance.
2160 %
2161 % o exception: return any errors or warnings in this structure.
2162 %
2163 */
MeanShiftImage(const Image * image,const size_t width,const size_t height,const double color_distance,ExceptionInfo * exception)2164 MagickExport Image *MeanShiftImage(const Image *image,const size_t width,
2165 const size_t height,const double color_distance,ExceptionInfo *exception)
2166 {
2167 #define MaxMeanShiftIterations 100
2168 #define MeanShiftImageTag "MeanShift/Image"
2169
2170 CacheView
2171 *image_view,
2172 *mean_view,
2173 *pixel_view;
2174
2175 Image
2176 *mean_image;
2177
2178 MagickBooleanType
2179 status;
2180
2181 MagickOffsetType
2182 progress;
2183
2184 ssize_t
2185 y;
2186
2187 assert(image != (const Image *) NULL);
2188 assert(image->signature == MagickCoreSignature);
2189 if (image->debug != MagickFalse)
2190 (void) LogMagickEvent(TraceEvent,GetMagickModule(),"%s",image->filename);
2191 assert(exception != (ExceptionInfo *) NULL);
2192 assert(exception->signature == MagickCoreSignature);
2193 mean_image=CloneImage(image,0,0,MagickTrue,exception);
2194 if (mean_image == (Image *) NULL)
2195 return((Image *) NULL);
2196 if (SetImageStorageClass(mean_image,DirectClass,exception) == MagickFalse)
2197 {
2198 mean_image=DestroyImage(mean_image);
2199 return((Image *) NULL);
2200 }
2201 status=MagickTrue;
2202 progress=0;
2203 image_view=AcquireVirtualCacheView(image,exception);
2204 pixel_view=AcquireVirtualCacheView(image,exception);
2205 mean_view=AcquireAuthenticCacheView(mean_image,exception);
2206 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2207 #pragma omp parallel for schedule(static) shared(status,progress) \
2208 magick_number_threads(mean_image,mean_image,mean_image->rows,1)
2209 #endif
2210 for (y=0; y < (ssize_t) mean_image->rows; y++)
2211 {
2212 register const Quantum
2213 *magick_restrict p;
2214
2215 register Quantum
2216 *magick_restrict q;
2217
2218 register ssize_t
2219 x;
2220
2221 if (status == MagickFalse)
2222 continue;
2223 p=GetCacheViewVirtualPixels(image_view,0,y,image->columns,1,exception);
2224 q=GetCacheViewAuthenticPixels(mean_view,0,y,mean_image->columns,1,
2225 exception);
2226 if ((p == (const Quantum *) NULL) || (q == (Quantum *) NULL))
2227 {
2228 status=MagickFalse;
2229 continue;
2230 }
2231 for (x=0; x < (ssize_t) mean_image->columns; x++)
2232 {
2233 PixelInfo
2234 mean_pixel,
2235 previous_pixel;
2236
2237 PointInfo
2238 mean_location,
2239 previous_location;
2240
2241 register ssize_t
2242 i;
2243
2244 GetPixelInfo(image,&mean_pixel);
2245 GetPixelInfoPixel(image,p,&mean_pixel);
2246 mean_location.x=(double) x;
2247 mean_location.y=(double) y;
2248 for (i=0; i < MaxMeanShiftIterations; i++)
2249 {
2250 double
2251 distance,
2252 gamma;
2253
2254 PixelInfo
2255 sum_pixel;
2256
2257 PointInfo
2258 sum_location;
2259
2260 ssize_t
2261 count,
2262 v;
2263
2264 sum_location.x=0.0;
2265 sum_location.y=0.0;
2266 GetPixelInfo(image,&sum_pixel);
2267 previous_location=mean_location;
2268 previous_pixel=mean_pixel;
2269 count=0;
2270 for (v=(-((ssize_t) height/2)); v <= (((ssize_t) height/2)); v++)
2271 {
2272 ssize_t
2273 u;
2274
2275 for (u=(-((ssize_t) width/2)); u <= (((ssize_t) width/2)); u++)
2276 {
2277 if ((v*v+u*u) <= (ssize_t) ((width/2)*(height/2)))
2278 {
2279 PixelInfo
2280 pixel;
2281
2282 status=GetOneCacheViewVirtualPixelInfo(pixel_view,(ssize_t)
2283 MagickRound(mean_location.x+u),(ssize_t) MagickRound(
2284 mean_location.y+v),&pixel,exception);
2285 distance=(mean_pixel.red-pixel.red)*(mean_pixel.red-pixel.red)+
2286 (mean_pixel.green-pixel.green)*(mean_pixel.green-pixel.green)+
2287 (mean_pixel.blue-pixel.blue)*(mean_pixel.blue-pixel.blue);
2288 if (distance <= (color_distance*color_distance))
2289 {
2290 sum_location.x+=mean_location.x+u;
2291 sum_location.y+=mean_location.y+v;
2292 sum_pixel.red+=pixel.red;
2293 sum_pixel.green+=pixel.green;
2294 sum_pixel.blue+=pixel.blue;
2295 sum_pixel.alpha+=pixel.alpha;
2296 count++;
2297 }
2298 }
2299 }
2300 }
2301 gamma=PerceptibleReciprocal(count);
2302 mean_location.x=gamma*sum_location.x;
2303 mean_location.y=gamma*sum_location.y;
2304 mean_pixel.red=gamma*sum_pixel.red;
2305 mean_pixel.green=gamma*sum_pixel.green;
2306 mean_pixel.blue=gamma*sum_pixel.blue;
2307 mean_pixel.alpha=gamma*sum_pixel.alpha;
2308 distance=(mean_location.x-previous_location.x)*
2309 (mean_location.x-previous_location.x)+
2310 (mean_location.y-previous_location.y)*
2311 (mean_location.y-previous_location.y)+
2312 255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)*
2313 255.0*QuantumScale*(mean_pixel.red-previous_pixel.red)+
2314 255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)*
2315 255.0*QuantumScale*(mean_pixel.green-previous_pixel.green)+
2316 255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue)*
2317 255.0*QuantumScale*(mean_pixel.blue-previous_pixel.blue);
2318 if (distance <= 3.0)
2319 break;
2320 }
2321 SetPixelRed(mean_image,ClampToQuantum(mean_pixel.red),q);
2322 SetPixelGreen(mean_image,ClampToQuantum(mean_pixel.green),q);
2323 SetPixelBlue(mean_image,ClampToQuantum(mean_pixel.blue),q);
2324 SetPixelAlpha(mean_image,ClampToQuantum(mean_pixel.alpha),q);
2325 p+=GetPixelChannels(image);
2326 q+=GetPixelChannels(mean_image);
2327 }
2328 if (SyncCacheViewAuthenticPixels(mean_view,exception) == MagickFalse)
2329 status=MagickFalse;
2330 if (image->progress_monitor != (MagickProgressMonitor) NULL)
2331 {
2332 MagickBooleanType
2333 proceed;
2334
2335 #if defined(MAGICKCORE_OPENMP_SUPPORT)
2336 #pragma omp atomic
2337 #endif
2338 progress++;
2339 proceed=SetImageProgress(image,MeanShiftImageTag,progress,image->rows);
2340 if (proceed == MagickFalse)
2341 status=MagickFalse;
2342 }
2343 }
2344 mean_view=DestroyCacheView(mean_view);
2345 pixel_view=DestroyCacheView(pixel_view);
2346 image_view=DestroyCacheView(image_view);
2347 return(mean_image);
2348 }
2349