1/* optional args to functions 2 */ 3 4get_option options defaults f 5 = error (_ "unknown parameter " ++ f), hits == [] 6 = hits?0 7{ 8 hits = [v :: [n, v] <- options ++ defaults; n == f]; 9} 10 11/* Various operators as functions. 12 */ 13 14logical_and a b = a && b; 15logical_or a b = a || b; 16bitwise_and a b = a & b; 17bitwise_or a b = a | b; 18eor a b = a ^ b; 19left_shift a b = a << b; 20right_shift a b = a >> b; 21not a = !a; 22 23less a b = a < b; 24more a b = a > b; 25less_equal a b = a <= b; 26more_equal a b = a >= b; 27equal a b = a == b; 28not_equal a b = a != b; 29pointer_equal a b = a === b; 30not_pointer_equal a b = a !== b; 31 32add a b = a + b; 33subtract a b = a - b; 34multiply a b = a * b; 35divide a b = a / b; 36idivide a b = (int) ((int) a / (int) b); 37power a b = a ** b; 38square x = x * x; 39remainder a b = a % b; 40 41cons a b = a : b; 42dot a b = a . ( b ); 43join a b = a ++ b; 44// 'difference' is defined in _list 45subscript a b = a ? b; 46 47generate s n f = [s, n .. f]; 48comma r i = (r, i); 49 50compose f g = f @ g; 51 52// our only trinary operator is actually a binary operator 53if_then_else a x = if a then x?0 else x?1; 54 55cast_unsigned_char x = (unsigned char) x; 56cast_signed_char x = (signed char) x; 57cast_unsigned_short x = (unsigned short) x; 58cast_signed_short x = (signed short) x; 59cast_unsigned_int x = (unsigned int) x; 60cast_signed_int x = (signed int) x; 61cast_float x = (float) x; 62cast_double x = (double) x; 63cast_complex x = (complex) x; 64cast_double_complex x = (double complex) x; 65 66unary_minus x = -x; 67negate x = !x; 68complement x = ~x; 69unary_plus x = +x; 70 71// the function we call for "a -> v" expressions 72mksvpair s v 73 = [s, v], is_string s 74 = error "not str on lhs of ->"; 75 76// the vector ops ... im is an image, vec is a real_list 77vec op_name im vec 78 = im_lintra_vec ones im vec, 79 op_name == "add" || op_name == "add'" 80 = im_lintra_vec ones (-1 * im) vec, 81 op_name == "subtract'" 82 = im_lintra_vec ones im inv, 83 op_name == "subtract" 84 = im_lintra_vec vec im zeros, 85 op_name == "multiply" || op_name == "multiply'" 86 = im_lintra_vec vec (1 / im) zeros, 87 op_name == "divide'" 88 = im_lintra_vec recip im zeros, 89 op_name == "divide" 90 = im_expntra_vec im vec, 91 op_name == "power'" 92 = im_powtra_vec im vec, 93 op_name == "power" 94 = im_remainderconst_vec im vec, 95 op_name == "remainder" 96 = im_andimage_vec im vec, 97 op_name == "bitwise_and" || op_name == "bitwise_and'" 98 = im_orimage_vec im vec, 99 op_name == "bitwise_or" || op_name == "bitwise_or'" 100 = im_eorimage_vec im vec, 101 op_name == "eor" || op_name == "eor'" 102 = im_equal_vec im vec, 103 op_name == "equal" || op_name == "equal'" 104 = im_notequal_vec im vec, 105 op_name == "not_equal" || op_name == "not_equal'" 106 = im_less_vec im vec, 107 op_name == "less" 108 = im_moreeq_vec im vec, 109 op_name == "less'" 110 = im_lesseq_vec im vec, 111 op_name == "less_equal" 112 = im_more_vec im vec, 113 op_name == "less_equal'" 114 = error ("unimplemented vector operation: " ++ op_name) 115{ 116 zeros = replicate (len vec) 0; 117 ones = replicate (len vec) 1; 118 recip = map (divide 1) vec; 119 inv = map (multiply (-1)) vec; 120} 121 122// make a name value pair 123mknvpair n v 124 = [n, v], is_string n 125 = error "not [char] on LHS of =>"; 126 127/* Macbeth chart patch names. 128 */ 129macbeth_names = [ 130 "Dark skin", 131 "Light skin", 132 "Blue sky", 133 "Foliage", 134 "Blue flower", 135 "Bluish green", 136 "Orange", 137 "Purplish blue", 138 "Moderate red", 139 "Purple", 140 "Yellow green", 141 "Orange yellow", 142 "Blue", 143 "Green", 144 "Red", 145 "Yellow", 146 "Magenta", 147 "Cyan", 148 "White (density 0.05)", 149 "Neutral 8 (density 0.23)", 150 "Neutral 6.5 (density 0.44)", 151 "Neutral 5 (density 0.70)", 152 "Neutral 3.5 (density 1.05)", 153 "Black (density 1.50)" 154]; 155 156bandsplit x 157 = oo_unary_function bandsplit_op x, is_class x 158 = map (subscript x) [0 .. bands - 1], is_image x 159 = error (_ "bad arguments to " ++ "bandsplit") 160{ 161 bands = get_header "Bands" x; 162 bandsplit_op = Operator "bandsplit" (map Image @ bandsplit) 163 Operator_type.COMPOUND false; 164} 165 166bandjoin l 167 = wrapper joined, has_wrapper 168 = joined, is_listof has_image l 169 = error (_ "bad arguments to " ++ "bandjoin") 170{ 171 has_wrapper = has_member_list (has_member "Image") l; 172 wrapper = get_member_list (has_member "Image") (get_member "Image") l; 173 joined = im_gbandjoin (map get_image l); 174} 175 176bandand x 177 = oo_unary_function bandand_op x, is_class x 178 = foldr1 bitwise_and (bandsplit x), is_image x 179 = error (_ "bad arguments to " ++ "bandand") 180{ 181 bandand_op = Operator "bandand" bandand Operator_type.COMPOUND_REWRAP false; 182} 183 184bandor x 185 = oo_unary_function bandor_op x, is_class x 186 = foldr1 bitwise_or (bandsplit x), is_image x 187 = error (_ "bad arguments to " ++ "bandor") 188{ 189 bandor_op = Operator "bandor" bandor Operator_type.COMPOUND_REWRAP false; 190} 191 192sum x 193 = oo_unary_function sum_op x, is_class x 194 = im_avg x * (get_width x) * (get_height x) * (get_bands x), is_image x 195 = sum_list x, is_list x 196 = error (_ "bad arguments (" ++ print x ++ ") to " ++ "sum") 197{ 198 sum_op = Operator "sum" sum Operator_type.COMPOUND false; 199 200 // add elements in a nested-list thing 201 sum_list l 202 = foldr acc 0 l 203 { 204 acc x total 205 = total + sum x, is_list x 206 = total + x; 207 } 208} 209 210product x 211 = oo_unary_function product_op x, is_class x 212 = product_list x, is_list x 213 // (product image) doesn't make much sense :( 214 = error (_ "bad arguments (" ++ print x ++ ") to " ++ "product") 215{ 216 product_op = Operator "product" product Operator_type.COMPOUND false; 217 218 product_list l 219 = foldr prod 1 l 220 { 221 prod x total 222 = total * product x, is_list x 223 = total * x; 224 } 225} 226 227mean x 228 = oo_unary_function mean_op x, is_class x 229 = im_avg x, is_image x 230 = mean_list x, is_list x 231 = error (_ "bad arguments (" ++ print x ++ ") to " ++ "mean") 232{ 233 mean_op = Operator "mean" mean Operator_type.COMPOUND false; 234 235 mean_list l = sum l / size l; 236 237 // number of elements in some sort of nested-list thing 238 size l 239 = foldr acc 0 l 240 { 241 acc x total 242 = total + size x, is_list x 243 = total + 1; 244 } 245} 246 247meang x 248 = (appl (power e) @ mean @ appl log) x 249{ 250 appl fn x 251 = map fn x, is_list x 252 = fn x; 253} 254 255skew x 256 = oo_unary_function skew_op x, is_class x 257 = sum ((x - m) ** 3) / ((N - 1) * s ** 3), is_image x 258 = sum ((Group x' - m) ** 3).value / ((N - 1) * s ** 3), is_list x 259 = error (_ "bad arguments (" ++ print x ++ ") to " ++ "skew") 260{ 261 skew_op = Operator "skew" skew Operator_type.COMPOUND false; 262 263 // squash any large matrix down to a flat list ... much simpler 264 x' 265 = x, is_image x; 266 = flatten x; 267 268 m = mean x'; 269 s = deviation x'; 270 w = get_width x'; 271 h = get_height x'; 272 b = get_bands x'; 273 274 N 275 = w * h * b, is_image x' 276 = len x'; 277} 278 279kurtosis x 280 = oo_unary_function kurtosis_op x, is_class x 281 = sum ((x - m) ** 4) / ((N - 1) * s ** 4), is_image x 282 = sum ((Group x' - m) ** 4).value / ((N - 1) * s ** 4), is_list x 283 = error (_ "bad arguments (" ++ print x ++ ") to " ++ "kurtosis") 284{ 285 kurtosis_op = Operator "kurtosis" kurtosis Operator_type.COMPOUND false; 286 287 // squash any large matrix down to a flat list ... much simpler 288 x' 289 = x, is_image x; 290 = flatten x; 291 292 m = mean x'; 293 s = deviation x'; 294 w = get_width x'; 295 h = get_height x'; 296 b = get_bands x'; 297 298 N 299 = len x', is_list x'; 300 = w * h * b; 301} 302 303// zero-excluding mean 304meanze x 305 = oo_unary_function meanze_op x, is_class x 306 = meanze_image_hist x, is_image x && 307 (fmt == Image_format.UCHAR || fmt == Image_format.USHORT) 308 = meanze_image x, is_image x 309 = meanze_list x, is_list x 310 = error (_ "bad arguments (" ++ print x ++ ") to " ++ "meanze") 311{ 312 fmt = get_format x; 313 314 meanze_op = Operator "meanze" meanze Operator_type.COMPOUND false; 315 316 meanze_list l = sum l / size l; 317 318 // number of non-zero elements in some sort of nested-list thing 319 size l 320 = foldr acc 0 l 321 { 322 acc x total 323 = total + size x, is_list x 324 = total + 1, x != 0; 325 = total; 326 } 327 328 // add elements in a nested-list thing 329 sum l 330 = foldr acc 0 l 331 { 332 acc x total 333 = total + sum x, is_list x 334 = total + x; 335 } 336 337 // image mean, for any image type 338 meanze_image i 339 = sum / N 340 { 341 w = get_width i; 342 h = get_height i; 343 b = get_bands i; 344 345 st = stats i; 346 sum = st.value?0?2; 347 348 // find non-zero pixels (not zero in all bands) 349 zp = im_notequal_vec i (replicate b 0); 350 351 // number of non-zero pixels 352 N = b * (mean zp * w * h) / 255; 353 } 354 355 // image mean for 8 and 16-bit unsigned images 356 // we can use a histogram, yay, and save a pass through the image 357 meanze_image_hist i 358 = sum / N 359 { 360 // histogram, knock out zeros 361 hist = hist_find i; 362 black = image_new 1 1 (get_bands hist) 363 (get_format hist) (get_coding hist) (get_type hist) 0 0 0; 364 histze = insert 0 0 black hist; 365 366 // matching identity 367 iden 368 = im_identity_ushort (get_bands hist) (get_width hist), 369 (get_width hist) > 256 370 = im_identity (get_bands hist); 371 372 // number of non-zero pixels 373 N = mean histze * 256; 374 375 // sum of pixels 376 sum = mean (hist * iden) * 256; 377 } 378} 379 380deviation x 381 = oo_unary_function deviation_op x, is_class x 382 = im_deviate x, is_image x 383 = deviation_list x, is_real_list x || is_matrix x 384 = error (_ "bad arguments to " ++ "deviation") 385{ 386 deviation_op = Operator "deviation" 387 deviation Operator_type.COMPOUND false; 388 389 deviation_list l 390 = (abs (s2 - (s * s / n)) / (n - 1)) ** 0.5 391 { 392 [n, s, s2] = sum_sum2_list l; 393 } 394 395 // return n, sum, sum of squares for a list of reals 396 sum_sum2_list x 397 = foldr accumulate [0, 0, 0] x 398 { 399 accumulate x sofar 400 = [n + 1, x + s, x * x + s2], is_real x 401 = [n + n', s + s', s2 + s2'], is_list x 402 = error "sum_sum2_list: not real or [real]" 403 { 404 [n, s, s2] = sofar; 405 [n', s', s2'] = sum_sum2_list x; 406 } 407 } 408} 409 410deviationze x 411 = oo_unary_function deviationze_op x, is_class x 412 = deviationze_image x, is_image x 413 = deviationze_list x, is_real_list x || is_matrix x 414 = error (_ "bad arguments to " ++ "deviationze") 415{ 416 deviationze_op = Operator "deviationze" 417 deviationze Operator_type.COMPOUND false; 418 419 deviationze_list l 420 = (abs (s2 - (s * s / n)) / (n - 1)) ** 0.5 421 { 422 [n, s, s2] = sum_sum2_list l; 423 } 424 425 // return number of non-zero elements, sum, sum of squares for a list of 426 // reals 427 sum_sum2_list x 428 = foldr accumulate [0, 0, 0] x 429 { 430 accumulate x sofar 431 = sofar, is_real x && x == 0 432 = [n + 1, x + s, x * x + s2], is_real x 433 = [n + n', s + s', s2 + s2'], is_list x 434 = error "sum_sum2_list: not real or [real]" 435 { 436 [n, s, s2] = sofar; 437 [n', s', s2'] = sum_sum2_list x; 438 } 439 } 440 441 deviationze_image i 442 = ((sum2 - sum * sum / N) / (N - 1)) ** 0.5 443 { 444 w = get_width i; 445 h = get_height i; 446 b = get_bands i; 447 448 st = stats i; 449 sum = st.value?0?2; 450 sum2 = st.value?0?3; 451 452 // find non-zero pixels (not zero in all bands) 453 zp = im_notequal_vec i (replicate b 0); 454 455 // number of non-zero pixels 456 N = b * (mean zp * w * h) / 255; 457 } 458} 459 460// find the centre of gravity of a histogram 461gravity x 462 = oo_unary_function gravity_op x, is_class x 463 = im_hist_gravity x, is_hist x 464 = gravity_list x, is_list x 465 = error (_ "bad arguments to " ++ "gravity") 466{ 467 gravity_op = Operator "gravity" gravity Operator_type.COMPOUND false; 468 469 // centre of gravity of a histogram... use the histogram to weight an 470 // identity, then sum, then find the mean element 471 im_hist_gravity h 472 = m 473 { 474 // make horizontal 475 h' 476 = rot270 h, get_width h == 1 477 = h, get_height h == 1 478 = error "width or height not 1"; 479 480 // number of elements 481 w = get_width h'; 482 483 // matching identity 484 i 485 = im_identity_ushort 1 w, w <= 2 ** 16 - 1 486 = make_xy w 1 ? 0; 487 488 // weight identity and sum 489 s = mean (i * h') * w; 490 491 // sum of original histogram 492 s' = mean h * w; 493 494 // weighted mean 495 m = s / s'; 496 } 497 498 gravity_list l 499 = m 500 { 501 w = len l; 502 503 // matching identity 504 i = [0, 1 .. w - 1]; 505 506 // weight identity and sum 507 s = sum (map2 multiply i l); 508 509 // sum of original histogram 510 s' = sum l; 511 512 // weighted mean 513 m = s / s'; 514 } 515} 516 517project x 518 = oo_unary_function project_op x, is_class x 519 = im_project x, is_image x 520 = error (_ "bad arguments to " ++ "project") 521{ 522 project_op = Operator "project" project Operator_type.COMPOUND false; 523} 524 525abs x 526 = oo_unary_function abs_op x, is_class x 527 = im_abs x, is_image x 528 = abs_cmplx x, is_complex x 529 = abs_num x, is_real x 530 = abs_list x, is_real_list x 531 = abs_list (map abs_list x), is_matrix x 532 = error (_ "bad arguments to " ++ "abs") 533{ 534 abs_op = Operator "abs" abs Operator_type.COMPOUND false; 535 536 abs_list l = (sum (map square l)) ** 0.5; 537 538 abs_num n 539 = n, n >= 0 540 = -n; 541 542 abs_cmplx c = ((re c)**2 + (im c)**2) ** 0.5; 543} 544 545copy x 546 = oo_unary_function copy_op x, is_class x 547 = im_copy x, is_image x 548 = x 549{ 550 copy_op = Operator "copy" copy Operator_type.COMPOUND_REWRAP false; 551} 552 553// like abs, but treat pixels as vectors ... ie. always get a 1-band image 554// back ... also treat matricies as lists of vectors 555// handy for dE from object difference 556abs_vec x 557 = oo_unary_function abs_vec_op x, is_class x 558 = abs_vec_image x, is_image x 559 = abs_vec_cmplx x, is_complex x 560 = abs_vec_num x, is_real x 561 = abs_vec_list x, is_real_list x 562 = mean (map abs_vec_list x), is_matrix x 563 = error (_ "bad arguments to " ++ "abs_vec") 564{ 565 abs_vec_op = Operator "abs_vec" 566 abs_vec Operator_type.COMPOUND false; 567 568 abs_vec_list l = (sum (map square l)) ** 0.5; 569 570 abs_vec_num n 571 = n, n >= 0 572 = -n; 573 574 abs_vec_cmplx c = ((re c)**2 + (im c)**2) ** 0.5; 575 576 abs_vec_image im 577 = (sum (map square (bandsplit im))) ** 0.5; 578} 579 580transpose x 581 = oo_unary_function transpose_op x, is_class x 582 = transpose_image x, is_image x 583 = transpose_list x, is_listof is_list x 584 = error (_ "bad arguments to " ++ "transpose") 585{ 586 transpose_op = Operator "transpose" 587 transpose Operator_type.COMPOUND_REWRAP false; 588 589 transpose_list l 590 = [], l' == [] 591 = (map hd l') : (transpose_list (map tl l')) 592 { 593 l' = takewhile (not_equal []) l; 594 } 595 596 transpose_image = im_flipver @ im_rot270; 597} 598 599rot45 x 600 = oo_unary_function rot45_op x, is_class x 601 = error "rot45 image: not implemented", is_image x 602 = error (_ "bad arguments to " ++ "rot45") 603{ 604 rot45_op = Operator "rot45" 605 rot45_object Operator_type.COMPOUND_REWRAP false; 606 607 rot45_object x 608 = rot45_matrix x, is_odd_square_matrix x 609 = error "rot45 image: not implemented", is_image x 610 = error (_ "bad arguments to " ++ "rot45"); 611 612 // slow, but what the heck 613 rot45_matrix l = (im_rotate_dmask45 (Matrix l)).value; 614} 615 616// apply an image function to a [[real]] ... matrix is converted to a 1 band 617// image for processing 618apply_matrix_as_image fn m 619 = (get_value @ im_vips2mask @ fn @ im_mask2vips @ Matrix) m; 620 621// a general image/matrix operation where the mat version is most easily done 622// by converting mat->image->mat 623apply_matim_operation name fn x 624 = oo_unary_function class_op x, is_class x 625 = fn x, is_image x 626 = apply_matrix_as_image fn x, is_matrix x 627 = error (_ "bad arguments to " ++ name) 628{ 629 class_op = Operator name 630 (apply_matim_operation name fn) Operator_type.COMPOUND_REWRAP false; 631} 632 633rot90 = apply_matim_operation "rot90" im_rot90; 634rot180 = apply_matim_operation "rot180" im_rot180; 635rot270 = apply_matim_operation "rot270" im_rot270; 636rotquad = apply_matim_operation "rotquad" im_rotquad; 637fliplr = apply_matim_operation "fliplr" im_fliphor; 638fliptb = apply_matim_operation "flipud" im_flipver; 639 640image_set_type type x 641 = oo_unary_function image_set_type_op x, is_class x 642 = im_copy_set x (to_real type) 643 (get_header "Xres" x) (get_header "Yres" x) 644 (get_header "Xoffset" x) (get_header "Yoffset" x), 645 is_image x 646 = error (_ "bad arguments to " ++ "image_set_type:" ++ 647 print type ++ " " ++ print x) 648{ 649 image_set_type_op = Operator "image_set_type" 650 (image_set_type type) Operator_type.COMPOUND_REWRAP false; 651} 652 653image_set_origin xoff yoff x 654 = oo_unary_function image_set_origin_op x, is_class x 655 = im_copy_set x 656 (get_header "Type" x) 657 (get_header "Xres" x) (get_header "Yres" x) 658 (to_real xoff) (to_real yoff), 659 is_image x 660 = error (_ "bad arguments to " ++ "image_set_origin") 661{ 662 image_set_origin_op = Operator "image_set_origin" 663 (image_set_origin xoff yoff) 664 Operator_type.COMPOUND_REWRAP false; 665} 666 667cache tile_width tile_height max_tiles x 668 = oo_unary_function cache_op x, is_class x 669 = im_tile_cache_random x (to_real tile_width) (to_real tile_height) 670 (to_real max_tiles), is_image x 671 = error (_ "bad arguments to " ++ "cache") 672{ 673 cache_op = Operator "cache" 674 (cache tile_width tile_height max_tiles) 675 Operator_type.COMPOUND_REWRAP false; 676} 677 678tile across down x 679 = oo_unary_function tile_op x, is_class x 680 = im_replicate x (to_real across) (to_real down), is_image x 681 = error (_ "bad arguments to " ++ "tile") 682{ 683 tile_op = Operator "tile" 684 (tile across down) Operator_type.COMPOUND_REWRAP false; 685} 686 687grid tile_height across down x 688 = oo_unary_function grid_op x, is_class x 689 = im_grid x (to_real tile_height) (to_real across) (to_real down), 690 is_image x 691 = error (_ "bad arguments to " ++ "grid") 692{ 693 grid_op = Operator "grid" 694 (grid tile_height across down) Operator_type.COMPOUND_REWRAP false; 695} 696 697max_pair a b 698 = a, a > b 699 = b; 700 701min_pair a b 702 = a, a < b 703 = b; 704 705range min value max = min_pair max (max_pair min value); 706 707max x 708 = oo_unary_function max_op x, is_class x 709 = im_max x, is_image x 710 = max_list x, is_list x 711 = x, is_number x 712 = error (_ "bad arguments to " ++ "max") 713{ 714 max_op = Operator "max" max Operator_type.COMPOUND false; 715 716 max_list x 717 = error "max []", x == [] 718 = foldr1 max_pair x, is_real_list x 719 = foldr1 max_pair (map max_list x), is_list x 720 = max x; 721} 722 723min x 724 = oo_unary_function min_op x, is_class x 725 = im_min x, is_image x 726 = min_list x, is_list x 727 = x, is_number x 728 = error (_ "bad arguments to " ++ "min") 729{ 730 min_op = Operator "min" min Operator_type.COMPOUND false; 731 732 min_list x 733 = error "min []", x == [] 734 = foldr1 min_pair x, is_real_list x 735 = foldr1 min_pair (map min_list x), is_list x 736 = min x; 737} 738 739maxpos x 740 = oo_unary_function maxpos_op x, is_class x 741 = im_maxpos x, is_image x 742 = maxpos_matrix x, is_matrix x 743 = maxpos_list x, is_list x 744 = error (_ "bad arguments to " ++ "maxpos") 745{ 746 maxpos_op = Operator "maxpos" maxpos Operator_type.COMPOUND false; 747 748 maxpos_matrix m 749 = (-1, -1), m == [[]] 750 = (indexes?row, row) 751 { 752 max_value = max (Matrix m); 753 indexes = map (index (equal max_value)) m; 754 row = index (not_equal (-1)) indexes; 755 } 756 757 maxpos_list l 758 = -1, l == [] 759 = index (equal (max l)) l; 760} 761 762minpos x 763 = oo_unary_function minpos_op x, is_class x 764 = im_minpos x, is_image x 765 = minpos_matrix x, is_matrix x 766 = minpos_list x, is_list x 767 = error (_ "bad arguments to " ++ "minpos") 768{ 769 minpos_op = Operator "minpos" minpos Operator_type.COMPOUND false; 770 771 minpos_matrix m 772 = (-1, -1), m == [[]] 773 = (indexes?row, row) 774 { 775 min_value = min (Matrix m); 776 indexes = map (index (equal min_value)) m; 777 row = index (not_equal (-1)) indexes; 778 } 779 780 minpos_list l 781 = -1, l == [] 782 = index (equal (min l)) l; 783} 784 785stats x 786 = oo_unary_function stats_op x, is_class x 787 = im_stats x, is_image x 788 = im_stats (to_image x).value, is_matrix x 789 = error (_ "bad arguments to " ++ "stats") 790{ 791 stats_op = Operator "stats" 792 stats Operator_type.COMPOUND false; 793} 794 795e = 2.7182818284590452354; 796 797pi = 3.14159265358979323846; 798 799rad d = 2 * pi * (d / 360); 800 801deg r = 360 * r / (2 * pi); 802 803sign x 804 = oo_unary_function sign_op x, is_class x 805 = im_sign x, is_image x 806 = sign_cmplx x, is_complex x 807 = sign_num x, is_real x 808 = error (_ "bad arguments to " ++ "sign") 809{ 810 sign_op = Operator "sign" sign Operator_type.COMPOUND_REWRAP false; 811 812 sign_num n 813 = 0, n == 0 814 = 1, n > 0 815 = -1; 816 817 sign_cmplx c 818 = (0, 0), mod == 0 819 = (re c / mod, im c / mod) 820 { 821 mod = abs c; 822 } 823} 824 825rint x 826 = oo_unary_function rint_op x, is_class x 827 = im_rint x, is_image x 828 = rint_value x, is_number x 829 = error (_ "bad arguments to " ++ "rint") 830{ 831 rint_op = Operator "rint" rint Operator_type.ARITHMETIC false; 832 833 rint_value x 834 = (int) (x + 0.5), x > 0 835 = (int) (x - 0.5); 836} 837 838scale x 839 = oo_unary_function scale_op x, is_class x 840 = (unsigned char) x, is_number x 841 = im_scale x, is_image x 842 = scale_list x, is_real_list x || is_matrix x 843 = error (_ "bad arguments to " ++ "scale") 844{ 845 scale_op = Operator "scale" scale Operator_type.COMPOUND_REWRAP false; 846 847 scale_list l 848 = apply_scale s o l 849 { 850 mn = find_limit min_pair l; 851 mx = find_limit max_pair l; 852 s = 255.0 / (mx - mn); 853 o = -(mn * s); 854 } 855 856 find_limit fn l 857 = find_limit fn (map (find_limit fn) l), is_listof is_list l 858 = foldr1 fn l; 859 860 apply_scale s o x 861 = x * s + o, is_number x 862 = map (apply_scale s o) x; 863} 864 865scaleps x 866 = oo_unary_function scale_op x, is_class x 867 = im_scaleps x, is_image x 868 = error (_ "bad arguments to " ++ "scale") 869{ 870 scale_op = Operator "scaleps" 871 scaleps Operator_type.COMPOUND_REWRAP false; 872} 873 874fwfft x 875 = oo_unary_function fwfft_op x, is_class x 876 = im_fwfft x, is_image x 877 = error (_ "bad arguments to " ++ "fwfft") 878{ 879 fwfft_op = Operator "fwfft" 880 fwfft Operator_type.COMPOUND_REWRAP false; 881} 882 883invfft x 884 = oo_unary_function invfft_op x, is_class x 885 = im_invfftr x, is_image x 886 = error (_ "bad arguments to " ++ "invfft") 887{ 888 invfft_op = Operator "invfft" 889 invfft Operator_type.COMPOUND_REWRAP false; 890} 891 892falsecolour x 893 = oo_unary_function falsecolour_op x, is_class x 894 = image_set_type Image_type.sRGB (im_falsecolour x), is_image x 895 = error (_ "bad arguments to " ++ "falsecolour") 896{ 897 falsecolour_op = Operator "falsecolour" 898 falsecolour Operator_type.COMPOUND_REWRAP false; 899} 900 901polar x 902 = oo_unary_function polar_op x, is_class x 903 = im_c2amph x, is_image x 904 = polar_cmplx x, is_complex x 905 = error (_ "bad arguments to " ++ "polar") 906{ 907 polar_op = Operator "polar" polar Operator_type.COMPOUND false; 908 909 polar_cmplx r 910 = (l, a) 911 { 912 a 913 = 270, x == 0 && y < 0 914 = 90, x == 0 && y >= 0 915 = 360 + atan (y / x), x > 0 && y < 0 916 = atan (y / x), x > 0 && y >= 0 917 = 180 + atan (y / x); 918 919 l = (x ** 2 + y ** 2) ** 0.5; 920 921 x = re r; 922 y = im r; 923 } 924} 925 926rectangular x 927 = oo_unary_function rectangular_op x, is_class x 928 = im_c2rect x, is_image x 929 = rectangular_cmplx x, is_complex x 930 = error (_ "bad arguments to " ++ "rectangular") 931{ 932 rectangular_op = Operator "rectangular" 933 rectangular Operator_type.COMPOUND false; 934 935 rectangular_cmplx p 936 = (x, y) 937 { 938 l = re p; 939 a = im p; 940 941 x = l * cos a; 942 y = l * sin a; 943 } 944} 945 946// we can't use colour_unary: that likes 3 band only 947recomb matrix x 948 = oo_unary_function recomb_op x, is_class x 949 = im_recomb x matrix, is_image x 950 = recomb_real_list x, is_real_list x 951 = map recomb_real_list x, is_matrix x 952 = error (_ "bad arguments to " ++ "recomb") 953{ 954 // COMPOUND_REWRAP ... signal to the colour class to go to image and 955 // back 956 recomb_op = Operator "recomb" 957 (recomb matrix) Operator_type.COMPOUND_REWRAP false; 958 959 // process [1,2,3 ..] as an image 960 recomb_real_list l 961 = (to_matrix im').value?0 962 { 963 im = (float) (to_image (Vector l)).value; 964 im' = recomb matrix im; 965 } 966} 967 968extract_area x y w h obj 969 = oo_unary_function extract_area_op obj, is_class obj 970 = im_extract_area obj x' y' w' h', is_image obj 971 = map (extract_range x' w') (extract_range y' h' obj), is_matrix obj 972 = error (_ "bad arguments to " ++ "extract_area") 973{ 974 x' = to_real x; 975 y' = to_real y; 976 w' = to_real w; 977 h' = to_real h; 978 979 extract_area_op = Operator "extract_area" (extract_area x y w h) 980 Operator_type.COMPOUND_REWRAP false; 981 982 extract_range from length list 983 = (take length @ drop from) list; 984} 985 986extract_band b obj = subscript obj b; 987 988extract_row y obj 989 = oo_unary_function extract_row_op obj, is_class obj 990 = extract_area 0 y' (get_width obj) 1 obj, is_image obj 991 = [obj?y'], is_matrix obj 992 = error (_ "bad arguments to " ++ "extract_row") 993{ 994 y' = to_real y; 995 996 extract_row_op = Operator "extract_row" (extract_row y) 997 Operator_type.COMPOUND_REWRAP false; 998} 999 1000extract_column x obj 1001 = oo_unary_function extract_column_op obj, is_class obj 1002 = extract_area x' 0 1 height obj, is_image obj 1003 = map (\row [row?x']) obj, is_matrix obj 1004 = error (_ "bad arguments to " ++ "extract_column") 1005{ 1006 x' = to_real x; 1007 height = get_header "Ysize" obj; 1008 1009 extract_column_op = Operator "extract_column" (extract_column x) 1010 Operator_type.COMPOUND_REWRAP false; 1011} 1012 1013blend cond in1 in2 1014 = oo_binary_function blend_op cond [in1,in2], is_class cond 1015 = im_blend (get_image cond) (get_image in1) (get_image in2), 1016 has_image cond && has_image in1 && has_image in2 1017 = error (_ "bad arguments to " ++ "blend" ++ ": " ++ 1018 join_sep ", " (map print [cond, in1, in2])) 1019{ 1020 blend_op = Operator "blend" 1021 blend_obj Operator_type.COMPOUND_REWRAP false; 1022 1023 blend_obj cond x 1024 = blend_result_image 1025 { 1026 [then_part, else_part] = x; 1027 1028 // get things about our output from inputs in this order 1029 objects = [then_part, else_part, cond]; 1030 1031 // properties of our output image 1032 target_width = get_member_list has_width get_width objects; 1033 target_height = get_member_list has_height get_height objects; 1034 target_bands = get_member_list has_bands get_bands objects; 1035 target_format = get_member_list has_format get_format objects; 1036 target_type = get_member_list has_type get_type objects; 1037 1038 to_image x 1039 = to_image_size target_width target_height 1040 target_bands target_format x; 1041 1042 [then_image, else_image] = map to_image [then_part, else_part]; 1043 1044 blend_result_image = image_set_type target_type 1045 (im_blend cond then_image else_image); 1046 } 1047} 1048 1049// do big first: we want to keep big's class, if possible 1050// eg. big is a Plot, small is a 1x1 Image 1051insert x y small big 1052 = oo_binary'_function insert_op small big, is_class big 1053 = oo_binary_function insert_op small big, is_class small 1054 = im_insert big small (to_real x) (to_real y), 1055 is_image small && is_image big 1056 = error (_ "bad arguments to " ++ "insert") 1057{ 1058 insert_op = Operator "insert" 1059 (insert x y) Operator_type.COMPOUND_REWRAP false; 1060} 1061 1062insert_noexpand x y small big 1063 = oo_binary_function insert_noexpand_op small big, is_class small 1064 = oo_binary'_function insert_noexpand_op small big, is_class big 1065 = im_insert_noexpand big small (to_real x) (to_real y), 1066 is_image small && is_image big 1067 = error (_ "bad arguments to " ++ "insert_noexpand") 1068{ 1069 insert_noexpand_op = Operator "insert_noexpand" 1070 (insert_noexpand x y) Operator_type.COMPOUND_REWRAP false; 1071} 1072 1073decode im 1074 = oo_unary_function decode_op im, is_class im 1075 = decode_im im, is_image im 1076 = error (_ "bad arguments to " ++ "decode") 1077{ 1078 decode_op = Operator "decode" 1079 decode Operator_type.COMPOUND_REWRAP false; 1080 1081 decode_im im 1082 = im_LabQ2Lab im, get_coding im == Image_coding.LABPACK 1083 = im_rad2float im, get_coding im == Image_coding.RAD 1084 = im; 1085} 1086 1087measure_draw across down measure image 1088 = mark 1089{ 1090 patch_width = image.width / across; 1091 sample_width = patch_width * (measure / 100); 1092 left_margin = (patch_width - sample_width) / 2; 1093 patch_height = image.height / down; 1094 sample_height = patch_height * (measure / 100); 1095 top_margin = (patch_height - sample_height) / 2; 1096 1097 cods = [[x * patch_width + left_margin, y * patch_height + top_margin] :: 1098 y <- [0 .. down - 1]; x <- [0 .. across - 1]]; 1099 x = map (extract 0) cods; 1100 y = map (extract 1) cods; 1101 1102 outer = mkim [$pixel => 255] sample_width sample_height 1; 1103 inner = mkim [] (sample_width - 4) (sample_height - 4) 1; 1104 patch = insert 2 2 inner outer; 1105 1106 bg = mkim [] image.width image.height 1; 1107 1108 mask = Image (im_insertset bg.value patch.value x y); 1109 1110 image' = colour_transform_to Image_type.sRGB image; 1111 1112 mark = if mask then Vector [0, 255, 0] else image'; 1113} 1114 1115measure_sample across down measure image 1116 = measures 1117{ 1118 patch_width = image.width / across; 1119 sample_width = patch_width * (measure / 100); 1120 left_margin = (patch_width - sample_width) / 2; 1121 patch_height = image.height / down; 1122 sample_height = patch_height * (measure / 100); 1123 top_margin = (patch_height - sample_height) / 2; 1124 1125 cods = [[x * patch_width + left_margin, y * patch_height + top_margin] :: 1126 y <- [0 .. down - 1]; x <- [0 .. across - 1]]; 1127 1128 image' = decode image; 1129 patches = map (\p extract_area p?0 p?1 sample_width sample_height image') 1130 cods; 1131 measures = Matrix (map (map mean) (map bandsplit patches)); 1132} 1133 1134extract_bands b n obj 1135 = oo_unary_function extract_bands_op obj, is_class obj 1136 = im_extract_bands obj (to_real b) (to_real n), is_image obj 1137 = error (_ "bad arguments to " ++ "extract_bands") 1138{ 1139 extract_bands_op = Operator "extract_bands" 1140 (extract_bands b n) Operator_type.COMPOUND_REWRAP false; 1141} 1142 1143invert x 1144 = oo_unary_function invert_op x, is_class x 1145 = im_invert x, is_image x 1146 = 255 - x, is_real x 1147 = error (_ "bad arguments to " ++ "invert") 1148{ 1149 invert_op = Operator "invert" invert Operator_type.COMPOUND false; 1150} 1151 1152transform ipol wrap params image 1153 = oo_unary_function transform_op image, is_class image 1154 = im_transform image 1155 (to_matrix params) (to_real ipol) (to_real wrap), is_image image 1156 = error (_ "bad arguments to " ++ "transform") 1157{ 1158 transform_op = Operator "transform" 1159 (transform ipol wrap params) 1160 Operator_type.COMPOUND_REWRAP false; 1161} 1162 1163transform_search max_error max_iterations order ipol wrap sample reference 1164 = oo_binary_function transform_search_op sample reference, is_class sample 1165 = oo_binary'_function transform_search_op sample reference, 1166 is_class reference 1167 = im_transform_search sample reference 1168 (to_real max_error) (to_real max_iterations) (to_real order) 1169 (to_real ipol) (to_real wrap), 1170 is_image sample && is_image reference 1171 = error (_ "bad arguments to " ++ "transform_search") 1172{ 1173 transform_search_op = Operator "transform_search" 1174 (transform_search max_error max_iterations order ipol wrap) 1175 Operator_type.COMPOUND false; 1176} 1177 1178rotate interp angle image 1179 = oo_binary_function rotate_op angle image, is_class angle 1180 = oo_binary'_function rotate_op angle image, is_class image 1181 = im_affinei_all image interp.value a (-b) b a 0 0, 1182 is_real angle && is_image image 1183 = error (_ "bad arguments to " ++ "rotate") 1184{ 1185 rotate_op = Operator "rotate" 1186 (rotate interp) Operator_type.COMPOUND_REWRAP false; 1187 a = cos angle; 1188 b = sin angle; 1189} 1190 1191matrix_binary fn a b 1192 = itom (fn (mtoi a) (mtoi b)) 1193{ 1194 mtoi x = im_mask2vips (Matrix x); 1195 itom x = (im_vips2mask x).value; 1196} 1197 1198join_lr a b 1199 = oo_binary_function join_lr_op a b, is_class a 1200 = oo_binary'_function join_lr_op a b, is_class b 1201 = join_im a b, is_image a && is_image b 1202 = matrix_binary join_im a b, is_matrix a && is_matrix b 1203 = error (_ "bad arguments to " ++ "join_lr") 1204{ 1205 join_lr_op = Operator "join_lr" 1206 join_lr Operator_type.COMPOUND_REWRAP false; 1207 1208 join_im a b = insert (get_width a) 0 b a; 1209} 1210 1211join_tb a b 1212 = oo_binary_function join_tb_op a b, is_class a 1213 = oo_binary'_function join_tb_op a b, is_class b 1214 = join_im a b, is_image a && is_image b 1215 = matrix_binary join_im a b, is_matrix a && is_matrix b 1216 = error (_ "bad arguments to " ++ "join_tb") 1217{ 1218 join_tb_op = Operator "join_tb" 1219 join_tb Operator_type.COMPOUND_REWRAP false; 1220 1221 join_im a b = insert 0 (get_height a) b a; 1222} 1223 1224conj x 1225 = oo_unary_function conj_op x, is_class x 1226 = (re x, -im x), 1227 is_complex x || 1228 (is_image x && format == Image_format.COMPLEX) || 1229 (is_image x && format == Image_format.DPCOMPLEX) 1230 // assume it's some sort of real 1231 = x 1232{ 1233 format = get_header "BandFmt" x; 1234 conj_op = Operator "conj" conj Operator_type.COMPOUND false; 1235} 1236 1237clip2fmt format image 1238 = oo_unary_function clip2fmt_op image, is_class image 1239 = im_clip2fmt image (to_real format), is_image image 1240 = error (_ "bad arguments to " ++ "clip2fmt") 1241{ 1242 clip2fmt_op = Operator "clip2fmt" 1243 (clip2fmt format) Operator_type.COMPOUND_REWRAP false; 1244} 1245 1246embed type x y w h im 1247 = oo_unary_function embed_op im, is_class im 1248 = im_embed im (to_real type) 1249 (to_real x) (to_real y) (to_real w) (to_real h), is_image im 1250 = error (_ "bad arguments to " ++ "embed") 1251{ 1252 embed_op = Operator "embed" 1253 (embed type x y w h) Operator_type.COMPOUND_REWRAP false; 1254} 1255 1256/* Morph a mask with a [[real]] matrix ... turn m2 into an image, morph it 1257 * with m1, turn it back to a matrix again. 1258 */ 1259_morph_2_masks fn m1 m2 1260 = m'' 1261{ 1262 // The [[real]] can contain 128 values ... squeeze them out! 1263 image = im_mask2vips (Matrix m2) == 255; 1264 m2_width = get_width image; 1265 m2_height = get_height image; 1266 1267 // need to embed m2 in an image large enough for us to be able to 1268 // position m1 all around the edges, with a 1 pixel overlap 1269 image' = embed 0 1270 (m1.width / 2) (m1.height / 2) 1271 (m2_width + (m1.width - 1)) (m2_height + (m1.height - 1)) 1272 image; 1273 1274 // morph! 1275 image'' = fn m1 image'; 1276 1277 // back to mask 1278 m' = im_vips2mask ((double) image''); 1279 1280 // Turn 0 in output to 128 (don't care). 1281 m'' 1282 = map (map fn) m'.value 1283 { 1284 fn a 1285 = 128, a == 0; 1286 = a; 1287 } 1288} 1289 1290dilate mask image 1291 = oo_unary_function dilate_op image, is_class image 1292 = im_dilate image (to_matrix mask), is_image image 1293 = error (_ "bad arguments to " ++ "dilate") 1294{ 1295 dilate_op = Operator "dilate" 1296 dilate_object Operator_type.COMPOUND_REWRAP false; 1297 1298 dilate_object x 1299 = _morph_2_masks dilate mask x, is_matrix x 1300 = dilate mask x; 1301} 1302 1303erode mask image 1304 = oo_unary_function erode_op image, is_class image 1305 = im_erode image (to_matrix mask), is_image image 1306 = error (_ "bad arguments to " ++ "erode") 1307{ 1308 erode_op = Operator "erode" 1309 erode_object Operator_type.COMPOUND_REWRAP false; 1310 1311 erode_object x 1312 = _morph_2_masks erode mask x, is_matrix x 1313 = erode mask x; 1314} 1315 1316conv mask image 1317 = oo_unary_function conv_op image, is_class image 1318 = im_conv image (to_matrix mask), is_image image 1319 = error (_ "bad arguments to " ++ "conv" ++ ": " ++ 1320 "conv (" ++ print mask ++ ") (" ++ print image ++ ")") 1321{ 1322 conv_op = Operator "conv" 1323 (conv mask) Operator_type.COMPOUND_REWRAP false; 1324} 1325 1326convf mask image 1327 = oo_unary_function convf_op image, is_class image 1328 = im_conv_f image (to_matrix mask), is_image image 1329 = error (_ "bad arguments to " ++ "convf" ++ ": " ++ 1330 "convf (" ++ print mask ++ ") (" ++ print image ++ ")") 1331{ 1332 convf_op = Operator "convf" 1333 (convf mask) Operator_type.COMPOUND_REWRAP false; 1334} 1335 1336convsep mask image 1337 = oo_unary_function convsep_op image, is_class image 1338 = im_convsep image (to_matrix mask), is_image image 1339 = error (_ "bad arguments to " ++ "convsep") 1340{ 1341 convsep_op = Operator "convsep" 1342 (convsep mask) Operator_type.COMPOUND_REWRAP false; 1343} 1344 1345aconvsep layers mask image 1346 = oo_unary_function aconvsep_op image, is_class image 1347 = im_aconvsep image (to_matrix mask) (to_real layers), is_image image 1348 = error (_ "bad arguments to " ++ "aconvsep") 1349{ 1350 aconvsep_op = Operator "aconvsep" 1351 (aconvsep layers mask) Operator_type.COMPOUND_REWRAP false; 1352} 1353 1354convsepf mask image 1355 = oo_unary_function convsepf_op image, is_class image 1356 = im_convsep_f image (to_matrix mask), is_image image 1357 = error (_ "bad arguments to " ++ "convsepf") 1358{ 1359 convsepf_op = Operator "convsepf" 1360 (convsepf mask) Operator_type.COMPOUND_REWRAP false; 1361} 1362 1363rank w h n image 1364 = oo_unary_function rank_op image, is_class image 1365 = im_rank image (to_real w) (to_real h) (to_real n), is_image image 1366 = error (_ "bad arguments to " ++ "rank") 1367{ 1368 rank_op = Operator "rank" 1369 (rank w h n) Operator_type.COMPOUND_REWRAP false; 1370} 1371 1372rank_image n x 1373 = rlist x.value, is_Group x 1374 = rlist x, is_list x 1375 = error (_ "bad arguments to " ++ "rank_image") 1376{ 1377 rlist l 1378 = wrapper ranked, has_wrapper 1379 = ranked 1380 { 1381 has_wrapper = has_member_list (has_member "Image") l; 1382 wrapper = get_member_list (has_member "Image") (get_member "Image") l; 1383 ranked = im_rank_image (map get_image l) (to_real n); 1384 } 1385} 1386 1387// find the correlation surface for a small image within a big one 1388correlate small big 1389 = oo_binary_function correlate_op small big, is_class small 1390 = oo_binary'_function correlate_op small big, is_class big 1391 = im_spcor big small, is_image small && is_image big 1392 = error (_ "bad arguments to " ++ "correlate") 1393{ 1394 correlate_op = Operator "correlate" 1395 correlate Operator_type.COMPOUND_REWRAP false; 1396} 1397 1398// just sum-of-squares-of-differences 1399correlate_fast small big 1400 = oo_binary_function correlate_fast_op small big, is_class small 1401 = oo_binary'_function correlate_fast_op small big, is_class big 1402 = im_fastcor big small, is_image small && is_image big 1403 = error (_ "bad arguments to " ++ "correlate_fast") 1404{ 1405 correlate_fast_op = Operator "correlate_fast" 1406 correlate_fast Operator_type.COMPOUND_REWRAP false; 1407} 1408 1409// set Type, wrap as Plot_hist if it's a class 1410hist_tag x 1411 = oo_unary_function hist_tag_op x, is_class x 1412 = image_set_type Image_type.HISTOGRAM x, is_image x 1413 = error (_ "bad arguments to " ++ "hist_tag") 1414{ 1415 hist_tag_op = Operator "hist_tag" 1416 (Plot_histogram @ hist_tag) Operator_type.COMPOUND false; 1417} 1418 1419hist_find x 1420 = oo_unary_function hist_find_op x, is_class x 1421 = im_histgr x (-1), is_image x 1422 = error (_ "bad arguments to " ++ "hist_find") 1423{ 1424 hist_find_op = Operator "hist_find" 1425 (Plot_histogram @ hist_find) Operator_type.COMPOUND false; 1426} 1427 1428hist_find_nD bins image 1429 = oo_unary_function hist_find_nD_op image, is_class image 1430 = im_histnD image (to_real bins), is_image image 1431 = error (_ "bad arguments to " ++ "hist_find_nD") 1432{ 1433 hist_find_nD_op = Operator "hist_find_nD" 1434 (hist_find_nD bins) Operator_type.COMPOUND_REWRAP false; 1435} 1436 1437hist_find_indexed mode index value 1438 = oo_binary_function hist_find_indexed_op index value, is_class index 1439 = oo_binary'_function hist_find_indexed_op index value, is_class value 1440 = indexed index value, is_image index && is_image value 1441 = error (_ "bad arguments to " ++ "hist_find_indexed") 1442{ 1443 hist_find_indexed_op = Operator "hist_find_indexed" 1444 (compose (compose Plot_histogram) (hist_find_indexed mode)) 1445 Operator_type.COMPOUND false; 1446 1447 indexed index value 1448 = out 1449 { 1450 [out] = vips_call "hist_find_indexed" [value, index] [ 1451 "combine" => mode 1452 ]; 1453 } 1454} 1455 1456hist_entropy x 1457 = oo_unary_function hist_entropy_op x, is_class x 1458 = entropy x, is_image x 1459 = error (_ "bad arguments to " ++ "hist_entropy") 1460{ 1461 hist_entropy_op = Operator "hist_entropy" 1462 hist_entropy Operator_type.COMPOUND_REWRAP false; 1463 1464 entropy x 1465 = out 1466 { 1467 [out] = vips_call "hist_entropy" [x] [ 1468 ]; 1469 } 1470} 1471 1472hist_map hist image 1473 = oo_binary_function hist_map_op hist image, is_class hist 1474 = oo_binary'_function hist_map_op hist image, is_class image 1475 = im_maplut image hist, is_image hist && is_image image 1476 = error (_ "bad arguments to " ++ "hist_map") 1477{ 1478 // can't use rewrap, as we want to always wrap as image 1479 hist_map_op = Operator "hist_map" 1480 (compose (compose Image) hist_map) Operator_type.COMPOUND false; 1481} 1482 1483hist_cum hist 1484 = oo_unary_function hist_cum_op hist, is_class hist 1485 = im_histcum hist, is_image hist 1486 = error (_ "bad arguments to " ++ "hist_cum") 1487{ 1488 hist_cum_op = Operator "hist_cum" 1489 hist_cum Operator_type.COMPOUND_REWRAP false; 1490} 1491 1492hist_diff hist 1493 = oo_unary_function hist_diff_op hist, is_class hist 1494 = im_histdiff hist, is_image hist 1495 = error (_ "bad arguments to " ++ "hist_diff") 1496{ 1497 hist_diff_op = Operator "hist_diff" 1498 hist_diff Operator_type.COMPOUND_REWRAP false; 1499 1500 im_histdiff h 1501 = (conv (Matrix [[-1, 1]]) @ clip2fmt (fmt (get_format h))) h 1502 { 1503 // up the format so it can represent the whole range of 1504 // possible values from this mask 1505 fmt x 1506 = Image_format.SHORT, 1507 x == Image_format.UCHAR || x == Image_format.CHAR 1508 = Image_format.INT, 1509 x == Image_format.USHORT || x == Image_format.SHORT || 1510 x == Image_format.UINT 1511 = x; 1512 } 1513} 1514 1515hist_norm hist 1516 = oo_unary_function hist_norm_op hist, is_class hist 1517 = im_histnorm hist, is_image hist 1518 = error (_ "bad arguments to " ++ "hist_norm") 1519{ 1520 hist_norm_op = Operator "hist_norm" 1521 hist_norm Operator_type.COMPOUND_REWRAP false; 1522} 1523 1524hist_inv hist 1525 = oo_unary_function hist_inv_op hist, is_class hist 1526 = inv hist, is_image hist 1527 = error (_ "bad arguments to " ++ "hist_inv") 1528{ 1529 hist_inv_op = Operator "hist_inv" 1530 hist_inv Operator_type.COMPOUND_REWRAP false; 1531 1532 inv im 1533 = im_invertlut (to_matrix im''') len 1534 { 1535 // need a vertical doublemask 1536 im' 1537 = rot90 im, get_width im > 1 && get_height im == 1 1538 = im, get_width im == 1 && get_height im > 1 1539 = error (_ "not a hist"); 1540 len = get_height im'; 1541 1542 // values must be scaled to 0 - 1 1543 im'' = im' / (max im'); 1544 1545 // add an index column on the left 1546 // again, must be in 0-1 1547 y = ((make_xy 1 len)?1) / len; 1548 im''' = y ++ im''; 1549 } 1550} 1551 1552hist_match in ref 1553 = oo_binary_function hist_match_op in ref, is_class in 1554 = oo_binary'_function hist_match_op in ref, is_class ref 1555 = im_histspec in ref, is_image in && is_image ref 1556 = error (_ "bad arguments to " ++ "hist_match") 1557{ 1558 hist_match_op = Operator "hist_match" 1559 hist_match Operator_type.COMPOUND_REWRAP false; 1560} 1561 1562hist_equalize x = hist_map ((hist_norm @ hist_cum @ hist_find) x) x; 1563 1564hist_equalize_local w h l image 1565 = oo_unary_function hist_equalize_local_op image, is_class image 1566 = out, is_image image 1567 = error (_ "bad arguments to " ++ "hist_equalize_local") 1568{ 1569 hist_equalize_local_op = Operator "hist_equalize_local" 1570 (hist_equalize_local w h l) Operator_type.COMPOUND_REWRAP false; 1571 1572 [out] = vips_call "hist_local" 1573 [image, to_real w, to_real h] [$max_slope => to_real l]; 1574} 1575 1576// find the threshold below which are percent of the image (percent in [0,1]) 1577// eg. hist_thresh 0.1 x == 12, then x < 12 will light up 10% of the pixels 1578hist_thresh percent image 1579 = x 1580{ 1581 // our own normaliser ... we don't want to norm channels separately 1582 // norm to [0,1] 1583 my_hist_norm h = h / max h; 1584 1585 // normalised cumulative hist 1586 // we sum the channels before we normalise, because we want to treat them 1587 // all the same 1588 h = (my_hist_norm @ sum @ bandsplit @ hist_cum @ hist_find) 1589 image.value; 1590 1591 // threshold that, then use im_profile to search for the x position in the 1592 // histogram 1593 x = mean (im_profile (h > percent) 1); 1594} 1595 1596/* Sometimes useful, despite plotting now being built in, for making 1597 * diagnostic images. 1598 */ 1599hist_plot hist 1600 = oo_unary_function hist_plot_op hist, is_class hist 1601 = im_histplot hist, is_image hist 1602 = error (_ "bad arguments to " ++ "hist_plot") 1603{ 1604 hist_plot_op = Operator "hist_plot" 1605 (Image @ hist_plot) Operator_type.COMPOUND false; 1606} 1607 1608zerox d x 1609 = oo_unary_function zerox_op x, is_class x 1610 = im_zerox x (to_real d), is_image x 1611 = error (_ "bad arguments to " ++ "zerox") 1612{ 1613 zerox_op = Operator "zerox" 1614 (zerox d) Operator_type.COMPOUND_REWRAP false; 1615} 1616 1617reduce kernel xshr yshr image 1618 = oo_unary_function reduce_op image, is_class image 1619 = reduce_im image, is_image image 1620 = error (_ "bad arguments to " ++ "reduce") 1621{ 1622 reduce_op = Operator "reduce" 1623 reduce_im Operator_type.COMPOUND_REWRAP false; 1624 1625 reduce_im im 1626 = out 1627 { 1628 [out] = vips_call "reduce" [im, xshr, yshr] [$kernel => kernel.value]; 1629 } 1630} 1631 1632similarity interpolate scale angle image 1633 = oo_unary_function similarity_op image, is_class image 1634 = similarity_im image, is_image image 1635 = error (_ "bad arguments to " ++ "similarity") 1636{ 1637 similarity_op = Operator "similarity" 1638 similarity_im Operator_type.COMPOUND_REWRAP false; 1639 1640 similarity_im im 1641 = out 1642 { 1643 [out] = vips_call "similarity" [im] [ 1644 $interpolate => interpolate.value, 1645 $scale => scale, 1646 $angle => angle 1647 ]; 1648 } 1649} 1650 1651resize kernel xfac yfac image 1652 = oo_unary_function resize_op image, is_class image 1653 = resize_im image, is_image image 1654 = error (_ "bad arguments to " ++ "resize") 1655{ 1656 resize_op = Operator "resize" 1657 resize_im Operator_type.COMPOUND_REWRAP false; 1658 1659 xfac' = to_real xfac; 1660 yfac' = to_real yfac; 1661 1662 rxfac' = 1 / xfac'; 1663 ryfac' = 1 / yfac'; 1664 1665 is_nn = kernel.type == Kernel_type.NEAREST_NEIGHBOUR; 1666 1667 resize_im im 1668 // upscale by integer factor, nearest neighbour 1669 = im_zoom im xfac' yfac', 1670 is_int xfac' && is_int yfac' && 1671 xfac' >= 1 && yfac' >= 1 && 1672 is_nn 1673 1674 // downscale by integer factor, nearest neighbour 1675 = im_subsample im rxfac' ryfac', 1676 is_int rxfac' && is_int ryfac' && 1677 rxfac' >= 1 && ryfac' >= 1 && 1678 is_nn 1679 1680 // everything else ... we just pass on to vips_resize() 1681 = vips_resize kernel xfac' yfac' im 1682 { 1683 vips_resize kernel hscale vscale im 1684 = out 1685 { 1686 [out] = vips_call "resize" [im, hscale] 1687 [$vscale => vscale, $kernel => kernel.value]; 1688 } 1689 } 1690} 1691 1692sharpen radius x1 y2 y3 m1 m2 in 1693 = oo_unary_function sharpen_op in, is_class in 1694 = im_sharpen in (to_real radius) 1695 (to_real x1) (to_real y2) (to_real y3) 1696 (to_real m1) (to_real m2), is_image in 1697 = error (_ "bad arguments to " ++ "sharpen") 1698{ 1699 sharpen_op = Operator "sharpen" 1700 (sharpen radius x1 y2 y3 m1 m2) 1701 Operator_type.COMPOUND_REWRAP false; 1702} 1703 1704tone_analyse s m h sa ma ha in 1705 = oo_unary_function tone_analyse_op in, is_class in 1706 = im_tone_analyse in 1707 (to_real s) (to_real m) (to_real h) 1708 (to_real sa) (to_real ma) (to_real ha), is_image in 1709 = error (_ "bad arguments to " ++ "tone_analyse") 1710{ 1711 tone_analyse_op = Operator "tone_analyse" 1712 (Plot_histogram @ tone_analyse s m h sa ma ha) 1713 Operator_type.COMPOUND false; 1714} 1715 1716tone_map hist image 1717 = oo_binary_function tone_map_op hist image, is_class hist 1718 = oo_binary'_function tone_map_op hist image, is_class image 1719 = im_tone_map image hist, is_image hist && is_image image 1720 = error (_ "bad arguments to " ++ "tone_map") 1721{ 1722 tone_map_op = Operator "tone_map" 1723 tone_map Operator_type.COMPOUND_REWRAP false; 1724} 1725 1726tone_build fmt b w s m h sa ma ha 1727 = (Plot_histogram @ clip2fmt fmt) 1728 (im_tone_build_range mx mx 1729 (to_real b) (to_real w) 1730 (to_real s) (to_real m) (to_real h) 1731 (to_real sa) (to_real ma) (to_real ha)) 1732{ 1733 mx = Image_format.maxval fmt; 1734} 1735 1736icc_export depth profile intent in 1737 = oo_unary_function icc_export_op in, is_class in 1738 = im_icc_export_depth in 1739 (to_real depth) (expand profile) (to_real intent), is_image in 1740 = error (_ "bad arguments to " ++ "icc_export") 1741{ 1742 icc_export_op = Operator "icc_export" 1743 (icc_export depth profile intent) 1744 Operator_type.COMPOUND_REWRAP false; 1745} 1746 1747icc_import profile intent in 1748 = oo_unary_function icc_import_op in, is_class in 1749 = im_icc_import in 1750 (expand profile) (to_real intent), is_image in 1751 = error (_ "bad arguments to " ++ "icc_import") 1752{ 1753 icc_import_op = Operator "icc_import" 1754 (icc_import profile intent) 1755 Operator_type.COMPOUND_REWRAP false; 1756} 1757 1758icc_import_embedded intent in 1759 = oo_unary_function icc_import_embedded_op in, is_class in 1760 = im_icc_import_embedded in (to_real intent), is_image in 1761 = error (_ "bad arguments to " ++ "icc_import_embedded") 1762{ 1763 icc_import_embedded_op = Operator "icc_import_embedded" 1764 (icc_import_embedded intent) 1765 Operator_type.COMPOUND_REWRAP false; 1766} 1767 1768icc_transform in_profile out_profile intent in 1769 = oo_unary_function icc_transform_op in, is_class in 1770 = im_icc_transform in 1771 (expand in_profile) (expand out_profile) 1772 (to_real intent), is_image in 1773 = error (_ "bad arguments to " ++ "icc_transform") 1774{ 1775 icc_transform_op = Operator "icc_transform" 1776 (icc_transform in_profile out_profile intent) 1777 Operator_type.COMPOUND_REWRAP false; 1778} 1779 1780icc_ac2rc profile in 1781 = oo_unary_function icc_ac2rc_op in, is_class in 1782 = im_icc_ac2rc in (expand profile), is_image in 1783 = error (_ "bad arguments to " ++ "icc_ac2rc") 1784{ 1785 icc_ac2rc_op = Operator "icc_ac2rc" 1786 (icc_ac2rc profile) 1787 Operator_type.COMPOUND_REWRAP false; 1788} 1789 1790 1791// Given a xywh rect, flip it around so wh are always positive 1792rect_normalise x y w h 1793 = [x', y', w', h'] 1794{ 1795 x' 1796 = x + w, w < 0 1797 = x; 1798 y' 1799 = y + h, h < 0 1800 = y; 1801 w' = abs w; 1802 h' = abs h; 1803} 1804 1805draw_flood_blob x y ink image 1806 = oo_unary_function draw_flood_blob_op image, is_class image 1807 = im_draw_flood_blob image (to_real x) (to_real y) ink, is_image image 1808 = error (_ "bad arguments to " ++ "draw_flood_blob") 1809{ 1810 draw_flood_blob_op = Operator "draw_flood_blob" 1811 (draw_flood_blob x y ink) Operator_type.COMPOUND_REWRAP false; 1812} 1813 1814draw_flood x y ink image 1815 = oo_unary_function draw_flood_op image, is_class image 1816 = im_draw_flood image (to_real x) (to_real y) ink, is_image image 1817 = error (_ "bad arguments to " ++ "draw_flood") 1818{ 1819 draw_flood_op = Operator "draw_flood" 1820 (draw_flood x y ink) Operator_type.COMPOUND_REWRAP false; 1821} 1822 1823/* This version of draw_rect uses insert_noexpand and will be fast, even for 1824 * huge images. 1825 */ 1826draw_rect_width x y w h f t ink image 1827 = oo_unary_function draw_rect_width_op image, is_class image 1828 = my_draw_rect_width image (to_int x) (to_int y) 1829 (to_int w) (to_int h) (to_int f) (to_int t) ink, 1830 is_image image 1831 = error (_ "bad arguments to " ++ "draw_rect_width") 1832{ 1833 draw_rect_width_op = Operator "draw_rect_width" 1834 (draw_rect_width x y w h f t ink) 1835 Operator_type.COMPOUND_REWRAP false; 1836 1837 my_draw_rect_width image x y w h f t ink 1838 = insert x' y' (block w' h') image, f == 1 1839 = (insert x' y' (block w' t) @ 1840 insert (x' + w' - t) y' (block t h') @ 1841 insert x' (y' + h' - t) (block w' t) @ 1842 insert x' y' (block t h')) image 1843 { 1844 insert = insert_noexpand; 1845 block w h = image_new w h (get_bands image) (get_format image) 1846 (get_coding image) (get_type image) ink' 0 0; 1847 ink' 1848 = Vector ink, is_list ink 1849 = ink; 1850 [x', y', w', h'] = rect_normalise x y w h; 1851 } 1852} 1853 1854/* Default to 1 pixel wide edges. 1855 */ 1856draw_rect x y w h f ink image 1857 = draw_rect_width x y w h f 1 ink image; 1858 1859/* This version of draw_rect uses the paintbox rect draw operation. It is an 1860 * inplace operation and will use bucketloads of memory. 1861 */ 1862draw_rect_paintbox x y w h f ink image 1863 = oo_unary_function draw_rect_op image, is_class image 1864 = im_draw_rect image (to_real x) (to_real y) 1865 (to_real w) (to_real h) (to_real f) ink, is_image image 1866 = error (_ "bad arguments to " ++ "draw_rect_paintbox") 1867{ 1868 draw_rect_op = Operator "draw_rect" 1869 (draw_rect x y w h f ink) Operator_type.COMPOUND_REWRAP false; 1870} 1871 1872draw_circle x y r f ink image 1873 = oo_unary_function draw_circle_op image, is_class image 1874 = im_draw_circle image (to_real x) (to_real y) 1875 (to_real r) (to_real f) ink, is_image image 1876 = error (_ "bad arguments to " ++ "draw_circle") 1877{ 1878 draw_circle_op = Operator "draw_circle" 1879 (draw_circle x y r f ink) Operator_type.COMPOUND_REWRAP false; 1880} 1881 1882draw_line x1 y1 x2 y2 ink image 1883 = oo_unary_function draw_line_op image, is_class image 1884 = im_draw_line image (to_real x1) (to_real y1) 1885 (to_real x2) (to_real y2) ink, is_image image 1886 = error (_ "bad arguments to " ++ "draw_line") 1887{ 1888 draw_line_op = Operator "draw_line" 1889 (draw_line x1 y1 x2 y2 ink) Operator_type.COMPOUND_REWRAP false; 1890} 1891 1892print_base base in 1893 = oo_unary_function print_base_op in, is_class in 1894 = map (print_base base) in, is_list in 1895 = print_base_real, is_real in 1896 = error (_ "bad arguments to " ++ "print_base") 1897{ 1898 print_base_op 1899 = Operator "print_base" (print_base base) Operator_type.COMPOUND false; 1900 1901 print_base_real 1902 = error "print_base: bad base", base < 2 || base > 16 1903 = "0", in < 0 || chars == [] 1904 = reverse chars 1905 { 1906 digits = map (\x x % base) 1907 (takewhile (not_equal 0) 1908 (iterate (\x idivide x base) in)); 1909 chars = map tohd digits; 1910 1911 tohd x 1912 = (char) ((int) '0' + x), x < 10 1913 = (char) ((int) 'A' + (x - 10)); 1914 } 1915} 1916 1917/* id x: the identity function 1918 * 1919 * id :: * -> * 1920 */ 1921id x = x; 1922 1923/* const x y: junk y, return x 1924 * 1925 * (const 3) is the function that always returns 3. 1926 * const :: * -> ** -> * 1927 */ 1928const x y = x; 1929 1930/* converse fn a b: swap order of args to fn 1931 * 1932 * converse fn a b == fn b a 1933 * converse :: (* -> ** -> ***) -> ** -> * -> *** 1934 */ 1935converse fn a b = fn b a; 1936 1937/* fix fn x: find the fixed point of a function 1938 */ 1939fix fn x = limit (iterate fn x); 1940 1941/* until pred fn n: apply fn to n until pred succeeds; return that value 1942 * 1943 * until (more 1000) (multiply 2) 1 = 1024 1944 * until :: (* -> bool) -> (* -> *) -> * -> * 1945 */ 1946until pred fn n 1947 = n, pred n 1948 = until pred fn (fn n); 1949 1950/* Infinite list of primes. 1951 */ 1952primes 1953 = 1 : (sieve [2 ..]) 1954{ 1955 sieve l = hd l : sieve (filter (nmultiple (hd l)) (tl l)); 1956 nmultiple n x = x / n != (int) (x / n); 1957} 1958 1959/* Map an n-ary function (pass the args as a list) over groups of objects. 1960 * The objects can be single objects or groups. If more than one 1961 * object is a group, we iterate for the length of the smallest group. 1962 * Don't loop over plain lists, since we want (eg.) "mean [1, 2, 3]" to work. 1963 * Treat [] as no-value --- ie. if any of the inputs are [] we put [] into the 1964 * output and don't apply the function. 1965 1966 copy-pasted into _types, keep in sync 1967 1968 */ 1969map_nary fn args 1970 = fn args, groups == [] 1971 = Group (map process [0, 1 .. shortest - 1]) 1972{ 1973 // find all the group arguments 1974 groups = filter is_Group args; 1975 1976 // what's the length of the shortest group arg? 1977 shortest = foldr1 min_pair (map (len @ get_value) groups); 1978 1979 // process index n ... pull that member from each argument 1980 // recurse to handle application, so we work for nested groups too 1981 process n 1982 = NULL, any (map (is_noval n) args) 1983 = map_nary fn (map (extract n) args) 1984 { 1985 extract n arg 1986 = arg.value?n, is_Group arg 1987 = arg; 1988 1989 is_noval n arg = is_Group arg && arg.value?n == NULL; 1990 } 1991} 1992 1993/* Map a 1-ary function over an object. 1994 */ 1995map_unary fn a = map_nary (list_1ary fn) [a]; 1996 1997/* Map a 2-ary function over a pair of objects. 1998 */ 1999map_binary fn a b = map_nary (list_2ary fn) [a, b]; 2000 2001/* Map a 3-ary function over three objects. 2002 */ 2003map_trinary fn a b c = map_nary (list_3ary fn) [a, b, c]; 2004 2005/* Map a 4-ary function over three objects. 2006 */ 2007map_quaternary fn a b c d = map_nary (list_4ary fn) [a, b, c, d]; 2008 2009/* Same as map_nary, but for lists. Handy for (eg.) implementing arith ops on 2010 * vectors and matricies. 2011 */ 2012map_nary_list fn args 2013 = fn args, lists == [] 2014 = map process [0, 1 .. shortest - 1] 2015{ 2016 // find all the list arguments 2017 lists = filter is_list args; 2018 2019 // what's the length of the shortest list arg? 2020 shortest = foldr1 min_pair (map len lists); 2021 2022 // process index n ... pull that member from each argument 2023 // recurse to handle application, so we work for nested lists too 2024 process n 2025 = map_nary_list fn (map (extract n) args) 2026 { 2027 extract n arg 2028 = arg?n, is_list arg 2029 = arg; 2030 } 2031} 2032 2033map_unaryl fn a = map_nary_list (list_1ary fn) [a]; 2034 2035map_binaryl fn a b = map_nary_list (list_2ary fn) [a, b]; 2036 2037/* Remove features smaller than x pixels across from an image. This used to be 2038 * rather complex ... convsep is now good enough to use. 2039 */ 2040smooth x image = convsep (matrix_gaussian_blur (to_real x * 2)) image; 2041 2042/* Chop up an image into a list of lists of smaller images. Pad edges with 2043 * black. 2044 */ 2045imagearray_chop tile_width tile_height hoverlap voverlap i 2046 = map chop' [0, vstep .. last_y] 2047{ 2048 width = get_width i; 2049 height = get_height i; 2050 bands = get_bands i; 2051 format = get_format i; 2052 type = get_type i; 2053 2054 tile_width' = to_real tile_width; 2055 tile_height' = to_real tile_height; 2056 hoverlap' = to_real hoverlap; 2057 voverlap' = to_real voverlap; 2058 2059 /* Unique pixels per tile. 2060 */ 2061 hstep = tile_width' - hoverlap'; 2062 vstep = tile_height' - voverlap'; 2063 2064 /* Number of tiles across and down. Remember the case where width == 2065 * hstep. 2066 */ 2067 across = (int) ((width - 1) / hstep) + 1; 2068 down = (int) ((height - 1) / vstep) + 1; 2069 2070 /* top/left of final tile. 2071 */ 2072 last_x = hstep * (across - 1); 2073 last_y = vstep * (down - 1); 2074 2075 /* How much do we need to pad by? 2076 */ 2077 sx = last_x + tile_width'; 2078 sy = last_y + tile_height'; 2079 2080 /* Expand image with black to pad size. 2081 */ 2082 pad = embed 0 0 0 sx sy i; 2083 2084 /* Chop up a row. 2085 */ 2086 chop' y 2087 = map chop'' [0, hstep .. last_x] 2088 { 2089 chop'' x = extract_area x y tile_width' tile_height' pad; 2090 } 2091} 2092 2093/* Reassemble image. 2094 */ 2095imagearray_assemble hoverlap voverlap il 2096 = (image_set_origin 0 0 @ foldl1 tbj @ map (foldl1 lrj)) il 2097{ 2098 lrj l r = insert (get_width l + hoverlap) 0 r l; 2099 tbj t b = insert 0 (get_height t + voverlap) b t; 2100} 2101 2102/* Generate an nxn identity matrix. 2103 */ 2104identity_matrix n 2105 = error "identity_matrix: n > 0", n < 1 2106 = map line [0 .. n - 1] 2107{ 2108 line p = take p [0, 0 ..] ++ [1] ++ take (n - p - 1) [0, 0 ..]; 2109} 2110 2111/* Incomplete gamma function Q(a, x) == 1 - P(a, x) 2112 2113 FIXME ... this is now a builtin, until we can get a nice List class 2114 (requires overloadable ':') 2115 2116gammq a x 2117 = error "bad args", x < 0 || a <= 0 2118 = 1 - gamser, x < a + 1 2119 = gammcf 2120{ 2121 gamser = (gser a x)?0; 2122 gammcf = (gcf a x)?0; 2123} 2124 */ 2125 2126/* Incomplete gamma function P(a, x) evaluated as series representation. Also 2127 * return ln(gamma(a)) ... nr in c, pp 218 2128 */ 2129gser a x 2130 = [gamser, gln] 2131{ 2132 gln = gammln a; 2133 gamser 2134 = error "bad args", x < 0 2135 = 0, x == 0 2136 = 1 // fix this 2137 { 2138 // maximum iterations 2139 maxit = 100; 2140 2141 ap = List [a + 1, a + 2 ...]; 2142 xoap = x / ap; 2143 2144 del = map product (prefixes xoap.value); 2145 2146 2147 2148 2149 2150 2151 2152/* 2153 del = map (multiply (1 / a)) (map product (prefixes xoap)) 2154 2155 del = 2156 2157 xap = iterate (multiply 2158 */ 2159 2160 /* Generate all prefixes of a list ... [1,2,3] -> [[1], [1, 2], [1, 2, 2161 * 3], [1, 2, 3, 4] ...] 2162 */ 2163 prefixes l = map (converse take l) [1..]; 2164 2165 } 2166} 2167 2168/* ln(gamma(xx)) ... nr in c, pp 214 2169 */ 2170gammln xx 2171 = gln 2172{ 2173 cof = [76.18009172947146, -86.50532032941677, 24.01409824083091, 2174 -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5]; 2175 y = take 6 (iterate (add 1) (xx + 1)); 2176 ser = 1.000000000190015 + sum (map2 divide cof y); 2177 tmp = xx + 0.5; 2178 tmp' = tmp - ((xx + 0.5) * log tmp); 2179 gln = -tmp + log (2.5066282746310005 * ser / xx); 2180} 2181 2182/* make a LUT from a scatter 2183 */ 2184buildlut x 2185 = Plot_histogram (im_buildlut x), is_Matrix x && x.width > 1 2186 = im_buildlut (Matrix x), is_matrix x && is_list_len_more 1 x?0 2187 = error (_ "bad arguments to " ++ "buildlut"); 2188 2189/* Linear regression. Return a class with the stuff we need in. 2190 * from s15.2, p 665 NR in C 2191 * 2192 * Also calculate R2, see eg.: 2193 * https://en.wikipedia.org/wiki/Coefficient_of_determination 2194 */ 2195linreg xes yes 2196 = obj 2197{ 2198 obj = class { 2199 // in case we ever get shown in the workspace 2200 _vislevel = 2; 2201 2202 slope = sum [t * y :: [t, y] <- zip2 tes yes] / st2; 2203 intercept = (sy - sx * slope) / ss; 2204 2205 chi2 = sum [(y - intercept - slope * x) ** 2 :: [x, y] <- zip2 xes yes]; 2206 2207 siga = (chi2 / (ss - 2)) ** 0.5 * 2208 ((1 + sx ** 2 / (ss * st2)) / ss) ** 0.5; 2209 sigb = (chi2 / (ss - 2)) ** 0.5 * (1 / st2) ** 0.5; 2210 2211 // for compat with linregw, see below 2212 q = 1.0; 2213 2214 R2 = 1 - chi2 / sum [(y - my) ** 2 :: y <- yes]; 2215 } 2216 2217 ss = len xes; 2218 sx = sum xes; 2219 sy = sum yes; 2220 my = sy / ss; 2221 sxoss = sx / ss; 2222 2223 tes = [x - sxoss :: x <- xes]; 2224 st2 = sum [t ** 2 :: t <- tes]; 2225} 2226 2227/* Weighted linear regression. Xes, yes and a list of deviations. 2228 */ 2229linregw xes yes devs 2230 = obj 2231{ 2232 obj = class { 2233 // in case we ever get shown in the workspace 2234 _vislevel = 2; 2235 2236 slope = sum [(t * y) / sd :: [t, y, sd] <- zip3 tes yes devs] / st2; 2237 intercept = (sy - sx * slope) / ss; 2238 2239 chi2 = sum [((y - intercept - slope * x) / sd) ** 2 :: 2240 [x, y, sd] <- zip3 xes yes devs]; 2241 2242 siga = ((1 + sx * sx / (ss * st2)) / ss) ** 0.5; 2243 sigb = (1 / st2) ** 0.5; 2244 2245 q = gammq (0.5 * (len xes - 2)) (0.5 * chi2); 2246 2247 R2 = 1 - chi2 / sum [(y - my) ** 2 :: y <- yes]; 2248 } 2249 2250 wt = [sd ** -0.5 :: sd <- devs]; 2251 2252 ss = sum wt; 2253 sx = sum [x * w :: [x, w] <- zip2 xes wt]; 2254 sy = sum [y * w :: [y, w] <- zip2 yes wt]; 2255 my = sy / len xes; 2256 sxoss = sx / ss; 2257 2258 tes = [(x - sxoss) / sd :: [x, sd] <- zip2 xes devs]; 2259 st2 = sum [t ** 2 :: t <- tes]; 2260} 2261 2262/* Clustering: pass in a list of points, repeatedly merge the 2263 * closest two points until no two points are closer than the threshold. 2264 * Return [merged-points, corresponding-weights]. A weight is a list of the 2265 * indexes we merged to make that point, ie. len weight == how significant 2266 * this point is. 2267 * 2268 * eg. 2269 * cluster 12 [152,154,155,42,159] == 2270 * [[155,42],[[1,2,0,4],[3]]] 2271 */ 2272cluster thresh points 2273 = oo_unary_function cluster_op points, is_class points 2274 // can't use [0..len points - 1], in case len points == 0 2275 = merge [points, map (converse cons []) (take (len points) [0 ..])], 2276 is_list points 2277 = error (_ "bad arguments to " ++ "cluster") 2278{ 2279 cluster_op = Operator "cluster" 2280 (cluster thresh) Operator_type.COMPOUND false; 2281 2282 merge x 2283 = x, m < 2 || d > thresh 2284 = merge [points', weights'] 2285 { 2286 [points, weights] = x; 2287 m = len points; 2288 2289 // generate indexes of all possible pairs, avoiding comparing a thing 2290 // to itself, and assuming that dist is reflexive 2291 // first index is always less than 2nd index 2292 // the +1,+2 makes sure we have an increasing generator, otherwise we 2293 // can get [3 .. 4] (for example), which will make a decreasing 2294 // sequence 2295 pairs = [[x, y] :: x <- [0 .. m - 1]; y <- [x + 1, x + 2 .. m - 1]]; 2296 2297 // distance function 2298 // arg is eg. [3,1], meaning get distance from point 3 to point 1 2299 dist x 2300 = abs (points?i - points?j) 2301 { 2302 [i, j] = x; 2303 } 2304 2305 // smallest distance, then the two points we merge 2306 p = minpos (map dist pairs); 2307 d = dist pairs?p; 2308 [i, j] = pairs?p; 2309 2310 // new point and new weight 2311 nw = weights?i ++ weights?j; 2312 np = (points?i * len weights?i + points?j * len weights?j) / len nw; 2313 2314 // remove element i from a list 2315 remove i l = take i l ++ drop (i + 1) l; 2316 2317 // remove two old points, add the new merged one 2318 // i < j (see "pairs", above) 2319 points' = np : remove i (remove j points); 2320 weights' = nw : remove i (remove j weights); 2321 } 2322} 2323 2324/* Extract the area of an image around an arrow. 2325 * Transform the image to make the arrow horizontal, then displace by hd and 2326 * vd pxels, then cut out a bit h pixels high centered on the arrow. 2327 */ 2328extract_arrow hd vd h arrow 2329 = extract_area (re p' + hd) (im p' - h / 2 + vd) (re pv) h im' 2330{ 2331 // the line as a polar vector 2332 pv = polar (arrow.width, arrow.height); 2333 a = im pv; 2334 2335 // smallest rotation that will make the line horizontal 2336 a' 2337 = 360 - a, a > 270 2338 = 180 - a, a > 90 2339 = -a; 2340 2341 im' = rotate Interpolate_bilinear a' arrow.image; 2342 2343 // look at the start and end of the arrow, pick the leftmost 2344 p 2345 = (arrow.left, arrow.top), arrow.left <= arrow.right 2346 = (arrow.right, arrow.bottom); 2347 2348 // transform that point to im' space 2349 p' = rectangular (polar p + (0, a')) + (im'.xoffset, im'.yoffset); 2350} 2351 2352/* You'd think these would go in _convert, but they are not really colour ops, 2353 * so put them here. 2354 */ 2355rad2float image 2356 = oo_unary_function rad2float_op image, is_class image 2357 = im_rad2float image, is_image image 2358 = error (_ "bad arguments to " ++ "rad2float") 2359{ 2360 rad2float_op = Operator "rad2float" 2361 rad2float Operator_type.COMPOUND_REWRAP false; 2362} 2363 2364float2rad image 2365 = oo_unary_function float2rad_op image, is_class image 2366 = im_float2rad image, is_image image 2367 = error (_ "bad arguments to " ++ "float2rad") 2368{ 2369 float2rad_op = Operator "float2rad" 2370 float2rad Operator_type.COMPOUND_REWRAP false; 2371} 2372 2373segment x 2374 = oo_unary_function segment_op x, is_class x 2375 = image', is_image x 2376 = error (_ "bad arguments to " ++ "segment") 2377{ 2378 segment_op = Operator "segment" 2379 segment Operator_type.COMPOUND_REWRAP false; 2380 2381 [image, nsegs] = im_segment x; 2382 image' = im_copy_set_meta image "n-segments" nsegs; 2383} 2384 2385point a b 2386 = oo_binary_function point_op a b, is_class a 2387 = oo_binary'_function point_op a b, is_class b 2388 = im_read_point b x y, is_image b 2389 = [b?x?y], is_matrix b 2390 = [b?x], is_real_list b && y == 0 2391 = [b?y], is_real_list b && x == 0 2392 = error (_ "bad arguments to " ++ "point") 2393{ 2394 point_op = Operator "point" 2395 (\a\b Vector (point a b)) Operator_type.COMPOUND false; 2396 2397 (x, y) 2398 = a, is_complex a; 2399 = (a?0, a?1), is_real_list a && is_list_len 2 a 2400 = error "bad position format"; 2401} 2402 2403/* One image in, one out. 2404 */ 2405system_image command x 2406 = oo_unary_function system_image_op x, is_class x 2407 = system x, is_image x 2408 = error (_ "bad arguments to " ++ "system_image") 2409{ 2410 system_image_op = Operator "system_image" 2411 (system_image command) Operator_type.COMPOUND_REWRAP false; 2412 2413 system im 2414 = image_out 2415 { 2416 [image_out, log] = 2417 im_system_image (get_image im) "%s.tif" "%s.tif" command; 2418 } 2419} 2420 2421/* Two images in, one out. 2422 */ 2423system_image2 command x1 x2 2424 = oo_binary_function system_image2_op x1 x2, is_class x1 2425 = oo_binary'_function system_image2_op x1 x2, is_class x2 2426 = system x1 x2, is_image x1 && is_image x2 2427 = error (_ "bad arguments to " ++ "system_image2") 2428{ 2429 system_image2_op = Operator "system_image2" 2430 (system_image2 command) Operator_type.COMPOUND_REWRAP false; 2431 2432 system x1 x2 2433 = image_out 2434 { 2435 [image_out] = vips_call "system" [command] [ 2436 $in => [x1, x2], 2437 $out => true, 2438 $out_format => "%s.tif", 2439 $in_format => "%s.tif" 2440 ]; 2441 } 2442} 2443 2444/* Three images in, one out. 2445 */ 2446system_image3 command x1 x2 x3 2447 = oo_binary_function system_image2_op x2 x3, is_class x2 2448 = oo_binary'_function system_image2_op x2 x3, is_class x3 2449 = system x1 x2 x3, is_image x1 && is_image x2 && is_image x3 2450 = error (_ "bad arguments to " ++ "system_image3") 2451{ 2452 system_image2_op = Operator "system_image2" 2453 (system_image3 command x1) Operator_type.COMPOUND_REWRAP false; 2454 2455 system x1 x2 x3 2456 = image_out 2457 { 2458 [image_out] = vips_call "system" [command] [ 2459 $in => [x1, x2, x3], 2460 $out => true, 2461 $out_format => "%s.tif", 2462 $in_format => "%s.tif" 2463 ]; 2464 } 2465} 2466 2467/* Zero images in, one out. 2468 */ 2469system_image0 command 2470 = Image image_out 2471{ 2472 [image_out] = vips_call "system" [command] [ 2473 $out => true, 2474 $out_format => "%s.tif" 2475 ]; 2476} 2477 2478hough_line w h x 2479 = oo_unary_function hough_line_op x, is_class x 2480 = hline (to_real w) (to_real h) x, is_image x 2481 = error (_ "bad arguments to " ++ "hough_line") 2482{ 2483 hough_line_op = Operator "hough_line" 2484 (hough_line w h) Operator_type.COMPOUND_REWRAP false; 2485 2486 hline w h x 2487 = pspace 2488 { 2489 [pspace] = vips_call "hough_line" [x] [ 2490 $width => w, 2491 $height => h 2492 ]; 2493 } 2494} 2495 2496hough_circle s mn mx x 2497 = oo_unary_function hough_circle_op x, is_class x 2498 = hcircle (to_real s) (to_real mn) (to_real mx) x, is_image x 2499 = error (_ "bad arguments to " ++ "hough_circle") 2500{ 2501 hough_circle_op = Operator "hough_circle" 2502 (hough_circle s mn mx) Operator_type.COMPOUND_REWRAP false; 2503 2504 hcircle s mn mx x 2505 = pspace 2506 { 2507 [pspace] = vips_call "hough_circle" [x] [ 2508 $scale => s, 2509 $min_radius => mn, 2510 $max_radius => mx 2511 ]; 2512 } 2513} 2514 2515mapim interp ind in 2516 = oo_binary_function mapim_op ind in, is_class ind 2517 = oo_binary'_function mapim_op ind in, is_class in 2518 = mapim_fn ind in, is_image ind && is_image in 2519 = error (_ "bad arguments to " ++ "mapim") 2520{ 2521 mapim_op = Operator "mapim" 2522 (mapim interp) Operator_type.COMPOUND_REWRAP false; 2523 2524 mapim_fn ind im 2525 = out 2526 { 2527 [out] = vips_call "mapim" [im, ind] [$interpolate => interp]; 2528 } 2529} 2530 2531perlin cell width height 2532 = Image im 2533{ 2534 [im] = vips_call "perlin" [to_real width, to_real height] [ 2535 $cell_size => to_real cell 2536 ]; 2537} 2538 2539worley cell width height 2540 = Image im 2541{ 2542 [im] = vips_call "worley" [to_real width, to_real height] [ 2543 $cell_size => to_real cell 2544 ]; 2545} 2546 2547gaussnoise width height mean sigma 2548 = im 2549{ 2550 [im] = vips_call "gaussnoise" [to_real width, to_real height] [ 2551 $mean => to_real mean, 2552 $sigma => to_real sigma 2553 ]; 2554} 2555 2556flattenimage bg x 2557 = oo_unary_function flatten_op x, is_class x 2558 = flt (to_vector bg) x, is_image x 2559 = error (_ "bad arguments to " ++ "flattenimage") 2560{ 2561 flatten_op = Operator "flatten" 2562 (flattenimage bg) Operator_type.COMPOUND_REWRAP false; 2563 2564 flt bg x 2565 = out 2566 { 2567 [out] = vips_call "flatten" [x] [ 2568 $background => bg.value 2569 ]; 2570 } 2571} 2572 2573premultiply x 2574 = oo_unary_function premultiply_op x, is_class x 2575 = prem x, is_image x 2576 = error (_ "bad arguments to " ++ "premultiply") 2577{ 2578 premultiply_op = Operator "premultiply" 2579 premultiply Operator_type.COMPOUND_REWRAP false; 2580 2581 prem x 2582 = out 2583 { 2584 [out] = vips_call "premultiply" [x] [ 2585 ]; 2586 } 2587} 2588 2589unpremultiply x 2590 = oo_unary_function unpremultiply_op x, is_class x 2591 = unprem x, is_image x 2592 = error (_ "bad arguments to " ++ "unpremultiply") 2593{ 2594 unpremultiply_op = Operator "unpremultiply" 2595 unpremultiply Operator_type.COMPOUND_REWRAP false; 2596 2597 unprem x 2598 = out 2599 { 2600 [out] = vips_call "unpremultiply" [x] [ 2601 ]; 2602 } 2603} 2604 2605hist_entropy x 2606 = oo_unary_function hist_entropy_op x, is_class x 2607 = entropy x, is_image x 2608 = error (_ "bad arguments to " ++ "hist_entropy") 2609{ 2610 hist_entropy_op = Operator "hist_entropy" 2611 hist_entropy Operator_type.COMPOUND_REWRAP false; 2612 2613 entropy x 2614 = out 2615 { 2616 [out] = vips_call "hist_entropy" [x] [ 2617 ]; 2618 } 2619} 2620 2621canny sigma precision x 2622 = oo_unary_function canny_op x, is_class x 2623 = canny_im (to_real sigma) (to_int precision) x, is_image x 2624 = error (_ "bad arguments to " ++ "canny") 2625{ 2626 canny_op = Operator "canny" 2627 (canny sigma precision) Operator_type.COMPOUND_REWRAP false; 2628 2629 canny_im sigma precision x 2630 = out 2631 { 2632 [out] = vips_call "canny" [x] [ 2633 $sigma => sigma, 2634 $precision => precision 2635 ]; 2636 } 2637 2638} 2639 2640sobel x 2641 = oo_unary_function sobel_op x, is_class x 2642 = sobel_im x, is_image x 2643 = error (_ "bad arguments to " ++ "sobel") 2644{ 2645 sobel_op = Operator "sobel" 2646 sobel Operator_type.COMPOUND_REWRAP false; 2647 2648 sobel_im x 2649 = out 2650 { 2651 [out] = vips_call "sobel" [x] [ 2652 ]; 2653 } 2654} 2655